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

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