42#if defined(RAD_DEBUG_BLOCKKEEP) && !defined(HAVE_SACADO_UNINIT)
43#undef RAD_DEBUG_BLOCKKEEP
54#undef RAD_ALLOW_WANTDERIV
55#ifndef RAD_DISALLOW_WANTDERIV
56#define RAD_ALLOW_WANTDERIV
62#define RAD_REINIT_0(x)
63#define RAD_REINIT_1(x)
64#define RAD_REINIT_2(x)
69#define RAD_REINIT_1(x) x
70#ifdef RAD_ALLOW_WANTDERIV
72#define IndepADvar IndepADvar_1n
75#define IndepADvar IndepADvar_1
79#define RAD_REINIT_2(x) x
81#define RAD_cvchk(x) if (x.gen != x.IndepADvar_root.gen) { x.cv = new ADvari<Double>(x.Val);\
82 x.gen = x.IndepADvar_root.gen; }
83#ifdef RAD_ALLOW_WANTDERIV
84#define IndepADvar IndepADvar_2n
87#define IndepADvar IndepADvar_2
91Botch! Expected
"RAD_REINIT" to be 0, 1, or 2.
94#undef RAD_ALLOW_WANTDERIV
96#define RAD_REINIT_0(x) x
99#ifdef RAD_ALLOW_WANTDERIV
100#define Allow_noderiv(x) x
102#define Allow_noderiv(x)
107#undef RAD_AUTO_AD_Const
108#undef RAD_DEBUG_BLOCKKEEP
115#ifndef RAD_AUTO_AD_Const
116#define RAD_AUTO_AD_Const
130#ifndef RAD_AUTO_AD_Const
131#ifdef RAD_DEBUG_BLOCKKEEP
139#ifdef RAD_AUTO_AD_Const
140#undef RAD_DEBUG_BLOCKKEEP
142#ifdef RAD_DEBUG_BLOCKKEEP
143#if !(RAD_DEBUG_BLOCKKEEP > 0)
144#undef RAD_DEBUG_BLOCKKEEP
146extern "C" void _uninit_f2c(
void *
x,
int type,
long len);
152struct UninitType<float> {
153 static const int utype = 4;
157struct UninitType<double> {
158 static const int utype = 5;
162struct UninitType<
std::complex<T> > {
163 static const int utype = UninitType<T>::utype + 2;
172 template<
typename T>
class
205#define Dtype typename DoubleAvoid<Double>::dtype
206#define Ltype typename DoubleAvoid<Double>::ltype
207#define Itype typename DoubleAvoid<Double>::itype
208#define Ttype typename DoubleAvoid<Double>::ttype
210 template<
typename Double>
class IndepADvar;
211 template<
typename Double>
class ConstADvar;
212 template<
typename Double>
class ConstADvari;
213 template<
typename Double>
class ADvar;
214 template<
typename Double>
class ADvari;
215 template<
typename Double>
class ADvar1;
216 template<
typename Double>
class ADvar1s;
217 template<
typename Double>
class ADvar2;
218 template<
typename Double>
class ADvar2q;
219 template<
typename Double>
class ADvarn;
220 template<
typename Double>
class Derp;
222 template<
typename Double>
struct
227 char memblk[1000*
sizeof(double)];
230 template<
typename Double>
class
243 size_t DMleft, nderps;
245#ifdef RAD_DEBUG_BLOCKKEEP
270 template<
typename Double>
class
273 bool fpval_implies_const;
275 friend class ADvar<Double>;
280 template<
typename Double>
class
283 friend class ADvarn<Double>;
287 inline void *
operator new(size_t,
Derp *
x) {
return x; }
288 inline void operator delete(
void*,
Derp *) {}
297 Derp(
const ADVari *);
298 Derp(
const Double *,
const ADVari *);
299 Derp(
const Double *,
const ADVari *,
const ADVari *);
300 inline void *
operator new(
size_t len) {
return (
Derp*)ADVari::adc.Memalloc(len); }
305 template<
typename Double> Derp<Double>*
306ADcontext<Double>::new_Derp(
const Double *
a,
const ADvari<Double> *b,
const ADvari<Double> *
c)
311 if (this->DMleft == 0) {
318 this->DMleft = nderps;
320 rv = &((DErp*)DBusy->
memblk)[--this->DMleft];
328#define Ai const Base< ADvari<Double> >&
329#define AI const Base< IndepADvar<Double> >&
330#define T template<typename Double>
357#define F ADvari<Double>&
392#ifdef HAVE_SACADO_CXX11
407 template<
typename Double>ADvari<Double>&
ADf1(Double f, Double g,
const IndepADvar<Double> &
x);
408 template<
typename Double>ADvari<Double>&
ADf2(Double f, Double gx, Double gy,
409 const IndepADvar<Double> &
x,
const IndepADvar<Double> &
y);
410 template<
typename Double>ADvari<Double>&
ADfn(Double f,
int n,
411 const IndepADvar<Double> *
x,
const Double *g);
413 template<
typename Double> IndepADvar<Double>&
ADvar_operatoreq(IndepADvar<Double>*,
414 const ADvari<Double>&);
415 template<
typename Double> ADvar<Double>&
ADvar_operatoreq(ADvar<Double>*,
const ADvari<Double>&);
416 template<
typename Double>
void AD_Const(
const IndepADvar<Double>&);
417 template<
typename Double>
void AD_Const1(Double*,
const IndepADvar<Double>&);
418 template<
typename Double> ADvari<Double>&
ADf1(Double, Double,
const ADvari<Double>&);
419 template<
typename Double> ADvari<Double>&
ADf2(Double, Double, Double,
420 const ADvari<Double>&,
const ADvari<Double>&);
421 template<
typename Double> ADvari<Double>&
ADf2(Double, Double, Double,
422 const IndepADvar<Double>&,
const ADvari<Double>&);
423 template<
typename Double> ADvari<Double>&
ADf2(Double, Double, Double,
424 const ADvari<Double>&,
const IndepADvar<Double>&);
425 template<
typename Double> Double
val(
const ADvari<Double>&);
427 template<
typename Double>
class
433#ifdef RAD_AUTO_AD_Const
443 static ADvari *First_ADvari, **Last_ADvari;
445 inline void ADvari_padv(
const IndepADVar *v) {
447 Last_ADvari = &this->Next;
454 static int gcgen_cur, last_opno, zap_gcgen, zap_gcgen1, zap_opno;
455 static FILE *debug_file;
460 void *
operator new(
size_t len) {
463 rv->gcgen = gcgen_cur;
464 rv->opno = ++last_opno;
465 if (last_opno == zap_opno && gcgen_cur == zap_gcgen)
466 printf(
"Got to opno %d\n", last_opno);
472 void operator delete(
void*) {}
473#ifdef RAD_ALLOW_WANTDERIV
474 inline ADvari(Double t): Val(t), aval(0.), wantderiv(1) {}
475 inline ADvari(): Val(0.), aval(0.), wantderiv(1) {}
476 inline void no_deriv() { wantderiv = 0; }
478 inline ADvari(Double t): Val(t), aval(0.) {}
481#ifdef RAD_AUTO_AD_Const
483 friend class ADvar<Double>;
484 friend class ADvar1<Double>;
486 friend class ADvar2<Double>;
489 ADvari(
const IndepADVar *, Double);
490 ADvari(
const IndepADVar *, Double, Double);
494#define Ai const Base< ADvari >&
496#define T1(r,f) F r f <>(Ai);
530#ifdef HAVE_SACADO_CXX11
550 inline operator Double() {
return this->Val; }
551 inline operator Double()
const {
return this->Val; }
556 template<
typename Double>
class
562 ADvar1(Double val1,
const Double *a1,
const ADVari *c1): ADVari(val1) {
563#ifdef RAD_ALLOW_WANTDERIV
564 if (!ADVari::adc.new_Derp(a1,
this,c1))
567 ADVari::adc.new_Derp(a1,
this,c1);
574 ADVari(val1), d(a1,this,c1) {}
575#ifdef RAD_AUTO_AD_Const
576 typedef typename ADVari::IndepADVar IndepADVar;
578 ADvar1(
const IndepADVar*,
const IndepADVar&);
579 ADvar1(
const IndepADVar*,
const ADVari&);
580 ADvar1(
const Double val1,
const Double *a1,
const ADVari *c1,
const ADVar *v):
581 ADVari(val1), d(a1,this,c1) {
583 this->ADvari_padv(v);
590 template<
typename Double>
class
605 inline void *
operator new(
size_t len) {
return ADVari::adc.Memalloc(len); }
606 inline ConstADvari(Double t): ADVari(t) {}
608 static void aval_reset(
void);
611 template<
typename Double>
class
618 typedef unsigned long ADGenType;
619 mutable ADGenType gen;
625 template<
typename Double>
class
635 IndepADvar_root.ADvprev = (this->ADvprev = IndepADvar_root.ADvprev)->ADvnext =
this;
636 this->ADvnext = &IndepADvar_root;
647 (this->ADvnext->ADvprev = this->ADvprev)->ADvnext = this->ADvnext;
661template<
typename Double> IndepADvar_base0<Double>
662 IndepADvar_base<Double>::IndepADvar_root(0.);
664 template<
typename Double>
void
665IndepADvar_base<Double>::Move_to_end() {
666 if (
this != IndepADvar_root.ADvprev) {
667 (this->ADvnext->ADvprev = this->ADvprev)->ADvnext = this->ADvnext;
668 IndepADvar_root.ADvprev = (this->ADvprev = IndepADvar_root.ADvprev)->ADvnext =
this;
669 this->ADvnext = &IndepADvar_root;
673template<
typename Double>
typename IndepADvar_base0<Double>::ADGenType
674 IndepADvar_base<Double>::gen0(1);
678 template<
typename Double>
class
687#ifdef RAD_AUTO_AD_Const
692#elif defined(RAD_EQ_ALIAS)
703 friend class ADvar<Double>;
705 friend class ADvar1<Double>;
706 friend class ADvarn<Double>;
714#ifdef RAD_ALLOW_WANTDERIV
715 inline int Wantderiv() {
return this->wantderiv; }
718 IndepADvar(
Ttype,
int);
719 IndepADvar(Double,
int);
724#ifdef RAD_AUTO_AD_Const
726 inline IndepADvar() { cv = 0; }
727 inline ~IndepADvar() {
742 inline operator ADVari&()
const {
748 inline operator ADVari*()
const {
749 ADVari *tcv = this->cv;
759 inline Double
val()
const {
766 inline Double
adj()
const {
793#define Ai const Base< ADVari >&
794#define AI const Base< IndepADvar >&
802#define T1(f) friend ADVari& f<> (AI);
804#define F friend ADVari&
843#ifdef HAVE_SACADO_CXX11
856 template<
typename Double>
class
866 void ADvar_ctr(Double d) {
869 if (ConstADVari::cadc.fpval_implies_const)
872#ifdef RAD_AUTO_AD_Const
874 x->ADvari_padv(
this);
888 friend class ADvar1<Double>;
892 inline ADvar(
double i) { ADvar_ctr(Double(
i)); }
894 inline ADvar(
long i) { ADvar_ctr(Double(
i)); }
897#ifdef RAD_AUTO_AD_Const
898 inline ADvar(IndepADVar &
x) {
899 this->cv =
x.cv ?
new ADVar1(
this,
x) : 0;
902 inline ADvar(ADVari &
x) { this->cv = &
x;
x.ADvari_padv(
this); }
903 inline ADvar& operator=(IndepADVar &
x) {
906 this->cv =
new ADVar1(
this,
x);
909 inline ADvar& operator=(ADVari &
x) {
912 this->cv =
new ADVar1(
this,
x);
916 friend ADvar& ADvar_operatoreq<>(ADvar*,
const ADVari&);
919 inline ADvar(
const IndepADVar &
x) {
921 this->cv = (ADVari*)
x.cv;
923 inline ADvar(
const ADVari &
x) { this->cv = (ADVari*)&
x; }
924 inline ADvar& operator=(IndepADVar &
x) {
926 this->cv = (ADVari*)
x.cv;
return *
this;
928 inline ADvar& operator=(
const ADVari &
x) { this->cv = (ADVari*)&
x;
return *
this; }
934 this->cv =
new ADVar1(
x.cv->Val, &this->cv->adc.One,
x.cv);
944 this->cv =
new ADVar1(
x.cv->Val, &this->cv->adc.One, (
ADVari*)
x.cv);
952 this->cv =
new ADVar1(
x.Val, &this->cv->adc.One, &
x);
969 {
return ConstADVari::cadc.fpval_implies_const; }
971 { ConstADVari::cadc.fpval_implies_const = newval; }
973 bool oldval = ConstADVari::cadc.fpval_implies_const;
974 ConstADVari::cadc.fpval_implies_const = newval;
978 inline static bool get_fpval_implies_const(
void) {
return true; }
979 inline static void set_fpval_implies_const(
bool newval) {}
980 inline static bool setget_fpval_implies_const(
bool newval) {
return newval; }
986 static inline void aval_reset() { ConstADVari::aval_reset(); }
994template<
typename Double>
998template<
typename Double>
1001template<
typename Double>
1002 inline void AD_Const(
const IndepADvar<Double>&v) {}
1005 template<
typename Double>
class
1023 void ConstADvar_ctr(Double);
1033#ifdef RAD_NO_CONST_UPDATE
1039 this->cv =
new ADVari(d);
1048 this->cv =
new ADVar1(
this,d);
1057#ifdef RAD_ALLOW_WANTDERIV
1058 template<
typename Double>
class
1062 typedef ADvar<Double>
ADVar;
1063 typedef IndepADvar<Double> IndepADVar;
1064 typedef typename IndepADVar::ADVari ADVari;
1065 typedef ADvar1<Double> ADVar1;
1067 void ADvar_ndctr(Double d) {
1068 ADVari *
x =
new ADVari(d);
1076 inline ADvar_nd(
double i):
ADVar((void*)0,0) { ADvar_ndctr(Double(
i)); }
1078 inline ADvar_nd(
long i):
ADVar((void*)0,0) { ADvar_ndctr(Double(
i)); }
1079 inline ADvar_nd(Double d):
ADVar((void*)0,0) { ADvar_ndctr(d); }
1080 inline ~ADvar_nd() {}
1083 this->cv =
new ADVar1(
x.cv->Val, &this->cv->adc.One,
x.cv);
1090 this->cv =
new ADVar1(
x.Val, &this->cv->adc.One, &
x);
1096#define ADvar_nd ADvar
1099 template<
typename Double>
class
1109#ifdef RAD_AUTO_AD_Const
1110 typedef typename ADVar1::ADVar
ADVar;
1111 ADvar1s(Double val1, Double a1,
const ADVari *c1,
const ADVar *v):
1112 ADVar1(val1,&
a,c1,v),
a(a1) {}
1117 template<
typename Double>
class
1124 ADvar2(Double val1,
const ADVari *Lcv,
const Double *Lc,
1125 const ADVari *Rcv,
const Double *Rc): ADVari(val1) {
1126#ifdef RAD_ALLOW_WANTDERIV
1128 Ld = ADVari::adc.new_Derp(Lc,
this, Lcv);
1129 Rd = ADVari::adc.new_Derp(Rc,
this, Rcv);
1133 ADVari::adc.new_Derp(Lc,
this, Lcv);
1134 ADVari::adc.new_Derp(Rc,
this, Rcv);
1140 const ADVari *Rcv,
const Double *Rc):
1142 dR.
next = DErp::LastDerp;
1144 DErp::LastDerp = &dL;
1151#ifdef RAD_AUTO_AD_Const
1153 ADvar2(Double val1,
const ADVari *Lcv,
const Double *Lc,
1154 const ADVari *Rcv,
const Double *Rc,
ADVar *v):
1156 dR.
next = DErp::LastDerp;
1158 DErp::LastDerp = &dL;
1165 this->ADvari_padv(v);
1171 template<
typename Double>
class
1180#ifdef RAD_ALLOW_WANTDERIV
1182 Ld = ADVari::adc.new_Derp(&Lp,
this, Lcv);
1183 Rd = ADVari::adc.new_Derp(&Rp,
this, Rcv);
1187 ADVari::adc.new_Derp(&Lp,
this, Lcv);
1188 ADVari::adc.new_Derp(&Rp,
this, Rcv);
1195 this->dR.next = DErp::LastDerp;
1196 this->dL.next = &this->dR;
1197 DErp::LastDerp = &this->dL;
1202 this->dL.b = this->dR.b =
this;
1204#ifdef RAD_AUTO_AD_Const
1205 typedef typename ADVar2::ADVar
ADVar;
1206 ADvar2q(Double val1, Double Lp, Double Rp,
const ADVari *Lcv,
1207 const ADVari *Rcv,
const ADVar *v):
1208 ADVar2(val1),
a(Lp), b(Rp) {
1209 this->dR.next = DErp::LastDerp;
1210 this->dL.next = &this->dR;
1211 DErp::LastDerp = &this->dL;
1216 this->dL.b = this->dR.b =
this;
1218 this->ADvari_padv(v);
1224 template<
typename Double>
class
1232#ifdef RAD_ALLOW_WANTDERIV
1234 for(
i = nd = 0;
i < n1;
i++)
1235 if (ADVari::adc.new_Derp(g+
i,
this,
x[
i].cv))
1240 for(
int i = 0;
i < n1;
i++)
1241 ADVari::adc.new_Derp(g+
i,
this,
x[
i].cv);
1254 a1 =
a = (Double*)ADVari::adc.Memalloc(n*
sizeof(*
a));
1255 d1 = Da = (
DErp*)ADVari::adc.Memalloc(n*
sizeof(
DErp));
1256 dlast = DErp::LastDerp;
1257 for(
i = 0;
i < n1;
i++, d1++) {
1265 DErp::LastDerp = dlast;
1270template<
typename Double>
1275template<
typename Double>
1279 return L.
Val <
R.Val;
1281template<
typename Double>
1286template<
typename Double>
1292template<
typename Double>
1296 return L.
Val <=
R.Val;
1298template<
typename Double>
1303template<
typename Double>
1309template<
typename Double>
1313 return L.
Val ==
R.Val;
1315template<
typename Double>
1320template<
typename Double>
1326template<
typename Double>
1330 return L.
Val !=
R.Val;
1332template<
typename Double>
1337template<
typename Double>
1343template<
typename Double>
1347 return L.
Val >=
R.Val;
1349template<
typename Double>
1354template<
typename Double>
1360template<
typename Double>
1364 return L.
Val >
R.Val;
1366template<
typename Double>
1371template<
typename Double>
1377template<
typename Double>
1380 return Mbase + (Mleft -= len);
1381 return new_ADmemblock(len);
1385template<
typename Double>
1388template<
typename Double>
1392template<
typename Double>
1394 a(*a1), b(b1),
c(c1) {}
1397template<
typename Double>
1403template<
typename Double>
1410template<
typename Double>
1429#ifdef RAD_AUTO_AD_Const
1435#ifndef RAD_DEBUG_gcgen1
1436#define RAD_DEBUG_gcgen1 -1
1452 Mbase = (
char*)First->
memblk;
1453 Mleft =
sizeof(First->
memblk);
1454 rad_need_reinit = 0;
1455#ifdef RAD_DEBUG_BLOCKKEEP
1456 rad_busy_blocks = 0;
1464 DMleft = nderps =
sizeof(DBusy->
memblk)/
sizeof(
DErp);
1473 for(mb = ADVari::adc.Busy; mb; mb = mb1) {
1477 for(mb = ADVari::adc.Free; mb; mb = mb1) {
1481 for(mb = ConstADVari::cadc.Busy; mb; mb = mb1) {
1485 ConstADVari::cadc.Busy = ADVari::adc.Busy = ADVari::adc.Free = 0;
1486 ConstADVari::cadc.Mleft = ADVari::adc.Mleft = 0;
1487 ConstADVari::cadc.Mbase = ADVari::adc.Mbase = 0;
1489 for(mb = ADVari::adc.DBusy; mb; mb = mb1) {
1493 for(mb = ADVari::adc.DFree; mb; mb = mb1) {
1497 ADVari::adc.DBusy = ADVari::adc.DFree = 0;
1498 ADVari::adc.DMleft = 0;
1499 ConstADVari::cadc.Mbase = ADVari::adc.Mbase = 0;
1501 ConstADVari::lastcad = 0;
1510 if (ConstADVari::cadc.Busy || ADVari::adc.Busy || ADVari::adc.Free
1512 || ADVari::adc.DBusy || ADVari::adc.DFree
1515 ADVari::adc.do_init();
1516 ConstADVari::cadc.do_init();
1519template<
typename Double>
void*
1523#ifdef RAD_AUTO_AD_Const
1526#ifdef RAD_Const_WARN
1538 if ((rad_need_reinit & 1) &&
this == &ADVari::adc) {
1539 rad_need_reinit &= ~1;
1541#ifdef RAD_DEBUG_BLOCKKEEP
1542 Mleft = rad_mleft_save;
1543 if (Mleft <
sizeof(First->
memblk))
1545 UninitType<Double>::utype,
1546 (
sizeof(First->
memblk) - Mleft)
1548 if ((mb = Busy->
next)) {
1550 for(;; mb = mb->
next) {
1552 UninitType<Double>::utype,
1559 rad_Oldcurmb = Busy;
1560 if (rad_busy_blocks >= RAD_DEBUG_BLOCKKEEP) {
1561 rad_busy_blocks = 0;
1565 for(mb = Busy; mb != mb0; mb = mb1) {
1572 Mbase = (
char*)First->
memblk;
1573 Mleft =
sizeof(First->
memblk);
1580 for(mb = Busy; mb != mb0; mb = mb1) {
1587 Mbase = (
char*)First->
memblk;
1588 Mleft =
sizeof(First->
memblk);
1590#ifdef RAD_AUTO_AD_Const
1591 *ADVari::Last_ADvari = 0;
1592 ADVari::Last_ADvari = &ADVari::First_ADvari;
1593 a = ADVari::First_ADvari;
1598#ifdef RAD_Const_WARN
1599 if ((
i =
a->opno) > 0)
1615 while((mb1 = mb->
next)) {
1626 while((vb = vb->ADvnext) != vb0)
1627 if ((tcv = ((IADv*)vb)->cv))
1629#elif RAD_REINIT == 2
1634 return Mbase + (Mleft -= len);
1641#ifdef RAD_DEBUG_BLOCKKEEP
1646 return (Mbase = (
char*)
x->memblk) +
1647 (Mleft =
sizeof(First->
memblk) - len);
1650template<
typename Double>
void
1657 if (ADVari::adc.rad_need_reinit && wantgrad) {
1658 mb = ADVari::adc.DBusy;
1659 d = ((
DErp*)mb->
memblk) + ADVari::adc.DMleft;
1660 de = ((
DErp*)mb->
memblk) + ADVari::adc.nderps;
1664 if (!(mb = mb->
next))
1667 de = d + ADVari::adc.nderps;
1673 if (ADVari::adc.rad_need_reinit && wantgrad) {
1674 for(d = DErp::LastDerp; d; d = d->
next)
1679 if (!(ADVari::adc.rad_need_reinit & 1)) {
1680 ADVari::adc.rad_need_reinit = 1;
1681 ADVari::adc.rad_mleft_save = ADVari::adc.Mleft;
1682 ADVari::adc.Mleft = 0;
1686 if (ADVari::gcgen_cur == ADVari::zap_gcgen1 && wantgrad) {
1688 if (!(fname = getenv(
"RAD_DEBUG_FILE")))
1689 fname =
"rad_debug.out";
1693 ADVari::debug_file = fopen(fname,
"w");
1694 ADVari::zap_gcgen1 = -1;
1698 if (ADVari::adc.DMleft < ADVari::adc.nderps && wantgrad) {
1699 mb = ADVari::adc.DBusy;
1700 d = ((
DErp*)mb->
memblk) + ADVari::adc.DMleft;
1701 de = ((
DErp*)mb->
memblk) + ADVari::adc.nderps;
1705 if (ADVari::debug_file) {
1706 for(; d < de; d++) {
1707 fprintf(ADVari::debug_file,
"%d\t%d\t%g + %g * %g",
1710 fprintf(ADVari::debug_file,
" = %g\n", d->
c->
aval);
1717 if (!(mb = mb->
next))
1720 de = d + ADVari::adc.nderps;
1724 if ((d = DErp::LastDerp) && wantgrad) {
1727 if (ADVari::debug_file)
1729 fprintf(ADVari::debug_file,
"%d\t%d\t%g + %g * %g",
1732 fprintf(ADVari::debug_file,
" = %g\n", d->
c->
aval);
1733 }
while((d = d->
next));
1737 while((d = d->
next));
1740 if (ADVari::debug_file) {
1741 fclose(ADVari::debug_file);
1742 ADVari::debug_file = 0;
1747 if (ADVari::debug_file)
1748 fflush(ADVari::debug_file);
1749 ADVari::gcgen_cur++;
1750 ADVari::last_opno = 0;
1754 template<
typename Double>
void
1762 if (ADVari::adc.rad_need_reinit) {
1763 mb = ADVari::adc.DBusy;
1764 d = ((
DErp*)mb->
memblk) + ADVari::adc.DMleft;
1765 de = ((
DErp*)mb->
memblk) + ADVari::adc.nderps;
1769 if (!(mb = mb->
next))
1772 de = d + ADVari::adc.nderps;
1778 if (ADVari::adc.rad_need_reinit) {
1779 for(d = DErp::LastDerp; d; d = d->
next)
1784 if (!(ADVari::adc.rad_need_reinit & 1)) {
1785 ADVari::adc.rad_need_reinit = 1;
1786 ADVari::adc.rad_mleft_save = ADVari::adc.Mleft;
1787 ADVari::adc.Mleft = 0;
1791 if (ADVari::gcgen_cur == ADVari::zap_gcgen1) {
1793 if (!(fname = getenv(
"RAD_DEBUG_FILE")))
1794 fname =
"rad_debug.out";
1798 ADVari::debug_file = fopen(fname,
"w");
1799 ADVari::zap_gcgen1 = -1;
1803 if (ADVari::adc.DMleft < ADVari::adc.nderps) {
1804 for(
i = 0;
i < n;
i++)
1806 mb = ADVari::adc.DBusy;
1807 d = ((
DErp*)mb->
memblk) + ADVari::adc.DMleft;
1808 de = ((
DErp*)mb->
memblk) + ADVari::adc.nderps;
1812 if (ADVari::debug_file) {
1813 for(; d < de; d++) {
1814 fprintf(ADVari::debug_file,
"%d\t%d\t%g + %g * %g",
1817 fprintf(ADVari::debug_file,
" = %g\n", d->
c->
aval);
1824 if (!(mb = mb->
next))
1827 de = d + ADVari::adc.nderps;
1831 if ((d = DErp::LastDerp) != 0) {
1832 for(
i = 0;
i < n;
i++)
1835 if (ADVari::debug_file)
1837 fprintf(ADVari::debug_file,
"%d\t%d\t%g + %g * %g",
1840 fprintf(ADVari::debug_file,
" = %g\n", d->
c->
aval);
1841 }
while((d = d->
next));
1845 while((d = d->
next));
1848 if (ADVari::debug_file) {
1849 fclose(ADVari::debug_file);
1850 ADVari::debug_file = 0;
1855 if (ADVari::debug_file)
1856 fflush(ADVari::debug_file);
1857 ADVari::gcgen_cur++;
1858 ADVari::last_opno = 0;
1862 template<
typename Double>
void
1870 template<
typename Double>
void
1873 for(
DErp *d = DErp::LastDerp; d; d = d->next) {
1875 *(
const_cast<Double*
>(d->a)) = Double(0.0);
1877 d->b->aval = Double(0.0);
1878 d->b->Val = Double(0.0);
1881 d->c->aval = Double(0.0);
1882 d->c->Val = Double(0.0);
1887 template<
typename Double>
1892 RAD_REINIT_2(Val = d; this->gen = this->IndepADvar_root.gen;)
1895 template<
typename Double>
1900 RAD_REINIT_2(Val =
i; this->gen = this->IndepADvar_root.gen;)
1903 template<
typename Double>
1908 RAD_REINIT_2(Val =
i; this->gen = this->IndepADvar_root.gen;)
1911 template<
typename Double>
1916 RAD_REINIT_2(Val =
i; this->gen = this->IndepADvar_root.gen;)
1919 template<
typename Double>
1924 RAD_REINIT_2(this->Val = 0.; this->gen = this->IndepADvar_root.gen;)
1927 template<
typename Double>
void
1932 RAD_REINIT_2(this->Val = d; this->gen = this->IndepADvar_root.gen;)
1935 template<
typename Double>
1942 x.cv->adc.new_Derp(&
x.adc.One,
y,
x.cv);
1947 RAD_REINIT_2(this->Val =
y->Val; this->gen = this->IndepADvar_root.gen;)
1950 template<
typename Double>
1957 x.cv->adc.new_Derp(&
x.cv->adc.One,
y,
x.cv);
1962 RAD_REINIT_2(this->Val =
y->Val; this->gen = this->IndepADvar_root.gen;)
1965 template<
typename Double>
1971 x.adc.new_Derp(&
x.adc.One,
y, &
x);
1976 RAD_REINIT_2(this->Val =
y->Val; this->gen = this->IndepADvar_root.gen;)
1979 template<
typename Double>
1985 ConstADVari *ncv =
new ConstADVari(v.
val());
1986#ifdef RAD_AUTO_AD_Const
1993 template<
typename Double>
1997#ifdef RAD_ALLOW_WANTDERIV
1999 if (this->gen != this->IndepADvar_root.gen) {
2001 this->gen = this->IndepADvar_root.gen;
2004 return this->wantderiv = cv->wantderiv = n;
2010 template<
typename Double>
2020#elif RAD_REINIT == 1
2025#ifdef RAD_AUTO_AD_Const
2027 template<
typename Double>
2030 this->ADvari_padv(
x);
2034 template<
typename Double>
2037 this->ADvari_padv(
x);
2041 template<
typename Double>
2043 ADVari(
y.cv->Val), d((const Double*)&ADcontext<Double>::One, (ADVari*)this,
y.cv)
2045 this->ADvari_padv(
x);
2048 template<
typename Double>
2050 ADVari(
y.Val), d((const Double*)&ADcontext<Double>::One, this, &
y)
2052 this->ADvari_padv(
x);
2057 template<
typename Double>
2063 RAD_REINIT_2(This->Val =
x.Val; This->gen = This->IndepADvar_root.gen;)
2067 template<
typename Double>
2073 RAD_REINIT_2(This->Val =
x.Val; This->gen = This->IndepADvar_root.gen;)
2080 template<
typename Double>
2084#ifdef RAD_AUTO_AD_Const
2087 this->cv =
new ADVari(
this,d);
2089 this->cv =
new ADVari(d);
2096 template<
typename Double>
2100#ifdef RAD_AUTO_AD_Const
2103 this->cv =
new ADVari(
this,d);
2105 this->cv =
RAD_REINIT_0(ConstADVari::cadc.fpval_implies_const
2109 RAD_REINIT_2(this->Val = d; this->gen = this->IndepADvar_root.gen;)
2114 template<
typename Double>
2121 template<
typename Double>
2129#ifdef RAD_AUTO_AD_Const
2130#define RAD_ACA ,this
2135 template<
typename Double>
2146 template<
typename Double>
2153 template<
typename Double>
2164 template<
typename Double>
2171 template<
typename Double>
2179 template<
typename Double>
2190 template<
typename Double>
2197 template<
typename Double>
2208 template<
typename Double>
2215 template<
typename Double>
2223 template<
typename Double>
2234 template<
typename Double>
2241 template<
typename Double>
2252 template<
typename Double>
2259 template<
typename Double>
2264 Double Lv = L.
Val, Rv =
R.Val, pL = 1. / Rv, q = Lv/Rv;
2268 template<
typename Double>
2272 Double Lv = Lcv->
Val, Rv =
R.Val, pL = 1. / Rv, q = Lv/Rv;
2280 template<
typename Double>
2287 template<
typename Double>
2291 Double recip = 1. /
R.Val;
2292 Double q = L * recip;
2296 template<
typename Double>
2307 template<
typename Double>
2312 return *(
new ADvar1s<Double>(std::acos(t), -1./std::sqrt(1. - t*t), &v));
2315 template<
typename Double>
2319 Double t = v.
Val, t1 = std::sqrt(t*t - 1.);
2323 template<
typename Double>
2328 return *(
new ADvar1s<Double>(std::asin(t), 1./std::sqrt(1. - t*t), &v));
2331 template<
typename Double>
2335 Double t = v.
Val, td = 1., t1 = std::sqrt(t*t + 1.);
2343 template<
typename Double>
2351 template<
typename Double>
2356 return *(
new ADvar1s<Double>(0.5*std::log((1.+t)/(1.-t)), 1./(1. - t*t), &v));
2359 template<
typename Double>
2368 template<
typename Double>
2372 Double
y =
R.Val, t =
x*
x +
y*
y;
2376 template<
typename Double>
2384 template<
typename Double>
2393 template<
typename Double>
2402 template<
typename Double>
2411 template<
typename Double>
2420 template<
typename Double>
2429 template<
typename Double>
2438 template<
typename Double>
2445 template<
typename Double>
2452 template<
typename Double>
2457 Double t = std::exp(v.
Val);
2461 rcv->
d.a = &rcv->
Val;
2467 template<
typename Double>
2475 template<
typename Double>
2479 static double num = 1. / std::log(10.);
2484 template<
typename Double>
2489 Double
x = L.
Val,
y =
R.Val, t = std::pow(
x,
y);
2493 template<
typename Double>
2497 Double t = std::pow(
x,
R.Val);
2501 template<
typename Double>
2505 Double
x = L.
Val, t = std::pow(
x,
y);
2509 template<
typename Double>
2516 template<
typename Double>
2523 template<
typename Double>
2527 Double t = std::sqrt(v.
Val);
2531 template<
typename Double>
2535 Double t = std::cos(v.
Val);
2539 template<
typename Double>
2543 Double t = 1. / std::cosh(v.
Val);
2547 template<
typename Double>
2553 if ((t = v.
Val) < 0) {
2560 template<
typename Double>
2569 if ((t = v.
Val) < 0) {
2576#ifdef HAVE_SACADO_CXX11
2577template<
typename Double>
2579cbrt(
const Base< ADvari<Double> > &vv) {
2580 const ADvari<Double>& v = vv.derived();
2581 Double t = std::cbrt(v.Val);
2582 return *(
new ADvar1s<Double>(t, 1.0/(3.0*t*t), &v));
2586 template<
typename Double>
2592 template<
typename Double>
2598 template<
typename Double>
2604 template<
typename Double>
2610 template<
typename Double>
2616 template<
typename Double>
2622 template<
typename Double>
2628 template<
typename Double>
2634 template<
typename Double>
2641#define A (ADvari<Double>*)
2642#ifdef RAD_Const_WARN
2643#define C(x) (((x)->opno < 0) ? RAD_Const_Warn(x) : 0, *A x)
2647#define T template<typename Double> inline
2648#define F ADvari<Double>&
2649#define Ai const Base< ADvari<Double> >&
2650#define AI const Base< IndepADvar<Double> >&
2652#define CAI(x,y) const IndepADvar<Double> & x = y.derived()
2653#define CAi(x,y) const ADvari<Double> & x = y.derived()
2655 T r f(Ai LL, AI RR) { CAi(L,LL); CAI(R,RR); RAD_cvchk(R) return f(L, C(R.cv)); } \
2656 T r f(AI LL, Ai RR) { CAI(L,LL); CAi(R,RR); RAD_cvchk(L) return f(C(L.cv), R); }\
2657 T r f(AI LL, AI RR) { CAI(L,LL); CAI(R,RR); RAD_cvchk(L) RAD_cvchk(R) return f(C(L.cv), C(R.cv)); }\
2658 T r f(AI LL, D R) { CAI(L,LL); RAD_cvchk(L) return f(C(L.cv), R); } \
2659 T r f(D L, AI RR) { CAI(R,RR); RAD_cvchk(R) return f(L, C(R.cv)); } \
2660 T r f(Ai L, Dtype R) { return f(L, (D)R); }\
2661 T r f(AI L, Dtype R) { return f(L, (D)R); }\
2662 T r f(Ai L, Ltype R) { return f(L, (D)R); }\
2663 T r f(AI L, Ltype R) { return f(L, (D)R); }\
2664 T r f(Ai L, Itype R) { return f(L, (D)R); }\
2665 T r f(AI L, Itype R) { return f(L, (D)R); }\
2666 T r f(Dtype L, Ai R) { return f((D)L, R); }\
2667 T r f(Dtype L, AI R) {return f((D)L, R); }\
2668 T r f(Ltype L, Ai R) { return f((D)L, R); }\
2669 T r f(Ltype L, AI R) { return f((D)L, R); }\
2670 T r f(Itype L, Ai R) { return f((D)L, R); }\
2671 T r f(Itype L, AI R) { return f((D)L, R); }
2692 T F f(AI xx) { CAI(x,xx); RAD_cvchk(x) return f(C(x.cv)); }
2714#ifdef HAVE_SACADO_CXX11
int RAD_Const_Warn(void *v)
expr expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c *expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr2 expr1 expr2 expr1 expr1 expr1 c
ADmemblock< Double > ADMemblock
static void Gradcomp(int wantgrad)
static void Weighted_Gradcomp(size_t, ADVar **, Double *)
void * Memalloc(size_t len)
static void Outvar_Gradcomp(ADVar &)
void * new_ADmemblock(size_t)
static const Double negOne
ADvar1(Double val1, const Double *a1, const ADVari *c1)
ADvar1(Double val1, const ADVari *c1)
ADvar1s(Double val1, Double a1, const ADVari *c1)
ADvar2(Double val1, const ADVari *Lcv, const Double *Lc, const ADVari *Rcv, const Double *Rc)
ADvar2q(Double val1, Double Lp, Double Rp, const ADVari *Lcv, const ADVari *Rcv)
ConstADvari< Double > ConstADVari
IndepADvar< Double > IndepADVar
static void Outvar_Gradcomp(ADvar &v)
ADvar & operator-=(const ADVari &)
static void set_fpval_implies_const(bool newval)
ADvar & operator+=(Double)
ADvar & operator/=(Double)
ADvar & operator-=(Double)
static void Gradcomp(int wantgrad)
ADvar & operator=(const ADVari &x)
static bool get_fpval_implies_const(void)
IndepADVar::ADVari ADVari
ADvar & operator*=(const ADVari &)
ADvar & operator+=(const ADVari &)
const ADVari & ADvar(const IndepADVar &x)
static bool setget_fpval_implies_const(bool newval)
Allow_noderiv(inline ADvar(void *v, int wd):IndepADVar(v, wd) {}) friend ADvar &ADvar_operatoreq<>(ADvar *
ADvar & operator=(Double)
static void Weighted_Gradcomp(size_t n, ADvar **v, Double *w)
ADvar & operator*=(Double)
ADvar & operator/=(const ADVari &)
Allow_noderiv(mutable int wantderiv;) void *operator new(size_t len)
IndepADvar< Double > IndepADVar
static ADcontext< Double > adc
ScalarType< value_type >::type scalar_type
ADVari::IndepADVar IndepADVar
ADvarn(Double val1, int n1, const IndepADVar *x, const Double *g)
ConstADvar & operator=(Double d)
ConstADvar & operator*=(ADVari &)
ADVar::ConstADVari ConstADVari
ConstADvar & operator*=(Double)
void ConstADvar_ctr(Double)
ConstADvar & operator+=(ADVari &)
ConstADvar & operator-=(ADVari &)
ConstADvar & operator+=(Double)
ConstADvar & operator=(ADVari &d)
ConstADvar & operator-=(Double)
ADVar::IndepADVar IndepADVar
ConstADvar & operator/=(Double)
ConstADvar & operator/=(ADVari &)
static CADcontext< Double > cadc
static void aval_reset(void)
IndepADvar_base(Allow_noderiv(int wd))
IndepADvar() Allow_noderiv(
static void AD_Const(const IndepADvar &)
static void Weighted_Gradcomp(size_t n, ADVar **v, Double *w)
static void Gradcomp(int wantgrad)
IndepADvar & operator=(IndepADvar &x)
static void Outvar_Gradcomp(ADVar &v)
ADvari< Double > & atan(const Base< ADvari< Double > > &vv)
ADvari< Double > & operator*(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
ADvari< Double > & tanh(const Base< ADvari< Double > > &vv)
int operator==(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
ADvari< Double > & ADf1(Double f, Double g, const IndepADvar< Double > &x)
ADvari< Double > & sin(const Base< ADvari< Double > > &vv)
int operator<(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
void AD_Const1(Double *, const IndepADvar< Double > &)
int operator!=(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
IndepADvar< Double > & ADvar_operatoreq(IndepADvar< Double > *, const ADvari< Double > &)
ADvari< Double > & tan(const Base< ADvari< Double > > &vv)
int operator>(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
ADvari< Double > & operator-(const Base< ADvari< Double > > &TT)
ADvari< Double > & operator+(const Base< ADvari< Double > > &TT)
ADvari< Double > & asin(const Base< ADvari< Double > > &vv)
ADvari< Double > & log10(const Base< ADvari< Double > > &vv)
ADvari< Double > & sinh(const Base< ADvari< Double > > &vv)
ADvari< Double > & cos(const Base< ADvari< Double > > &vv)
ADvari< Double > & fabs(const Base< ADvari< Double > > &vv)
ADvari< Double > & max(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
void AD_Const(const IndepADvar< Double > &)
ADvari< Double > & ADf2(Double f, Double gx, Double gy, const IndepADvar< Double > &x, const IndepADvar< Double > &y)
ADvari< Double > & operator/(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
ADvari< Double > & log(const Base< ADvari< Double > > &vv)
ADvari< Double > & pow(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
ADvari< Double > & acos(const Base< ADvari< Double > > &vv)
int operator<=(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
ADvari< Double > & ADfn(Double f, int n, const IndepADvar< Double > *x, const Double *g)
ADvari< Double > & atanh(const Base< ADvari< Double > > &vv)
ADvari< Double > & min(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
ADvari< Double > & atan2(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
ADvari< Double > & cosh(const Base< ADvari< Double > > &vv)
ADvari< Double > & asinh(const Base< ADvari< Double > > &vv)
ADvari< Double > & acosh(const Base< ADvari< Double > > &vv)
int operator>=(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
ADvari< Double > & sqrt(const Base< ADvari< Double > > &vv)
ADvari< Double > & abs(const Base< ADvari< Double > > &vv)
ADvari< Double > & exp(const Base< ADvari< Double > > &vv)
Base class for Sacado types to control overload resolution.
const derived_type & derived() const
char memblk[1000 *sizeof(double)]
Turn ADvar into a meta-function class usable with mpl::apply.
ADT_RAD IndepADvar< double > AI
ADT_RAD ADvari< double > Ai
Sacado::RadVec::ADvar< double > ADVar
void _uninit_f2c(void *x, int type, long len)