source: trunk/example/atomic/eigen_mat_inv.cpp @ 3810

Last change on this file since 3810 was 3810, checked in by bradbell, 4 years ago

merge to branch: trunk
from repository: https://github.com/coin-or/CppAD
start hash code: dd967ef41b8d6731d90ebb8b3e7e8b863565289c
end hash code: f77cd2e443b8d14f354f97efc4c26369227dfe81

commit f77cd2e443b8d14f354f97efc4c26369227dfe81
Author: Brad Bell <bradbell@…>
Date: Sat Mar 26 06:47:33 2016 -0700

  1. Add reverse mode to atomic_eigen_mat_inv.
  2. Advance version to cppad-20160326.


eigen_mat_mul.cpp: impove reverse mode comments.

commit 94c81399191d06ae9cef81c271f520b0e78f079d
Author: Brad Bell <bradbell@…>
Date: Fri Mar 25 12:16:03 2016 -0700

whats_new_16.omh: add comment about changes to eigen_mat_mul.

commit d647288c9e022f16ef94229db612a774808a5de3
Author: Brad Bell <bradbell@…>
Date: Fri Mar 25 11:59:45 2016 -0700

Replace pack and unpack by member function that inverts matrices.

commit 6fba95d90e50590abd438946f47025e0a6e34c2c
Author: Brad Bell <bradbell@…>
Date: Fri Mar 25 11:45:43 2016 -0700

Replace pack and unpack by member function that multiplies matrices.

commit cac97c115abad9fba0d0c3706f4bc0e5a7f617f0
Author: Brad Bell <bradbell@…>
Date: Fri Mar 25 11:11:16 2016 -0700

Change atomic_eigen_mat_div -> atomic_eigen_mat_inv.

commit 1f7e52a4d96986f7be2ad2d6abeafc086c4f7622
Author: Brad Bell <bradbell@…>
Date: Fri Mar 25 08:09:39 2016 -0700

  1. Implement and test forward for atomic_eigen_mat_div.
  2. Advance version to 20160325.


eigen_mat_mul.hpp: edit comments and documentation.

commit e953f669d8e9e122dc2f182c28c051c9ba75c705
Author: Brad Bell <bradbell@…>
Date: Fri Mar 25 04:32:45 2016 -0700

eigen_mat_mul.cpp: No longer under construction.

commit 3a62044eaebb6393676276e5b843fbe5db4f45a2
Author: Brad Bell <bradbell@…>
Date: Fri Mar 25 04:07:18 2016 -0700

eigen_mat_mul.hpp: add assert for f_left_ and f_right_ in reverse.

commit 721d64783547487fe69100c7190628efab7149c3
Author: Brad Bell <bradbell@…>
Date: Fri Mar 25 03:59:05 2016 -0700

eigen_mat_mul.hpp: typo Publice -> Public.

File size: 6.3 KB
Line 
1// $Id$
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-16 Bradley M. Bell
4
5CppAD is distributed under multiple licenses. This distribution is under
6the terms of the
7                    Eclipse Public License Version 1.0.
8
9A copy of this license is included in the COPYING file of this distribution.
10Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
11-------------------------------------------------------------------------- */
12
13/*
14$begin atomic_eigen_mat_inv.cpp$$
15$spell
16        mul
17        Eigen
18$$
19
20$section  Atomic Eigen Matrix Division: Example and Test$$
21
22$head Description$$
23The $cref ADFun$$ function object $icode f$$ for this example is
24$latex \[
25f(x) =
26\left( \begin{array}{cc}
27        x_0   & 0 \\
28        0   & x_1
29\end{array} \right)^{-1}
30\left( \begin{array}{c}
31        0   \\
32        x_2
33\end{array} \right)
34=
35\left( \begin{array}{c}
36        0 \\
37        x_2 / x_1 )
38\end{array} \right)
39\] $$
40
41$children%
42        cppad/example/eigen_mat_inv.hpp
43%$$
44
45$head Class Definition$$
46This example uses the file $cref atomic_eigen_mat_inv.hpp$$
47which defines matrix multiply as a $cref atomic_base$$ operation.
48
49$nospell
50
51$head Use Atomic Function$$
52$srccode%cpp% */
53# include <cppad/cppad.hpp>
54# include <cppad/example/eigen_mat_inv.hpp>
55# include <cppad/example/eigen_mat_mul.hpp>
56
57
58bool eigen_mat_inv(void)
59{
60        typedef double                                            scalar;
61        typedef CppAD::AD<scalar>                                 ad_scalar;
62        typedef typename atomic_eigen_mat_inv<scalar>::ad_matrix  ad_matrix;
63        //
64        bool ok    = true;
65        scalar eps = 10. * std::numeric_limits<scalar>::epsilon();
66        using CppAD::NearEqual;
67        //
68/* %$$
69$subhead Constructor$$
70$srccode%cpp% */
71        // -------------------------------------------------------------------
72        // object that multiplies a 2x2 matrix times a 2x1 matrix
73        size_t nr_left  = 2;
74        size_t n_middle = 2;
75        size_t nc_right = 1;
76        atomic_eigen_mat_mul<scalar> mat_mul(nr_left, n_middle, nc_right);
77        // -------------------------------------------------------------------
78        // object that computes invers of a 2x2 matrix
79        size_t nr  = 2;
80        atomic_eigen_mat_inv<scalar> mat_inv(nr);
81        // -------------------------------------------------------------------
82        // declare independent variable vector x
83        size_t n = 3;
84        CPPAD_TESTVECTOR(ad_scalar) ad_x(n);
85        for(size_t j = 0; j < n; j++)
86                ad_x[j] = ad_scalar(j);
87        CppAD::Independent(ad_x);
88        // -------------------------------------------------------------------
89        // left = [ x[0]  0    ]
90        //        [ 0     x[1] ]
91        ad_matrix ad_left(nr_left, nr_left);
92        ad_left(0, 0) = ad_x[0];
93        ad_left(0, 1) = ad_scalar(0.0);
94        ad_left(1, 0) = ad_scalar(0.0);
95        ad_left(1, 1) = ad_x[1];
96        // -------------------------------------------------------------------
97        // right = [ 0 , x[2] ]^T
98        ad_matrix ad_right(nr_left, nc_right);
99        ad_right(0, 0) = ad_scalar(0.0);
100        ad_right(1, 0) = ad_x[2];
101        // -------------------------------------------------------------------
102        // use atomic operation to compute left^{-1}
103        ad_matrix ad_left_inv = mat_inv.op(ad_left);
104        // use atomic operation to multiply left^{-1} * right
105        ad_matrix ad_result   = mat_mul.op(ad_left_inv, ad_right);
106        // -------------------------------------------------------------------
107        // check that first component of result is a parameter
108        // and the second component is a varaible.
109        ok &= Parameter( ad_result(0, 0) );
110        ok &= Variable(  ad_result(1, 0) );
111        // -------------------------------------------------------------------
112        // declare the dependent variable vector y
113        size_t m = 2;
114        CPPAD_TESTVECTOR(ad_scalar) ad_y(2);
115        for(size_t i = 0; i < m; i++)
116                ad_y[i] = ad_result(i, 0);
117        CppAD::ADFun<scalar> f(ad_x, ad_y);
118        // -------------------------------------------------------------------
119        // check zero order forward mode
120        CPPAD_TESTVECTOR(scalar) x(n), y(m);
121        for(size_t i = 0; i < n; i++)
122                x[i] = scalar(i + 2);
123        y   = f.Forward(0, x);
124        ok &= NearEqual(y[0], 0.0,          eps, eps);
125        ok &= NearEqual(y[1], x[2] / x[1],  eps, eps);
126        // -------------------------------------------------------------------
127        // check first order forward mode
128        CPPAD_TESTVECTOR(scalar) x1(n), y1(m);
129        x1[0] = 1.0;
130        x1[1] = 0.0;
131        x1[2] = 0.0;
132        y1    = f.Forward(1, x1);
133        ok   &= NearEqual(y1[0], 0.0,        eps, eps);
134        ok   &= NearEqual(y1[1], 0.0,        eps, eps);
135        x1[0] = 0.0;
136        x1[1] = 0.0;
137        x1[2] = 1.0;
138        y1    = f.Forward(1, x1);
139        ok   &= NearEqual(y1[0], 0.0,        eps, eps);
140        ok   &= NearEqual(y1[1], 1.0 / x[1], eps, eps);
141        x1[0] = 0.0;
142        x1[1] = 1.0;
143        x1[2] = 0.0;
144        y1    = f.Forward(1, x1);
145        ok   &= NearEqual(y1[0], 0.0,                  eps, eps);
146        ok   &= NearEqual(y1[1], - x[2] / (x[1]*x[1]), eps, eps);
147        // -------------------------------------------------------------------
148        // check second order forward mode
149        CPPAD_TESTVECTOR(scalar) x2(n), y2(m);
150        x2[0] = 0.0;
151        x2[1] = 0.0;
152        x2[2] = 0.0;
153        scalar  f1_x1_x1 = 2.0 * x[2] / (x[1] * x[1] * x[1] );
154        y2    = f.Forward(2, x2);
155        ok   &= NearEqual(y2[0], 0.0,            eps, eps);
156        ok   &= NearEqual(y2[1], f1_x1_x1 / 2.0, eps, eps);
157        // -------------------------------------------------------------------
158        // check first order reverse
159        CPPAD_TESTVECTOR(scalar) w(m), d1w(n);
160        w[0] = 1.0;
161        w[1] = 0.0;
162        d1w  = f.Reverse(1, w);
163        ok  &= NearEqual(d1w[0], 0.0, eps, eps);
164        ok  &= NearEqual(d1w[1], 0.0, eps, eps);
165        ok  &= NearEqual(d1w[2], 0.0, eps, eps);
166        w[0] = 0.0;
167        w[1] = 1.0;
168        d1w  = f.Reverse(1, w);
169        ok  &= NearEqual(d1w[0], 0.0,                  eps, eps);
170        ok  &= NearEqual(d1w[1], - x[2] / (x[1]*x[1]), eps, eps);
171        ok  &= NearEqual(d1w[2], 1.0 / x[1],           eps, eps);
172        // -------------------------------------------------------------------
173        // check second order reverse
174        CPPAD_TESTVECTOR(scalar) d2w(2 * n);
175        d2w  = f.Reverse(2, w);
176        // partial f_1 w.r.t x_0
177        ok  &= NearEqual(d2w[0 * 2 + 0], 0.0,                  eps, eps);
178        // partial f_1 w.r.t x_1
179        ok  &= NearEqual(d2w[1 * 2 + 0], - x[2] / (x[1]*x[1]), eps, eps);
180        // partial f_1 w.r.t x_2
181        ok  &= NearEqual(d2w[2 * 2 + 0], 1.0 / x[1],           eps, eps);
182        // partial f_1 w.r.t x_1, x_0
183        ok  &= NearEqual(d2w[0 * 2 + 1], 0.0,                  eps, eps);
184        // partial f_1 w.r.t x_1, x_1
185        ok  &= NearEqual(d2w[1 * 2 + 1], f1_x1_x1,             eps, eps);
186        // partial f_1 w.r.t x_1, x_2
187        ok  &= NearEqual(d2w[2 * 2 + 1], - 1.0 / (x[1]*x[1]),  eps, eps);
188        // -------------------------------------------------------------------
189        return ok;
190}
191/* %$$
192$$ $comment end nospell$$
193$end
194*/
Note: See TracBrowser for help on using the repository browser.