CppADCodeGen  HEAD
A C++ Algorithmic Differentiation Package with Source Code Generation
model_c_source_gen_jac.hpp
1 #ifndef CPPAD_CG_MODEL_C_SOURCE_GEN_JAC_INCLUDED
2 #define CPPAD_CG_MODEL_C_SOURCE_GEN_JAC_INCLUDED
3 /* --------------------------------------------------------------------------
4  * CppADCodeGen: C++ Algorithmic Differentiation with Source Code Generation:
5  * Copyright (C) 2012 Ciengis
6  *
7  * CppADCodeGen is distributed under multiple licenses:
8  *
9  * - Eclipse Public License Version 1.0 (EPL1), and
10  * - GNU General Public License Version 3 (GPL3).
11  *
12  * EPL1 terms and conditions can be found in the file "epl-v10.txt", while
13  * terms and conditions for the GPL3 can be found in the file "gpl3.txt".
14  * ----------------------------------------------------------------------------
15  * Author: Joao Leal
16  */
17 
18 namespace CppAD {
19 namespace cg {
20 
21 template<class Base>
22 void ModelCSourceGen<Base>::generateJacobianSource() {
23  using std::vector;
24 
25  const std::string jobName = "Jacobian";
26 
27  startingJob("'" + jobName + "'", JobTimer::GRAPH);
28 
29  CodeHandler<Base> handler;
30  handler.setJobTimer(_jobTimer);
31 
32  vector<CGBase> indVars(_fun.Domain());
33  handler.makeVariables(indVars);
34  if (_x.size() > 0) {
35  for (size_t i = 0; i < indVars.size(); i++) {
36  indVars[i].setValue(_x[i]);
37  }
38  }
39 
40  size_t m = _fun.Range();
41  size_t n = _fun.Domain();
42 
43  vector<CGBase> jac(n * m);
44  if (_jacMode == JacobianADMode::Automatic) {
45  jac = _fun.Jacobian(indVars);
46  } else if (_jacMode == JacobianADMode::Forward) {
47  JacobianFor(_fun, indVars, jac);
48  } else {
49  JacobianRev(_fun, indVars, jac);
50  }
51 
52  finishedJob();
53 
54  LanguageC<Base> langC(_baseTypeName);
55  langC.setMaxAssignmentsPerFunction(_maxAssignPerFunc, &_sources);
56  langC.setMaxOperationsPerAssignment(_maxOperationsPerAssignment);
57  langC.setParameterPrecision(_parameterPrecision);
58  langC.setGenerateFunction(_name + "_" + FUNCTION_JACOBIAN);
59 
60  std::ostringstream code;
61  std::unique_ptr<VariableNameGenerator<Base> > nameGen(createVariableNameGenerator("jac"));
62 
63  handler.generateCode(code, langC, jac, *nameGen, _atomicFunctions, jobName);
64 }
65 
66 template<class Base>
67 void ModelCSourceGen<Base>::generateSparseJacobianSource(MultiThreadingType multiThreadingType) {
68  size_t m = _fun.Range();
69  size_t n = _fun.Domain();
70 
74  determineJacobianSparsity();
75 
76  bool forwardMode;
77 
78  if (_jacMode == JacobianADMode::Automatic) {
79  if (_custom_jac.defined) {
80  forwardMode = estimateBestJacobianADMode(_jacSparsity.rows, _jacSparsity.cols);
81  } else {
82  forwardMode = n <= m;
83  }
84  } else {
85  forwardMode = _jacMode == JacobianADMode::Forward;
86  }
87 
91  if (_sparseJacobianReusesOne && _forwardOne && forwardMode) {
92  generateSparseJacobianForRevSource(true, multiThreadingType);
93  } else if (_sparseJacobianReusesOne && _reverseOne && !forwardMode) {
94  generateSparseJacobianForRevSource(false, multiThreadingType);
95  } else {
96  generateSparseJacobianSource(forwardMode);
97  }
98 }
99 
100 template<class Base>
101 void ModelCSourceGen<Base>::generateSparseJacobianSource(bool forward) {
102  using std::vector;
103 
104  const std::string jobName = "sparse Jacobian";
105 
106  //size_t m = _fun.Range();
107  size_t n = _fun.Domain();
108 
109  startingJob("'" + jobName + "'", JobTimer::GRAPH);
110 
111  CodeHandler<Base> handler;
112  handler.setJobTimer(_jobTimer);
113 
114  vector<CGBase> indVars(n);
115  handler.makeVariables(indVars);
116  if (_x.size() > 0) {
117  for (size_t i = 0; i < n; i++) {
118  indVars[i].setValue(_x[i]);
119  }
120  }
121 
122  vector<CGBase> jac(_jacSparsity.rows.size());
123  if (_loopTapes.empty()) {
124  //printSparsityPattern(_jacSparsity.sparsity, "jac sparsity");
125  CppAD::sparse_jacobian_work work;
126  if (forward) {
127  _fun.SparseJacobianForward(indVars, _jacSparsity.sparsity, _jacSparsity.rows, _jacSparsity.cols, jac, work);
128  } else {
129  _fun.SparseJacobianReverse(indVars, _jacSparsity.sparsity, _jacSparsity.rows, _jacSparsity.cols, jac, work);
130  }
131 
132  } else {
133  jac = prepareSparseJacobianWithLoops(handler, indVars, forward);
134  }
135 
136  finishedJob();
137 
138  LanguageC<Base> langC(_baseTypeName);
139  langC.setMaxAssignmentsPerFunction(_maxAssignPerFunc, &_sources);
140  langC.setMaxOperationsPerAssignment(_maxOperationsPerAssignment);
141  langC.setParameterPrecision(_parameterPrecision);
142  langC.setGenerateFunction(_name + "_" + FUNCTION_SPARSE_JACOBIAN);
143 
144  std::ostringstream code;
145  std::unique_ptr<VariableNameGenerator<Base> > nameGen(createVariableNameGenerator("jac"));
146 
147  handler.generateCode(code, langC, jac, *nameGen, _atomicFunctions, jobName);
148 }
149 
150 template<class Base>
151 void ModelCSourceGen<Base>::generateSparseJacobianForRevSource(bool forward,
152  MultiThreadingType multiThreadingType) {
153  //size_t m = _fun.Range();
154  //size_t n = _fun.Domain();
155  using namespace std;
156 
157  std::map<size_t, CompressedVectorInfo> jacInfo;
158  string functionRevFor, revForSuffix;
159  if (forward) {
160  // jacInfo[var].index{equations}
161  for (size_t e = 0; e < _jacSparsity.rows.size(); e++) {
162  jacInfo[_jacSparsity.cols[e]].indexes.push_back(_jacSparsity.rows[e]);
163  }
164  for (auto& it : jacInfo) {
165  size_t col = it.first;
166  it.second.locations = determineOrderByCol(col, it.second.indexes, _jacSparsity.rows, _jacSparsity.cols);
167  }
168  _cache.str("");
169  _cache << _name << "_" << FUNCTION_SPARSE_FORWARD_ONE;
170  functionRevFor = _cache.str();
171  revForSuffix = "indep";
172  } else {
173  // jacInfo[equation].index{vars}
174  for (size_t e = 0; e < _jacSparsity.rows.size(); e++) {
175  jacInfo[_jacSparsity.rows[e]].indexes.push_back(_jacSparsity.cols[e]);
176  }
177  for (auto& it : jacInfo) {
178  size_t row = it.first;
179  it.second.locations = determineOrderByRow(row, it.second.indexes, _jacSparsity.rows, _jacSparsity.cols);
180  }
181  _cache.str("");
182  _cache << _name << "_" << FUNCTION_SPARSE_REVERSE_ONE;
183  functionRevFor = _cache.str();
184  revForSuffix = "dep";
185  }
186 
187 
192  for (auto& it : jacInfo) {
193  const std::vector<size_t>& els = it.second.indexes;
194  const std::vector<set<size_t> >& location = it.second.locations;
195  CPPADCG_ASSERT_UNKNOWN(els.size() == location.size());
196  CPPADCG_ASSERT_UNKNOWN(els.size() > 0);
197 
198  bool passed = true;
199  size_t jacArrayStart = *location[0].begin();
200  for (size_t e = 0; e < els.size(); e++) {
201  if (location[e].size() > 1) {
202  passed = false; // too many elements
203  break;
204  }
205  if (*location[e].begin() != jacArrayStart + e) {
206  passed = false; // wrong order
207  break;
208  }
209  }
210  it.second.ordered = passed;
211  }
212 
213  size_t maxCompressedSize = 0;
214  map<size_t, bool>::const_iterator itOrd;
215  map<size_t, std::vector<size_t> >::const_iterator it;
216  for (const auto& it : jacInfo) {
217  if (it.second.indexes.size() > maxCompressedSize && !it.second.ordered)
218  maxCompressedSize = it.second.indexes.size();
219  }
220 
221  if (!_loopTapes.empty()) {
225  if (forward) {
226  generateSparseJacobianWithLoopsSourceFromForRev(jacInfo, maxCompressedSize,
227  FUNCTION_SPARSE_FORWARD_ONE, "indep", "jcol",
228  _nonLoopFor1Elements, _loopFor1Groups,
229  generateFunctionNameLoopFor1);
230  } else {
231  generateSparseJacobianWithLoopsSourceFromForRev(jacInfo, maxCompressedSize,
232  FUNCTION_SPARSE_REVERSE_ONE, "dep", "jrow",
233  _nonLoopRev1Elements, _loopRev1Groups,
234  generateFunctionNameLoopRev1);
235  }
236  return;
237  }
238 
239  _cache.str("");
240  _cache << _name << "_" << FUNCTION_SPARSE_JACOBIAN;
241  string functionName(_cache.str());
242 
243  if(!_multiThreading || multiThreadingType == MultiThreadingType::NONE) {
244  _sources[functionName + ".c"] = generateSparseJacobianForRevSingleThreadSource(functionName, jacInfo, maxCompressedSize, functionRevFor, revForSuffix, forward);
245  } else {
246  _sources[functionName + ".c"] = generateSparseJacobianForRevMultiThreadSource(functionName, jacInfo, maxCompressedSize, functionRevFor, revForSuffix, forward, multiThreadingType);
247  }
248 
249  _cache.str("");
250 }
251 
252 template<class Base>
253 std::string ModelCSourceGen<Base>::generateSparseJacobianForRevSingleThreadSource(const std::string& functionName,
254  std::map<size_t, CompressedVectorInfo> jacInfo,
255  size_t maxCompressedSize,
256  const std::string& functionRevFor,
257  const std::string& revForSuffix,
258  bool forward) {
259 
260  LanguageC<Base> langC(_baseTypeName);
261  std::string argsDcl = langC.generateDefaultFunctionArgumentsDcl();
262  std::vector<std::string> argsDcl2 = langC.generateDefaultFunctionArgumentsDcl2();
263 
264  _cache.str("");
265  _cache << "#include <stdlib.h>\n"
266  "\n"
267  << LanguageC<Base>::ATOMICFUN_STRUCT_DEFINITION << "\n\n";
268  generateFunctionDeclarationSource(_cache, functionRevFor, revForSuffix, jacInfo, argsDcl);
269  _cache << "\n";
270  LanguageC<Base>::printFunctionDeclaration(_cache, "void", functionName, argsDcl2);
271  _cache << " {\n"
272  " " << _baseTypeName << " const * inLocal[2];\n"
273  " " << _baseTypeName << " inLocal1 = 1;\n"
274  " " << _baseTypeName << " * outLocal[1];\n"
275  " " << _baseTypeName << " compressed[" << maxCompressedSize << "];\n"
276  " " << _baseTypeName << " * jac = out[0];\n"
277  "\n"
278  " inLocal[0] = in[0];\n"
279  " inLocal[1] = &inLocal1;\n"
280  " outLocal[0] = compressed;\n";
281 
282  langC.setArgumentIn("inLocal");
283  langC.setArgumentOut("outLocal");
284  std::string argsLocal = langC.generateDefaultFunctionArguments();
285 
286  bool previousCompressed = true;
287  for (const auto& it : jacInfo) {
288  size_t index = it.first;
289  const std::vector<size_t>& els = it.second.indexes;
290  const std::vector<std::set<size_t> >& location = it.second.locations;
291  CPPADCG_ASSERT_UNKNOWN(els.size() == location.size());
292 
293  _cache << "\n";
294  bool compressed = !it.second.ordered;
295  if (!compressed) {
296  _cache << " outLocal[0] = &jac[" << *location[0].begin() << "];\n";
297  } else if (!previousCompressed) {
298  _cache << " outLocal[0] = compressed;\n";
299  }
300  _cache << " " << functionRevFor << "_" << revForSuffix << index << "(" << argsLocal << ");\n";
301  if (compressed) {
302  for (size_t e = 0; e < els.size(); e++) {
303  _cache << " ";
304  for (size_t itl : location[e]) {
305  _cache << "jac[" << (itl) << "] = ";
306  }
307  _cache << "compressed[" << e << "];\n";
308  }
309  }
310  previousCompressed = compressed;
311  }
312 
313  _cache << "\n"
314  "}\n";
315 
316  return _cache.str();
317 }
318 
319 template<class Base>
320 std::string ModelCSourceGen<Base>::generateSparseJacobianForRevMultiThreadSource(const std::string& functionName,
321  std::map<size_t, CompressedVectorInfo> jacInfo,
322  size_t maxCompressedSize,
323  const std::string& functionRevFor,
324  const std::string& revForSuffix,
325  bool forward,
326  MultiThreadingType multiThreadingType) {
327  LanguageC<Base> langC(_baseTypeName);
328  std::string argsDcl = langC.generateDefaultFunctionArgumentsDcl();
329  std::vector<std::string> argsDcl2 = langC.generateDefaultFunctionArgumentsDcl2();
330 
331  _cache.str("");
332  _cache << "#include <stdlib.h>\n"
333  "\n"
334  << LanguageC<Base>::ATOMICFUN_STRUCT_DEFINITION << "\n\n";
335  generateFunctionDeclarationSource(_cache, functionRevFor, revForSuffix, jacInfo, argsDcl);
336 
337  langC.setArgumentIn("inLocal");
338  langC.setArgumentOut("outLocal");
339  std::string argsLocal = langC.generateDefaultFunctionArguments();
340 
344  for (const auto& it : jacInfo) {
345  size_t index = it.first;
346  const std::vector<size_t>& els = it.second.indexes;
347  const std::vector<std::set<size_t> >& location = it.second.locations;
348  CPPADCG_ASSERT_UNKNOWN(els.size() == location.size());
349 
350  bool compressed = !it.second.ordered;
351  if (!compressed) {
352  continue;
353  }
354 
355  std::string functionNameWrap = functionRevFor + "_" + revForSuffix + std::to_string(index) + "_wrap";
356  LanguageC<Base>::printFunctionDeclaration(_cache, "void", functionNameWrap, argsDcl2);
357  _cache << " {\n"
358  " " << _baseTypeName << " const * inLocal[2];\n"
359  " " << _baseTypeName << " inLocal1 = 1;\n"
360  " " << _baseTypeName << " * outLocal[1];\n"
361  " " << _baseTypeName << " compressed[" << it.second.indexes.size() << "];\n"
362  " " << _baseTypeName << " * jac = out[0];\n"
363  "\n"
364  " inLocal[0] = in[0];\n"
365  " inLocal[1] = &inLocal1;\n"
366  " outLocal[0] = compressed;\n";
367 
368  _cache << " " << functionRevFor << "_" << revForSuffix << index << "(" << argsLocal << ");\n";
369  for (size_t e = 0; e < els.size(); e++) {
370  _cache << " ";
371  for (size_t itl : location[e]) {
372  _cache << "jac[" << (itl) << "] = ";
373  }
374  _cache << "compressed[" << e << "];\n";
375  }
376  _cache << "}\n";
377  }
378 
379  _cache << "\n"
380  "typedef void (*cppadcg_function_type) (" << argsDcl << ");\n";
381 
385  if(multiThreadingType == MultiThreadingType::OPENMP) {
386  _cache << "\n";
387  printFileStartOpenMP(_cache);
388  _cache << "\n";
389 
390  } else {
391  assert(multiThreadingType == MultiThreadingType::PTHREADS);
392 
393  printFileStartPThreads(_cache, _baseTypeName);
394  }
395 
399  _cache << "\n"
400  "void " << functionName << "(" << argsDcl << ") {\n"
401  " static const cppadcg_function_type p[" << jacInfo.size() << "] = {";
402  for (const auto& it : jacInfo) {
403  size_t index = it.first;
404  if (index != jacInfo.begin()->first) _cache << ", ";
405  if (it.second.ordered) {
406  _cache << functionRevFor << "_" << revForSuffix << index;
407  } else {
408  _cache << functionRevFor << "_" << revForSuffix << index << "_wrap";
409  }
410  }
411  _cache << "};\n"
412  " static const long offset["<< jacInfo.size() <<"] = {";
413  for (const auto& it : jacInfo) {
414  if (it.first != jacInfo.begin()->first) _cache << ", ";
415  if (it.second.ordered) {
416  _cache << *it.second.locations[0].begin();
417  } else {
418  _cache << "0";
419  }
420  }
421  _cache << "};\n"
422  " " << _baseTypeName << " inLocal1 = 1;\n"
423  " " << _baseTypeName << " const * inLocal[2] = {in[0], &inLocal1};\n"
424  " " << _baseTypeName << " * outLocal[1];\n"
425  " " << _baseTypeName << " * jac = out[0];\n"
426  " long i;\n"
427  "\n";
428 
429  if(multiThreadingType == MultiThreadingType::OPENMP) {
430  printFunctionStartOpenMP(_cache, jacInfo.size());
431  _cache << "\n";
432  printLoopStartOpenMP(_cache, jacInfo.size());
433  _cache << " outLocal[0] = &jac[offset[i]];\n"
434  " (*p[i])(" << argsLocal << ");\n";
435  printLoopEndOpenMP(_cache, jacInfo.size());
436  _cache << "\n";
437 
438  } else {
439  assert(multiThreadingType == MultiThreadingType::PTHREADS);
440 
441  printFunctionStartPThreads(_cache, jacInfo.size());
442  _cache << "\n"
443  " for(i = 0; i < " << jacInfo.size() << "; ++i) {\n"
444  " args[i] = (ExecArgStruct*) malloc(sizeof(ExecArgStruct));\n"
445  " args[i]->func = p[i];\n"
446  " args[i]->in = inLocal;\n"
447  " args[i]->out[0] = &jac[offset[i]];\n"
448  " args[i]->atomicFun = " << langC.getArgumentAtomic() << ";\n"
449  " }\n"
450  "\n";
451  printFunctionEndPThreads(_cache, jacInfo.size());
452  }
453 
454  _cache << "\n"
455  "}\n";
456 
457  return _cache.str();
458 }
459 
460 template<class Base>
461 void ModelCSourceGen<Base>::determineJacobianSparsity() {
462  if (_jacSparsity.sparsity.size() > 0) {
463  return;
464  }
465 
469  _jacSparsity.sparsity = jacobianSparsitySet<SparsitySetType, CGBase> (_fun);
470 
471  if (!_custom_jac.defined) {
472  generateSparsityIndexes(_jacSparsity.sparsity, _jacSparsity.rows, _jacSparsity.cols);
473 
474  } else {
475  _jacSparsity.rows = _custom_jac.row;
476  _jacSparsity.cols = _custom_jac.col;
477  }
478 }
479 
480 template<class Base>
481 void ModelCSourceGen<Base>::generateJacobianSparsitySource() {
482  determineJacobianSparsity();
483 
484  generateSparsity2DSource(_name + "_" + FUNCTION_JACOBIAN_SPARSITY, _jacSparsity);
485  _sources[_name + "_" + FUNCTION_JACOBIAN_SPARSITY + ".c"] = _cache.str();
486  _cache.str("");
487 }
488 
489 } // END cg namespace
490 } // END CppAD namespace
491 
492 #endif
CppAD::cg::ModelCSourceGen::generateSparseJacobianSource
virtual void generateSparseJacobianSource(MultiThreadingType multiThreadingType)
Definition: model_c_source_gen_jac.hpp:67
std
STL namespace.
CppAD::cg::LanguageC::setMaxOperationsPerAssignment
void setMaxOperationsPerAssignment(size_t maxOperationsPerAssignment)
Definition: language_c.hpp:273
CppAD::cg::LanguageC::printFunctionDeclaration
static void printFunctionDeclaration(std::ostringstream &out, const std::string &returnType, const std::string &functionName, const std::vector< std::string > &arguments, const std::vector< std::string > &arguments2={})
Definition: language_c.hpp:538
CppAD::cg::ModelCSourceGen::generateSparseJacobianForRevMultiThreadSource
virtual std::string generateSparseJacobianForRevMultiThreadSource(const std::string &functionName, std::map< size_t, CompressedVectorInfo > jacInfo, size_t maxCompressedSize, const std::string &functionRevFor, const std::string &revForSuffix, bool forward, MultiThreadingType multiThreadingType)
Definition: model_c_source_gen_jac.hpp:320
CppAD::cg::CodeHandler::makeVariables
void makeVariables(VectorCG &variables)
Definition: code_handler.hpp:207
CppAD::cg::LanguageC::setMaxAssignmentsPerFunction
virtual void setMaxAssignmentsPerFunction(size_t maxAssignmentsPerFunction, std::map< std::string, std::string > *sources)
Definition: language_c.hpp:252
CppAD::cg::LanguageC::setParameterPrecision
virtual void setParameterPrecision(size_t p)
Definition: language_c.hpp:237
CppAD::cg::CodeHandler::generateCode
virtual void generateCode(std::ostream &out, Language< Base > &lang, CppAD::vector< CGB > &dependent, VariableNameGenerator< Base > &nameGen, const std::string &jobName="source")
Definition: code_handler_impl.hpp:214
CppAD::cg::ModelCSourceGen::generateSparseJacobianForRevSource
virtual void generateSparseJacobianForRevSource(bool forward, MultiThreadingType multiThreadingType)
Definition: model_c_source_gen_jac.hpp:151