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