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