source: trunk/test_more/eigen_mat_inv.cpp @ 3811

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

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

commit 33d1f5b8d837d83c8baf656045d85d8b45f0b297
Author: Brad Bell <bradbell@…>
Date: Sun Mar 27 04:01:14 2016 -0700

  1. Test non-symmetric reverse mode in eigen_mat_inv (and fix this example).
  2. Advance version to cppad-20160327.

commit 238e214ff36ca835efca615c609576dc8bf5038d
Author: Brad Bell <bradbell@…>
Date: Sat Mar 26 18:10:59 2016 -0700

eigen_mat_inv.hpp: use inverse matrix (since we are computing it).
eigen_mat_mul.hpp: comment seperator.
eigen_mat_inv.cpp: matrix during recording is non-singular.
eigen_mat_inv.cpp: test first and second order forward.

commit fc918b0476cc8ea66abdf2904a71fc93472d279d
Author: Brad Bell <bradbell@…>
Date: Sat Mar 26 16:26:51 2016 -0700

Add test_more/eigen_mat_inv.cpp.


eigen_mat_inv.hpp: fix calculation of vy.
eigen_mat_inv.cpp: fix section title.
CMakeLists.txt: no longer necessary to have special library for eigen tests.
eigen_mat_inv.cpp: Test with non-symmetric matrix.

File size: 6.0 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/*
13Test atomic_eigen_mat_inv using a non-symetric matrix
14
15f(x) = [ x[0]   -1  ]^{-1}
16       [ 2     x[1] ]
17
18     = [ x[1]  1    ] / (x[0] * x[1] + 2)
19       [ -2    x[0] ]
20
21y[0] = x[1] / (x[0] * x[1] + 2)
22y[1] = 1.0  / (x[0] * x[1] + 2)
23y[2] = -2.0 / (x[0] * x[1] + 2)
24y[3] = x[0] / (x[0] * x[1] + 2)
25*/
26# include <cppad/cppad.hpp>
27# include <cppad/example/eigen_mat_inv.hpp>
28
29
30bool eigen_mat_inv(void)
31{
32        typedef double                                            scalar;
33        typedef CppAD::AD<scalar>                                 ad_scalar;
34        typedef typename atomic_eigen_mat_inv<scalar>::ad_matrix  ad_matrix;
35        //
36        bool ok    = true;
37        scalar eps = 10. * std::numeric_limits<scalar>::epsilon();
38        using CppAD::NearEqual;
39        // -------------------------------------------------------------------
40        // object that computes invers of a 2x2 matrix
41        size_t nr  = 2;
42        atomic_eigen_mat_inv<scalar> mat_inv(nr);
43        // -------------------------------------------------------------------
44        // declare independent variable vector x
45        size_t n = 2;
46        CPPAD_TESTVECTOR(ad_scalar) ad_x(n);
47        for(size_t j = 0; j < n; j++)
48                ad_x[j] = ad_scalar(j);
49        CppAD::Independent(ad_x);
50        // -------------------------------------------------------------------
51        // arg = [ x[0]  -1   ]
52        //       [ 2     x[1] ]
53        ad_matrix ad_arg(nr, nr);
54        ad_arg(0, 0) = ad_x[0];
55        ad_arg(0, 1) = ad_scalar(-1.0);
56        ad_arg(1, 0) = ad_scalar(2.0);
57        ad_arg(1, 1) = ad_x[1];
58        // -------------------------------------------------------------------
59        // use atomic operation to compute arg^{-1}
60        ad_matrix ad_result = mat_inv.op(ad_arg);
61        // -------------------------------------------------------------------
62        // declare the dependent variable vector y
63        size_t m = 4;
64        CPPAD_TESTVECTOR(ad_scalar) ad_y(4);
65        for(size_t i = 0; i < nr; i++)
66                for(size_t j = 0; j < nr; j++)
67                        ad_y[ i * nr + j ] = ad_result(i, j);
68        /* Used to test hand calculated derivaives
69        CppAD::AD<scalar> ad_dinv = 1.0 / (ad_x[0] * ad_x[1] + 2.0);
70        ad_y[0] = ad_x[1] * ad_dinv;
71        ad_y[1] = 1.0  * ad_dinv;
72        ad_y[2] = -2.0 * ad_dinv;
73        ad_y[3] = ad_x[0] * ad_dinv;
74        */
75        CppAD::ADFun<scalar> f(ad_x, ad_y);
76        // -------------------------------------------------------------------
77        // check zero order forward mode
78        CPPAD_TESTVECTOR(scalar) x(n), y(m);
79        for(size_t i = 0; i < n; i++)
80                x[i] = scalar(i + 2);
81        scalar dinv = 1.0 / (x[0] * x[1] + 2.0);
82        y          = f.Forward(0, x);
83        ok        &= NearEqual(y[0], x[1] * dinv,  eps, eps);
84        ok        &= NearEqual(y[1], 1.0  * dinv,  eps, eps);
85        ok        &= NearEqual(y[2], -2.0 * dinv,  eps, eps);
86        ok        &= NearEqual(y[3], x[0] * dinv,  eps, eps);
87        // -------------------------------------------------------------------
88        // check first order forward mode
89        CPPAD_TESTVECTOR(scalar) x1(n), y1(m);
90        scalar dinv_x0 = - x[1] * dinv * dinv;
91        x1[0] = 1.0;
92        x1[1] = 0.0;
93        y1    = f.Forward(1, x1);
94        ok   &= NearEqual(y1[0], x[1] * dinv_x0,        eps, eps);
95        ok   &= NearEqual(y1[1], 1.0  * dinv_x0,        eps, eps);
96        ok   &= NearEqual(y1[2], -2.0 * dinv_x0,        eps, eps);
97        ok   &= NearEqual(y1[3], dinv + x[0] * dinv_x0, eps, eps);
98        //
99        scalar dinv_x1 = - x[0] * dinv * dinv;
100        x1[0] = 0.0;
101        x1[1] = 1.0;
102        y1    = f.Forward(1, x1);
103        ok   &= NearEqual(y1[0], dinv + x[1] * dinv_x1, eps, eps);
104        ok   &= NearEqual(y1[1], 1.0  * dinv_x1,        eps, eps);
105        ok   &= NearEqual(y1[2], -2.0 * dinv_x1,        eps, eps);
106        ok   &= NearEqual(y1[3], x[0] * dinv_x1,        eps, eps);
107        // -------------------------------------------------------------------
108        // check second order forward mode
109        CPPAD_TESTVECTOR(scalar) x2(n), y2(m);
110        scalar dinv_x1_x1 = 2.0 * x[0] * x[0] * dinv * dinv * dinv;
111        x2[0] = 0.0;
112        x2[1] = 0.0;
113        y2    = f.Forward(2, x2);
114        ok   &= NearEqual(2.0*y2[0], 2.0*dinv_x1 + x[1]*dinv_x1_x1, eps, eps);
115        ok   &= NearEqual(2.0*y2[1], 1.0  * dinv_x1_x1,             eps, eps);
116        ok   &= NearEqual(2.0*y2[2], -2.0 * dinv_x1_x1,             eps, eps);
117        ok   &= NearEqual(2.0*y2[3], x[0] * dinv_x1_x1,             eps, eps);
118        // -------------------------------------------------------------------
119        // check first order reverse
120        CPPAD_TESTVECTOR(scalar) w(m), d1w(n);
121        for(size_t i = 0; i < m; i++)
122                w[i] = 0.0;
123        w[0] = 1.0;
124        d1w  = f.Reverse(1, w);
125        ok  &= NearEqual(d1w[0], x[1] * dinv_x0,        eps, eps);
126        ok  &= NearEqual(d1w[1], dinv + x[1] * dinv_x1, eps, eps);
127        w[0] = 0.0;
128        w[1] = 1.0;
129        d1w  = f.Reverse(1, w);
130        ok  &= NearEqual(d1w[0], 1.0 * dinv_x0,         eps, eps);
131        ok  &= NearEqual(d1w[1], 1.0 * dinv_x1,         eps, eps);
132        w[1] = 0.0;
133        w[2] = 1.0;
134        d1w  = f.Reverse(1, w);
135        ok  &= NearEqual(d1w[0], -2.0 * dinv_x0,        eps, eps);
136        ok  &= NearEqual(d1w[1], -2.0 * dinv_x1,        eps, eps);
137        w[2] = 0.0;
138        w[3] = 1.0;
139        d1w  = f.Reverse(1, w);
140        ok  &= NearEqual(d1w[0], dinv + x[0] * dinv_x0, eps, eps);
141        ok  &= NearEqual(d1w[1], x[0] * dinv_x1,        eps, eps);
142        // -------------------------------------------------------------------
143        // check second order reverse
144        CPPAD_TESTVECTOR(scalar) d2w(2 * n);
145        // dinv_x1 = - x[0] * dinv * dinv;
146        scalar dinv_x1_x0 = 2.0 * x[0] * x[1] * dinv * dinv * dinv - dinv * dinv;
147        d2w  = f.Reverse(2, w);
148        // partial f_3 w.r.t x_0
149        ok  &= NearEqual(d2w[0 * 2 + 0], dinv + x[0] * dinv_x0, eps, eps);
150        // partial f_3 w.r.t x_1
151        ok  &= NearEqual(d2w[1 * 2 + 0], x[0] * dinv_x1,              eps, eps);
152        // partial f_3 w.r.t. x_1, x_0
153        ok  &= NearEqual(d2w[0 * 2 + 1], dinv_x1 + x[0] * dinv_x1_x0, eps, eps);
154        // partial f_3 w.r.t. x_1, x_1
155        ok  &= NearEqual(d2w[1 * 2 + 1], x[0] * dinv_x1_x1,           eps, eps);
156        // -------------------------------------------------------------------
157        return ok;
158}
Note: See TracBrowser for help on using the repository browser.