2 # ifndef CPPAD_EXAMPLE_EIGEN_MAT_INV_HPP
3 # define CPPAD_EXAMPLE_EIGEN_MAT_INV_HPP
137 # include <Eigen/Core>
149 template <
class Base>
158 typedef Eigen::Matrix<
161 typedef Eigen::Matrix<
168 "atom_eigen_mat_inv" ,
169 CppAD::atomic_base<Base>::set_sparsity_enum
177 {
size_t nr = size_t( arg.rows() );
180 assert( nr ==
size_t( arg.cols() ) );
183 CPPAD_TESTVECTOR(
ad_scalar) packed_arg(nx);
185 for(
size_t i = 0; i < ny; i++)
186 packed_arg[1 + i] = arg.data()[i];
191 CPPAD_TESTVECTOR(
ad_scalar) packed_result(ny);
192 (*this)(packed_arg, packed_result);
196 for(
size_t i = 0; i < ny; i++)
197 result.data()[i] = packed_result[i];
230 {
size_t n_order = q + 1;
236 assert( vx.
size() == 0 || nx == vx.
size() );
237 assert( vx.
size() == 0 || ny == vy.
size() );
238 assert( nx * n_order == tx.
size() );
239 assert( ny * n_order == ty.
size() );
243 assert( f_arg_.size() == f_result_.size() );
244 if( f_arg_.size() < n_order )
245 { f_arg_.resize(n_order);
246 f_result_.resize(n_order);
248 for(
size_t k = 0; k < n_order; k++)
249 { f_arg_[k].resize(nr, nr);
250 f_result_[k].resize(nr, nr);
255 for(
size_t k = 0; k < n_order; k++)
257 for(
size_t i = 0; i < ny; i++)
258 f_arg_[k].data()[i] = tx[ (1 + i) * n_order + k ];
264 f_result_[0] = f_arg_[0].inverse();
265 for(
size_t k = 1; k < n_order; k++)
267 matrix f_sum = matrix::Zero(nr, nr);
269 for(
size_t ell = 1; ell <= k; ell++)
270 f_sum -= f_arg_[ell] * f_result_[k-ell];
272 f_result_[k] = f_result_[0] * f_sum;
276 for(
size_t k = 0; k < n_order; k++)
277 {
for(
size_t i = 0; i < ny; i++)
278 ty[ i * n_order + k ] = f_result_[k].data()[i];
288 for(
size_t i = 0; i < ny; i++)
290 for(
size_t i = 0; i < ny; i++)
310 {
size_t n_order = q + 1;
317 assert( nx * n_order == tx.
size() );
318 assert( ny * n_order == ty.
size() );
323 assert( f_arg_.size() == f_result_.size() );
325 assert( f_arg_.size() >= n_order );
328 assert( r_arg_.size() == r_result_.size() );
329 if( r_arg_.size() < n_order )
330 { r_arg_.resize(n_order);
331 r_result_.resize(n_order);
333 for(
size_t k = 0; k < n_order; k++)
334 { r_arg_[k].resize(nr, nr);
335 r_result_[k].resize(nr, nr);
340 for(
size_t k = 0; k < n_order; k++)
342 for(
size_t i = 0; i < ny; i++)
343 f_arg_[k].data()[i] = tx[ (1 + i) * n_order + k ];
347 for(
size_t k = 0; k < n_order; k++)
348 {
for(
size_t i = 0; i < ny; i++)
349 r_result_[k].data()[i] = py[ i * n_order + k ];
353 for(
size_t k = 0; k < n_order; k++)
354 r_arg_[k] = matrix::Zero(nr, nr);
358 for(
size_t k1 = n_order; k1 > 1; k1--)
362 r_result_[k] * f_result_[k].transpose() * f_arg_[0].transpose();
364 for(
size_t ell = 1; ell <= k; ell++)
366 r_arg_[ell] -= f_result_[0].transpose()
367 * r_result_[k] * f_result_[k-ell].transpose();
369 r_result_[k-ell] -= f_arg_[ell].transpose()
370 * f_result_[0].transpose() * r_result_[k];
374 f_result_[0].transpose() * r_result_[0] * f_result_[0].transpose();
377 for(
size_t k = 0; k < n_order; k++)
378 {
for(
size_t i = 0; i < ny; i++)
379 px[ (1 + i) * n_order + k ] = r_arg_[k].
data()[i];
includes the entire CppAD package in the necessary order.
CppAD::vector< matrix > r_result_
Eigen::Matrix< scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > matrix
Type * data(void)
raw pointer to the data
ad_matrix op(const ad_matrix &arg)
CppAD::AD< scalar > ad_scalar
Eigen::Matrix< ad_scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > ad_matrix
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)
size_t size(void) const
number of elements currently in this vector.
int Integer(const std::complex< double > &x)
atomic_eigen_mat_inv(void)
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)
CppAD::vector< matrix > f_result_