source: trunk/test_more/old_mat_mul.hpp @ 3757

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

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

commit aec9ede58726ce1a40825eac63d4d8f7cb40bfba
Author: Brad Bell <bradbell@…>
Date: Sun Nov 29 21:24:17 2015 -0800

Move omh/libaray to omh/utility becasue it corresponds to cppad/utility.

commit 17b6e02f29dcc2a1ab2d2e6980a979374e73e19e
Author: Brad Bell <bradbell@…>
Date: Sun Nov 29 19:03:22 2015 -0800

Batch edit the entire souce directory using the script

bin/batch_edit.sh

  1. With a few exceptions, move cppad/*.hpp to cppad/utility/*.hpp
  2. Change include gaurd from CPPAD_<NAME>_INCLUDED to CPPAD_<NAME>_HPP
  3. Make the source code control $Id$ commands more uniform.
  4. Remove invisible white space
  5. Create utility.hpp (includes entire utilty directory)

commit 4a8fbe249f97116cac86583a7d9bf596dbb06473
Author: Brad Bell <bradbell@…>
Date: Sun Nov 29 18:12:50 2015 -0800

new version that passes bin/check_all.sh

commit f041069fe3e24d68133f829353f7200cf2c299ed
Author: Brad Bell <bradbell@…>
Date: Sun Nov 29 14:27:08 2015 -0800

pases bin/check_all.sh

commit 32ac61696305b7625938d4c488203cdcc19f36dc
Author: Brad Bell <bradbell@…>
Date: Sun Nov 29 06:13:56 2015 -0800

passed bin/check_all.sh

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