source: trunk/cppad/example/matrix_mul.hpp @ 2903

Last change on this file since 2903 was 2903, checked in by bradbell, 7 years ago

Add lines from atomic_base documentation to corresponding examples.

package.sh: ignore the bug/build directory.

  • Property svn:keywords set to Id
File size: 14.9 KB
Line 
1// $Id: matrix_mul.hpp 2903 2013-09-19 15:53:13Z bradbell $
2# ifndef CPPAD_MATRIX_MUL_INCLUDED
3# define CPPAD_MATRIX_MUL_INCLUDED
4
5/* --------------------------------------------------------------------------
6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 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_matrix_mul.hpp$$
18$spell
19$$
20
21$section Matrix Multiply as an Atomic Operation$$
22$index multiply, matrix atomic operation$$
23$index atomic, matrix multiply operation$$
24$index matrix, multiply atomic operation$$
25
26$nospell
27
28$head Start Class Definition$$
29$codep */
30# include <cppad/cppad.hpp>
31namespace { // Begin empty namespace
32using CppAD::vector;
33//
34void my_union(
35        std::set<size_t>&         result  ,
36        const std::set<size_t>&   left    ,
37        const std::set<size_t>&   right   )
38{       std::set<size_t> temp;
39        std::set_union(
40                left.begin()              ,
41                left.end()                ,
42                right.begin()             ,
43                right.end()               ,
44                std::inserter(temp, temp.begin())
45        );
46        result.swap(temp);
47}
48//
49// matrix result = left * right
50class matrix_mul : public CppAD::atomic_base<double> {
51/* $$
52$head Constructor$$
53$codep */
54        private:
55        // number of rows in left operand and in the result
56        const size_t nr_result_;
57        // number of columns in left operand and rows in right operand
58        const size_t n_middle_;
59        // number of columns in right operand and in the result
60        const size_t nc_result_;
61        // dimension of the domain space
62        const size_t n_;
63        // dimension of the range space
64        const size_t m_;
65        public:
66        // ---------------------------------------------------------------------
67        // constructor
68        matrix_mul(size_t nr_result, size_t n_middle, size_t nc_result)
69        : CppAD::atomic_base<double>("matrix_mul"),
70        nr_result_(nr_result) ,
71        n_middle_(n_middle)    ,
72        nc_result_(nc_result) ,
73        n_( nr_result * n_middle + n_middle * nc_result ) ,
74        m_( n_middle * nc_result ) 
75        { }
76        private:
77/* $$
78$head Left Operand Element Index$$
79$codep */
80        // left matrix element index in the taylor coefficient vector tx.
81        size_t left(
82                size_t i  , // left matrix row index
83                size_t j  , // left matrix column index
84                size_t k  , // Taylor coeffocient order
85                size_t nk ) // number of Taylor coefficients in tx
86        {       assert( i < nr_result_ );
87                assert( j < n_middle_ );
88                return (i * n_middle_ + j) * nk + k;
89        }
90/* $$
91$head Right Operand Element Index$$
92$codep */
93        // right matrix element index in the taylor coefficient vector tx.
94        size_t right(
95                size_t i  , // right matrix row index
96                size_t j  , // right matrix column index
97                size_t k  , // Taylor coeffocient order
98                size_t nk ) // number of Taylor coefficients in tx
99        {       assert( i < n_middle_  );
100                assert( j < nc_result_ );
101                size_t offset = nr_result_ * n_middle_;
102                return (offset + i * nc_result_ + j) * nk + k; 
103        }
104/* $$
105$head Result Element Index$$
106$codep */
107        // result matrix element index in the taylor coefficient vector ty.
108        size_t result(
109                size_t i  , // result matrix row index
110                size_t j  , // result matrix column index
111                size_t k  , // Taylor coeffocient order
112                size_t nk ) // number of Taylor coefficients in ty
113        {       assert( i < nr_result_  );
114                assert( j < nc_result_ );
115                return (i * nc_result_ + j) * nk + k;
116        }
117/* $$
118$head Forward Matrix Multipliy$$
119$codep */
120        // Forward mode multiply Taylor coefficients in tx and sum into ty
121        // (for one pair of left and right orders)
122        void forward_multiply(
123                size_t                 k_left  , // order for left coefficients
124                size_t                 k_right , // order for right coefficients
125                const vector<double>&  tx      , // domain space Taylor coefficients
126                      vector<double>&  ty      ) // range space Taylor coefficients
127        {       size_t nk       = tx.size() / n_;
128                assert( nk == ty.size() / m_ );
129                //
130                size_t k_result = k_left + k_right;
131                assert( k_result < nk );
132                //
133                for(size_t i = 0; i < nr_result_; i++)
134                {       for(size_t j = 0; j < nc_result_; j++)
135                        {       double sum = 0.0;
136                                for(size_t ell = 0; ell < n_middle_; ell++)
137                                {       size_t i_left  = left(i, ell, k_left, nk);
138                                        size_t i_right = right(ell, j,  k_right, nk);
139                                        sum           += tx[i_left] * tx[i_right];
140                                }
141                                size_t i_result = result(i, j, k_result, nk);
142                                ty[i_result]   += sum;
143                        }
144                }
145        }
146/* $$
147$head Reverse Matrix Multipliy$$
148$codep */
149        // Reverse mode partials of Taylor coefficients and sum into px
150        // (for one pair of left and right orders)
151        void reverse_multiply(
152                size_t                 k_left  , // order for left coefficients
153                size_t                 k_right , // order for right coefficients
154                const vector<double>&  tx      , // domain space Taylor coefficients
155                const vector<double>&  ty      , // range space Taylor coefficients
156                      vector<double>&  px      , // partials w.r.t. tx
157                const vector<double>&  py      ) // partials w.r.t. ty
158        {       size_t nk       = tx.size() / n_;
159                assert( nk == ty.size() / m_ );
160                assert( tx.size() == px.size() );
161                assert( ty.size() == py.size() );
162                //
163                size_t k_result = k_left + k_right;
164                assert( k_result < nk );
165                //
166                for(size_t i = 0; i < nr_result_; i++)
167                {       for(size_t j = 0; j < nc_result_; j++)
168                        {       size_t i_result = result(i, j, k_result, nk);
169                                for(size_t ell = 0; ell < n_middle_; ell++)
170                                {       size_t i_left  = left(i, ell, k_left, nk);
171                                        size_t i_right = right(ell, j,  k_right, nk);
172                                        // sum        += tx[i_left] * tx[i_right];
173                                        px[i_left]    += tx[i_right] * py[i_result];
174                                        px[i_right]   += tx[i_left]  * py[i_result];
175                                }
176                        }
177                }
178                return;
179        }
180/* $$
181$head forward$$
182$codep */
183        // forward mode routine called by CppAD
184        bool forward(
185                size_t                    q ,
186                size_t                    p ,
187                const vector<bool>&      vx ,
188                      vector<bool>&      vy ,
189                const vector<double>&    tx ,
190                      vector<double>&    ty
191        )
192        {       size_t p1 = p + 1;
193                assert( vx.size() == 0 || n_ == vx.size() );
194                assert( vx.size() == 0 || m_ == vy.size() );
195                assert( n_ * p1 == tx.size() );
196                assert( m_ * p1 == ty.size() );
197                size_t i, j, ell;
198
199                // check if we are computing vy information
200                if( vx.size() > 0 )
201                {       size_t nk = 1;
202                        size_t k  = 0;
203                        for(i = 0; i < nr_result_; i++)
204                        {       for(j = 0; j < nc_result_; j++)
205                                {       bool var = false;
206                                        for(ell = 0; ell < n_middle_; ell++)
207                                        {       size_t i_left  = left(i, ell, k, nk);
208                                                size_t i_right = right(ell, j, k, nk);
209                                                bool   nz_left = vx[i_left] |(tx[i_left]  != 0.);
210                                                bool  nz_right = vx[i_right]|(tx[i_right] != 0.);
211                                                // if not multiplying by the constant zero
212                                                if( nz_left & nz_right )
213                                                                var |= vx[i_left] | vx[i_right];
214                                        }
215                                        size_t i_result = result(i, j, k, nk);
216                                        vy[i_result] = var;
217                                }
218                        }
219                }
220
221                // initialize result as zero
222                size_t k;
223                for(i = 0; i < nr_result_; i++)
224                {       for(j = 0; j < nc_result_; j++)
225                        {       for(k = q; k <= p; k++)
226                                        ty[ result(i, j, k, p1) ] = 0.0;
227                        }
228                }
229                for(k = q; k <= p; k++)
230                {       // sum the produces that result in order k
231                        for(ell = 0; ell <= k; ell++)
232                                forward_multiply(ell, k - ell, tx, ty);
233                }
234
235                // all orders are implented, so always return true
236                return true;
237        }
238/* $$
239$head reverse$$
240$codep */
241        // reverse mode routine called by CppAD
242        virtual bool reverse(
243                size_t                     p ,
244                const vector<double>&     tx ,
245                const vector<double>&     ty ,
246                      vector<double>&     px ,
247                const vector<double>&     py
248        )
249        {       size_t p1 = p + 1;
250                assert( n_ * p1 == tx.size() );
251                assert( m_ * p1 == ty.size() );
252                assert( px.size() == tx.size() );
253                assert( py.size() == ty.size() );
254
255                // initialize summation
256                for(size_t i = 0; i < px.size(); i++)
257                        px[i] = 0.0;
258
259                // number of orders to differentiate
260                size_t k = p1;
261                while(k--)
262                {       // differentiate the produces that result in order k
263                        for(size_t ell = 0; ell <= k; ell++)
264                                reverse_multiply(ell, k - ell, tx, ty, px, py);
265                }
266
267                // all orders are implented, so always return true
268                return true;
269        }
270/* $$
271$head for_sparse_jac$$
272$codep */
273        // forward Jacobian sparsity routine called by CppAD
274        virtual bool for_sparse_jac(
275                size_t                                q ,
276                const vector<bool>&                   r ,
277                      vector<bool>&                   s )
278        {       assert( n_ * q == r.size() );
279                assert( m_ * q == s.size() );
280                size_t p;
281
282                // sparsity for S(x) = f'(x) * R
283                size_t nk = 1;
284                size_t k  = 0;
285                for(size_t i = 0; i < nr_result_; i++)
286                {       for(size_t j = 0; j < nc_result_; j++)
287                        {       size_t i_result = result(i, j, k, nk);
288                                for(p = 0; p < q; p++)
289                                        s[i_result * q + p] = false; 
290                                for(size_t ell = 0; ell < n_middle_; ell++)
291                                {       size_t i_left  = left(i, ell, k, nk);
292                                        size_t i_right = right(ell, j, k, nk);
293                                        for(p = 0; p < q; p++)
294                                        {       s[i_result * q + p] |= r[i_left * q + p ]; 
295                                                s[i_result * q + p] |= r[i_right * q + p ]; 
296                                        }
297                                }
298                        }
299                }
300                return true;
301        }
302        virtual bool for_sparse_jac(
303                size_t                                q ,
304                const vector< std::set<size_t> >&     r ,
305                      vector< std::set<size_t> >&     s )
306        {       assert( n_ == r.size() );
307                assert( m_ == s.size() );
308
309                // sparsity for S(x) = f'(x) * R
310                size_t nk = 1;
311                size_t k  = 0;
312                for(size_t i = 0; i < nr_result_; i++)
313                {       for(size_t j = 0; j < nc_result_; j++)
314                        {       size_t i_result = result(i, j, k, nk);
315                                s[i_result].clear();
316                                for(size_t ell = 0; ell < n_middle_; ell++)
317                                {       size_t i_left  = left(i, ell, k, nk);
318                                        size_t i_right = right(ell, j, k, nk);
319                                        //
320                                        my_union( s[i_result], s[i_result], r[i_left] );
321                                        my_union( s[i_result], s[i_result], r[i_right] );
322                                }
323                        }
324                }
325                return true;
326        }
327/* $$
328$head rev_sparse_jac$$
329$codep */
330        // reverse Jacobian sparsity routine called by CppAD
331        virtual bool rev_sparse_jac(
332                size_t                                q ,
333                const vector<bool>&                  rt ,
334                      vector<bool>&                  st )
335        {       assert( n_ * q == st.size() );
336                assert( m_ * q == rt.size() );
337                size_t i, j, p;
338
339                // initialize
340                for(i = 0; i < n_; i++)
341                {       for(p = 0; p < q; p++)
342                                st[ i * q + p ] = false;
343                }
344
345                // sparsity for S(x)^T = f'(x)^T * R^T
346                size_t nk = 1;
347                size_t k  = 0;
348                for(i = 0; i < nr_result_; i++)
349                {       for(j = 0; j < nc_result_; j++)
350                        {       size_t i_result = result(i, j, k, nk);
351                                for(size_t ell = 0; ell < n_middle_; ell++)
352                                {       size_t i_left  = left(i, ell, k, nk);
353                                        size_t i_right = right(ell, j, k, nk);
354                                        for(p = 0; p < q; p++)
355                                        {       st[i_left * q + p] |= rt[i_result * q + p];
356                                                st[i_right* q + p] |= rt[i_result * q + p];
357                                        }
358                                }
359                        }
360                }
361                return true;
362        }
363        virtual bool rev_sparse_jac(
364                size_t                                q ,
365                const vector< std::set<size_t> >&    rt ,
366                      vector< std::set<size_t> >&    st )
367        {       assert( n_ == st.size() );
368                assert( m_ == rt.size() );
369                size_t i, j;
370
371                // initialize
372                for(i = 0; i < n_; i++)
373                        st[i].clear();
374
375                // sparsity for S(x)^T = f'(x)^T * R^T
376                size_t nk = 1;
377                size_t k  = 0;
378                for(i = 0; i < nr_result_; i++)
379                {       for(j = 0; j < nc_result_; j++)
380                        {       size_t i_result = result(i, j, k, nk);
381                                for(size_t ell = 0; ell < n_middle_; ell++)
382                                {       size_t i_left  = left(i, ell, k, nk);
383                                        size_t i_right = right(ell, j, k, nk);
384                                        //
385                                        my_union(st[i_left],  st[i_left],  rt[i_result]);
386                                        my_union(st[i_right], st[i_right], rt[i_result]);
387                                }
388                        }
389                }
390                return true;
391        }
392/* $$
393$head rev_sparse_hes$$
394$codep */
395        // reverse Hessian sparsity routine called by CppAD
396        virtual bool rev_sparse_hes(
397                const vector<bool>&                   vx,
398                const vector<bool>&                   s ,
399                      vector<bool>&                   t ,
400                size_t                                q ,
401                const vector< std::set<size_t> >&     r ,
402                const vector< std::set<size_t> >&     u ,
403                      vector< std::set<size_t> >&     v )
404        {       size_t n = vx.size();   
405                assert( t.size() == n );
406                assert( r.size() == n );
407                assert( v.size() == n );
408# ifndef NDEBUG
409                size_t m = s.size();
410                assert( u.size() == m );
411# endif
412                size_t i, j;
413                //
414                // initilaize sparsity patterns as false
415                for(j = 0; j < n; j++)
416                {       t[j] = false;
417                        v[j].clear();
418                }
419                size_t nk = 1;
420                size_t k  = 0;
421                for(i = 0; i < nr_result_; i++)
422                {       for(j = 0; j < nc_result_; j++)
423                        {       size_t i_result = result(i, j, k, nk);
424                                for(size_t ell = 0; ell < n_middle_; ell++)
425                                {       size_t i_left  = left(i, ell, k, nk);
426                                        size_t i_right = right(ell, j, k, nk);
427                                        //
428                                        // Compute sparsity for T(x) = S(x) * f'(x).
429                                        // We need not use vx with f'(x) back propagation.
430                                        t[i_left]  |= s[i_result];
431                                        t[i_right] |= s[i_result];
432
433                                        // V(x) = f'(x)^T * U(x) +  S(x) * f''(x) * R
434                                        // U(x) = g''(y) * f'(x) * R
435                                        // S(x) = g'(y)
436
437                                        // back propagate f'(x)^T * U(x)
438                                        // (no need to use vx with f'(x) propogation)
439                                        my_union(v[i_left],  v[i_left],  u[i_result] );
440                                        my_union(v[i_right], v[i_right], u[i_result] );
441
442                                        // back propagate S(x) * f''(x) * R
443                                        // (here is where we must check for cross terms)
444                                        if( s[i_result] & vx[i_left] & vx[i_right] )
445                                        {       my_union(v[i_left],  v[i_left],  r[i_right] );
446                                                my_union(v[i_right], v[i_right], r[i_left]  );
447                                        }
448                                }
449                        }
450                }
451                return true;
452        }
453        virtual bool rev_sparse_hes(
454                const vector<bool>&                   vx,
455                const vector<bool>&                   s ,
456                      vector<bool>&                   t ,
457                size_t                                q ,
458                const vector<bool>&                   r ,
459                const vector<bool>&                   u ,
460                      vector<bool>&                   v )
461        {       size_t n = vx.size();
462                assert( t.size() == n );
463                assert( r.size() == n * q );
464                assert( v.size() == n * q );
465# ifndef NDEBUG
466                size_t m = s.size();
467                assert( u.size() == m * q );
468# endif
469                size_t i, j, p;
470                //
471                // initilaize sparsity patterns as false
472                for(j = 0; j < n; j++)
473                {       t[j] = false;
474                        for(p = 0; p < q; p++)
475                                v[j * q + p] = false;
476                }
477                size_t nk = 1;
478                size_t k  = 0;
479                for(i = 0; i < nr_result_; i++)
480                {       for(j = 0; j < nc_result_; j++)
481                        {       size_t i_result = result(i, j, k, nk);
482                                for(size_t ell = 0; ell < n_middle_; ell++)
483                                {       size_t i_left  = left(i, ell, k, nk);
484                                        size_t i_right = right(ell, j, k, nk);
485                                        //
486                                        // Compute sparsity for T(x) = S(x) * f'(x).
487                                        // We so not need to use vx with f'(x) propagation.
488                                        t[i_left]  |= s[i_result];
489                                        t[i_right] |= s[i_result];
490
491                                        // V(x) = f'(x)^T * U(x) +  S(x) * f''(x) * R
492                                        // U(x) = g''(y) * f'(x) * R
493                                        // S(x) = g'(y)
494
495                                        // back propagate f'(x)^T * U(x)
496                                        // (no need to use vx with f'(x) propogation)
497                                        for(p = 0; p < q; p++)
498                                        {       v[ i_left  * q + p] |= u[ i_result * q + p];
499                                                v[ i_right * q + p] |= u[ i_result * q + p];
500                                        }
501
502                                        // back propagate S(x) * f''(x) * R
503                                        // (here is where we must check for cross terms)
504                                        if( s[i_result] & vx[i_left] & vx[i_right] )
505                                        {       for(p = 0; p < q; p++)
506                                                {       v[i_left * q + p]  |= r[i_right * q + p];
507                                                        v[i_right * q + p] |= r[i_left * q + p];
508                                                }
509                                        }
510                                }
511                        }
512                }
513                return true;
514        }
515
516/* $$
517$head End Class Definition$$
518$codep */
519}; // End of matrix_mul class
520}  // End empty namespace
521/* $$
522$$ $comment end nospell$$
523$end
524*/
525
526
527# endif
Note: See TracBrowser for help on using the repository browser.