Sacado_rad.hpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                           Sacado Package
00005 //                 Copyright (2006) 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 // RAD package (Reverse Automatic Differentiation) --
00031 // a package specialized for function and gradient evaluations.
00032 // Written in 2004 by David M. Gay at Sandia National Labs, Albuquerque, NM.
00033 
00034 #ifndef SACADO_RAD_H
00035 #define SACADO_RAD_H
00036 
00037 #include <stddef.h>
00038 #include <cmath>
00039 #include <math.h>
00040 
00041 namespace Sacado {
00042 namespace Radnt { // nontemplated RAD
00043 
00044 // -DRAD_NO_USING_STDCC is needed, e.g., with Sun CC 5.7
00045 #ifndef RAD_NO_USING_STDCC
00046   // Bring math functions into scope
00047   using std::exp;
00048   using std::log;
00049   using std::log10;
00050   using std::sqrt;
00051   using std::cos;
00052   using std::sin;
00053   using std::tan;
00054   using std::acos;
00055   using std::asin;
00056   using std::atan;
00057   using std::cosh;
00058   using std::sinh;
00059   using std::tanh;
00060   using std::abs;
00061   using std::fabs;
00062   using std::atan2;
00063   using std::pow;
00064 #endif //!RAD_NO_USING_STDCC
00065 
00066  class ADvar;
00067  class ADvari;
00068  class ADvar1;
00069  class ADvar2;
00070  class ConstADvar;
00071  class Derp;
00072  class IndepADvar;
00073 
00074  extern ADvari& ADf1(double f, double g, const ADvari &x);
00075  extern ADvari& ADf2(double f, double gx, double gy, const ADvari &x, const ADvari &y);
00076  extern ADvari& ADfn(double f, int n, const ADvar *x, const double *g);
00077 
00078  struct
00079 ADmemblock {  // We get memory in ADmemblock chunks and never give it back,
00080     // but reuse it once computations start anew after call(s) on
00081     // ADcontext::Gradcomp() or ADcontext::Weighted_Gradcomp().
00082   ADmemblock *next;
00083   double memblk[1000];
00084   };
00085 
00086  class
00087 ADcontext { // A singleton class: one instance in radops.c
00088   ADmemblock *Busy, *Free;
00089   char *Mbase;
00090   size_t Mleft;
00091   ADmemblock First;
00092   void *new_ADmemblock(size_t);
00093  public:
00094   ADcontext();
00095   void *Memalloc(size_t len);
00096   static void Gradcomp(int);
00097   static inline void Gradcomp() { Gradcomp(1); }
00098   static void Weighted_Gradcomp(int, ADvar**, double*);
00099   };
00100 
00101 inline void *ADcontext::Memalloc(size_t len) {
00102     if (Mleft >= len)
00103       return Mbase + (Mleft -= len);
00104     return new_ADmemblock(len);
00105     }
00106 
00107  class
00108 CADcontext: public ADcontext {
00109  protected:
00110   bool fpval_implies_const;
00111  public:
00112   friend class ADvar;
00113   CADcontext(): ADcontext() { fpval_implies_const = false; }
00114   static const double One, negOne;
00115   };
00116 
00117  class
00118 Derp {    // one derivative-propagation operation
00119  public:
00120   friend class ADvarn;
00121   static Derp *LastDerp;
00122   Derp *next;
00123   const double *a;
00124   const ADvari *b;
00125   const ADvari *c;
00126   void *operator new(size_t);
00127   void operator delete(void*) {} /*Should never be called.*/
00128   inline Derp(){};
00129   Derp(const ADvari *);
00130   Derp(const double *, const ADvari *);
00131   Derp(const double *, const ADvari *, const ADvari *);
00132   /* c->aval += a * b->aval; */
00133   };
00134 
00135 inline Derp::Derp(const ADvari *c1): c(c1) {
00136     next = LastDerp;
00137     LastDerp = this;
00138     }
00139 
00140 inline Derp::Derp(const double *a1, const ADvari *c1): a(a1), c(c1) {
00141     next = LastDerp;
00142     LastDerp = this;
00143     }
00144 
00145 inline Derp::Derp(const double *a1, const ADvari *b1, const ADvari *c1): a(a1), b(b1), c(c1) {
00146     next = LastDerp;
00147     LastDerp = this;
00148     }
00149 
00150  class
00151 ADvari {  // implementation of an ADvar
00152  public:
00153   double Val;   // result of this operation
00154   mutable double aval;  // adjoint -- partial of final result w.r.t. this Val
00155   void *operator new(size_t len) { return ADvari::adc.Memalloc(len); }
00156   void operator delete(void*) {} /*Should never be called.*/
00157   inline ADvari(double t): Val(t), aval(0.) {}
00158   inline ADvari(double t, double ta): Val(t), aval(ta) {}
00159   inline ADvari(): Val(0.), aval(0.) {}
00160   static ADcontext adc;
00161 #ifdef RAD_AUTO_AD_Const
00162   friend class ADcontext;
00163   friend class ADvar1;
00164   friend class ADvar;
00165   friend class ConstADvar;
00166   friend class IndepADvar;
00167  private:
00168   ADvari *Next;
00169   static ADvari *First_ADvari, **Last_ADvari;
00170  protected:
00171   IndepADvar *padv;
00172  public:
00173   ADvari(const IndepADvar *, double);
00174 #endif
00175 #define F friend
00176 #define R ADvari&
00177 #define Ai const ADvari&
00178 #define T1(r,f) F r f(Ai);
00179 #define T2(r,f) \
00180 F r f(Ai,Ai); \
00181 F r f(double,Ai); \
00182 F r f(Ai,double);
00183   T1(R,operator+)
00184   T2(R,operator+)
00185   T1(R,operator-)
00186   T2(R,operator-)
00187   T2(R,operator*)
00188   T2(R,operator/)
00189   T1(R,abs)
00190   T1(R,acos)
00191   T1(R,acosh)
00192   T1(R,asin)
00193   T1(R,asinh)
00194   T1(R,atan)
00195   T1(R,atanh)
00196   T2(R,atan2)
00197   T2(R,max)
00198   T2(R,min)
00199   T1(R,cos)
00200   T1(R,cosh)
00201   T1(R,exp)
00202   T1(R,log)
00203   T1(R,log10)
00204   T2(R,pow)
00205   T1(R,sin)
00206   T1(R,sinh)
00207   T1(R,sqrt)
00208   T1(R,tan)
00209   T1(R,tanh)
00210   T1(R,fabs)
00211   T1(R,copy)
00212   T2(int,operator<)
00213   T2(int,operator<=)
00214   T2(int,operator==)
00215   T2(int,operator!=)
00216   T2(int,operator>=)
00217   T2(int,operator>)
00218 #undef T2
00219 #undef T1
00220 #undef Ai
00221 #undef R
00222 #undef F
00223 
00224   friend ADvari& ADf1(double f, double g, const ADvari &x);
00225   friend ADvari& ADf2(double f, double gx, double gy, const ADvari &x, const ADvari &y);
00226   friend ADvari& ADfn(double f, int n, const ADvar *x, const double *g);
00227   };
00228 
00229  inline void* Derp::operator new(size_t len) { return ADvari::adc.Memalloc(len); }
00230 
00231 
00232  class
00233 ADvar1: public ADvari { // simplest unary ops
00234  public:
00235   Derp d;
00236   ADvar1(double val1): ADvari(val1) {}
00237   ADvar1(double val1, const ADvari *c1): d(c1) { Val = val1; }
00238   ADvar1(double val1, const double *a1, const ADvari *c1): d(a1,this,c1) { Val = val1; }
00239 #ifdef RAD_AUTO_AD_Const
00240   ADvar1(const IndepADvar *, const IndepADvar &);
00241   ADvar1(const IndepADvar *, const ADvari &);
00242 #endif
00243   };
00244 
00245  class
00246 ConstADvari: public ADvari {
00247  private:
00248   ConstADvari *prevcad;
00249   ConstADvari() {}; // prevent construction without value (?)
00250   static ConstADvari *lastcad;
00251  public:
00252   static CADcontext cadc;
00253   inline void *operator new(size_t len) { return ConstADvari::cadc.Memalloc(len); }
00254   inline ConstADvari(double t): ADvari(t) { prevcad = lastcad; lastcad = this; }
00255   static void aval_reset(void);
00256   };
00257 
00258  class
00259 IndepADvar
00260 {
00261  private:
00262   inline IndepADvar& operator=(const IndepADvar &x) {
00263     /* private to prevent assignment */
00264 #ifdef RAD_AUTO_AD_Const
00265     if (cv)
00266       cv->padv = 0;
00267     cv = new ADvar1(this,x);
00268     return *this;
00269 #else
00270 #ifdef RAD_EQ_ALIAS
00271     cv = x.cv;
00272     return *this;
00273 #else
00274     return ADvar_operatoreq(this,*x.cv);
00275 #endif
00276 #endif /* RAD_AUTO_AD_Const */
00277     }
00278  protected:
00279   static void AD_Const(const IndepADvar&);
00280   ADvari *cv;
00281  public:
00282   typedef double value_type;
00283   friend class ADvar;
00284   friend class ADvar1;
00285   friend class ADvarn;
00286   friend class ADcontext;
00287   IndepADvar(double);
00288   IndepADvar(int);
00289   IndepADvar(long);
00290   IndepADvar& operator=(double);
00291 #ifdef RAD_AUTO_AD_Const
00292   inline IndepADvar(const IndepADvar &x) { cv = x.cv ? new ADvar1(this, x) : 0; };
00293   inline IndepADvar() { cv = 0; }
00294   inline ~IndepADvar() {
00295       if (cv)
00296         cv->padv = 0;
00297     }
00298 #else
00299   inline IndepADvar() {
00300 #ifndef RAD_EQ_ALIAS
00301     cv = 0;
00302 #endif
00303     }
00304   inline ~IndepADvar() {}
00305   friend IndepADvar& ADvar_operatoreq(IndepADvar*, const ADvari&);
00306 #endif
00307 
00308   friend void AD_Const(const IndepADvar&);
00309 
00310   inline operator ADvari&() { return *cv; }
00311   inline operator ADvari&() const { return *(ADvari*)cv; }
00312 
00313   inline double val() const { return cv->Val; }
00314   inline double adj() const { return cv->aval; }
00315   static inline void Gradcomp(int wantgrad)
00316     { ADcontext::Gradcomp(wantgrad); }
00317   static inline void Gradcomp()
00318     { ADcontext::Gradcomp(1); }
00319   static inline void aval_reset() { ConstADvari::aval_reset(); }
00320   static inline void Weighted_Gradcomp(int n, ADvar **v, double *w)
00321         { ADcontext::Weighted_Gradcomp(n, v, w); }
00322 
00323 
00324 #define Ai const ADvari&
00325 #define AI const IndepADvar&
00326 #define T2(r,f) \
00327  r f(AI,AI);\
00328  r f(Ai,AI);\
00329  r f(AI,Ai);\
00330  r f(double,AI);\
00331  r f(AI,double);
00332 #define T1(f) friend ADvari& f(AI);
00333 
00334 #define F friend ADvari&
00335 T2(F, operator+)
00336 T2(F, operator-)
00337 T2(F, operator*)
00338 T2(F, operator/)
00339 T2(F, atan2)
00340 T2(F, max)
00341 T2(F, min)
00342 T2(F, pow)
00343 #undef F
00344 #define F friend int
00345 T2(F, operator<)
00346 T2(F, operator<=)
00347 T2(F, operator==)
00348 T2(F, operator!=)
00349 T2(F, operator>=)
00350 T2(F, operator>)
00351 
00352 T1(operator+)
00353 T1(operator-)
00354 T1(abs)
00355 T1(acos)
00356 T1(acosh)
00357 T1(asin)
00358 T1(asinh)
00359 T1(atan)
00360 T1(atanh)
00361 T1(cos)
00362 T1(cosh)
00363 T1(exp)
00364 T1(log)
00365 T1(log10)
00366 T1(sin)
00367 T1(sinh)
00368 T1(sqrt)
00369 T1(tan)
00370 T1(tanh)
00371 T1(fabs)
00372 T1(copy)
00373 
00374 #undef T1
00375 #undef T2
00376 #undef F
00377 #undef Ai
00378 #undef AI
00379 
00380   };
00381 
00382  class
00383 ADvar: public IndepADvar {    // an "active" variable
00384   void ADvar_ctr(double d);
00385  public:
00386   inline ADvar() { /* cv = 0; */ }
00387   inline ADvar(double d) { ADvar_ctr(d); }
00388   inline ADvar(int i) { ADvar_ctr((double)i); }
00389   inline ADvar(long i)  { ADvar_ctr((double)i); }
00390   inline ~ADvar() {}
00391 #ifdef RAD_AUTO_AD_Const
00392   friend class ADvar1;
00393   inline ADvar(const IndepADvar &x) { cv = x.cv ? new ADvar1(this, x) : 0; }
00394   inline ADvar(ADvari &x) { cv = &x; x.padv = this; }
00395   inline ADvar& operator=(const IndepADvar &x) {
00396     if (cv)
00397       cv->padv = 0;
00398     cv = new ADvar1(this,x);
00399     return *this;
00400     }
00401   inline ADvar& operator=(const ADvari &x) {
00402     if (cv)
00403       cv->padv = 0;
00404     cv = new ADvar1(this, x);
00405     return *this;
00406     }
00407 #else
00408   friend ADvar& ADvar_operatoreq(ADvar*, const ADvari&);
00409 #ifdef RAD_EQ_ALIAS
00410   /* allow aliasing v and w after "v = w;" */
00411   inline ADvar(const IndepADvar &x) { cv = x.cv; }
00412   inline ADvar(const ADvari &x) { cv = (ADvari*)&x; }
00413   inline ADvar& operator=(const ADvari &x) { cv = (ADvari*)&x;   return *this; }
00414   inline ADvar& operator=(const IndepADvar &x)  { cv = (ADvari*)x.cv; return *this; }
00415 #else
00416   ADvar(const IndepADvar &x) { cv = x.cv ? new ADvar1(x.cv->Val, &CADcontext::One, x.cv) : 0; }
00417   ADvar(const ADvari &x) { cv = new ADvar1(x.Val, &CADcontext::One, &x); }
00418   inline ADvar& operator=(const ADvari &x) { return ADvar_operatoreq(this,x); }
00419   inline ADvar& operator=(const IndepADvar &x)    { return ADvar_operatoreq(this,*x.cv); }
00420 #endif /* RAD_EQ_ALIAS */
00421 #endif /* RAD_AUTO_AD_Const */
00422   ADvar& operator=(double);
00423   ADvar& operator+=(const ADvari&);
00424   ADvar& operator+=(double);
00425   ADvar& operator-=(const ADvari&);
00426   ADvar& operator-=(double);
00427   ADvar& operator*=(const ADvari&);
00428   ADvar& operator*=(double);
00429   ADvar& operator/=(const ADvari&);
00430   ADvar& operator/=(double);
00431   inline static bool get_fpval_implies_const(void)
00432     { return ConstADvari::cadc.fpval_implies_const; }
00433   inline static void set_fpval_implies_const(bool newval)
00434     { ConstADvari::cadc.fpval_implies_const = newval; }
00435   inline static bool setget_fpval_implies_const(bool newval) {
00436     bool oldval = ConstADvari::cadc.fpval_implies_const;
00437     ConstADvari::cadc.fpval_implies_const = newval;
00438     return oldval;
00439     }
00440   static inline void Gradcomp(int wantgrad)
00441         { ADcontext::Gradcomp(wantgrad); }
00442   static inline void Gradcomp()
00443         { ADcontext::Gradcomp(1); }
00444   static inline void aval_reset() { ConstADvari::aval_reset(); }
00445   static inline void Weighted_Gradcomp(int n, ADvar **v, double *w)
00446         { ADcontext::Weighted_Gradcomp(n, v, w); }
00447   };
00448 
00449  inline void AD_Const(const IndepADvar&v) { IndepADvar::AD_Const(v); }
00450 
00451  class
00452 ConstADvar: public ADvar {
00453  private: // disable op=
00454   ConstADvar& operator+=(const ADvari&);
00455   ConstADvar& operator+=(double);
00456   ConstADvar& operator-=(const ADvari&);
00457   ConstADvar& operator-=(double);
00458   ConstADvar& operator*=(const ADvari&);
00459   ConstADvar& operator*=(double);
00460   ConstADvar& operator/=(const ADvari&);
00461   ConstADvar& operator/=(double);
00462   void ConstADvar_ctr(double);
00463  public:
00464   inline ConstADvar(double d) { ConstADvar_ctr(d); }
00465   inline ConstADvar(int i)  { ConstADvar_ctr((double)i); }
00466   inline ConstADvar(long i) { ConstADvar_ctr((double)i); }
00467   ConstADvar(const ADvari &x);
00468 #ifdef RAD_AUTO_AD_Const
00469   ConstADvar(const IndepADvar &x) { cv = new ADvar1(this,x); }
00470 #endif
00471   inline ~ConstADvar(){}
00472 #ifdef RAD_NO_CONST_UPDATE
00473  private:
00474 #endif
00475   ConstADvar();
00476   inline ConstADvar& operator=(double d) { cv->Val = d; return *this; }
00477   inline ConstADvar& operator=(const IndepADvar& d) { cv->Val = d.val(); return *this; }
00478  };
00479 
00480  class
00481 ADvar1s: public ADvar1 { // unary ops with partial "a"
00482  public:
00483   double a;
00484   ADvar1s(double val1, double a1, const ADvari *c1): ADvar1(val1,&a,c1), a(a1) {}
00485   };
00486 
00487  class
00488 ADvar2: public ADvari { // basic binary ops
00489  public:
00490   Derp dL, dR;
00491   ADvar2(double val1): ADvari(val1) {}
00492   ADvar2(double val1, const ADvari *Lcv, const double *Lc, const ADvari *Rcv,
00493       const double *Rc): ADvari(val1) {
00494     dR.next = Derp::LastDerp;
00495     dL.next = &dR;
00496     Derp::LastDerp = &dL;
00497     dL.a = Lc;
00498     dL.c = Lcv;
00499     dR.a = Rc;
00500     dR.c = Rcv;
00501     dL.b = dR.b = this;
00502     }
00503   };
00504 
00505  class
00506 ADvar2q: public ADvar2 { // binary ops with partials "a", "b"
00507  public:
00508   double a, b;
00509   ADvar2q(double val1, double Lp, double Rp, const ADvari *Lcv, const ADvari *Rcv):
00510       ADvar2(val1), a(Lp), b(Rp) {
00511     dR.next = Derp::LastDerp;
00512     dL.next = &dR;
00513     Derp::LastDerp = &dL;
00514     dL.a = &a;
00515     dL.c = Lcv;
00516     dR.a = &b;
00517     dR.c = Rcv;
00518     dL.b = dR.b = this;
00519     }
00520   };
00521 
00522  class
00523 ADvarn: public ADvari { // n-ary ops with partials "a"
00524  public:
00525   int n;
00526   double *a;
00527   Derp *Da;
00528   ADvarn(double val1, int n1, const ADvar *x, const double *g);
00529   };
00530 
00531 inline ADvari &operator+(ADvari &T) { return T; }
00532 inline ADvari &operator+(const ADvari &T) { return (ADvari&) T; }
00533 
00534 inline int operator<(const ADvari &L, const ADvari &R) { return L.Val < R.Val; }
00535 inline int operator<(const ADvari &L, double R) { return L.Val < R; }
00536 inline int operator<(double L, const ADvari &R) { return L < R.Val; }
00537 
00538 inline int operator<=(const ADvari &L, const ADvari &R) { return L.Val <= R.Val; }
00539 inline int operator<=(const ADvari &L, double R) { return L.Val <= R; }
00540 inline int operator<=(double L, const ADvari &R) { return L <= R.Val; }
00541 
00542 inline int operator==(const ADvari &L, const ADvari &R) { return L.Val == R.Val; }
00543 inline int operator==(const ADvari &L, double R) { return L.Val == R; }
00544 inline int operator==(double L, const ADvari &R) { return L == R.Val; }
00545 
00546 inline int operator!=(const ADvari &L, const ADvari &R) { return L.Val != R.Val; }
00547 inline int operator!=(const ADvari &L, double R) { return L.Val != R; }
00548 inline int operator!=(double L, const ADvari &R) { return L != R.Val; }
00549 
00550 inline int operator>=(const ADvari &L, const ADvari &R) { return L.Val >= R.Val; }
00551 inline int operator>=(const ADvari &L, double R) { return L.Val >= R; }
00552 inline int operator>=(double L, const ADvari &R) { return L >= R.Val; }
00553 
00554 inline int operator>(const ADvari &L, const ADvari &R) { return L.Val > R.Val; }
00555 inline int operator>(const ADvari &L, double R) { return L.Val > R; }
00556 inline int operator>(double L, const ADvari &R) { return L > R.Val; }
00557 
00558 inline ADvari& copy(const IndepADvar &x)
00559 { return *(new ADvar1(x.cv->Val, &CADcontext::One, x.cv)); }
00560 
00561 inline ADvari& copy(const ADvari &x)
00562 { return *(new ADvar1(x.Val, &CADcontext::One, &x)); }
00563 
00564 inline ADvari& abs(const ADvari &x)
00565 { return fabs(x); }
00566 
00567 #define A (ADvari*)
00568 #define T inline
00569 #define F ADvari&
00570 #define Ai const ADvari&
00571 #define AI const IndepADvar&
00572 #define D double
00573 #define T2(r,f) \
00574  T r f(Ai L, AI R) { return f(L, *A R.cv); }\
00575  T r f(AI L, Ai R) { return f(*A L.cv, R); }\
00576  T r f(AI L, AI R) { return f(*A L.cv, *A R.cv); }\
00577  T r f(AI L, D R) { return f(*A L.cv, R); }\
00578  T r f(D L, AI R) { return f(L, *A R.cv); }
00579 
00580 T2(F, operator+)
00581 T2(F, operator-)
00582 T2(F, operator*)
00583 T2(F, operator/)
00584 T2(F, atan2)
00585 T2(F, pow)
00586 T2(F, max)
00587 T2(F, min)
00588 T2(int, operator<)
00589 T2(int, operator<=)
00590 T2(int, operator==)
00591 T2(int, operator!=)
00592 T2(int, operator>=)
00593 T2(int, operator>)
00594 
00595 #undef T2
00596 #undef D
00597 
00598 #define T1(f)\
00599  T F f(AI x) { return f(*A x.cv); }
00600 
00601 T1(operator+)
00602 T1(operator-)
00603 T1(abs)
00604 T1(acos)
00605 T1(acosh)
00606 T1(asin)
00607 T1(asinh)
00608 T1(atan)
00609 T1(atanh)
00610 T1(cos)
00611 T1(cosh)
00612 T1(exp)
00613 T1(log)
00614 T1(log10)
00615 T1(sin)
00616 T1(sinh)
00617 T1(sqrt)
00618 T1(tan)
00619 T1(tanh)
00620 T1(fabs)
00621 
00622 #undef T1
00623 #undef AI
00624 #undef Ai
00625 #undef F
00626 #undef T
00627 #undef A
00628 
00629 } // namespace Radnt
00630 } // namespace Sacado
00631 #endif /* SACADO_RAD_H */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:19:34 2011 for Sacado Package Browser (Single Doxygen Collection) by  doxygen 1.6.3