Sacado Package Browser (Single Doxygen Collection) Version of the Day
Sacado_radops2.cpp
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 // Support routines for Hessian-vector products via the RAD package
00031 // (Reverse Automatic Differentiation) -- a package specialized for
00032 // function and gradient evaluations, and herewith extended for Hessian-
00033 // vector products.  This specialization is for several Hessian-vector
00034 // products at one point, where the vectors are determined on the fly,
00035 // e.g., by the conjugate-gradient algorithm.  Where the vectors are known
00036 // in advance, nestings of RAD and FAD might be more efficient.
00037 // RAD Hessian-vector product extension written in 2007 by David M. Gay at
00038 // Sandia National Labs, Albuquerque, NM.
00039 
00040 #include "Sacado_rad2.hpp"
00041 
00042 #ifndef SACADO_NO_NAMESPACE
00043 namespace Sacado {
00044 namespace Rad2d { // "2" for 2nd derivatives, "d" for "double"
00045 #endif
00046 
00047 #ifdef RAD_AUTO_AD_Const
00048 ADvari *ADvari::First_ADvari, **ADvari::Last_ADvari = &ADvari::First_ADvari;
00049 #undef RAD_DEBUG_BLOCKKEEP
00050 #else
00051 #ifdef RAD_DEBUG_BLOCKKEEP
00052 #if !(RAD_DEBUG_BLOCKKEEP > 0)
00053 #undef RAD_DEBUG_BLOCKKEEP
00054 #else
00055 extern "C" void _uninit_f2c(void *x, int type, long len);
00056 #define TYDREAL 5
00057 static ADmemblock *rad_Oldcurmb;
00058 static int rad_busy_blocks;
00059 #endif /*RAD_DEBUG_BLOCKKEEP > 0*/
00060 #endif /*RAD_DEBUG_BLOCKKEEP*/
00061 #endif /*RAD_AUTO_AD_Const*/
00062 
00063 Derp *Derp::LastDerp = 0;
00064 ADcontext ADvari::adc;
00065 CADcontext ConstADvari::cadc;
00066 ConstADvari *ConstADvari::lastcad;
00067 static int rad_need_reinit;
00068 #ifdef RAD_DEBUG_BLOCKKEEP
00069 static size_t rad_mleft_save;
00070 #endif
00071 
00072 const double CADcontext::One = 1., CADcontext::negOne = -1.;
00073 
00074 ADcontext::ADcontext()
00075 {
00076   First.next = 0;
00077   Busy = &First;
00078   Free = 0;
00079   Mbase = (char*)First.memblk;
00080   Mleft = sizeof(First.memblk);
00081   Aibusy = &AiFirst;
00082   Aifree = 0;
00083   Ainext = AiFirst.pADvari;
00084   AiFirst.next = AiFirst.prev = 0;
00085   AiFirst.limit = Ailimit = AiFirst.pADvari + ADvari_block::Gulp;
00086   }
00087 
00088  void*
00089 ADcontext::new_ADmemblock(size_t len)
00090 {
00091   ADmemblock *mb, *mb0, *mb1, *mbf, *x;
00092   ADvari_block *b;
00093 #ifdef RAD_AUTO_AD_Const
00094   ADvari *a, *anext;
00095   IndepADvar *v;
00096 #endif /* RAD_AUTO_AD_Const */
00097 
00098   if (rad_need_reinit && this == &ADvari::adc) {
00099     rad_need_reinit = 0;
00100     Derp::LastDerp = 0;
00101     Aibusy = b = &AiFirst;
00102     Aifree = b->next;
00103     b->next = b->prev = 0;
00104     Ailimit = b->limit = (Ainext = b->pADvari) + ADvari_block::Gulp;
00105 #ifdef RAD_DEBUG_BLOCKKEEP
00106     Mleft = rad_mleft_save;
00107     if (Mleft < sizeof(First.memblk))
00108       _uninit_f2c(Mbase + Mleft, TYDREAL,
00109        (sizeof(First.memblk) - Mleft)/sizeof(double));
00110     if ((mb = Busy->next)) {
00111       if (!(mb0 = rad_Oldcurmb))
00112         mb0 = &First;
00113       for(;; mb = mb->next) {
00114         _uninit_f2c(mb->memblk, TYDREAL,
00115           sizeof(First.memblk)/sizeof(double));
00116         if (mb == mb0)
00117           break;
00118         }
00119       }
00120     rad_Oldcurmb = Busy;
00121     if (rad_busy_blocks >= RAD_DEBUG_BLOCKKEEP) {
00122       rad_busy_blocks = 0;
00123       rad_Oldcurmb = 0;
00124       mb0 = &First;
00125       mbf =  Free;
00126       for(mb = Busy; mb != mb0; mb = mb1) {
00127         mb1 = mb->next;
00128         mb->next = mbf;
00129         mbf = mb;
00130         }
00131       Free = mbf;
00132       Busy = mb;
00133       Mbase = (char*)First.memblk;
00134       Mleft = sizeof(First.memblk);
00135       }
00136 
00137 #else /* !RAD_DEBUG_BLOCKKEEP */
00138 
00139     mb0 = &ADvari::adc.First;
00140     mbf =  ADvari::adc.Free;
00141     for(mb = ADvari::adc.Busy; mb != mb0; mb = mb1) {
00142       mb1 = mb->next;
00143       mb->next = mbf;
00144       mbf = mb;
00145       }
00146     ADvari::adc.Free = mbf;
00147     ADvari::adc.Busy = mb;
00148     ADvari::adc.Mbase = (char*)ADvari::adc.First.memblk;
00149     ADvari::adc.Mleft = sizeof(ADvari::adc.First.memblk);
00150 #ifdef RAD_AUTO_AD_Const
00151     *ADvari::Last_ADvari = 0;
00152     ADvari::Last_ADvari = &ADvari::First_ADvari;
00153     if ((anext = ADvari::First_ADvari)) {
00154       while((a = anext)) {
00155         anext = a->Next;
00156         if ((v = a->padv))
00157           v->cv = new ADvari(v, a->Val);
00158         }
00159       }
00160 #endif /*RAD_AUTO_AD_Const*/
00161 #endif /* RAD_DEBUG_BLOCKKEEP */
00162     if (Mleft >= len)
00163       return Mbase + (Mleft -= len);
00164     }
00165 
00166   if ((x = Free))
00167     Free = x->next;
00168   else
00169     x = new ADmemblock;
00170 #ifdef RAD_DEBUG_BLOCKKEEP
00171   rad_busy_blocks++;
00172 #endif
00173   x->next = Busy;
00174   Busy = x;
00175   return (Mbase = (char*)x->memblk) +
00176     (Mleft = sizeof(First.memblk) - len);
00177   }
00178 
00179  void
00180 ADcontext::new_ADvari_block()
00181 {
00182   ADvari_block *ob, *nb;
00183   ob = Aibusy;
00184   ob->limit = Ailimit;  // should be unnecessary, but harmless
00185   if ((nb = Aifree))
00186     Aifree = nb->next;
00187   else
00188     nb = new ADvari_block;
00189   Aibusy = Aibusy->next = nb;
00190   nb->limit = Ailimit = (Ainext = nb->pADvari) + ADvari_block::Gulp;
00191   ob->next = nb;
00192   nb->prev = ob;
00193   nb->next = 0;
00194   }
00195 
00196  void
00197 ADcontext::Gradcomp(int wantgrad)
00198 {
00199   Derp *d;
00200 
00201   if (rad_need_reinit && wantgrad) {
00202     for(d = Derp::LastDerp; d; d = d->next)
00203       d->c->aval = 0;
00204     }
00205   else {
00206     rad_need_reinit = 1;
00207 #ifdef RAD_DEBUG_BLOCKKEEP
00208     rad_mleft_save = ADvari::adc.Mleft;
00209 #endif
00210     ADvari::adc.Mleft = 0;
00211     }
00212 
00213   if ((d = Derp::LastDerp) && wantgrad) {
00214     d->b->aval = 1;
00215     do d->c->aval += *d->a * d->b->aval;
00216     while((d = d->next));
00217     }
00218   }
00219 
00220  void
00221 ADcontext::Weighted_Gradcomp(int n, ADvar **v, double *w)
00222 {
00223   Derp *d;
00224   int i;
00225 
00226   if (rad_need_reinit) {
00227     for(d = Derp::LastDerp; d; d = d->next)
00228       d->c->aval = 0;
00229     }
00230   else {
00231     rad_need_reinit = 1;
00232 #ifdef RAD_DEBUG_BLOCKKEEP
00233     rad_mleft_save = ADvari::adc.Mleft;
00234 #endif
00235     ADvari::adc.Mleft = 0;
00236     }
00237   if ((d = Derp::LastDerp)) {
00238     for(i = 0; i < n; i++)
00239       v[i]->cv->aval = w[i];
00240     do d->c->aval += *d->a * d->b->aval;
00241     while((d = d->next));
00242     }
00243   }
00244 
00245 #ifdef RAD_AUTO_AD_Const
00246 
00247 IndepADvar::IndepADvar(double d)
00248 {
00249   cv = new ADvari(this, d);
00250   }
00251 
00252 IndepADvar::IndepADvar(int d)
00253 {
00254   cv = new ADvari(this, (double)d);
00255   }
00256 
00257 IndepADvar::IndepADvar(long d)
00258 {
00259   cv = new ADvari(this, (double)d);
00260   }
00261 
00262 #else 
00264 IndepADvar::IndepADvar(double d)
00265 {
00266   ADvari *x = new ADvari(Hv_const, d);
00267   cv = x;
00268   }
00269 
00270 IndepADvar::IndepADvar(int d)
00271 {
00272   ADvari *x = new ADvari(Hv_const, (double)d);
00273   cv = x;
00274   }
00275 
00276 IndepADvar::IndepADvar(long d)
00277 {
00278   ADvari *x = new ADvari(Hv_const, (double)d);
00279   cv = x;
00280   }
00281 
00282 #endif /*RAD_AUTO_AD_Const*/
00283 
00284  void
00285 ADvar::ADvar_ctr(double d)
00286 {
00287 #ifdef RAD_AUTO_AD_Const
00288   ADvari *x = new ADvari(this, d);
00289 #else
00290   ADvari *x = ConstADvari::cadc.fpval_implies_const
00291     ? new ConstADvari(d)
00292     : new ADvari(Hv_const, d);
00293 #endif
00294   cv = x;
00295   }
00296 
00297 ConstADvar::ConstADvar()
00298 {
00299   ConstADvari *x = new ConstADvari(0.);
00300   cv = x;
00301   }
00302 
00303  void
00304 ConstADvar::ConstADvar_ctr(double d)
00305 {
00306   ConstADvari *x = new ConstADvari(d);
00307   cv = x;
00308   }
00309 ConstADvar::ConstADvar(const ADvari &x)
00310 {
00311   ConstADvari *y = new ConstADvari(x.Val);
00312   new Derp(&CADcontext::One, y, &x); /* for side effect; value not used */
00313   cv = y;
00314   }
00315 
00316  void
00317 IndepADvar::AD_Const(const IndepADvar &v)
00318 {
00319   ConstADvari *ncv = new ConstADvari(v.val());
00320   ((IndepADvar*)&v)->cv = ncv;
00321   }
00322 
00323  void
00324 ConstADvari::aval_reset()
00325 {
00326   ConstADvari *x = ConstADvari::lastcad;
00327   while(x) {
00328     x->aval = 0;
00329     x = x->prevcad;
00330     }
00331   }
00332 
00333 #ifdef RAD_AUTO_AD_Const
00334 
00335 ADvari::ADvari(const IndepADvar *x, double d): Val(d), aval(0.)
00336 {
00337   opclass = Hv_const;
00338   *Last_ADvari = this;
00339   Last_ADvari = &Next;
00340   padv = (IndepADvar*)x;
00341   }
00342 
00343 ADvar1::ADvar1(const IndepADvar *x, const IndepADvar &y):
00344   ADvari(Hv_copy, y.cv->Val), d(&CADcontext::One, this, y.cv)
00345 {
00346   *ADvari::Last_ADvari = this;
00347   ADvari::Last_ADvari = &Next;
00348   padv = (IndepADvar*)x;
00349   }
00350 
00351 ADvar1::ADvar1(const IndepADvar *x, const ADvari &y):
00352   ADvari(Hv_copy, y.Val), d(&CADcontext::One, this, &y)
00353 {
00354   *ADvari::Last_ADvari = this;
00355   ADvari::Last_ADvari = &Next;
00356   padv = (IndepADvar*)x;
00357   }
00358 
00359 #else
00360 
00361  ADvar&
00362 ADvar_operatoreq(ADvar *This, const ADvari &x)
00363 { This->cv = new ADvar1(Hv_copy, x.Val, &CADcontext::One, (ADvari*)&x); return *This; }
00364 
00365  IndepADvar&
00366 ADvar_operatoreq(IndepADvar *This, const ADvari &x)
00367 { This->cv = new ADvar1(Hv_copy, x.Val, &CADcontext::One, (ADvari*)&x); return *This; }
00368 
00369 #endif /*RAD_AUTO_AD_Const*/
00370 
00371 ADvar2q::ADvar2q(double val1, double Lp, double Rp,
00372     double LR, double R2, const ADvari *Lcv, const ADvari *Rcv):
00373       ADvar2(Hv_quotLR,val1,Lcv,&pL,Rcv,&pR),
00374       pL(Lp), pR(Rp), pLR(LR), pR2(R2) { }
00375 
00376 ADvar2g::ADvar2g(double val1, double Lp, double Rp,
00377     double L2, double LR, double R2, const ADvari *Lcv, const ADvari *Rcv):
00378       ADvar2(Hv_binary,val1,Lcv,&pL,Rcv,&pR),
00379       pL(Lp), pR(Rp), pL2(L2), pLR(LR), pR2(R2) { }
00380 
00381  IndepADvar&
00382 IndepADvar::operator=(double d)
00383 {
00384 #ifdef RAD_AUTO_AD_Const
00385   if (cv)
00386     cv->padv = 0;
00387   cv = new ADvari(this,d);
00388 #else
00389   cv = new ADvari(Hv_const, d);
00390 #endif
00391   return *this;
00392   }
00393 
00394  ADvar&
00395 ADvar::operator=(double d)
00396 {
00397 #ifdef RAD_AUTO_AD_Const
00398   if (cv)
00399     cv->padv = 0;
00400   cv = new ADvari(this, d);
00401 #else
00402   cv = ConstADvari::cadc.fpval_implies_const
00403     ? new ConstADvari(d)
00404     : new ADvari(Hv_const, d);
00405 #endif
00406   return *this;
00407   }
00408 
00409  ADvari&
00410 operator-(const ADvari &T) {
00411   return *(new ADvar1(Hv_negate, -T.Val, &CADcontext::negOne, &T));
00412   }
00413 
00414  ADvari&
00415 operator+(const ADvari &L, const ADvari &R) {
00416   return *(new ADvar2(Hv_plusLR,L.Val + R.Val, &L, &CADcontext::One, &R, &CADcontext::One));
00417   }
00418 
00419  ADvar&
00420 ADvar::operator+=(const ADvari &R) {
00421   ADvari *Lcv = cv;
00422 #ifdef RAD_AUTO_AD_Const
00423   Lcv->padv = 0;
00424 #endif
00425   cv = new ADvar2(Hv_plusLR, Lcv->Val + R.Val, Lcv, &CADcontext::One, &R, &CADcontext::One);
00426   return *this;
00427   }
00428 
00429  ADvari&
00430 operator+(const ADvari &L, double R) {
00431   return *(new ADvar1(Hv_copy, L.Val + R, &CADcontext::One, &L));
00432   }
00433 
00434  ADvar&
00435 ADvar::operator+=(double R) {
00436   ADvari *tcv = cv;
00437 #ifdef RAD_AUTO_AD_Const
00438   tcv->padv = 0;
00439 #endif
00440   cv = new ADvar1(Hv_copy, tcv->Val + R, &CADcontext::One, tcv);
00441   return *this;
00442   }
00443 
00444  ADvari&
00445 operator+(double L, const ADvari &R) {
00446   return *(new ADvar1(Hv_copy, L + R.Val, &CADcontext::One, &R));
00447   }
00448 
00449  ADvari&
00450 operator-(const ADvari &L, const ADvari &R) {
00451   return *(new ADvar2(Hv_minusLR,L.Val - R.Val, &L, &CADcontext::One, &R, &CADcontext::negOne));
00452   }
00453 
00454  ADvar&
00455 ADvar::operator-=(const ADvari &R) {
00456   ADvari *Lcv = cv;
00457 #ifdef RAD_AUTO_AD_Const
00458   Lcv->padv = 0;
00459 #endif
00460   cv = new ADvar2(Hv_minusLR, Lcv->Val - R.Val, Lcv, &CADcontext::One, &R, &CADcontext::negOne);
00461   return *this;
00462   }
00463 
00464  ADvari&
00465 operator-(const ADvari &L, double R) {
00466   return *(new ADvar1(Hv_copy, L.Val - R, &CADcontext::One, &L));
00467   }
00468 
00469  ADvar&
00470 ADvar::operator-=(double R) {
00471   ADvari *tcv = cv;
00472 #ifdef RAD_AUTO_AD_Const
00473   tcv->padv = 0;
00474 #endif
00475   cv = new ADvar1(Hv_copy, tcv->Val - R, &CADcontext::One, tcv);
00476   return *this;
00477   }
00478 
00479  ADvari&
00480 operator-(double L, const ADvari &R) {
00481   return *(new ADvar1(Hv_negate, L - R.Val, &CADcontext::negOne, &R));
00482   }
00483 
00484  ADvari&
00485 operator*(const ADvari &L, const ADvari &R) {
00486   return *(new ADvar2(Hv_timesLR, L.Val * R.Val, &L, &R.Val, &R, &L.Val));
00487   }
00488 
00489  ADvar&
00490 ADvar::operator*=(const ADvari &R) {
00491   ADvari *Lcv = cv;
00492 #ifdef RAD_AUTO_AD_Const
00493   Lcv->padv = 0;
00494 #endif
00495   cv = new ADvar2(Hv_timesLR, Lcv->Val * R.Val, Lcv, &R.Val, &R, &Lcv->Val);
00496   return *this;
00497   }
00498 
00499  ADvari&
00500 operator*(const ADvari &L, double R) {
00501   return *(new ADvar1s(L.Val * R, R, &L));
00502   }
00503 
00504  ADvar&
00505 ADvar::operator*=(double R) {
00506   ADvari *Lcv = cv;
00507 #ifdef RAD_AUTO_AD_Const
00508   Lcv->padv = 0;
00509 #endif
00510   cv = new ADvar1s(Lcv->Val * R, R, Lcv);
00511   return *this;
00512   }
00513 
00514  ADvari&
00515 operator*(double L, const ADvari &R) {
00516   return *(new ADvar1s(L * R.Val, L, &R));
00517   }
00518 
00519  ADvari&
00520 operator/(const ADvari &L, const ADvari &R) {
00521   double Lv = L.Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv, qpL = q*pL;
00522   return *(new ADvar2q(q, pL, -qpL, -pL*pL, 2.*pL*qpL, &L, &R));
00523   }
00524 
00525  ADvar&
00526 ADvar::operator/=(const ADvari &R) {
00527   ADvari *Lcv = cv;
00528 #ifdef RAD_AUTO_AD_Const
00529   Lcv->padv = 0;
00530 #endif
00531   double Lv = Lcv->Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv, qpL = q*pL;
00532   cv = new ADvar2q(q, pL, -qpL, -pL*pL, 2.*pL*qpL, Lcv, &R);
00533   return *this;
00534   }
00535 
00536  ADvari&
00537 operator/(const ADvari &L, double R) {
00538   return *(new ADvar1s(L.Val / R, 1./R, &L));
00539   }
00540 
00541  ADvari&
00542 operator/(double L, const ADvari &R) {
00543   double recip = 1. / R.Val;
00544   double q = L * recip;
00545   double d1 = -q*recip;
00546   return *(new ADvar1g(q, d1, -q*d1, &R));
00547   }
00548 
00549  ADvar&
00550 ADvar::operator/=(double R) {
00551   ADvari *Lcv = cv;
00552 #ifdef RAD_AUTO_AD_Const
00553   Lcv->padv = 0;
00554 #endif
00555   cv = new ADvar1s(Lcv->Val / R, 1./R, Lcv);
00556   return *this;
00557   }
00558 
00559  ADvari&
00560 acos(const ADvari &v) {
00561   double t = v.Val;
00562   double t1 = 1. - t*t;
00563   double d1 = -1. / sqrt(t1);
00564   return *(new ADvar1g(acos(t), d1, t*d1/t1, &v));
00565   }
00566 
00567  ADvari&
00568 acosh(const ADvari &v) {
00569   double d1, t, t1, t2;
00570   t = v.Val;
00571   t1 = sqrt(t2 = t*t - 1.);
00572   d1 = 1./t1;
00573   return *(new ADvar1g(log(t + t1), d1, -t*d1/t2, &v));
00574   }
00575 
00576  ADvari&
00577 asin(const ADvari &v) {
00578   double d1, t, t1;
00579   t = v.Val;
00580   d1 = 1. / sqrt(t1 = 1. - t*t);
00581   return *(new ADvar1g(asin(t), d1, t*d1/t1, &v));
00582   }
00583 
00584  ADvari&
00585 asinh(const ADvari &v) {
00586   double d1, t, t1, t2, td;
00587   t = v.Val;
00588   t1 = sqrt(t2 = t*t + 1.);
00589   d1 = 1. / t1;
00590   td = 1.;
00591   if (t < 0.)
00592     td = -1.;
00593   return *(new ADvar1g(td*log(t*td + t1), d1, -(t/t2)*d1, &v));
00594   }
00595 
00596  ADvari&
00597 atan(const ADvari &v) {
00598   double d1, t;
00599   t = v.Val;
00600   d1 = 1./(1. + t*t);
00601   return *(new ADvar1g(atan(t), d1, -(t+t)*d1*d1, &v));
00602   }
00603 
00604  ADvari&
00605 atanh(const ADvari &v) {
00606   double t = v.Val;
00607   double d1 = 1./(1. - t*t);
00608   return *(new ADvar1g(0.5*log((1.+t)/(1.-t)), d1, (t+t)*d1*d1, &v));
00609   }
00610 
00611  ADvari&
00612 max(const ADvari &L, const ADvari &R) {
00613   const ADvari &x = L.Val >= R.Val ? L : R;
00614   return *(new ADvar1(Hv_copy, x.Val, &CADcontext::One, &x));
00615   }
00616 
00617  ADvari&
00618 max(double L, const ADvari &R) {
00619   if (L >= R.Val)
00620     return *(new ADvari(Hv_const, L));
00621   return *(new ADvar1(Hv_copy, R.Val, &CADcontext::One, &R));
00622   }
00623 
00624  ADvari&
00625 max(const ADvari &L, double R) {
00626   if (L.Val >= R)
00627     return *(new ADvar1(Hv_copy, L.Val, &CADcontext::One, &L));
00628   return *(new ADvari(Hv_const, R));
00629   }
00630 
00631  ADvari&
00632 min(const ADvari &L, const ADvari &R) {
00633   const ADvari &x = L.Val <= R.Val ? L : R;
00634   return *(new ADvar1(Hv_copy, x.Val, &CADcontext::One, &x));
00635   }
00636 
00637  ADvari&
00638 min(double L, const ADvari &R) {
00639   if (L <= R.Val)
00640     return *(new ADvari(Hv_const, L));
00641   return *(new ADvar1(Hv_copy, R.Val, &CADcontext::One, &R));
00642   }
00643 
00644  ADvari&
00645 min(const ADvari &L, double R) {
00646   if (L.Val <= R)
00647     return *(new ADvar1(Hv_copy, L.Val, &CADcontext::One, &L));
00648   return *(new ADvari(Hv_const, R));
00649   }
00650 
00651  ADvari&
00652 atan2(const ADvari &L, const ADvari &R) {
00653   double x = L.Val, y = R.Val;
00654   double x2 = x*x, y2 = y*y;
00655   double t = 1./(x2 + y2);
00656   double t2 = t*t;
00657   double R2 = 2.*t2*x*y;
00658   return *(new ADvar2g(atan2(x,y), y*t, -x*t, -R2, t2*(x2 - y2), R2, &L, &R));
00659   }
00660 
00661  ADvari&
00662 atan2(double x, const ADvari &R) {
00663   double y = R.Val;
00664   double x2 = x*x, y2 = y*y;
00665   double t = 1./(x2 + y2);
00666   return *(new ADvar1g(atan2(x,y), -x*t, 2.*t*t*x*y, &R));
00667   }
00668 
00669  ADvari&
00670 atan2(const ADvari &L, double y) {
00671   double x = L.Val;
00672   double x2 = x*x, y2 = y*y;
00673   double t = 1./(x2 + y2);
00674   return *(new ADvar1g(atan2(x,y), y*t, -2.*t*t*x*y, &L));
00675   }
00676 
00677  ADvari&
00678 cos(const ADvari &v) {
00679   double t = cos(v.Val);
00680   return *(new ADvar1g(t, -sin(v.Val), -t, &v));
00681   }
00682 
00683  ADvari&
00684 cosh(const ADvari &v) {
00685   double t = cosh(v.Val);
00686   return *(new ADvar1g(t, sinh(v.Val), t, &v));
00687   }
00688 
00689  ADvari&
00690 exp(const ADvari &v) {
00691   double t = exp(v.Val);
00692   return *(new ADvar1g(t, t, t, &v));
00693   }
00694 
00695  ADvari&
00696 log(const ADvari &v) {
00697   double x = v.Val;
00698   double d1 = 1. / x;
00699   return *(new ADvar1g(log(x), d1, -d1*d1, &v));
00700   }
00701 
00702  ADvari&
00703 log10(const ADvari &v) {
00704   static double num = 1. / log(10.);
00705   double x = v.Val, t = 1. / x;
00706   double d1 = num * t;
00707   return *(new ADvar1g(log10(x), d1, -d1*t, &v));
00708   }
00709 
00710  ADvari&
00711 pow(const ADvari &L, const ADvari &R) {
00712   double x = L.Val, y = R.Val, t = pow(x,y);
00713   double xym1 = t/x;
00714   double xlog = log(x);
00715   double dx = y*xym1;
00716   double dy = t * xlog;
00717   return *(new ADvar2g(t, dx, dy, (y-1.)*dx/x, xym1*(1. + y*xlog), dy*xlog, &L, &R));
00718   }
00719 
00720  ADvari&
00721 pow(double x, const ADvari &R) {
00722   double y = R.Val, t = pow(x,y);
00723   double xlog = log(x);
00724   double dy = t * xlog;
00725   return *(new ADvar1g(t, dy, dy*xlog, &R));
00726   }
00727 
00728  ADvari&
00729 pow(const ADvari &L, double y) {
00730   double x = L.Val, t = pow(x,y);
00731   double dx = y*t/x;
00732   return *(new ADvar1g(t, dx, (y-1.)*dx/x, &L));
00733   }
00734 
00735  ADvari&
00736 sin(const ADvari &v) {
00737   double t = sin(v.Val);
00738   return *(new ADvar1g(t, cos(v.Val), -t, &v));
00739   }
00740 
00741  ADvari&
00742 sinh(const ADvari &v) {
00743   double t = sinh(v.Val);
00744   return *(new ADvar1g(t, cosh(v.Val), t, &v));
00745   }
00746 
00747  ADvari&
00748 sqrt(const ADvari &v) {
00749   double t = sqrt(v.Val);
00750   double d1 = 0.5/t;
00751   return *(new ADvar1g(t, d1, -0.5*d1/v.Val, &v));
00752   }
00753 
00754  ADvari&
00755 tan(const ADvari &v) {
00756   double d1, rv, t;
00757   rv = tan(v.Val);
00758   t = 1. / cos(v.Val);
00759   d1 = t*t;
00760   return *(new ADvar1g(rv, d1, (rv+rv)*d1, &v));
00761   }
00762 
00763  ADvari&
00764 tanh(const ADvari &v) {
00765   double d1, rv, t;
00766   rv = tanh(v.Val);
00767   t = 1. / cosh(v.Val);
00768   d1 = t*t;
00769   return *(new ADvar1g(rv, d1, -(rv+rv)*d1, &v));
00770   }
00771 
00772  ADvari&
00773 fabs(const ADvari &v) { // "fabs" is not the best choice of name,
00774       // but this name is used at Sandia.
00775   double t, p;
00776   p = 1;
00777   if ((t = v.Val) < 0) {
00778     t = -t;
00779     p = -p;
00780     }
00781   return *(new ADvar1g(t, p, 0., &v));
00782   }
00783 
00784  ADvari&
00785 ADf1(double f, double g, double h, const ADvari &x) {
00786   return *(new ADvar1g(f, g, h, &x));
00787   }
00788 
00789  ADvari&
00790 ADf2(double f, double gx, double gy, double hxx, double hxy, double hyy,
00791     const ADvari &x, const ADvari &y) {
00792   return *(new ADvar2g(f, gx, gy, hxx, hxy, hyy, &x, &y));
00793   }
00794 
00795  ADvarn::ADvarn(double val1, int n1, const ADvar *x, const double *g, const double *h):
00796     ADvari(Hv_nary,val1), n(n1)
00797 {
00798   Derp *d1, *dlast;
00799   double *a1;
00800   int i, nh;
00801 
00802   a1 = G = (double*)ADvari::adc.Memalloc(n1*sizeof(*g));
00803   d1 = D =  (Derp*)ADvari::adc.Memalloc(n1*sizeof(Derp));
00804   dlast = Derp::LastDerp;
00805   for(i = 0; i < n1; i++, d1++) {
00806     d1->next = dlast;
00807     dlast = d1;
00808     a1[i] = g[i];
00809     d1->a = &a1[i];
00810     d1->b = this;
00811     d1->c = x[i].cv;
00812     }
00813   Derp::LastDerp = dlast;
00814   nh = (n1*(n1+1)) >> 1;
00815   a1 = H = (double*)ADvari::adc.Memalloc(nh * sizeof(*h));
00816   for(i = 0; i < nh; i++)
00817     a1[i] = h[i];
00818   }
00819 
00820  ADvari&
00821 ADfn(double f, int n, const ADvar *x, const double *g, const double *h) {
00822   return *(new ADvarn(f, n, x, g, h));
00823   }
00824 
00825  void
00826 ADcontext::Hvprod(int n, ADvar **x, double *v, double *hv)
00827 {
00828   ADvari *a, *aL, *aR, **ap, **ape;
00829   ADvari_block *b, *b0;
00830   Derp *d;
00831   double aO, adO, *g, *h, *h0, t, tL, tR;
00832   int i, j, k, m;
00833   for(i = 0; i < n; i++) {
00834     a = x[i]->cv;
00835     a->dO = v[i];
00836     a->aO = a->adO = 0.;
00837     }
00838   ADvari::adc.Aibusy->limit = ADvari::adc.Ainext;
00839   a = 0;
00840   for(b0 = 0, b = &ADvari::adc.AiFirst; b; b0 = b, b = b->next) {
00841     ap = b->pADvari;
00842     ape = b->limit;
00843     while(ap < ape) {
00844       a = *ap++;
00845       a->aO = a->adO = 0.;
00846       switch(a->opclass) {
00847        case Hv_copy:
00848         a->dO = ((ADvar1*)a)->d.c->dO;
00849         break;
00850        case Hv_binary:
00851         a->dO =   ((ADvar2g*)a)->pL * ((ADvar2g*)a)->dL.c->dO
00852           + ((ADvar2g*)a)->pR * ((ADvar2g*)a)->dR.c->dO;
00853         break;
00854        case Hv_unary:
00855         a->dO = ((ADvar1g*)a)->pL * ((ADvar1g*)a)->d.c->dO;
00856         break;
00857        case Hv_negate:
00858         a->dO = -((ADvar1*)a)->d.c->dO;
00859         break;
00860        case Hv_plusLR:
00861         a->dO = ((ADvar2*)a)->dL.c->dO + ((ADvar2*)a)->dR.c->dO;
00862         break;
00863        case Hv_minusLR:
00864         a->dO = ((ADvar2*)a)->dL.c->dO - ((ADvar2*)a)->dR.c->dO;
00865         break;
00866        case Hv_timesL:
00867         a->dO = ((ADvar1s*)a)->pL * ((ADvar1s*)a)->d.c->dO;
00868         break;
00869        case Hv_timesLR:
00870         a->dO =   ((ADvar2*)a)->dR.c->Val * ((ADvar2*)a)->dL.c->dO
00871           + ((ADvar2*)a)->dL.c->Val * ((ADvar2*)a)->dR.c->dO;
00872         break;
00873        case Hv_quotLR:
00874         a->dO =   ((ADvar2q*)a)->pL * ((ADvar2q*)a)->dL.c->dO
00875           + ((ADvar2q*)a)->pR * ((ADvar2q*)a)->dR.c->dO;
00876        case Hv_const: // This case does not arise; the intent here
00877           // is to eliminate a red-herring compiler warning.
00878         break;
00879        case Hv_nary:
00880         d = ((ADvarn*)a)->D;
00881         m = ((ADvarn*)a)->n;
00882         g = ((ADvarn*)a)->G;
00883         t = 0.;
00884         for(i = 0; i < m; i++)
00885           t += g[i] * d[i].c->dO;
00886         a->dO = t;
00887        }
00888       }
00889     }
00890   if (a)
00891     a->adO = 1.;
00892   for(b = b0; b; b = b->prev) {
00893     ape = b->pADvari;
00894     ap = b->limit;
00895     while(ap > ape) {
00896       a = *--ap;
00897       aO = a->aO;
00898       adO = a->adO;
00899       switch(a->opclass) {
00900        case Hv_copy:
00901         aL = ((ADvar1*)a)->d.c;
00902         aL->aO += aO;
00903         aL->adO += adO;
00904         break;
00905        case Hv_binary:
00906         aL = ((ADvar2g*)a)->dL.c;
00907         aR = ((ADvar2g*)a)->dR.c;
00908         tL = adO*aL->dO;
00909         tR = adO*aR->dO;
00910         aL->aO += aO*((ADvar2g*)a)->pL
00911           + tL*((ADvar2g*)a)->pL2
00912           + tR*((ADvar2g*)a)->pLR;
00913         aR->aO += aO*((ADvar2g*)a)->pR
00914           + tL*((ADvar2g*)a)->pLR
00915           + tR*((ADvar2g*)a)->pR2;
00916         aL->adO += adO * ((ADvar2g*)a)->pL;
00917         aR->adO += adO * ((ADvar2g*)a)->pR;
00918         break;
00919        case Hv_unary:
00920         aL = ((ADvar1g*)a)->d.c;
00921         aL->aO += aO * ((ADvar1g*)a)->pL
00922           + adO * aL->dO * ((ADvar1g*)a)->pL2;
00923         aL->adO += adO * ((ADvar1g*)a)->pL;
00924         break;
00925        case Hv_negate:
00926         aL = ((ADvar1*)a)->d.c;
00927         aL->aO -= aO;
00928         aL->adO -= adO;
00929         break;
00930        case Hv_plusLR:
00931         aL = ((ADvar2*)a)->dL.c;
00932         aR = ((ADvar2*)a)->dR.c;
00933         aL->aO += aO;
00934         aL->adO += adO;
00935         aR->aO += aO;
00936         aR->adO += adO;
00937         break;
00938        case Hv_minusLR:
00939         aL = ((ADvar2*)a)->dL.c;
00940         aR = ((ADvar2*)a)->dR.c;
00941         aL->aO += aO;
00942         aL->adO += adO;
00943         aR->aO -= aO;
00944         aR->adO -= adO;
00945         break;
00946        case Hv_timesL:
00947         aL = ((ADvar1s*)a)->d.c;
00948         aL->aO += aO * (tL = ((ADvar1s*)a)->pL);
00949         aL->adO += adO * tL;
00950         break;
00951        case Hv_timesLR:
00952         aL = ((ADvar2*)a)->dL.c;
00953         aR = ((ADvar2*)a)->dR.c;
00954         aL->aO += aO * (tL = aR->Val) + adO*aR->dO;
00955         aR->aO += aO * (tR = aL->Val) + adO*aL->dO;
00956         aL->adO += adO * tL;
00957         aR->adO += adO * tR;
00958         break;
00959        case Hv_quotLR:
00960         aL = ((ADvar2q*)a)->dL.c;
00961         aR = ((ADvar2q*)a)->dR.c;
00962         tL = adO*aL->dO;
00963         tR = adO*aR->dO;
00964         aL->aO += aO*((ADvar2q*)a)->pL
00965           + tR*((ADvar2q*)a)->pLR;
00966         aR->aO += aO*((ADvar2q*)a)->pR
00967           + tL*((ADvar2q*)a)->pLR
00968           + tR*((ADvar2q*)a)->pR2;
00969         aL->adO += adO * ((ADvar2q*)a)->pL;
00970         aR->adO += adO * ((ADvar2q*)a)->pR;
00971        case Hv_const: // This case does not arise; the intent here
00972           // is to eliminate a red-herring compiler warning.
00973         break;
00974        case Hv_nary:
00975         d  = ((ADvarn*)a)->D;
00976         m  = ((ADvarn*)a)->n;
00977         g  = ((ADvarn*)a)->G;
00978         h0 = ((ADvarn*)a)->H;
00979         for(i = 0; i < m; i++) {
00980           aL = d[i].c;
00981           aL->adO += adO * (t = g[i]);
00982           aL->aO += t*aO;
00983           t = adO * aL->dO;
00984           for(h = h0, j = 0; j <= i; j++)
00985             d[j].c->aO += t * *h++;
00986           h0 = h--;
00987           for(k = j; j < m; j++)
00988             d[j].c->aO += t * *(h += k++);
00989           }
00990        }
00991       }
00992     }
00993   for(i = 0; i < n; i++) {
00994     a = x[i]->cv;
00995     a->dO = 0.;
00996     hv[i] = a->aO;
00997     }
00998   }
00999 
01000 #ifndef SACADO_NO_NAMESPACE
01001 } // namespace Rad2d
01002 } // namespace Sacado
01003 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines