source: trunk/test_more/old_mat_mul.hpp @ 3788

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

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

commit 431e0a227dbfe6172d265d9d79a2b5b258c5bc18
Author: Brad Bell <bradbell@…>
Date: Tue Feb 9 07:41:29 2016 -0700

  1. Change package.sh to automatically update version on (master branch only).
  2. Change version.sh copy to check and change.

commit a14455414810cfe3c3e4bca90090defc2528a353
Author: Brad Bell <bradbell@…>
Date: Tue Feb 9 06:19:54 2016 -0700

Change check_verbatim to check_srcfile
because all but one of the verbatim commands were changed to srcfile commands.


check_include_omh.sh: old check that file names did not change case (for cygwin development).

commit 4ce45b796b57629332ab46d8ae6df94e0a1ed998
Author: Brad Bell <bradbell@…>
Date: Tue Feb 9 06:04:57 2016 -0700

batch_edit.sh to change $verbatim and $codep to $srcfile and $srccode.


det_by_minor.c: remove some invisible white space.

commit 56553b88c9623c30d2222425a9640b95ce4c8281
Author: Brad Bell <bradbell@…>
Date: Mon Feb 8 18:01:49 2016 -0700

check_jenkins.sh: jenkins.sh no longer takes an argument.
jenkins.sh: fix name of script in error message.

commit 3b8a208cfc7e8ef3c928c17eb291aa3b90ff0050
Author: Brad Bell <bradbell@…>
Date: Mon Feb 8 07:57:02 2016 -0700

new_release.sh: track branches in comments, back to master at OK end, first check of response.

commit 442b7cbc45c022776e8257d3c3404dccdd06c420
Author: Brad Bell <bradbell@…>
Date: Mon Feb 8 06:01:11 2016 -0700

  1. Advance to release 20160000.1.
  2. Check using master version of new_release.sh.
  3. Make sure auto-tools version up to date.
  4. Ask user if doing further testing before commiting new release.

commit f7bdd1f48e72feb05d604da63914022809f45c28
Author: Brad Bell <bradbell@…>
Date: Sun Feb 7 07:59:41 2016 -0700

Add shared library version number to cppad_lib; i.e., cppad_lib.yyyy.mmdd.rel

commit a4c716552e3ad05b337aea58b643c9ad1cbcd4ac
Author: Brad Bell <bradbell@…>
Date: Sun Feb 7 05:25:39 2016 -0700

Make cppad_lib libarary conditional on colpack_prefix being specified.

commit 5e8890eb8de8b0cde146a6ed59c391d7c355ff24
Author: Brad Bell <bradbell@…>
Date: Tue Jan 26 10:49:37 2016 -0700

vector.hpp: fix The -> This.

commit e4e5442b069d7b00e197c31616da32eee20460b3
Merge: c4559d5 ed28b89
Author: Brad Bell <bradbell@…>
Date: Tue Jan 26 09:47:58 2016 -0700

Merge pull request #14 from barak/master


fix spelling in description of cppad_profile_flag on cmake command line

commit ed28b899c9fedab52a578aa7dd73818638081fe6
Author: Barak A. Pearlmutter <barak+git@…>
Date: Tue Jan 26 16:24:32 2016 +0000

typo

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-16 Bradley M. Bell
7
8CppAD is distributed under multiple licenses. This distribution is under
9the terms of the
10                    Eclipse Public License Version 1.0.
11
12A copy of this license is included in the COPYING file of this distribution.
13Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
14-------------------------------------------------------------------------- */
15
16/*
17$begin 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$srccode%cpp% */
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$srccode%cpp% */
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$srccode%cpp% */
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$srccode%cpp% */
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$srccode%cpp% */
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$srccode%cpp% */
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$srccode%cpp% */
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$srccode%cpp% */
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.