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 //!RAD_NO_USING_STDCC
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   inline void ADvari_padv(const IndepADVar *v) {
00414     *Last_ADvari = this;
00415     Last_ADvari = &this->Next;
00416     this->padv = v;
00417     }
00418 #endif //RAD_AUTO_AD_Const
00419 #ifdef RAD_DEBUG
00420   int gcgen;
00421   int opno;
00422   static int gcgen_cur, last_opno, zap_gcgen, zap_gcgen1, zap_opno;
00423   static FILE *debug_file;
00424 #endif
00425   static ADcontext<Double> adc;
00426   Advari_Opclass opclass;
00427   Double Val;   // result of this operation
00428   mutable Double aval;  // adjoint -- partial of final result w.r.t. this Val
00429   mutable Double dO;  // deriv of op w.r.t. t in x + t*p
00430   mutable Double aO;  // adjoint (in Hv computation) of op
00431   mutable Double adO; // adjoint (in Hv computation) of dO
00432   void *operator new(size_t len) {
00433 #ifdef RAD_DEBUG
00434     ADvari *rv = (ADvari*)ADvari::adc.Memalloc(len);
00435     rv->gcgen = gcgen_cur;
00436     rv->opno = ++last_opno;
00437     if (last_opno == zap_opno && gcgen_cur == zap_gcgen)
00438       printf("");
00439     return rv;
00440 #else
00441     return ADvari::adc.Memalloc(len);
00442 #endif
00443     }
00444   void operator delete(void*) {} /*Should never be called.*/
00445   inline ADvari(Advari_Opclass oc, Double t):
00446     opclass(oc), Val(t), aval(0.), dO(0.)
00447     { if (oc != Hv_const) ADvari::adc.ADvari_record(this); }
00448   inline ADvari(Advari_Opclass oc, Double t, Double ta):
00449     opclass(oc), Val(t), aval(ta), dO(0.)
00450     { if (oc != Hv_const) ADvari::adc.ADvari_record(this); }
00451  private:
00452   inline ADvari(): Val(0.), aval(0.), dO(0.) {} // prevent construction without value (?)
00453  public:
00454   friend class ConstADvari<Double>;
00455 #ifdef RAD_AUTO_AD_Const
00456   friend class ADcontext<Double>;
00457   friend class ADvar<Double>;
00458   friend class ADvar1<Double>;
00459   friend class ADvar1s<Double>;
00460   friend class ADvar2<Double>;
00461   friend class ADvar2q<Double>;
00462   friend class ConstADvar<Double>;
00463   ADvari(const IndepADVar *, Double);
00464   ADvari(const IndepADVar *, Double, Double);
00465   ADvari(const IndepADVar *, Double, Double, int);
00466 #endif
00467 #define F friend
00468 #define R ADvari&
00469 #define Ai const ADvari&
00470 #define T1(r,f) F r f <>(Ai);
00471 #define T2(r,f) \
00472 F r f <>(Ai,Ai); \
00473 F r f <>(Ttype,Ai); \
00474 F r f <>(Ai,Ttype); \
00475 F r f <>(double,Ai); \
00476 F r f <>(Ai,double); \
00477 F r f <>(long,Ai); \
00478 F r f <>(Ai,long); \
00479 F r f <>(int,Ai); \
00480 F r f <>(Ai,int);
00481   T1(R,operator+)
00482   T2(R,operator+)
00483   T1(R,operator-)
00484   T2(R,operator-)
00485   T2(R,operator*)
00486   T2(R,operator/)
00487   T1(R,abs)
00488   T1(R,acos)
00489   T1(R,acosh)
00490   T1(R,asin)
00491   T1(R,asinh)
00492   T1(R,atan)
00493   T1(R,atanh)
00494   T2(R,atan2)
00495   T2(R,max)
00496   T2(R,min)
00497   T1(R,cos)
00498   T1(R,cosh)
00499   T1(R,exp)
00500   T1(R,log)
00501   T1(R,log10)
00502   T2(R,pow)
00503   T1(R,sin)
00504   T1(R,sinh)
00505   T1(R,sqrt)
00506   T1(R,tan)
00507   T1(R,tanh)
00508   T1(R,fabs)
00509   T1(R,copy)
00510   T2(int,operator<)
00511   T2(int,operator<=)
00512   T2(int,operator==)
00513   T2(int,operator!=)
00514   T2(int,operator>=)
00515   T2(int,operator>)
00516 #undef T2
00517 #undef T1
00518 #undef Ai
00519 #undef R
00520 #undef F
00521 
00522   friend ADvari& ADf1<>(Double f, Double g, Double h, const ADvari &x);
00523   friend ADvari& ADf2<>(Double f, Double gx, Double gy, Double hxx,
00524       Double hxy, Double hyy, const ADvari &x, const ADvari &y);
00525   friend ADvari& ADfn<>(Double f, int n, const IndepADVar *x,
00526       const Double *g, const Double *h);
00527 
00528   inline operator Double() { return this->Val; }
00529   inline operator Double() const { return this->Val; }
00530   };
00531 
00532  template<typename Double> class
00533 ADvar1: public ADvari<Double> { // simplest unary ops
00534  public:
00535   typedef ADvari<Double> ADVari;
00536   Derp<Double> d;
00537   ADvar1(Advari_Opclass oc, Double val1, const Double *a1, const ADVari *c1):
00538     ADVari(oc,val1), d(a1,this,c1) {}
00539 #ifdef RAD_AUTO_AD_Const
00540   typedef typename ADVari::IndepADVar IndepADVar;
00541   typedef ADvar<Double> ADVar;
00542   ADvar1(const IndepADVar*, const IndepADVar&);
00543   ADvar1(const IndepADVar*, const ADVari&);
00544   ADvar1(Advari_Opclass oc, const Double val1, const Double *a1,
00545     const ADVari *c1, const ADVar *v):
00546       ADVari(oc,val1), d(a1,this,c1) {
00547     c1->padv = 0;
00548     this->ADvari_padv(v);
00549     }
00550 #endif
00551   };
00552 
00553 
00554  template<typename Double> class
00555 ConstADvari: public ADvari<Double> {
00556  private:
00557   ConstADvari *prevcad;
00558   ConstADvari() {}; // prevent construction without value (?)
00559   static ConstADvari *lastcad;
00560  public:
00561   typedef ADvari<Double> ADVari;
00562   typedef Derp<Double> DErp;
00563   static CADcontext<Double> cadc;
00564   inline void *operator new(size_t len) { return ConstADvari::cadc.Memalloc(len); }
00565   inline ConstADvari(Double t): ADVari(Hv_copy, t) { prevcad = lastcad; lastcad = this; }
00566   static void aval_reset(void);
00567   };
00568 
00569 
00570  template<typename Double> class
00571 IndepADvar {    // an independent ADvar
00572  private:
00573   IndepADvar& operator=(IndepADvar&x) {
00574     /* private to prevent assignment */
00575 #ifdef RAD_AUTO_AD_Const
00576     if (cv)
00577       cv->padv = 0;
00578     return new ADvar1<Double>(this,x);
00579 #elif defined(RAD_EQ_ALIAS)
00580     this->cv = x.cv;
00581     return *this;
00582 #else
00583     return ADvar_operatoreq(this,*x.cv);
00584 #endif //RAD_AUTO_AD_Const
00585     }
00586  protected:
00587   static void AD_Const(const IndepADvar&);
00588   mutable ADvari<Double> *cv;
00589  public:
00590   typedef Double value_type;
00591   friend class ADvar<Double>;
00592   friend class ADcontext<Double>;
00593   friend class ADvar1<Double>;
00594   friend class ADvarn<Double>;
00595   typedef ADvari<Double> ADVari;
00596   typedef ADvar<Double> ADVar;
00597   IndepADvar(Ttype);
00598   IndepADvar(double);
00599   IndepADvar(int);
00600   IndepADvar(long);
00601   IndepADvar& operator= (Double);
00602 #ifdef RAD_AUTO_AD_Const
00603   friend IndepADvar& ADvar_operatoreq<>(IndepADvar*, const ADVari&);
00604   inline IndepADvar(const IndepADvar &x) { cv = x.cv ? new ADvar1<Double>(this, x) : 0; };
00605   inline IndepADvar(ADVari *ncv) { cv = ncv; }
00606   inline IndepADvar() { cv = 0; }
00607   inline ~IndepADvar() {
00608           if (cv)
00609             cv->padv = 0;
00610           }
00611 #else
00612   inline IndepADvar() {
00613 #ifndef RAD_EQ_ALIAS
00614     cv = 0;
00615 #endif
00616     }
00617   inline ~IndepADvar() {}
00618   friend IndepADvar& ADvar_operatoreq<>(IndepADvar*, const ADVari&);
00619 #endif
00620 
00621 #ifdef RAD_Const_WARN
00622   inline operator ADVari&() const {
00623     ADVari *tcv = this->cv;
00624     if (tcv->opno < 0)
00625       RAD_Const_Warn(tcv);
00626     return *tcv;
00627     }
00628   inline operator ADVari*() const {
00629     ADVari *tcv = this->cv;
00630     if (tcv->opno < 0)
00631       RAD_Const_Warn(tcv);
00632     return tcv;
00633     }
00634 #else //RAD_Const_WARN
00635   inline operator ADVari&() const { return *this->cv; }
00636   inline operator ADVari*() const { return this->cv; }
00637 #endif //RAD_Const_WARN
00638 
00639   Double val() const { return cv->Val; }
00640   Double adj() const { return cv->aval; }
00641 
00642   friend void AD_Const1<>(Double*, const IndepADvar&);
00643 
00644   friend ADVari& ADf1<>(Double, Double, const IndepADvar&);
00645   friend ADVari& ADf2<>(Double, Double, Double, const IndepADvar&, const IndepADvar&);
00646   friend ADVari& ADf2<>(Double, Double, Double, const ADVari&, const IndepADvar&);
00647   friend ADVari& ADf2<>(Double, Double, Double, const IndepADvar&, const ADVari&);
00648 
00649   static inline void Gradcomp(int wantgrad)
00650         { ADcontext<Double>::Gradcomp(wantgrad); }
00651   static inline void Gradcomp()
00652         { ADcontext<Double>::Gradcomp(1); }
00653   static inline void Hvprod(int n, ADVar **vp, Double *v, Double *hv)
00654         { ADcontext<Double>::Hvprod(n, vp, v, hv); }
00655   static inline void aval_reset() { ConstADvari<Double>::aval_reset(); }
00656   static inline void Weighted_Gradcomp(int n, ADVar **v, Double *w)
00657         { ADcontext<Double>::Weighted_Gradcomp(n, v, w); }
00658 
00659   /* We use #define to deal with bizarre templating rules that apparently */
00660   /* require us to spell the some conversion explicitly */
00661 
00662 
00663 #define Ai const ADVari&
00664 #define AI const IndepADvar&
00665 #define D Double
00666 #define T2(r,f) \
00667  r f <>(AI,AI);\
00668  r f <>(Ai,AI);\
00669  r f <>(AI,Ai);\
00670  r f <>(Ttype,AI);\
00671  r f <>(double,AI);\
00672  r f <>(long,AI);\
00673  r f <>(int,AI);\
00674  r f <>(AI,Ttype);\
00675  r f <>(AI,double);\
00676  r f <>(AI,long);\
00677  r f <>(AI,int);
00678 #define T1(f) friend ADVari& f<> (AI);
00679 
00680 #define F friend ADVari&
00681 T2(F, operator+)
00682 T2(F, operator-)
00683 T2(F, operator*)
00684 T2(F, operator/)
00685 T2(F, atan2)
00686 T2(F, max)
00687 T2(F, min)
00688 T2(F, pow)
00689 #undef F
00690 #define F friend int
00691 T2(F, operator<)
00692 T2(F, operator<=)
00693 T2(F, operator==)
00694 T2(F, operator!=)
00695 T2(F, operator>=)
00696 T2(F, operator>)
00697 
00698 T1(operator+)
00699 T1(operator-)
00700 T1(abs)
00701 T1(acos)
00702 T1(acosh)
00703 T1(asin)
00704 T1(asinh)
00705 T1(atan)
00706 T1(atanh)
00707 T1(cos)
00708 T1(cosh)
00709 T1(exp)
00710 T1(log)
00711 T1(log10)
00712 T1(sin)
00713 T1(sinh)
00714 T1(sqrt)
00715 T1(tan)
00716 T1(tanh)
00717 T1(fabs)
00718 T1(copy)
00719 
00720 #undef F
00721 #undef T1
00722 #undef T2
00723 #undef D
00724 #undef AI
00725 #undef Ai
00726 
00727   };
00728 
00729  template<typename Double> class
00730 ADvar: public IndepADvar<Double> {  // an "active" variable
00731  public:
00733         template <typename U> struct apply { typedef ADvar<U> type; };
00734 
00735   typedef IndepADvar<Double> IndepADVar;
00736   typedef ADvari<Double> ADVari;
00737   typedef ConstADvari<Double> ConstADVari;
00738  private:
00739   void ADvar_ctr(Double d) {
00740     ADVari *x;
00741     if (ConstADVari::cadc.fpval_implies_const)
00742       x = new ConstADVari(d);
00743     else {
00744 #ifdef RAD_AUTO_AD_Const
00745       x = new ADVari((IndepADVar*)this, d);
00746       x->ADvari_padv(this);
00747 #else
00748       x = new ADVari(Hv_const, d);
00749 #endif
00750       }
00751     this->cv = x;
00752     }
00753  public:
00754   friend class ADvar1<Double>;
00755   typedef ADvar1<Double> ADVar1;
00756   ADvar() { /* cv = 0; */ }
00757   ADvar(Ttype d)  { ADvar_ctr(d); }
00758   ADvar(double i) { ADvar_ctr(Double(i)); }
00759   ADvar(int i)  { ADvar_ctr(Double(i)); }
00760   ADvar(long i) { ADvar_ctr(Double(i)); }
00761   inline ~ADvar() {}
00762 #ifdef RAD_AUTO_AD_Const
00763   ADvar(IndepADVar &x) {
00764     this->cv = x.cv ? new ADVar1(this, x) : 0;
00765     }
00766   ADvar(ADVari &x) :IndepADVar(&x) { x.ADvari_padv((IndepADVar*)this);}
00767   inline ADvar& operator=(IndepADVar &x) {
00768     if (this->cv)
00769       this->cv->padv = 0;
00770     this->cv = new ADVar1(this,x);
00771     return *this;
00772     }
00773   inline ADvar& operator=(ADVari &x) {
00774     if (this->cv)
00775       this->cv->padv = 0;
00776     this->cv = new ADVar1(this, x);
00777     return *this;
00778     }
00779 #else 
00780   friend ADvar& ADvar_operatoreq<>(ADvar*, const ADVari&);
00781 #ifdef RAD_EQ_ALIAS
00782   /* allow aliasing v and w after "v = w;" */
00783   inline ADvar(const IndepADVar &x) { this->cv = (ADVari*)x.cv; }
00784   inline ADvar(const ADVari &x) { this->cv = (ADVari*)&x; }
00785   inline ADvar& operator=(IndepADVar &x) { this->cv = (ADVari*)x.cv; return *this; }
00786   inline ADvar& operator=(const ADVari &x) { this->cv = (ADVari*)&x; return *this; }
00787 #else 
00788   friend ADvar& ADvar_operatoreq<>(ADvar*, const ADVari&);
00789   ADvar(const IndepADVar &x) { this->cv = x.cv ? new ADVar1(x.cv->Val, &this->cv->adc.One, x.cv) : 0; }
00790   ADvar(const ADvar&x) { this->cv = x.cv ?
00791     new ADVar1(Hv_copy, x.cv->Val, &this->cv->adc.One, (ADVari*)x.cv) : 0; }
00792   ADvar(const ADVari &x) { this->cv = new ADVar1(Hv_copy, x.Val, &this->cv->adc.One, &x); }
00793   inline ADvar& operator=(IndepADVar &x) { return ADvar_operatoreq(this,*x.cv); };
00794   inline ADvar& operator=(const ADVari &x) { return ADvar_operatoreq(this,x); };
00795 #endif /* RAD_EQ_ALIAS */
00796 #endif /* RAD_AUTO_AD_Const */
00797   ADvar& operator=(Double);
00798   ADvar& operator+=(const ADVari&);
00799   ADvar& operator+=(Double);
00800   ADvar& operator-=(const ADVari&);
00801   ADvar& operator-=(Double);
00802   ADvar& operator*=(const ADVari&);
00803   ADvar& operator*=(Double);
00804   ADvar& operator/=(const ADVari&);
00805   ADvar& operator/=(Double);
00806   inline static bool get_fpval_implies_const(void)
00807     { return ConstADVari::cadc.fpval_implies_const; }
00808   inline static void set_fpval_implies_const(bool newval)
00809     { ConstADVari::cadc.fpval_implies_const = newval; }
00810   inline static bool setget_fpval_implies_const(bool newval) {
00811     bool oldval = ConstADVari::cadc.fpval_implies_const;
00812     ConstADVari::cadc.fpval_implies_const = newval;
00813     return oldval;
00814     }
00815   static inline void Gradcomp(int wantgrad)
00816         { ADcontext<Double>::Gradcomp(wantgrad); }
00817   static inline void Gradcomp()
00818         { ADcontext<Double>::Gradcomp(1); }
00819   static inline void aval_reset() { ConstADVari::aval_reset(); }
00820   static inline void Weighted_Gradcomp(int n, ADvar **v, Double *w)
00821         { ADcontext<Double>::Weighted_Gradcomp(n, v, w); }
00822   };
00823 
00824 template<typename Double>
00825  inline void AD_Const1(Double *notused, const IndepADvar<Double>&v)
00826 { IndepADvar<Double>::AD_Const(v); }
00827 
00828 template<typename Double>
00829  inline void AD_Const(const IndepADvar<Double>&v) { AD_Const1((Double*)0, v); }
00830 
00831  template<typename Double> class
00832 ConstADvar: public ADvar<Double> {
00833  public:
00834   typedef ADvar<Double> ADVar;
00835   typedef typename ADVar::ADVari ADVari;
00836   typedef typename ADVar::ConstADVari ConstADVari;
00837   typedef Derp<Double> DErp;
00838   typedef typename ADVar::IndepADVar IndepADVar;
00839  private: // disable op=
00840   ConstADvar& operator+=(ADVari&);
00841   ConstADvar& operator+=(Double);
00842   ConstADvar& operator-=(ADVari&);
00843   ConstADvar& operator-=(Double);
00844   ConstADvar& operator*=(ADVari&);
00845   ConstADvar& operator*=(Double);
00846   ConstADvar& operator/=(ADVari&);
00847   ConstADvar& operator/=(Double);
00848   void ConstADvar_ctr(Double);
00849  public:
00850   ConstADvar(Ttype d) { ConstADvar_ctr(d); }
00851   ConstADvar(double i)  { ConstADvar_ctr(Double(i)); }
00852   ConstADvar(int i) { ConstADvar_ctr(Double(i)); }
00853   ConstADvar(long i)  { ConstADvar_ctr(Double(i)); }
00854   ConstADvar(const IndepADVar&);
00855   ConstADvar(const ConstADvar&);
00856   ConstADvar(const ADVari&);
00857   inline ~ConstADvar() {}
00858 #ifdef RAD_NO_CONST_UPDATE
00859  private:
00860 #endif
00861   ConstADvar();
00862   inline ConstADvar& operator=(Double d) { this->cv->Val = d; return *this; }
00863   inline ConstADvar& operator=(ADVari& d) { this->cv->Val = d.Val; return *this; }
00864  };
00865 
00866  template<typename Double> class
00867 ADvar1s: public ADvar1<Double> { // unary ops with partials
00868  public:
00869   typedef ADvar1<Double> ADVar1;
00870   typedef typename ADVar1::ADVari ADVari;
00871   Double pL;  // deriv of op w.r.t. left operand L
00872   ADvar1s(Double val1, Double a1, const ADVari *c1):
00873     ADVar1(Hv_timesL,val1,&pL,c1), pL(a1) {}
00874 #ifdef RAD_AUTO_AD_Const
00875   Double a;
00876   typedef typename ADVar1::ADVar ADVar;
00877   ADvar1s(Double val1, Double a1, const ADVari *c1, const ADVar *v):
00878     ADVar1(Hv_timesL,val1,&a,c1,v), a(a1) {}
00879 #endif
00880   };
00881 
00882  template<typename Double> class
00883 ADvar1g: public ADvar1<Double> { // unary ops with partials
00884  public:
00885   typedef ADvar1<Double> ADVar1;
00886   typedef typename ADVar1::ADVari ADVari;
00887   Double pL;  // deriv of op w.r.t. left operand L
00888   Double pL2; // partial of op w.r.t. L,L
00889   ADvar1g(Double val1, Double d1, Double d2, const ADVari *c1):
00890     ADVar1(Hv_unary,val1,&pL,c1), pL(d1), pL2(d2) {}
00891   };
00892 
00893  template<typename Double> class
00894 ADvar2: public ADvari<Double> { // basic binary ops
00895  public:
00896   typedef ADvari<Double> ADVari;
00897   typedef Derp<Double> DErp;
00898   DErp dL, dR;
00899   ADvar2(Advari_Opclass oc, Double val1, const ADVari *Lcv, const Double *Lc,
00900     const ADVari *Rcv, const Double *Rc):
00901       ADVari(oc,val1) {
00902     dR.next = DErp::LastDerp;
00903     dL.next = &dR;
00904     DErp::LastDerp = &dL;
00905     dL.a = Lc;
00906     dL.c = (ADVari*)Lcv;
00907     dR.a = Rc;
00908     dR.c = (ADVari*)Rcv;
00909     dL.b = dR.b = this;
00910     }
00911 #ifdef RAD_AUTO_AD_Const
00912   typedef ADvar<Double> ADVar;
00913   ADvar2(Advari_Opclass oc, Double val1, const ADVari *Lcv, const Double *Lc,
00914     const ADVari *Rcv, const Double *Rc, ADVar *v):
00915       ADVari(oc,val1) {
00916     dR.next = DErp::LastDerp;
00917     dL.next = &dR;
00918     DErp::LastDerp = &dL;
00919     dL.a = Lc;
00920     dL.c = Lcv;
00921     dR.a = Rc;
00922     dR.c = Rcv;
00923     dL.b = dR.b = this;
00924     Lcv->padv = 0;
00925     this->ADvari_padv(v);
00926     }
00927 #endif
00928   };
00929 
00930  template<typename Double> class
00931 ADvar2q: public ADvar2<Double> { // binary ops with partials
00932  public:
00933   typedef ADvar2<Double> ADVar2;
00934   typedef typename ADVar2::ADVari ADVari;
00935   typedef typename ADVar2::DErp DErp;
00936   Double pL;  // deriv of op w.r.t. left operand L
00937   Double pR;  // deriv of op w.r.t. right operand R
00938   Double pLR; // second partial w.r.t. L,R
00939   Double pR2; // second partial w.r.t. R,R
00940   ADvar2q(Double val1, Double Lp, Double Rp, Double LR, Double R2,
00941     const ADVari *Lcv, const ADVari *Rcv):
00942       ADVar2(Hv_quotLR,val1,Lcv,&pL,Rcv,&pR),
00943       pL(Lp), pR(Rp), pLR(LR), pR2(R2) {}
00944 #ifdef RAD_AUTO_AD_Const
00945   typedef typename ADVar2::ADVar ADVar;
00946   ADvar2q(Double val1, Double Lp, Double Rp, Double LR, Double R2,
00947     const ADVari *Lcv, const ADVari *Rcv, const ADVar *v):
00948       ADVar2(Hv_quotLR,val1,Lcv,&pL,Rcv,&pR,v),
00949       pL(Lp), pR(Rp), pLR(LR), pR2(R2) {}
00950 #endif
00951   };
00952 
00953  template<typename Double> class
00954 ADvar2g: public ADvar2<Double> { // general binary ops with partials
00955  public:
00956   typedef ADvar2<Double> ADVar2;
00957   typedef typename ADVar2::ADVari ADVari;
00958   Double pL;  // deriv of op w.r.t. left operand L
00959   Double pR;  // deriv of op w.r.t. right operand R
00960   Double pL2; // second partial w.r.t. L,L
00961   Double pLR; // second partial w.r.t. L,R
00962   Double pR2; // second partial w.r.t. R,R
00963   ADvar2g(Double val1, Double Lp, Double Rp, Double L2, Double LR, Double R2,
00964     const ADVari *Lcv, const ADVari *Rcv):
00965       ADVar2(Hv_binary,val1,Lcv,&pL,Rcv,&pR),
00966       pL(Lp), pR(Rp), pL2(L2), pLR(LR), pR2(R2) { }
00967   };
00968 
00969  template<typename Double> class
00970 ADvarn: public ADvari<Double> { // n-ary ops with partials g and
00971         // 2nd partials h (lower triangle, rowwise)
00972  public:
00973   typedef ADvari<Double> ADVari;
00974   typedef ADvar<Double> ADVar;
00975   typedef typename ADVari::IndepADVar IndepADVar;
00976   typedef Derp<Double> DErp;
00977   int n;
00978   Double *G, *H;
00979   DErp *D;
00980   ADvarn(Double val1, int n1, const IndepADVar *x, const Double *g, const Double *h):
00981       ADVari(Hv_nary,val1), n(n1) {
00982     DErp *d1, *dlast;
00983     Double *a1;
00984     int i, nh;
00985 
00986     a1 = G = (Double*)ADVari::adc.Memalloc(n1*sizeof(*g));
00987     d1 = D = (DErp*)ADVari::adc.Memalloc(n1*sizeof(DErp));
00988     dlast = DErp::LastDerp;
00989     for(i = 0; i < n1; i++, d1++) {
00990       d1->next = dlast;
00991       dlast = d1;
00992       a1[i] = g[i];
00993       d1->a = &a1[i];
00994       d1->b = this;
00995       d1->c = x[i].cv;
00996       }
00997     DErp::LastDerp = dlast;
00998     nh = (n1*(n1+1)) >> 1;
00999     a1 = H = (double*)ADVari::adc.Memalloc(nh * sizeof(*h));
01000     for(i = 0; i < nh; i++)
01001       a1[i] = h[i];
01002     }
01003   };
01004 
01005 template<typename Double>
01006  inline ADvari<Double>& operator+(const ADvari<Double> &T) { return *(ADvari<Double>*)&T; }
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 int operator>(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val > R.Val; }
01045 template<typename Double>
01046  inline int operator>(const ADvari<Double> &L, Double R) { return L.Val > R; }
01047 template<typename Double>
01048  inline int operator>(Double L, const ADvari<Double> &R) { return L > R.Val; }
01049 
01050 template<typename Double>
01051  inline void *ADcontext<Double>::Memalloc(size_t len) {
01052     if (Mleft >= len)
01053       return Mbase + (Mleft -= len);
01054     return new_ADmemblock(len);
01055     }
01056 
01057 template<typename Double>
01058  inline Derp<Double>::Derp(const ADVari *c1): c((ADVari*)c1) {
01059     next = LastDerp;
01060     LastDerp = this;
01061     }
01062 
01063 template<typename Double>
01064  inline Derp<Double>::Derp(const Double *a1, const ADVari *c1): a(a1), c((ADVari*)c1) {
01065     next = LastDerp;
01066     LastDerp = this;
01067     }
01068 
01069 template<typename Double>
01070  inline Derp<Double>::Derp(const Double *a1, const ADVari *b1, const ADVari *c1):
01071   a(a1), b(b1), c((ADVari*)c1) {
01072     next = LastDerp;
01073     LastDerp = this;
01074     }
01075 
01076 /**** radops ****/
01077 
01078 template<typename Double> Derp<Double> *Derp<Double>::LastDerp = 0;
01079 template<typename Double> ADcontext<Double> ADvari<Double>::adc;
01080 template<typename Double> const Double ADcontext<Double>::One = 1.;
01081 template<typename Double> const Double ADcontext<Double>::negOne = -1.;
01082 template<typename Double> CADcontext<Double> ConstADvari<Double>::cadc;
01083 template<typename Double> ConstADvari<Double> *ConstADvari<Double>::lastcad;
01084 
01085 #ifdef RAD_AUTO_AD_Const
01086 template<typename Double> ADvari<Double>*   ADvari<Double>::First_ADvari;
01087 template<typename Double> ADvari<Double>**  ADvari<Double>::Last_ADvari = &ADvari<Double>::First_ADvari;
01088 #endif
01089 
01090 #ifdef RAD_DEBUG
01091 #ifndef RAD_DEBUG_gcgen1
01092 #define RAD_DEBUG_gcgen1 -1
01093 #endif
01094 template<typename Double> int ADvari<Double>::gcgen_cur;
01095 template<typename Double> int ADvari<Double>::last_opno;
01096 template<typename Double> int ADvari<Double>::zap_gcgen;
01097 template<typename Double> int ADvari<Double>::zap_gcgen1 = RAD_DEBUG_gcgen1;
01098 template<typename Double> int ADvari<Double>::zap_opno;
01099 template<typename Double> FILE *ADvari<Double>::debug_file;
01100 #endif
01101 
01102 
01103  template<typename Double>
01104 ADcontext<Double>::ADcontext()
01105 {
01106   ADVari_block *fb;
01107 
01108   First = (ADMemblock*)First0;
01109   First->next = 0;
01110   Busy = First;
01111   Free = 0;
01112   Mbase = (char*)First->memblk;
01113   Mleft = sizeof(First->memblk);
01114   AiFirst = Aibusy = fb = (ADVari_block*)First1;
01115   Aifree = 0;
01116   Ainext = fb->pADvari;
01117   fb->next = fb->prev = 0;
01118   fb->limit = Ailimit = fb->pADvari + ADVari_block::Gulp;
01119   rad_need_reinit = 0;
01120 #ifdef RAD_DEBUG_BLOCKKEEP
01121   rad_busy_blocks = 0;
01122   rad_mleft_save = 0;
01123   rad_Oldcurmb = 0;
01124 #endif
01125   }
01126 
01127 template<typename Double> void*
01128 ADcontext<Double>::new_ADmemblock(size_t len)
01129 {
01130   ADMemblock *mb, *mb0, *mb1, *mbf, *x;
01131   ADVari_block *b;
01132 #ifdef RAD_AUTO_AD_Const
01133   ADVari *a, *anext;
01134   IndepADvar<Double> *v;
01135 #ifdef RAD_Const_WARN
01136   ADVari *cv;
01137   int i, j;
01138 #endif
01139 #endif /*RAD_AUTO_AD_Const*/
01140 
01141   if ((rad_need_reinit & 1) && this == &ADVari::adc) {
01142     rad_need_reinit &= ~1;
01143     DErp::LastDerp = 0;
01144     Aibusy = b = AiFirst;
01145     Aifree = b->next;
01146     b->next = b->prev = 0;
01147     Ailimit = b->limit = (Ainext = b->pADvari) + ADVari_block::Gulp;
01148 #ifdef RAD_DEBUG_BLOCKKEEP
01149     Mleft = rad_mleft_save;
01150     if (Mleft < sizeof(First->memblk))
01151       _uninit_f2c(Mbase + Mleft,
01152         UninitType<Double>::utype,
01153         (sizeof(First->memblk) - Mleft)
01154         /sizeof(typename Sacado::ValueType<Double>::type));
01155     if ((mb = Busy->next)) {
01156       if (!(mb0 = rad_Oldcurmb))
01157         mb0 = (ADMemblock*)First0;
01158       for(;; mb = mb->next) {
01159         _uninit_f2c(mb->memblk,
01160           UninitType<Double>::utype,
01161           sizeof(First->memblk)
01162           /sizeof(typename Sacado::ValueType<Double>::type));
01163         if (mb == mb0)
01164           break;
01165         }
01166       }
01167     rad_Oldcurmb = Busy;
01168     if (rad_busy_blocks >= RAD_DEBUG_BLOCKKEEP) {
01169       rad_busy_blocks = 0;
01170       rad_Oldcurmb = 0;
01171       mb0 = (ADMemblock*)First0;
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       }
01183 
01184 #else /* !RAD_DEBUG_BLOCKKEEP */
01185 
01186     mb0 = First;
01187     mbf =  Free;
01188     for(mb = Busy; mb != mb0; mb = mb1) {
01189       mb1 = mb->next;
01190       mb->next = mbf;
01191       mbf = mb;
01192       }
01193     Free = mbf;
01194     Busy = mb;
01195     Mbase = (char*)First->memblk;
01196     Mleft = sizeof(First->memblk);
01197 #ifdef RAD_AUTO_AD_Const
01198     if (ADVari::adc.rad_need_reinit & 2) {
01199       ADVari::adc.rad_need_reinit &= ~2;
01200       *ADVari::Last_ADvari = 0;
01201       ADVari::Last_ADvari = &ADVari::First_ADvari;
01202       if ((anext = ADVari::First_ADvari)) {
01203         while((a = anext)) {
01204           anext = a->Next;
01205           if ((v = (IndepADvar<Double> *)a->padv)) {
01206 #ifdef RAD_Const_WARN
01207             if ((i = a->opno) > 0)
01208               i = -i;
01209             j = a->gcgen;
01210             v->cv = cv = new ADVari(v, a->Val, a->aval);
01211             cv->opno = i;
01212             cv->gcgen = j;
01213 #else
01214             v->cv = new ADVari(v, a->Val, a->aval);
01215 #endif
01216             }
01217           }
01218         }
01219       }
01220 #endif /*RAD_AUTO_AD_Const*/
01221 #endif /*RAD_DEBUG_BLOCKKEEP*/
01222     if (Mleft >= len)
01223       return Mbase + (Mleft -= len);
01224     }
01225 
01226   if ((x = Free))
01227     Free = x->next;
01228   else
01229     x = new ADMemblock;
01230 #ifdef RAD_DEBUG_BLOCKKEEP
01231   rad_busy_blocks++;
01232 #endif
01233   x->next = Busy;
01234   Busy = x;
01235   return (Mbase = (char*)x->memblk) +
01236     (Mleft = sizeof(First->memblk) - len);
01237   }
01238 
01239 template<typename Double> void
01240 ADcontext<Double>::new_ADvari_block()
01241 {
01242   ADVari_block *ob, *nb;
01243   ob = Aibusy;
01244   ob->limit = Ailimit;  // should be unnecessary, but harmless
01245   if ((nb = Aifree))
01246     Aifree = nb->next;
01247   else
01248     nb = new ADVari_block;
01249   Aibusy = Aibusy->next = nb;
01250   nb->limit = Ailimit = (Ainext = nb->pADvari) + ADVari_block::Gulp;
01251   ob->next = nb;
01252   nb->prev = ob;
01253   nb->next = 0;
01254   }
01255 
01256 template<typename Double> void
01257 ADcontext<Double>::Gradcomp(int wantgrad)
01258 {
01259   DErp *d;
01260 
01261   if (ADVari::adc.rad_need_reinit) {
01262     for(d = DErp::LastDerp; d; d = d->next)
01263       d->c->aval = 0;
01264     }
01265   if (!(ADVari::adc.rad_need_reinit & 1)) {
01266     ADVari::adc.rad_need_reinit = 1;
01267     ADVari::adc.rad_mleft_save = ADVari::adc.Mleft;
01268     ADVari::adc.Mleft = 0;
01269     }
01270 #ifdef RAD_DEBUG
01271   if (ADVari::gcgen_cur == ADVari::zap_gcgen1) {
01272     const char *fname;
01273     if (!(fname = getenv("RAD_DEBUG_FILE")))
01274       fname = "rad_debug.out";
01275     else if (!*fname)
01276       fname = 0;
01277     if (fname)
01278       ADVari::debug_file = fopen(fname, "w");
01279     ADVari::zap_gcgen1 = -1;
01280     }
01281 #endif
01282   if ((d = DErp::LastDerp)) {
01283 #ifdef RAD_AUTO_AD_Const
01284      ADVari::adc.rad_need_reinit |= 2;
01285 #endif /*RAD_AUTO_AD_Const*/
01286     if (wantgrad) {
01287       d->b->aval = 1;
01288 #ifdef RAD_DEBUG
01289       if (ADVari::debug_file)
01290         do {
01291           fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
01292             d->c->opno, d->b->opno, d->c->aval,
01293             *d->a, d->b->aval);
01294           d->c->aval += *d->a * d->b->aval;
01295           fprintf(ADVari::debug_file, " = %g\n", d->c->aval);
01296           fflush(ADVari::debug_file);
01297           } while((d = d->next));
01298       else
01299 #endif
01300       do d->c->aval += *d->a * d->b->aval;
01301       while((d = d->next));
01302       DErp::LastDerp = 0;
01303       }
01304     }
01305 #ifdef RAD_DEBUG
01306   if (ADVari::debug_file) {
01307     fclose(ADVari::debug_file);
01308     ADVari::debug_file = 0;
01309     }
01310 #endif //RAD_DEBUG
01311 #ifdef RAD_DEBUG
01312   ADVari::gcgen_cur++;
01313   ADVari::last_opno = 0;
01314 #endif
01315   }
01316 
01317 template<typename Double> void
01318 ADcontext<Double>::Weighted_Gradcomp(int n, ADVar **V, Double *w)
01319 {
01320   DErp *d;
01321   int i;
01322 #ifdef RAD_Const_WARN
01323   ADVari *cv;
01324   int j;
01325 #endif
01326 #ifdef RAD_AUTO_AD_Const
01327   ADVari *a, *anext;
01328   IndepADvar<Double> *v;
01329 #endif /*RAD_AUTO_AD_Const*/
01330 
01331   if (ADVari::adc.rad_need_reinit) {
01332     for(d = DErp::LastDerp; d; d = d->next)
01333       d->c->aval = 0;
01334     }
01335   if (!(ADVari::adc.rad_need_reinit & 1)) {
01336     ADVari::adc.rad_need_reinit = 1;
01337     ADVari::adc.rad_mleft_save = ADVari::adc.Mleft;
01338     ADVari::adc.Mleft = 0;
01339     }
01340 #ifdef RAD_DEBUG
01341   if (ADVari::gcgen_cur == ADVari::zap_gcgen1) {
01342     const char *fname;
01343     if (!(fname = getenv("RAD_DEBUG_FILE")))
01344       fname = "rad_debug.out";
01345     else if (!*fname)
01346       fname = 0;
01347     if (fname)
01348       ADVari::debug_file = fopen(fname, "w");
01349     ADVari::zap_gcgen1 = -1;
01350     }
01351 #endif
01352   if ((d = DErp::LastDerp)) {
01353     for(i = 0; i < n; i++)
01354       V[i]->cv->aval = w[i];
01355 #ifdef RAD_DEBUG
01356     if (ADVari::debug_file)
01357       do {
01358         fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
01359           d->c->opno, d->b->opno, d->c->aval, *d->a, d->b->aval);
01360         d->c->aval += *d->a * d->b->aval;
01361         fprintf(ADVari::debug_file, " = %g\n", d->c->aval);
01362         fflush(ADVari::debug_file);
01363         } while((d = d->next));
01364     else
01365 #endif
01366     do d->c->aval += *d->a * d->b->aval;
01367     while((d = d->next));
01368     }
01369 #ifdef RAD_DEBUG
01370   if (ADVari::debug_file) {
01371     fclose(ADVari::debug_file);
01372     ADVari::debug_file = 0;
01373     }
01374 #endif //RAD_DEBUG
01375 #ifdef RAD_AUTO_AD_Const
01376   *ADVari::Last_ADvari = 0;
01377   ADVari::Last_ADvari = &ADVari::First_ADvari;
01378   if ((anext = ADVari::First_ADvari) && !(ADVari::adc.rad_need_reinit & 2)) {
01379     ADVari::adc.rad_need_reinit = 3;
01380     while((a = anext)) {
01381       anext = a->Next;
01382       if ((v = (IndepADvar<Double> *)a->padv)) {
01383 #ifdef RAD_Const_WARN
01384         if ((i = a->opno) > 0)
01385           i = -i;
01386         j = a->gcgen;
01387         v->cv = cv = new ADVari(v, a->Val, a->aval);
01388         cv->opno = i;
01389         cv->gcgen = j;
01390 #else
01391         v->cv = new ADVari(v, a->Val, a->aval);
01392 #endif
01393         }
01394       }
01395     DErp::LastDerp = 0;
01396     }
01397 #endif /*RAD_AUTO_AD_Const*/
01398 #ifdef RAD_DEBUG
01399   ADVari::gcgen_cur++;
01400   ADVari::last_opno = 0;
01401 #endif
01402     }
01403 
01404  template<typename Double>
01405 IndepADvar<Double>::IndepADvar(Ttype d)
01406 {
01407 
01408   ADVari *x = new ADVari(Hv_const, d);
01409   cv = x;
01410   }
01411 
01412  template<typename Double>
01413 IndepADvar<Double>::IndepADvar(double i)
01414 {
01415 
01416   ADVari *x = new ADVari(Hv_const, Double(i));
01417   cv = x;
01418   }
01419 
01420  template<typename Double>
01421 IndepADvar<Double>::IndepADvar(int i)
01422 {
01423 
01424   ADVari *x = new ADVari(Hv_const, Double(i));
01425   cv = x;
01426   }
01427 
01428  template<typename Double>
01429 IndepADvar<Double>::IndepADvar(long i)
01430 {
01431 
01432   ADVari *x = new ADVari(Hv_const, Double(i));
01433   cv = x;
01434   }
01435 
01436  template<typename Double>
01437 ConstADvar<Double>::ConstADvar()
01438 {
01439   ConstADVari *x = new ConstADVari(0.);
01440   this->cv = x;
01441   }
01442 
01443  template<typename Double> void
01444 ConstADvar<Double>::ConstADvar_ctr(Double d)
01445 {
01446   ConstADVari *x = new ConstADVari(d);
01447   this->cv = x;
01448   }
01449 
01450  template<typename Double>
01451 ConstADvar<Double>::ConstADvar(const IndepADVar &x)
01452 {
01453   ConstADVari *y = new ConstADVari(x.cv->Val);
01454   DErp *d = new DErp(&x.adc.One, y, x.cv);
01455   this->cv = y;
01456   }
01457 
01458  template<typename Double>
01459 ConstADvar<Double>::ConstADvar(const ConstADvar &x)
01460 {
01461   ConstADVari *y = new ConstADVari(x.cv->Val);
01462   DErp *d = new DErp(&x.cv->adc.One, y, (ADVari*)x.cv);
01463   this->cv = y;
01464   }
01465 
01466  template<typename Double>
01467 ConstADvar<Double>::ConstADvar(const ADVari &x)
01468 {
01469   ConstADVari *y = new ConstADVari(x.Val);
01470   DErp *d = new DErp(&x.adc.One, y, &x);
01471   this->cv = y;
01472   }
01473 
01474  template<typename Double>
01475  void
01476 IndepADvar<Double>::AD_Const(const IndepADvar &v)
01477 {
01478   typedef ConstADvari<Double> ConstADVari;
01479 
01480   ConstADVari *ncv = new ConstADVari(v.val());
01481 #ifdef RAD_AUTO_AD_Const
01482   if (v.cv)
01483     v.cv->padv = 0;
01484 #endif
01485   v.cv = ncv;
01486   }
01487 
01488  template<typename Double>
01489  void
01490 ConstADvari<Double>::aval_reset()
01491 {
01492   ConstADvari *x = ConstADvari::lastcad;
01493   while(x) {
01494     x->aval = 0;
01495     x = x->prevcad;
01496     }
01497   }
01498 
01499 #ifdef RAD_AUTO_AD_Const
01500 
01501  template<typename Double>
01502 ADvari<Double>::ADvari(const IndepADVar *x, Double d): Val(d), aval(0.)
01503 {
01504   opclass = Hv_const;
01505   this->ADvari_padv(x);
01506   }
01507 
01508  template<typename Double>
01509 ADvari<Double>::ADvari(const IndepADVar *x, Double d, Double g): Val(d), aval(g)
01510 {
01511   opclass = Hv_const;
01512   this->ADvari_padv(x);
01513   }
01514 
01515  template<typename Double>
01516 ADvar1<Double>::ADvar1(const IndepADVar *x, const IndepADVar &y):
01517   ADVari(Hv_copy, y.cv->Val), d((const Double*)&ADcontext<Double>::One, (ADVari*)this, y.cv)
01518 {
01519   this->ADvari_padv(x);
01520   }
01521 
01522  template<typename Double>
01523 ADvar1<Double>::ADvar1(const IndepADVar *x, const ADVari &y):
01524   ADVari(Hv_copy, y.Val), d((const Double*)&ADcontext<Double>::One, this, &y)
01525 {
01526   this->ADvari_padv(x);
01527   }
01528 
01529 #endif /* RAD_AUTO_AD_Const */
01530 
01531  template<typename Double>
01532  IndepADvar<Double>&
01533 ADvar_operatoreq(IndepADvar<Double> *This, const ADvari<Double> &x)
01534 { This->cv = new ADvar1<Double>(Hv_copy, x.Val, &x.adc.One, &x);
01535   return *(IndepADvar<Double>*) This; }
01536 
01537  template<typename Double>
01538  ADvar<Double>&
01539 ADvar_operatoreq(ADvar<Double> *This, const ADvari<Double> &x)
01540 { This->cv = new ADvar1<Double>(Hv_copy, x.Val, &x.adc.One, &x); return *(ADvar<Double>*) This; }
01541 
01542  template<typename Double>
01543  IndepADvar<Double>&
01544 IndepADvar<Double>::operator=(Double d)
01545 {
01546 #ifdef RAD_AUTO_AD_Const
01547   ADVari *ncv = new ADVari(this, d);
01548   if (this->cv)
01549     this->cv->padv = 0;
01550   this->cv = ncv;
01551 #else
01552   this->cv = new ADVari(Hv_const, d);
01553 #endif
01554   return *this;
01555   }
01556 
01557  template<typename Double>
01558  ADvar<Double>&
01559 ADvar<Double>::operator=(Double d)
01560 {
01561 #ifdef RAD_AUTO_AD_Const
01562   ADVari *nv = new ADVari(this, d);
01563   if (this->cv)
01564     this->cv->padv = 0;
01565   this->cv = nv;
01566 #else
01567   this->cv = ConstADVari::cadc.fpval_implies_const
01568     ? new ConstADVari(d)
01569     : new ADVari(Hv_const, d);
01570 #endif
01571   return *this;
01572   }
01573 
01574  template<typename Double>
01575  ADvari<Double>&
01576 operator-(const ADvari<Double> &T) {
01577   return *(new ADvar1<Double>(Hv_negate, -T.Val, &T.adc.negOne, &T));
01578   }
01579 
01580  template<typename Double>
01581  ADvari<Double>&
01582 operator+(const ADvari<Double> &L, const ADvari<Double> &R) {
01583   return *(new ADvar2<Double>(Hv_plusLR, L.Val + R.Val, &L, &L.adc.One, &R, &L.adc.One));
01584   }
01585 
01586 #ifdef RAD_AUTO_AD_Const
01587 #define RAD_ACA ,this
01588 #else
01589 #define RAD_ACA /*nothing*/
01590 #endif
01591 
01592  template<typename Double>
01593  ADvar<Double>&
01594 ADvar<Double>::operator+=(const ADVari &R) {
01595   ADVari *Lcv = this->cv;
01596   this->cv = new ADvar2<Double>(Hv_plusLR, Lcv->Val + R.Val, Lcv,
01597           &R.adc.One, &R, &R.adc.One RAD_ACA);
01598   return *this;
01599   }
01600 
01601  template<typename Double>
01602  ADvari<Double>&
01603 operator+(const ADvari<Double> &L, Double R) {
01604   return *(new ADvar1<Double>(Hv_copy, L.Val + R, &L.adc.One, &L));
01605   }
01606 
01607  template<typename Double>
01608  ADvar<Double>&
01609 ADvar<Double>::operator+=(Double R) {
01610   ADVari *tcv = this->cv;
01611   this->cv = new ADVar1(Hv_copy, tcv->Val + R, &tcv->adc.One, tcv RAD_ACA);
01612   return *this;
01613   }
01614 
01615  template<typename Double>
01616  ADvari<Double>&
01617 operator+(Double L, const ADvari<Double> &R) {
01618   return *(new ADvar1<Double>(Hv_copy, L + R.Val, &R.adc.One, &R));
01619   }
01620 
01621  template<typename Double>
01622  ADvari<Double>&
01623 operator-(const ADvari<Double> &L, const ADvari<Double> &R) {
01624   return *(new ADvar2<Double>(Hv_minusLR, L.Val - R.Val, &L, &L.adc.One, &R, &L.adc.negOne));
01625   }
01626 
01627  template<typename Double>
01628  ADvar<Double>&
01629 ADvar<Double>::operator-=(const ADVari &R) {
01630   ADVari *Lcv = this->cv;
01631   this->cv = new ADvar2<Double>(Hv_minusLR,Lcv->Val - R.Val, Lcv,
01632           &R.adc.One, &R, &R.adc.negOne RAD_ACA);
01633   return *this;
01634   }
01635 
01636  template<typename Double>
01637  ADvari<Double>&
01638 operator-(const ADvari<Double> &L, Double R) {
01639   return *(new ADvar1<Double>(Hv_copy, L.Val - R, &L.adc.One, &L));
01640   }
01641 
01642  template<typename Double>
01643  ADvar<Double>&
01644 ADvar<Double>::operator-=(Double R) {
01645   ADVari *tcv = this->cv;
01646   this->cv = new ADVar1(Hv_copy, tcv->Val - R, &tcv->adc.One, tcv RAD_ACA);
01647   return *this;
01648   }
01649 
01650  template<typename Double>
01651  ADvari<Double>&
01652 operator-(Double L, const ADvari<Double> &R) {
01653   return *(new ADvar1<Double>(Hv_negate, L - R.Val, &R.adc.negOne, &R));
01654   }
01655 
01656  template<typename Double>
01657  ADvari<Double>&
01658 operator*(const ADvari<Double> &L, const ADvari<Double> &R) {
01659   return *(new ADvar2<Double>(Hv_timesLR, L.Val * R.Val, &L, &R.Val, &R, &L.Val));
01660   }
01661 
01662  template<typename Double>
01663  ADvar<Double>&
01664 ADvar<Double>::operator*=(const ADVari &R) {
01665   ADVari *Lcv = this->cv;
01666   this->cv = new ADvar2<Double>(Hv_timesLR, Lcv->Val * R.Val, Lcv,
01667           &R.Val, &R, &Lcv->Val RAD_ACA);
01668   return *this;
01669   }
01670 
01671  template<typename Double>
01672  ADvari<Double>&
01673 operator*(const ADvari<Double> &L, Double R) {
01674   return *(new ADvar1s<Double>(L.Val * R, R, &L));
01675   }
01676 
01677  template<typename Double>
01678  ADvar<Double>&
01679 ADvar<Double>::operator*=(Double R) {
01680   ADVari *Lcv = this->cv;
01681   this->cv = new ADvar1s<Double>(Lcv->Val * R, R, Lcv RAD_ACA);
01682   return *this;
01683   }
01684 
01685  template<typename Double>
01686  ADvari<Double>&
01687 operator*(Double L, const ADvari<Double> &R) {
01688   return *(new ADvar1s<Double>(L * R.Val, L, &R));
01689   }
01690 
01691  template<typename Double>
01692  ADvari<Double>&
01693 operator/(const ADvari<Double> &L, const ADvari<Double> &R) {
01694   Double Lv = L.Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv, qpL = q*pL;
01695   return *(new ADvar2q<Double>(q, pL, -qpL, -pL*pL, 2.*pL*qpL, &L, &R));
01696   }
01697 
01698  template<typename Double>
01699  ADvar<Double>&
01700 ADvar<Double>::operator/=(const ADVari &R) {
01701   ADVari *Lcv = this->cv;
01702   Double Lv = Lcv->Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv, qpL = q*pL;
01703   this->cv = new ADvar2q<Double>(q, pL, -qpL, -pL*pL, 2.*pL*qpL, Lcv, &R RAD_ACA);
01704   return *this;
01705   }
01706 
01707  template<typename Double>
01708  ADvari<Double>&
01709 operator/(const ADvari<Double> &L, Double R) {
01710   return *(new ADvar1s<Double>(L.Val / R, 1./R, &L));
01711   }
01712 
01713  template<typename Double>
01714  ADvari<Double>&
01715 operator/(Double L, const ADvari<Double> &R) {
01716   Double recip = 1. / R.Val;
01717   Double q = L * recip;
01718   Double d1 = -q*recip;
01719   return *(new ADvar1g<Double>(q, d1, -q*d1, &R));
01720   }
01721 
01722  template<typename Double>
01723  ADvar<Double>&
01724 ADvar<Double>::operator/=(Double R) {
01725   ADVari *Lcv = this->cv;
01726   this->cv = new ADvar1s<Double>(Lcv->Val / R, 1./R, Lcv RAD_ACA);
01727   return *this;
01728   }
01729 
01730  template<typename Double>
01731  ADvari<Double>&
01732 acos(const ADvari<Double> &v) {
01733   Double t = v.Val, t1 = 1. - t*t, d1 = -1./std::sqrt(t1);
01734   return *(new ADvar1g<Double>(std::acos(t), d1, t*d1/t1, &v));
01735   }
01736 
01737  template<typename Double>
01738  ADvari<Double>&
01739 acosh(const ADvari<Double> &v) {
01740   Double d1, t, t1, t2;
01741   t = v.Val;
01742   t1 = std::sqrt(t2 = t*t - 1.);
01743   d1 = 1. / t1;
01744   return *(new ADvar1g<Double>(std::log(t + t1), d1, -t*d1/t2, &v));
01745   }
01746 
01747  template<typename Double>
01748  ADvari<Double>&
01749 asin(const ADvari<Double> &v) {
01750   Double d1, t, t1;
01751   t = v.Val;
01752   d1 = 1. / std::sqrt(t1 = 1. - t*t);
01753   return *(new ADvar1g<Double>(std::asin(t), d1, t*d1/t1, &v));
01754   }
01755 
01756  template<typename Double>
01757  ADvari<Double>&
01758 asinh(const ADvari<Double> &v) {
01759   Double d1, t, t1, t2, td;
01760   t = v.Val;
01761   t1 = std::sqrt(t2 = t*t + 1.);
01762   d1 = 1. / t1;
01763   td = 1.;
01764   if (t < 0.)
01765     td = -1.;
01766   return *(new ADvar1g<Double>(td*std::log(t*td + t1), d1, -(t/t2)*d1, &v));
01767   }
01768 
01769  template<typename Double>
01770  ADvari<Double>&
01771 atan(const ADvari<Double> &v) {
01772   Double t = v.Val, d1 = 1./(1. + t*t);
01773   return *(new ADvar1g<Double>(std::atan(t), d1, -(t+t)*d1*d1, &v));
01774   }
01775 
01776  template<typename Double>
01777  ADvari<Double>&
01778 atanh(const ADvari<Double> &v) {
01779   Double t = v.Val, d1 = 1./(1. - t*t);
01780   return *(new ADvar1g<Double>(0.5*std::log((1.+t)/(1.-t)), d1, (t+t)*d1*d1, &v));
01781   }
01782 
01783  template<typename Double>
01784  ADvari<Double>&
01785 atan2(const ADvari<Double> &L, const ADvari<Double> &R) {
01786   Double R2, t, t2, x, x2, y, y2;
01787   x = L.Val;
01788   y = R.Val;
01789   x2 = x*x;
01790   y2 = y*y;
01791   t = 1./(x2 + y2);
01792   t2 = t*t;
01793   R2 = 2.*t2*x*y;
01794   return *(new ADvar2g<Double>(std::atan2(x,y), y*t, -x*t, -R2, t2*(x2 - y2), R2, &L, &R));
01795   }
01796 
01797  template<typename Double>
01798  ADvari<Double>&
01799 atan2(Double x, const ADvari<Double> &R) {
01800   Double t, x2, y, y2;
01801   y = R.Val;
01802   x2 = x*x;
01803   y2 = y*y;
01804   t = 1./(x2 + y2);
01805   return *(new ADvar1g<Double>(std::atan2(x,y), -x*t, 2.*t*t*x*y, &R));
01806   }
01807 
01808  template<typename Double>
01809  ADvari<Double>&
01810 atan2(const ADvari<Double> &L, Double y) {
01811   Double t, x, x2, y2;
01812   x = L.Val;
01813   x2 = x*x;
01814   y2 = y*y;
01815   t = 1./(x2 + y2);
01816   return *(new ADvar1g<Double>(std::atan2(x,y), y*t, -2.*t*t*x*y, &L));
01817   }
01818 
01819  template<typename Double>
01820  ADvari<Double>&
01821 max(const ADvari<Double> &L, const ADvari<Double> &R) {
01822   const ADvari<Double> &x = L.Val >= R.Val ? L : R;
01823   return *(new ADvar1<Double>(Hv_copy, x.Val, &x.adc.One, &x));
01824   }
01825 
01826  template<typename Double>
01827  ADvari<Double>&
01828 max(Double L, const ADvari<Double> &R) {
01829   if (L >= R.Val)
01830     return *(new ADvari<Double>(Hv_const, L));
01831   return *(new ADvar1<Double>(Hv_copy, R.Val, &R.adc.One, &R));
01832   }
01833 
01834  template<typename Double>
01835  ADvari<Double>&
01836 max(const ADvari<Double> &L, Double R) {
01837   if (L.Val >= R)
01838     return *(new ADvar1<Double>(Hv_copy, L.Val, &L.adc.One, &L));
01839   return *(new ADvari<Double>(Hv_const, R));
01840   }
01841 
01842  template<typename Double>
01843  ADvari<Double>&
01844 min(const ADvari<Double> &L, const ADvari<Double> &R) {
01845   const ADvari<Double> &x = L.Val <= R.Val ? L : R;
01846   return *(new ADvar1<Double>(Hv_copy, x.Val, &x.adc.One, &x));
01847   }
01848 
01849  template<typename Double>
01850  ADvari<Double>&
01851 min(Double L, const ADvari<Double> &R) {
01852   if (L <= R.Val)
01853     return *(new ADvari<Double>(Hv_const, L));
01854   return *(new ADvar1<Double>(Hv_copy, R.Val, &R.adc.One, &R));
01855   }
01856 
01857  template<typename Double>
01858  ADvari<Double>&
01859 min(const ADvari<Double> &L, Double R) {
01860   if (L.Val <= R)
01861     return *(new ADvar1<Double>(Hv_copy, L.Val, &L.adc.One, &L));
01862   return *(new ADvari<Double>(Hv_const, R));
01863   }
01864 
01865  template<typename Double>
01866  ADvari<Double>&
01867 cos(const ADvari<Double> &v) {
01868   Double t = std::cos(v.Val);
01869   return *(new ADvar1g<Double>(t, -std::sin(v.Val), -t, &v));
01870   }
01871 
01872  template<typename Double>
01873  ADvari<Double>&
01874 cosh(const ADvari<Double> &v) {
01875   Double t = std::cosh(v.Val);
01876   return *(new ADvar1g<Double>(t, std::sinh(v.Val), t, &v));
01877   }
01878 
01879  template<typename Double>
01880  ADvari<Double>&
01881 exp(const ADvari<Double> &v) {
01882   Double t = std::exp(v.Val);
01883   return *(new ADvar1g<Double>(t, t, t, &v));
01884   }
01885 
01886  template<typename Double>
01887  ADvari<Double>&
01888 log(const ADvari<Double> &v) {
01889   Double x = v.Val, d1 = 1. / x;
01890   return *(new ADvar1g<Double>(std::log(x), d1, -d1*d1, &v));
01891   }
01892 
01893  template<typename Double>
01894  ADvari<Double>&
01895 log10(const ADvari<Double> &v) {
01896   static double num = 1. / std::log(10.);
01897   Double d1, t, x;
01898   x = v.Val;
01899   t = 1. / x;
01900   d1 = num * t;
01901   return *(new ADvar1g<Double>(std::log10(x), d1, -d1*t, &v));
01902   }
01903 
01904  template<typename Double>
01905  ADvari<Double>&
01906 pow(const ADvari<Double> &L, const ADvari<Double> &R) {
01907   Double dx, dy, t, x, xlog, xym1, y;
01908   x = L.Val;
01909   y = R.Val;
01910   t = std::pow(x,y);
01911   xym1 = t / x;
01912   xlog = std::log(x);
01913   dx = y*xym1;
01914   dy = t * xlog;
01915   return *(new ADvar2g<Double>(t, dx, dy, (y-1.)*dx/x, xym1*(1. + y*xlog), dy*xlog, &L, &R));
01916   }
01917 
01918  template<typename Double>
01919  ADvari<Double>&
01920 pow(Double x, const ADvari<Double> &R) {
01921   Double dy, t, xlog, y;
01922   y = R.Val;
01923   t = std::pow(x,y);
01924   xlog = std::log(x);
01925   dy = t * xlog;
01926   return *(new ADvar1g<Double>(t, dy, dy*xlog, &R));
01927   }
01928 
01929  template<typename Double>
01930  ADvari<Double>&
01931 pow(const ADvari<Double> &L, Double y) {
01932   Double dx, t, x;
01933   x = L.Val;
01934   t = std::pow(x,y);
01935   dx = y*t/x;
01936   return *(new ADvar1g<Double>(t, dx, (y-1.)*dx/x, &L));
01937   }
01938 
01939  template<typename Double>
01940  ADvari<Double>&
01941 sin(const ADvari<Double> &v) {
01942   Double t = std::sin(v.Val);
01943   return *(new ADvar1g<Double>(t, std::cos(v.Val), -t, &v));
01944   }
01945 
01946  template<typename Double>
01947  ADvari<Double>&
01948 sinh(const ADvari<Double> &v) {
01949   Double t = std::sinh(v.Val);
01950   return *(new ADvar1g<Double>(t, std::cosh(v.Val), t, &v));
01951   }
01952 
01953  template<typename Double>
01954  ADvari<Double>&
01955 sqrt(const ADvari<Double> &v) {
01956   Double t = std::sqrt(v.Val);
01957   Double d1 = 0.5 / t;
01958   return *(new ADvar1g<Double>(t, d1, -0.5*d1/v.Val, &v));
01959   }
01960 
01961  template<typename Double>
01962  ADvari<Double>&
01963 tan(const ADvari<Double> &v) {
01964   Double d1, rv, t;
01965   rv = std::tan(v.Val);
01966   t = 1. / std::cos(v.Val);
01967   d1 = t*t;
01968   return *(new ADvar1g<Double>(rv, d1, (rv+rv)*d1, &v));
01969   }
01970 
01971  template<typename Double>
01972  ADvari<Double>&
01973 tanh(const ADvari<Double> &v) {
01974   Double d1, rv, t;
01975   rv = std::tanh(v.Val);
01976   t = 1. / std::cosh(v.Val);
01977   d1 = t*t;
01978   return *(new ADvar1g<Double>(rv, d1, -(rv+rv)*d1, &v));
01979   }
01980 
01981  template<typename Double>
01982  ADvari<Double>&
01983 abs(const ADvari<Double> &v) {
01984   Double t, p;
01985   p = 1.;
01986   if ((t = v.Val) < 0) {
01987     t = -t;
01988     p = -p;
01989     }
01990   return *(new ADvar1g<Double>(t, p, 0., &v));
01991   }
01992 
01993  template<typename Double>
01994  ADvari<Double>&
01995 fabs(const ADvari<Double> &v) { // Synonym for "abs"
01996         // "fabs" is not the best choice of name,
01997         // but this name is used at Sandia.
01998   Double t, p;
01999   p = 1.;
02000   if ((t = v.Val) < 0) {
02001     t = -t;
02002     p = -p;
02003     }
02004   return *(new ADvar1g<Double>(t, p, 0., &v));
02005   }
02006 
02007  template<typename Double>
02008  ADvari<Double>&
02009 ADf1(Double f, Double g, Double h, const ADvari<Double> &x) {
02010   return *(new ADvar1g<Double>(f, g, h, &x));
02011   }
02012 
02013  template<typename Double>
02014  inline ADvari<Double>&
02015 ADf1(Double f, Double g, Double h, const IndepADvar<Double> &x) {
02016   return *(new ADvar1g<Double>(f, g, h, x.cv));
02017   }
02018 
02019  template<typename Double>
02020  ADvari<Double>&
02021 ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
02022     const ADvari<Double> &x, const ADvari<Double> &y) {
02023   return *(new ADvar2g<Double>(f, gx, gy, hxx, hxy, hyy, &x, &y));
02024   }
02025 
02026  template<typename Double>
02027  ADvari<Double>&
02028 ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
02029     const ADvari<Double> &x, const IndepADvar<Double> &y) {
02030   return *(new ADvar2g<Double>(f, gx, gy, hxx, hxy, hyy, &x, y.cv));
02031   }
02032 
02033  template<typename Double>
02034  ADvari<Double>&
02035 ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
02036     const IndepADvar<Double> &x, const ADvari<Double> &y) {
02037   return *(new ADvar2g<Double>(f, gx, gy, hxx, hxy, hyy, x.cv, &y));
02038   }
02039 
02040  template<typename Double>
02041  ADvari<Double>&
02042 ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
02043     const IndepADvar<Double> &x, const IndepADvar<Double> &y) {
02044   return *(new ADvar2g<Double>(f, gx, gy, hxx, hxy, hyy, x.cv, y.cv));
02045   }
02046 
02047  template<typename Double>
02048  ADvari<Double>&
02049 ADfn(Double f, int n, const IndepADvar<Double> *x, const Double *g, const Double *h) {
02050   return *(new ADvarn<Double>(f, n, x, g, h));
02051   }
02052 
02053  template<typename Double>
02054  inline ADvari<Double>&
02055 ADfn(Double f, int n, const ADvar<Double> *x, const Double *g, const Double *h) {
02056   return ADfn<Double>(f, n, (IndepADvar<Double>*)x, g, h);
02057   }
02058 
02059  template<typename Double>
02060  void
02061 ADcontext<Double>::Hvprod(int n, ADvar<Double> **x, Double *v, Double *hv)
02062 {
02063   ADVari *a, *aL, *aR, **ap, **ape;
02064   ADVari_block *b, *b0;
02065   DErp *d;
02066   Double aO, adO, *g, *h, *h0, t, tL, tR;
02067   int i, j, k, m;
02068   for(i = 0; i < n; i++) {
02069     a = x[i]->cv;
02070     a->dO = v[i];
02071     a->aO = a->adO = 0.;
02072     }
02073   ADVari::adc.Aibusy->limit = ADVari::adc.Ainext;
02074   a = 0;
02075   for(b0 = 0, b = ADVari::adc.AiFirst; b; b0 = b, b = b->next) {
02076     ap = b->pADvari;
02077     ape = b->limit;
02078     while(ap < ape) {
02079       a = *ap++;
02080       a->aO = a->adO = 0.;
02081       switch(a->opclass) {
02082        case Hv_copy:
02083         a->dO = ((ADVar1*)a)->d.c->dO;
02084         break;
02085        case Hv_binary:
02086         a->dO =   ((ADVar2g*)a)->pL * ((ADVar2g*)a)->dL.c->dO
02087           + ((ADVar2g*)a)->pR * ((ADVar2g*)a)->dR.c->dO;
02088         break;
02089        case Hv_unary:
02090         a->dO = ((ADVar1g*)a)->pL * ((ADVar1g*)a)->d.c->dO;
02091         break;
02092        case Hv_negate:
02093         a->dO = -((ADVar1*)a)->d.c->dO;
02094         break;
02095        case Hv_plusLR:
02096         a->dO = ((ADVar2*)a)->dL.c->dO + ((ADVar2*)a)->dR.c->dO;
02097         break;
02098        case Hv_minusLR:
02099         a->dO = ((ADVar2*)a)->dL.c->dO - ((ADVar2*)a)->dR.c->dO;
02100         break;
02101        case Hv_timesL:
02102         a->dO = ((ADVar1s*)a)->pL * ((ADVar1s*)a)->d.c->dO;
02103         break;
02104        case Hv_timesLR:
02105         a->dO =   ((ADVar2*)a)->dR.c->Val * ((ADVar2*)a)->dL.c->dO
02106           + ((ADVar2*)a)->dL.c->Val * ((ADVar2*)a)->dR.c->dO;
02107         break;
02108        case Hv_quotLR:
02109         a->dO =   ((ADVar2q*)a)->pL * ((ADVar2q*)a)->dL.c->dO
02110           + ((ADVar2q*)a)->pR * ((ADVar2q*)a)->dR.c->dO;
02111         break;
02112        case Hv_nary:
02113         d = ((ADVarn*)a)->D;
02114         m = ((ADVarn*)a)->n;
02115         g = ((ADVarn*)a)->G;
02116         t = 0.;
02117         for(i = 0; i < m; i++)
02118           t += g[i] * d[i].c->dO;
02119         a->dO = t;
02120        }
02121       }
02122     }
02123   if (a)
02124     a->adO = 1.;
02125   for(b = b0; b; b = b->prev) {
02126     ape = b->pADvari;
02127     ap = b->limit;
02128     while(ap > ape) {
02129       a = *--ap;
02130       aO = a->aO;
02131       adO = a->adO;
02132       switch(a->opclass) {
02133        case Hv_copy:
02134         aL = ((ADVar1*)a)->d.c;
02135         aL->aO += aO;
02136         aL->adO += adO;
02137         break;
02138        case Hv_binary:
02139         aL = ((ADVar2g*)a)->dL.c;
02140         aR = ((ADVar2g*)a)->dR.c;
02141         tL = adO*aL->dO;
02142         tR = adO*aR->dO;
02143         aL->aO += aO*((ADVar2g*)a)->pL
02144           + tL*((ADVar2g*)a)->pL2
02145           + tR*((ADVar2g*)a)->pLR;
02146         aR->aO += aO*((ADVar2g*)a)->pR
02147           + tL*((ADVar2g*)a)->pLR
02148           + tR*((ADVar2g*)a)->pR2;
02149         aL->adO += adO * ((ADVar2g*)a)->pL;
02150         aR->adO += adO * ((ADVar2g*)a)->pR;
02151         break;
02152        case Hv_unary:
02153         aL = ((ADVar1g*)a)->d.c;
02154         aL->aO += aO * ((ADVar1g*)a)->pL
02155           + adO * aL->dO * ((ADVar1g*)a)->pL2;
02156         aL->adO += adO * ((ADVar1g*)a)->pL;
02157         break;
02158        case Hv_negate:
02159         aL = ((ADVar1*)a)->d.c;
02160         aL->aO -= aO;
02161         aL->adO -= adO;
02162         break;
02163        case Hv_plusLR:
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_minusLR:
02172         aL = ((ADVar2*)a)->dL.c;
02173         aR = ((ADVar2*)a)->dR.c;
02174         aL->aO += aO;
02175         aL->adO += adO;
02176         aR->aO -= aO;
02177         aR->adO -= adO;
02178         break;
02179        case Hv_timesL:
02180         aL = ((ADVar1s*)a)->d.c;
02181         aL->aO += aO * (tL = ((ADVar1s*)a)->pL);
02182         aL->adO += adO * tL;
02183         break;
02184        case Hv_timesLR:
02185         aL = ((ADVar2*)a)->dL.c;
02186         aR = ((ADVar2*)a)->dR.c;
02187         aL->aO += aO * (tL = aR->Val) + adO*aR->dO;
02188         aR->aO += aO * (tR = aL->Val) + adO*aL->dO;
02189         aL->adO += adO * tL;
02190         aR->adO += adO * tR;
02191         break;
02192        case Hv_quotLR:
02193         aL = ((ADVar2q*)a)->dL.c;
02194         aR = ((ADVar2q*)a)->dR.c;
02195         tL = adO*aL->dO;
02196         tR = adO*aR->dO;
02197         aL->aO += aO*((ADVar2q*)a)->pL
02198           + tR*((ADVar2q*)a)->pLR;
02199         aR->aO += aO*((ADVar2q*)a)->pR
02200           + tL*((ADVar2q*)a)->pLR
02201           + tR*((ADVar2q*)a)->pR2;
02202         aL->adO += adO * ((ADVar2q*)a)->pL;
02203         aR->adO += adO * ((ADVar2q*)a)->pR;
02204         break;
02205        case Hv_nary:
02206         d  = ((ADVarn*)a)->D;
02207         m  = ((ADVarn*)a)->n;
02208         g  = ((ADVarn*)a)->G;
02209         h0 = ((ADVarn*)a)->H;
02210         for(i = 0; i < m; i++) {
02211           aL = d[i].c;
02212           aL->adO += adO * (t = g[i]);
02213           aL->aO += t*aO;
02214           t = adO * aL->dO;
02215           for(h = h0, j = 0; j <= i; j++)
02216             d[j].c->aO += t * *h++;
02217           h0 = h--;
02218           for(k = j; j < m; j++)
02219             d[j].c->aO += t * *(h += k++);
02220           }
02221        }
02222       }
02223     }
02224   for(i = 0; i < n; i++) {
02225     a = x[i]->cv;
02226     a->dO = 0.;
02227     hv[i] = a->aO;
02228     }
02229   }
02230 
02231  template<typename Double>
02232  inline Double
02233 val(const ADvari<Double> &x) {
02234   return x.Val;
02235   }
02236 
02237 #undef RAD_ACA
02238 #define A (ADvari<Double>*)
02239 #ifdef RAD_Const_WARN
02240 #define C(x) (((x)->opno < 0) ? RAD_Const_Warn(x) : 0, *A x)
02241 #else
02242 #define C(x) *A x
02243 #endif
02244 #define T template<typename Double> inline
02245 #define F ADvari<Double>&
02246 #define Ai const ADvari<Double>&
02247 #define AI const IndepADvar<Double>&
02248 #define D Double
02249 #define T2(r,f) \
02250  T r f(Ai L, AI R) { return f(L, C(R.cv)); }\
02251  T r f(AI L, Ai R) { return f(C(L.cv), R); }\
02252  T r f(AI L, AI R) { return f(C(L.cv), C(R.cv)); }\
02253  T r f(AI L, D R) { return f(C(L.cv), R); }\
02254  T r f(Ai L, Dtype R) { return f(L, (D)R); }\
02255  T r f(AI L, Dtype R) { return f(C(L.cv), (D)R); }\
02256  T r f(Ai L, long R) { return f(L, (D)R); }\
02257  T r f(AI L, long R) { return f(C(L.cv), (D)R); }\
02258  T r f(Ai L, int R) { return f(L, (D)R); }\
02259  T r f(AI L, int R) { return f(C(L.cv), (D)R); }\
02260  T r f(D L, AI R) { return f(L, C(R.cv)); }\
02261  T r f(Dtype L, Ai R) { return f((D)L, R); }\
02262  T r f(Dtype L, AI R) { return f((D)L, C(R.cv)); }\
02263  T r f(long L, Ai R) { return f((D)L, R); }\
02264  T r f(long L, AI R) { return f((D)L, C(R.cv)); }\
02265  T r f(int L, Ai R) { return f((D)L, R); }\
02266  T r f(int L, AI R) { return f((D)L, C(R.cv)); }
02267 
02268 T2(F, operator+)
02269 T2(F, operator-)
02270 T2(F, operator*)
02271 T2(F, operator/)
02272 T2(F, atan2)
02273 T2(F, pow)
02274 T2(F, max)
02275 T2(F, min)
02276 T2(int, operator<)
02277 T2(int, operator<=)
02278 T2(int, operator==)
02279 T2(int, operator!=)
02280 T2(int, operator>=)
02281 T2(int, operator>)
02282 
02283 #undef T2
02284 #undef D
02285 
02286 #define T1(f)\
02287  T F f(AI x) { return f(C(x.cv)); }
02288 
02289 T1(operator+)
02290 T1(operator-)
02291 T1(abs)
02292 T1(acos)
02293 T1(acosh)
02294 T1(asin)
02295 T1(asinh)
02296 T1(atan)
02297 T1(atanh)
02298 T1(cos)
02299 T1(cosh)
02300 T1(exp)
02301 T1(log)
02302 T1(log10)
02303 T1(sin)
02304 T1(sinh)
02305 T1(sqrt)
02306 T1(tan)
02307 T1(tanh)
02308 T1(fabs)
02309 
02310 T F copy(AI x)
02311 { return *(new ADvar1<Double>(Hv_copy, x.cv->Val, &ADcontext<Double>::One, (ADvari<Double>*)x.cv)); }
02312 
02313 T F copy(Ai x)
02314 { return *(new ADvar1<Double>(Hv_copy, x.Val, &ADcontext<Double>::One, (ADvari<Double>*)&x)); }
02315 
02316 #undef T1
02317 #undef AI
02318 #undef Ai
02319 #undef F
02320 #undef T
02321 #undef A
02322 #undef C
02323 #undef Ttype
02324 #undef Dtype
02325 
02326 } /* namespace Rad2 */
02327 } /* namespace Sacado */
02328 
02329 #endif /* SACADO_TRAD2_H */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:19:37 2011 for Sacado Package Browser (Single Doxygen Collection) by  doxygen 1.6.3