00001 # ifndef CPPAD_MULTI_NEWTON_INCLUDED 00002 # define CPPAD_MULTI_NEWTON_INCLUDED 00003 00004 /* -------------------------------------------------------------------------- 00005 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 Bradley M. Bell 00006 00007 CppAD is distributed under multiple licenses. This distribution is under 00008 the terms of the 00009 Common Public License Version 1.0. 00010 00011 A copy of this license is included in the COPYING file of this distribution. 00012 Please visit http://www.coin-or.org/CppAD/ for information on other licenses. 00013 -------------------------------------------------------------------------- */ 00014 00015 /* 00016 $begin multi_newton$$ 00017 $spell 00018 xout 00019 xlow 00020 xup 00021 itr 00022 CppAD 00023 const 00024 $$ 00025 00026 $index OpenMP, Newton's method$$ 00027 $index multi-thread, Newton's method$$ 00028 $index example, OpenMP Newton's method$$ 00029 00030 $section Multi-Threaded Newton's Method Routine$$ 00031 00032 $head Syntax$$ 00033 $syntax%multi_newton(% 00034 xout%, %fun%, %n_grid%, %xlow%, %xup%, %epsilon%, %max_itr%)%$$ 00035 00036 00037 $head Purpose$$ 00038 Determine the argument values $latex x \in [a, b]$$ (where $latex a < b$$) 00039 such that $latex f(x) = 0$$. 00040 00041 $head Method$$ 00042 For $latex i = 0 , \ldots , n$$, 00043 we define the $th i$$ grid point $latex g_i$$ 00044 and the $th i$$ interval $latex I_i$$ by 00045 $latex \[ 00046 \begin{array}{rcl} 00047 g_i & = & a \frac{n - i}{n} + b \frac{i}{n} 00048 \\ 00049 I_i & = & [ g_i , g_{i+1} ] 00050 \end{array} 00051 \] $$ 00052 Newton's method is applied starting 00053 at the center of each of the intervals $latex I_i$$ for 00054 $latex i = 0 , \ldots , n-1$$ 00055 and at most one zero is found for each interval. 00056 00057 00058 $head xout$$ 00059 The argument $italic xout$$ has the prototype 00060 $syntax% 00061 CppAD::vector<double> &%xout% 00062 %$$ 00063 The input size and value of the elements of $italic xout$$ do not matter. 00064 Upon return from $code multi_newton$$, 00065 the size of $italic xout$$ is less than $latex n$$ and 00066 $latex \[ 00067 | f( xout[i] ) | \leq epsilon 00068 \] $$ 00069 for each valid index $italic i$$. 00070 Two $latex x$$ solutions are considered equal (and joined as one) if 00071 the absolute difference between the solutions is less than 00072 $latex (b - a) / n$$. 00073 00074 $head fun$$ 00075 The argument $italic fun$$ has prototype 00076 $syntax% 00077 %Fun% &%fun% 00078 %$$ 00079 This argument must evaluate the function $latex f(x)$$ 00080 using the syntax 00081 $syntax% 00082 %f% = %fun%(%x%) 00083 %$$ 00084 where the argument $italic x$$ and the result $italic f$$ 00085 have the prototypes 00086 $syntax% 00087 const AD<double> &%x% 00088 AD<double> %f% 00089 %$$. 00090 00091 00092 $head n_grid$$ 00093 The argument $italic n_grid$$ has prototype 00094 $syntax% 00095 size_t %n_grid% 00096 %$$ 00097 It specifies the number of grid points; i.e., $latex n$$ 00098 in the $cref/method/multi_newton/Method/$$ above. 00099 00100 $head xlow$$ 00101 The argument $italic xlow$$ has prototype 00102 $syntax% 00103 double %xlow% 00104 %$$ 00105 It specifies the lower limit for the entire search; i.e., $latex a$$ 00106 in the $cref/method/multi_newton/Method/$$ above. 00107 00108 $head xup$$ 00109 The argument $italic xup$$ has prototype 00110 $syntax% 00111 double %xup% 00112 %$$ 00113 It specifies the upper limit for the entire search; i.e., $latex b$$ 00114 in the $cref/method/multi_newton/Method/$$ above. 00115 00116 $head epsilon$$ 00117 The argument $italic epsilon$$ has prototype 00118 $syntax% 00119 double %epsilon% 00120 %$$ 00121 It specifies the convergence criteria for Newton's method in terms 00122 of how small the function value must be. 00123 00124 $head max_itr$$ 00125 The argument $italic max_itr$$ has prototype 00126 $syntax% 00127 size_t %max_itr% 00128 %$$ 00129 It specifies the maximum number of iterations of Newton's method to try 00130 before giving up on convergence. 00131 00132 $end 00133 --------------------------------------------------------------------------- 00134 $begin multi_newton.hpp$$ 00135 00136 $index multi_newton, source$$ 00137 $index source, multi_newton$$ 00138 $index example, OpenMP$$ 00139 $index example, multi-thread$$ 00140 $index OpenMP, example$$ 00141 $index multi-thread, example$$ 00142 00143 00144 $section OpenMP Multi-Threading Newton's Method Source Code$$ 00145 00146 $code 00147 $verbatim%openmp/multi_newton.hpp%0%// BEGIN PROGRAM%// END PROGRAM%1%$$ 00148 $$ 00149 $end 00150 --------------------------------------------------------------------------- 00151 */ 00152 // BEGIN PROGRAM 00153 00154 # include <cppad/cppad.hpp> 00155 # include <cassert> 00156 00157 # ifdef _OPENMP 00158 # include <omp.h> 00159 # endif 00160 00161 namespace { // BEGIN CppAD namespace 00162 00163 template <class Fun> 00164 void one_newton(double &fcur, double &xcur, Fun &fun, 00165 double xlow, double xin, double xup, double epsilon, size_t max_itr) 00166 { using CppAD::AD; 00167 using CppAD::vector; 00168 using CppAD::abs; 00169 00170 // domain space vector 00171 size_t n = 1; 00172 vector< AD<double> > X(n); 00173 // range space vector 00174 size_t m = 1; 00175 vector< AD<double> > Y(m); 00176 // domain and range differentials 00177 vector<double> dx(n), dy(m); 00178 00179 size_t itr; 00180 xcur = xin; 00181 for(itr = 0; itr < max_itr; itr++) 00182 { // domain space vector 00183 X[0] = xcur; 00184 CppAD::Independent(X); 00185 // range space vector 00186 Y[0] = fun(X[0]); 00187 // F : X -> Y 00188 CppAD::ADFun<double> F(X, Y); 00189 // fcur = F(xcur) 00190 fcur = Value(Y[0]); 00191 // evaluate dfcur = F'(xcur) 00192 dx[0] = 1; 00193 dy = F.Forward(1, dx); 00194 double dfcur = dy[0]; 00195 // check end of iterations 00196 if( abs(fcur) <= epsilon ) 00197 return; 00198 if( (xcur == xlow) & (fcur * dfcur > 0.) ) 00199 return; 00200 if( (xcur == xup) & (fcur * dfcur < 0.) ) 00201 return; 00202 if( dfcur == 0. ) 00203 return; 00204 // next Newton iterate 00205 double delta_x = - fcur / dfcur; 00206 if( xlow - xcur >= delta_x ) 00207 xcur = xlow; 00208 else if( xup - xcur <= delta_x ) 00209 xcur = xup; 00210 else xcur = xcur + delta_x; 00211 } 00212 return; 00213 } 00214 00215 template <class Fun> 00216 void multi_newton( 00217 CppAD::vector<double> &xout , 00218 Fun &fun , 00219 size_t n_grid , 00220 double xlow , 00221 double xup , 00222 double epsilon , 00223 size_t max_itr ) 00224 { using CppAD::AD; 00225 using CppAD::vector; 00226 using CppAD::abs; 00227 00228 // check argument values 00229 assert( xlow < xup ); 00230 assert( n_grid > 0 ); 00231 00232 // OpenMP uses integers in place of size_t 00233 int i, n = int(n_grid); 00234 00235 // set up grid 00236 vector<double> grid(n_grid + 1); 00237 vector<double> fcur(n_grid), xcur(n_grid), xmid(n_grid); 00238 double dx = (xup - xlow) / double(n_grid); 00239 for(i = 0; size_t(i) < n_grid; i++) 00240 { grid[i] = xlow + i * dx; 00241 xmid[i] = xlow + (i + .5) * dx; 00242 } 00243 grid[n_grid] = xup; 00244 00245 # ifdef _OPENMP 00246 # pragma omp parallel for 00247 # endif 00248 for(i = 0; i < n; i++) 00249 { one_newton( 00250 fcur[i] , 00251 xcur[i] , 00252 fun , 00253 grid[i] , 00254 xmid[i] , 00255 grid[i+1] , 00256 epsilon , 00257 max_itr 00258 ); 00259 } 00260 // end omp parallel for 00261 00262 // remove duplicates and points that are not solutions 00263 double xlast = xlow; 00264 size_t ilast = 0; 00265 size_t n_zero = 0; 00266 for(i = 0; size_t(i) < n_grid; i++) 00267 { if( abs( fcur[i] ) <= epsilon ) 00268 { if( n_zero == 0 ) 00269 { xcur[n_zero++] = xlast = xcur[i]; 00270 ilast = i; 00271 } 00272 else if( fabs( xcur[i] - xlast ) > dx ) 00273 { xcur[n_zero++] = xlast = xcur[i]; 00274 ilast = i; 00275 } 00276 else if( fabs( fcur[i] ) < fabs( fcur[ilast] ) ) 00277 { xcur[n_zero - 1] = xlast = xcur[i]; 00278 ilast = i; 00279 } 00280 } 00281 } 00282 00283 // resize output vector and set its values 00284 xout.resize(n_zero); 00285 for(i = 0; size_t(i) < n_zero; i++) 00286 xout[i] = xcur[i]; 00287 00288 return; 00289 } 00290 00291 } // END CppAD namespace 00292 00293 // END PROGRAM 00294 00295 # endif