00001 /*************************************************************************** 00002 newtonroot.h - description 00003 ------------------- 00004 begin : Mon May 8 2000 00005 copyright : (C) 2000 by Michael Peeters 00006 email : Michael.Peeters@vub.ac.be 00007 ***************************************************************************/ 00008 00009 /*************************************************************************** 00010 * * 00011 * This program is free software; you can redistribute it and/or modify * 00012 * it under the terms of the GNU General Public License as published by * 00013 * the Free Software Foundation; either version 2 of the License, or * 00014 * (at your option) any later version. * 00015 * * 00016 ***************************************************************************/ 00017 00018 #ifndef NEWTONROOT_H 00019 #define NEWTONROOT_H 00020 00021 #include "numerictypes.h" 00022 #include "numerictraits.h" 00023 #include <stdexcept> 00024 #include "jacobian.h" 00025 #include "linesearch.h" 00026 #include "lusolve.h" 00027 #include "utility.h" 00028 00029 // #include "debugmacro.h" 00030 00031 namespace MODEL { 00032 00036 // Previously: template <integer dims, class NT = NumericTraits<number,dims> > 00037 template<integer dims, typename nelem=number, class NT = NumericTraits<nelem,dims> > 00038 class NewtonRoot : public VectorFunction<dims,number,NT> 00039 { 00040 public: 00041 typedef typename NT::number numT; 00042 typedef typename NT::vect vect; 00043 typedef typename NT::matrix matrix; 00044 typedef typename NT::vf vf; 00045 typedef VectorFunction<dims,number,NT> base; 00046 00047 NewtonRoot(vf& f, bool own=false) 00048 : func( (own)? (f.clone()) : (&f) ), owned(own), 00049 wrongmin(false), ls(*func), fdjac(*func) {} 00050 00051 ~NewtonRoot(){ if (owned) delete func;} 00052 00054 NewtonRoot(const NewtonRoot& nr) 00055 : func(nr.func), owned(false), 00056 wrongmin(nr.wrongmin), ls(*func), fdjac(*func) {} 00057 00058 virtual NewtonRoot* clone () const {return new NewtonRoot(*this);} 00059 00064 virtual const vect& function(vect& fu,const vect& startu); 00065 00067 bool wrong_min(void){return wrongmin;} 00068 00070 bool no_root(void) 00071 { 00072 return noroot; 00073 } 00074 // no assignment allowed 00075 PRIVATE_ASSIGN(NewtonRoot); 00076 00077 00078 private: 00079 vf* func; 00080 bool owned; 00081 bool wrongmin; 00082 bool noroot; 00083 00084 LineSearch<dims,NT> ls; 00085 Jacobian<dims,NT> fdjac; 00086 00087 public: 00088 static const integer maxiterations=10000; // NRC: MAXITS 00089 static const numT tolerancef=1.E-8; // NRC: TOLF 00090 static const numT tolerancemin=1.E-6; // NRC: TOLMIN 00091 static const numT maxstep=100.; // NRC: STPMX 00092 static const numT tolerancex=1.E-8; // NRC: TOLX 00093 }; 00094 00095 // --- NON INLINED --- 00096 // template <integer dims, typename nelem, class NT > 00097 // NewtonRoot<dims,nelem,NT>::NewtonRoot(const NewtonRoot& nr) : ls(nr.ls),fdjac(nr.fdjac) {cout 00098 // << "NewtonRoot : BOOM" << endl; throw std::logic_error("Copying NewtonRoot is evil");} 00099 00100 // Prev: template <integer dims, class NT > 00101 template <integer dims, typename nelem, class NT > 00102 const NewtonRoot<dims,nelem,NT>::vect& 00103 NewtonRoot<dims,nelem,NT>::function(vect& fu,const vect& startu) 00104 { 00105 // We are looking 00106 noroot=false; 00107 00108 // Calculate the initial function value and norm 00109 numT f=ls.norm()(startu); 00110 vect fvec=ls.norm().function_value(); 00111 00112 // alias for clarity 00113 vect& u(fu); // The return vector IS the final point where we end 00114 u=startu; // And we start fro, the initial point 00115 00116 00117 // Test if we are too close to a zero 00118 numT testtol=norm(fvec); // More efficient with square 00119 // for (integer i=0;i<dims;i++) 00120 // { numT tm(0.); if ( (tm=abs(fvec[i]))> testtol) testtol=tm; } 00121 00122 if (testtol<0.01*tolerancef*tolerancef) return fu; // close enough to zero already 00123 00124 // Calculate maximal step 00125 numT sum=norm(u); 00126 numT stpmax=maxstep*max(sqrt(sum),numT(dims)); 00127 00128 for (integer its=0;its<maxiterations;its++) 00129 { 00130 // DEBUG_Cellular << "newtonroot:" << u << endl; 00131 00132 // Calculate Jacobian 00133 matrix j=fdjac.calculate(u,fvec); 00134 00135 vect gradient(0.); 00136 00137 // Calculate Gradient 00138 for (integer i=0;i<dims;i++) 00139 { 00140 numT sum(0.); 00141 for (integer k=0;k<dims;k++) sum += j[k][i]*fvec[k]; 00142 gradient[i]=sum; 00143 } 00144 00145 // Keep previous point 00146 vect uold(u); 00147 numT fold(f); 00148 00149 // Descent direction 00150 vect p(-fvec); 00151 00152 // Solve system using LU decomposition 00153 // ludcmp(fjac,n,indx,&d); 00154 LUSolve<dims,NT> lus(j); // init 00155 00156 lus.solve(p); // solve by backsubstitution 00157 00158 ls(uold,fold,gradient,p,u,f,stpmax); // Linesearch with this value 00159 00160 fvec=ls.norm().function_value(); // Get value (without recalculating) 00161 00162 wrongmin=ls.converged(); // Error coming from ls 00163 00164 // Test if we are too close to a zero 00165 //testtol=0.; 00166 //for (integer i=0;i<dims;i++) 00167 // { numT tm(0.); if ( (tm=abs(fvec[i]) ) > testtol) testtol=tm; } 00168 00169 testtol=norm(fvec); 00170 if (testtol<tolerancef*tolerancef) 00171 { 00172 // DEBUG_Grainy << "newtonroot: exited normally" << endl; 00173 wrongmin=false; 00174 return fu; 00175 } // close enough to zero already 00176 00177 // Make sure we are not at a spurious convergence (grad f=0) 00178 if (wrongmin) 00179 { 00180 // Besides, whenever we arrive here, something boogery has happened 00181 testtol=0.0; 00182 numT den=max(f,0.5*number(dims)); 00183 numT tm(0.); 00184 for (integer i=0;i<dims;i++) 00185 { 00186 tm=abs(gradient[i])*max(abs(u[i]),number(1.0))/den; 00187 if (tm > testtol) testtol=tm; 00188 } 00189 wrongmin = (testtol < tolerancemin)||(!(norm(fvec)<tolerancef*tolerancef)); 00190 // DEBUG_Grainy << "newtonroot: gradient?" << gradient << 00191 // "zero:" << fvec <<endl; 00192 // DEBUG_Grainy<< "newtonroot: exited thru spurious 00193 // convergence"<<endl; 00194 return fu; 00195 } 00196 00197 testtol =0.0; // check for converence on u 00198 numT tm(0.); 00199 for (integer i=0;i<dims;i++) { 00200 tm=(abs(u[i]-uold[i]))/max(abs(u[i]),number(1.0)); 00201 if (tm > testtol) testtol=tm; 00202 } 00203 if (testtol < tolerancex) 00204 { 00205 // DEBUG_Grainy << "newtonroot: converged on u"<<endl; 00206 return fu; 00207 } 00208 } 00209 noroot=true; 00210 throw std::logic_error("NewtonRoot: maxiteration exceeded!"); 00211 return fu=vect(0.); // BAD BAD BAD we should never come here 00212 } 00213 00214 00215 00216 } // end MODEL 00217 #endif 00218 00219 /********************************************************************* 00220 $Id: newtonroot.h,v 1.1.2.1 2002/02/06 12:48:04 mpeeters Exp $ 00221 ********************************************************************** 00222 00223 $Log: newtonroot.h,v $ 00224 Revision 1.1.2.1 2002/02/06 12:48:04 mpeeters 00225 Changes minor type error. 00226 00227 Revision 1.1 2001/05/22 10:54:55 mpeeters 00228 Moved sources and headers for libModel to model/ subdirectory, in an attempt to rationalize the source tree. This should make things "netter". 00229 00230 Revision 1.4 2001/05/21 11:53:16 mpeeters 00231 Removed Makefile.in, which is automatically generated anyway. 00232 00233 Revision 1.3 2000/09/22 08:54:45 mpeeters 00234 Change of default tolerances when looking for roots. 00235 00236 Revision 1.2 2000/09/15 10:26:31 mpeeters 00237 Cleaned out header and added CVS tails to files, removed superfuous 00238 @author comments, inserted dates 00239 00240 *********************************************************************/
More Info? Michael Peeters. Also, check our research website: www.alna.vub.ac.be
Last update: June 2002.
Looking for Open Source? Check out
!