2 # ifndef CPPAD_EXAMPLE_EIGEN_CHOLESKY_HPP
3 # define CPPAD_EXAMPLE_EIGEN_CHOLESKY_HPP
36 # include <Eigen/Dense>
56 typedef Eigen::Matrix<
59 typedef Eigen::Matrix<
63 typedef Eigen::TriangularView<matrix, Eigen::Lower>
lower_view;
69 "atom_eigen_cholesky" ,
70 CppAD::atomic_base<Base>::set_sparsity_enum
78 {
size_t nr = size_t( arg.rows() );
79 size_t ny = ( (nr + 1 ) * nr ) / 2;
81 assert( nr ==
size_t( arg.cols() ) );
84 CPPAD_TESTVECTOR(
ad_scalar) packed_arg(nx);
88 for(
size_t i = 0; i < nr; i++)
89 {
for(
size_t j = 0; j <= i; j++)
90 packed_arg[index++] = arg(i, j);
92 assert( index == nx );
97 CPPAD_TESTVECTOR(
ad_scalar) packed_result(ny);
98 (*this)(packed_arg, packed_result);
101 ad_matrix result = ad_matrix::Zero(nr, nr);
103 for(
size_t i = 0; i < nr; i++)
104 {
for(
size_t j = 0; j <= i; j++)
105 result(i, j) = packed_result[index++];
139 {
size_t n_order = q + 1;
141 size_t ny = ((nr + 1) * nr) / 2;
145 assert( vx.
size() == 0 || nx == vx.
size() );
146 assert( vx.
size() == 0 || ny == vy.
size() );
147 assert( nx * n_order == tx.
size() );
148 assert( ny * n_order == ty.
size() );
152 assert( f_arg_.size() == f_result_.size() );
153 if( f_arg_.size() < n_order )
154 { f_arg_.resize(n_order);
155 f_result_.resize(n_order);
157 for(
size_t k = 0; k < n_order; k++)
158 { f_arg_[k].resize(nr, nr);
159 f_result_[k].resize(nr, nr);
164 for(
size_t k = 0; k < n_order; k++)
167 for(
size_t i = 0; i < nr; i++)
168 {
for(
size_t j = 0; j <= i; j++)
169 { f_arg_[k](i, j) = tx[ index * n_order + k ];
170 f_arg_[k](j, i) = f_arg_[k](i, j);
179 Eigen::LLT<matrix> cholesky(f_arg_[0]);
180 f_result_[0] = cholesky.matrixL();
181 lower_view L_0 = f_result_[0].template triangularView<Eigen::Lower>();
182 for(
size_t k = 1; k < n_order; k++)
186 for(
size_t ell = 1; ell < k; ell++)
187 f_sum -= f_result_[ell] * f_result_[k-ell].transpose();
189 matrix temp = L_0.template solve<Eigen::OnTheLeft>(f_sum);
190 temp = L_0.transpose().template solve<Eigen::OnTheRight>(temp);
192 for(
size_t i = 0; i < nr; i++)
193 temp(i, i) /=
scalar(2.0);
195 lower_view view = temp.template triangularView<Eigen::Lower>();
196 f_result_[k] = f_result_[0] * view;
200 for(
size_t k = 0; k < n_order; k++)
202 for(
size_t i = 0; i < nr; i++)
203 {
for(
size_t j = 0; j <= i; j++)
204 { ty[ index * n_order + k ] = f_result_[k](i, j);
217 for(
size_t i = 0; i < ny; i++)
219 for(
size_t i = 0; i < ny; i++)
240 {
size_t n_order = q + 1;
243 size_t ny = ( (nr + 1 ) * nr ) / 2;
247 assert( nx * n_order == tx.
size() );
248 assert( ny * n_order == ty.
size() );
253 assert( f_arg_.size() == f_result_.size() );
255 assert( f_arg_.size() >= n_order );
258 assert( r_arg_.size() == r_result_.size() );
259 if( r_arg_.size() < n_order )
260 { r_arg_.resize(n_order);
261 r_result_.resize(n_order);
263 for(
size_t k = 0; k < n_order; k++)
264 { r_arg_[k].resize(nr, nr);
265 r_result_[k].resize(nr, nr);
270 for(
size_t k = 0; k < n_order; k++)
273 for(
size_t i = 0; i < nr; i++)
274 {
for(
size_t j = 0; j <= i; j++)
275 { f_arg_[k](i, j) = tx[ index * n_order + k ];
276 f_arg_[k](j, i) = f_arg_[k](i, j);
283 for(
size_t k = 0; k < n_order; k++)
284 { r_result_[k] = matrix::Zero(nr, nr);
286 for(
size_t i = 0; i < nr; i++)
287 {
for(
size_t j = 0; j <= i; j++)
288 { r_result_[k](i, j) = py[ index * n_order + k ];
295 for(
size_t k = 0; k < n_order; k++)
296 r_arg_[k] = matrix::Zero(nr, nr);
299 lower_view L_0 = f_result_[0].template triangularView<Eigen::Lower>();
301 for(
size_t k1 = n_order; k1 > 1; k1--)
305 matrix tmp1 = L_0.transpose() * r_result_[k];
308 for(
size_t i = 0; i < nr; i++)
309 tmp1(i, i) /=
scalar(2.0);
310 matrix tmp2 = tmp1.template triangularView<Eigen::Lower>();
313 tmp1 = L_0.transpose().template solve<Eigen::OnTheLeft>( tmp2 );
316 matrix M_k = L_0.transpose().template
317 solve<Eigen::OnTheLeft>( tmp1.transpose() );
320 matrix barB_k =
scalar(0.5) * ( M_k + M_k.transpose() );
322 barB_k =
scalar(-1.0) * barB_k;
326 temp = temp.template triangularView<Eigen::Lower>();
329 r_result_[0] += temp;
332 for(
size_t ell = 1; ell < k; ell++)
334 temp =
scalar(2.0) * barB_k * f_result_[k-ell];
335 r_result_[ell] += temp.template triangularView<Eigen::Lower>();
339 matrix M_0 = L_0.transpose() * r_result_[0];
340 for(
size_t i = 0; i < nr; i++)
342 M_0 = M_0.template triangularView<Eigen::Lower>();
343 M_0 = L_0.template solve<Eigen::OnTheRight>( M_0 );
344 M_0 = L_0.transpose().template solve<Eigen::OnTheLeft>( M_0 );
346 r_arg_[0] +=
scalar(0.5) * ( M_0 + M_0.transpose() );
350 for(
size_t k = 0; k < n_order; k++)
352 px[ index * n_order + k ] = 0.0;
354 for(
size_t i = 0; i < nr; i++)
355 {
for(
size_t j = 0; j < i; j++)
356 { px[ index * n_order + k ] = 2.0 * r_arg_[k](i, j);
359 px[ index * n_order + k] = r_arg_[k](i, i);
includes the entire CppAD package in the necessary order.
Eigen::TriangularView< matrix, Eigen::Lower > lower_view
atomic_eigen_cholesky(void)
Eigen::Matrix< ad_scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > ad_matrix
CppAD::vector< matrix > f_result_
ad_matrix op(const ad_matrix &arg)
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.
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 > r_result_
int Integer(const std::complex< double > &x)
Eigen::Matrix< scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > matrix
CppAD::AD< scalar > ad_scalar