1 /* --------------------------------------------------------------------------
2 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-20 Bradley M. Bell
3
4 CppAD is distributed under the terms of the
5 Eclipse Public License Version 2.0.
6
7 This Source Code may also be made available under the following
8 Secondary License when the conditions for such availability set forth
9 in the Eclipse Public License, Version 2.0 are satisfied:
10 GNU General Public License, Version 2.0 or later.
11 ---------------------------------------------------------------------------- */
12 /*
13 $begin adolc_det_minor.cpp$$
14 $spell
15 onetape
16 cppad
17 zos
18 fos
19 adouble
20 CppAD
21 typedef
22 adolc
23 Lu
24 Adolc
25 det
26 hpp
27 const
28 bool
29 srand
30 $$
31
32 $section Adolc Speed: Gradient of Determinant by Minor Expansion$$
33
34
35 $head Specifications$$
36 See $cref link_det_minor$$.
37
38 $head Implementation$$
39 $srccode%cpp% */
40 // suppress conversion warnings before other includes
41 # include <cppad/wno_conversion.hpp>
42 //
43 # include <adolc/adolc.h>
44 # include <cppad/utility/vector.hpp>
45 # include <cppad/speed/det_by_minor.hpp>
46 # include <cppad/speed/uniform_01.hpp>
47
48 // list of possible options
49 # include <map>
50 extern std::map<std::string, bool> global_option;
51
52 namespace {
setup(int tag,size_t size,const CppAD::vector<double> & matrix)53 void setup(int tag, size_t size, const CppAD::vector<double>& matrix)
54 { // number of independent variables
55 int n = size * size;
56
57 // object for computing determinant
58 CppAD::det_by_minor<adouble> a_det(size);
59
60 // declare independent variables
61 int keep = 1; // keep forward mode results
62 trace_on(tag, keep);
63 CppAD::vector<adouble> a_A(n);
64 for(int j = 0; j < n; ++j)
65 a_A[j] <<= matrix[j];
66
67 // AD computation of the determinant
68 adouble a_detA = a_det(a_A);
69
70 // create function object f : A -> detA
71 double f;
72 a_detA >>= f;
73 trace_off();
74 }
75 }
76
link_det_minor(const std::string & job,size_t size,size_t repeat,CppAD::vector<double> & matrix,CppAD::vector<double> & gradient)77 bool link_det_minor(
78 const std::string& job ,
79 size_t size ,
80 size_t repeat ,
81 CppAD::vector<double> &matrix ,
82 CppAD::vector<double> &gradient )
83 {
84 // --------------------------------------------------------------------
85 // check global options
86 // Allow colpack true even though it is not used below because it is
87 // true durng the adolc correctness tests.
88 const char* valid[] = { "onetape", "optimize", "colpack"};
89 size_t n_valid = sizeof(valid) / sizeof(valid[0]);
90 typedef std::map<std::string, bool>::iterator iterator;
91 //
92 for(iterator itr=global_option.begin(); itr!=global_option.end(); ++itr)
93 { if( itr->second )
94 { bool ok = false;
95 for(size_t i = 0; i < n_valid; i++)
96 ok |= itr->first == valid[i];
97 if( ! ok )
98 return false;
99 }
100 }
101 // -----------------------------------------------------
102 // size corresponding to current tape
103 static size_t static_size = 0;
104 //
105 // number of independent variables
106 int n = size * size;
107 //
108 // tape identifier
109 int tag = 0;
110 //
111 bool onetape = global_option["onetape"];
112 // ----------------------------------------------------------------------
113 if( job == "setup" )
114 { if( onetape )
115 { // get a matrix
116 CppAD::uniform_01(size_t(n), matrix);
117 //
118 // recrod the tape
119 setup(tag, size, matrix);
120 static_size = size;
121 }
122 else
123 { static_size = 0;
124 }
125 return true;
126 }
127 if( job == "teardown" )
128 { // 2DO: How does one free an adolc tape ?
129 return true;
130 }
131 // ----------------------------------------------------------------------
132 CPPAD_ASSERT_UNKNOWN( job == "run" );
133 //
134 // number of dependent variables
135 int m = 1;
136 //
137 // vectors of reverse mode weights
138 CppAD::vector<double> u(m);
139 u[0] = 1.;
140 //
141 if( onetape ) while(repeat--)
142 { if( size != static_size )
143 { CPPAD_ASSERT_UNKNOWN( size == static_size );
144 }
145
146 // choose a matrix
147 CppAD::uniform_01(n, matrix);
148
149 // evaluate the determinant at the new matrix value
150 int keep = 1; // keep this forward mode result
151 double f; // function result
152 zos_forward(tag, m, n, keep, matrix.data(), &f);
153
154 // evaluate and return gradient using reverse mode
155 fos_reverse(tag, m, n, u.data(), gradient.data());
156 }
157 else while(repeat--)
158 {
159 // choose a matrix
160 CppAD::uniform_01(n, matrix);
161
162 // record the tape
163 setup(tag, size, matrix);
164
165 // evaluate and return gradient using reverse mode
166 fos_reverse(tag, m, n, u.data(), gradient.data());
167 }
168 // --------------------------------------------------------------------
169 return true;
170 }
171 /* %$$
172 $end
173 */
174