Sacado_trad2.hpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                           Sacado Package
00005 //                 Copyright (2007) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
00025 // (etphipp@sandia.gov).
00026 //
00027 // ***********************************************************************
00028 // @HEADER
00029 
00030 // Extension of the RAD package (Reverse Automatic Differentiation) --
00031 // a package specialized for function and gradient evaluations -- to
00032 // Hessian-vector products.
00033 // Written in 2007 by David M. Gay at Sandia National Labs, Albuquerque, NM.
00034 
00035 #ifndef SACADO_TRAD2_H
00036 #define SACADO_TRAD2_H
00037 
00038 #include "Sacado_ConfigDefs.h"
00039 #include "Sacado_trad2_Traits.hpp"
00040 
00041 #include <stddef.h>
00042 #include <cmath>
00043 #include <math.h>
00044 
00045 #ifdef RAD_Const_WARN // ==> RAD_AUTO_AD_Const and RAD_DEBUG
00046 #ifndef RAD_AUTO_AD_Const
00047 #define RAD_AUTO_AD_Const
00048 #endif
00049 #ifndef RAD_DEBUG
00050 #define RAD_DEBUG
00051 #endif
00052 extern "C" int RAD_Const_Warn(const void*);// outside any namespace for
00053           // ease in setting breakpoints
00054 #endif // RAD_Const_WARN
00055 
00056 #ifdef RAD_DEBUG
00057 #include <cstdio>
00058 #include <stdlib.h>
00059 #endif
00060 
00061 #ifndef RAD_AUTO_AD_Const
00062 #ifdef RAD_DEBUG_BLOCKKEEP
00063 #include <complex>  // must be here when SACADO_NAMESPACE is #defined
00064 #endif
00065 #endif
00066 
00067 namespace Sacado {
00068 namespace Rad2 {
00069 
00070 // -DRAD_NO_USING_STDCC is needed, e.g., with Sun CC 5.7
00071 #ifndef RAD_NO_USING_STDCC
00072   // Bring math functions into scope
00073   using std::exp;
00074   using std::log;
00075   using std::log10;
00076   using std::sqrt;
00077   using std::cos;
00078   using std::sin;
00079   using std::tan;
00080   using std::acos;
00081   using std::asin;
00082   using std::atan;
00083   using std::cosh;
00084   using std::sinh;
00085   using std::tanh;
00086   using std::abs;
00087   using std::fabs;
00088   using std::atan2;
00089   using std::pow;
00090 #endif 
00091 
00092 #ifdef RAD_AUTO_AD_Const
00093 #undef RAD_DEBUG_BLOCKKEEP
00094 #else 
00095 #ifdef RAD_DEBUG_BLOCKKEEP
00096 #if !(RAD_DEBUG_BLOCKKEEP > 0)
00097 #undef RAD_DEBUG_BLOCKKEEP
00098 #else
00099 extern "C" void _uninit_f2c(void *x, int type, long len);
00100 
00101 template <typename T>
00102 struct UninitType {};
00103 
00104 template <>
00105 struct UninitType<float> {
00106   static const int utype = 4;
00107 };
00108 
00109 template <>
00110 struct UninitType<double> {
00111   static const int utype = 5;
00112 };
00113 
00114 template <typename T>
00115 struct UninitType< std::complex<T> > {
00116   static const int utype = UninitType<T>::utype + 2;
00117 };
00118 
00119 #endif /*RAD_DEBUG_BLOCKKEEP > 0*/
00120 #endif /*RAD_DEBUG_BLOCKKEEP*/
00121 #endif /*RAD_AUTO_AD_Const*/
00122 
00123  class RAD_DoubleIgnore {};
00124 
00125  template<typename T> class
00126 DoubleAvoid {
00127  public:
00128   typedef double  dtype;
00129   typedef T ttype;
00130   };
00131  template<> class
00132 DoubleAvoid<double> {
00133  public:
00134   typedef RAD_DoubleIgnore &dtype;
00135   typedef RAD_DoubleIgnore &ttype;
00136   };
00137 
00138 #define Dtype typename DoubleAvoid<Double>::dtype
00139 #define Ttype typename DoubleAvoid<Double>::ttype
00140 
00141  template<typename Double> class IndepADvar;
00142  template<typename Double> class ConstADvar;
00143  template<typename Double> class ConstADvari;
00144  template<typename Double> class ADvar;
00145  template<typename Double> class ADvar1;
00146  template<typename Double> class ADvar1g;
00147  template<typename Double> class ADvar1s;
00148  template<typename Double> class ADvar2;
00149  template<typename Double> class ADvar2g;
00150  template<typename Double> class ADvar2q;
00151  template<typename Double> class ADvari;
00152  template<typename Double> class ADvari_block;
00153  template<typename Double> class ADvarn;
00154  template<typename Double> class Derp;
00155 
00156  template<typename Double> struct
00157 ADmemblock {  // We get memory in ADmemblock chunks and never give it back,
00158     // but reuse it once computations start anew after call(s) on
00159     // ADcontext::Gradcomp() or ADcontext::Weighted_Gradcomp().
00160   ADmemblock *next;
00161   Double memblk[2000];
00162   };
00163 
00164  template<typename Double> class
00165 ADvari_block {
00166  public:
00167   typedef ADvari<Double> ADVari;
00168   enum { Gulp = 1021 };
00169   ADvari_block *next, *prev;
00170   ADVari **limit;
00171   ADVari *pADvari[Gulp];
00172   };
00173 
00174  template<typename Double> class
00175 ADcontext { // A singleton class: one instance in radops.c
00176   typedef ADmemblock<Double> ADMemblock;
00177   typedef ADvari <Double> ADVari;
00178   typedef ADvar1 <Double> ADVar1;
00179   typedef ADvar2 <Double> ADVar2;
00180   typedef ADvarn <Double> ADVarn;
00181   typedef ADvar1g<Double> ADVar1g;
00182   typedef ADvar1s<Double> ADVar1s;
00183   typedef ADvar2g<Double> ADVar2g;
00184   typedef ADvar2q<Double> ADVar2q;
00185   typedef ADvari_block<Double> ADVari_block;
00186   ADMemblock *Busy, *Free;
00187   char *Mbase;
00188   size_t Mleft;
00189   ADVari **Ailimit, **Ainext;
00190   ADVari_block *Aibusy, *Aifree;
00191   ADMemblock *First;
00192   ADVari_block *AiFirst;
00193   double First0[(sizeof(ADMemblock) + sizeof(double) - 1) / sizeof(double)];
00194   double First1[(sizeof(ADVari_block) + sizeof(double) - 1) / sizeof(double)];
00195   void *new_ADmemblock(size_t);
00196   void new_ADvari_block();
00197   typedef ADvar<Double> ADVar;
00198   typedef Derp<Double> DErp;
00199   int rad_need_reinit;
00200   size_t rad_mleft_save;
00201 #ifdef RAD_DEBUG_BLOCKKEEP
00202   int rad_busy_blocks;
00203   ADMemblock *rad_Oldcurmb;
00204 #endif
00205  public:
00206   static const Double One, negOne;
00207   ADcontext();
00208   void *Memalloc(size_t len);
00209   static void Gradcomp(int wantgrad);
00210   static void Hvprod(int, ADVar**, Double*, Double*);
00211   static void aval_reset(void);
00212   static void Weighted_Gradcomp(int, ADVar**, Double*);
00213   static inline void Gradcomp() { Gradcomp(1); }
00214   inline void ADvari_record(ADVari *x) {
00215     if (Ainext >= Ailimit)
00216       new_ADvari_block();
00217     *Ainext++ = x;
00218     }
00219   };
00220 
00221  template<typename Double> class
00222 CADcontext: public ADcontext<Double> {  // for possibly constant ADvar values
00223  protected:
00224   bool fpval_implies_const;
00225  public:
00226   friend class ADvar<Double>;
00227   CADcontext(): ADcontext<Double>() { fpval_implies_const = false; }
00228   };
00229 
00230  template<typename Double> class
00231 Derp {    // one derivative-propagation operation
00232  public:
00233   friend class ADvarn<Double>;
00234   typedef ADvari<Double> ADVari;
00235   static Derp *LastDerp;
00236   Derp *next;
00237   const Double *a;
00238   const ADVari *b;
00239   mutable ADVari *c;
00240   Derp(){};
00241   Derp(const ADVari *);
00242   Derp(const Double *, const ADVari *);
00243   Derp(const Double *, const ADVari *, const ADVari *);
00244   /* c->aval += a * b->aval; */
00245   };
00246 
00247 
00248 // Now we use #define to overcome bad design in the C++ templating system
00249 
00250 #define Ai const ADvari<Double>&
00251 #define AI const IndepADvar<Double>&
00252 #define T template<typename Double>
00253 #define D Double
00254 #define T1(f) \
00255 T F f (AI); \
00256 T F f (Ai);
00257 #define T2(r,f) \
00258  T r f(Ai,Ai); \
00259  T r f(Ai,D); \
00260  T r f(Ai,Dtype); \
00261  T r f(Ai,long); \
00262  T r f(Ai,int); \
00263  T r f(D,Ai); \
00264  T r f(Dtype,Ai); \
00265  T r f(long,Ai); \
00266  T r f(int,Ai); \
00267  T r f(AI,D); \
00268  T r f(AI,Dtype); \
00269  T r f(AI,long); \
00270  T r f(AI,int); \
00271  T r f(D,AI); \
00272  T r f(Dtype,AI); \
00273  T r f(long,AI); \
00274  T r f(int,AI); \
00275  T r f(Ai,AI);\
00276  T r f(AI,Ai);\
00277  T r f(AI,AI);
00278 
00279 #define F ADvari<Double>&
00280 T2(F, operator+)
00281 T2(F, operator-)
00282 T2(F, operator*)
00283 T2(F, operator/)
00284 T2(F, atan2)
00285 T2(F, pow)
00286 T2(F, max)
00287 T2(F, min)
00288 T2(int, operator<)
00289 T2(int, operator<=)
00290 T2(int, operator==)
00291 T2(int, operator!=)
00292 T2(int, operator>=)
00293 T2(int, operator>)
00294 T1(operator+)
00295 T1(operator-)
00296 T1(abs)
00297 T1(acos)
00298 T1(acosh)
00299 T1(asin)
00300 T1(asinh)
00301 T1(atan)
00302 T1(atanh)
00303 T1(cos)
00304 T1(cosh)
00305 T1(exp)
00306 T1(log)
00307 T1(log10)
00308 T1(sin)
00309 T1(sinh)
00310 T1(sqrt)
00311 T1(tan)
00312 T1(tanh)
00313 T1(fabs)
00314 
00315 T F copy(AI);
00316 T F copy(Ai);
00317 
00318 #undef F
00319 #undef T2
00320 #undef T1
00321 #undef D
00322 #undef T
00323 #undef AI
00324 #undef Ai
00325 
00326 } /* namespace Rad2 */
00327 } /* namespace Sacado */
00328 #define SNS Sacado::Rad2
00329 namespace std { // Moved here from bottom for use in testing nesting of Rad with itself.
00330   using SNS::exp;
00331   using SNS::log;
00332   using SNS::log10;
00333   using SNS::sqrt;
00334   using SNS::cos;
00335   using SNS::sin;
00336   using SNS::tan;
00337   using SNS::acos;
00338   using SNS::asin;
00339   using SNS::atan;
00340   using SNS::cosh;
00341   using SNS::sinh;
00342   using SNS::tanh;
00343   using SNS::abs;
00344   using SNS::fabs;
00345   using SNS::atan2;
00346   using SNS::pow;
00347 }
00348 #undef SNS
00349 namespace Sacado {
00350 namespace Rad2 {
00351 
00352  template<typename Double>ADvari<Double>& ADf1(Double f, Double g, const IndepADvar<Double> &x);
00353  template<typename Double>ADvari<Double>& ADf2(Double f, Double gx, Double gy,
00354   const IndepADvar<Double> &x, const IndepADvar<Double> &y);
00355  template<typename Double>ADvari<Double>& ADfn(Double f, int n,
00356   const IndepADvar<Double> *x, const Double *g);
00357 
00358  template<typename Double> IndepADvar<Double>& ADvar_operatoreq(IndepADvar<Double>*,
00359   const ADvari<Double>&);
00360  template<typename Double> ADvar<Double>& ADvar_operatoreq(ADvar<Double>*, const ADvari<Double>&);
00361  template<typename Double> void AD_Const(const IndepADvar<Double>&);
00362  template<typename Double> void AD_Const1(Double*, const IndepADvar<Double>&);
00363  template<typename Double> ADvari<Double>& ADf1(Double, Double, const ADvari<Double>&);
00364  template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
00365   const ADvari<Double>&, const ADvari<Double>&);
00366  template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
00367   const IndepADvar<Double>&, const ADvari<Double>&);
00368  template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
00369   const ADvari<Double>&, const IndepADvar<Double>&);
00370  template<typename Double> Double val(const ADvari<Double>&);
00371 
00372  enum Advari_Opclass {
00373   Hv_const,
00374   Hv_copy,
00375   Hv_binary,
00376   Hv_unary,
00377   Hv_negate,
00378   Hv_plusLR,
00379   Hv_minusLR,
00380   Hv_timesL,
00381   Hv_timesLR,
00382   Hv_quotLR,
00383   Hv_nary
00384   };
00385 
00386  template<typename Double> ADvari<Double>&
00387 ADf1(Double f, Double g, Double h, const ADvari<Double> &x);
00388 
00389  template<typename Double> ADvari<Double>&
00390 ADf2(Double f, Double gx, Double gy, Double hxx,
00391       Double hxy, Double hyy, const ADvari<Double> &x, const ADvari<Double> &y);
00392 
00393  template<typename Double> ADvari<Double>&
00394 ADfn(Double f, int n, const IndepADvar<Double> *x, const Double *g, const Double *h);
00395 
00396  template<typename Double> class
00397 ADvari {  // implementation of an ADvar
00398  public:
00399   typedef Double value_type;
00400   typedef IndepADvar<Double> IndepADVar;
00401 #ifdef RAD_AUTO_AD_Const
00402   friend class IndepADvar<Double>;
00403 #ifdef RAD_Const_WARN
00404   mutable const IndepADVar *padv;
00405 #else
00406  protected:
00407   mutable const IndepADVar *padv;
00408 #endif //RAD_Const_WARN
00409  private:
00410   ADvari *Next;
00411   static ADvari *First_ADvari, **Last_ADvari;
00412  public:
00413 #endif //RAD_AUTO_AD_Const
00414 #ifdef RAD_DEBUG
00415   int gcgen;
00416   int opno;
00417   static int gcgen_cur, last_opno, zap_gcgen, zap_gcgen1, zap_opno;
00418   static FILE *debug_file;
00419 #endif
00420   static ADcontext<Double> adc;
00421   Advari_Opclass opclass;
00422   Double Val;   // result of this operation
00423   mutable Double aval;  // adjoint -- partial of final result w.r.t. this Val
00424   mutable Double dO;  // deriv of op w.r.t. t in x + t*p
00425   mutable Double aO;  // adjoint (in Hv computation) of op
00426   mutable Double adO; // adjoint (in Hv computation) of dO
00427   void *operator new(size_t len) {
00428 #ifdef RAD_DEBUG
00429     ADvari *rv = (ADvari*)ADvari::adc.Memalloc(len);
00430     rv->gcgen = gcgen_cur;
00431     rv->opno = ++last_opno;
00432     if (last_opno == zap_opno && gcgen_cur == zap_gcgen)
00433       printf("");
00434     return rv;
00435 #else
00436     return ADvari::adc.Memalloc(len);
00437 #endif
00438     }
00439   void operator delete(void*) {} /*Should never be called.*/
00440   inline ADvari(Advari_Opclass oc, Double t):
00441     opclass(oc), Val(t), aval(0.), dO(0.)
00442     { if (oc != Hv_const) ADvari::adc.ADvari_record(this); }
00443   inline ADvari(Advari_Opclass oc, Double t, Double ta):
00444     opclass(oc), Val(t), aval(ta), dO(0.)
00445     { if (oc != Hv_const) ADvari::adc.ADvari_record(this); }
00446  private:
00447   inline ADvari(): Val(0.), aval(0.), dO(0.) {} // prevent construction without value (?)
00448  public:
00449   friend class ConstADvari<Double>;
00450 #ifdef RAD_AUTO_AD_Const
00451   friend class ADcontext<Double>;
00452   friend class ADvar<Double>;
00453   friend class ADvar1<Double>;
00454   friend class ADvar1s<Double>;
00455   friend class ADvar2<Double>;
00456   friend class ADvar2q<Double>;
00457   friend class ConstADvar<Double>;
00458   ADvari(const IndepADVar *, Double);
00459   ADvari(const IndepADVar *, Double, Double);
00460   ADvari(const IndepADVar *, Double, Double, int);
00461 #endif
00462 #define F friend
00463 #define R ADvari&
00464 #define Ai const ADvari&
00465 #define T1(r,f) F r f <>(Ai);
00466 #define T2(r,f) \
00467 F r f <>(Ai,Ai); \
00468 F r f <>(Ttype,Ai); \
00469 F r f <>(Ai,Ttype); \
00470 F r f <>(double,Ai); \
00471 F r f <>(Ai,double); \
00472 F r f <>(long,Ai); \
00473 F r f <>(Ai,long); \
00474 F r f <>(int,Ai); \
00475 F r f <>(Ai,int);
00476   T1(R,operator+)
00477   T2(R,operator+)
00478   T1(R,operator-)
00479   T2(R,operator-)
00480   T2(R,operator*)
00481   T2(R,operator/)
00482   T1(R,abs)
00483   T1(R,acos)
00484   T1(R,acosh)
00485   T1(R,asin)
00486   T1(R,asinh)
00487   T1(R,atan)
00488   T1(R,atanh)
00489   T2(R,atan2)
00490   T2(R,max)
00491   T2(R,min)
00492   T1(R,cos)
00493   T1(R,cosh)
00494   T1(R,exp)
00495   T1(R,log)
00496   T1(R,log10)
00497   T2(R,pow)
00498   T1(R,sin)
00499   T1(R,sinh)
00500   T1(R,sqrt)
00501   T1(R,tan)
00502   T1(R,tanh)
00503   T1(R,fabs)
00504   T1(R,copy)
00505   T2(int,operator<)
00506   T2(int,operator<=)
00507   T2(int,operator==)
00508   T2(int,operator!=)
00509   T2(int,operator>=)
00510   T2(int,operator>)
00511 #undef T2
00512 #undef T1
00513 #undef Ai
00514 #undef R
00515 #undef F
00516 
00517   friend ADvari& ADf1<>(Double f, Double g, Double h, const ADvari &x);
00518   friend ADvari& ADf2<>(Double f, Double gx, Double gy, Double hxx,
00519       Double hxy, Double hyy, const ADvari &x, const ADvari &y);
00520   friend ADvari& ADfn<>(Double f, int n, const IndepADVar *x,
00521       const Double *g, const Double *h);
00522 
00523   inline operator Double() { return this->Val; }
00524   inline operator Double() const { return this->Val; }
00525   };
00526 
00527  template<typename Double> class
00528 ADvar1: public ADvari<Double> { // simplest unary ops
00529  public:
00530   typedef ADvari<Double> ADVari;
00531   Derp<Double> d;
00532   ADvar1(Advari_Opclass oc, Double val1, const Double *a1, const ADVari *c1):
00533     ADVari(oc,val1), d(a1,this,c1) {}
00534 #ifdef RAD_AUTO_AD_Const
00535   typedef typename ADVari::IndepADVar IndepADVar;
00536   typedef ADvar<Double> ADVar;
00537   ADvar1(const IndepADVar*, const IndepADVar&);
00538   ADvar1(const IndepADVar*, const ADVari&);
00539   ADvar1(const Double val1, const Double *a1, const ADVari *c1, const ADVar *v):
00540       ADVari(val1), d(a1,this,c1) {
00541     c1->padv = 0;
00542     *ADVari::Last_ADvari = this;
00543     ADVari::Last_ADvari = &this->Next;
00544     this->padv = v;
00545     }
00546 #endif
00547   };
00548 
00549 
00550  template<typename Double> class
00551 ConstADvari: public ADvari<Double> {
00552  private:
00553   ConstADvari *prevcad;
00554   ConstADvari() {}; // prevent construction without value (?)
00555   static ConstADvari *lastcad;
00556  public:
00557   typedef ADvari<Double> ADVari;
00558   typedef Derp<Double> DErp;
00559   static CADcontext<Double> cadc;
00560   inline void *operator new(size_t len) { return ConstADvari::cadc.Memalloc(len); }
00561   inline ConstADvari(Double t): ADVari(Hv_copy, t) { prevcad = lastcad; lastcad = this; }
00562   static void aval_reset(void);
00563   };
00564 
00565 
00566  template<typename Double> class
00567 IndepADvar {    // an independent ADvar
00568  private:
00569   IndepADvar& operator=(IndepADvar&x) {
00570     /* private to prevent assignment */
00571 #ifdef RAD_AUTO_AD_Const
00572     if (cv)
00573       cv->padv = 0;
00574     cv = new ADvar1<Double>(this,x);
00575     return *this;
00576 #else
00577 #ifdef RAD_EQ_ALIAS
00578     this->cv = x.cv;
00579     return *this;
00580 #else
00581     return ADvar_operatoreq(this,*x.cv);
00582 #endif
00583 #endif
00584     }
00585  protected:
00586   static void AD_Const(const IndepADvar&);
00587   mutable ADvari<Double> *cv;
00588  public:
00589   typedef Double value_type;
00590   friend class ADvar<Double>;
00591   friend class ADcontext<Double>;
00592   friend class ADvar1<Double>;
00593   friend class ADvarn<Double>;
00594   typedef ADvari<Double> ADVari;
00595   typedef ADvar<Double> ADVar;
00596   IndepADvar(Ttype);
00597   IndepADvar(double);
00598   IndepADvar(int);
00599   IndepADvar(long);
00600   IndepADvar& operator= (Double);
00601 #ifdef RAD_AUTO_AD_Const
00602   inline IndepADvar(const IndepADvar &x) { cv = x.cv ? new ADvar1<Double>(this, x) : 0; };
00603   inline IndepADvar() { cv = 0; }
00604   inline ~IndepADvar() {
00605           if (cv)
00606             cv->padv = 0;
00607           }
00608 #else
00609   inline IndepADvar() {
00610 #ifndef RAD_EQ_ALIAS
00611     cv = 0;
00612 #endif
00613     }
00614   inline ~IndepADvar() {}
00615   friend IndepADvar& ADvar_operatoreq<>(IndepADvar*, const ADVari&);
00616 #endif
00617 
00618 #ifdef RAD_Const_WARN
00619   inline operator ADVari&() const {
00620     ADVari *tcv = this->cv;
00621     if (tcv->opno < 0)
00622       RAD_Const_Warn(tcv);
00623     return *tcv;
00624     }
00625   inline operator ADVari*() const {
00626     ADVari *tcv = this->cv;
00627     if (tcv->opno < 0)
00628       RAD_Const_Warn(tcv);
00629     return tcv;
00630     }
00631 #else //RAD_Const_WARN
00632   inline operator ADVari&() const { return *this->cv; }
00633   inline operator ADVari*() const { return this->cv; }
00634 #endif //RAD_Const_WARN
00635 
00636   Double val() const { return cv->Val; }
00637   Double adj() const { return cv->aval; }
00638 
00639   friend void AD_Const1<>(Double*, const IndepADvar&);
00640 
00641   friend ADVari& ADf1<>(Double, Double, const IndepADvar&);
00642   friend ADVari& ADf2<>(Double, Double, Double, const IndepADvar&, const IndepADvar&);
00643   friend ADVari& ADf2<>(Double, Double, Double, const ADVari&, const IndepADvar&);
00644   friend ADVari& ADf2<>(Double, Double, Double, const IndepADvar&, const ADVari&);
00645 
00646   static inline void Gradcomp(int wantgrad)
00647         { ADcontext<Double>::Gradcomp(wantgrad); }
00648   static inline void Gradcomp()
00649         { ADcontext<Double>::Gradcomp(1); }
00650   static inline void Hvprod(int n, ADVar **vp, Double *v, Double *hv)
00651         { ADcontext<Double>::Hvprod(n, vp, v, hv); }
00652   static inline void aval_reset() { ConstADvari<Double>::aval_reset(); }
00653   static inline void Weighted_Gradcomp(int n, ADVar **v, Double *w)
00654         { ADcontext<Double>::Weighted_Gradcomp(n, v, w); }
00655 
00656   /* We use #define to deal with bizarre templating rules that apparently */
00657   /* require us to spell the some conversion explicitly */
00658 
00659 
00660 #define Ai const ADVari&
00661 #define AI const IndepADvar&
00662 #define D Double
00663 #define T2(r,f) \
00664  r f <>(AI,AI);\
00665  r f <>(Ai,AI);\
00666  r f <>(AI,Ai);\
00667  r f <>(Ttype,AI);\
00668  r f <>(double,AI);\
00669  r f <>(long,AI);\
00670  r f <>(int,AI);\
00671  r f <>(AI,Ttype);\
00672  r f <>(AI,double);\
00673  r f <>(AI,long);\
00674  r f <>(AI,int);
00675 #define T1(f) friend ADVari& f<> (AI);
00676 
00677 #define F friend ADVari&
00678 T2(F, operator+)
00679 T2(F, operator-)
00680 T2(F, operator*)
00681 T2(F, operator/)
00682 T2(F, atan2)
00683 T2(F, max)
00684 T2(F, min)
00685 T2(F, pow)
00686 #undef F
00687 #define F friend int
00688 T2(F, operator<)
00689 T2(F, operator<=)
00690 T2(F, operator==)
00691 T2(F, operator!=)
00692 T2(F, operator>=)
00693 T2(F, operator>)
00694 
00695 T1(operator+)
00696 T1(operator-)
00697 T1(abs)
00698 T1(acos)
00699 T1(acosh)
00700 T1(asin)
00701 T1(asinh)
00702 T1(atan)
00703 T1(atanh)
00704 T1(cos)
00705 T1(cosh)
00706 T1(exp)
00707 T1(log)
00708 T1(log10)
00709 T1(sin)
00710 T1(sinh)
00711 T1(sqrt)
00712 T1(tan)
00713 T1(tanh)
00714 T1(fabs)
00715 T1(copy)
00716 
00717 #undef F
00718 #undef T1
00719 #undef T2
00720 #undef D
00721 #undef AI
00722 #undef Ai
00723 
00724   };
00725 
00726  template<typename Double> class
00727 ADvar: public IndepADvar<Double> {  // an "active" variable
00728  public:
00730         template <typename U> struct apply { typedef ADvar<U> type; };
00731 
00732   typedef IndepADvar<Double> IndepADVar;
00733   typedef typename IndepADVar::ADVari ADVari;
00734   typedef ConstADvari<Double> ConstADVari;
00735  private:
00736   void ADvar_ctr(Double d) {
00737     ADVari *x = ConstADVari::cadc.fpval_implies_const
00738       ? new ConstADVari(d)
00739 #ifdef RAD_AUTO_AD_Const
00740       : new ADVari((IndepADVar*)this, d);
00741 #else
00742       : new ADVari(Hv_const, d);
00743 #endif
00744     this->cv = x;
00745     }
00746  public:
00747   friend class ADvar1<Double>;
00748   typedef ADvar1<Double> ADVar1;
00749   ADvar() { /* cv = 0; */ }
00750   ADvar(Ttype d)  { ADvar_ctr(d); }
00751   ADvar(double i) { ADvar_ctr(Double(i)); }
00752   ADvar(int i)  { ADvar_ctr(Double(i)); }
00753   ADvar(long i) { ADvar_ctr(Double(i)); }
00754   inline ~ADvar() {}
00755 #ifdef RAD_AUTO_AD_Const
00756   inline ADvar(IndepADVar &x) {
00757     this->cv = x.cv ? new ADVar1(this, x) : 0;
00758     }
00759   inline ADvar(ADVari &x) { this->cv = &x; x.padv = this; }
00760   inline ADvar& operator=(IndepADVar &x) {
00761     if (this->cv)
00762       this->cv->padv = 0;
00763     this->cv = new ADVar1(this,x);
00764     return *this;
00765     }
00766   inline ADvar& operator=(ADVari &x) {
00767     if (this->cv)
00768       this->cv->padv = 0;
00769     this->cv = new ADVar1(this, x);
00770     return *this;
00771     }
00772 #else 
00773   friend ADvar& ADvar_operatoreq<>(ADvar*, const ADVari&);
00774 #ifdef RAD_EQ_ALIAS
00775   /* allow aliasing v and w after "v = w;" */
00776   inline ADvar(const IndepADVar &x) { this->cv = (ADVari*)x.cv; }
00777   inline ADvar(const ADVari &x) { this->cv = (ADVari*)&x; }
00778   inline ADvar& operator=(IndepADVar &x) { this->cv = (ADVari*)x.cv; return *this; }
00779   inline ADvar& operator=(const ADVari &x) { this->cv = (ADVari*)&x; return *this; }
00780 #else 
00781   ADvar(const IndepADVar &x) { this->cv = x.cv ? new ADVar1(x.cv->Val, &this->cv->adc.One, x.cv) : 0; }
00782   ADvar(const ADvar&x) { this->cv = x.cv ?
00783     new ADVar1(Hv_copy, x.cv->Val, &this->cv->adc.One, (ADVari*)x.cv) : 0; }
00784   ADvar(const ADVari &x) { this->cv = new ADVar1(Hv_copy, x.Val, &this->cv->adc.One, &x); }
00785   inline ADvar& operator=(IndepADVar &x) { return ADvar_operatoreq(this,*x.cv); };
00786   inline ADvar& operator=(const ADVari &x) { return ADvar_operatoreq(this,x); };
00787 #endif /* RAD_EQ_ALIAS */
00788 #endif /* RAD_AUTO_AD_Const */
00789   ADvar& operator=(Double);
00790   ADvar& operator+=(const ADVari&);
00791   ADvar& operator+=(Double);
00792   ADvar& operator-=(const ADVari&);
00793   ADvar& operator-=(Double);
00794   ADvar& operator*=(const ADVari&);
00795   ADvar& operator*=(Double);
00796   ADvar& operator/=(const ADVari&);
00797   ADvar& operator/=(Double);
00798   inline static bool get_fpval_implies_const(void)
00799     { return ConstADVari::cadc.fpval_implies_const; }
00800   inline static void set_fpval_implies_const(bool newval)
00801     { ConstADVari::cadc.fpval_implies_const = newval; }
00802   inline static bool setget_fpval_implies_const(bool newval) {
00803     bool oldval = ConstADVari::cadc.fpval_implies_const;
00804     ConstADVari::cadc.fpval_implies_const = newval;
00805     return oldval;
00806     }
00807   static inline void Gradcomp(int wantgrad)
00808         { ADcontext<Double>::Gradcomp(wantgrad); }
00809   static inline void Gradcomp()
00810         { ADcontext<Double>::Gradcomp(1); }
00811   static inline void aval_reset() { ConstADVari::aval_reset(); }
00812   static inline void Weighted_Gradcomp(int n, ADvar **v, Double *w)
00813         { ADcontext<Double>::Weighted_Gradcomp(n, v, w); }
00814   };
00815 
00816 template<typename Double>
00817  inline void AD_Const1(Double *notused, const IndepADvar<Double>&v)
00818 { IndepADvar<Double>::AD_Const(v); }
00819 
00820 template<typename Double>
00821  inline void AD_Const(const IndepADvar<Double>&v) { AD_Const1((Double*)0, v); }
00822 
00823  template<typename Double> class
00824 ConstADvar: public ADvar<Double> {
00825  public:
00826   typedef ADvar<Double> ADVar;
00827   typedef typename ADVar::ADVari ADVari;
00828   typedef typename ADVar::ConstADVari ConstADVari;
00829   typedef Derp<Double> DErp;
00830   typedef typename ADVar::IndepADVar IndepADVar;
00831  private: // disable op=
00832   ConstADvar& operator+=(ADVari&);
00833   ConstADvar& operator+=(Double);
00834   ConstADvar& operator-=(ADVari&);
00835   ConstADvar& operator-=(Double);
00836   ConstADvar& operator*=(ADVari&);
00837   ConstADvar& operator*=(Double);
00838   ConstADvar& operator/=(ADVari&);
00839   ConstADvar& operator/=(Double);
00840   void ConstADvar_ctr(Double);
00841  public:
00842   ConstADvar(Ttype d) { ConstADvar_ctr(d); }
00843   ConstADvar(double i)  { ConstADvar_ctr(Double(i)); }
00844   ConstADvar(int i) { ConstADvar_ctr(Double(i)); }
00845   ConstADvar(long i)  { ConstADvar_ctr(Double(i)); }
00846   ConstADvar(const IndepADVar&);
00847   ConstADvar(const ConstADvar&);
00848   ConstADvar(const ADVari&);
00849   inline ~ConstADvar() {}
00850 #ifdef RAD_NO_CONST_UPDATE
00851  private:
00852 #endif
00853   ConstADvar();
00854   inline ConstADvar& operator=(Double d) { this->cv->Val = d; return *this; }
00855   inline ConstADvar& operator=(ADVari& d) { this->cv->Val = d.Val; return *this; }
00856  };
00857 
00858  template<typename Double> class
00859 ADvar1s: public ADvar1<Double> { // unary ops with partials
00860  public:
00861   typedef ADvar1<Double> ADVar1;
00862   typedef typename ADVar1::ADVari ADVari;
00863   Double pL;  // deriv of op w.r.t. left operand L
00864   ADvar1s(Double val1, Double a1, const ADVari *c1):
00865     ADVar1(Hv_timesL,val1,&pL,c1), pL(a1) {}
00866 #ifdef RAD_AUTO_AD_Const
00867   typedef typename ADVar1::ADVar ADVar;
00868   ADvar1s(Double val1, Double a1, const ADVari *c1, const ADVar *v):
00869     ADVar1(val1,&a,c1,v), a(a1) {}
00870 #endif
00871   };
00872 
00873  template<typename Double> class
00874 ADvar1g: public ADvar1<Double> { // unary ops with partials
00875  public:
00876   typedef ADvar1<Double> ADVar1;
00877   typedef typename ADVar1::ADVari ADVari;
00878   Double pL;  // deriv of op w.r.t. left operand L
00879   Double pL2; // partial of op w.r.t. L,L
00880   ADvar1g(Double val1, Double d1, Double d2, const ADVari *c1):
00881     ADVar1(Hv_unary,val1,&pL,c1), pL(d1), pL2(d2) {}
00882   };
00883 
00884  template<typename Double> class
00885 ADvar2: public ADvari<Double> { // basic binary ops
00886  public:
00887   typedef ADvari<Double> ADVari;
00888   typedef Derp<Double> DErp;
00889   DErp dL, dR;
00890   ADvar2(Advari_Opclass oc, Double val1, const ADVari *Lcv, const Double *Lc,
00891     const ADVari *Rcv, const Double *Rc):
00892       ADVari(oc,val1) {
00893     dR.next = DErp::LastDerp;
00894     dL.next = &dR;
00895     DErp::LastDerp = &dL;
00896     dL.a = Lc;
00897     dL.c = (ADVari*)Lcv;
00898     dR.a = Rc;
00899     dR.c = (ADVari*)Rcv;
00900     dL.b = dR.b = this;
00901     }
00902 #ifdef RAD_AUTO_AD_Const
00903   typedef ADvar<Double> ADVar;
00904   ADvar2(Double val1, const ADVari *Lcv, const Double *Lc,
00905     const ADVari *Rcv, const Double *Rc, ADVar *v):
00906       ADVari(val1) {
00907     dR.next = DErp::LastDerp;
00908     dL.next = &dR;
00909     DErp::LastDerp = &dL;
00910     dL.a = Lc;
00911     dL.c = Lcv;
00912     dR.a = Rc;
00913     dR.c = Rcv;
00914     dL.b = dR.b = this;
00915     Lcv->padv = 0;
00916     *ADVari::Last_ADvari = this;
00917     ADVari::Last_ADvari = &this->Next;
00918     this->padv = v;
00919     }
00920 #endif
00921   };
00922 
00923  template<typename Double> class
00924 ADvar2q: public ADvar2<Double> { // binary ops with partials
00925  public:
00926   typedef ADvar2<Double> ADVar2;
00927   typedef typename ADVar2::ADVari ADVari;
00928   typedef typename ADVar2::DErp DErp;
00929   Double pL;  // deriv of op w.r.t. left operand L
00930   Double pR;  // deriv of op w.r.t. right operand R
00931   Double pLR; // second partial w.r.t. L,R
00932   Double pR2; // second partial w.r.t. R,R
00933   ADvar2q(Double val1, Double Lp, Double Rp, Double LR, Double R2,
00934     const ADVari *Lcv, const ADVari *Rcv):
00935       ADVar2(Hv_quotLR,val1,Lcv,&pL,Rcv,&pR),
00936       pL(Lp), pR(Rp), pLR(LR), pR2(R2) {}
00937 #ifdef RAD_AUTO_AD_Const
00938   typedef typename ADVar2::ADVar ADVar;
00939   ADvar2q(Double val1, Double Lp, Double Rp, Double LR, Double R2,
00940     const ADVari *Lcv, const ADVari *Rcv, const ADVar *v):
00941       ADVar2(Hv_quotLR,val1,Lcv,&pL,Rcv,&pR,v),
00942       pL(Lp), pR(Rp), pLR(LR), pR2(R2) {}
00943 #endif
00944   };
00945 
00946  template<typename Double> class
00947 ADvar2g: public ADvar2<Double> { // general binary ops with partials
00948  public:
00949   typedef ADvar2<Double> ADVar2;
00950   typedef typename ADVar2::ADVari ADVari;
00951   Double pL;  // deriv of op w.r.t. left operand L
00952   Double pR;  // deriv of op w.r.t. right operand R
00953   Double pL2; // second partial w.r.t. L,L
00954   Double pLR; // second partial w.r.t. L,R
00955   Double pR2; // second partial w.r.t. R,R
00956   ADvar2g(Double val1, Double Lp, Double Rp, Double L2, Double LR, Double R2,
00957     const ADVari *Lcv, const ADVari *Rcv):
00958       ADVar2(Hv_binary,val1,Lcv,&pL,Rcv,&pR),
00959       pL(Lp), pR(Rp), pL2(L2), pLR(LR), pR2(R2) { }
00960   };
00961 
00962  template<typename Double> class
00963 ADvarn: public ADvari<Double> { // n-ary ops with partials g and
00964         // 2nd partials h (lower triangle, rowwise)
00965  public:
00966   typedef ADvari<Double> ADVari;
00967   typedef ADvar<Double> ADVar;
00968   typedef typename ADVari::IndepADVar IndepADVar;
00969   typedef Derp<Double> DErp;
00970   int n;
00971   Double *G, *H;
00972   DErp *D;
00973   ADvarn(Double val1, int n1, const IndepADVar *x, const Double *g, const Double *h):
00974       ADVari(Hv_nary,val1), n(n1) {
00975     DErp *d1, *dlast;
00976     Double *a1;
00977     int i, nh;
00978 
00979     a1 = G = (Double*)ADVari::adc.Memalloc(n1*sizeof(*g));
00980     d1 = D = (DErp*)ADVari::adc.Memalloc(n1*sizeof(DErp));
00981     dlast = DErp::LastDerp;
00982     for(i = 0; i < n1; i++, d1++) {
00983       d1->next = dlast;
00984       dlast = d1;
00985       a1[i] = g[i];
00986       d1->a = &a1[i];
00987       d1->b = this;
00988       d1->c = x[i].cv;
00989       }
00990     DErp::LastDerp = dlast;
00991     nh = (n1*(n1+1)) >> 1;
00992     a1 = H = (double*)ADVari::adc.Memalloc(nh * sizeof(*h));
00993     for(i = 0; i < nh; i++)
00994       a1[i] = h[i];
00995     }
00996   };
00997 
00998 template<typename Double>
00999  inline ADvari<Double>& operator+(const ADvari<Double> &T) { return *(ADvari<Double>*)&T; }
01000 
01001 template<typename Double>
01002  inline int operator<(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val < R.Val; }
01003 template<typename Double>
01004  inline int operator<(const ADvari<Double> &L, Double R) { return L.Val < R; }
01005 template<typename Double>
01006  inline int operator<(Double L, const ADvari<Double> &R) { return L < R.Val; }
01007 
01008 template<typename Double>
01009  inline int operator<=(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val <= R.Val; }
01010 template<typename Double>
01011  inline int operator<=(const ADvari<Double> &L, Double R) { return L.Val <= R; }
01012 template<typename Double>
01013  inline int operator<=(Double L, const ADvari<Double> &R) { return L <= R.Val; }
01014 
01015 template<typename Double>
01016  inline int operator==(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val == R.Val; }
01017 template<typename Double>
01018  inline int operator==(const ADvari<Double> &L, Double R) { return L.Val == R; }
01019 template<typename Double>
01020  inline int operator==(Double L, const ADvari<Double> &R) { return L == R.Val; }
01021 
01022 template<typename Double>
01023  inline int operator!=(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val != R.Val; }
01024 template<typename Double>
01025  inline int operator!=(const ADvari<Double> &L, Double R) { return L.Val != R; }
01026 template<typename Double>
01027  inline int operator!=(Double L, const ADvari<Double> &R) { return L != R.Val; }
01028 
01029 template<typename Double>
01030  inline int operator>=(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val >= R.Val; }
01031 template<typename Double>
01032  inline int operator>=(const ADvari<Double> &L, Double R) { return L.Val >= R; }
01033 template<typename Double>
01034  inline int operator>=(Double L, const ADvari<Double> &R) { return L >= R.Val; }
01035 
01036 template<typename Double>
01037  inline int operator>(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val > R.Val; }
01038 template<typename Double>
01039  inline int operator>(const ADvari<Double> &L, Double R) { return L.Val > R; }
01040 template<typename Double>
01041  inline int operator>(Double L, const ADvari<Double> &R) { return L > R.Val; }
01042 
01043 template<typename Double>
01044  inline void *ADcontext<Double>::Memalloc(size_t len) {
01045     if (Mleft >= len)
01046       return Mbase + (Mleft -= len);
01047     return new_ADmemblock(len);
01048     }
01049 
01050 template<typename Double>
01051  inline Derp<Double>::Derp(const ADVari *c1): c((ADVari*)c1) {
01052     next = LastDerp;
01053     LastDerp = this;
01054     }
01055 
01056 template<typename Double>
01057  inline Derp<Double>::Derp(const Double *a1, const ADVari *c1): a(a1), c((ADVari*)c1) {
01058     next = LastDerp;
01059     LastDerp = this;
01060     }
01061 
01062 template<typename Double>
01063  inline Derp<Double>::Derp(const Double *a1, const ADVari *b1, const ADVari *c1):
01064   a(a1), b(b1), c((ADVari*)c1) {
01065     next = LastDerp;
01066     LastDerp = this;
01067     }
01068 
01069 /**** radops ****/
01070 
01071 template<typename Double> Derp<Double> *Derp<Double>::LastDerp = 0;
01072 template<typename Double> ADcontext<Double> ADvari<Double>::adc;
01073 template<typename Double> const Double ADcontext<Double>::One = 1.;
01074 template<typename Double> const Double ADcontext<Double>::negOne = -1.;
01075 template<typename Double> CADcontext<Double> ConstADvari<Double>::cadc;
01076 template<typename Double> ConstADvari<Double> *ConstADvari<Double>::lastcad;
01077 
01078 #ifdef RAD_AUTO_AD_Const
01079 template<typename Double> ADvari<Double>*   ADvari<Double>::First_ADvari;
01080 template<typename Double> ADvari<Double>**  ADvari<Double>::Last_ADvari = &ADvari<Double>::First_ADvari;
01081 #endif
01082 
01083 #ifdef RAD_DEBUG
01084 #ifndef RAD_DEBUG_gcgen1
01085 #define RAD_DEBUG_gcgen1 -1
01086 #endif
01087 template<typename Double> int ADvari<Double>::gcgen_cur;
01088 template<typename Double> int ADvari<Double>::last_opno;
01089 template<typename Double> int ADvari<Double>::zap_gcgen;
01090 template<typename Double> int ADvari<Double>::zap_gcgen1 = RAD_DEBUG_gcgen1;
01091 template<typename Double> int ADvari<Double>::zap_opno;
01092 template<typename Double> FILE *ADvari<Double>::debug_file;
01093 #endif
01094 
01095 
01096  template<typename Double>
01097 ADcontext<Double>::ADcontext()
01098 {
01099   ADVari_block *fb;
01100 
01101   First = (ADMemblock*)First0;
01102   First->next = 0;
01103   Busy = First;
01104   Free = 0;
01105   Mbase = (char*)First->memblk;
01106   Mleft = sizeof(First->memblk);
01107   AiFirst = Aibusy = fb = (ADVari_block*)First1;
01108   Aifree = 0;
01109   Ainext = fb->pADvari;
01110   fb->next = fb->prev = 0;
01111   fb->limit = Ailimit = fb->pADvari + ADVari_block::Gulp;
01112   rad_need_reinit = 0;
01113 #ifdef RAD_DEBUG_BLOCKKEEP
01114   rad_busy_blocks = 0;
01115   rad_mleft_save = 0;
01116   rad_Oldcurmb = 0;
01117 #endif
01118   }
01119 
01120 template<typename Double> void*
01121 ADcontext<Double>::new_ADmemblock(size_t len)
01122 {
01123   ADMemblock *mb, *mb0, *mb1, *mbf, *x;
01124   ADVari_block *b;
01125 
01126   if ((rad_need_reinit & 1) && this == &ADVari::adc) {
01127     rad_need_reinit &= ~1;
01128     DErp::LastDerp = 0;
01129     Aibusy = b = AiFirst;
01130     Aifree = b->next;
01131     b->next = b->prev = 0;
01132     Ailimit = b->limit = (Ainext = b->pADvari) + ADVari_block::Gulp;
01133 #ifdef RAD_DEBUG_BLOCKKEEP
01134     Mleft = rad_mleft_save;
01135     if (Mleft < sizeof(First->memblk))
01136       _uninit_f2c(Mbase + Mleft,
01137         UninitType<Double>::utype,
01138         (sizeof(First->memblk) - Mleft)
01139         /sizeof(typename Sacado::ValueType<Double>::type));
01140     if ((mb = Busy->next)) {
01141       if (!(mb0 = rad_Oldcurmb))
01142         mb0 = (ADMemblock*)First0;
01143       for(;; mb = mb->next) {
01144         _uninit_f2c(mb->memblk,
01145           UninitType<Double>::utype,
01146           sizeof(First->memblk)
01147           /sizeof(typename Sacado::ValueType<Double>::type));
01148         if (mb == mb0)
01149           break;
01150         }
01151       }
01152     rad_Oldcurmb = Busy;
01153     if (rad_busy_blocks >= RAD_DEBUG_BLOCKKEEP) {
01154       rad_busy_blocks = 0;
01155       rad_Oldcurmb = 0;
01156       mb0 = (ADMemblock*)First0;
01157       mbf =  Free;
01158       for(mb = Busy; mb != mb0; mb = mb1) {
01159         mb1 = mb->next;
01160         mb->next = mbf;
01161         mbf = mb;
01162         }
01163       Free = mbf;
01164       Busy = mb;
01165       Mbase = (char*)First->memblk;
01166       Mleft = sizeof(First->memblk);
01167       }
01168 
01169 #else /* !RAD_DEBUG_BLOCKKEEP */
01170 
01171     mb0 = First;
01172     mbf =  Free;
01173     for(mb = Busy; mb != mb0; mb = mb1) {
01174       mb1 = mb->next;
01175       mb->next = mbf;
01176       mbf = mb;
01177       }
01178     Free = mbf;
01179     Busy = mb;
01180     Mbase = (char*)First->memblk;
01181     Mleft = sizeof(First->memblk);
01182 #endif /*RAD_DEBUG_BLOCKKEEP*/
01183     if (Mleft >= len)
01184       return Mbase + (Mleft -= len);
01185     }
01186 
01187   if ((x = Free))
01188     Free = x->next;
01189   else
01190     x = new ADMemblock;
01191 #ifdef RAD_DEBUG_BLOCKKEEP
01192   rad_busy_blocks++;
01193 #endif
01194   x->next = Busy;
01195   Busy = x;
01196   return (Mbase = (char*)x->memblk) +
01197     (Mleft = sizeof(First->memblk) - len);
01198   }
01199 
01200 template<typename Double> void
01201 ADcontext<Double>::new_ADvari_block()
01202 {
01203   ADVari_block *ob, *nb;
01204   ob = Aibusy;
01205   ob->limit = Ailimit;  // should be unnecessary, but harmless
01206   if ((nb = Aifree))
01207     Aifree = nb->next;
01208   else
01209     nb = new ADVari_block;
01210   Aibusy = Aibusy->next = nb;
01211   nb->limit = Ailimit = (Ainext = nb->pADvari) + ADVari_block::Gulp;
01212   ob->next = nb;
01213   nb->prev = ob;
01214   nb->next = 0;
01215   }
01216 
01217 template<typename Double> void
01218 ADcontext<Double>::Gradcomp(int wantgrad)
01219 {
01220   DErp *d;
01221 #ifdef RAD_AUTO_AD_Const
01222   ADVari *a, *anext;
01223   IndepADvar<Double> *v;
01224 #ifdef RAD_Const_WARN
01225   ADVari *cv;
01226   int i, j;
01227 #endif
01228 #endif /*RAD_AUTO_AD_Const*/
01229 
01230   if (ADVari::adc.rad_need_reinit) {
01231     for(d = DErp::LastDerp; d; d = d->next)
01232       d->c->aval = 0;
01233     }
01234   if (!(ADVari::adc.rad_need_reinit & 1)) {
01235     ADVari::adc.rad_need_reinit = 1;
01236     ADVari::adc.rad_mleft_save = ADVari::adc.Mleft;
01237     ADVari::adc.Mleft = 0;
01238     }
01239 #ifdef RAD_DEBUG
01240   if (ADVari::gcgen_cur == ADVari::zap_gcgen1) {
01241     const char *fname;
01242     if (!(fname = getenv("RAD_DEBUG_FILE")))
01243       fname = "rad_debug.out";
01244     else if (!*fname)
01245       fname = 0;
01246     if (fname)
01247       ADVari::debug_file = fopen(fname, "w");
01248     ADVari::zap_gcgen1 = -1;
01249     }
01250 #endif
01251   if ((d = DErp::LastDerp) && wantgrad) {
01252     d->b->aval = 1;
01253 #ifdef RAD_DEBUG
01254     if (ADVari::debug_file)
01255       do {
01256         fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
01257           d->c->opno, d->b->opno, d->c->aval, *d->a, d->b->aval);
01258         d->c->aval += *d->a * d->b->aval;
01259         fprintf(ADVari::debug_file, " = %g\n", d->c->aval);
01260         fflush(ADVari::debug_file);
01261         } while((d = d->next));
01262     else
01263 #endif
01264     do d->c->aval += *d->a * d->b->aval;
01265     while((d = d->next));
01266     }
01267 #ifdef RAD_DEBUG
01268   if (ADVari::debug_file) {
01269     fclose(ADVari::debug_file);
01270     ADVari::debug_file = 0;
01271     }
01272 #endif //RAD_DEBUG
01273 #ifdef RAD_AUTO_AD_Const
01274   *ADVari::Last_ADvari = 0;
01275   ADVari::Last_ADvari = &ADVari::First_ADvari;
01276   if ((anext = ADVari::First_ADvari) && !(ADVari::adc.rad_need_reinit & 2)) {
01277     ADVari::adc.rad_need_reinit = 3;
01278     while((a = anext)) {
01279       anext = a->Next;
01280       if ((v = (IndepADvar<Double> *)a->padv)) {
01281 #ifdef RAD_Const_WARN
01282         if ((i = a->opno) > 0)
01283           i = -i;
01284         j = a->gcgen;
01285         v->cv = cv = new ADVari(v, a->Val, a->aval);
01286         cv->opno = i;
01287         cv->gcgen = j;
01288 #else
01289         v->cv = new ADVari(v, a->Val, a->aval);
01290 #endif
01291         }
01292       }
01293     DErp::LastDerp = 0;
01294     }
01295 #endif /*RAD_AUTO_AD_Const*/
01296 #ifdef RAD_DEBUG
01297   ADVari::gcgen_cur++;
01298   ADVari::last_opno = 0;
01299 #endif
01300   }
01301 
01302 template<typename Double> void
01303 ADcontext<Double>::Weighted_Gradcomp(int n, ADVar **V, Double *w)
01304 {
01305   DErp *d;
01306   int i;
01307 #ifdef RAD_Const_WARN
01308   ADVari *cv;
01309   int j;
01310 #endif
01311 #ifdef RAD_AUTO_AD_Const
01312   ADVari *a, *anext;
01313   IndepADvar<Double> *v;
01314 #endif /*RAD_AUTO_AD_Const*/
01315 
01316   if (ADVari::adc.rad_need_reinit) {
01317     for(d = DErp::LastDerp; d; d = d->next)
01318       d->c->aval = 0;
01319     }
01320   if (!(ADVari::adc.rad_need_reinit & 1)) {
01321     ADVari::adc.rad_need_reinit = 1;
01322     ADVari::adc.rad_mleft_save = ADVari::adc.Mleft;
01323     ADVari::adc.Mleft = 0;
01324     }
01325 #ifdef RAD_DEBUG
01326   if (ADVari::gcgen_cur == ADVari::zap_gcgen1) {
01327     const char *fname;
01328     if (!(fname = getenv("RAD_DEBUG_FILE")))
01329       fname = "rad_debug.out";
01330     else if (!*fname)
01331       fname = 0;
01332     if (fname)
01333       ADVari::debug_file = fopen(fname, "w");
01334     ADVari::zap_gcgen1 = -1;
01335     }
01336 #endif
01337   if ((d = DErp::LastDerp)) {
01338     for(i = 0; i < n; i++)
01339       V[i]->cv->aval = w[i];
01340 #ifdef RAD_DEBUG
01341     if (ADVari::debug_file)
01342       do {
01343         fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
01344           d->c->opno, d->b->opno, d->c->aval, *d->a, d->b->aval);
01345         d->c->aval += *d->a * d->b->aval;
01346         fprintf(ADVari::debug_file, " = %g\n", d->c->aval);
01347         fflush(ADVari::debug_file);
01348         } while((d = d->next));
01349     else
01350 #endif
01351     do d->c->aval += *d->a * d->b->aval;
01352     while((d = d->next));
01353     }
01354 #ifdef RAD_DEBUG
01355   if (ADVari::debug_file) {
01356     fclose(ADVari::debug_file);
01357     ADVari::debug_file = 0;
01358     }
01359 #endif //RAD_DEBUG
01360 #ifdef RAD_AUTO_AD_Const
01361   *ADVari::Last_ADvari = 0;
01362   ADVari::Last_ADvari = &ADVari::First_ADvari;
01363   if ((anext = ADVari::First_ADvari) && !(ADVari::adc.rad_need_reinit & 2)) {
01364     ADVari::adc.rad_need_reinit = 3;
01365     while((a = anext)) {
01366       anext = a->Next;
01367       if ((v = (IndepADvar<Double> *)a->padv)) {
01368 #ifdef RAD_Const_WARN
01369         if ((i = a->opno) > 0)
01370           i = -i;
01371         j = a->gcgen;
01372         v->cv = cv = new ADVari(v, a->Val, a->aval);
01373         cv->opno = i;
01374         cv->gcgen = j;
01375 #else
01376         v->cv = new ADVari(v, a->Val, a->aval);
01377 #endif
01378         }
01379       }
01380     DErp::LastDerp = 0;
01381     }
01382 #endif /*RAD_AUTO_AD_Const*/
01383 #ifdef RAD_DEBUG
01384   ADVari::gcgen_cur++;
01385   ADVari::last_opno = 0;
01386 #endif
01387     }
01388 
01389  template<typename Double>
01390 IndepADvar<Double>::IndepADvar(Ttype d)
01391 {
01392 
01393   ADVari *x = new ADVari(Hv_const, d);
01394   cv = x;
01395   }
01396 
01397  template<typename Double>
01398 IndepADvar<Double>::IndepADvar(double i)
01399 {
01400 
01401   ADVari *x = new ADVari(Hv_const, Double(i));
01402   cv = x;
01403   }
01404 
01405  template<typename Double>
01406 IndepADvar<Double>::IndepADvar(int i)
01407 {
01408 
01409   ADVari *x = new ADVari(Hv_const, Double(i));
01410   cv = x;
01411   }
01412 
01413  template<typename Double>
01414 IndepADvar<Double>::IndepADvar(long i)
01415 {
01416 
01417   ADVari *x = new ADVari(Hv_const, Double(i));
01418   cv = x;
01419   }
01420 
01421  template<typename Double>
01422 ConstADvar<Double>::ConstADvar()
01423 {
01424   ConstADVari *x = new ConstADVari(0.);
01425   this->cv = x;
01426   }
01427 
01428  template<typename Double> void
01429 ConstADvar<Double>::ConstADvar_ctr(Double d)
01430 {
01431   ConstADVari *x = new ConstADVari(d);
01432   this->cv = x;
01433   }
01434 
01435  template<typename Double>
01436 ConstADvar<Double>::ConstADvar(const IndepADVar &x)
01437 {
01438   ConstADVari *y = new ConstADVari(x.cv->Val);
01439   DErp *d = new DErp(&x.adc.One, y, x.cv);
01440   this->cv = y;
01441   }
01442 
01443  template<typename Double>
01444 ConstADvar<Double>::ConstADvar(const ConstADvar &x)
01445 {
01446   ConstADVari *y = new ConstADVari(x.cv->Val);
01447   DErp *d = new DErp(&x.cv->adc.One, y, (ADVari*)x.cv);
01448   this->cv = y;
01449   }
01450 
01451  template<typename Double>
01452 ConstADvar<Double>::ConstADvar(const ADVari &x)
01453 {
01454   ConstADVari *y = new ConstADVari(x.Val);
01455   DErp *d = new DErp(&x.adc.One, y, &x);
01456   this->cv = y;
01457   }
01458 
01459  template<typename Double>
01460  void
01461 IndepADvar<Double>::AD_Const(const IndepADvar &v)
01462 {
01463   typedef ConstADvari<Double> ConstADVari;
01464 
01465   ConstADVari *ncv = new ConstADVari(v.val());
01466 #ifdef RAD_AUTO_AD_Const
01467   v.cv->padv = 0;
01468 #endif
01469   v.cv = ncv;
01470   }
01471 
01472  template<typename Double>
01473  void
01474 ConstADvari<Double>::aval_reset()
01475 {
01476   ConstADvari *x = ConstADvari::lastcad;
01477   while(x) {
01478     x->aval = 0;
01479     x = x->prevcad;
01480     }
01481   }
01482 
01483 #ifdef RAD_AUTO_AD_Const
01484 
01485  template<typename Double>
01486 ADvari<Double>::ADvari(const IndepADVar *x, Double d): Val(d), aval(0.)
01487 {
01488   opclass = Hv_const;
01489   *Last_ADvari = this;
01490   Last_ADvari = &Next;
01491   this->padv = (IndepADVar*)x;
01492   }
01493 
01494  template<typename Double>
01495 ADvari<Double>::ADvari(const IndepADVar *x, Double d, Double g): Val(d), aval(g)
01496 {
01497   *Last_ADvari = this;
01498   Last_ADvari = &Next;
01499   this->padv = (IndepADVar*)x;
01500   }
01501 
01502  template<typename Double>
01503 ADvar1<Double>::ADvar1(const IndepADVar *x, const IndepADVar &y):
01504   ADVari(Hv_copy, y.cv->Val), d((const Double*)&ADcontext<Double>::One, (ADVari*)this, y.cv)
01505 {
01506   *ADVari::Last_ADvari = this;
01507   ADVari::Last_ADvari = &this->Next;
01508   this->padv = x;
01509   }
01510 
01511  template<typename Double>
01512 ADvar1<Double>::ADvar1(const IndepADVar *x, const ADVari &y):
01513   ADVari(Hv_copy, y.Val), d((const Double*)&ADcontext<Double>::One, this, &y)
01514 {
01515   *ADVari::Last_ADvari = this;
01516   ADVari::Last_ADvari = &this->Next;
01517   this->padv = x;
01518   }
01519 
01520 #else /* !RAD_AUTO_AD_Const */
01521 
01522  template<typename Double>
01523  IndepADvar<Double>&
01524 ADvar_operatoreq(IndepADvar<Double> *This, const ADvari<Double> &x)
01525 { This->cv = new ADvar1<Double>(Hv_copy, x.Val, &x.adc.One, &x);
01526   return *(IndepADvar<Double>*) This; }
01527 
01528  template<typename Double>
01529  ADvar<Double>&
01530 ADvar_operatoreq(ADvar<Double> *This, const ADvari<Double> &x)
01531 { This->cv = new ADvar1<Double>(Hv_copy, x.Val, &x.adc.One, &x); return *(ADvar<Double>*) This; }
01532 
01533 #endif /* RAD_AUTO_AD_Const */
01534 
01535 
01536  template<typename Double>
01537  IndepADvar<Double>&
01538 IndepADvar<Double>::operator=(Double d)
01539 {
01540 #ifdef RAD_AUTO_AD_Const
01541   if (this->cv)
01542     this->cv->padv = 0;
01543   this->cv = new ADVari(Hv_const, this,d);
01544 #else
01545   this->cv = new ADVari(Hv_const, d);
01546 #endif
01547   return *this;
01548   }
01549 
01550  template<typename Double>
01551  ADvar<Double>&
01552 ADvar<Double>::operator=(Double d)
01553 {
01554 #ifdef RAD_AUTO_AD_Const
01555   if (this->cv)
01556     this->cv->padv = 0;
01557   this->cv = new ADVari(Hv_const, this,d);
01558 #else
01559   this->cv = ConstADVari::cadc.fpval_implies_const
01560     ? new ConstADVari(d)
01561     : new ADVari(Hv_const, d);
01562 #endif
01563   return *this;
01564   }
01565 
01566  template<typename Double>
01567  ADvari<Double>&
01568 operator-(const ADvari<Double> &T) {
01569   return *(new ADvar1<Double>(Hv_negate, -T.Val, &T.adc.negOne, &T));
01570   }
01571 
01572  template<typename Double>
01573  ADvari<Double>&
01574 operator+(const ADvari<Double> &L, const ADvari<Double> &R) {
01575   return *(new ADvar2<Double>(Hv_plusLR, L.Val + R.Val, &L, &L.adc.One, &R, &L.adc.One));
01576   }
01577 
01578 #ifdef RAD_AUTO_AD_Const
01579 #define RAD_ACA ,this
01580 #else
01581 #define RAD_ACA /*nothing*/
01582 #endif
01583 
01584  template<typename Double>
01585  ADvar<Double>&
01586 ADvar<Double>::operator+=(const ADVari &R) {
01587   ADVari *Lcv = this->cv;
01588   this->cv = new ADvar2<Double>(Hv_plusLR, Lcv->Val + R.Val, Lcv,
01589           &R.adc.One, &R, &R.adc.One RAD_ACA);
01590   return *this;
01591   }
01592 
01593  template<typename Double>
01594  ADvari<Double>&
01595 operator+(const ADvari<Double> &L, Double R) {
01596   return *(new ADvar1<Double>(Hv_copy, L.Val + R, &L.adc.One, &L));
01597   }
01598 
01599  template<typename Double>
01600  ADvar<Double>&
01601 ADvar<Double>::operator+=(Double R) {
01602   ADVari *tcv = this->cv;
01603   this->cv = new ADVar1(Hv_copy, tcv->Val + R, &tcv->adc.One, tcv RAD_ACA);
01604   return *this;
01605   }
01606 
01607  template<typename Double>
01608  ADvari<Double>&
01609 operator+(Double L, const ADvari<Double> &R) {
01610   return *(new ADvar1<Double>(Hv_copy, L + R.Val, &R.adc.One, &R));
01611   }
01612 
01613  template<typename Double>
01614  ADvari<Double>&
01615 operator-(const ADvari<Double> &L, const ADvari<Double> &R) {
01616   return *(new ADvar2<Double>(Hv_minusLR, L.Val - R.Val, &L, &L.adc.One, &R, &L.adc.negOne));
01617   }
01618 
01619  template<typename Double>
01620  ADvar<Double>&
01621 ADvar<Double>::operator-=(const ADVari &R) {
01622   ADVari *Lcv = this->cv;
01623   this->cv = new ADvar2<Double>(Hv_minusLR,Lcv->Val - R.Val, Lcv,
01624           &R.adc.One, &R, &R.adc.negOne RAD_ACA);
01625   return *this;
01626   }
01627 
01628  template<typename Double>
01629  ADvari<Double>&
01630 operator-(const ADvari<Double> &L, Double R) {
01631   return *(new ADvar1<Double>(Hv_copy, L.Val - R, &L.adc.One, &L));
01632   }
01633 
01634  template<typename Double>
01635  ADvar<Double>&
01636 ADvar<Double>::operator-=(Double R) {
01637   ADVari *tcv = this->cv;
01638   this->cv = new ADVar1(Hv_copy, tcv->Val - R, &tcv->adc.One, tcv RAD_ACA);
01639   return *this;
01640   }
01641 
01642  template<typename Double>
01643  ADvari<Double>&
01644 operator-(Double L, const ADvari<Double> &R) {
01645   return *(new ADvar1<Double>(Hv_negate, L - R.Val, &R.adc.negOne, &R));
01646   }
01647 
01648  template<typename Double>
01649  ADvari<Double>&
01650 operator*(const ADvari<Double> &L, const ADvari<Double> &R) {
01651   return *(new ADvar2<Double>(Hv_timesLR, L.Val * R.Val, &L, &R.Val, &R, &L.Val));
01652   }
01653 
01654  template<typename Double>
01655  ADvar<Double>&
01656 ADvar<Double>::operator*=(const ADVari &R) {
01657   ADVari *Lcv = this->cv;
01658   this->cv = new ADvar2<Double>(Hv_timesLR, Lcv->Val * R.Val, Lcv,
01659           &R.Val, &R, &Lcv->Val RAD_ACA);
01660   return *this;
01661   }
01662 
01663  template<typename Double>
01664  ADvari<Double>&
01665 operator*(const ADvari<Double> &L, Double R) {
01666   return *(new ADvar1s<Double>(L.Val * R, R, &L));
01667   }
01668 
01669  template<typename Double>
01670  ADvar<Double>&
01671 ADvar<Double>::operator*=(Double R) {
01672   ADVari *Lcv = this->cv;
01673   this->cv = new ADvar1s<Double>(Lcv->Val * R, R, Lcv RAD_ACA);
01674   return *this;
01675   }
01676 
01677  template<typename Double>
01678  ADvari<Double>&
01679 operator*(Double L, const ADvari<Double> &R) {
01680   return *(new ADvar1s<Double>(L * R.Val, L, &R));
01681   }
01682 
01683  template<typename Double>
01684  ADvari<Double>&
01685 operator/(const ADvari<Double> &L, const ADvari<Double> &R) {
01686   Double Lv = L.Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv, qpL = q*pL;
01687   return *(new ADvar2q<Double>(q, pL, -qpL, -pL*pL, 2.*pL*qpL, &L, &R));
01688   }
01689 
01690  template<typename Double>
01691  ADvar<Double>&
01692 ADvar<Double>::operator/=(const ADVari &R) {
01693   ADVari *Lcv = this->cv;
01694   Double Lv = Lcv->Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv, qpL = q*pL;
01695   this->cv = new ADvar2q<Double>(q, pL, -qpL, -pL*pL, 2.*pL*qpL, Lcv, &R RAD_ACA);
01696   return *this;
01697   }
01698 
01699  template<typename Double>
01700  ADvari<Double>&
01701 operator/(const ADvari<Double> &L, Double R) {
01702   return *(new ADvar1s<Double>(L.Val / R, 1./R, &L));
01703   }
01704 
01705  template<typename Double>
01706  ADvari<Double>&
01707 operator/(Double L, const ADvari<Double> &R) {
01708   Double recip = 1. / R.Val;
01709   Double q = L * recip;
01710   Double d1 = -q*recip;
01711   return *(new ADvar1g<Double>(q, d1, -q*d1, &R));
01712   }
01713 
01714  template<typename Double>
01715  ADvar<Double>&
01716 ADvar<Double>::operator/=(Double R) {
01717   ADVari *Lcv = this->cv;
01718   this->cv = new ADvar1s<Double>(Lcv->Val / R, 1./R, Lcv RAD_ACA);
01719   return *this;
01720   }
01721 
01722  template<typename Double>
01723  ADvari<Double>&
01724 acos(const ADvari<Double> &v) {
01725   Double t = v.Val, t1 = 1. - t*t, d1 = -1./std::sqrt(t1);
01726   return *(new ADvar1g<Double>(std::acos(t), d1, t*d1/t1, &v));
01727   }
01728 
01729  template<typename Double>
01730  ADvari<Double>&
01731 acosh(const ADvari<Double> &v) {
01732   Double d1, t, t1, t2;
01733   t = v.Val;
01734   t1 = std::sqrt(t2 = t*t - 1.);
01735   d1 = 1. / t1;
01736   return *(new ADvar1g<Double>(std::log(t + t1), d1, -t*d1/t2, &v));
01737   }
01738 
01739  template<typename Double>
01740  ADvari<Double>&
01741 asin(const ADvari<Double> &v) {
01742   Double d1, t, t1;
01743   t = v.Val;
01744   d1 = 1. / std::sqrt(t1 = 1. - t*t);
01745   return *(new ADvar1g<Double>(std::asin(t), d1, t*d1/t1, &v));
01746   }
01747 
01748  template<typename Double>
01749  ADvari<Double>&
01750 asinh(const ADvari<Double> &v) {
01751   Double d1, t, t1, t2, td;
01752   t = v.Val;
01753   t1 = std::sqrt(t2 = t*t + 1.);
01754   d1 = 1. / t1;
01755   td = 1.;
01756   if (t < 0.)
01757     td = -1.;
01758   return *(new ADvar1g<Double>(td*std::log(t*td + t1), d1, -(t/t2)*d1, &v));
01759   }
01760 
01761  template<typename Double>
01762  ADvari<Double>&
01763 atan(const ADvari<Double> &v) {
01764   Double t = v.Val, d1 = 1./(1. + t*t);
01765   return *(new ADvar1g<Double>(std::atan(t), d1, -(t+t)*d1*d1, &v));
01766   }
01767 
01768  template<typename Double>
01769  ADvari<Double>&
01770 atanh(const ADvari<Double> &v) {
01771   Double t = v.Val, d1 = 1./(1. - t*t);
01772   return *(new ADvar1g<Double>(0.5*std::log((1.+t)/(1.-t)), d1, (t+t)*d1*d1, &v));
01773   }
01774 
01775  template<typename Double>
01776  ADvari<Double>&
01777 atan2(const ADvari<Double> &L, const ADvari<Double> &R) {
01778   Double R2, t, t2, x, x2, y, y2;
01779   x = L.Val;
01780   y = R.Val;
01781   x2 = x*x;
01782   y2 = y*y;
01783   t = 1./(x2 + y2);
01784   t2 = t*t;
01785   R2 = 2.*t2*x*y;
01786   return *(new ADvar2g<Double>(std::atan2(x,y), y*t, -x*t, -R2, t2*(x2 - y2), R2, &L, &R));
01787   }
01788 
01789  template<typename Double>
01790  ADvari<Double>&
01791 atan2(Double x, const ADvari<Double> &R) {
01792   Double t, x2, y, y2;
01793   y = R.Val;
01794   x2 = x*x;
01795   y2 = y*y;
01796   t = 1./(x2 + y2);
01797   return *(new ADvar1g<Double>(std::atan2(x,y), -x*t, 2.*t*t*x*y, &R));
01798   }
01799 
01800  template<typename Double>
01801  ADvari<Double>&
01802 atan2(const ADvari<Double> &L, Double y) {
01803   Double t, x, x2, y2;
01804   x = L.Val;
01805   x2 = x*x;
01806   y2 = y*y;
01807   t = 1./(x2 + y2);
01808   return *(new ADvar1g<Double>(std::atan2(x,y), y*t, -2.*t*t*x*y, &L));
01809   }
01810 
01811  template<typename Double>
01812  ADvari<Double>&
01813 max(const ADvari<Double> &L, const ADvari<Double> &R) {
01814   const ADvari<Double> &x = L.Val >= R.Val ? L : R;
01815   return *(new ADvar1<Double>(Hv_copy, x.Val, &x.adc.One, &x));
01816   }
01817 
01818  template<typename Double>
01819  ADvari<Double>&
01820 max(Double L, const ADvari<Double> &R) {
01821   if (L >= R.Val)
01822     return *(new ADvari<Double>(Hv_const, L));
01823   return *(new ADvar1<Double>(Hv_copy, R.Val, &R.adc.One, &R));
01824   }
01825 
01826  template<typename Double>
01827  ADvari<Double>&
01828 max(const ADvari<Double> &L, Double R) {
01829   if (L.Val >= R)
01830     return *(new ADvar1<Double>(Hv_copy, L.Val, &L.adc.One, &L));
01831   return *(new ADvari<Double>(Hv_const, R));
01832   }
01833 
01834  template<typename Double>
01835  ADvari<Double>&
01836 min(const ADvari<Double> &L, const ADvari<Double> &R) {
01837   const ADvari<Double> &x = L.Val <= R.Val ? L : R;
01838   return *(new ADvar1<Double>(Hv_copy, x.Val, &x.adc.One, &x));
01839   }
01840 
01841  template<typename Double>
01842  ADvari<Double>&
01843 min(Double L, const ADvari<Double> &R) {
01844   if (L <= R.Val)
01845     return *(new ADvari<Double>(Hv_const, L));
01846   return *(new ADvar1<Double>(Hv_copy, R.Val, &R.adc.One, &R));
01847   }
01848 
01849  template<typename Double>
01850  ADvari<Double>&
01851 min(const ADvari<Double> &L, Double R) {
01852   if (L.Val <= R)
01853     return *(new ADvar1<Double>(Hv_copy, L.Val, &L.adc.One, &L));
01854   return *(new ADvari<Double>(Hv_const, R));
01855   }
01856 
01857  template<typename Double>
01858  ADvari<Double>&
01859 cos(const ADvari<Double> &v) {
01860   Double t = std::cos(v.Val);
01861   return *(new ADvar1g<Double>(t, -std::sin(v.Val), -t, &v));
01862   }
01863 
01864  template<typename Double>
01865  ADvari<Double>&
01866 cosh(const ADvari<Double> &v) {
01867   Double t = std::cosh(v.Val);
01868   return *(new ADvar1g<Double>(t, std::sinh(v.Val), t, &v));
01869   }
01870 
01871  template<typename Double>
01872  ADvari<Double>&
01873 exp(const ADvari<Double> &v) {
01874   Double t = std::exp(v.Val);
01875   return *(new ADvar1g<Double>(t, t, t, &v));
01876   }
01877 
01878  template<typename Double>
01879  ADvari<Double>&
01880 log(const ADvari<Double> &v) {
01881   Double x = v.Val, d1 = 1. / x;
01882   return *(new ADvar1g<Double>(std::log(x), d1, -d1*d1, &v));
01883   }
01884 
01885  template<typename Double>
01886  ADvari<Double>&
01887 log10(const ADvari<Double> &v) {
01888   static double num = 1. / std::log(10.);
01889   Double d1, t, x;
01890   x = v.Val;
01891   t = 1. / x;
01892   d1 = num * t;
01893   return *(new ADvar1g<Double>(std::log10(x), d1, -d1*t, &v));
01894   }
01895 
01896  template<typename Double>
01897  ADvari<Double>&
01898 pow(const ADvari<Double> &L, const ADvari<Double> &R) {
01899   Double dx, dy, t, x, xlog, xym1, y;
01900   x = L.Val;
01901   y = R.Val;
01902   t = std::pow(x,y);
01903   xym1 = t / x;
01904   xlog = std::log(x);
01905   dx = y*xym1;
01906   dy = t * xlog;
01907   return *(new ADvar2g<Double>(t, dx, dy, (y-1.)*dx/x, xym1*(1. + y*xlog), dy*xlog, &L, &R));
01908   }
01909 
01910  template<typename Double>
01911  ADvari<Double>&
01912 pow(Double x, const ADvari<Double> &R) {
01913   Double dy, t, xlog, y;
01914   y = R.Val;
01915   t = std::pow(x,y);
01916   xlog = std::log(x);
01917   dy = t * xlog;
01918   return *(new ADvar1g<Double>(t, dy, dy*xlog, &R));
01919   }
01920 
01921  template<typename Double>
01922  ADvari<Double>&
01923 pow(const ADvari<Double> &L, Double y) {
01924   Double dx, t, x;
01925   x = L.Val;
01926   t = std::pow(x,y);
01927   dx = y*t/x;
01928   return *(new ADvar1g<Double>(t, dx, (y-1.)*dx/x, &L));
01929   }
01930 
01931  template<typename Double>
01932  ADvari<Double>&
01933 sin(const ADvari<Double> &v) {
01934   Double t = std::sin(v.Val);
01935   return *(new ADvar1g<Double>(t, std::cos(v.Val), -t, &v));
01936   }
01937 
01938  template<typename Double>
01939  ADvari<Double>&
01940 sinh(const ADvari<Double> &v) {
01941   Double t = std::sinh(v.Val);
01942   return *(new ADvar1g<Double>(t, std::cosh(v.Val), t, &v));
01943   }
01944 
01945  template<typename Double>
01946  ADvari<Double>&
01947 sqrt(const ADvari<Double> &v) {
01948   Double t = std::sqrt(v.Val);
01949   Double d1 = 0.5 / t;
01950   return *(new ADvar1g<Double>(t, d1, -0.5*d1/v.Val, &v));
01951   }
01952 
01953  template<typename Double>
01954  ADvari<Double>&
01955 tan(const ADvari<Double> &v) {
01956   Double d1, rv, t;
01957   rv = std::tan(v.Val);
01958   t = 1. / std::cos(v.Val);
01959   d1 = t*t;
01960   return *(new ADvar1g<Double>(rv, d1, (rv+rv)*d1, &v));
01961   }
01962 
01963  template<typename Double>
01964  ADvari<Double>&
01965 tanh(const ADvari<Double> &v) {
01966   Double d1, rv, t;
01967   rv = std::tanh(v.Val);
01968   t = 1. / std::cosh(v.Val);
01969   d1 = t*t;
01970   return *(new ADvar1g<Double>(rv, d1, -(rv+rv)*d1, &v));
01971   }
01972 
01973  template<typename Double>
01974  ADvari<Double>&
01975 abs(const ADvari<Double> &v) {
01976   Double t, p;
01977   p = 1.;
01978   if ((t = v.Val) < 0) {
01979     t = -t;
01980     p = -p;
01981     }
01982   return *(new ADvar1g<Double>(t, p, 0., &v));
01983   }
01984 
01985  template<typename Double>
01986  ADvari<Double>&
01987 fabs(const ADvari<Double> &v) { // Synonym for "abs"
01988         // "fabs" is not the best choice of name,
01989         // but this name is used at Sandia.
01990   Double t, p;
01991   p = 1.;
01992   if ((t = v.Val) < 0) {
01993     t = -t;
01994     p = -p;
01995     }
01996   return *(new ADvar1g<Double>(t, p, 0., &v));
01997   }
01998 
01999  template<typename Double>
02000  ADvari<Double>&
02001 ADf1(Double f, Double g, Double h, const ADvari<Double> &x) {
02002   return *(new ADvar1g<Double>(f, g, h, &x));
02003   }
02004 
02005  template<typename Double>
02006  inline ADvari<Double>&
02007 ADf1(Double f, Double g, Double h, const IndepADvar<Double> &x) {
02008   return *(new ADvar1g<Double>(f, g, h, x.cv));
02009   }
02010 
02011  template<typename Double>
02012  ADvari<Double>&
02013 ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
02014     const ADvari<Double> &x, const ADvari<Double> &y) {
02015   return *(new ADvar2g<Double>(f, gx, gy, hxx, hxy, hyy, &x, &y));
02016   }
02017 
02018  template<typename Double>
02019  ADvari<Double>&
02020 ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
02021     const ADvari<Double> &x, const IndepADvar<Double> &y) {
02022   return *(new ADvar2g<Double>(f, gx, gy, hxx, hxy, hyy, &x, y.cv));
02023   }
02024 
02025  template<typename Double>
02026  ADvari<Double>&
02027 ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
02028     const IndepADvar<Double> &x, const ADvari<Double> &y) {
02029   return *(new ADvar2g<Double>(f, gx, gy, hxx, hxy, hyy, x.cv, &y));
02030   }
02031 
02032  template<typename Double>
02033  ADvari<Double>&
02034 ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
02035     const IndepADvar<Double> &x, const IndepADvar<Double> &y) {
02036   return *(new ADvar2g<Double>(f, gx, gy, hxx, hxy, hyy, x.cv, y.cv));
02037   }
02038 
02039  template<typename Double>
02040  ADvari<Double>&
02041 ADfn(Double f, int n, const IndepADvar<Double> *x, const Double *g, const Double *h) {
02042   return *(new ADvarn<Double>(f, n, x, g, h));
02043   }
02044 
02045  template<typename Double>
02046  inline ADvari<Double>&
02047 ADfn(Double f, int n, const ADvar<Double> *x, const Double *g, const Double *h) {
02048   return ADfn<Double>(f, n, (IndepADvar<Double>*)x, g, h);
02049   }
02050 
02051  template<typename Double>
02052  void
02053 ADcontext<Double>::Hvprod(int n, ADvar<Double> **x, Double *v, Double *hv)
02054 {
02055   ADVari *a, *aL, *aR, **ap, **ape;
02056   ADVari_block *b, *b0;
02057   DErp *d;
02058   Double aO, adO, *g, *h, *h0, t, tL, tR;
02059   int i, j, k, m;
02060   for(i = 0; i < n; i++) {
02061     a = x[i]->cv;
02062     a->dO = v[i];
02063     a->aO = a->adO = 0.;
02064     }
02065   ADVari::adc.Aibusy->limit = ADVari::adc.Ainext;
02066   a = 0;
02067   for(b0 = 0, b = ADVari::adc.AiFirst; b; b0 = b, b = b->next) {
02068     ap = b->pADvari;
02069     ape = b->limit;
02070     while(ap < ape) {
02071       a = *ap++;
02072       a->aO = a->adO = 0.;
02073       switch(a->opclass) {
02074        case Hv_copy:
02075         a->dO = ((ADVar1*)a)->d.c->dO;
02076         break;
02077        case Hv_binary:
02078         a->dO =   ((ADVar2g*)a)->pL * ((ADVar2g*)a)->dL.c->dO
02079           + ((ADVar2g*)a)->pR * ((ADVar2g*)a)->dR.c->dO;
02080         break;
02081        case Hv_unary:
02082         a->dO = ((ADVar1g*)a)->pL * ((ADVar1g*)a)->d.c->dO;
02083         break;
02084        case Hv_negate:
02085         a->dO = -((ADVar1*)a)->d.c->dO;
02086         break;
02087        case Hv_plusLR:
02088         a->dO = ((ADVar2*)a)->dL.c->dO + ((ADVar2*)a)->dR.c->dO;
02089         break;
02090        case Hv_minusLR:
02091         a->dO = ((ADVar2*)a)->dL.c->dO - ((ADVar2*)a)->dR.c->dO;
02092         break;
02093        case Hv_timesL:
02094         a->dO = ((ADVar1s*)a)->pL * ((ADVar1s*)a)->d.c->dO;
02095         break;
02096        case Hv_timesLR:
02097         a->dO =   ((ADVar2*)a)->dR.c->Val * ((ADVar2*)a)->dL.c->dO
02098           + ((ADVar2*)a)->dL.c->Val * ((ADVar2*)a)->dR.c->dO;
02099         break;
02100        case Hv_quotLR:
02101         a->dO =   ((ADVar2q*)a)->pL * ((ADVar2q*)a)->dL.c->dO
02102           + ((ADVar2q*)a)->pR * ((ADVar2q*)a)->dR.c->dO;
02103         break;
02104        case Hv_nary:
02105         d = ((ADVarn*)a)->D;
02106         m = ((ADVarn*)a)->n;
02107         g = ((ADVarn*)a)->G;
02108         t = 0.;
02109         for(i = 0; i < m; i++)
02110           t += g[i] * d[i].c->dO;
02111         a->dO = t;
02112        }
02113       }
02114     }
02115   if (a)
02116     a->adO = 1.;
02117   for(b = b0; b; b = b->prev) {
02118     ape = b->pADvari;
02119     ap = b->limit;
02120     while(ap > ape) {
02121       a = *--ap;
02122       aO = a->aO;
02123       adO = a->adO;
02124       switch(a->opclass) {
02125        case Hv_copy:
02126         aL = ((ADVar1*)a)->d.c;
02127         aL->aO += aO;
02128         aL->adO += adO;
02129         break;
02130        case Hv_binary:
02131         aL = ((ADVar2g*)a)->dL.c;
02132         aR = ((ADVar2g*)a)->dR.c;
02133         tL = adO*aL->dO;
02134         tR = adO*aR->dO;
02135         aL->aO += aO*((ADVar2g*)a)->pL
02136           + tL*((ADVar2g*)a)->pL2
02137           + tR*((ADVar2g*)a)->pLR;
02138         aR->aO += aO*((ADVar2g*)a)->pR
02139           + tL*((ADVar2g*)a)->pLR
02140           + tR*((ADVar2g*)a)->pR2;
02141         aL->adO += adO * ((ADVar2g*)a)->pL;
02142         aR->adO += adO * ((ADVar2g*)a)->pR;
02143         break;
02144        case Hv_unary:
02145         aL = ((ADVar1g*)a)->d.c;
02146         aL->aO += aO * ((ADVar1g*)a)->pL
02147           + adO * aL->dO * ((ADVar1g*)a)->pL2;
02148         aL->adO += adO * ((ADVar1g*)a)->pL;
02149         break;
02150        case Hv_negate:
02151         aL = ((ADVar1*)a)->d.c;
02152         aL->aO -= aO;
02153         aL->adO -= adO;
02154         break;
02155        case Hv_plusLR:
02156         aL = ((ADVar2*)a)->dL.c;
02157         aR = ((ADVar2*)a)->dR.c;
02158         aL->aO += aO;
02159         aL->adO += adO;
02160         aR->aO += aO;
02161         aR->adO += adO;
02162         break;
02163        case Hv_minusLR:
02164         aL = ((ADVar2*)a)->dL.c;
02165         aR = ((ADVar2*)a)->dR.c;
02166         aL->aO += aO;
02167         aL->adO += adO;
02168         aR->aO -= aO;
02169         aR->adO -= adO;
02170         break;
02171        case Hv_timesL:
02172         aL = ((ADVar1s*)a)->d.c;
02173         aL->aO += aO * (tL = ((ADVar1s*)a)->pL);
02174         aL->adO += adO * tL;
02175         break;
02176        case Hv_timesLR:
02177         aL = ((ADVar2*)a)->dL.c;
02178         aR = ((ADVar2*)a)->dR.c;
02179         aL->aO += aO * (tL = aR->Val) + adO*aR->dO;
02180         aR->aO += aO * (tR = aL->Val) + adO*aL->dO;
02181         aL->adO += adO * tL;
02182         aR->adO += adO * tR;
02183         break;
02184        case Hv_quotLR:
02185         aL = ((ADVar2q*)a)->dL.c;
02186         aR = ((ADVar2q*)a)->dR.c;
02187         tL = adO*aL->dO;
02188         tR = adO*aR->dO;
02189         aL->aO += aO*((ADVar2q*)a)->pL
02190           + tR*((ADVar2q*)a)->pLR;
02191         aR->aO += aO*((ADVar2q*)a)->pR
02192           + tL*((ADVar2q*)a)->pLR
02193           + tR*((ADVar2q*)a)->pR2;
02194         aL->adO += adO * ((ADVar2q*)a)->pL;
02195         aR->adO += adO * ((ADVar2q*)a)->pR;
02196         break;
02197        case Hv_nary:
02198         d  = ((ADVarn*)a)->D;
02199         m  = ((ADVarn*)a)->n;
02200         g  = ((ADVarn*)a)->G;
02201         h0 = ((ADVarn*)a)->H;
02202         for(i = 0; i < m; i++) {
02203           aL = d[i].c;
02204           aL->adO += adO * (t = g[i]);
02205           aL->aO += t*aO;
02206           t = adO * aL->dO;
02207           for(h = h0, j = 0; j <= i; j++)
02208             d[j].c->aO += t * *h++;
02209           h0 = h--;
02210           for(k = j; j < m; j++)
02211             d[j].c->aO += t * *(h += k++);
02212           }
02213        }
02214       }
02215     }
02216   for(i = 0; i < n; i++) {
02217     a = x[i]->cv;
02218     a->dO = 0.;
02219     hv[i] = a->aO;
02220     }
02221   }
02222 
02223  template<typename Double>
02224  inline Double
02225 val(const ADvari<Double> &x) {
02226   return x.Val;
02227   }
02228 
02229 #undef RAD_ACA
02230 #define A (ADvari<Double>*)
02231 #ifdef RAD_Const_WARN
02232 #define C(x) (((x)->opno < 0) ? RAD_Const_Warn(x) : 0, *A x)
02233 #else
02234 #define C(x) *A x
02235 #endif
02236 #define T template<typename Double> inline
02237 #define F ADvari<Double>&
02238 #define Ai const ADvari<Double>&
02239 #define AI const IndepADvar<Double>&
02240 #define D Double
02241 #define T2(r,f) \
02242  T r f(Ai L, AI R) { return f(L, C(R.cv)); }\
02243  T r f(AI L, Ai R) { return f(C(L.cv), R); }\
02244  T r f(AI L, AI R) { return f(C(L.cv), C(R.cv)); }\
02245  T r f(AI L, D R) { return f(C(L.cv), R); }\
02246  T r f(Ai L, Dtype R) { return f(L, (D)R); }\
02247  T r f(AI L, Dtype R) { return f(C(L.cv), (D)R); }\
02248  T r f(Ai L, long R) { return f(L, (D)R); }\
02249  T r f(AI L, long R) { return f(C(L.cv), (D)R); }\
02250  T r f(Ai L, int R) { return f(L, (D)R); }\
02251  T r f(AI L, int R) { return f(C(L.cv), (D)R); }\
02252  T r f(D L, AI R) { return f(L, C(R.cv)); }\
02253  T r f(Dtype L, Ai R) { return f((D)L, R); }\
02254  T r f(Dtype L, AI R) { return f((D)L, C(R.cv)); }\
02255  T r f(long L, Ai R) { return f((D)L, R); }\
02256  T r f(long L, AI R) { return f((D)L, C(R.cv)); }\
02257  T r f(int L, Ai R) { return f((D)L, R); }\
02258  T r f(int L, AI R) { return f((D)L, C(R.cv)); }
02259 
02260 T2(F, operator+)
02261 T2(F, operator-)
02262 T2(F, operator*)
02263 T2(F, operator/)
02264 T2(F, atan2)
02265 T2(F, pow)
02266 T2(F, max)
02267 T2(F, min)
02268 T2(int, operator<)
02269 T2(int, operator<=)
02270 T2(int, operator==)
02271 T2(int, operator!=)
02272 T2(int, operator>=)
02273 T2(int, operator>)
02274 
02275 #undef T2
02276 #undef D
02277 
02278 #define T1(f)\
02279  T F f(AI x) { return f(C(x.cv)); }
02280 
02281 T1(operator+)
02282 T1(operator-)
02283 T1(abs)
02284 T1(acos)
02285 T1(acosh)
02286 T1(asin)
02287 T1(asinh)
02288 T1(atan)
02289 T1(atanh)
02290 T1(cos)
02291 T1(cosh)
02292 T1(exp)
02293 T1(log)
02294 T1(log10)
02295 T1(sin)
02296 T1(sinh)
02297 T1(sqrt)
02298 T1(tan)
02299 T1(tanh)
02300 T1(fabs)
02301 
02302 T F copy(AI x)
02303 { return *(new ADvar1<Double>(Hv_copy, x.cv->Val, &ADcontext<Double>::One, (ADvari<Double>*)x.cv)); }
02304 
02305 T F copy(Ai x)
02306 { return *(new ADvar1<Double>(Hv_copy, x.Val, &ADcontext<Double>::One, (ADvari<Double>*)&x)); }
02307 
02308 #undef T1
02309 #undef AI
02310 #undef Ai
02311 #undef F
02312 #undef T
02313 #undef A
02314 #undef C
02315 #undef Ttype
02316 #undef Dtype
02317 
02318 } /* namespace Rad2 */
02319 } /* namespace Sacado */
02320 
02321 #endif /* SACADO_TRAD2_H */

Generated on Wed May 12 21:39:38 2010 for Sacado Package Browser (Single Doxygen Collection) by  doxygen 1.4.7