35 #ifndef SACADO_TRAD2_H
36 #define SACADO_TRAD2_H
44 #if defined(RAD_DEBUG_BLOCKKEEP) && !defined(HAVE_SACADO_UNINIT)
45 #undef RAD_DEBUG_BLOCKKEEP
48 #ifdef RAD_Const_WARN // ==> RAD_AUTO_AD_Const and RAD_DEBUG
49 #ifndef RAD_AUTO_AD_Const
50 #define RAD_AUTO_AD_Const
57 #endif // RAD_Const_WARN
64 #ifndef RAD_AUTO_AD_Const
65 #ifdef RAD_DEBUG_BLOCKKEEP
74 #ifndef RAD_NO_USING_STDCC
95 #ifdef RAD_AUTO_AD_Const
96 #undef RAD_DEBUG_BLOCKKEEP
98 #ifdef RAD_DEBUG_BLOCKKEEP
99 #if !(RAD_DEBUG_BLOCKKEEP > 0)
100 #undef RAD_DEBUG_BLOCKKEEP
102 extern "C" void _uninit_f2c(
void *x,
int type,
long len);
104 template <
typename T>
105 struct UninitType {};
108 struct UninitType<float> {
109 static const int utype = 4;
113 struct UninitType<double> {
114 static const int utype = 5;
117 template <
typename T>
118 struct UninitType<
std::complex<T> > {
119 static const int utype = UninitType<T>::utype + 2;
128 template<
typename T>
class
141 #define Dtype typename DoubleAvoid<Double>::dtype
142 #define Ttype typename DoubleAvoid<Double>::ttype
147 template<
typename Double>
class ADvar;
157 template<
typename Double>
class Derp;
159 template<
typename Double>
struct
167 template<
typename Double>
class
171 enum { Gulp = 1021 };
177 template<
typename Double>
class
196 double First0[(
sizeof(
ADMemblock) +
sizeof(
double) - 1) /
sizeof(
double)];
197 double First1[(
sizeof(
ADVari_block) +
sizeof(
double) - 1) /
sizeof(
double)];
198 void *new_ADmemblock(
size_t);
199 void new_ADvari_block();
204 #ifdef RAD_DEBUG_BLOCKKEEP
209 static const Double
One, negOne;
211 void *Memalloc(
size_t len);
212 static void Gradcomp(
int wantgrad);
213 static void Hvprod(
int,
ADVar**, Double*, Double*);
214 static void aval_reset(
void);
215 static void Weighted_Gradcomp(
int,
ADVar**, Double*);
218 if (Ainext >= Ailimit)
224 template<
typename Double>
class
233 template<
typename Double>
class
244 Derp(
const ADVari *);
245 Derp(
const Double *,
const ADVari *);
246 Derp(
const Double *,
const ADVari *,
const ADVari *);
253 #define Ai const ADvari<Double>&
254 #define AI const IndepADvar<Double>&
255 #define T template<typename Double>
282 #define F ADvari<Double>&
331 #define SNS Sacado::Rad2
355 template<
typename Double>ADvari<Double>&
ADf1(Double f, Double g,
const IndepADvar<Double> &x);
356 template<
typename Double>ADvari<Double>&
ADf2(Double f, Double gx, Double gy,
357 const IndepADvar<Double> &x,
const IndepADvar<Double> &y);
358 template<
typename Double>ADvari<Double>&
ADfn(Double f,
int n,
359 const IndepADvar<Double> *x,
const Double *g);
361 template<
typename Double> IndepADvar<Double>&
ADvar_operatoreq(IndepADvar<Double>*,
362 const ADvari<Double>&);
363 template<
typename Double> ADvar<Double>&
ADvar_operatoreq(ADvar<Double>*,
const ADvari<Double>&);
364 template<
typename Double>
void AD_Const(
const IndepADvar<Double>&);
365 template<
typename Double>
void AD_Const1(Double*,
const IndepADvar<Double>&);
366 template<
typename Double> ADvari<Double>&
ADf1(Double, Double,
const ADvari<Double>&);
367 template<
typename Double> ADvari<Double>&
ADf2(Double, Double, Double,
368 const ADvari<Double>&,
const ADvari<Double>&);
369 template<
typename Double> ADvari<Double>&
ADf2(Double, Double, Double,
370 const IndepADvar<Double>&,
const ADvari<Double>&);
371 template<
typename Double> ADvari<Double>&
ADf2(Double, Double, Double,
372 const ADvari<Double>&,
const IndepADvar<Double>&);
373 template<
typename Double> Double
val(
const ADvari<Double>&);
389 template<
typename Double> ADvari<Double>&
390 ADf1(Double f, Double g, Double h,
const ADvari<Double> &x);
392 template<
typename Double> ADvari<Double>&
393 ADf2(Double f, Double gx, Double gy, Double hxx,
394 Double hxy, Double hyy,
const ADvari<Double> &x,
const ADvari<Double> &y);
396 template<
typename Double> ADvari<Double>&
397 ADfn(Double f,
int n,
const IndepADvar<Double> *x,
const Double *g,
const Double *h);
399 template<
typename Double>
class
405 #ifdef RAD_AUTO_AD_Const
407 #ifdef RAD_Const_WARN
412 #endif //RAD_Const_WARN
415 static ADvari *First_ADvari, **Last_ADvari;
417 inline void ADvari_padv(
const IndepADVar *v) {
419 Last_ADvari = &this->Next;
422 #endif //RAD_AUTO_AD_Const
426 static int gcgen_cur, last_opno, zap_gcgen, zap_gcgen1, zap_opno;
427 static FILE *debug_file;
436 void *
operator new(
size_t len) {
439 rv->gcgen = gcgen_cur;
440 rv->opno = ++last_opno;
441 if (last_opno == zap_opno && gcgen_cur == zap_gcgen)
448 void operator delete(
void*) {}
450 opclass(oc), Val(t), aval(0.), dO(0.)
453 opclass(oc), Val(t), aval(ta), dO(0.)
456 inline ADvari(): Val(0.), aval(0.), dO(0.) {}
459 #ifdef RAD_AUTO_AD_Const
461 friend class ADvar<Double>;
462 friend class ADvar1<Double>;
464 friend class ADvar2<Double>;
473 #define Ai const ADvari&
474 #define T1(r,f) F r f <>(Ai);
477 F r f <>(Ttype,Ai); \
478 F r f <>(Ai,Ttype); \
479 F r f <>(double,Ai); \
480 F r f <>(Ai,double); \
526 friend ADvari& ADf1<>(Double f, Double g, Double h,
const ADvari &x);
527 friend ADvari& ADf2<>(Double f, Double gx, Double gy, Double hxx,
528 Double hxy, Double hyy,
const ADvari &x,
const ADvari &y);
529 friend ADvari& ADfn<>(Double f,
int n,
const IndepADVar *x,
530 const Double *g,
const Double *h);
532 inline operator Double() {
return this->Val; }
533 inline operator Double()
const {
return this->Val; }
536 template<
typename Double>
class
537 ADvar1:
public ADvari<Double> {
542 ADVari(oc,val1), d(a1,this,c1) {}
543 #ifdef RAD_AUTO_AD_Const
544 typedef typename ADVari::IndepADVar IndepADVar;
546 ADvar1(
const IndepADVar*,
const IndepADVar&);
547 ADvar1(
const IndepADVar*,
const ADVari&);
549 const ADVari *c1,
const ADVar *v):
550 ADVari(oc,val1), d(a1,this,c1) {
552 this->ADvari_padv(v);
558 template<
typename Double>
class
559 ConstADvari:
public ADvari<Double> {
570 static void aval_reset(
void);
574 template<
typename Double>
class
579 #ifdef RAD_AUTO_AD_Const
583 #elif defined(RAD_EQ_ALIAS)
588 #endif //RAD_AUTO_AD_Const
606 #ifdef RAD_AUTO_AD_Const
610 inline IndepADvar() { cv = 0; }
611 inline ~IndepADvar() {
625 #ifdef RAD_Const_WARN
626 inline operator ADVari&()
const {
627 ADVari *tcv = this->cv;
632 inline operator ADVari*()
const {
633 ADVari *tcv = this->cv;
638 #else //RAD_Const_WARN
639 inline operator ADVari&()
const {
return *this->cv; }
640 inline operator ADVari*()
const {
return this->cv; }
641 #endif //RAD_Const_WARN
646 friend void AD_Const1<>(Double*,
const IndepADvar&);
648 friend ADVari& ADf1<>(Double, Double,
const IndepADvar&);
650 friend ADVari& ADf2<>(Double, Double, Double,
const ADVari&,
const IndepADvar&);
651 friend ADVari& ADf2<>(Double, Double, Double,
const IndepADvar&,
const ADVari&);
657 static inline void Hvprod(
int n,
ADVar **vp, Double *v, Double *hv)
667 #define Ai const ADVari&
668 #define AI const IndepADvar&
682 #define T1(f) friend ADVari& f<> (AI);
684 #define F friend ADVari&
733 template<
typename Double>
class
734 ADvar:
public IndepADvar<Double> {
745 if (ConstADVari::cadc.fpval_implies_const)
748 #ifdef RAD_AUTO_AD_Const
750 x->ADvari_padv(
this);
762 ADvar(
double i) { ADvar_ctr(Double(i)); }
763 ADvar(
int i) { ADvar_ctr(Double(i)); }
764 ADvar(
long i) { ADvar_ctr(Double(i)); }
766 #ifdef RAD_AUTO_AD_Const
767 ADvar(IndepADVar &x) {
768 this->cv = x.cv ?
new ADVar1(
this, x) : 0;
770 ADvar(ADVari &x) :IndepADVar(&x) { x.ADvari_padv((IndepADVar*)
this);}
771 inline ADvar& operator=(IndepADVar &x) {
774 this->cv =
new ADVar1(
this,x);
777 inline ADvar& operator=(ADVari &x) {
780 this->cv =
new ADVar1(
this, x);
784 friend ADvar& ADvar_operatoreq<>(ADvar*, const ADVari&);
787 inline ADvar(
const IndepADVar &x) { this->cv = (ADVari*)x.cv; }
788 inline ADvar(
const ADVari &x) { this->cv = (ADVari*)&x; }
789 inline ADvar& operator=(IndepADVar &x) { this->cv = (ADVari*)x.cv;
return *
this; }
790 inline ADvar& operator=(
const ADVari &x) { this->cv = (ADVari*)&x;
return *
this; }
792 friend ADvar& ADvar_operatoreq<>(ADvar*, const ADVari&);
801 ADvar& operator=(Double);
802 ADvar& operator+=(
const ADVari&);
803 ADvar& operator+=(Double);
804 ADvar& operator-=(
const ADVari&);
805 ADvar& operator-=(Double);
806 ADvar& operator*=(
const ADVari&);
807 ADvar& operator*=(Double);
808 ADvar& operator/=(
const ADVari&);
809 ADvar& operator/=(Double);
811 {
return ConstADVari::cadc.fpval_implies_const; }
813 { ConstADVari::cadc.fpval_implies_const = newval; }
815 bool oldval = ConstADVari::cadc.fpval_implies_const;
816 ConstADVari::cadc.fpval_implies_const = newval;
823 static inline void aval_reset() { ConstADVari::aval_reset(); }
828 template<
typename Double>
832 template<
typename Double>
835 template<
typename Double>
class
836 ConstADvar:
public ADvar<Double> {
852 void ConstADvar_ctr(Double);
862 #ifdef RAD_NO_CONST_UPDATE
870 template<
typename Double>
class
871 ADvar1s:
public ADvar1<Double> {
878 #ifdef RAD_AUTO_AD_Const
881 ADvar1s(Double val1, Double a1,
const ADVari *c1,
const ADVar *v):
886 template<
typename Double>
class
887 ADvar1g:
public ADvar1<Double> {
897 template<
typename Double>
class
898 ADvar2:
public ADvari<Double> {
904 const ADVari *Rcv,
const Double *Rc):
906 dR.
next = DErp::LastDerp;
908 DErp::LastDerp = &dL;
915 #ifdef RAD_AUTO_AD_Const
918 const ADVari *Rcv,
const Double *Rc,
ADVar *v):
920 dR.
next = DErp::LastDerp;
922 DErp::LastDerp = &dL;
929 this->ADvari_padv(v);
934 template<
typename Double>
class
935 ADvar2q:
public ADvar2<Double> {
944 ADvar2q(Double val1, Double Lp, Double Rp, Double LR, Double R2,
947 pL(Lp), pR(Rp), pLR(LR), pR2(R2) {}
948 #ifdef RAD_AUTO_AD_Const
950 ADvar2q(Double val1, Double Lp, Double Rp, Double LR, Double R2,
951 const ADVari *Lcv,
const ADVari *Rcv,
const ADVar *v):
952 ADVar2(
Hv_quotLR,val1,Lcv,&pL,Rcv,&pR,v),
953 pL(Lp), pR(Rp), pLR(LR), pR2(R2) {}
957 template<
typename Double>
class
958 ADvar2g:
public ADvar2<Double> {
967 ADvar2g(Double val1, Double Lp, Double Rp, Double L2, Double LR, Double R2,
970 pL(Lp), pR(Rp), pL2(L2), pLR(LR), pR2(R2) { }
973 template<
typename Double>
class
974 ADvarn:
public ADvari<Double> {
990 a1 = G = (Double*)ADVari::adc.Memalloc(n1*
sizeof(*
g));
991 d1 =
D = (
DErp*)ADVari::adc.Memalloc(n1*
sizeof(
DErp));
992 dlast = DErp::LastDerp;
993 for(i = 0; i < n1; i++, d1++) {
1001 DErp::LastDerp = dlast;
1002 nh = (n1*(n1+1)) >> 1;
1003 a1 = H = (
double*)ADVari::adc.Memalloc(nh *
sizeof(*h));
1004 for(i = 0; i < nh; i++)
1009 template<
typename Double>
1012 template<
typename Double>
1014 template<
typename Double>
1016 template<
typename Double>
1019 template<
typename Double>
1021 template<
typename Double>
1023 template<
typename Double>
1026 template<
typename Double>
1028 template<
typename Double>
1030 template<
typename Double>
1033 template<
typename Double>
1035 template<
typename Double>
1037 template<
typename Double>
1040 template<
typename Double>
1042 template<
typename Double>
1044 template<
typename Double>
1047 template<
typename Double>
1049 template<
typename Double>
1051 template<
typename Double>
1054 template<
typename Double>
1057 return Mbase + (Mleft -= len);
1058 return new_ADmemblock(len);
1061 template<
typename Double>
1067 template<
typename Double>
1073 template<
typename Double>
1089 #ifdef RAD_AUTO_AD_Const
1095 #ifndef RAD_DEBUG_gcgen1
1096 #define RAD_DEBUG_gcgen1 -1
1107 template<
typename Double>
1116 Mbase = (
char*)First->memblk;
1117 Mleft =
sizeof(First->memblk);
1122 fb->
limit = Ailimit = fb->
pADvari + ADVari_block::Gulp;
1124 #ifdef RAD_DEBUG_BLOCKKEEP
1125 rad_busy_blocks = 0;
1131 template<
typename Double>
void*
1136 #ifdef RAD_AUTO_AD_Const
1139 #ifdef RAD_Const_WARN
1148 Aibusy = b = AiFirst;
1151 Ailimit = b->
limit = (Ainext = b->
pADvari) + ADVari_block::Gulp;
1152 #ifdef RAD_DEBUG_BLOCKKEEP
1153 Mleft = rad_mleft_save;
1154 if (Mleft <
sizeof(First->memblk))
1156 UninitType<Double>::utype,
1157 (
sizeof(First->memblk) - Mleft)
1159 if ((mb = Busy->
next)) {
1160 if (!(mb0 = rad_Oldcurmb))
1162 for(;; mb = mb->
next) {
1164 UninitType<Double>::utype,
1165 sizeof(First->memblk)
1171 rad_Oldcurmb = Busy;
1172 if (rad_busy_blocks >= RAD_DEBUG_BLOCKKEEP) {
1173 rad_busy_blocks = 0;
1177 for(mb = Busy; mb != mb0; mb = mb1) {
1184 Mbase = (
char*)First->
memblk;
1185 Mleft =
sizeof(First->memblk);
1192 for(mb = Busy; mb != mb0; mb = mb1) {
1199 Mbase = (
char*)First->
memblk;
1200 Mleft =
sizeof(First->memblk);
1201 #ifdef RAD_AUTO_AD_Const
1203 ADVari::adc.rad_need_reinit &= ~2;
1204 *ADVari::Last_ADvari = 0;
1205 ADVari::Last_ADvari = &ADVari::First_ADvari;
1206 if ((anext = ADVari::First_ADvari)) {
1207 while((
a = anext)) {
1210 #ifdef RAD_Const_WARN
1211 if ((i =
a->opno) > 0)
1227 return Mbase + (Mleft -= len);
1234 #ifdef RAD_DEBUG_BLOCKKEEP
1239 return (Mbase = (
char*)x->
memblk) +
1240 (Mleft =
sizeof(First->memblk) - len);
1243 template<
typename Double>
void
1248 ob->
limit = Ailimit;
1253 Aibusy = Aibusy->
next = nb;
1254 nb->
limit = Ailimit = (Ainext = nb->
pADvari) + ADVari_block::Gulp;
1260 template<
typename Double>
void
1266 for(d = DErp::LastDerp; d; d = d->
next)
1270 ADVari::adc.rad_need_reinit = 1;
1271 ADVari::adc.rad_mleft_save = ADVari::adc.Mleft;
1272 ADVari::adc.Mleft = 0;
1275 if (ADVari::gcgen_cur == ADVari::zap_gcgen1) {
1277 if (!(fname = getenv(
"RAD_DEBUG_FILE")))
1278 fname =
"rad_debug.out";
1282 ADVari::debug_file = fopen(fname,
"w");
1283 ADVari::zap_gcgen1 = -1;
1286 if ((d = DErp::LastDerp) != 0) {
1287 #ifdef RAD_AUTO_AD_Const
1288 ADVari::adc.rad_need_reinit |= 2;
1293 if (ADVari::debug_file)
1295 fprintf(ADVari::debug_file,
"%d\t%d\t%g + %g * %g",
1296 d->
c->opno, d->
b->opno, d->
c->
aval,
1299 fprintf(ADVari::debug_file,
" = %g\n", d->
c->
aval);
1300 fflush(ADVari::debug_file);
1301 }
while((d = d->
next));
1305 while((d = d->
next));
1310 if (ADVari::debug_file) {
1311 fclose(ADVari::debug_file);
1312 ADVari::debug_file = 0;
1316 ADVari::gcgen_cur++;
1317 ADVari::last_opno = 0;
1321 template<
typename Double>
void
1326 #ifdef RAD_Const_WARN
1330 #ifdef RAD_AUTO_AD_Const
1336 for(d = DErp::LastDerp; d; d = d->
next)
1340 ADVari::adc.rad_need_reinit = 1;
1341 ADVari::adc.rad_mleft_save = ADVari::adc.Mleft;
1342 ADVari::adc.Mleft = 0;
1345 if (ADVari::gcgen_cur == ADVari::zap_gcgen1) {
1347 if (!(fname = getenv(
"RAD_DEBUG_FILE")))
1348 fname =
"rad_debug.out";
1352 ADVari::debug_file = fopen(fname,
"w");
1353 ADVari::zap_gcgen1 = -1;
1356 if ((d = DErp::LastDerp) != 0) {
1357 for(i = 0; i <
n; i++)
1358 V[i]->cv->
aval = w[i];
1360 if (ADVari::debug_file)
1362 fprintf(ADVari::debug_file,
"%d\t%d\t%g + %g * %g",
1365 fprintf(ADVari::debug_file,
" = %g\n", d->
c->
aval);
1366 fflush(ADVari::debug_file);
1367 }
while((d = d->
next));
1371 while((d = d->
next));
1374 if (ADVari::debug_file) {
1375 fclose(ADVari::debug_file);
1376 ADVari::debug_file = 0;
1379 #ifdef RAD_AUTO_AD_Const
1380 *ADVari::Last_ADvari = 0;
1381 ADVari::Last_ADvari = &ADVari::First_ADvari;
1382 if ((anext = ADVari::First_ADvari) && !(ADVari::adc.
rad_need_reinit & 2)) {
1383 ADVari::adc.rad_need_reinit = 3;
1384 while((
a = anext)) {
1387 #ifdef RAD_Const_WARN
1388 if ((i =
a->opno) > 0)
1403 ADVari::gcgen_cur++;
1404 ADVari::last_opno = 0;
1408 template<
typename Double>
1416 template<
typename Double>
1424 template<
typename Double>
1432 template<
typename Double>
1440 template<
typename Double>
1447 template<
typename Double>
void
1454 template<
typename Double>
1462 template<
typename Double>
1470 template<
typename Double>
1478 template<
typename Double>
1484 ConstADVari *ncv =
new ConstADVari(v.
val());
1485 #ifdef RAD_AUTO_AD_Const
1492 template<
typename Double>
1503 #ifdef RAD_AUTO_AD_Const
1505 template<
typename Double>
1509 this->ADvari_padv(x);
1512 template<
typename Double>
1516 this->ADvari_padv(x);
1519 template<
typename Double>
1521 ADVari(
Hv_copy, y.cv->Val), d((const Double*)&ADcontext<Double>::One, (ADVari*)this, y.cv)
1523 this->ADvari_padv(x);
1526 template<
typename Double>
1528 ADVari(
Hv_copy, y.Val), d((const Double*)&ADcontext<Double>::One, this, &y)
1530 this->ADvari_padv(x);
1535 template<
typename Double>
1541 template<
typename Double>
1546 template<
typename Double>
1550 #ifdef RAD_AUTO_AD_Const
1561 template<
typename Double>
1565 #ifdef RAD_AUTO_AD_Const
1571 this->cv = ConstADVari::cadc.fpval_implies_const
1578 template<
typename Double>
1584 template<
typename Double>
1590 #ifdef RAD_AUTO_AD_Const
1591 #define RAD_ACA ,this
1596 template<
typename Double>
1605 template<
typename Double>
1611 template<
typename Double>
1619 template<
typename Double>
1625 template<
typename Double>
1631 template<
typename Double>
1640 template<
typename Double>
1646 template<
typename Double>
1654 template<
typename Double>
1660 template<
typename Double>
1666 template<
typename Double>
1675 template<
typename Double>
1681 template<
typename Double>
1689 template<
typename Double>
1695 template<
typename Double>
1698 Double Lv = L.
Val, Rv =
R.Val, pL = 1. / Rv, q = Lv/Rv, qpL = q*pL;
1702 template<
typename Double>
1706 Double Lv = Lcv->
Val, Rv =
R.Val, pL = 1. / Rv, q = Lv/Rv, qpL = q*pL;
1711 template<
typename Double>
1717 template<
typename Double>
1720 Double recip = 1. /
R.Val;
1721 Double q = L * recip;
1722 Double d1 = -q*recip;
1726 template<
typename Double>
1734 template<
typename Double>
1737 Double t = v.
Val, t1 = 1. - t*t, d1 = -1./
std::sqrt(t1);
1741 template<
typename Double>
1744 Double d1, t, t1, t2;
1751 template<
typename Double>
1760 template<
typename Double>
1763 Double d1, t, t1, t2, td;
1773 template<
typename Double>
1776 Double t = v.
Val, d1 = 1./(1. + t*t);
1780 template<
typename Double>
1783 Double t = v.
Val, d1 = 1./(1. - t*t);
1787 template<
typename Double>
1790 Double R2, t, t2, x, x2, y, y2;
1801 template<
typename Double>
1804 Double t, x2, y, y2;
1812 template<
typename Double>
1815 Double t, x, x2, y2;
1823 template<
typename Double>
1830 template<
typename Double>
1838 template<
typename Double>
1846 template<
typename Double>
1853 template<
typename Double>
1861 template<
typename Double>
1869 template<
typename Double>
1876 template<
typename Double>
1883 template<
typename Double>
1890 template<
typename Double>
1893 Double x = v.
Val, d1 = 1. / x;
1897 template<
typename Double>
1900 static double num = 1. /
std::log(10.);
1908 template<
typename Double>
1911 Double
dx, dy, t, x, xlog, xym1, y;
1922 template<
typename Double>
1925 Double dy, t, xlog, y;
1933 template<
typename Double>
1943 template<
typename Double>
1950 template<
typename Double>
1957 template<
typename Double>
1961 Double d1 = 0.5 / t;
1965 template<
typename Double>
1975 template<
typename Double>
1985 template<
typename Double>
1990 if ((t = v.
Val) < 0) {
1997 template<
typename Double>
2004 if ((t = v.
Val) < 0) {
2011 template<
typename Double>
2017 template<
typename Double>
2018 inline ADvari<Double>&
2023 template<
typename Double>
2025 ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
2030 template<
typename Double>
2032 ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
2037 template<
typename Double>
2039 ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
2044 template<
typename Double>
2046 ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
2051 template<
typename Double>
2057 template<
typename Double>
2058 inline ADvari<Double>&
2063 template<
typename Double>
2067 ADVari *
a, *aL, *aR, **ap, **ape;
2070 Double aO, adO, *
g, *h, *h0, t, tL, tR;
2072 for(i = 0; i <
n; i++) {
2075 a->aO =
a->adO = 0.;
2077 ADVari::adc.Aibusy->limit = ADVari::adc.Ainext;
2079 for(b0 = 0, b = ADVari::adc.AiFirst; b; b0 = b, b = b->
next) {
2084 a->aO =
a->adO = 0.;
2085 switch(
a->opclass) {
2121 for(i = 0; i < m; i++)
2122 t +=
g[i] * d[i].
c->dO;
2132 for(b = b0; b; b = b->
prev) {
2139 switch(
a->opclass) {
2189 aL->
adO += adO * tL;
2194 aL->
aO += aO * (tL = aR->
Val) + adO*aR->
dO;
2195 aR->
aO += aO * (tR = aL->
Val) + adO*aL->
dO;
2196 aL->
adO += adO * tL;
2197 aR->
adO += adO * tR;
2217 for(i = 0; i < m; i++) {
2219 aL->
adO += adO * (t =
g[i]);
2222 for(h = h0, j = 0; j <= i; j++)
2223 d[j].
c->aO += t * *h++;
2225 for(k = j; j < m; j++)
2226 d[j].
c->
aO += t * *(h += k++);
2233 for(i = 0; i <
n; i++) {
2240 template<
typename Double>
2247 #define A (ADvari<Double>*)
2248 #ifdef RAD_Const_WARN
2249 #define C(x) (((x)->opno < 0) ? RAD_Const_Warn(x) : 0, *A x)
2253 #define T template<typename Double> inline
2254 #define F ADvari<Double>&
2255 #define Ai const ADvari<Double>&
2256 #define AI const IndepADvar<Double>&
2259 T r f(Ai L, AI R) { return f(L, C(R.cv)); }\
2260 T r f(AI L, Ai R) { return f(C(L.cv), R); }\
2261 T r f(AI L, AI R) { return f(C(L.cv), C(R.cv)); }\
2262 T r f(AI L, D R) { return f(C(L.cv), R); }\
2263 T r f(Ai L, Dtype R) { return f(L, (D)R); }\
2264 T r f(AI L, Dtype R) { return f(C(L.cv), (D)R); }\
2265 T r f(Ai L, long R) { return f(L, (D)R); }\
2266 T r f(AI L, long R) { return f(C(L.cv), (D)R); }\
2267 T r f(Ai L, int R) { return f(L, (D)R); }\
2268 T r f(AI L, int R) { return f(C(L.cv), (D)R); }\
2269 T r f(D L, AI R) { return f(L, C(R.cv)); }\
2270 T r f(Dtype L, Ai R) { return f((D)L, R); }\
2271 T r f(Dtype L, AI R) { return f((D)L, C(R.cv)); }\
2272 T r f(long L, Ai R) { return f((D)L, R); }\
2273 T r f(long L, AI R) { return f((D)L, C(R.cv)); }\
2274 T r f(int L, Ai R) { return f((D)L, R); }\
2275 T r f(int L, AI R) { return f((D)L, C(R.cv)); }
2296 T F f(AI x) { return f(C(x.cv)); }