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) 67 #if RAD_REINIT == 1 //{{ 69 #define RAD_REINIT_1(x) x 70 #ifdef RAD_ALLOW_WANTDERIV 71 #define ADvar ADvar_1n 72 #define IndepADvar IndepADvar_1n 75 #define IndepADvar IndepADvar_1 76 #endif //RAD_ALLOW_WANTDERIV 77 #elif RAD_REINIT == 2 //}{ 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 85 #define ADvar ADvar_2n 87 #define IndepADvar IndepADvar_2 89 #endif //RAD_ALLOW_WANTDERIV 90 #elif RAD_REINIT != 0 //}{ 91 Botch! 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) 106 #undef RAD_Const_WARN 107 #undef RAD_AUTO_AD_Const 108 #undef RAD_DEBUG_BLOCKKEEP 114 #ifdef RAD_Const_WARN // ==> RAD_AUTO_AD_Const and RAD_DEBUG 115 #ifndef RAD_AUTO_AD_Const 116 #define RAD_AUTO_AD_Const 123 #endif // RAD_Const_WARN 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 146 extern "C" void _uninit_f2c(
void *
x,
int type,
long len);
148 template <
typename T>
149 struct UninitType {};
152 struct UninitType<float> {
153 static const int utype = 4;
157 struct UninitType<double> {
158 static const int utype = 5;
161 template <
typename T>
162 struct 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 213 template<
typename Double>
class ADvar;
220 template<
typename Double>
class Derp;
222 template<
typename Double>
struct 230 template<
typename Double>
class 243 size_t DMleft, nderps;
245 #ifdef RAD_DEBUG_BLOCKKEEP 249 void *new_ADmemblock(
size_t);
254 void *Memalloc(
size_t len);
255 static void Gradcomp(
int wantgrad);
257 static void aval_reset();
258 static void free_all();
259 static void re_init();
260 static void Weighted_Gradcomp(
size_t,
ADVar**,
Double*);
261 static void Outvar_Gradcomp(
ADVar&);
263 DErp *new_Derp(
const Double *,
const ADVari*,
const ADVari*);
269 template<
typename Double>
class 279 template<
typename Double>
class 286 inline void *
operator new(size_t,
Derp *
x) {
return x; }
287 inline void operator delete(
void*,
Derp *) {}
296 Derp(
const ADVari *);
298 Derp(
const Double *,
const ADVari *,
const ADVari *);
299 inline void *
operator new(
size_t len) {
return (
Derp*)ADVari::adc.Memalloc(len); }
303 #if RAD_REINIT > 0 //{ 304 template<
typename Double> Derp<Double>*
305 ADcontext<Double>::new_Derp(
const Double *
a,
const ADvari<Double> *b,
const ADvari<Double> *
c)
310 if (this->DMleft == 0) {
317 this->DMleft = nderps;
319 rv = &((DErp*)DBusy->
memblk)[--this->DMleft];
327 #define Ai const Base< ADvari<Double> >& 328 #define AI const Base< IndepADvar<Double> >& 329 #define T template<typename Double> 356 #define F ADvari<Double>& 391 #ifdef HAVE_SACADO_CXX11 406 template<
typename Double>ADvari<Double>&
ADf1(
Double f,
Double g,
const IndepADvar<Double> &
x);
408 const IndepADvar<Double> &
x,
const IndepADvar<Double> &
y);
409 template<
typename Double>ADvari<Double>&
ADfn(
Double f,
int n,
410 const IndepADvar<Double> *
x,
const Double *g);
412 template<
typename Double> IndepADvar<Double>&
ADvar_operatoreq(IndepADvar<Double>*,
413 const ADvari<Double>&);
414 template<
typename Double> ADvar<Double>&
ADvar_operatoreq(ADvar<Double>*,
const ADvari<Double>&);
415 template<
typename Double>
void AD_Const(
const IndepADvar<Double>&);
416 template<
typename Double>
void AD_Const1(
Double*,
const IndepADvar<Double>&);
417 template<
typename Double> ADvari<Double>&
ADf1(
Double,
Double,
const ADvari<Double>&);
419 const ADvari<Double>&,
const ADvari<Double>&);
421 const IndepADvar<Double>&,
const ADvari<Double>&);
423 const ADvari<Double>&,
const IndepADvar<Double>&);
424 template<
typename Double>
Double val(
const ADvari<Double>&);
426 template<
typename Double>
class 427 ADvari :
public Base< ADvari<Double> > {
432 #ifdef RAD_AUTO_AD_Const 434 #ifdef RAD_Const_WARN 439 #endif //RAD_Const_WARN 442 static ADvari *First_ADvari, **Last_ADvari;
444 inline void ADvari_padv(
const IndepADVar *v) {
446 Last_ADvari = &this->Next;
449 #endif //RAD_AUTO_AD_Const 453 static int gcgen_cur, last_opno, zap_gcgen, zap_gcgen1, zap_opno;
454 static FILE *debug_file;
459 void *
operator new(
size_t len) {
462 rv->gcgen = gcgen_cur;
463 rv->opno = ++last_opno;
464 if (last_opno == zap_opno && gcgen_cur == zap_gcgen)
465 printf(
"Got to opno %d\n", last_opno);
471 void operator delete(
void*) {}
472 #ifdef RAD_ALLOW_WANTDERIV 473 inline ADvari(
Double t): Val(t), aval(0.), wantderiv(1) {}
474 inline ADvari(): Val(0.), aval(0.), wantderiv(1) {}
475 inline void no_deriv() { wantderiv = 0; }
480 #ifdef RAD_AUTO_AD_Const 493 #define Ai const Base< ADvari >& 495 #define T1(r,f) F r f <>(Ai); 529 #ifdef HAVE_SACADO_CXX11 547 friend ADvari& ADfn<>(
Double f,
int n,
const IndepADVar *
x,
const Double *g);
549 inline operator Double() {
return this->Val; }
550 inline operator Double()
const {
return this->Val; }
555 template<
typename Double>
class 562 #ifdef RAD_ALLOW_WANTDERIV 563 if (!ADVari::adc.new_Derp(a1,
this,c1))
566 ADVari::adc.new_Derp(a1,
this,c1);
569 #else // RAD_REINIT == 0 573 ADVari(val1), d(a1,this,c1) {}
574 #ifdef RAD_AUTO_AD_Const 575 typedef typename ADVari::IndepADVar IndepADVar;
577 ADvar1(
const IndepADVar*,
const IndepADVar&);
578 ADvar1(
const IndepADVar*,
const ADVari&);
580 ADVari(val1), d(a1,this,c1) {
582 this->ADvari_padv(v);
585 #endif // RAD_REININT > 0 589 template<
typename Double>
class 590 ConstADvari:
public ADvari<Double> {
604 inline void *
operator new(
size_t len) {
return ADVari::adc.Memalloc(len); }
605 inline ConstADvari(
Double t): ADVari(t) {}
607 static void aval_reset(
void);
610 template<
typename Double>
class 616 #elif RAD_REINIT == 2 617 typedef unsigned long ADGenType;
618 mutable ADGenType gen;
624 template<
typename Double>
class 634 IndepADvar_root.ADvprev = (this->ADvprev = IndepADvar_root.ADvprev)->ADvnext =
this;
635 this->ADvnext = &IndepADvar_root;
646 (this->ADvnext->ADvprev = this->ADvprev)->ADvnext = this->ADvnext;
659 #if RAD_REINIT > 0 //{ 660 template<
typename Double> IndepADvar_base0<Double>
661 IndepADvar_base<Double>::IndepADvar_root(0.);
663 template<
typename Double>
void 664 IndepADvar_base<Double>::Move_to_end() {
665 if (
this != IndepADvar_root.ADvprev) {
666 (this->ADvnext->ADvprev = this->ADvprev)->ADvnext = this->ADvnext;
667 IndepADvar_root.ADvprev = (this->ADvprev = IndepADvar_root.ADvprev)->ADvnext =
this;
668 this->ADvnext = &IndepADvar_root;
671 #elif RAD_REINIT == 2 672 template<
typename Double>
typename IndepADvar_base0<Double>::ADGenType
673 IndepADvar_base<Double>::gen0(1);
677 template<
typename Double>
class 678 IndepADvar:
protected IndepADvar_base<Double>,
public Base< IndepADvar<Double> > {
680 static void AD_Const(
const IndepADvar&);
686 #ifdef RAD_AUTO_AD_Const 691 #elif defined(RAD_EQ_ALIAS) 696 #endif //RAD_AUTO_AD_Const 713 #ifdef RAD_ALLOW_WANTDERIV 714 inline int Wantderiv() {
return this->wantderiv; }
717 IndepADvar(
Ttype,
int);
723 #ifdef RAD_AUTO_AD_Const 725 inline IndepADvar() { cv = 0; }
726 inline ~IndepADvar() {
740 #ifdef RAD_Const_WARN 741 inline operator ADVari&()
const {
742 ADVari *tcv = this->cv;
747 inline operator ADVari*()
const {
748 ADVari *tcv = this->cv;
753 #else //RAD_Const_WARN 756 #endif //RAD_Const_WARN 792 #define Ai const Base< ADVari >& 793 #define AI const Base< IndepADvar >& 801 #define T1(f) friend ADVari& f<> (AI); 803 #define F friend ADVari& 842 #ifdef HAVE_SACADO_CXX11 855 template<
typename Double>
class 856 ADvar:
public IndepADvar<Double> {
866 #if RAD_REINIT == 0 //{ 868 if (ConstADVari::cadc.fpval_implies_const)
871 #ifdef RAD_AUTO_AD_Const //{ 873 x->ADvari_padv(
this);
896 #ifdef RAD_AUTO_AD_Const
897 inline ADvar(IndepADVar &
x) {
898 this->cv =
x.cv ?
new ADVar1(
this,
x) : 0;
901 inline ADvar(ADVari &
x) { this->cv = &
x;
x.ADvari_padv(
this); }
902 inline ADvar& operator=(IndepADVar &
x) {
905 this->cv =
new ADVar1(
this,
x);
908 inline ADvar& operator=(ADVari &
x) {
911 this->cv =
new ADVar1(
this,
x);
915 friend ADvar& ADvar_operatoreq<>(ADvar*, const ADVari&); 918 inline ADvar(
const IndepADVar &
x) {
920 this->cv = (ADVari*)
x.cv;
922 inline ADvar(
const ADVari &
x) { this->cv = (ADVari*)&
x; }
923 inline ADvar& operator=(IndepADVar &
x) {
925 this->cv = (ADVari*)
x.cv;
return *
this;
927 inline ADvar& operator=(
const ADVari &
x) { this->cv = (ADVari*)&
x;
return *
this; }
929 ADvar(const IndepADVar &x) { 933 this->cv =
new ADVar1(
x.cv->Val, &this->cv->adc.One,
x.cv);
943 this->cv =
new ADVar1(
x.cv->Val, &this->cv->adc.One, (
ADVari*)
x.cv);
951 this->cv =
new ADVar1(
x.Val, &this->cv->adc.One, &
x);
958 ADvar& operator+=(
const ADVari&);
960 ADvar& operator-=(
const ADVari&);
962 ADvar& operator*=(
const ADVari&);
964 ADvar& operator/=(
const ADVari&);
968 {
return ConstADVari::cadc.fpval_implies_const; }
970 { ConstADVari::cadc.fpval_implies_const = newval; }
972 bool oldval = ConstADVari::cadc.fpval_implies_const;
973 ConstADVari::cadc.fpval_implies_const = newval;
977 inline static bool get_fpval_implies_const(
void) {
return true; }
978 inline static void set_fpval_implies_const(
bool newval) {}
979 inline static bool setget_fpval_implies_const(
bool newval) {
return newval; }
985 static inline void aval_reset() { ConstADVari::aval_reset(); }
993 template<
typename Double>
997 template<
typename Double>
1000 template<
typename Double>
1001 inline void AD_Const(
const IndepADvar<Double>&v) {}
1002 #endif //RAD_REINIT == 0 1004 template<
typename Double>
class 1005 ConstADvar:
public ADvar<Double> {
1022 void ConstADvar_ctr(
Double);
1032 #ifdef RAD_NO_CONST_UPDATE 1038 this->cv =
new ADVari(d);
1047 this->cv =
new ADVar1(
this,d);
1056 #ifdef RAD_ALLOW_WANTDERIV 1057 template<
typename Double>
class 1061 typedef ADvar<Double>
ADVar;
1062 typedef IndepADvar<Double> IndepADVar;
1063 typedef typename IndepADVar::ADVari ADVari;
1064 typedef ADvar1<Double> ADVar1;
1066 void ADvar_ndctr(
Double d) {
1067 ADVari *
x =
new ADVari(d);
1082 this->cv =
new ADVar1(
x.cv->Val, &this->cv->adc.One,
x.cv);
1089 this->cv =
new ADVar1(
x.Val, &this->cv->adc.One, &
x);
1095 #define ADvar_nd ADvar 1096 #endif //RAD_ALLOW_WANTDERIV 1098 template<
typename Double>
class 1099 ADvar1s:
public ADvar1<Double> {
1103 #if RAD_REINIT > 0 //{{ 1108 #ifdef RAD_AUTO_AD_Const 1111 ADVar1(val1,&
a,c1,v),
a(a1) {}
1113 #endif // }} RAD_REINIT > 0 1116 template<
typename Double>
class 1117 ADvar2:
public ADvari<Double> {
1122 #if RAD_REINIT > 0 //{{ 1124 const ADVari *Rcv,
const Double *Rc): ADVari(val1) {
1125 #ifdef RAD_ALLOW_WANTDERIV 1127 Ld = ADVari::adc.new_Derp(Lc,
this, Lcv);
1128 Rd = ADVari::adc.new_Derp(Rc,
this, Rcv);
1132 ADVari::adc.new_Derp(Lc, this, Lcv); 1133 ADVari::adc.new_Derp(Rc,
this, Rcv);
1134 #endif //RAD_ALLOW_WANTDERIV 1136 #else //}{ RAD_REINIT == 0 1141 dR.
next = DErp::LastDerp;
1143 DErp::LastDerp = &dL;
1150 #ifdef RAD_AUTO_AD_Const 1155 dR.
next = DErp::LastDerp;
1157 DErp::LastDerp = &dL;
1164 this->ADvari_padv(v);
1167 #endif // }} RAD_REINIT 1170 template<
typename Double>
class 1171 ADvar2q:
public ADvar2<Double> {
1176 #if RAD_REINIT > 0 //{{ 1179 #ifdef RAD_ALLOW_WANTDERIV 1181 Ld = ADVari::adc.new_Derp(&Lp,
this, Lcv);
1182 Rd = ADVari::adc.new_Derp(&Rp,
this, Rcv);
1186 ADVari::adc.new_Derp(&Lp, this, Lcv); 1187 ADVari::adc.new_Derp(&Rp,
this, Rcv);
1188 #endif //RAD_ALLOW_WANTDERIV 1190 #else //}{ RADINIT == 0 1194 this->dR.next = DErp::LastDerp;
1195 this->dL.next = &this->dR;
1196 DErp::LastDerp = &this->dL;
1201 this->dL.b = this->dR.b =
this;
1203 #ifdef RAD_AUTO_AD_Const 1206 const ADVari *Rcv,
const ADVar *v):
1207 ADVar2(val1),
a(Lp), b(Rp) {
1208 this->dR.next = DErp::LastDerp;
1209 this->dL.next = &this->dR;
1210 DErp::LastDerp = &this->dL;
1215 this->dL.b = this->dR.b =
this;
1217 this->ADvari_padv(v);
1220 #endif //}} RAD_REINIT > 0 1223 template<
typename Double>
class 1224 ADvarn:
public ADvari<Double> {
1229 #if RAD_REINIT > 0 // {{ 1231 #ifdef RAD_ALLOW_WANTDERIV 1233 for(
i = nd = 0;
i < n1;
i++)
1234 if (ADVari::adc.new_Derp(g+
i,
this,
x[
i].cv))
1239 for(
int i = 0;
i < n1;
i++)
1240 ADVari::adc.new_Derp(g+
i,
this,
x[
i].cv);
1241 #endif // RAD_ALLOW_WANTDERIV 1253 a1 =
a = (
Double*)ADVari::adc.Memalloc(n*
sizeof(*
a));
1254 d1 = Da = (
DErp*)ADVari::adc.Memalloc(n*
sizeof(
DErp));
1255 dlast = DErp::LastDerp;
1256 for(
i = 0;
i < n1;
i++, d1++) {
1264 DErp::LastDerp = dlast;
1266 #endif //}} RAD_REINIT > 0 1269 template<
typename Double>
1274 template<
typename Double>
1278 return L.
Val <
R.Val;
1280 template<
typename Double>
1281 inline int operator<(const Base< ADvari<Double> > &LL,
Double R) {
1285 template<
typename Double>
1286 inline int operator<(Double L, const Base< ADvari<Double> > &RR) {
1291 template<
typename Double>
1295 return L.
Val <=
R.Val;
1297 template<
typename Double>
1298 inline int operator<=(const Base< ADvari<Double> > &LL,
Double R) {
1302 template<
typename Double>
1303 inline int operator<=(Double L, const Base< ADvari<Double> > &RR) {
1308 template<
typename Double>
1312 return L.
Val ==
R.Val;
1314 template<
typename Double>
1319 template<
typename Double>
1325 template<
typename Double>
1329 return L.
Val !=
R.Val;
1331 template<
typename Double>
1336 template<
typename Double>
1342 template<
typename Double>
1346 return L.
Val >=
R.Val;
1348 template<
typename Double>
1353 template<
typename Double>
1359 template<
typename Double>
1363 return L.
Val >
R.Val;
1365 template<
typename Double>
1370 template<
typename Double>
1376 template<
typename Double>
1379 return Mbase + (Mleft -= len);
1380 return new_ADmemblock(len);
1382 #if RAD_REINIT > 0 //{{ 1384 template<
typename Double>
1387 template<
typename Double>
1391 template<
typename Double>
1393 a(*a1), b(b1),
c(c1) {}
1394 #else //}{ RAD_REINIT == 0 1396 template<
typename Double>
1402 template<
typename Double>
1409 template<
typename Double>
1414 #endif //}} RAD_REINIT 1428 #ifdef RAD_AUTO_AD_Const
1434 #ifndef RAD_DEBUG_gcgen1 1435 #define RAD_DEBUG_gcgen1 -1 1451 Mbase = (
char*)First->
memblk;
1452 Mleft =
sizeof(First->
memblk);
1454 #ifdef RAD_DEBUG_BLOCKKEEP 1455 rad_busy_blocks = 0;
1463 DMleft = nderps =
sizeof(DBusy->
memblk)/
sizeof(
DErp);
1472 for(mb = ADVari::adc.Busy; mb; mb = mb1) {
1476 for(mb = ADVari::adc.Free; mb; mb = mb1) {
1480 for(mb = ConstADVari::cadc.Busy; mb; mb = mb1) {
1484 ConstADVari::cadc.Busy = ADVari::adc.Busy = ADVari::adc.Free = 0;
1485 ConstADVari::cadc.Mleft = ADVari::adc.Mleft = 0;
1486 ConstADVari::cadc.Mbase = ADVari::adc.Mbase = 0;
1488 for(mb = ADVari::adc.DBusy; mb; mb = mb1) {
1492 for(mb = ADVari::adc.DFree; mb; mb = mb1) {
1496 ADVari::adc.DBusy = ADVari::adc.DFree = 0;
1497 ADVari::adc.DMleft = 0;
1498 ConstADVari::cadc.Mbase = ADVari::adc.Mbase = 0;
1500 ConstADVari::lastcad = 0;
1509 if (ConstADVari::cadc.Busy || ADVari::adc.Busy || ADVari::adc.Free
1511 || ADVari::adc.DBusy || ADVari::adc.DFree
1514 ADVari::adc.do_init();
1515 ConstADVari::cadc.do_init();
1518 template<
typename Double>
void*
1522 #ifdef RAD_AUTO_AD_Const 1525 #ifdef RAD_Const_WARN 1540 #ifdef RAD_DEBUG_BLOCKKEEP
1541 Mleft = rad_mleft_save;
1542 if (Mleft <
sizeof(First->
memblk))
1544 UninitType<Double>::utype,
1545 (
sizeof(First->
memblk) - Mleft)
1547 if ((mb = Busy->
next)) {
1549 for(;; mb = mb->
next) {
1551 UninitType<Double>::utype,
1558 rad_Oldcurmb = Busy;
1559 if (rad_busy_blocks >= RAD_DEBUG_BLOCKKEEP) {
1560 rad_busy_blocks = 0;
1564 for(mb = Busy; mb != mb0; mb = mb1) {
1571 Mbase = (
char*)First->
memblk;
1572 Mleft =
sizeof(First->
memblk);
1579 for(mb = Busy; mb != mb0; mb = mb1) {
1586 Mbase = (
char*)First->
memblk;
1587 Mleft =
sizeof(First->
memblk);
1589 #ifdef RAD_AUTO_AD_Const // { 1590 *ADVari::Last_ADvari = 0;
1591 ADVari::Last_ADvari = &ADVari::First_ADvari;
1592 a = ADVari::First_ADvari;
1597 #ifdef RAD_Const_WARN 1598 if ((
i =
a->opno) > 0)
1611 #endif // } RAD_AUTO_AD_Const 1612 #if RAD_REINIT > 0 //{ 1614 while((mb1 = mb->
next)) {
1625 while((vb = vb->ADvnext) != vb0)
1626 if ((tcv = ((IADv*)vb)->cv))
1628 #elif RAD_REINIT == 2 1633 return Mbase + (Mleft -= len);
1640 #ifdef RAD_DEBUG_BLOCKKEEP
1645 return (Mbase = (
char*)
x->memblk) +
1646 (Mleft =
sizeof(First->
memblk) - len);
1649 template<
typename Double>
void 1652 #if RAD_REINIT > 0 //{{ 1657 mb = ADVari::adc.DBusy;
1658 d = ((
DErp*)mb->
memblk) + ADVari::adc.DMleft;
1659 de = ((
DErp*)mb->
memblk) + ADVari::adc.nderps;
1663 if (!(mb = mb->
next))
1666 de = d + ADVari::adc.nderps;
1669 #else //}{ RAD_REINIT == 0 1673 for(d = DErp::LastDerp; d; d = d->
next)
1676 #endif //}} RAD_REINIT 1679 ADVari::adc.rad_need_reinit = 1;
1680 ADVari::adc.rad_mleft_save = ADVari::adc.Mleft;
1681 ADVari::adc.Mleft = 0;
1685 if (ADVari::gcgen_cur == ADVari::zap_gcgen1 && wantgrad) {
1687 if (!(fname = getenv(
"RAD_DEBUG_FILE")))
1688 fname =
"rad_debug.out";
1692 ADVari::debug_file = fopen(fname,
"w");
1693 ADVari::zap_gcgen1 = -1;
1696 #if RAD_REINIT > 0 //{{ 1697 if (ADVari::adc.DMleft < ADVari::adc.nderps && wantgrad) {
1698 mb = ADVari::adc.DBusy;
1699 d = ((
DErp*)mb->
memblk) + ADVari::adc.DMleft;
1700 de = ((
DErp*)mb->
memblk) + ADVari::adc.nderps;
1704 if (ADVari::debug_file) {
1705 for(; d < de; d++) {
1706 fprintf(ADVari::debug_file,
"%d\t%d\t%g + %g * %g",
1709 fprintf(ADVari::debug_file,
" = %g\n", d->
c->
aval);
1716 if (!(mb = mb->
next))
1719 de = d + ADVari::adc.nderps;
1722 #else //}{ RAD_REINIT == 0 1723 if ((d = DErp::LastDerp) && wantgrad) {
1726 if (ADVari::debug_file)
1728 fprintf(ADVari::debug_file,
"%d\t%d\t%g + %g * %g",
1731 fprintf(ADVari::debug_file,
" = %g\n", d->
c->
aval);
1732 }
while((d = d->
next));
1736 while((d = d->
next));
1739 if (ADVari::debug_file) {
1740 fclose(ADVari::debug_file);
1741 ADVari::debug_file = 0;
1744 #endif // }} RAD_REINIT 1746 if (ADVari::debug_file)
1747 fflush(ADVari::debug_file);
1748 ADVari::gcgen_cur++;
1749 ADVari::last_opno = 0;
1753 template<
typename Double>
void 1757 #if RAD_REINIT > 0 //{{ 1762 mb = ADVari::adc.DBusy;
1763 d = ((
DErp*)mb->
memblk) + ADVari::adc.DMleft;
1764 de = ((
DErp*)mb->
memblk) + ADVari::adc.nderps;
1768 if (!(mb = mb->
next))
1771 de = d + ADVari::adc.nderps;
1774 #else //}{ RAD_REINIT == 0 1778 for(d = DErp::LastDerp; d; d = d->
next)
1781 #endif //}} RAD_REINIT 1784 ADVari::adc.rad_need_reinit = 1;
1785 ADVari::adc.rad_mleft_save = ADVari::adc.Mleft;
1786 ADVari::adc.Mleft = 0;
1790 if (ADVari::gcgen_cur == ADVari::zap_gcgen1) {
1792 if (!(fname = getenv(
"RAD_DEBUG_FILE")))
1793 fname =
"rad_debug.out";
1797 ADVari::debug_file = fopen(fname,
"w");
1798 ADVari::zap_gcgen1 = -1;
1801 #if RAD_REINIT > 0 //{{ 1802 if (ADVari::adc.DMleft < ADVari::adc.nderps) {
1803 for(
i = 0;
i < n;
i++)
1805 mb = ADVari::adc.DBusy;
1806 d = ((
DErp*)mb->
memblk) + ADVari::adc.DMleft;
1807 de = ((
DErp*)mb->
memblk) + ADVari::adc.nderps;
1811 if (ADVari::debug_file) {
1812 for(; d < de; d++) {
1813 fprintf(ADVari::debug_file,
"%d\t%d\t%g + %g * %g",
1816 fprintf(ADVari::debug_file,
" = %g\n", d->
c->
aval);
1823 if (!(mb = mb->
next))
1826 de = d + ADVari::adc.nderps;
1829 #else //}{ RAD_REINIT == 0 1830 if ((d = DErp::LastDerp) != 0) {
1831 for(
i = 0;
i < n;
i++)
1834 if (ADVari::debug_file)
1836 fprintf(ADVari::debug_file,
"%d\t%d\t%g + %g * %g",
1839 fprintf(ADVari::debug_file,
" = %g\n", d->
c->
aval);
1840 }
while((d = d->
next));
1844 while((d = d->
next));
1847 if (ADVari::debug_file) {
1848 fclose(ADVari::debug_file);
1849 ADVari::debug_file = 0;
1852 #endif // }} RAD_REINIT 1854 if (ADVari::debug_file)
1855 fflush(ADVari::debug_file);
1856 ADVari::gcgen_cur++;
1857 ADVari::last_opno = 0;
1861 template<
typename Double>
void 1869 template<
typename Double>
1874 RAD_REINIT_2(Val = d; this->gen = this->IndepADvar_root.gen;)
1877 template<
typename Double>
1882 RAD_REINIT_2(Val =
i; this->gen = this->IndepADvar_root.gen;)
1885 template<
typename Double>
1890 RAD_REINIT_2(Val =
i; this->gen = this->IndepADvar_root.gen;)
1893 template<
typename Double>
1898 RAD_REINIT_2(Val =
i; this->gen = this->IndepADvar_root.gen;)
1901 template<
typename Double>
1906 RAD_REINIT_2(this->Val = 0.; this->gen = this->IndepADvar_root.gen;)
1909 template<
typename Double>
void 1914 RAD_REINIT_2(this->Val = d; this->gen = this->IndepADvar_root.gen;)
1917 template<
typename Double>
1924 x.cv->adc.new_Derp(&
x.adc.One,
y,
x.cv);
1929 RAD_REINIT_2(this->Val =
y->Val; this->gen = this->IndepADvar_root.gen;)
1932 template<
typename Double>
1939 x.cv->adc.new_Derp(&
x.cv->adc.One,
y,
x.cv);
1944 RAD_REINIT_2(this->Val =
y->Val; this->gen = this->IndepADvar_root.gen;)
1947 template<
typename Double>
1953 x.adc.new_Derp(&
x.adc.One,
y, &
x);
1958 RAD_REINIT_2(this->Val =
y->Val; this->gen = this->IndepADvar_root.gen;)
1961 template<
typename Double>
1967 ConstADVari *ncv =
new ConstADVari(v.
val());
1968 #ifdef RAD_AUTO_AD_Const 1975 template<
typename Double>
1979 #ifdef RAD_ALLOW_WANTDERIV 1981 if (this->gen != this->IndepADvar_root.gen) {
1983 this->gen = this->IndepADvar_root.gen;
1986 return this->wantderiv = cv->wantderiv = n;
1989 #endif // RAD_ALLOW_WANTDERIV 1992 template<
typename Double>
2002 #elif RAD_REINIT == 1 2007 #ifdef RAD_AUTO_AD_Const 2009 template<
typename Double>
2012 this->ADvari_padv(
x);
2016 template<
typename Double>
2019 this->ADvari_padv(
x);
2023 template<
typename Double>
2025 ADVari(
y.cv->Val), d((const
Double*)&ADcontext<
Double>::One, (ADVari*)this,
y.cv)
2027 this->ADvari_padv(
x);
2030 template<
typename Double>
2032 ADVari(
y.Val), d((const
Double*)&ADcontext<
Double>::One, this, &
y)
2034 this->ADvari_padv(
x);
2039 template<
typename Double>
2045 RAD_REINIT_2(This->Val =
x.Val; This->gen = This->IndepADvar_root.gen;)
2049 template<
typename Double>
2055 RAD_REINIT_2(This->Val =
x.Val; This->gen = This->IndepADvar_root.gen;)
2062 template<
typename Double>
2066 #ifdef RAD_AUTO_AD_Const 2069 this->cv =
new ADVari(
this,d);
2071 this->cv =
new ADVari(d);
2078 template<
typename Double>
2082 #ifdef RAD_AUTO_AD_Const 2085 this->cv =
new ADVari(
this,d);
2087 this->cv =
RAD_REINIT_0(ConstADVari::cadc.fpval_implies_const
2091 RAD_REINIT_2(this->Val = d; this->gen = this->IndepADvar_root.gen;)
2096 template<
typename Double>
2103 template<
typename Double>
2111 #ifdef RAD_AUTO_AD_Const 2112 #define RAD_ACA ,this 2117 template<
typename Double>
2128 template<
typename Double>
2135 template<
typename Double>
2146 template<
typename Double>
2153 template<
typename Double>
2161 template<
typename Double>
2172 template<
typename Double>
2179 template<
typename Double>
2190 template<
typename Double>
2197 template<
typename Double>
2205 template<
typename Double>
2216 template<
typename Double>
2223 template<
typename Double>
2234 template<
typename Double>
2241 template<
typename Double>
2246 Double Lv = L.
Val, Rv =
R.Val, pL = 1. / Rv, q = Lv/Rv;
2250 template<
typename Double>
2254 Double Lv = Lcv->
Val, Rv =
R.Val, pL = 1. / Rv, q = Lv/Rv;
2262 template<
typename Double>
2269 template<
typename Double>
2278 template<
typename Double>
2289 template<
typename Double>
2297 template<
typename Double>
2305 template<
typename Double>
2313 template<
typename Double>
2325 template<
typename Double>
2333 template<
typename Double>
2341 template<
typename Double>
2350 template<
typename Double>
2358 template<
typename Double>
2366 template<
typename Double>
2375 template<
typename Double>
2384 template<
typename Double>
2393 template<
typename Double>
2402 template<
typename Double>
2411 template<
typename Double>
2420 template<
typename Double>
2427 template<
typename Double>
2434 template<
typename Double>
2443 rcv->
d.a = &rcv->
Val;
2449 template<
typename Double>
2457 template<
typename Double>
2461 static double num = 1. /
std::log(10.);
2466 template<
typename Double>
2475 template<
typename Double>
2483 template<
typename Double>
2491 template<
typename Double>
2498 template<
typename Double>
2505 template<
typename Double>
2513 template<
typename Double>
2521 template<
typename Double>
2529 template<
typename Double>
2535 if ((t = v.
Val) < 0) {
2542 template<
typename Double>
2551 if ((t = v.
Val) < 0) {
2558 #ifdef HAVE_SACADO_CXX11 2559 template<
typename Double>
2561 cbrt(
const Base< ADvari<Double> > &vv) {
2562 const ADvari<Double>& v = vv.derived();
2564 return *(
new ADvar1s<Double>(t, 1.0/(3.0*t*t), &v));
2568 template<
typename Double>
2574 template<
typename Double>
2575 inline ADvari<Double>&
2580 template<
typename Double>
2586 template<
typename Double>
2592 template<
typename Double>
2598 template<
typename Double>
2604 template<
typename Double>
2610 template<
typename Double>
2611 inline ADvari<Double>&
2616 template<
typename Double>
2623 #define A (ADvari<Double>*) 2624 #ifdef RAD_Const_WARN 2625 #define C(x) (((x)->opno < 0) ? RAD_Const_Warn(x) : 0, *A x) 2629 #define T template<typename Double> inline 2630 #define F ADvari<Double>& 2631 #define Ai const Base< ADvari<Double> >& 2632 #define AI const Base< IndepADvar<Double> >& 2634 #define CAI(x,y) const IndepADvar<Double> & x = y.derived() 2635 #define CAi(x,y) const ADvari<Double> & x = y.derived() 2637 T r f(Ai LL, AI RR) { CAi(L,LL); CAI(R,RR); RAD_cvchk(R) return f(L, C(R.cv)); } \ 2638 T r f(AI LL, Ai RR) { CAI(L,LL); CAi(R,RR); RAD_cvchk(L) return f(C(L.cv), R); }\ 2639 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)); }\ 2640 T r f(AI LL, D R) { CAI(L,LL); RAD_cvchk(L) return f(C(L.cv), R); } \ 2641 T r f(D L, AI RR) { CAI(R,RR); RAD_cvchk(R) return f(L, C(R.cv)); } \ 2642 T r f(Ai L, Dtype R) { return f(L, (D)R); }\ 2643 T r f(AI L, Dtype R) { return f(L, (D)R); }\ 2644 T r f(Ai L, Ltype R) { return f(L, (D)R); }\ 2645 T r f(AI L, Ltype R) { return f(L, (D)R); }\ 2646 T r f(Ai L, Itype R) { return f(L, (D)R); }\ 2647 T r f(AI L, Itype R) { return f(L, (D)R); }\ 2648 T r f(Dtype L, Ai R) { return f((D)L, R); }\ 2649 T r f(Dtype L, AI R) {return f((D)L, R); }\ 2650 T r f(Ltype L, Ai R) { return f((D)L, R); }\ 2651 T r f(Ltype L, AI R) { return f((D)L, R); }\ 2652 T r f(Itype L, Ai R) { return f((D)L, R); }\ 2653 T r f(Itype L, AI R) { return f((D)L, R); } 2674 T F f(AI xx) { CAI(x,xx); RAD_cvchk(x) return f(C(x.cv)); } 2696 #ifdef HAVE_SACADO_CXX11 2731 #undef Allow_noderiv
static void Outvar_Gradcomp(ADVar &v)
ADvari< Double > & asinh(const Base< ADvari< Double > > &vv)
ADT_RAD ADvari< double > Ai
ADvari< Double > & pow(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
static ADcontext< Double > adc
ADvari< Double > & pow(const Base< ADvari< Double > > &LL, Double y)
ADvari< Double > & operator/(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
static void aval_reset(void)
ADvari< Double > & exp(const Base< ADvari< Double > > &vv)
const derived_type & derived() const
ADvar2q(Double val1, Double Lp, Double Rp, const ADVari *Lcv, const ADVari *Rcv)
ADvari< Double > & atan2(const Base< ADvari< Double > > &LL, Double y)
void AD_Const1(Double *, const IndepADvar< Double > &)
ADT_RAD IndepADvar< double > AI
static CADcontext< Double > cadc
ADvar & operator+=(const ADVari &)
Sacado::RadVec::ADvar< double > ADVar
ADvari< Double > & abs(const Base< ADvari< Double > > &vv)
ADvari< Double > & ADf2(Double f, Double gx, Double gy, const IndepADvar< Double > &x, const IndepADvar< Double > &y)
void AD_Const(const IndepADvar< Double > &)
ADvar & operator*=(const ADVari &)
void ConstADvar_ctr(Double)
static bool setget_fpval_implies_const(bool newval)
int operator>(Double L, const Base< ADvari< Double > > &RR)
ADvari< Double > & log(const Base< ADvari< Double > > &vv)
Double val(const ADvari< Double > &)
IndepADvar & operator=(IndepADvar &x)
Turn ADvar into a meta-function class usable with mpl::apply.
ADvari< Double > & atan2(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
ADvar1(Double val1, const ADVari *c1)
ConstADvar & operator=(ADVari &d)
Base class for Sacado types to control overload resolution.
IndepADvar< Double > & ADvar_operatoreq(IndepADvar< Double > *, const ADvari< Double > &)
ADvari< Double > & cosh(const Base< ADvari< Double > > &vv)
void _uninit_f2c(void *x, int type, long len)
ADvari< Double > & sqrt(const Base< ADvari< Double > > &vv)
ADvari< Double > & max(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
ADvar1s(Double val1, Double a1, const ADVari *c1)
static void AD_Const(const IndepADvar &)
IndepADvar() Allow_noderiv(
IndepADVar::ADVari ADVari
static void Gradcomp(int wantgrad)
ADvari< Double > & asin(const Base< ADvari< Double > > &vv)
static void Weighted_Gradcomp(size_t, ADVar **, Double *)
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
FloatingPoint< double > Double
ADvar2(Double val1, const ADVari *Lcv, const Double *Lc, const ADVari *Rcv, const Double *Rc)
static void set_fpval_implies_const(bool newval)
static int rad_need_reinit
ADvari< Double > & atanh(const Base< ADvari< Double > > &vv)
ConstADvar & operator=(Double d)
int RAD_Const_Warn(void *v)
ADvar & operator-=(const ADVari &)
ADvar< Double > & ADvar_operatoreq(ADvar< Double > *, const ADvari< Double > &)
ADvari< Double > & operator+(Double L, const Base< ADvari< Double > > &RR)
ADvar & operator=(const ADVari &x)
ADvar1(Double val1, const Double *a1, const ADVari *c1)
static void Weighted_Gradcomp(size_t n, ADVar **v, Double *w)
ScalarType< value_type >::type scalar_type
ADvari< Double > & acosh(const Base< ADvari< Double > > &vv)
static void Gradcomp(int wantgrad)
ADvari< Double > & min(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
ADvar & operator/=(const ADVari &)
ADvari< Double > & tan(const Base< ADvari< Double > > &vv)
ADvari< Double > & log10(const Base< ADvari< Double > > &vv)
ADvari< Double > & acos(const Base< ADvari< Double > > &vv)
ADvari< Double > & operator-(const Base< ADvari< Double > > &TT)
ADvari< Double > & tanh(const Base< ADvari< Double > > &vv)
ADvari< Double > & fabs(const Base< ADvari< Double > > &vv)
bool operator==(const Handle< T > &h1, const Handle< T > &h2)
Compare two handles.
ADmemblock< Double > ADMemblock
static void Outvar_Gradcomp(ADvar &v)
ADvari< Double > & operator*(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
IndepADvar_base(Allow_noderiv(int wd))
ADvari< Double > & operator+(const Base< ADvari< Double > > &TT)
ADvarn(Double val1, int n1, const IndepADVar *x, const Double *g)
static void Outvar_Gradcomp(ADVar &)
Allow_noderiv(mutable int wantderiv;) void *operator new(size_t len)
ADvari< Double > & sinh(const Base< ADvari< Double > > &vv)
int operator>=(Double L, const Base< ADvari< Double > > &RR)
IndepADvar< Double > IndepADVar
int operator!=(Double L, const Base< ADvari< Double > > &RR)
ADVari::IndepADVar IndepADVar
ADvari< Double > & sin(const Base< ADvari< Double > > &vv)
void * new_ADmemblock(size_t)
ADvari< Double > & cos(const Base< ADvari< Double > > &vv)
ADvari< Double > & ADf1(Double f, Double g, const IndepADvar< Double > &x)
static void Weighted_Gradcomp(size_t n, ADvar **v, Double *w)
IndepADvar< Double > IndepADVar
ADvari< Double > & ADfn(Double f, int n, const IndepADvar< Double > *x, const Double *g)
static bool get_fpval_implies_const(void)
ConstADvari< Double > ConstADVari
ADVar::ConstADVari ConstADVari
ADvari< Double > & atan(const Base< ADvari< Double > > &vv)
ADVar::IndepADVar IndepADVar