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);
252 static const Double
One, negOne;
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; }
295 Derp(
const ADVari *);
296 Derp(
const Double *,
const ADVari *);
297 Derp(
const Double *,
const ADVari *,
const ADVari *);
298 inline void *
operator new(
size_t len) {
return (
Derp*)ADVari::adc.Memalloc(len); }
302 #if RAD_REINIT > 0 //{ 303 template<
typename Double> Derp<Double>*
304 ADcontext<Double>::new_Derp(
const Double *
a,
const ADvari<Double> *b,
const ADvari<Double> *
c)
309 if (this->DMleft == 0) {
316 this->DMleft = nderps;
318 rv = &((DErp*)DBusy->
memblk)[--this->DMleft];
326 #define Ai const Base< ADvari<Double> >& 327 #define AI const Base< IndepADvar<Double> >& 328 #define T template<typename Double> 355 #define F ADvari<Double>& 390 #ifdef HAVE_SACADO_CXX11 405 template<
typename Double>ADvari<Double>&
ADf1(Double f, Double g,
const IndepADvar<Double> &x);
406 template<
typename Double>ADvari<Double>&
ADf2(Double f, Double gx, Double gy,
407 const IndepADvar<Double> &x,
const IndepADvar<Double> &y);
408 template<
typename Double>ADvari<Double>&
ADfn(Double f,
int n,
409 const IndepADvar<Double> *x,
const Double *g);
411 template<
typename Double> IndepADvar<Double>&
ADvar_operatoreq(IndepADvar<Double>*,
412 const ADvari<Double>&);
413 template<
typename Double> ADvar<Double>&
ADvar_operatoreq(ADvar<Double>*,
const ADvari<Double>&);
414 template<
typename Double>
void AD_Const(
const IndepADvar<Double>&);
415 template<
typename Double>
void AD_Const1(Double*,
const IndepADvar<Double>&);
416 template<
typename Double> ADvari<Double>&
ADf1(Double, Double,
const ADvari<Double>&);
417 template<
typename Double> ADvari<Double>&
ADf2(Double, Double, Double,
418 const ADvari<Double>&,
const ADvari<Double>&);
419 template<
typename Double> ADvari<Double>&
ADf2(Double, Double, Double,
420 const IndepADvar<Double>&,
const ADvari<Double>&);
421 template<
typename Double> ADvari<Double>&
ADf2(Double, Double, Double,
422 const ADvari<Double>&,
const IndepADvar<Double>&);
423 template<
typename Double> Double
val(
const ADvari<Double>&);
425 template<
typename Double>
class 426 ADvari :
public Base< ADvari<Double> > {
431 #ifdef RAD_AUTO_AD_Const 433 #ifdef RAD_Const_WARN 438 #endif //RAD_Const_WARN 441 static ADvari *First_ADvari, **Last_ADvari;
443 inline void ADvari_padv(
const IndepADVar *v) {
445 Last_ADvari = &this->Next;
448 #endif //RAD_AUTO_AD_Const 452 static int gcgen_cur, last_opno, zap_gcgen, zap_gcgen1, zap_opno;
453 static FILE *debug_file;
458 void *
operator new(
size_t len) {
461 rv->gcgen = gcgen_cur;
462 rv->opno = ++last_opno;
463 if (last_opno == zap_opno && gcgen_cur == zap_gcgen)
464 printf(
"Got to opno %d\n", last_opno);
470 void operator delete(
void*) {}
471 #ifdef RAD_ALLOW_WANTDERIV 472 inline ADvari(Double t): Val(t), aval(0.), wantderiv(1) {}
473 inline ADvari(): Val(0.), aval(0.), wantderiv(1) {}
474 inline void no_deriv() { wantderiv = 0; }
476 inline ADvari(Double t): Val(t), aval(0.) {}
479 #ifdef RAD_AUTO_AD_Const 481 friend class ADvar<Double>;
482 friend class ADvar1<Double>;
484 friend class ADvar2<Double>;
487 ADvari(
const IndepADVar *, Double);
488 ADvari(
const IndepADVar *, Double, Double);
492 #define Ai const Base< ADvari >& 494 #define T1(r,f) F r f <>(Ai); 528 #ifdef HAVE_SACADO_CXX11 544 friend ADvari& ADf1<>(Double
f, Double
g,
const ADvari &x);
545 friend ADvari& ADf2<>(Double
f, Double gx, Double gy,
const ADvari &x,
const ADvari &y);
546 friend ADvari& ADfn<>(Double
f,
int n,
const IndepADVar *x,
const Double *
g);
548 inline operator Double() {
return this->Val; }
549 inline operator Double()
const {
return this->Val; }
554 template<
typename Double>
class 560 ADvar1(Double val1,
const Double *a1,
const ADVari *c1): ADVari(val1) {
561 #ifdef RAD_ALLOW_WANTDERIV 562 if (!ADVari::adc.new_Derp(a1,
this,c1))
565 ADVari::adc.new_Derp(a1,
this,c1);
568 #else // RAD_REINIT == 0 572 ADVari(val1), d(a1,this,c1) {}
573 #ifdef RAD_AUTO_AD_Const 574 typedef typename ADVari::IndepADVar IndepADVar;
576 ADvar1(
const IndepADVar*,
const IndepADVar&);
577 ADvar1(
const IndepADVar*,
const ADVari&);
578 ADvar1(
const Double val1,
const Double *a1,
const ADVari *c1,
const ADVar *v):
579 ADVari(val1), d(a1,this,c1) {
581 this->ADvari_padv(v);
584 #endif // RAD_REININT > 0 588 template<
typename Double>
class 589 ConstADvari:
public ADvari<Double> {
603 inline void *
operator new(
size_t len) {
return ADVari::adc.Memalloc(len); }
604 inline ConstADvari(Double t): ADVari(t) {}
606 static void aval_reset(
void);
609 template<
typename Double>
class 615 #elif RAD_REINIT == 2 616 typedef unsigned long ADGenType;
617 mutable ADGenType gen;
623 template<
typename Double>
class 633 IndepADvar_root.ADvprev = (this->ADvprev = IndepADvar_root.ADvprev)->ADvnext =
this;
634 this->ADvnext = &IndepADvar_root;
645 (this->ADvnext->ADvprev = this->ADvprev)->ADvnext = this->ADvnext;
658 #if RAD_REINIT > 0 //{ 659 template<
typename Double> IndepADvar_base0<Double>
660 IndepADvar_base<Double>::IndepADvar_root(0.);
662 template<
typename Double>
void 663 IndepADvar_base<Double>::Move_to_end() {
664 if (
this != IndepADvar_root.ADvprev) {
665 (this->ADvnext->ADvprev = this->ADvprev)->ADvnext = this->ADvnext;
666 IndepADvar_root.ADvprev = (this->ADvprev = IndepADvar_root.ADvprev)->ADvnext =
this;
667 this->ADvnext = &IndepADvar_root;
670 #elif RAD_REINIT == 2 671 template<
typename Double>
typename IndepADvar_base0<Double>::ADGenType
672 IndepADvar_base<Double>::gen0(1);
676 template<
typename Double>
class 677 IndepADvar:
protected IndepADvar_base<Double>,
public Base< IndepADvar<Double> > {
679 static void AD_Const(
const IndepADvar&);
685 #ifdef RAD_AUTO_AD_Const 690 #elif defined(RAD_EQ_ALIAS) 695 #endif //RAD_AUTO_AD_Const 712 #ifdef RAD_ALLOW_WANTDERIV 713 inline int Wantderiv() {
return this->wantderiv; }
716 IndepADvar(
Ttype,
int);
717 IndepADvar(Double,
int);
722 #ifdef RAD_AUTO_AD_Const 724 inline IndepADvar() { cv = 0; }
725 inline ~IndepADvar() {
739 #ifdef RAD_Const_WARN 740 inline operator ADVari&()
const {
741 ADVari *tcv = this->cv;
746 inline operator ADVari*()
const {
747 ADVari *tcv = this->cv;
752 #else //RAD_Const_WARN 755 #endif //RAD_Const_WARN 757 inline Double
val()
const {
764 inline Double
adj()
const {
770 friend void AD_Const1<>(Double*,
const IndepADvar&);
772 friend ADVari& ADf1<>(Double, Double,
const IndepADvar&);
774 friend ADVari& ADf2<>(Double, Double, Double,
const ADVari&,
const IndepADvar&);
775 friend ADVari& ADf2<>(Double, Double, Double,
const IndepADvar&,
const ADVari&);
791 #define Ai const Base< ADVari >& 792 #define AI const Base< IndepADvar >& 800 #define T1(f) friend ADVari& f<> (AI); 802 #define F friend ADVari& 841 #ifdef HAVE_SACADO_CXX11 854 template<
typename Double>
class 855 ADvar:
public IndepADvar<Double> {
865 #if RAD_REINIT == 0 //{ 867 if (ConstADVari::cadc.fpval_implies_const)
870 #ifdef RAD_AUTO_AD_Const //{ 872 x->ADvari_padv(
this);
890 inline ADvar(
double i) { ADvar_ctr(Double(i)); }
891 inline ADvar(
int i) { ADvar_ctr(Double(i)); }
892 inline ADvar(
long i) { ADvar_ctr(Double(i)); }
895 #ifdef RAD_AUTO_AD_Const
896 inline ADvar(IndepADVar &x) {
897 this->cv = x.cv ?
new ADVar1(
this, x) : 0;
900 inline ADvar(ADVari &x) { this->cv = &x; x.ADvari_padv(
this); }
901 inline ADvar& operator=(IndepADVar &x) {
904 this->cv =
new ADVar1(
this,x);
907 inline ADvar& operator=(ADVari &x) {
910 this->cv =
new ADVar1(
this, x);
914 friend ADvar& ADvar_operatoreq<>(ADvar*, const ADVari&); 917 inline ADvar(
const IndepADVar &x) {
919 this->cv = (ADVari*)x.cv;
921 inline ADvar(
const ADVari &x) { this->cv = (ADVari*)&x; }
922 inline ADvar& operator=(IndepADVar &x) {
924 this->cv = (ADVari*)x.cv;
return *
this;
926 inline ADvar& operator=(
const ADVari &x) { this->cv = (ADVari*)&x;
return *
this; }
928 ADvar(const IndepADVar &x) { 932 this->cv =
new ADVar1(x.cv->Val, &this->cv->adc.One, x.cv);
950 this->cv =
new ADVar1(x.
Val, &this->cv->adc.One, &x);
956 ADvar& operator=(Double);
957 ADvar& operator+=(
const ADVari&);
958 ADvar& operator+=(Double);
959 ADvar& operator-=(
const ADVari&);
960 ADvar& operator-=(Double);
961 ADvar& operator*=(
const ADVari&);
962 ADvar& operator*=(Double);
963 ADvar& operator/=(
const ADVari&);
964 ADvar& operator/=(Double);
967 {
return ConstADVari::cadc.fpval_implies_const; }
969 { ConstADVari::cadc.fpval_implies_const = newval; }
971 bool oldval = ConstADVari::cadc.fpval_implies_const;
972 ConstADVari::cadc.fpval_implies_const = newval;
976 inline static bool get_fpval_implies_const(
void) {
return true; }
977 inline static void set_fpval_implies_const(
bool newval) {}
978 inline static bool setget_fpval_implies_const(
bool newval) {
return newval; }
984 static inline void aval_reset() { ConstADVari::aval_reset(); }
992 template<
typename Double>
996 template<
typename Double>
999 template<
typename Double>
1000 inline void AD_Const(
const IndepADvar<Double>&v) {}
1001 #endif //RAD_REINIT == 0 1003 template<
typename Double>
class 1004 ConstADvar:
public ADvar<Double> {
1021 void ConstADvar_ctr(Double);
1031 #ifdef RAD_NO_CONST_UPDATE 1037 this->cv =
new ADVari(d);
1046 this->cv =
new ADVar1(
this,d);
1055 #ifdef RAD_ALLOW_WANTDERIV 1056 template<
typename Double>
class 1060 typedef ADvar<Double>
ADVar;
1061 typedef IndepADvar<Double> IndepADVar;
1062 typedef typename IndepADVar::ADVari ADVari;
1063 typedef ADvar1<Double> ADVar1;
1065 void ADvar_ndctr(Double d) {
1066 ADVari *x =
new ADVari(d);
1074 inline ADvar_nd(
double i):
ADVar((void*)0,0) { ADvar_ndctr(Double(i)); }
1075 inline ADvar_nd(
int i):
ADVar((void*)0,0) { ADvar_ndctr(Double(i)); }
1076 inline ADvar_nd(
long i):
ADVar((void*)0,0) { ADvar_ndctr(Double(i)); }
1077 inline ADvar_nd(Double d):
ADVar((void*)0,0) { ADvar_ndctr(d); }
1081 this->cv =
new ADVar1(x.cv->Val, &this->cv->adc.One, x.cv);
1088 this->cv =
new ADVar1(x.Val, &this->cv->adc.One, &x);
1094 #define ADvar_nd ADvar 1095 #endif //RAD_ALLOW_WANTDERIV 1097 template<
typename Double>
class 1098 ADvar1s:
public ADvar1<Double> {
1102 #if RAD_REINIT > 0 //{{ 1107 #ifdef RAD_AUTO_AD_Const 1109 ADvar1s(Double val1, Double a1,
const ADVari *c1,
const ADVar *v):
1110 ADVar1(val1,&
a,c1,v),
a(a1) {}
1112 #endif // }} RAD_REINIT > 0 1115 template<
typename Double>
class 1116 ADvar2:
public ADvari<Double> {
1121 #if RAD_REINIT > 0 //{{ 1122 ADvar2(Double val1,
const ADVari *Lcv,
const Double *Lc,
1123 const ADVari *Rcv,
const Double *Rc): ADVari(val1) {
1124 #ifdef RAD_ALLOW_WANTDERIV 1126 Ld = ADVari::adc.new_Derp(Lc,
this, Lcv);
1127 Rd = ADVari::adc.new_Derp(Rc,
this, Rcv);
1131 ADVari::adc.new_Derp(Lc, this, Lcv); 1132 ADVari::adc.new_Derp(Rc,
this, Rcv);
1133 #endif //RAD_ALLOW_WANTDERIV 1135 #else //}{ RAD_REINIT == 0 1138 const ADVari *Rcv,
const Double *Rc):
1140 dR.
next = DErp::LastDerp;
1142 DErp::LastDerp = &dL;
1149 #ifdef RAD_AUTO_AD_Const 1151 ADvar2(Double val1,
const ADVari *Lcv,
const Double *Lc,
1152 const ADVari *Rcv,
const Double *Rc,
ADVar *v):
1154 dR.
next = DErp::LastDerp;
1156 DErp::LastDerp = &dL;
1163 this->ADvari_padv(v);
1166 #endif // }} RAD_REINIT 1169 template<
typename Double>
class 1170 ADvar2q:
public ADvar2<Double> {
1175 #if RAD_REINIT > 0 //{{ 1178 #ifdef RAD_ALLOW_WANTDERIV 1180 Ld = ADVari::adc.new_Derp(&Lp,
this, Lcv);
1181 Rd = ADVari::adc.new_Derp(&Rp,
this, Rcv);
1185 ADVari::adc.new_Derp(&Lp, this, Lcv); 1186 ADVari::adc.new_Derp(&Rp,
this, Rcv);
1187 #endif //RAD_ALLOW_WANTDERIV 1189 #else //}{ RADINIT == 0 1193 this->dR.next = DErp::LastDerp;
1194 this->dL.next = &this->dR;
1195 DErp::LastDerp = &this->dL;
1200 this->dL.b = this->dR.b =
this;
1202 #ifdef RAD_AUTO_AD_Const 1204 ADvar2q(Double val1, Double Lp, Double Rp,
const ADVari *Lcv,
1205 const ADVari *Rcv,
const ADVar *v):
1206 ADVar2(val1),
a(Lp), b(Rp) {
1207 this->dR.next = DErp::LastDerp;
1208 this->dL.next = &this->dR;
1209 DErp::LastDerp = &this->dL;
1214 this->dL.b = this->dR.b =
this;
1216 this->ADvari_padv(v);
1219 #endif //}} RAD_REINIT > 0 1222 template<
typename Double>
class 1223 ADvarn:
public ADvari<Double> {
1228 #if RAD_REINIT > 0 // {{ 1230 #ifdef RAD_ALLOW_WANTDERIV 1232 for(i = nd = 0; i < n1; i++)
1233 if (ADVari::adc.new_Derp(
g+i,
this, x[i].cv))
1238 for(
int i = 0; i < n1; i++)
1239 ADVari::adc.new_Derp(
g+i,
this, x[i].cv);
1240 #endif // RAD_ALLOW_WANTDERIV 1252 a1 =
a = (Double*)ADVari::adc.Memalloc(
n*
sizeof(*
a));
1253 d1 = Da = (
DErp*)ADVari::adc.Memalloc(
n*
sizeof(
DErp));
1254 dlast = DErp::LastDerp;
1255 for(i = 0; i < n1; i++, d1++) {
1263 DErp::LastDerp = dlast;
1265 #endif //}} RAD_REINIT > 0 1268 template<
typename Double>
1273 template<
typename Double>
1277 return L.
Val <
R.Val;
1279 template<
typename Double>
1280 inline int operator<(const Base< ADvari<Double> > &LL, Double
R) {
1284 template<
typename Double>
1285 inline int operator<(Double L, const Base< ADvari<Double> > &RR) {
1290 template<
typename Double>
1294 return L.
Val <=
R.Val;
1296 template<
typename Double>
1297 inline int operator<=(const Base< ADvari<Double> > &LL, Double
R) {
1301 template<
typename Double>
1302 inline int operator<=(Double L, const Base< ADvari<Double> > &RR) {
1307 template<
typename Double>
1311 return L.
Val ==
R.Val;
1313 template<
typename Double>
1318 template<
typename Double>
1324 template<
typename Double>
1328 return L.
Val !=
R.Val;
1330 template<
typename Double>
1335 template<
typename Double>
1341 template<
typename Double>
1345 return L.
Val >=
R.Val;
1347 template<
typename Double>
1352 template<
typename Double>
1358 template<
typename Double>
1362 return L.
Val >
R.Val;
1364 template<
typename Double>
1369 template<
typename Double>
1375 template<
typename Double>
1378 return Mbase + (Mleft -= len);
1379 return new_ADmemblock(len);
1381 #if RAD_REINIT > 0 //{{ 1383 template<
typename Double>
1386 template<
typename Double>
1390 template<
typename Double>
1392 a(*a1), b(b1),
c(c1) {}
1393 #else //}{ RAD_REINIT == 0 1395 template<
typename Double>
1401 template<
typename Double>
1408 template<
typename Double>
1413 #endif //}} RAD_REINIT 1427 #ifdef RAD_AUTO_AD_Const
1433 #ifndef RAD_DEBUG_gcgen1 1434 #define RAD_DEBUG_gcgen1 -1 1450 Mbase = (
char*)First->
memblk;
1451 Mleft =
sizeof(First->
memblk);
1453 #ifdef RAD_DEBUG_BLOCKKEEP 1454 rad_busy_blocks = 0;
1462 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;
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();
1519 template<
typename Double>
void*
1523 #ifdef RAD_AUTO_AD_Const 1526 #ifdef RAD_Const_WARN 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 if ((
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++)
1804 V[i]->cv->
aval = w[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++)
1832 V[i]->cv->
aval = w[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>
1881 cv =
new ADVari(Double(i));
1882 RAD_REINIT_2(Val = i; this->gen = this->IndepADvar_root.gen;)
1885 template<
typename Double>
1889 cv =
new ADVari(Double(i));
1890 RAD_REINIT_2(Val = i; this->gen = this->IndepADvar_root.gen;)
1893 template<
typename Double>
1897 cv =
new ADVari(Double(i));
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>
2273 Double recip = 1. /
R.Val;
2274 Double q = L * recip;
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>
2346 Double x = L.
Val, y =
R.Val, t = x*x + y*y;
2350 template<
typename Double>
2354 Double y =
R.Val, t = x*x + y*y;
2358 template<
typename Double>
2362 Double x = L.
Val, t = x*x + y*y;
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)
bool operator!=(const Allocator< T > &a_t, const Allocator< U > &a_u)
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
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
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