CppAD: A C++ Algorithmic Differentiation Package  20171217
mat_mul.hpp
Go to the documentation of this file.
1 // \$Id\$
4
5 /* --------------------------------------------------------------------------
7
9 the terms of the
10  Eclipse Public License Version 1.0.
11
12 A copy of this license is included in the COPYING file of this distribution.
14 -------------------------------------------------------------------------- */
15
16 /*
17 \$begin atomic_mat_mul.hpp\$\$
18 \$spell
19  Taylor
20  ty
21  px
23  jac
24  hes
25  nr
26  nc
27 \$\$
28
29 \$section Matrix Multiply as an Atomic Operation\$\$
30
32 \$cref atomic_eigen_mat_mul.hpp\$\$
33
35 This example puts the matrix dimensions in the atomic function arguments,
36 instead of the \$cref/constructor/atomic_ctor/\$\$, so that they can
37 be different for different calls to the atomic function.
38 These dimensions are:
39 \$table
40 \$icode nr_left\$\$ \$cnext number of rows in the left matrix \$rend
41 \$icode n_middle\$\$ \$cnext rows in the left matrix and columns in right \$rend
42 \$icode nc_right\$\$ \$cnext number of columns in the right matrix
43 \$tend
44
46 \$srccode%cpp% */
48 namespace { // Begin empty namespace
50 //
52 //
53 // matrix result = left * right
54 class atomic_mat_mul : public CppAD::atomic_base<double> {
55 /* %\$\$
57 \$srccode%cpp% */
58 public:
59  // ---------------------------------------------------------------------
60  // constructor
62  { }
63 private:
64 /* %\$\$
65 \$head Left Operand Element Index\$\$
66 Index in the Taylor coefficient matrix \$icode tx\$\$ of a left matrix element.
67 \$srccode%cpp% */
68  size_t left(
69  size_t i , // left matrix row index
70  size_t j , // left matrix column index
71  size_t k , // Taylor coeffocient order
72  size_t nk , // number of Taylor coefficients in tx
73  size_t nr_left , // rows in left matrix
74  size_t n_middle , // rows in left and columns in right
75  size_t nc_right ) // columns in right matrix
76  { assert( i < nr_left );
77  assert( j < n_middle );
78  return (3 + i * n_middle + j) * nk + k;
79  }
80 /* %\$\$
81 \$head Right Operand Element Index\$\$
82 Index in the Taylor coefficient matrix \$icode tx\$\$ of a right matrix element.
83 \$srccode%cpp% */
84  size_t right(
85  size_t i , // right matrix row index
86  size_t j , // right matrix column index
87  size_t k , // Taylor coeffocient order
88  size_t nk , // number of Taylor coefficients in tx
89  size_t nr_left , // rows in left matrix
90  size_t n_middle , // rows in left and columns in right
91  size_t nc_right ) // columns in right matrix
92  { assert( i < n_middle );
93  assert( j < nc_right );
94  size_t offset = 3 + nr_left * n_middle;
95  return (offset + i * nc_right + j) * nk + k;
96  }
97 /* %\$\$
99 Index in the Taylor coefficient matrix \$icode ty\$\$ of a result matrix element.
100 \$srccode%cpp% */
101  size_t result(
102  size_t i , // result matrix row index
103  size_t j , // result matrix column index
104  size_t k , // Taylor coeffocient order
105  size_t nk , // number of Taylor coefficients in ty
106  size_t nr_left , // rows in left matrix
107  size_t n_middle , // rows in left and columns in right
108  size_t nc_right ) // columns in right matrix
109  { assert( i < nr_left );
110  assert( j < nc_right );
111  return (i * nc_right + j) * nk + k;
112  }
113 /* %\$\$
115 Forward mode multiply Taylor coefficients in \$icode tx\$\$ and sum into
116 \$icode ty\$\$ (for one pair of left and right orders)
117 \$srccode%cpp% */
119  size_t k_left , // order for left coefficients
120  size_t k_right , // order for right coefficients
121  const vector<double>& tx , // domain space Taylor coefficients
122  vector<double>& ty , // range space Taylor coefficients
123  size_t nr_left , // rows in left matrix
124  size_t n_middle , // rows in left and columns in right
125  size_t nc_right ) // columns in right matrix
126  {
127  size_t nx = 3 + (nr_left + nc_right) * n_middle;
128  size_t nk = tx.size() / nx;
129 # ifndef NDEBUG
130  size_t ny = nr_left * nc_right;
131  assert( nk == ty.size() / ny );
132 # endif
133  //
134  size_t k_result = k_left + k_right;
135  assert( k_result < nk );
136  //
137  for(size_t i = 0; i < nr_left; i++)
138  { for(size_t j = 0; j < nc_right; j++)
139  { double sum = 0.0;
140  for(size_t ell = 0; ell < n_middle; ell++)
141  { size_t i_left = left(
142  i, ell, k_left, nk, nr_left, n_middle, nc_right
143  );
144  size_t i_right = right(
145  ell, j, k_right, nk, nr_left, n_middle, nc_right
146  );
147  sum += tx[i_left] * tx[i_right];
148  }
149  size_t i_result = result(
150  i, j, k_result, nk, nr_left, n_middle, nc_right
151  );
152  ty[i_result] += sum;
153  }
154  }
155  }
156 /* %\$\$
158 Reverse mode partials of Taylor coefficients and sum into \$icode px\$\$
159 (for one pair of left and right orders)
160 \$srccode%cpp% */
162  size_t k_left , // order for left coefficients
163  size_t k_right , // order for right coefficients
164  const vector<double>& tx , // domain space Taylor coefficients
165  const vector<double>& ty , // range space Taylor coefficients
166  vector<double>& px , // partials w.r.t. tx
167  const vector<double>& py , // partials w.r.t. ty
168  size_t nr_left , // rows in left matrix
169  size_t n_middle , // rows in left and columns in right
170  size_t nc_right ) // columns in right matrix
171  {
172  size_t nx = 3 + (nr_left + nc_right) * n_middle;
173  size_t nk = tx.size() / nx;
174 # ifndef NDEBUG
175  size_t ny = nr_left * nc_right;
176  assert( nk == ty.size() / ny );
177 # endif
178  assert( tx.size() == px.size() );
179  assert( ty.size() == py.size() );
180  //
181  size_t k_result = k_left + k_right;
182  assert( k_result < nk );
183  //
184  for(size_t i = 0; i < nr_left; i++)
185  { for(size_t j = 0; j < nc_right; j++)
186  { size_t i_result = result(
187  i, j, k_result, nk, nr_left, n_middle, nc_right
188  );
189  for(size_t ell = 0; ell < n_middle; ell++)
190  { size_t i_left = left(
191  i, ell, k_left, nk, nr_left, n_middle, nc_right
192  );
193  size_t i_right = right(
194  ell, j, k_right, nk, nr_left, n_middle, nc_right
195  );
196  // sum += tx[i_left] * tx[i_right];
197  px[i_left] += tx[i_right] * py[i_result];
198  px[i_right] += tx[i_left] * py[i_result];
199  }
200  }
201  }
202  return;
203  }
204 /* %\$\$
206 Routine called by CppAD during \$cref Forward\$\$ mode.
207 \$srccode%cpp% */
208  virtual bool forward(
209  size_t q ,
210  size_t p ,
211  const vector<bool>& vx ,
212  vector<bool>& vy ,
213  const vector<double>& tx ,
214  vector<double>& ty
215  )
216  { size_t n_order = p + 1;
217  size_t nr_left = size_t( tx[ 0 * n_order + 0 ] );
218  size_t n_middle = size_t( tx[ 1 * n_order + 0 ] );
219  size_t nc_right = size_t( tx[ 2 * n_order + 0 ] );
220 # ifndef NDEBUG
221  size_t nx = 3 + (nr_left + nc_right) * n_middle;
222  size_t ny = nr_left * nc_right;
223 # endif
224  assert( vx.size() == 0 || nx == vx.size() );
225  assert( vx.size() == 0 || ny == vy.size() );
226  assert( nx * n_order == tx.size() );
227  assert( ny * n_order == ty.size() );
228  size_t i, j, ell;
229
230  // check if we are computing vy information
231  if( vx.size() > 0 )
232  { size_t nk = 1;
233  size_t k = 0;
234  for(i = 0; i < nr_left; i++)
235  { for(j = 0; j < nc_right; j++)
236  { bool var = false;
237  for(ell = 0; ell < n_middle; ell++)
238  { size_t i_left = left(
239  i, ell, k, nk, nr_left, n_middle, nc_right
240  );
241  size_t i_right = right(
242  ell, j, k, nk, nr_left, n_middle, nc_right
243  );
244  bool nz_left = vx[i_left] |(tx[i_left] != 0.);
245  bool nz_right = vx[i_right]|(tx[i_right] != 0.);
246  // if not multiplying by the constant zero
247  if( nz_left & nz_right )
248  var |= bool(vx[i_left]) | bool(vx[i_right]);
249  }
250  size_t i_result = result(
251  i, j, k, nk, nr_left, n_middle, nc_right
252  );
253  vy[i_result] = var;
254  }
255  }
256  }
257
258  // initialize result as zero
259  size_t k;
260  for(i = 0; i < nr_left; i++)
261  { for(j = 0; j < nc_right; j++)
262  { for(k = q; k <= p; k++)
263  { size_t i_result = result(
264  i, j, k, n_order, nr_left, n_middle, nc_right
265  );
266  ty[i_result] = 0.0;
267  }
268  }
269  }
270  for(k = q; k <= p; k++)
271  { // sum the produces that result in order k
272  for(ell = 0; ell <= k; ell++)
273  forward_multiply(
274  ell, k - ell, tx, ty, nr_left, n_middle, nc_right
275  );
276  }
277
278  // all orders are implented, so always return true
279  return true;
280  }
281 /* %\$\$
283 Routine called by CppAD during \$cref Reverse\$\$ mode.
284 \$srccode%cpp% */
285  virtual bool reverse(
286  size_t p ,
287  const vector<double>& tx ,
288  const vector<double>& ty ,
289  vector<double>& px ,
290  const vector<double>& py
291  )
292  { size_t n_order = p + 1;
293  size_t nr_left = size_t( tx[ 0 * n_order + 0 ] );
294  size_t n_middle = size_t( tx[ 1 * n_order + 0 ] );
295  size_t nc_right = size_t( tx[ 2 * n_order + 0 ] );
296 # ifndef NDEBUG
297  size_t nx = 3 + (nr_left + nc_right) * n_middle;
298  size_t ny = nr_left * nc_right;
299 # endif
300  assert( nx * n_order == tx.size() );
301  assert( ny * n_order == ty.size() );
302  assert( px.size() == tx.size() );
303  assert( py.size() == ty.size() );
304
305  // initialize summation
306  for(size_t i = 0; i < px.size(); i++)
307  px[i] = 0.0;
308
309  // number of orders to differentiate
310  size_t k = n_order;
311  while(k--)
312  { // differentiate the produces that result in order k
313  for(size_t ell = 0; ell <= k; ell++)
314  reverse_multiply(
315  ell, k - ell, tx, ty, px, py, nr_left, n_middle, nc_right
316  );
317  }
318
319  // all orders are implented, so always return true
320  return true;
321  }
322 /* %\$\$
324 Routines called by CppAD during \$cref ForSparseJac\$\$.
325 \$srccode%cpp% */
326  // boolean sparsity patterns
327  virtual bool for_sparse_jac(
328  size_t q ,
329  const vector<bool>& r ,
330  vector<bool>& s ,
331  const vector<double>& x )
332  {
333  size_t nr_left = size_t( CppAD::Integer( x[0] ) );
334  size_t n_middle = size_t( CppAD::Integer( x[1] ) );
335  size_t nc_right = size_t( CppAD::Integer( x[2] ) );
336 # ifndef NDEBUG
337  size_t nx = 3 + (nr_left + nc_right) * n_middle;
338  size_t ny = nr_left * nc_right;
339 # endif
340  assert( nx == x.size() );
341  assert( nx * q == r.size() );
342  assert( ny * q == s.size() );
343  size_t p;
344
345  // sparsity for S(x) = f'(x) * R
346  size_t nk = 1;
347  size_t k = 0;
348  for(size_t i = 0; i < nr_left; i++)
349  { for(size_t j = 0; j < nc_right; j++)
350  { size_t i_result = result(
351  i, j, k, nk, nr_left, n_middle, nc_right
352  );
353  for(p = 0; p < q; p++)
354  s[i_result * q + p] = false;
355  for(size_t ell = 0; ell < n_middle; ell++)
356  { size_t i_left = left(
357  i, ell, k, nk, nr_left, n_middle, nc_right
358  );
359  size_t i_right = right(
360  ell, j, k, nk, nr_left, n_middle, nc_right
361  );
362  for(p = 0; p < q; p++)
363  { // cast avoids Microsoft warning (should not be needed)
364  s[i_result * q + p] |= bool( r[i_left * q + p ] );
365  s[i_result * q + p] |= bool( r[i_right * q + p ] );
366  }
367  }
368  }
369  }
370  return true;
371  }
372  // set sparsity patterns
373  virtual bool for_sparse_jac(
374  size_t q ,
375  const vector< std::set<size_t> >& r ,
376  vector< std::set<size_t> >& s ,
377  const vector<double>& x )
378  {
379  size_t nr_left = size_t( CppAD::Integer( x[0] ) );
380  size_t n_middle = size_t( CppAD::Integer( x[1] ) );
381  size_t nc_right = size_t( CppAD::Integer( x[2] ) );
382 # ifndef NDEBUG
383  size_t nx = 3 + (nr_left + nc_right) * n_middle;
384  size_t ny = nr_left * nc_right;
385 # endif
386  assert( nx == x.size() );
387  assert( nx == r.size() );
388  assert( ny == s.size() );
389
390  // sparsity for S(x) = f'(x) * R
391  size_t nk = 1;
392  size_t k = 0;
393  for(size_t i = 0; i < nr_left; i++)
394  { for(size_t j = 0; j < nc_right; j++)
395  { size_t i_result = result(
396  i, j, k, nk, nr_left, n_middle, nc_right
397  );
398  s[i_result].clear();
399  for(size_t ell = 0; ell < n_middle; ell++)
400  { size_t i_left = left(
401  i, ell, k, nk, nr_left, n_middle, nc_right
402  );
403  size_t i_right = right(
404  ell, j, k, nk, nr_left, n_middle, nc_right
405  );
406  //
407  s[i_result] = set_union(s[i_result], r[i_left] );
408  s[i_result] = set_union(s[i_result], r[i_right] );
409  }
410  }
411  }
412  return true;
413  }
414 /* %\$\$
416 Routines called by CppAD during \$cref RevSparseJac\$\$.
417 \$srccode%cpp% */
418  // boolean sparsity patterns
419  virtual bool rev_sparse_jac(
420  size_t q ,
421  const vector<bool>& rt ,
422  vector<bool>& st ,
423  const vector<double>& x )
424  {
425  size_t nr_left = size_t( CppAD::Integer( x[0] ) );
426  size_t n_middle = size_t( CppAD::Integer( x[1] ) );
427  size_t nc_right = size_t( CppAD::Integer( x[2] ) );
428  size_t nx = 3 + (nr_left + nc_right) * n_middle;
429 # ifndef NDEBUG
430  size_t ny = nr_left * nc_right;
431 # endif
432  assert( nx == x.size() );
433  assert( nx * q == st.size() );
434  assert( ny * q == rt.size() );
435  size_t i, j, p;
436
437  // initialize
438  for(i = 0; i < nx; i++)
439  { for(p = 0; p < q; p++)
440  st[ i * q + p ] = false;
441  }
442
443  // sparsity for S(x)^T = f'(x)^T * R^T
444  size_t nk = 1;
445  size_t k = 0;
446  for(i = 0; i < nr_left; i++)
447  { for(j = 0; j < nc_right; j++)
448  { size_t i_result = result(
449  i, j, k, nk, nr_left, n_middle, nc_right
450  );
451  for(size_t ell = 0; ell < n_middle; ell++)
452  { size_t i_left = left(
453  i, ell, k, nk, nr_left, n_middle, nc_right
454  );
455  size_t i_right = right(
456  ell, j, k, nk, nr_left, n_middle, nc_right
457  );
458  for(p = 0; p < q; p++)
459  { st[i_left * q + p] |= bool( rt[i_result * q + p] );
460  st[i_right* q + p] |= bool( rt[i_result * q + p] );
461  }
462  }
463  }
464  }
465  return true;
466  }
467  // set sparsity patterns
468  virtual bool rev_sparse_jac(
469  size_t q ,
470  const vector< std::set<size_t> >& rt ,
471  vector< std::set<size_t> >& st ,
472  const vector<double>& x )
473  {
474  size_t nr_left = size_t( CppAD::Integer( x[0] ) );
475  size_t n_middle = size_t( CppAD::Integer( x[1] ) );
476  size_t nc_right = size_t( CppAD::Integer( x[2] ) );
477  size_t nx = 3 + (nr_left + nc_right) * n_middle;
478 # ifndef NDEBUG
479  size_t ny = nr_left * nc_right;
480 # endif
481  assert( nx == x.size() );
482  assert( nx == st.size() );
483  assert( ny == rt.size() );
484  size_t i, j;
485
486  // initialize
487  for(i = 0; i < nx; i++)
488  st[i].clear();
489
490  // sparsity for S(x)^T = f'(x)^T * R^T
491  size_t nk = 1;
492  size_t k = 0;
493  for(i = 0; i < nr_left; i++)
494  { for(j = 0; j < nc_right; j++)
495  { size_t i_result = result(
496  i, j, k, nk, nr_left, n_middle, nc_right
497  );
498  for(size_t ell = 0; ell < n_middle; ell++)
499  { size_t i_left = left(
500  i, ell, k, nk, nr_left, n_middle, nc_right
501  );
502  size_t i_right = right(
503  ell, j, k, nk, nr_left, n_middle, nc_right
504  );
505  //
506  st[i_left] = set_union(st[i_left], rt[i_result]);
507  st[i_right] = set_union(st[i_right], rt[i_result]);
508  }
509  }
510  }
511  return true;
512  }
513 /* %\$\$
515 Routines called by \$cref RevSparseHes\$\$.
516 \$srccode%cpp% */
517  // set sparsity patterns
518  virtual bool rev_sparse_hes(
519  const vector<bool>& vx,
520  const vector<bool>& s ,
521  vector<bool>& t ,
522  size_t q ,
523  const vector< std::set<size_t> >& r ,
524  const vector< std::set<size_t> >& u ,
525  vector< std::set<size_t> >& v ,
526  const vector<double>& x )
527  {
528  size_t nr_left = size_t( CppAD::Integer( x[0] ) );
529  size_t n_middle = size_t( CppAD::Integer( x[1] ) );
530  size_t nc_right = size_t( CppAD::Integer( x[2] ) );
531  size_t nx = 3 + (nr_left + nc_right) * n_middle;
532 # ifndef NDEBUG
533  size_t ny = nr_left * nc_right;
534 # endif
535  assert( x.size() == nx );
536  assert( vx.size() == nx );
537  assert( t.size() == nx );
538  assert( r.size() == nx );
539  assert( v.size() == nx );
540  assert( s.size() == ny );
541  assert( u.size() == ny );
542  //
543  size_t i, j;
544  //
545  // initilaize sparsity patterns as false
546  for(j = 0; j < nx; j++)
547  { t[j] = false;
548  v[j].clear();
549  }
550  size_t nk = 1;
551  size_t k = 0;
552  for(i = 0; i < nr_left; i++)
553  { for(j = 0; j < nc_right; j++)
554  { size_t i_result = result(
555  i, j, k, nk, nr_left, n_middle, nc_right
556  );
557  for(size_t ell = 0; ell < n_middle; ell++)
558  { size_t i_left = left(
559  i, ell, k, nk, nr_left, n_middle, nc_right
560  );
561  size_t i_right = right(
562  ell, j, k, nk, nr_left, n_middle, nc_right
563  );
564  //
565  // Compute sparsity for T(x) = S(x) * f'(x).
566  // We need not use vx with f'(x) back propagation.
567  t[i_left] |= bool( s[i_result] );
568  t[i_right] |= bool( s[i_result] );
569
570  // V(x) = f'(x)^T * U(x) + S(x) * f''(x) * R
571  // U(x) = g''(y) * f'(x) * R
572  // S(x) = g'(y)
573
574  // back propagate f'(x)^T * U(x)
575  // (no need to use vx with f'(x) propogation)
576  v[i_left] = set_union(v[i_left], u[i_result] );
577  v[i_right] = set_union(v[i_right], u[i_result] );
578
579  // back propagate S(x) * f''(x) * R
580  // (here is where we must check for cross terms)
581  if( s[i_result] & vx[i_left] & vx[i_right] )
582  { v[i_left] = set_union(v[i_left], r[i_right] );
583  v[i_right] = set_union(v[i_right], r[i_left] );
584  }
585  }
586  }
587  }
588  return true;
589  }
590  // bool sparsity
591  virtual bool rev_sparse_hes(
592  const vector<bool>& vx,
593  const vector<bool>& s ,
594  vector<bool>& t ,
595  size_t q ,
596  const vector<bool>& r ,
597  const vector<bool>& u ,
598  vector<bool>& v ,
599  const vector<double>& x )
600  {
601  size_t nr_left = size_t( CppAD::Integer( x[0] ) );
602  size_t n_middle = size_t( CppAD::Integer( x[1] ) );
603  size_t nc_right = size_t( CppAD::Integer( x[2] ) );
604  size_t nx = 3 + (nr_left + nc_right) * n_middle;
605 # ifndef NDEBUG
606  size_t ny = nr_left * nc_right;
607 # endif
608  assert( x.size() == nx );
609  assert( vx.size() == nx );
610  assert( t.size() == nx );
611  assert( r.size() == nx * q );
612  assert( v.size() == nx * q );
613  assert( s.size() == ny );
614  assert( u.size() == ny * q );
615  size_t i, j, p;
616  //
617  // initilaize sparsity patterns as false
618  for(j = 0; j < nx; j++)
619  { t[j] = false;
620  for(p = 0; p < q; p++)
621  v[j * q + p] = false;
622  }
623  size_t nk = 1;
624  size_t k = 0;
625  for(i = 0; i < nr_left; i++)
626  { for(j = 0; j < nc_right; j++)
627  { size_t i_result = result(
628  i, j, k, nk, nr_left, n_middle, nc_right
629  );
630  for(size_t ell = 0; ell < n_middle; ell++)
631  { size_t i_left = left(
632  i, ell, k, nk, nr_left, n_middle, nc_right
633  );
634  size_t i_right = right(
635  ell, j, k, nk, nr_left, n_middle, nc_right
636  );
637  //
638  // Compute sparsity for T(x) = S(x) * f'(x).
639  // We so not need to use vx with f'(x) propagation.
640  t[i_left] |= bool( s[i_result] );
641  t[i_right] |= bool( s[i_result] );
642
643  // V(x) = f'(x)^T * U(x) + S(x) * f''(x) * R
644  // U(x) = g''(y) * f'(x) * R
645  // S(x) = g'(y)
646
647  // back propagate f'(x)^T * U(x)
648  // (no need to use vx with f'(x) propogation)
649  for(p = 0; p < q; p++)
650  { v[ i_left * q + p] |= bool( u[ i_result * q + p] );
651  v[ i_right * q + p] |= bool( u[ i_result * q + p] );
652  }
653
654  // back propagate S(x) * f''(x) * R
655  // (here is where we must check for cross terms)
656  if( s[i_result] & vx[i_left] & vx[i_right] )
657  { for(p = 0; p < q; p++)
658  { v[i_left * q + p] |= bool( r[i_right * q + p] );
659  v[i_right * q + p] |= bool( r[i_left * q + p] );
660  }
661  }
662  }
663  }
664  }
665  return true;
666  }
667 /* %\$\$
669 \$srccode%cpp% */
670 }; // End of mat_mul class
671 } // End empty namespace
672 /* %\$\$
673 \$comment end nospell\$\$
674 \$end
675 */
676
677
678 # endif
includes the entire CppAD package in the necessary order.
The CppAD Simple Vector template class.
Definition: vector.hpp:338
virtual bool for_sparse_jac(size_t q, const vector< std::set< size_t > > &r, vector< std::set< size_t > > &s, const vector< double > &x)
Definition: mat_mul.hpp:373
size_t right(size_t i, size_t j, size_t k, size_t nk, size_t nr_left, size_t n_middle, size_t nc_right)
Definition: mat_mul.hpp:84
std::set< Element > set_union(const std::set< Element > &left, const std::set< Element > &right)
Definition: set_union.hpp:73
void clear(void)
free memory and set number of elements to zero
Definition: vector.hpp:418
size_t size(void) const
number of elements currently in this vector.
Definition: vector.hpp:387
virtual bool rev_sparse_hes(const vector< bool > &vx, const vector< bool > &s, vector< bool > &t, size_t q, const vector< std::set< size_t > > &r, const vector< std::set< size_t > > &u, vector< std::set< size_t > > &v, const vector< double > &x)
Definition: mat_mul.hpp:518
virtual bool rev_sparse_jac(size_t q, const vector< bool > &rt, vector< bool > &st, const vector< double > &x)
Definition: mat_mul.hpp:419
virtual bool for_sparse_jac(size_t q, const vector< bool > &r, vector< bool > &s, const vector< double > &x)
Definition: mat_mul.hpp:327
size_t left(size_t i, size_t j, size_t k, size_t nk, size_t nr_left, size_t n_middle, size_t nc_right)
Definition: mat_mul.hpp:68
void reverse_multiply(size_t k_left, size_t k_right, const vector< double > &tx, const vector< double > &ty, vector< double > &px, const vector< double > &py, size_t nr_left, size_t n_middle, size_t nc_right)
Definition: mat_mul.hpp:161
virtual bool rev_sparse_jac(size_t q, const vector< std::set< size_t > > &rt, vector< std::set< size_t > > &st, const vector< double > &x)
Definition: mat_mul.hpp:468
virtual bool reverse(size_t p, const vector< double > &tx, const vector< double > &ty, vector< double > &px, const vector< double > &py)
Definition: mat_mul.hpp:285
virtual bool rev_sparse_hes(const vector< bool > &vx, const vector< bool > &s, vector< bool > &t, size_t q, const vector< bool > &r, const vector< bool > &u, vector< bool > &v, const vector< double > &x)
Definition: mat_mul.hpp:591
int Integer(const std::complex< double > &x)
size_t result(size_t i, size_t j, size_t k, size_t nk, size_t nr_left, size_t n_middle, size_t nc_right)
Definition: mat_mul.hpp:101
void forward_multiply(size_t k_left, size_t k_right, const vector< double > &tx, vector< double > &ty, size_t nr_left, size_t n_middle, size_t nc_right)
Definition: mat_mul.hpp:118
virtual bool forward(size_t q, size_t p, const vector< bool > &vx, vector< bool > &vy, const vector< double > &tx, vector< double > &ty)
Definition: mat_mul.hpp:208