CppAD: A C++ Algorithmic Differentiation Package
20130918
|
00001 /* $Id$ */ 00002 /* -------------------------------------------------------------------------- 00003 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell 00004 00005 CppAD is distributed under multiple licenses. This distribution is under 00006 the terms of the 00007 Eclipse Public License Version 1.0. 00008 00009 A copy of this license is included in the COPYING file of this distribution. 00010 Please visit http://www.coin-or.org/CppAD/ for information on other licenses. 00011 -------------------------------------------------------------------------- */ 00012 00013 /* 00014 $begin link_ode$$ 00015 $spell 00016 Jacobian 00017 fp 00018 bool 00019 CppAD 00020 $$ 00021 00022 $index link_ode$$ 00023 $index ode, speed test$$ 00024 $index speed, test ode$$ 00025 $index test, ode speed$$ 00026 00027 $section Speed Testing the Jacobian of Ode Solution$$ 00028 00029 $head Prototype$$ 00030 $codei%extern bool link_ode( 00031 size_t %size% , 00032 size_t %repeat% , 00033 CppAD::vector<double> &%x% , 00034 CppAD::vector<double> &%jacobian% 00035 ); 00036 %$$ 00037 00038 $head Purpose$$ 00039 Each $cref/package/speed_main/package/$$ 00040 must define a version of this routine as specified below. 00041 This is used by the $cref speed_main$$ program 00042 to run the corresponding speed and correctness tests. 00043 00044 $head Method$$ 00045 The same template routine $cref ode_evaluate$$ is used 00046 by th different AD packages. 00047 00048 $head f$$ 00049 The function 00050 $latex f : \B{R}^n \rightarrow \B{R}^n$$ that is defined and computed by 00051 evaluating $cref ode_evaluate$$ with a call of the form 00052 $codei% 00053 ode_evaluate(%x%, %p%, %fp%) 00054 %$$ 00055 with $icode p$$ equal to zero. 00056 Calls with the value $icode p$$ equal to one are used to check 00057 the derivative values. 00058 00059 $head Return Value$$ 00060 If this speed test is not yet 00061 supported by a particular $icode package$$, 00062 the corresponding return value for $code link_ode$$ 00063 should be $code false$$. 00064 00065 $head size$$ 00066 The argument $icode size$$ 00067 is the number of variables in the ordinary differential equations 00068 which is also equal to $latex n$$. 00069 00070 $head repeat$$ 00071 The argument $icode repeat$$ is the number of times the 00072 Jacobian is computed. 00073 00074 $head x$$ 00075 The argument $icode x$$ is a vector with $latex n$$ elements. 00076 The input value of the elements of $icode x$$ does not matter. 00077 On output, it has been set to the 00078 argument value for which the function, 00079 or its derivative, is being evaluated. 00080 The value of this vector must change with each repetition. 00081 00082 $head jacobian$$ 00083 The argument $icode jacobian$$ is a vector with $latex n^2$$ elements. 00084 The input value of its elements does not matter. 00085 The output value of its elements is the Jacobian of the function $latex f(x)$$ 00086 that corresponds to output value of $icode x$$. 00087 To be more specific, for 00088 $latex i = 0 , \ldots , n-1$$ and $latex j = 0 , \ldots , n-1$$, 00089 $latex \[ 00090 \D{f[i]}{x[j]} (x) = jacobian [ i \cdot n + j ] 00091 \] $$ 00092 00093 $subhead double$$ 00094 In the case where $icode package$$ is $code double$$, 00095 only the first $latex n$$ element of $icode jacobian$$ 00096 are modified and they are to the function value 00097 $latex f(x)$$ corresponding to the output value of $icode x$$. 00098 00099 $end 00100 ----------------------------------------------------------------------------- 00101 */ 00102 # include <cppad/vector.hpp> 00103 # include <cppad/speed/ode_evaluate.hpp> 00104 # include <cppad/near_equal.hpp> 00105 00106 extern bool link_ode( 00107 size_t size , 00108 size_t repeat , 00109 CppAD::vector<double> &x , 00110 CppAD::vector<double> &jacobian 00111 ); 00112 bool available_ode(void) 00113 { size_t n = 1; 00114 size_t repeat = 1; 00115 CppAD::vector<double> x(n); 00116 CppAD::vector<double> jacobian(n * n); 00117 00118 return link_ode(n, repeat, x, jacobian); 00119 } 00120 bool correct_ode(bool is_package_double) 00121 { bool ok = true; 00122 00123 size_t n = 5; 00124 size_t repeat = 1; 00125 CppAD::vector<double> x(n); 00126 CppAD::vector<double> jacobian(n * n); 00127 00128 link_ode(n, repeat, x, jacobian); 00129 00130 size_t size = n * n; 00131 size_t p = 1; 00132 if( is_package_double ) 00133 { p = 0; // check function value 00134 size = n; 00135 } 00136 CppAD::vector<double> check(size); 00137 CppAD::ode_evaluate(x, p, check); 00138 size_t k; 00139 for(k = 0; k < size; k++) 00140 ok &= CppAD::NearEqual(check[k], jacobian[k], 1e-6, 1e-6); 00141 00142 return ok; 00143 } 00144 void speed_ode(size_t n, size_t repeat) 00145 { 00146 CppAD::vector<double> x(n); 00147 CppAD::vector<double> jacobian(n * n); 00148 link_ode(n, repeat, x, jacobian); 00149 return; 00150 }