Prev Next Index-> contents reference index search external Up-> CppAD AD ADValued atomic atomic_base atomic_eigen_mat_inv.cpp ADValued-> Arithmetic unary_standard_math binary_math CondExp Discrete numeric_limits atomic atomic-> checkpoint atomic_base atomic_base-> atomic_ctor atomic_option atomic_afun atomic_forward atomic_reverse atomic_for_sparse_jac atomic_rev_sparse_jac atomic_for_sparse_hes atomic_rev_sparse_hes atomic_base_clear atomic_get_started.cpp atomic_norm_sq.cpp atomic_reciprocal.cpp atomic_set_sparsity.cpp atomic_tangent.cpp atomic_eigen_mat_mul.cpp atomic_eigen_mat_inv.cpp atomic_eigen_cholesky.cpp atomic_mat_mul.cpp atomic_eigen_mat_inv.cpp-> atomic_eigen_mat_inv.hpp Headings-> Description Class Definition Use Atomic Function ---..Constructor

$\newcommand{\W}[1]{ \; #1 \; } \newcommand{\R}[1]{ {\rm #1} } \newcommand{\B}[1]{ {\bf #1} } \newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} } \newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} } \newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} } \newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }$
Atomic Eigen Matrix Inverse: Example and Test

Description
The ADFun function object f for this example is $$f(x) = \left( \begin{array}{cc} x_0 & 0 \\ 0 & x_1 \end{array} \right)^{-1} \left( \begin{array}{c} 0 \\ x_2 \end{array} \right) = \left( \begin{array}{c} 0 \\ x_2 / x_1 ) \end{array} \right)$$

Class Definition
This example uses the file atomic_eigen_mat_inv.hpp which defines matrix multiply as a atomic_base operation.

Use Atomic Function
# include <cppad/cppad.hpp>

bool eigen_mat_inv(void)
{
typedef double                                            scalar;
//
bool ok    = true;
scalar eps = 10. * std::numeric_limits<scalar>::epsilon();
//

Constructor
     // -------------------------------------------------------------------
// object that multiplies matrices
atomic_eigen_mat_mul<scalar> mat_mul;
// -------------------------------------------------------------------
// object that computes inverse of a square matrix
atomic_eigen_mat_inv<scalar> mat_inv;
// -------------------------------------------------------------------
// declare independent variable vector x
size_t n = 3;
for(size_t j = 0; j < n; j++)
// -------------------------------------------------------------------
// left = [ x[0]  0    ]
//        [ 0     x[1] ]
size_t nr_left  = 2;
// -------------------------------------------------------------------
// right = [ 0 , x[2] ]^T
size_t nc_right = 1;
// -------------------------------------------------------------------
// use atomic operation to compute left^{-1}
// use atomic operation to multiply left^{-1} * right
// -------------------------------------------------------------------
// declare the dependent variable vector y
size_t m = 2;
for(size_t i = 0; i < m; i++)
// -------------------------------------------------------------------
// check zero order forward mode
for(size_t i = 0; i < n; i++)
x[i] = scalar(i + 2);
y   = f.Forward(0, x);
ok &= NearEqual(y[0], 0.0,          eps, eps);
ok &= NearEqual(y[1], x[2] / x[1],  eps, eps);
// -------------------------------------------------------------------
// check first order forward mode
x1[0] = 1.0;
x1[1] = 0.0;
x1[2] = 0.0;
y1    = f.Forward(1, x1);
ok   &= NearEqual(y1[0], 0.0,        eps, eps);
ok   &= NearEqual(y1[1], 0.0,        eps, eps);
x1[0] = 0.0;
x1[1] = 0.0;
x1[2] = 1.0;
y1    = f.Forward(1, x1);
ok   &= NearEqual(y1[0], 0.0,        eps, eps);
ok   &= NearEqual(y1[1], 1.0 / x[1], eps, eps);
x1[0] = 0.0;
x1[1] = 1.0;
x1[2] = 0.0;
y1    = f.Forward(1, x1);
ok   &= NearEqual(y1[0], 0.0,                  eps, eps);
ok   &= NearEqual(y1[1], - x[2] / (x[1]*x[1]), eps, eps);
// -------------------------------------------------------------------
// check second order forward mode
x2[0] = 0.0;
x2[1] = 0.0;
x2[2] = 0.0;
scalar  f1_x1_x1 = 2.0 * x[2] / (x[1] * x[1] * x[1] );
y2    = f.Forward(2, x2);
ok   &= NearEqual(y2[0], 0.0,            eps, eps);
ok   &= NearEqual(y2[1], f1_x1_x1 / 2.0, eps, eps);
// -------------------------------------------------------------------
// check first order reverse
w[0] = 1.0;
w[1] = 0.0;
d1w  = f.Reverse(1, w);
ok  &= NearEqual(d1w[0], 0.0, eps, eps);
ok  &= NearEqual(d1w[1], 0.0, eps, eps);
ok  &= NearEqual(d1w[2], 0.0, eps, eps);
w[0] = 0.0;
w[1] = 1.0;
d1w  = f.Reverse(1, w);
ok  &= NearEqual(d1w[0], 0.0,                  eps, eps);
ok  &= NearEqual(d1w[1], - x[2] / (x[1]*x[1]), eps, eps);
ok  &= NearEqual(d1w[2], 1.0 / x[1],           eps, eps);
// -------------------------------------------------------------------
// check second order reverse
d2w  = f.Reverse(2, w);
// partial f_1 w.r.t x_0
ok  &= NearEqual(d2w[0 * 2 + 0], 0.0,                  eps, eps);
// partial f_1 w.r.t x_1
ok  &= NearEqual(d2w[1 * 2 + 0], - x[2] / (x[1]*x[1]), eps, eps);
// partial f_1 w.r.t x_2
ok  &= NearEqual(d2w[2 * 2 + 0], 1.0 / x[1],           eps, eps);
// partial f_1 w.r.t x_1, x_0
ok  &= NearEqual(d2w[0 * 2 + 1], 0.0,                  eps, eps);
// partial f_1 w.r.t x_1, x_1
ok  &= NearEqual(d2w[1 * 2 + 1], f1_x1_x1,             eps, eps);
// partial f_1 w.r.t x_1, x_2
ok  &= NearEqual(d2w[2 * 2 + 1], - 1.0 / (x[1]*x[1]),  eps, eps);
// -------------------------------------------------------------------
return ok;
}

Input File: example/atomic/eigen_mat_inv.cpp