1 # ifndef CPPAD_EXAMPLE_EIGEN_MAT_MUL_HPP
2 # define CPPAD_EXAMPLE_EIGEN_MAT_MUL_HPP
110 # include <Eigen/Core>
120 template <
class Base>
129 typedef Eigen::Matrix<
132 typedef Eigen::Matrix<
139 "atom_eigen_mat_mul" ,
140 CppAD::atomic_base<Base>::set_sparsity_enum
150 {
size_t nr_left = size_t( left.rows() );
151 size_t n_middle = size_t( left.cols() );
152 size_t nc_right = size_t( right.cols() );
153 assert( n_middle ==
size_t( right.rows() ) );
154 size_t nx = 3 + (nr_left + nc_right) * n_middle;
155 size_t ny = nr_left * nc_right;
156 size_t n_left = nr_left * n_middle;
157 size_t n_right = n_middle * nc_right;
158 size_t n_result = nr_left * nc_right;
160 assert( 3 + n_left + n_right == nx );
161 assert( n_result == ny );
164 CPPAD_TESTVECTOR(
ad_scalar) packed_arg(nx);
169 for(
size_t i = 0; i < n_left; i++)
170 packed_arg[3 + i] = left.data()[i];
171 for(
size_t i = 0; i < n_right; i++)
172 packed_arg[ 3 + n_left + i ] = right.data()[i];
177 CPPAD_TESTVECTOR(
ad_scalar) packed_result(ny);
178 (*this)(packed_arg, packed_result);
182 for(
size_t i = 0; i < n_result; i++)
183 result.data()[i] = packed_result[ i ];
217 {
size_t n_order = q + 1;
218 size_t nr_left = size_t(
CppAD::Integer( tx[ 0 * n_order + 0 ] ) );
219 size_t n_middle = size_t(
CppAD::Integer( tx[ 1 * n_order + 0 ] ) );
220 size_t nc_right = size_t(
CppAD::Integer( tx[ 2 * n_order + 0 ] ) );
222 size_t nx = 3 + (nr_left + nc_right) * n_middle;
223 size_t ny = nr_left * nc_right;
226 assert( vx.
size() == 0 || nx == vx.
size() );
227 assert( vx.
size() == 0 || ny == vy.
size() );
228 assert( nx * n_order == tx.
size() );
229 assert( ny * n_order == ty.
size() );
231 size_t n_left = nr_left * n_middle;
232 size_t n_right = n_middle * nc_right;
233 size_t n_result = nr_left * nc_right;
234 assert( 3 + n_left + n_right == nx );
235 assert( n_result == ny );
239 assert( f_left_.size() == f_right_.size() );
240 assert( f_left_.size() == f_result_.size() );
241 if( f_left_.size() < n_order )
242 { f_left_.resize(n_order);
243 f_right_.resize(n_order);
244 f_result_.resize(n_order);
246 for(
size_t k = 0; k < n_order; k++)
247 { f_left_[k].resize(nr_left, n_middle);
248 f_right_[k].resize(n_middle, nc_right);
249 f_result_[k].resize(nr_left, nc_right);
254 for(
size_t k = 0; k < n_order; k++)
256 for(
size_t i = 0; i < n_left; i++)
257 f_left_[k].data()[i] = tx[ (3 + i) * n_order + k ];
260 for(
size_t i = 0; i < n_right; i++)
261 f_right_[k].data()[i] = tx[ ( 3 + n_left + i) * n_order + k ];
266 for(
size_t k = 0; k < n_order; k++)
268 f_result_[k] = matrix::Zero(nr_left, nc_right);
269 for(
size_t ell = 0; ell <= k; ell++)
270 f_result_[k] += f_left_[ell] * f_right_[k-ell];
274 for(
size_t k = 0; k < n_order; k++)
275 {
for(
size_t i = 0; i < n_result; i++)
276 ty[ i * n_order + k ] = f_result_[k].data()[i];
286 assert( n_order == 1 );
287 for(
size_t i = 0; i < nr_left; i++)
288 {
for(
size_t j = 0; j < nc_right; j++)
290 for(
size_t ell = 0; ell < n_middle; ell++)
292 size_t index = 3 + i * n_middle + ell;
293 bool var_left = vx[index];
294 bool nz_left = var_left | (f_left_[0](i, ell) != zero);
296 index = 3 + n_left + ell * nc_right + j;
297 bool var_right = vx[index];
298 bool nz_right = var_right | (f_right_[0](ell, j) != zero);
300 var |= var_left & nz_right;
301 var |= nz_left & var_right;
303 size_t index = i * nc_right + j;
325 {
size_t n_order = q + 1;
326 size_t nr_left = size_t(
CppAD::Integer( tx[ 0 * n_order + 0 ] ) );
327 size_t n_middle = size_t(
CppAD::Integer( tx[ 1 * n_order + 0 ] ) );
328 size_t nc_right = size_t(
CppAD::Integer( tx[ 2 * n_order + 0 ] ) );
330 size_t nx = 3 + (nr_left + nc_right) * n_middle;
331 size_t ny = nr_left * nc_right;
334 assert( nx * n_order == tx.
size() );
335 assert( ny * n_order == ty.
size() );
339 size_t n_left = nr_left * n_middle;
340 size_t n_right = n_middle * nc_right;
341 size_t n_result = nr_left * nc_right;
342 assert( 3 + n_left + n_right == nx );
343 assert( n_result == ny );
346 assert( f_left_.size() == f_right_.size() );
347 assert( f_left_.size() == f_result_.size() );
349 assert( f_left_.size() >= n_order );
352 assert( r_left_.size() == r_right_.size() );
353 assert( r_left_.size() == r_result_.size() );
354 if( r_left_.size() < n_order )
355 { r_left_.resize(n_order);
356 r_right_.resize(n_order);
357 r_result_.resize(n_order);
359 for(
size_t k = 0; k < n_order; k++)
360 { r_left_[k].resize(nr_left, n_middle);
361 r_right_[k].resize(n_middle, nc_right);
362 r_result_[k].resize(nr_left, nc_right);
367 for(
size_t k = 0; k < n_order; k++)
369 for(
size_t i = 0; i < n_left; i++)
370 f_left_[k].data()[i] = tx[ (3 + i) * n_order + k ];
373 for(
size_t i = 0; i < n_right; i++)
374 f_right_[k].data()[i] = tx[ (3 + n_left + i) * n_order + k ];
378 for(
size_t k = 0; k < n_order; k++)
379 {
for(
size_t i = 0; i < n_result; i++)
380 r_result_[k].data()[i] = py[ i * n_order + k ];
384 for(
size_t k = 0; k < n_order; k++)
385 { r_left_[k] = matrix::Zero(nr_left, n_middle);
386 r_right_[k] = matrix::Zero(n_middle, nc_right);
390 for(
size_t k1 = n_order; k1 > 0; k1--)
392 for(
size_t ell = 0; ell <= k; ell++)
394 r_left_[ell] += r_result_[k] * f_right_[k-ell].transpose();
396 r_right_[k-ell] += f_left_[ell].transpose() * r_result_[k];
401 for(
size_t k = 0; k < n_order; k++)
403 px[ 0 * n_order + k ] = 0.0;
404 px[ 1 * n_order + k ] = 0.0;
405 px[ 2 * n_order + k ] = 0.0;
408 for(
size_t i = 0; i < n_left; i++)
409 px[ (3 + i) * n_order + k ] = r_left_[k].
data()[i];
412 for(
size_t i = 0; i < n_right; i++)
413 px[ (3 + i + n_left) * n_order + k] = r_right_[k].
data()[i];
435 size_t nx = 3 + (nr_left + nc_right) * n_middle;
436 size_t ny = nr_left * nc_right;
439 assert( nx == r.size() );
440 assert( ny == s.size() );
442 size_t n_left = nr_left * n_middle;
443 for(
size_t i = 0; i < nr_left; i++)
444 {
for(
size_t j = 0; j < nc_right; j++)
446 size_t i_result = i * nc_right + j;
448 for(
size_t ell = 0; ell < n_middle; ell++)
450 size_t i_left = 3 + i * n_middle + ell;
452 size_t i_right = 3 + n_left + ell * nc_right + j;
455 bool zero = x[i_left] == Base(0.0) || x[i_right] == Base(0);
483 size_t nx = 3 + (nr_left + nc_right) * n_middle;
485 size_t ny = nr_left * nc_right;
488 assert( nx == st.size() );
489 assert( ny == rt.size() );
492 for(
size_t i = 0; i < nx; i++)
496 size_t n_left = nr_left * n_middle;
497 for(
size_t i = 0; i < nr_left; i++)
498 {
for(
size_t j = 0; j < nc_right; j++)
500 size_t i_result = i * nc_right + j;
501 st[i_result].clear();
502 for(
size_t ell = 0; ell < n_middle; ell++)
504 size_t i_left = 3 + i * n_middle + ell;
506 size_t i_right = 3 + n_left + ell * nc_right + j;
532 size_t nx = 3 + (nr_left + nc_right) * n_middle;
534 size_t ny = nr_left * nc_right;
537 assert( vx.
size() == nx );
538 assert( r.
size() == nx );
539 assert( s.
size() == ny );
540 assert( h.size() == nx );
543 for(
size_t i = 0; i < nx; i++)
546 size_t n_left = nr_left * n_middle;
547 for(
size_t i = 0; i < nr_left; i++)
548 {
for(
size_t j = 0; j < nc_right; j++)
550 size_t i_result = i * nc_right + j;
552 {
for(
size_t ell = 0; ell < n_middle; ell++)
554 size_t i_left = 3 + i * n_middle + ell;
556 size_t i_right = 3 + n_left + ell * nc_right + j;
557 if( r[i_left] & r[i_right] )
558 { h[i_left].insert(i_right);
559 h[i_right].insert(i_left);
593 size_t nx = 3 + (nr_left + nc_right) * n_middle;
595 size_t ny = nr_left * nc_right;
598 assert( vx.
size() == nx );
599 assert( s.
size() == ny );
600 assert( t.
size() == nx );
601 assert( r.size() == nx );
602 assert( v.size() == nx );
605 for(
size_t j = 0; j < nx; j++)
610 size_t n_left = nr_left * n_middle;
611 for(
size_t i = 0; i < nr_left; i++)
612 {
for(
size_t j = 0; j < nc_right; j++)
614 size_t i_result = i * nc_right + j;
615 for(
size_t ell = 0; ell < n_middle; ell++)
617 size_t i_left = 3 + i * n_middle + ell;
619 size_t i_right = 3 + n_left + ell * nc_right + j;
622 t[i_left] |= bool( s[i_result] );
623 t[i_right] |= bool( s[i_result] );
635 if( s[i_result] & vx[i_left] & vx[i_right] )
includes the entire CppAD package in the necessary order.
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
Type * data(void)
raw pointer to the data
atomic_eigen_mat_mul(void)
ad_matrix op(const ad_matrix &left, const ad_matrix &right)
size_t size(void) const
number of elements currently in this vector.
CppAD::vector< matrix > r_right_
virtual bool forward(size_t p, size_t q, const CppAD::vector< bool > &vx, CppAD::vector< bool > &vy, const CppAD::vector< scalar > &tx, CppAD::vector< scalar > &ty)
Eigen::Matrix< scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > matrix
Eigen::Matrix< ad_scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > ad_matrix
virtual bool for_sparse_hes(const CppAD::vector< bool > &vx, const CppAD::vector< bool > &r, const CppAD::vector< bool > &s, CppAD::vector< std::set< size_t > > &h, const CppAD::vector< Base > &x)
virtual bool rev_sparse_jac(size_t q, const CppAD::vector< std::set< size_t > > &rt, CppAD::vector< std::set< size_t > > &st, const CppAD::vector< Base > &x)
int Integer(const std::complex< double > &x)
virtual bool for_sparse_jac(size_t q, const CppAD::vector< std::set< size_t > > &r, CppAD::vector< std::set< size_t > > &s, const CppAD::vector< Base > &x)
virtual bool reverse(size_t q, const CppAD::vector< double > &tx, const CppAD::vector< double > &ty, CppAD::vector< double > &px, const CppAD::vector< double > &py)
virtual bool rev_sparse_hes(const CppAD::vector< bool > &vx, const CppAD::vector< bool > &s, CppAD::vector< bool > &t, size_t q, const CppAD::vector< std::set< size_t > > &r, const CppAD::vector< std::set< size_t > > &u, CppAD::vector< std::set< size_t > > &v, const CppAD::vector< Base > &x)
CppAD::AD< scalar > ad_scalar
CppAD::vector< matrix > f_right_