source: trunk/example/mat_mul.hpp @ 2794

Last change on this file since 2794 was 2794, checked in by bradbell, 7 years ago
  1. Use CPPAD_NULL, intead of 0, for null pointer.

check_if_0.sh: Ignore subdirectories of new directories.
jenkins.sh: output logs when an error occurs.
acos_op.hpp: avoid use of q (will use it for an order index).
asin_op.hpp: avoid use of q (will use it for an order index).
forward_sweep.hpp: chnage d to p, use const in prototype.
div_op.hpp: minor edit of white space.
atom_usead_2.cpp: use ADFUN to compute variable/parameter information.

  • Property svn:keywords set to Id
File size: 12.6 KB
Line 
1// $Id: mat_mul.hpp 2794 2013-05-02 08:20:30Z bradbell $
2# ifndef CPPAD_MAT_MUL_INCLUDED
3# define CPPAD_MAT_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 mat_mul.hpp$$
18$spell
19        mat_mul.hpp
20        cppad
21        CppAD
22        namespace
23        struct
24        nr
25        nc
26        bool
27        vx
28        const
29        im
30        mj
31        ij
32        px
33        py
34        std
35        tx
36        ty
37        resize
38        nz
39        var
40        jac
41        Jacobian
42        hes
43        vy
44$$
45
46
47$section Define Matrix Multiply as a User Atomic Operation$$
48
49$index mat_mul, define$$
50$index matrix, multiply$$
51$index multiply, matrix$$
52$index user_atomic, test$$
53$index test, user_atomic$$
54$index user_atomic, example$$
55$index example, user_atomic$$
56
57$head Syntax$$
58This file is located in the $code example$$ directory.
59It can be copied to the current working directory and included
60with the syntax
61$codei%
62        # include "mat_mul.hpp"
63%$$
64
65$head Example$$
66The file $cref mat_mul.cpp$$ contains an example use of
67$code mat_mul.hpp$$.
68It returns true if it succeeds and false otherwise.
69
70$head Begin Source$$
71$codep */
72# include <cppad/cppad.hpp>      // Include CppAD definitions
73namespace {                      // Begin empty namespace
74        using CppAD::vector;        // Let vector denote CppAD::vector
75/* $$
76
77$head Extra Call Information$$
78$codep */
79        // Information we will attach to each mat_mul call
80        struct call_info {
81                size_t nr_result;
82                size_t n_middle;
83                size_t nc_result;
84                vector<bool>  vx;
85        };
86        vector<call_info> info_; // vector of call information
87
88        // number of orders for this operation (k + 1)
89        size_t n_order_ = 0;
90        // number of rows in the result matrix
91        size_t nr_result_ = 0;
92        // number of columns in left matrix and number of rows in right matrix
93        size_t n_middle_ = 0;
94        // number of columns in the result matrix
95        size_t nc_result_ = 0;
96        // which components of x are variables
97        vector<bool>* vx_ = CPPAD_NULL;
98
99        // get the information corresponding to this call
100        void get_info(size_t id, size_t k, size_t n, size_t m)
101        {       n_order_   = k + 1;     
102                nr_result_ = info_[id].nr_result; 
103                n_middle_  = info_[id].n_middle;
104                nc_result_ = info_[id].nc_result;
105                vx_        = &(info_[id].vx);
106
107                assert(n == nr_result_ * n_middle_ + n_middle_ * nc_result_);
108                assert(m ==  nr_result_ * nc_result_);
109        }
110
111/* $$
112$head Matrix Indexing$$
113$codep */ 
114        // Convert left matrix index pair and order to a single argument index
115        size_t left(size_t i, size_t j, size_t ell)
116        {       assert( i < nr_result_ );
117                assert( j < n_middle_ );
118                return (i * n_middle_ + j) * n_order_ + ell;
119        }
120        // Convert right matrix index pair and order to a single argument index
121        size_t right(size_t i, size_t j, size_t ell)
122        {       assert( i < n_middle_ );
123                assert( j < nc_result_ );
124                size_t offset = nr_result_ * n_middle_;
125                return (offset + i * nc_result_ + j) * n_order_ + ell;
126        }
127        // Convert result matrix index pair and order to a single result index
128        size_t result(size_t i, size_t j, size_t ell)
129        {       assert( i < nr_result_ );
130                assert( j < nc_result_ );
131                return (i * nc_result_ + j) * n_order_ + ell;
132        }
133/* $$
134
135$head One Matrix Multiply$$
136Forward mode matrix multiply left times right and sum into result:
137$codep */
138        void multiply_and_sum(
139                size_t                order_left , 
140                size_t                order_right, 
141                const vector<double>&         tx ,
142                vector<double>&               ty ) 
143        {       size_t i, j;
144                size_t order_result = order_left + order_right; 
145                for(i = 0; i < nr_result_; i++)
146                {       for(j = 0; j < nc_result_; j++)
147                        {       double sum = 0.;
148                                size_t middle, im_left, mj_right, ij_result;
149                                for(middle = 0; middle < n_middle_; middle++)
150                                {       im_left  = left(i, middle, order_left);
151                                        mj_right = right(middle, j, order_right);
152                                        sum     += tx[im_left] * tx[mj_right];
153                                }
154                                ij_result = result(i, j, order_result);
155                                ty[ ij_result ] += sum;
156                        }
157                }
158                return;
159        }
160/* $$
161
162$head Reverse Partials One Order$$
163Compute reverse mode partials for one order and sum into px:
164$codep */
165        void reverse_multiply(
166                size_t                order_left , 
167                size_t                order_right, 
168                const vector<double>&         tx ,
169                const vector<double>&         ty ,
170                vector<double>&               px ,
171                const vector<double>          py ) 
172        {       size_t i, j;
173                size_t order_result = order_left + order_right; 
174                for(i = 0; i < nr_result_; i++)
175                {       for(j = 0; j < nc_result_; j++)
176                        {       size_t middle, im_left, mj_right, ij_result;
177                                for(middle = 0; middle < n_middle_; middle++)
178                                {       ij_result = result(i, j, order_result);
179                                        im_left   = left(i, middle, order_left);
180                                        mj_right  = right(middle, j, order_right);
181                                        // sum       += tx[im_left]  * tx[mj_right];
182                                        px[im_left]  += tx[mj_right] * py[ij_result];
183                                        px[mj_right] += tx[im_left]  * py[ij_result];
184                                }
185                        }
186                }
187                return;
188        }
189/* $$
190$head Set Union$$
191$codep */
192        void my_union(
193                std::set<size_t>&         result  ,
194                const std::set<size_t>&   left    ,
195                const std::set<size_t>&   right   )
196        {       std::set<size_t> temp;
197                std::set_union(
198                        left.begin()              ,
199                        left.end()                ,
200                        right.begin()             ,
201                        right.end()               ,
202                        std::inserter(temp, temp.begin())
203                );
204                result.swap(temp);
205        }
206/* $$
207
208$head CppAD User Atomic Callback Functions$$
209$codep */
210        // ----------------------------------------------------------------------
211        // forward mode routine called by CppAD
212        bool mat_mul_forward(
213                size_t                   id ,
214                size_t                    k ,
215                size_t                    n ,
216                size_t                    m ,
217                const vector<bool>&      vx ,
218                vector<bool>&            vy ,
219                const vector<double>&    tx ,
220                vector<double>&          ty
221        )
222        {       size_t i, j, ell;
223                get_info(id, k, n, m);
224
225                // check if this is during the call to mat_mul(id, ax, ay)
226                if( vx.size() > 0 )
227                {       assert( k == 0 && vx.size() > 0 );
228
229                        // store the vx information in info_
230                        assert( vx_->size() == 0 );
231                        info_[id].vx.resize(n);
232                        for(j = 0; j < n; j++)
233                                info_[id].vx[j] = vx[j];
234                        assert( vx_->size() == n );
235                       
236                        // now compute vy
237                        for(i = 0; i < nr_result_; i++)
238                        {       for(j = 0; j < nc_result_; j++)
239                                {       // compute vy[ result(i, j, 0) ]
240                                        bool   var = false;
241                                        bool   nz_left, nz_right;
242                                        size_t middle, im_left, mj_right, ij_result;
243                                        for(middle = 0; middle < n_middle_; middle++)
244                                        {       im_left  = left(i, middle, k);
245                                                mj_right = right(middle, j, k);
246                                                nz_left  = vx[im_left]  | (tx[im_left] != 0.);
247                                                nz_right = vx[mj_right] | (tx[mj_right]!= 0.);
248                                                // if not multiplying by the constant zero
249                                                if( nz_left & nz_right )
250                                                        var |= (vx[im_left] | vx[mj_right]);
251                                        }
252                                        ij_result     = result(i, j, k);
253                                        vy[ij_result] = var;
254                                }
255                        }
256                }
257
258                // initialize result as zero
259                for(i = 0; i < nr_result_; i++)
260                {       for(j = 0; j < nc_result_; j++)
261                                ty[ result(i, j, k) ] = 0.;
262                }
263                // sum the product of proper orders
264                for(ell = 0; ell <=k; ell++)
265                        multiply_and_sum(ell, k-ell, tx, ty);
266
267                // All orders are implemented and there are no possible error
268                // conditions, so always return true.
269                return true;
270        }
271        // ----------------------------------------------------------------------
272        // reverse mode routine called by CppAD
273        bool mat_mul_reverse(
274                size_t                   id ,
275                size_t                    k ,
276                size_t                    n ,
277                size_t                    m ,
278                const vector<double>&    tx ,
279                const vector<double>&    ty ,
280                vector<double>&          px ,
281                const vector<double>&    py
282        )
283        {       get_info(id, k, n, m);
284
285                size_t ell = n * n_order_;
286                while(ell--)
287                        px[ell] = 0.;
288
289                size_t order = n_order_;
290                while(order--)
291                {       // reverse sum the products for specified order
292                        for(ell = 0; ell <=order; ell++)
293                                reverse_multiply(ell, order-ell, tx, ty, px, py);
294                }
295
296                // All orders are implemented and there are no possible error
297                // conditions, so always return true.
298                return true;
299        }
300
301        // ----------------------------------------------------------------------
302        // forward Jacobian sparsity routine called by CppAD
303        bool mat_mul_for_jac_sparse(
304                size_t                               id ,             
305                size_t                                n ,
306                size_t                                m ,
307                size_t                                q ,
308                const vector< std::set<size_t> >&     r ,
309                vector< std::set<size_t> >&           s )
310        {       size_t i, j, k, im_left, middle, mj_right, ij_result;
311                k = 0;
312                get_info(id, k, n, m);
313       
314                for(i = 0; i < nr_result_; i++)
315                {       for(j = 0; j < nc_result_; j++)
316                        {       ij_result = result(i, j, k);
317                                s[ij_result].clear();
318                                for(middle = 0; middle < n_middle_; middle++)
319                                {       im_left   = left(i, middle, k);
320                                        mj_right  = right(middle, j, k);
321
322                                        // s[ij_result] = union( s[ij_result], r[im_left] )
323                                        my_union(s[ij_result], s[ij_result], r[im_left]);
324
325                                        // s[ij_result] = union( s[ij_result], r[mj_right] )
326                                        my_union(s[ij_result], s[ij_result], r[mj_right]);
327                                }
328                        }
329                }
330                return true;
331        }
332        // ----------------------------------------------------------------------
333        // reverse Jacobian sparsity routine called by CppAD
334        bool mat_mul_rev_jac_sparse(
335                size_t                               id ,             
336                size_t                                n ,
337                size_t                                m ,
338                size_t                                q ,
339                vector< std::set<size_t> >&           r ,
340                const vector< std::set<size_t> >&     s )
341        {       size_t i, j, k, im_left, middle, mj_right, ij_result;
342                k = 0;
343                get_info(id, k, n, m);
344       
345                for(j = 0; j < n; j++)
346                        r[j].clear();
347
348                for(i = 0; i < nr_result_; i++)
349                {       for(j = 0; j < nc_result_; j++)
350                        {       ij_result = result(i, j, k);
351                                for(middle = 0; middle < n_middle_; middle++)
352                                {       im_left   = left(i, middle, k);
353                                        mj_right  = right(middle, j, k);
354
355                                        // r[im_left] = union( r[im_left], s[ij_result] )
356                                        my_union(r[im_left], r[im_left], s[ij_result]);
357
358                                        // r[mj_right] = union( r[mj_right], s[ij_result] )
359                                        my_union(r[mj_right], r[mj_right], s[ij_result]);
360                                }
361                        }
362                }
363                return true;
364        }
365        // ----------------------------------------------------------------------
366        // reverse Hessian sparsity routine called by CppAD
367        bool mat_mul_rev_hes_sparse(
368                size_t                               id ,             
369                size_t                                n ,
370                size_t                                m ,
371                size_t                                q ,
372                const vector< std::set<size_t> >&     r ,
373                const vector<bool>&                   s ,
374                vector<bool>&                         t ,
375                const vector< std::set<size_t> >&     u ,
376                vector< std::set<size_t> >&           v )
377        {       size_t i, j, k, im_left, middle, mj_right, ij_result;
378                k = 0;
379                get_info(id, k, n, m);
380       
381                for(j = 0; j < n; j++)
382                {       t[j] = false;   
383                        v[j].clear();
384                }
385
386                assert( vx_->size() == n );
387                for(i = 0; i < nr_result_; i++)
388                {       for(j = 0; j < nc_result_; j++)
389                        {       ij_result = result(i, j, k);
390                                for(middle = 0; middle < n_middle_; middle++)
391                                {       im_left   = left(i, middle, k);
392                                        mj_right  = right(middle, j, k);
393
394                                        // back propagate Jacobian sparsity
395                                        t[im_left]   = (t[im_left] | s[ij_result]);
396                                        t[mj_right]  = (t[mj_right] | s[ij_result]);
397                                        // Visual Studio C++ 2008 warns unsafe mix of int and
398                                        // bool if we use the following code directly above:
399                                        // t[im_left]  |= s[ij_result];
400                                        // t[mj_right] |= s[ij_result];
401
402                                        // back propagate Hessian sparsity
403                                        // v[im_left]  = union( v[im_left],  u[ij_result] )
404                                        // v[mj_right] = union( v[mj_right], u[ij_result] )
405                                        my_union(v[im_left],  v[im_left],  u[ij_result] );
406                                        my_union(v[mj_right], v[mj_right], u[ij_result] );
407
408                                        // Check for case where the (i,j) result element
409                                        // is in reverse Jacobian and both left and right
410                                        // operands in multiplication are variables
411                                        if(s[ij_result] & (*vx_)[im_left] & (*vx_)[mj_right])
412                                        {       // v[im_left] = union( v[im_left], r[mj_right] )
413                                                my_union(v[im_left], v[im_left], r[mj_right] );
414                                                // v[mj_right] = union( v[mj_right], r[im_left] )
415                                                my_union(v[mj_right], v[mj_right], r[im_left] );
416                                        }
417                                }
418                        }
419                }
420                return true;
421        }
422/* $$
423
424$head Declare mat_mul Function$$
425Declare the $code AD<double>$$ routine $codei%mat_mul(%id%, %ax%, %ay%)%$$
426and end empty namespace
427(we could use any $cref/simple vector template class/SimpleVector/$$
428instead of $code CppAD::vector$$):
429$codep */
430        CPPAD_USER_ATOMIC(
431                mat_mul                 , 
432                CppAD::vector           ,
433                double                  , 
434                mat_mul_forward         , 
435                mat_mul_reverse         ,
436                mat_mul_for_jac_sparse  ,
437                mat_mul_rev_jac_sparse  ,
438                mat_mul_rev_hes_sparse 
439        )
440} // End empty namespace
441/* $$
442$end
443*/
444
445# endif
Note: See TracBrowser for help on using the repository browser.