2 # ifndef CPPAD_EXAMPLE_MAT_MUL_HPP
3 # define CPPAD_EXAMPLE_MAT_MUL_HPP
76 { assert( i < nr_left );
77 assert( j < n_middle );
78 return (3 + i * n_middle + j) * nk + k;
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;
109 { assert( i < nr_left );
110 assert( j < nc_right );
111 return (i * nc_right + j) * nk + k;
127 size_t nx = 3 + (nr_left + nc_right) * n_middle;
128 size_t nk = tx.
size() / nx;
130 size_t ny = nr_left * nc_right;
131 assert( nk == ty.
size() / ny );
134 size_t k_result = k_left + k_right;
135 assert( k_result < nk );
137 for(
size_t i = 0; i < nr_left; i++)
138 {
for(
size_t j = 0; j < nc_right; j++)
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
144 size_t i_right = right(
145 ell, j, k_right, nk, nr_left, n_middle, nc_right
147 sum += tx[i_left] * tx[i_right];
149 size_t i_result = result(
150 i, j, k_result, nk, nr_left, n_middle, nc_right
172 size_t nx = 3 + (nr_left + nc_right) * n_middle;
173 size_t nk = tx.
size() / nx;
175 size_t ny = nr_left * nc_right;
176 assert( nk == ty.
size() / ny );
181 size_t k_result = k_left + k_right;
182 assert( k_result < nk );
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
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
193 size_t i_right = right(
194 ell, j, k_right, nk, nr_left, n_middle, nc_right
197 px[i_left] += tx[i_right] * py[i_result];
198 px[i_right] += tx[i_left] * py[i_result];
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 ] );
221 size_t nx = 3 + (nr_left + nc_right) * n_middle;
222 size_t ny = nr_left * nc_right;
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() );
234 for(i = 0; i < nr_left; i++)
235 {
for(j = 0; j < nc_right; j++)
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
241 size_t i_right = right(
242 ell, j, k, nk, nr_left, n_middle, nc_right
244 bool nz_left = vx[i_left] |(tx[i_left] != 0.);
245 bool nz_right = vx[i_right]|(tx[i_right] != 0.);
247 if( nz_left & nz_right )
248 var |= bool(vx[i_left]) | bool(vx[i_right]);
250 size_t i_result = result(
251 i, j, k, nk, nr_left, n_middle, nc_right
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
270 for(k = q; k <= p; k++)
272 for(ell = 0; ell <= k; ell++)
274 ell, k - ell, tx, ty, nr_left, n_middle, nc_right
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 ] );
297 size_t nx = 3 + (nr_left + nc_right) * n_middle;
298 size_t ny = nr_left * nc_right;
300 assert( nx * n_order == tx.
size() );
301 assert( ny * n_order == ty.
size() );
306 for(
size_t i = 0; i < px.
size(); i++)
313 for(
size_t ell = 0; ell <= k; ell++)
315 ell, k - ell, tx, ty, px, py, nr_left, n_middle, nc_right
337 size_t nx = 3 + (nr_left + nc_right) * n_middle;
338 size_t ny = nr_left * nc_right;
340 assert( nx == x.
size() );
341 assert( nx * q == r.
size() );
342 assert( ny * q == s.
size() );
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
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
359 size_t i_right = right(
360 ell, j, k, nk, nr_left, n_middle, nc_right
362 for(p = 0; p < q; p++)
364 s[i_result * q + p] |= bool( r[i_left * q + p ] );
365 s[i_result * q + p] |= bool( r[i_right * q + p ] );
375 const vector< std::set<size_t> >& r ,
376 vector< std::set<size_t> >& s ,
383 size_t nx = 3 + (nr_left + nc_right) * n_middle;
384 size_t ny = nr_left * nc_right;
386 assert( nx == x.
size() );
387 assert( nx == r.size() );
388 assert( ny == s.size() );
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
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
403 size_t i_right = right(
404 ell, j, k, nk, nr_left, n_middle, nc_right
407 s[i_result] =
set_union(s[i_result], r[i_left] );
408 s[i_result] =
set_union(s[i_result], r[i_right] );
428 size_t nx = 3 + (nr_left + nc_right) * n_middle;
430 size_t ny = nr_left * nc_right;
432 assert( nx == x.
size() );
433 assert( nx * q == st.
size() );
434 assert( ny * q == rt.
size() );
438 for(i = 0; i < nx; i++)
439 {
for(p = 0; p < q; p++)
440 st[ i * q + p ] =
false;
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
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
455 size_t i_right = right(
456 ell, j, k, nk, nr_left, n_middle, nc_right
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] );
470 const vector< std::set<size_t> >& rt ,
471 vector< std::set<size_t> >& st ,
477 size_t nx = 3 + (nr_left + nc_right) * n_middle;
479 size_t ny = nr_left * nc_right;
481 assert( nx == x.
size() );
482 assert( nx == st.size() );
483 assert( ny == rt.size() );
487 for(i = 0; i < nx; i++)
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
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
502 size_t i_right = right(
503 ell, j, k, nk, nr_left, n_middle, nc_right
506 st[i_left] =
set_union(st[i_left], rt[i_result]);
507 st[i_right] =
set_union(st[i_right], rt[i_result]);
523 const vector< std::set<size_t> >& r ,
524 const vector< std::set<size_t> >& u ,
525 vector< std::set<size_t> >& v ,
531 size_t nx = 3 + (nr_left + nc_right) * n_middle;
533 size_t ny = nr_left * nc_right;
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 );
546 for(j = 0; j < nx; j++)
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
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
561 size_t i_right = right(
562 ell, j, k, nk, nr_left, n_middle, nc_right
567 t[i_left] |= bool( s[i_result] );
568 t[i_right] |= bool( s[i_result] );
576 v[i_left] =
set_union(v[i_left], u[i_result] );
577 v[i_right] =
set_union(v[i_right], u[i_result] );
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] );
604 size_t nx = 3 + (nr_left + nc_right) * n_middle;
606 size_t ny = nr_left * nc_right;
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 );
618 for(j = 0; j < nx; j++)
620 for(p = 0; p < q; p++)
621 v[j * q + p] =
false;
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
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
634 size_t i_right = right(
635 ell, j, k, nk, nr_left, n_middle, nc_right
640 t[i_left] |= bool( s[i_result] );
641 t[i_right] |= bool( s[i_result] );
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] );
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] );
includes the entire CppAD package in the necessary order.
The CppAD Simple Vector template class.
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)
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)
std::set< Element > set_union(const std::set< Element > &left, const std::set< Element > &right)
void clear(void)
free memory and set number of elements to zero
size_t size(void) const
number of elements currently in this vector.
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)
virtual bool rev_sparse_jac(size_t q, const vector< bool > &rt, vector< bool > &st, const vector< double > &x)
virtual bool for_sparse_jac(size_t q, const vector< bool > &r, vector< bool > &s, const vector< double > &x)
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)
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)
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)
virtual bool reverse(size_t p, const vector< double > &tx, const vector< double > &ty, vector< double > &px, const vector< double > &py)
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)
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)
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)
virtual bool forward(size_t q, size_t p, const vector< bool > &vx, vector< bool > &vy, const vector< double > &tx, vector< double > &ty)