source: trunk/test_more/old_mat_mul.cpp @ 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: 8.1 KB
Line 
1// $Id$
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-16 Bradley M. Bell
4
5CppAD is distributed under multiple licenses. This distribution is under
6the terms of the
7                    Eclipse Public License Version 1.0.
8
9A copy of this license is included in the COPYING file of this distribution.
10Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
11-------------------------------------------------------------------------- */
12
13/*
14$begin old_mat_mul.cpp$$
15$spell
16        mul
17$$
18
19$section Old Matrix Multiply as a User Atomic Operation: Example and Test$$
20
21$head Deprecated 2013-05-27$$
22This example has been deprecated;
23see $cref atomic_mat_mul.cpp$$.
24
25$children%
26        test_more/old_mat_mul.hpp
27%$$
28$head Include File$$
29This routine uses the include file $cref old_mat_mul.hpp$$.
30
31$code
32$srcfile%test_more/old_mat_mul.cpp%0%// BEGIN C++%// END C++%1%$$
33$$
34
35$end
36*/
37// BEGIN C++
38# include <cppad/cppad.hpp>
39# include "old_mat_mul.hpp"
40
41bool old_mat_mul(void)
42{       bool ok = true;
43        using CppAD::AD;
44
45        // matrix sizes for this test
46        size_t nr_result = 2;
47        size_t n_middle  = 2;
48        size_t nc_result = 2;
49
50        // declare the AD<double> vectors ax and ay and X
51        size_t n = nr_result * n_middle + n_middle * nc_result;
52        size_t m = nr_result * nc_result;
53        CppAD::vector< AD<double> > X(4), ax(n), ay(m);
54        size_t i, j;
55        for(j = 0; j < X.size(); j++)
56                X[j] = (j + 1);
57
58        // X is the vector of independent variables
59        CppAD::Independent(X);
60        // left matrix
61        ax[0]  = X[0];  // left[0,0]   = x[0] = 1
62        ax[1]  = X[1];  // left[0,1]   = x[1] = 2
63        ax[2]  = 5.;    // left[1,0]   = 5
64        ax[3]  = 6.;    // left[1,1]   = 6
65        // right matrix
66        ax[4]  = X[2];  // right[0,0]  = x[2] = 3
67        ax[5]  = 7.;    // right[0,1]  = 7
68        ax[6]  = X[3];  // right[1,0]  = x[3] = 4
69        ax[7]  = 8.;    // right[1,1]  = 8
70        /*
71        [ x0 , x1 ] * [ x2 , 7 ] = [ x0*x2 + x1*x3 , x0*7 + x1*8 ]
72        [ 5  , 6 ]    [ x3 , 8 ]   [ 5*x2  + 6*x3  , 5*7 + 6*8 ]
73        */
74
75        // The call back routines need to know the dimensions of the matrices.
76        // Store information about the matrix multiply for this call to mat_mul.
77        call_info info;
78        info.nr_result = nr_result;
79        info.n_middle  = n_middle;
80        info.nc_result = nc_result;
81        // info.vx gets set by forward during call to mat_mul below
82        assert( info.vx.size() == 0 );
83        size_t id      = info_.size();
84        info_.push_back(info);
85
86        // user defined AD<double> version of matrix multiply
87        mat_mul(id, ax, ay);
88        //----------------------------------------------------------------------
89        // check AD<double>  results
90        ok &= ay[0] == (1*3 + 2*4); ok &= Variable( ay[0] );
91        ok &= ay[1] == (1*7 + 2*8); ok &= Variable( ay[1] );
92        ok &= ay[2] == (5*3 + 6*4); ok &= Variable( ay[2] );
93        ok &= ay[3] == (5*7 + 6*8); ok &= Parameter( ay[3] );
94        //----------------------------------------------------------------------
95        // use mat_mul to define a function g : X -> ay
96        CppAD::ADFun<double> G;
97        G.Dependent(X, ay);
98        // g(x) = [ x0*x2 + x1*x3 , x0*7 + x1*8 , 5*x2  + 6*x3  , 5*7 + 6*8 ]^T
99        //----------------------------------------------------------------------
100        // Test zero order forward mode evaluation of g(x)
101        CppAD::vector<double> x( X.size() ), y(m);
102        for(j = 0; j <  X.size() ; j++)
103                x[j] = j + 2;
104        y = G.Forward(0, x);
105        ok &= y[0] == x[0] * x[2] + x[1] * x[3];
106        ok &= y[1] == x[0] * 7.   + x[1] * 8.;
107        ok &= y[2] == 5. * x[2]   + 6. * x[3];
108        ok &= y[3] == 5. * 7.     + 6. * 8.;
109
110        //----------------------------------------------------------------------
111        // Test first order forward mode evaluation of g'(x) * [1, 2, 3, 4]^T
112        // g'(x) = [ x2, x3, x0, x1 ]
113        //         [ 7 ,  8,  0, 0  ]
114        //         [ 0 ,  0,  5, 6  ]
115        //         [ 0 ,  0,  0, 0  ]
116        CppAD::vector<double> dx( X.size() ), dy(m);
117        for(j = 0; j <  X.size() ; j++)
118                dx[j] = j + 1;
119        dy = G.Forward(1, dx);
120        ok &= dy[0] == 1. * x[2] + 2. * x[3] + 3. * x[0] + 4. * x[1];
121        ok &= dy[1] == 1. * 7.   + 2. * 8.   + 3. * 0.   + 4. * 0.;
122        ok &= dy[2] == 1. * 0.   + 2. * 0.   + 3. * 5.   + 4. * 6.;
123        ok &= dy[3] == 1. * 0.   + 2. * 0.   + 3. * 0.   + 4. * 0.;
124
125        //----------------------------------------------------------------------
126        // Test second order forward mode
127        // g_0^2 (x) = [ 0, 0, 1, 0 ], g_0^2 (x) * [1] = [3]
128        //             [ 0, 0, 0, 1 ]              [2]   [4]
129        //             [ 1, 0, 0, 0 ]              [3]   [1]
130        //             [ 0, 1, 0, 0 ]              [4]   [2]
131        CppAD::vector<double> ddx( X.size() ), ddy(m);
132        for(j = 0; j <  X.size() ; j++)
133                ddx[j] = 0.;
134        ddy = G.Forward(2, ddx);
135        // [1, 2, 3, 4] * g_0^2 (x) * [1, 2, 3, 4]^T = 1*3 + 2*4 + 3*1 + 4*2
136        ok &= 2. * ddy[0] == 1. * 3. + 2. * 4. + 3. * 1. + 4. * 2.;
137        // for i > 0, [1, 2, 3, 4] * g_i^2 (x) * [1, 2, 3, 4]^T = 0
138        ok &= ddy[1] == 0.;
139        ok &= ddy[2] == 0.;
140        ok &= ddy[3] == 0.;
141
142        //----------------------------------------------------------------------
143        // Test second order reverse mode
144        CppAD::vector<double> w(m), dw(2 *  X.size() );
145        for(i = 0; i < m; i++)
146                w[i] = 0.;
147        w[0] = 1.;
148        dw = G.Reverse(2, w);
149        // g_0'(x) = [ x2, x3, x0, x1 ]
150        ok &= dw[0*2 + 0] == x[2];
151        ok &= dw[1*2 + 0] == x[3];
152        ok &= dw[2*2 + 0] == x[0];
153        ok &= dw[3*2 + 0] == x[1];
154        // g_0'(x)   * [1, 2, 3, 4]  = 1 * x2 + 2 * x3 + 3 * x0 + 4 * x1
155        // g_0^2 (x) * [1, 2, 3, 4]  = [3, 4, 1, 2]
156        ok &= dw[0*2 + 1] == 3.;
157        ok &= dw[1*2 + 1] == 4.;
158        ok &= dw[2*2 + 1] == 1.;
159        ok &= dw[3*2 + 1] == 2.;
160
161        //----------------------------------------------------------------------
162        // Test forward and reverse Jacobian sparsity pattern
163        /*
164        [ x0 , x1 ] * [ x2 , 7 ] = [ x0*x2 + x1*x3 , x0*7 + x1*8 ]
165        [ 5  , 6 ]    [ x3 , 8 ]   [ 5*x2  + 6*x3  , 5*7 + 6*8 ]
166        so the sparsity pattern should be
167        s[0] = {0, 1, 2, 3}
168        s[1] = {0, 1}
169        s[2] = {2, 3}
170        s[3] = {}
171        */
172        CppAD::vector< std::set<size_t> > r( X.size() ), s(m);
173        for(j = 0; j <  X.size() ; j++)
174        {       assert( r[j].empty() );
175                r[j].insert(j);
176        }
177        s = G.ForSparseJac( X.size() , r);
178        for(j = 0; j <  X.size() ; j++)
179        {       // s[0] = {0, 1, 2, 3}
180                ok &= s[0].find(j) != s[0].end();
181                // s[1] = {0, 1}
182                if( j == 0 || j == 1 )
183                        ok &= s[1].find(j) != s[1].end();
184                else    ok &= s[1].find(j) == s[1].end();
185                // s[2] = {2, 3}
186                if( j == 2 || j == 3 )
187                        ok &= s[2].find(j) != s[2].end();
188                else    ok &= s[2].find(j) == s[2].end();
189        }
190        // s[3] == {}
191        ok &= s[3].empty();
192
193        //----------------------------------------------------------------------
194        // Test reverse Jacobian sparsity pattern
195        /*
196        [ x0 , x1 ] * [ x2 , 7 ] = [ x0*x2 + x1*x3 , x0*7 + x1*8 ]
197        [ 5  , 6 ]    [ x3 , 8 ]   [ 5*x2  + 6*x3  , 5*7 + 6*8 ]
198        so the sparsity pattern should be
199        r[0] = {0, 1, 2, 3}
200        r[1] = {0, 1}
201        r[2] = {2, 3}
202        r[3] = {}
203        */
204        for(i = 0; i <  m; i++)
205        {       s[i].clear();
206                s[i].insert(i);
207        }
208        r = G.RevSparseJac(m, s);
209        for(j = 0; j <  X.size() ; j++)
210        {       // r[0] = {0, 1, 2, 3}
211                ok &= r[0].find(j) != r[0].end();
212                // r[1] = {0, 1}
213                if( j == 0 || j == 1 )
214                        ok &= r[1].find(j) != r[1].end();
215                else    ok &= r[1].find(j) == r[1].end();
216                // r[2] = {2, 3}
217                if( j == 2 || j == 3 )
218                        ok &= r[2].find(j) != r[2].end();
219                else    ok &= r[2].find(j) == r[2].end();
220        }
221        // r[3] == {}
222        ok &= r[3].empty();
223
224        //----------------------------------------------------------------------
225        /* Test reverse Hessian sparsity pattern
226        g_0^2 (x) = [ 0, 0, 1, 0 ] and for i > 0, g_i^2 = 0
227                    [ 0, 0, 0, 1 ]
228                    [ 1, 0, 0, 0 ]
229                    [ 0, 1, 0, 0 ]
230        so for the sparsity pattern for the first component of g is
231        h[0] = {2}
232        h[1] = {3}
233        h[2] = {0}
234        h[3] = {1}
235        */
236        CppAD::vector< std::set<size_t> > h( X.size() ), t(1);
237        t[0].clear();
238        t[0].insert(0);
239        h = G.RevSparseHes(X.size() , t);
240        size_t check[] = {2, 3, 0, 1};
241        for(j = 0; j <  X.size() ; j++)
242        {       // h[j] = { check[j] }
243                for(i = 0; i < n; i++)
244                {       if( i == check[j] )
245                                ok &= h[j].find(i) != h[j].end();
246                        else    ok &= h[j].find(i) == h[j].end();
247                }
248        }
249        t[0].clear();
250        for( j = 1; j < X.size(); j++)
251                        t[0].insert(j);
252        h = G.RevSparseHes(X.size() , t);
253        for(j = 0; j <  X.size() ; j++)
254        {       // h[j] = { }
255                for(i = 0; i < X.size(); i++)
256                        ok &= h[j].find(i) == h[j].end();
257        }
258
259        // --------------------------------------------------------------------
260        // Free temporary work space. (If there are future calls to
261        // old_mat_mul they would create new temporary work space.)
262        CppAD::user_atomic<double>::clear();
263        info_.clear();
264
265        return ok;
266}
267// END C++
Note: See TracBrowser for help on using the repository browser.