source: trunk/example/atomic/eigen_mat_mul.cpp @ 3808

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

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

commit f8fc499c5720e9bead5d00b92db3d775c10e7043
Author: Brad Bell <bradbell@…>
Date: Thu Mar 24 05:05:58 2016 -0700

  1. Fix bug building example/atomic (conditional inclusion of eigen_mat_mul).
  2. Advance version to cppad-20160324.


run_cmake.sh: add conditional inclusion of eigen.
CMakeLists.txt: eigen_libs no longer necessary (using system incldue for eigen).
eigen_prefix.omh: update list of eigen examples.

commit 0021492a9b151fa056a98b1b9cbf16dfae39f6be
Author: Brad Bell <bradbell@…>
Date: Thu Mar 24 03:55:50 2016 -0700

eigen_mat_mul.hpp: move resize in reverse to beginning, fix reverse heading.

commit b6a96073abd0355b7cab1beb18cd3b2720ea08d5
Author: Brad Bell <bradbell@…>
Date: Thu Mar 24 03:47:14 2016 -0700

eigen_mat_mul.hpp: use data() and Zero to simplify reverse.

commit f7ff6b20162d88f39361441c38ab224f3fd2b929
Author: Brad Bell <bradbell@…>
Date: Thu Mar 24 03:38:12 2016 -0700

eigen_mat_mul.hpp: use data() to simplify forward.

commit 13fef46531f8fce5551e678955c0b0c63e564497
Author: Brad Bell <bradbell@…>
Date: Thu Mar 24 03:26:57 2016 -0700

eigen_mat_mul.hpp: use data() to simplify public pack and unpack.

commit 4cbd54b0775431a1b7ab2b4ff907593ad895ab02
Author: Brad Bell <bradbell@…>
Date: Wed Mar 23 18:26:32 2016 -0700

eigen_mat_mul.cpp: implement reverse, test first order.

commit 038c681ad17930e62faff713acf85993fdbd2921
Author: Brad Bell <bradbell@…>
Date: Wed Mar 23 17:08:08 2016 -0700

eigen_mat_mul.hpp: working on reverse mode.

commit b17b981270cab200abab07c55783d3fb02278e70
Author: Brad Bell <bradbell@…>
Date: Wed Mar 23 09:35:30 2016 -0700

eigen_mat_mul.hpp: more preparation for reverse mode.

commit 2f9f4633dde11cd004760aea36f341f45cafff37
Author: Brad Bell <bradbell@…>
Date: Wed Mar 23 08:43:39 2016 -0700

eigen_mat_mul.hpp: more preparation for reverse mode.

commit c75d7be73c89c06286931eb27904ae541415d4a4
Author: Brad Bell <bradbell@…>
Date: Wed Mar 23 08:06:58 2016 -0700

eigen_mat_mul.hpp: prepare for reverse mode.

commit c5321490d8cb8bbd1563e847a97c0813fd631296
Author: Brad Bell <bradbell@…>
Date: Wed Mar 23 07:07:08 2016 -0700

eigen_mat_mul.cpp: improve this example.

File size: 5.8 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_mul.cpp$$
15$spell
16        mul
17        Eigen
18$$
19
20$section  Atomic Eigen Matrix Multiply: Example and Test$$
21
22$head Under Construction$$
23This example is under construction. So far only forward mode has been
24implemented.
25
26$head Description$$
27The $cref ADFun$$ function object $icode f$$ for this example is
28$latex \[
29f(x) =
30\left( \begin{array}{cc}
31        0 & 0 \\
32        1 & 2 \\
33        2 & 4
34\end{array} \right)
35\left( \begin{array}{c}
36        x_0 \\
37        x_1
38\end{array} \right)
39=
40\left( \begin{array}{c}
41        0 \\
42        x_0 + 2 x_1 \\
43        2 x_0 + 4 x_1 )
44\end{array} \right)
45\] $$
46
47$children%
48        cppad/example/eigen_mat_mul.hpp
49%$$
50
51$head Class Definition$$
52This example uses the file $cref atomic_eigen_mat_mul.hpp$$
53which defines matrix multiply as a $cref atomic_base$$ operation.
54
55$nospell
56
57$head Use Atomic Function$$
58$srccode%cpp% */
59# include <cppad/cppad.hpp>
60# include <cppad/example/eigen_mat_mul.hpp>
61
62namespace {
63        typedef double            scalar;
64        typedef CppAD::AD<scalar> ad_scalar;
65        typedef typename atomic_eigen_mat_mul<scalar>::matrix     matrix;
66        typedef typename atomic_eigen_mat_mul<scalar>::ad_matrix  ad_matrix;
67
68        // use atomic operation to multiply two AD matrices
69        ad_matrix matrix_multiply(
70                atomic_eigen_mat_mul<scalar>& mat_mul ,
71                const ad_matrix&              left    ,
72                const ad_matrix&              right   )
73        {       size_t nr_left   = size_t( left.rows() );
74                size_t n_middle    = size_t( left.cols() );
75                size_t nc_right  = size_t ( right.cols() );
76                assert( size_t( right.rows() ) == n_middle );
77
78                // packed version of left and right
79                size_t nx = (nr_left + nc_right) * n_middle;
80                CPPAD_TESTVECTOR(ad_scalar) packed_arg(nx);
81                mat_mul.pack(packed_arg, left, right);
82
83                // packed version of result = left * right
84                size_t ny = nr_left * nc_right;
85                CPPAD_TESTVECTOR(ad_scalar) packed_result(ny);
86                mat_mul(packed_arg, packed_result);
87
88                // result matrix
89                ad_matrix result(nr_left, nc_right);
90                mat_mul.unpack(packed_result, result);
91
92                return result;
93        }
94
95}
96
97bool eigen_mat_mul(void)
98{       bool ok    = true;
99        scalar eps = 10. * std::numeric_limits<scalar>::epsilon();
100        using CppAD::NearEqual;
101        //
102/* %$$
103$subhead Constructor$$
104$srccode%cpp% */
105        // -------------------------------------------------------------------
106        // object that multiplies a 3x2 matrix times a 2x1 matrix
107        size_t nr_left  = 3;
108        size_t n_middle   = 2;
109        size_t nc_right = 1;
110        atomic_eigen_mat_mul<scalar> mat_mul(nr_left, n_middle, nc_right);
111        // -------------------------------------------------------------------
112        //        [ 0  0 ]
113        // left = [ 1  1 ]
114        //        [ 2  2 ]
115        ad_matrix ad_left(nr_left, n_middle);
116        for(size_t i = 0; i < nr_left; i++)
117        {       for(size_t j = 0; j < n_middle; j++)
118                        ad_left(i, j) = scalar( (j + 1) * i );
119        }
120        // -------------------------------------------------------------------
121        // declare independent variable vector x
122        size_t n = 2;
123        CPPAD_TESTVECTOR(ad_scalar) ad_x(n);
124        for(size_t j = 0; j < n; j++)
125                ad_x[j] = ad_scalar(j);
126        CppAD::Independent(ad_x);
127        // -------------------------------------------------------------------
128        // right = [ x[0] , x[1] ]^T
129        ad_matrix ad_right(n_middle, nc_right);
130        for(size_t i = 0; i < n_middle; i++)
131        {       for(size_t j = 0; j < nc_right; j++)
132                        ad_right(i, j) = ad_x[i];
133        }
134        // -------------------------------------------------------------------
135        // use atomic operation to multiply left * right
136        ad_matrix ad_result = matrix_multiply(mat_mul, ad_left, ad_right);
137        // -------------------------------------------------------------------
138        // check that first component of result is a parameter
139        // and the other components are varaibles.
140        ok &= Parameter( ad_result(0, 0) );
141        ok &= Variable(  ad_result(1, 0) );
142        ok &= Variable(  ad_result(2, 0) );
143        // -------------------------------------------------------------------
144        // declare the dependent variable vector y
145        size_t m = 3;
146        CPPAD_TESTVECTOR(ad_scalar) ad_y(m);
147        for(size_t i = 0; i < m; i++)
148                ad_y[i] = ad_result(i, 0);
149        CppAD::ADFun<scalar> f(ad_x, ad_y);
150        // -------------------------------------------------------------------
151        // check zero order forward mode
152        CPPAD_TESTVECTOR(scalar) x(n), y(m);
153        for(size_t i = 0; i < n; i++)
154                x[i] = scalar(i + 2);
155        y   = f.Forward(0, x);
156        ok &= NearEqual(y[0], 0.0,                     eps, eps);
157        ok &= NearEqual(y[1], x[0] + 2.0 * x[1],       eps, eps);
158        ok &= NearEqual(y[2], 2.0 * x[0] + 4.0 * x[1], eps, eps);
159        // -------------------------------------------------------------------
160        // check first order forward mode
161        CPPAD_TESTVECTOR(scalar) x1(n), y1(m);
162        x1[0] = 1.0;
163        x1[1] = 0.0;
164        y1    = f.Forward(1, x1);
165        ok   &= NearEqual(y1[0], 0.0, eps, eps);
166        ok   &= NearEqual(y1[1], 1.0, eps, eps);
167        ok   &= NearEqual(y1[2], 2.0, eps, eps);
168        x1[0] = 0.0;
169        x1[1] = 1.0;
170        y1    = f.Forward(1, x1);
171        ok   &= NearEqual(y1[0], 0.0, eps, eps);
172        ok   &= NearEqual(y1[1], 2.0, eps, eps);
173        ok   &= NearEqual(y1[2], 4.0, eps, eps);
174        // -------------------------------------------------------------------
175        // check first order reverse mode
176        CPPAD_TESTVECTOR(scalar) w(m), dw(n);
177        w[0]  = 0.0;
178        w[1]  = 1.0;
179        w[2]  = 0.0;
180        dw    = f.Reverse(1, w);
181        ok   &= NearEqual(dw[0], 1.0, eps, eps);
182        ok   &= NearEqual(dw[1], 2.0, eps, eps);
183        w[0]  = 0.0;
184        w[1]  = 0.0;
185        w[2]  = 1.0;
186        dw    = f.Reverse(1, w);
187        ok   &= NearEqual(dw[0], 2.0, eps, eps);
188        ok   &= NearEqual(dw[1], 4.0, eps, eps);
189        return ok;
190}
191/* %$$
192$$ $comment end nospell$$
193$end
194*/
Note: See TracBrowser for help on using the repository browser.