Sacado Package Browser (Single Doxygen Collection) Version of the Day
Sacado_radops.cpp
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 // Support routines for the 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 #include "Sacado_rad.hpp"
00035 
00036 namespace Sacado {
00037 namespace Radnt {
00038 
00039 #ifdef RAD_AUTO_AD_Const
00040 ADvari *ADvari::First_ADvari, **ADvari::Last_ADvari = &ADvari::First_ADvari;
00041 #undef RAD_DEBUG_BLOCKKEEP
00042 #else
00043 #ifdef RAD_DEBUG_BLOCKKEEP
00044 #if !(RAD_DEBUG_BLOCKKEEP > 0)
00045 #undef RAD_DEBUG_BLOCKKEEP
00046 #else
00047 extern "C" void _uninit_f2c(void *x, int type, long len);
00048 #define TYDREAL 5
00049 static ADmemblock *rad_Oldcurmb;
00050 static int rad_busy_blocks;
00051 #endif /*RAD_DEBUG_BLOCKKEEP > 0*/
00052 #endif /*RAD_DEBUG_BLOCKKEEP*/
00053 #endif /*RAD_AUTO_AD_Const*/
00054 
00055 Derp *Derp::LastDerp = 0;
00056 ADcontext ADvari::adc;
00057 CADcontext ConstADvari::cadc;
00058 ConstADvari *ConstADvari::lastcad;
00059 static int rad_need_reinit;
00060 #ifdef RAD_DEBUG_BLOCKKEEP
00061 static size_t rad_mleft_save;
00062 #endif
00063 
00064 const double CADcontext::One = 1., CADcontext::negOne = -1.;
00065 
00066 ADcontext::ADcontext()
00067 {
00068   First.next = 0;
00069   Busy = &First;
00070   Free = 0;
00071   Mbase = (char*)First.memblk;
00072   Mleft = sizeof(First.memblk);
00073   }
00074 
00075  void*
00076 ADcontext::new_ADmemblock(size_t len)
00077 {
00078   ADmemblock *mb, *mb0, *mb1, *mbf, *x;
00079 #ifdef RAD_AUTO_AD_Const
00080   ADvari *a, *anext;
00081   IndepADvar *v;
00082 #endif /* RAD_AUTO_AD_Const */
00083 
00084   if (rad_need_reinit && this == &ADvari::adc) {
00085     rad_need_reinit = 0;
00086     Derp::LastDerp = 0;
00087 #ifdef RAD_DEBUG_BLOCKKEEP
00088     Mleft = rad_mleft_save;
00089     if (Mleft < sizeof(First.memblk))
00090       _uninit_f2c(Mbase + Mleft, TYDREAL,
00091        (sizeof(First.memblk) - Mleft)/sizeof(double));
00092     if ((mb = Busy->next)) {
00093       if (!(mb0 = rad_Oldcurmb))
00094         mb0 = &First;
00095       for(;; mb = mb->next) {
00096         _uninit_f2c(mb->memblk, TYDREAL,
00097           sizeof(First.memblk)/sizeof(double));
00098         if (mb == mb0)
00099           break;
00100         }
00101       }
00102     rad_Oldcurmb = Busy;
00103     if (rad_busy_blocks >= RAD_DEBUG_BLOCKKEEP) {
00104       rad_busy_blocks = 0;
00105       rad_Oldcurmb = 0;
00106       mb0 = &First;
00107       mbf =  Free;
00108       for(mb = Busy; mb != mb0; mb = mb1) {
00109         mb1 = mb->next;
00110         mb->next = mbf;
00111         mbf = mb;
00112         }
00113       Free = mbf;
00114       Busy = mb;
00115       Mbase = (char*)First.memblk;
00116       Mleft = sizeof(First.memblk);
00117       }
00118 
00119 #else /* !RAD_DEBUG_BLOCKKEEP */
00120 
00121     mb0 = &ADvari::adc.First;
00122     mbf =  ADvari::adc.Free;
00123     for(mb = ADvari::adc.Busy; mb != mb0; mb = mb1) {
00124       mb1 = mb->next;
00125       mb->next = mbf;
00126       mbf = mb;
00127       }
00128     ADvari::adc.Free = mbf;
00129     ADvari::adc.Busy = mb;
00130     ADvari::adc.Mbase = (char*)ADvari::adc.First.memblk;
00131     ADvari::adc.Mleft = sizeof(ADvari::adc.First.memblk);
00132 #ifdef RAD_AUTO_AD_Const
00133     *ADvari::Last_ADvari = 0;
00134     ADvari::Last_ADvari = &ADvari::First_ADvari;
00135     if ((anext = ADvari::First_ADvari)) {
00136       while((a = anext)) {
00137         anext = a->Next;
00138         if ((v = a->padv))
00139           v->cv = new ADvari(v, a->Val);
00140         }
00141       }
00142 #endif /*RAD_AUTO_AD_Const*/
00143 #endif /* RAD_DEBUG_BLOCKKEEP */
00144     if (Mleft >= len)
00145       return Mbase + (Mleft -= len);
00146     }
00147 
00148   if ((x = Free))
00149     Free = x->next;
00150   else
00151     x = new ADmemblock;
00152 #ifdef RAD_DEBUG_BLOCKKEEP
00153   rad_busy_blocks++;
00154 #endif
00155   x->next = Busy;
00156   Busy = x;
00157   return (Mbase = (char*)x->memblk) +
00158     (Mleft = sizeof(First.memblk) - len);
00159   }
00160 
00161  void
00162 ADcontext::Gradcomp(int wantgrad)
00163 {
00164   Derp *d;
00165 
00166   if (rad_need_reinit && wantgrad) {
00167     for(d = Derp::LastDerp; d; d = d->next)
00168       d->c->aval = 0;
00169     }
00170   else {
00171     rad_need_reinit = 1;
00172 #ifdef RAD_DEBUG_BLOCKKEEP
00173     rad_mleft_save = ADvari::adc.Mleft;
00174 #endif
00175     ADvari::adc.Mleft = 0;
00176     }
00177 
00178   if ((d = Derp::LastDerp) && wantgrad) {
00179     d->b->aval = 1;
00180     do d->c->aval += *d->a * d->b->aval;
00181     while((d = d->next));
00182     }
00183   }
00184 
00185  void
00186 ADcontext::Weighted_Gradcomp(int n, ADvar **v, double *w)
00187 {
00188   Derp *d;
00189   int i;
00190 
00191   if (rad_need_reinit) {
00192     for(d = Derp::LastDerp; d; d = d->next)
00193       d->c->aval = 0;
00194     }
00195   else {
00196     rad_need_reinit = 1;
00197 #ifdef RAD_DEBUG_BLOCKKEEP
00198     rad_mleft_save = ADvari::adc.Mleft;
00199 #endif
00200     ADvari::adc.Mleft = 0;
00201     }
00202   if ((d = Derp::LastDerp)) {
00203     for(i = 0; i < n; i++)
00204       v[i]->cv->aval = w[i];
00205     do d->c->aval += *d->a * d->b->aval;
00206     while((d = d->next));
00207     }
00208   }
00209 
00210 #ifdef RAD_AUTO_AD_Const
00211 
00212 IndepADvar::IndepADvar(double d)
00213 {
00214   cv = new ADvari(this,d);
00215   }
00216 
00217 IndepADvar::IndepADvar(int d)
00218 {
00219   cv = new ADvari(this,d);
00220   }
00221 
00222 IndepADvar::IndepADvar(long d)
00223 {
00224   cv = new ADvari(this,d);
00225   }
00226 
00227 #else 
00229 IndepADvar::IndepADvar(double d)
00230 {
00231   ADvari *x = new ADvari(d);
00232   cv = x;
00233   }
00234 
00235 IndepADvar::IndepADvar(int d)
00236 {
00237   ADvari *x = new ADvari((double)d);
00238   cv = x;
00239   }
00240 
00241 IndepADvar::IndepADvar(long d)
00242 {
00243   ADvari *x = new ADvari((double)d);
00244   cv = x;
00245   }
00246 
00247 #endif /*RAD_AUTO_AD_Const*/
00248 
00249  void
00250 ADvar::ADvar_ctr(double d)
00251 {
00252 #ifdef RAD_AUTO_AD_Const
00253   ADvari *x = new ADvari(this,d);
00254 #else
00255   ADvari *x = ConstADvari::cadc.fpval_implies_const
00256     ? new ConstADvari(d)
00257     : new ADvari(d);
00258 #endif
00259   cv = x;
00260   }
00261 
00262 ConstADvar::ConstADvar()
00263 {
00264   ConstADvari *x = new ConstADvari(0.);
00265   cv = x;
00266   }
00267 
00268  void
00269 ConstADvar::ConstADvar_ctr(double d)
00270 {
00271   ConstADvari *x = new ConstADvari(d);
00272   cv = x;
00273   }
00274 ConstADvar::ConstADvar(const ADvari &x)
00275 {
00276   ConstADvari *y = new ConstADvari(x.Val);
00277   new Derp(&CADcontext::One, y, &x); /*for side effect; value not used */
00278   cv = y;
00279   }
00280 
00281  void
00282 IndepADvar::AD_Const(const IndepADvar &v)
00283 {
00284   ConstADvari *ncv = new ConstADvari(v.val());
00285 #ifdef RAD_AUTO_AD_Const
00286   v.cv->padv = 0;
00287 #endif
00288   ((IndepADvar*)&v)->cv = ncv;
00289   }
00290 
00291  void
00292 ConstADvari::aval_reset()
00293 {
00294   ConstADvari *x = ConstADvari::lastcad;
00295   while(x) {
00296     x->aval = 0;
00297     x = x->prevcad;
00298     }
00299   }
00300 
00301 #ifdef RAD_AUTO_AD_Const
00302 
00303 ADvari::ADvari(const IndepADvar *x, double d): Val(d), aval(0.)
00304 {
00305   *Last_ADvari = this;
00306   Last_ADvari = &Next;
00307   padv = (IndepADvar*)x;
00308   }
00309 
00310 ADvar1::ADvar1(const IndepADvar *x, const IndepADvar &y):
00311   ADvari(y.cv->Val), d(&CADcontext::One, this, y.cv)
00312 {
00313   *ADvari::Last_ADvari = this;
00314   ADvari::Last_ADvari = &Next;
00315   padv = (IndepADvar*)x;
00316   }
00317 
00318 ADvar1::ADvar1(const IndepADvar *x, const ADvari &y):
00319   ADvari(y.Val), d(&CADcontext::One, this, &y)
00320 {
00321   *ADvari::Last_ADvari = this;
00322   ADvari::Last_ADvari = &Next;
00323   padv = (IndepADvar*)x;
00324   }
00325 
00326 #else
00327 
00328  ADvar&
00329 ADvar_operatoreq(ADvar *This, const ADvari &x)
00330 { This->cv = new ADvar1(x.Val, &CADcontext::One, (ADvari*)&x); return *This; }
00331 
00332  IndepADvar&
00333 ADvar_operatoreq(IndepADvar *This, const ADvari &x)
00334 { This->cv = new ADvar1(x.Val, &CADcontext::One, (ADvari*)&x); return *This; }
00335 
00336 #endif /*RAD_AUTO_AD_Const*/
00337 
00338  IndepADvar&
00339 IndepADvar::operator=(double d)
00340 {
00341 #ifdef RAD_AUTO_AD_Const
00342   if (cv)
00343     cv->padv = 0;
00344   cv = new ADvari(this,d);
00345 #else
00346   cv = new ADvari(d);
00347 #endif
00348   return *this;
00349   }
00350 
00351  ADvar&
00352 ADvar::operator=(double d)
00353 {
00354 #ifdef RAD_AUTO_AD_Const
00355   if (cv)
00356     cv->padv = 0;
00357   cv = new ADvari(this,d);
00358 #else
00359   cv = ConstADvari::cadc.fpval_implies_const
00360     ? new ConstADvari(d)
00361     : new ADvari(d);
00362 #endif
00363   return *this;
00364   }
00365 
00366  ADvari&
00367 operator-(const ADvari &T) {
00368   return *(new ADvar1(-T.Val, &CADcontext::negOne, &T));
00369   }
00370 
00371  ADvari&
00372 operator+(const ADvari &L, const ADvari &R) {
00373   return *(new ADvar2(L.Val + R.Val, &L, &CADcontext::One, &R, &CADcontext::One));
00374   }
00375 
00376  ADvar&
00377 ADvar::operator+=(const ADvari &R) {
00378   ADvari *Lcv = cv;
00379 #ifdef RAD_AUTO_AD_Const
00380   Lcv->padv = 0;
00381 #endif
00382   cv = new ADvar2(Lcv->Val + R.Val, Lcv, &CADcontext::One, &R, &CADcontext::One);
00383   return *this;
00384   }
00385 
00386  ADvari&
00387 operator+(const ADvari &L, double R) {
00388   return *(new ADvar1(L.Val + R, &CADcontext::One, &L));
00389   }
00390 
00391  ADvar&
00392 ADvar::operator+=(double R) {
00393   ADvari *tcv = cv;
00394 #ifdef RAD_AUTO_AD_Const
00395   tcv->padv = 0;
00396 #endif
00397   cv = new ADvar1(tcv->Val + R, &CADcontext::One, tcv);
00398   return *this;
00399   }
00400 
00401  ADvari&
00402 operator+(double L, const ADvari &R) {
00403   return *(new ADvar1(L + R.Val, &CADcontext::One, &R));
00404   }
00405 
00406  ADvari&
00407 operator-(const ADvari &L, const ADvari &R) {
00408   return *(new ADvar2(L.Val - R.Val, &L, &CADcontext::One, &R, &CADcontext::negOne));
00409   }
00410 
00411  ADvar&
00412 ADvar::operator-=(const ADvari &R) {
00413   ADvari *Lcv = cv;
00414 #ifdef RAD_AUTO_AD_Const
00415   Lcv->padv = 0;
00416 #endif
00417   cv = new ADvar2(Lcv->Val - R.Val, Lcv, &CADcontext::One, &R, &CADcontext::negOne);
00418   return *this;
00419   }
00420 
00421  ADvari&
00422 operator-(const ADvari &L, double R) {
00423   return *(new ADvar1(L.Val - R, &CADcontext::One, &L));
00424   }
00425 
00426  ADvar&
00427 ADvar::operator-=(double R) {
00428   ADvari *tcv = cv;
00429 #ifdef RAD_AUTO_AD_Const
00430   tcv->padv = 0;
00431 #endif
00432   cv = new ADvar1(tcv->Val - R, &CADcontext::One, tcv);
00433   return *this;
00434   }
00435 
00436  ADvari&
00437 operator-(double L, const ADvari &R) {
00438   return *(new ADvar1(L - R.Val, &CADcontext::negOne, &R));
00439   }
00440 
00441  ADvari&
00442 operator*(const ADvari &L, const ADvari &R) {
00443   return *(new ADvar2(L.Val * R.Val, &L, &R.Val, &R, &L.Val));
00444   }
00445 
00446  ADvar&
00447 ADvar::operator*=(const ADvari &R) {
00448   ADvari *Lcv = cv;
00449 #ifdef RAD_AUTO_AD_Const
00450   Lcv->padv = 0;
00451 #endif
00452   cv = new ADvar2(Lcv->Val * R.Val, Lcv, &R.Val, &R, &Lcv->Val);
00453   return *this;
00454   }
00455 
00456  ADvari&
00457 operator*(const ADvari &L, double R) {
00458   return *(new ADvar1s(L.Val * R, R, &L));
00459   }
00460 
00461  ADvar&
00462 ADvar::operator*=(double R) {
00463   ADvari *Lcv = cv;
00464 #ifdef RAD_AUTO_AD_Const
00465   Lcv->padv = 0;
00466 #endif
00467   cv = new ADvar1s(Lcv->Val * R, R, Lcv);
00468   return *this;
00469   }
00470 
00471  ADvari&
00472 operator*(double L, const ADvari &R) {
00473   return *(new ADvar1s(L * R.Val, L, &R));
00474   }
00475 
00476  ADvari&
00477 operator/(const ADvari &L, const ADvari &R) {
00478   double Lv = L.Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv;
00479   return *(new ADvar2q(q, pL, -q*pL, &L, &R));
00480   }
00481 
00482  ADvar&
00483 ADvar::operator/=(const ADvari &R) {
00484   ADvari *Lcv = cv;
00485 #ifdef RAD_AUTO_AD_Const
00486   Lcv->padv = 0;
00487 #endif
00488   double Lv = Lcv->Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv;
00489   cv = new ADvar2q(q, pL, -q*pL, Lcv, &R);
00490   return *this;
00491   }
00492 
00493  ADvari&
00494 operator/(const ADvari &L, double R) {
00495   return *(new ADvar1s(L.Val / R, 1./R, &L));
00496   }
00497 
00498  ADvari&
00499 operator/(double L, const ADvari &R) {
00500   double recip = 1. / R.Val;
00501   double q = L * recip;
00502   return *(new ADvar1s(q, -q*recip, &R));
00503   }
00504 
00505  ADvar&
00506 ADvar::operator/=(double R) {
00507   ADvari *Lcv = cv;
00508 #ifdef RAD_AUTO_AD_Const
00509   Lcv->padv = 0;
00510 #endif
00511   cv = new ADvar1s(Lcv->Val / R, 1./R, Lcv);
00512   return *this;
00513   }
00514 
00515  ADvari&
00516 acos(const ADvari &v) {
00517   double t = v.Val;
00518   return *(new ADvar1s(acos(t), -1./sqrt(1. - t*t), &v));
00519   }
00520 
00521  ADvari&
00522 acosh(const ADvari &v) {
00523   double t = v.Val, t1 = sqrt(t*t - 1.);
00524   return *(new ADvar1s(log(t + t1), 1./t1, &v));
00525   }
00526 
00527  ADvari&
00528 asin(const ADvari &v) {
00529   double t = v.Val;
00530   return *(new ADvar1s(asin(t), 1./sqrt(1. - t*t), &v));
00531   }
00532 
00533  ADvari&
00534 asinh(const ADvari &v) {
00535   double t = v.Val, td = 1., t1 = sqrt(t*t + 1.);
00536   if (t < 0.) {
00537     t = -t;
00538     td = -1.;
00539     }
00540   return *(new ADvar1s(td*log(t + t1), 1./t1, &v));
00541   }
00542 
00543  ADvari&
00544 atan(const ADvari &v) {
00545   double t = v.Val;
00546   return *(new ADvar1s(atan(t), 1./(1. + t*t), &v));
00547   }
00548 
00549  ADvari&
00550 atanh(const ADvari &v) {
00551   double t = v.Val;
00552   return *(new ADvar1s(0.5*log((1.+t)/(1.-t)), 1./(1. - t*t), &v));
00553   }
00554 
00555  ADvari&
00556 max(const ADvari &L, const ADvari &R) {
00557   const ADvari &x = L.Val >= R.Val ? L : R;
00558   return *(new ADvar1(x.Val, &CADcontext::One, &x));
00559   }
00560 
00561  ADvari&
00562 max(double L, const ADvari &R) {
00563   if (L >= R.Val)
00564     return *(new ADvari(L));
00565   return *(new ADvar1(R.Val, &CADcontext::One, &R));
00566   }
00567 
00568  ADvari&
00569 max(const ADvari &L, double R) {
00570   if (L.Val >= R)
00571     return *(new ADvar1(L.Val, &CADcontext::One, &L));
00572   return *(new ADvari(R));
00573   }
00574 
00575  ADvari&
00576 min(const ADvari &L, const ADvari &R) {
00577   const ADvari &x = L.Val <= R.Val ? L : R;
00578   return *(new ADvar1(x.Val, &CADcontext::One, &x));
00579   }
00580 
00581  ADvari&
00582 min(double L, const ADvari &R) {
00583   if (L <= R.Val)
00584     return *(new ADvari(L));
00585   return *(new ADvar1(R.Val, &CADcontext::One, &R));
00586   }
00587 
00588  ADvari&
00589 min(const ADvari &L, double R) {
00590   if (L.Val <= R)
00591     return *(new ADvar1(L.Val, &CADcontext::One, &L));
00592   return *(new ADvari(R));
00593   }
00594 
00595  ADvari&
00596 atan2(const ADvari &L, const ADvari &R) {
00597   double x = L.Val, y = R.Val, t = x*x + y*y;
00598   return *(new ADvar2q(atan2(x,y), y/t, -x/t, &L, &R));
00599   }
00600 
00601  ADvari&
00602 atan2(double x, const ADvari &R) {
00603   double y = R.Val, t = x*x + y*y;
00604   return *(new ADvar1s(atan2(x,y), -x/t, &R));
00605   }
00606 
00607  ADvari&
00608 atan2(const ADvari &L, double y) {
00609   double x = L.Val, t = x*x + y*y;
00610   return *(new ADvar1s(atan2(x,y), y/t, &L));
00611   }
00612 
00613  ADvari&
00614 cos(const ADvari &v) {
00615   return *(new ADvar1s(cos(v.Val), -sin(v.Val), &v));
00616   }
00617 
00618  ADvari&
00619 cosh(const ADvari &v) {
00620   return *(new ADvar1s(cosh(v.Val), sinh(v.Val), &v));
00621   }
00622 
00623  ADvari&
00624 exp(const ADvari &v) {
00625   ADvar1* rcv = new ADvar1(exp(v.Val), &v);
00626   rcv->d.a = &rcv->Val;
00627   rcv->d.b = rcv;
00628   return *rcv;
00629   }
00630 
00631  ADvari&
00632 log(const ADvari &v) {
00633   double x = v.Val;
00634   return *(new ADvar1s(log(x), 1. / x, &v));
00635   }
00636 
00637  ADvari&
00638 log10(const ADvari &v) {
00639   static double num = 1. / log(10.);
00640   double x = v.Val;
00641   return *(new ADvar1s(log10(x), num / x, &v));
00642   }
00643 
00644  ADvari&
00645 pow(const ADvari &L, const ADvari &R) {
00646   double x = L.Val, y = R.Val, t = pow(x,y);
00647   return *(new ADvar2q(t, y*t/x, t*log(x), &L, &R));
00648   }
00649 
00650  ADvari&
00651 pow(double x, const ADvari &R) {
00652   double t = pow(x,R.Val);
00653   return *(new ADvar1s(t, t*log(x), &R));
00654   }
00655 
00656  ADvari&
00657 pow(const ADvari &L, double y) {
00658   double x = L.Val, t = pow(x,y);
00659   return *(new ADvar1s(t, y*t/x, &L));
00660   }
00661 
00662  ADvari&
00663 sin(const ADvari &v) {
00664   return *(new ADvar1s(sin(v.Val), cos(v.Val), &v));
00665   }
00666 
00667  ADvari&
00668 sinh(const ADvari &v) {
00669   return *(new ADvar1s(sinh(v.Val), cosh(v.Val), &v));
00670   }
00671 
00672  ADvari&
00673 sqrt(const ADvari &v) {
00674   double t = sqrt(v.Val);
00675   return *(new ADvar1s(t, 0.5/t, &v));
00676   }
00677 
00678  ADvari&
00679 tan(const ADvari &v) {
00680   double t = cos(v.Val);
00681   return *(new ADvar1s(tan(v.Val), 1./(t*t), &v));
00682   }
00683 
00684  ADvari&
00685 tanh(const ADvari &v) {
00686   double t = 1. / cosh(v.Val);
00687   return *(new ADvar1s(tanh(v.Val), t*t, &v));
00688   }
00689 
00690  ADvari&
00691 fabs(const ADvari &v) { // "fabs" is not the best choice of name,
00692       // but this name is used at Sandia.
00693   double t, p;
00694   p = 1;
00695   if ((t = v.Val) < 0) {
00696     t = -t;
00697     p = -p;
00698     }
00699   return *(new ADvar1s(t, p, &v));
00700   }
00701 
00702  ADvari&
00703 ADf1(double f, double g, const ADvari &x) {
00704   return *(new ADvar1s(f, g, &x));
00705   }
00706 
00707  ADvari&
00708 ADf2(double f, double gx, double gy, const ADvari &x, const ADvari &y) {
00709   return *(new ADvar2q(f, gx, gy, &x, &y));
00710   }
00711 
00712 ADvarn::ADvarn(double val1, int n1, const ADvar *x, const double *g): ADvari(val1), n(n1)
00713 {
00714   Derp *d1, *dlast;
00715   double *a1;
00716   int i;
00717 
00718   a1 = a = (double*)ADvari::adc.Memalloc(n*sizeof(*a));
00719   d1 = Da =  (Derp*)ADvari::adc.Memalloc(n*sizeof(Derp));
00720   dlast = Derp::LastDerp;
00721   for(i = 0; i < n1; i++, d1++) {
00722     d1->next = dlast;
00723     dlast = d1;
00724     a1[i] = g[i];
00725     d1->a = &a1[i];
00726     d1->b = this;
00727     d1->c = x[i].cv;
00728     }
00729   Derp::LastDerp = dlast;
00730   }
00731 
00732  ADvari&
00733 ADfn(double f, int n, const ADvar *x, const double *g) {
00734   return *(new ADvarn(f, n, x, g));
00735   }
00736 
00737 } // namespace Radnt
00738 } // namespace Sacado
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines