source: trunk/cppad/example/eigen_mat_mul.hpp @ 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: 11.3 KB
Line 
1// $Id$
2# ifndef CPPAD_EXAMPLE_EIGEN_MAT_MUL_HPP
3# define CPPAD_EXAMPLE_EIGEN_MAT_MUL_HPP
4
5/* --------------------------------------------------------------------------
6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-16 Bradley M. Bell
7
8CppAD is distributed under multiple licenses. This distribution is under
9the terms of the
10                    Eclipse Public License Version 1.0.
11
12A copy of this license is included in the COPYING file of this distribution.
13Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
14-------------------------------------------------------------------------- */
15
16/*
17$begin atomic_eigen_mat_mul.hpp$$
18$spell
19        Eigen
20$$
21
22$section Atomic Eigen Matrix Multiply Class$$
23
24$nospell
25
26$head Start Class Definition$$
27$srccode%cpp% */
28# include <cppad/cppad.hpp>
29# include <Eigen/Core>
30
31
32/* %$$
33$head Publice Types$$
34$srccode%cpp% */
35namespace { // BEGIN_EMPTY_NAMESPACE
36
37template <class Base>
38class atomic_eigen_mat_mul : public CppAD::atomic_base<Base> {
39public:
40        // -----------------------------------------------------------
41        // type of elements during calculation of derivatives
42        typedef Base              scalar;
43        // type of elements during taping
44        typedef CppAD::AD<scalar> ad_scalar;
45        // type of matrix during calculation of derivatives
46        typedef Eigen::Matrix<
47                scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>     matrix;
48        // type of matrix during taping
49        typedef Eigen::Matrix<
50                ad_scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > ad_matrix;
51/* %$$
52$head Public Constructor$$
53$srccode%cpp% */
54        // constructor
55        atomic_eigen_mat_mul(
56                // number of rows in left operand
57                const size_t nr_left  ,
58                // number of rows in left and columns in right operand
59                const size_t n_middle   ,
60                // number of columns in right operand
61                const size_t nc_right
62        ) :
63        CppAD::atomic_base<Base>(
64                "atom_eigen_mat_mul"                             ,
65                CppAD::atomic_base<Base>::set_sparsity_enum
66        ) ,
67        nr_left_(  nr_left  )                   ,
68        n_middle_(   n_middle   )                   ,
69        nc_right_( nc_right )                   ,
70        nx_( (nr_left + nc_right) * n_middle_ ) ,
71        ny_( nr_left * nc_right )
72        { }
73/* %$$
74$head Public Pack$$
75$srccode%cpp% */
76        template <class Matrix, class Vector>
77        void pack(
78                Vector&        packed  ,
79                const Matrix&  left    ,
80                const Matrix&  right   )
81        {       assert( packed.size() == nx_      );
82                assert( rows( left ) == nr_left_ );
83                assert( cols( left ) == n_middle_ );
84                assert( rows( right ) == n_middle_ );
85                assert( cols( right ) == nc_right_ );
86                //
87                size_t n_left = nr_left_ * n_middle_;
88                size_t n_right = n_middle_ * nc_right_;
89                assert( n_left + n_right == nx_ );
90                //
91                for(size_t i = 0; i < n_left; i++)
92                        packed[i] = left.data()[i];
93                for(size_t i = 0; i < n_right; i++)
94                        packed[ i + n_left ] = right.data()[i];
95                //
96                return;
97        }
98/* %$$
99$head Public Unpack$$
100$srccode%cpp% */
101        template <class Matrix, class Vector>
102        void unpack(
103                const Vector&   packed  ,
104                Matrix&         result  )
105        {       assert( packed.size() == ny_      );
106                assert( rows( result ) == nr_left_ );
107                assert( cols( result ) == nc_right_ );
108                //
109                size_t n_result = nr_left_ * nc_right_;
110                assert( n_result == ny_ );
111                //
112                for(size_t i = 0; i < n_result; i++)
113                        result.data()[i] = packed[ i ];
114                //
115                return;
116        }
117/* %$$
118$head Private Variables$$
119$srccode%cpp% */
120private:
121        // -------------------------------------------------------------
122        // number of of rows in left operand and result
123        const size_t nr_left_;
124        // number of of columns in left operand and rows in right operand
125        const size_t n_middle_;
126        // number of columns in right operand and result
127        const size_t nc_right_;
128        // size of the domain space
129        const size_t nx_;
130        // size of the range space
131        const size_t ny_;
132        // -------------------------------------------------------------
133        // one forward mode vector of matrices for left, right, and result
134        CppAD::vector<matrix> f_left_, f_right_, f_result_;
135        // one reverse mode vector of matrices for left, right, and result
136        CppAD::vector<matrix> r_left_, r_right_, r_result_;
137        // -------------------------------------------------------------
138/* %$$
139$head Private rows$$
140$srccode%cpp% */
141        // convert from int to size_t
142        static size_t rows(const matrix& x)
143        {       return size_t( x.rows() ); }
144        static size_t rows(const ad_matrix& x)
145        {       return size_t( x.rows() ); }
146/* %$$
147$head Private cols$$
148$srccode%cpp% */
149        // convert from int to size_t
150        static size_t cols(const matrix& x)
151        {       return size_t( x.cols() ); }
152        static size_t cols(const ad_matrix& x)
153        {       return size_t( x.cols() ); }
154/* %$$
155$head Private forward$$
156$srccode%cpp% */
157        // forward mode routine called by CppAD
158        virtual bool forward(
159                // lowest order Taylor coefficient we are evaluating
160                size_t                          p ,
161                // highest order Taylor coefficient we are evaluating
162                size_t                          q ,
163                // which components of x are variables
164                const CppAD::vector<bool>&      vx ,
165                // which components of y are variables
166                CppAD::vector<bool>&            vy ,
167                // tx [ j * (q+1) + k ] is x_j^k
168                const CppAD::vector<scalar>&    tx ,
169                // ty [ i * (q+1) + k ] is y_i^k
170                CppAD::vector<scalar>&          ty
171        )
172        {       size_t n_order = q + 1;
173                assert( vx.size() == 0 || nx_ == vx.size() );
174                assert( vx.size() == 0 || ny_ == vy.size() );
175                assert( nx_ * n_order == tx.size() );
176                assert( ny_ * n_order == ty.size() );
177                //
178                size_t n_left   = nr_left_ * n_middle_;
179                size_t n_right  = n_middle_ * nc_right_;
180                size_t n_result = nr_left_ * nc_right_;
181                assert( n_left + n_right == nx_ );
182                assert( n_result == ny_ );
183                //
184                // -------------------------------------------------------------------
185                // make sure f_left_, f_right_, and f_result_ are large enough
186                assert( f_left_.size() == f_right_.size() );
187                assert( f_left_.size() == f_result_.size() );
188                if( f_left_.size() < n_order )
189                {       f_left_.resize(n_order);
190                        f_right_.resize(n_order);
191                        f_result_.resize(n_order);
192                        //
193                        for(size_t k = 0; k < n_order; k++)
194                        {       f_left_[k].resize(nr_left_, n_middle_);
195                                f_right_[k].resize(n_middle_, nc_right_);
196                                f_result_[k].resize(nr_left_, nc_right_);
197                        }
198                }
199                // -------------------------------------------------------------------
200                // unpack tx into f_left and f_right
201                for(size_t k = 0; k < n_order; k++)
202                {       // unpack left values for this order
203                        for(size_t i = 0; i < n_left; i++)
204                                f_left_[k].data()[i] = tx[ i * n_order + k ];
205                        //
206                        // unpack right values for this order
207                        for(size_t i = 0; i < n_right; i++)
208                                f_right_[k].data()[i] = tx[ (i + n_left) * n_order + k ];
209                }
210                // -------------------------------------------------------------------
211                // result for each order
212                for(size_t k = 0; k < n_order; k++)
213                {       // result[k] = sum_ell left[ell] * right[k-ell]
214                        f_result_[k] = matrix::Zero(nr_left_, nc_right_);
215                        for(size_t ell = 0; ell <= k; ell++)
216                                f_result_[k] += f_left_[ell] * f_right_[k-ell];
217                }
218                // -------------------------------------------------------------------
219                // pack result_ into ty
220                for(size_t k = 0; k < n_order; k++)
221                {       for(size_t i = 0; i < n_result; i++)
222                                ty[ i * n_order + k ] = f_result_[k].data()[i];
223                }
224
225                // check if we are compute vy
226                if( vx.size() == 0 )
227                        return true;
228                // ------------------------------------------------------------------
229                // compute variable information for y; i.e., vy
230                // (note that the constant zero times a variable is a constant)
231                scalar zero(0.0);
232                assert( n_order == 1 );
233                for(size_t i = 0; i < nr_left_; i++)
234                {       for(size_t j = 0; j < nc_right_; j++)
235                        {       bool var = false;
236                                for(size_t ell = 0; ell < n_middle_; ell++)
237                                {       size_t index   = i * n_middle_ + ell;
238                                        bool var_left  = vx[index];
239                                        bool nz_left   = var_left | (f_left_[0](i, ell) != zero);
240                                        index          = nr_left_ * n_middle_;
241                                        index         += ell * nc_right_ + j;
242                                        bool var_right = vx[index];
243                                        bool nz_right  = var_right | (f_right_[0](ell, j) != zero);
244                                        var |= var_left & nz_right;
245                                        var |= nz_left  & var_right;
246                                }
247                                size_t index = i * nc_right_ + j;
248                                vy[index]    = var;
249                        }
250                }
251                return true;
252        }
253/* %$$
254$head Private reverse$$
255$srccode%cpp% */
256        // reverse mode routine called by CppAD
257        virtual bool reverse(
258                // highest order Taylor coefficient that we are computing deritive of
259                size_t                     q ,
260                // forward mode Taylor coefficients for x variables
261                const CppAD::vector<double>&     tx ,
262                // forward mode Taylor coefficients for y variables
263                const CppAD::vector<double>&     ty ,
264                // upon return, derivative of G[ F[ {x_j^k} ] ] w.r.t {x_j^k}
265                CppAD::vector<double>&           px ,
266                // derivative of G[ {y_i^k} ] w.r.t. {y_i^k}
267                const CppAD::vector<double>&     py
268        )
269        {       size_t n_order = q + 1;
270                assert( nx_ * n_order == tx.size() );
271                assert( ny_ * n_order == ty.size() );
272                assert( px.size() == tx.size() );
273                assert( py.size() == ty.size() );
274                //
275                size_t n_left   = nr_left_ * n_middle_;
276                size_t n_right  = n_middle_ * nc_right_;
277                size_t n_result = nr_left_ * nc_right_;
278                assert( n_left + n_right == nx_ );
279                assert( n_result == ny_ );
280                //
281                // -------------------------------------------------------------------
282                // make sure r_left_, r_right_, and r_result_ are large enough
283                assert( r_left_.size() == r_right_.size() );
284                assert( r_left_.size() == r_result_.size() );
285                if( r_left_.size() < n_order )
286                {       r_left_.resize(n_order);
287                        r_right_.resize(n_order);
288                        r_result_.resize(n_order);
289                        //
290                        for(size_t k = 0; k < n_order; k++)
291                        {       r_left_[k].resize(nr_left_, n_middle_);
292                                r_right_[k].resize(n_middle_, nc_right_);
293                                r_result_[k].resize(nr_left_, nc_right_);
294                        }
295                }
296                // -------------------------------------------------------------------
297                // unpack tx into f_left and f_right
298                for(size_t k = 0; k < n_order; k++)
299                {       // unpack left values for this order
300                        for(size_t i = 0; i < n_left; i++)
301                                f_left_[k].data()[i] = tx[ i * n_order + k ];
302                        //
303                        // unpack right values for this order
304                        for(size_t i = 0; i < n_right; i++)
305                                f_right_[k].data()[i] = tx[ (i + n_left) * n_order + k ];
306                }
307                // -------------------------------------------------------------------
308                // unpack result_ from py
309                for(size_t k = 0; k < n_order; k++)
310                {       for(size_t i = 0; i < n_result; i++)
311                                r_result_[k].data()[i] = py[ i * n_order + k ];
312                }
313                // -------------------------------------------------------------------
314                // initialize r_left_ and r_right_ as zero
315                for(size_t k = 0; k < n_order; k++)
316                {       r_left_[k]   = matrix::Zero(nr_left_, n_middle_);
317                        r_right_[k]  = matrix::Zero(n_middle_, nc_right_);
318                }
319                // -------------------------------------------------------------------
320                // matrix reverse mode calculation
321                for(size_t k1 = n_order; k1 > 0; k1--)
322                {       size_t k = k1 - 1;
323                        for(size_t ell = 0; ell <= k; ell++)
324                        {       // nr x nm       = nr x nc      * nc * nm
325                                r_left_[ell]    += r_result_[k] * f_right_[k-ell].transpose();
326                                // nm x nc       = nm x nr * nr * nc
327                                r_right_[k-ell] += f_left_[ell].transpose() * r_result_[k];
328                        }
329                }
330                // -------------------------------------------------------------------
331                // pack r_left and r_right int px
332                for(size_t k = 0; k < n_order; k++)
333                {       // pack left values for this order
334                        for(size_t i = 0; i < n_left; i++)
335                                px[ i * n_order + k ] = r_left_[k].data()[i];
336                        //
337                        // pack right values for this order
338                        for(size_t i = 0; i < n_right; i++)
339                                px[ (i + n_left) * n_order + k] = r_right_[k].data()[i];
340                }
341                //
342                return true;
343        }
344/* %$$
345$head End Class Definition$$
346$srccode%cpp% */
347}; // End of atomic_eigen_mat_mul class
348
349}  // END_EMPTY_NAMESPACE
350/* %$$
351$$ $comment end nospell$$
352$end
353*/
354
355
356# endif
Note: See TracBrowser for help on using the repository browser.