Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Sacado_tradvec.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Sacado Package
5// Copyright (2006) Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// This library is free software; you can redistribute it and/or modify
11// it under the terms of the GNU Lesser General Public License as
12// published by the Free Software Foundation; either version 2.1 of the
13// License, or (at your option) any later version.
14//
15// This library is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18// Lesser General Public License for more details.
19//
20// You should have received a copy of the GNU Lesser General Public
21// License along with this library; if not, write to the Free Software
22// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23// USA
24// Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
25// (etphipp@sandia.gov).
26//
27// ***********************************************************************
28// @HEADER
29
30// TRAD package (Templated Reverse Automatic Differentiation) --
31// a package specialized for function and gradient evaluations.
32// Written in 2004 and 2005 by David M. Gay at Sandia National Labs,
33// Albuquerque, NM.
34
35#ifndef SACADO_TRADVEC_H
36#define SACADO_TRADVEC_H
37
38#include "Sacado_ConfigDefs.h"
40
41#include <stddef.h>
42#include <Sacado_cmath.hpp>
43
44#if defined(RAD_DEBUG_BLOCKKEEP) && !defined(HAVE_SACADO_UNINIT)
45#undef RAD_DEBUG_BLOCKKEEP
46#endif
47
48#ifdef RAD_Const_WARN //{ ==> RAD_AUTO_AD_Const and RAD_DEBUG
49#ifndef RAD_AUTO_AD_Const
50#define RAD_AUTO_AD_Const
51#endif
52#ifndef RAD_DEBUG
53#define RAD_DEBUG
54#endif
55extern "C" int RAD_Const_Warn(const void*);// outside any namespace for
56 // ease in setting breakpoints
57#endif //} RAD_Const_WARN
58
59#ifdef RAD_DEBUG
60#include <cstdio>
61#include <stdlib.h>
62#endif
63
64#ifndef RAD_AUTO_AD_Const
65#ifdef RAD_DEBUG_BLOCKKEEP
66#include <complex> // must come before namespace Sacado
67#endif
68#endif
69
70namespace Sacado {
71namespace RadVec {
72
73#ifdef RAD_AUTO_AD_Const
74#undef RAD_DEBUG_BLOCKKEEP
75#define Padvinit ,padv(0)
76#else
77#define Padvinit /*nothing*/
78#ifdef RAD_DEBUG_BLOCKKEEP
79#if !(RAD_DEBUG_BLOCKKEEP > 0)
80#undef RAD_DEBUG_BLOCKKEEP
81#else
82extern "C" void _uninit_f2c(void *x, int type, long len);
83
84template <typename T>
85struct UninitType {};
86
87template <>
88struct UninitType<float> {
89 static const int utype = 4;
90};
91
92template <>
93struct UninitType<double> {
94 static const int utype = 5;
95};
96
97template <typename T>
98struct UninitType< std::complex<T> > {
99 static const int utype = UninitType<T>::utype + 2;
100};
101
102#endif /*RAD_DEBUG_BLOCKKEEP > 0*/
103#endif /*RAD_DEBUG_BLOCKKEEP*/
104#endif /*RAD_AUTO_AD_Const*/
105
107
108 template<typename T> class
110 public:
111 typedef double dtype;
112 typedef T ttype;
113 };
114 template<> class
116 public:
119 };
120
121#define Dtype typename DoubleAvoid<Double>::dtype
122#define Ttype typename DoubleAvoid<Double>::ttype
123
124 template<typename Double> class IndepADvar;
125 template<typename Double> class ConstADvar;
126 template<typename Double> class ConstADvari;
127 template<typename Double> class ADvar;
128 template<typename Double> class ADvari;
129 template<typename Double> class ADvar1;
130 template<typename Double> class ADvar1s;
131 template<typename Double> class ADvar2;
132 template<typename Double> class ADvar2q;
133 template<typename Double> class ADvarn;
134 template<typename Double> class Derp;
135
136 template<typename Double> struct
137ADmemblock { // We get memory in ADmemblock chunks and never give it back,
138 // but reuse it once computations start anew after call(s) on
139 // ADcontext::Gradcomp() or ADcontext::Weighted_Gradcomp().
141 Double memblk[1000];
142 };
143
144 template<typename Double> class
145ADcontext {
150
151 ADMemblock *Busy, *First, *Free, *Oldbusy;
152 char *Mbase;
153 size_t Mleft;
155 size_t ncur, nmax, rad_mleft_save;
156#ifdef RAD_DEBUG_BLOCKKEEP
157 int rad_busy_blocks;
158 ADMemblock *rad_Oldcurmb;
159#endif
160 void *new_ADmemblock(size_t);
161 void derp_init(size_t);
162 public:
163 static const Double One, negOne;
164 ADcontext();
165 void *Memalloc(size_t len);
166 static void Gradcomp();
167 static void aval_reset(void);
168 static void Weighted_Gradcomp(size_t, ADVar**, Double*);
169 static void Weighted_GradcompVec(size_t, size_t*, ADVar***, Double**);
170 static void Outvar_Gradcomp(ADVar&);
171 };
172
173 template<typename Double> class
174CADcontext: public ADcontext<Double> { // for possibly constant ADvar values
175 protected:
176 bool fpval_implies_const;
177 public:
178 friend class ADvar<Double>;
179 CADcontext(): ADcontext<Double>() { fpval_implies_const = false; }
180 };
181
182 template<typename Double> class
183Derp { // one derivative-propagation operation
184 public:
185 friend class ADvarn<Double>;
187 static Derp *LastDerp;
189 const Double *a;
190 const ADVari *b;
191 const ADVari *c;
192 Derp(){};
193 Derp(const ADVari *);
194 Derp(const Double *, const ADVari *);
195 Derp(const Double *, const ADVari *, const ADVari *);
196 inline void *operator new(size_t len) { return (Derp*)ADVari::adc.Memalloc(len); }
197 /* c->aval += a * b->aval; */
198 };
199
200// Now we use #define to overcome bad design in the C++ templating system
201
202#define Ai const ADvari<Double>&
203#define AI const IndepADvar<Double>&
204#define T template<typename Double>
205#define D Double
206#define T1(f) \
207T F f (AI); \
208T F f (Ai);
209#define T2(r,f) \
210 T r f(Ai,Ai); \
211 T r f(Ai,D); \
212 T r f(Ai,Dtype); \
213 T r f(Ai,long); \
214 T r f(Ai,int); \
215 T r f(D,Ai); \
216 T r f(Dtype,Ai); \
217 T r f(long,Ai); \
218 T r f(int,Ai); \
219 T r f(AI,D); \
220 T r f(AI,Dtype); \
221 T r f(AI,long); \
222 T r f(AI,int); \
223 T r f(D,AI); \
224 T r f(Dtype,AI); \
225 T r f(long,AI); \
226 T r f(int,AI); \
227 T r f(Ai,AI);\
228 T r f(AI,Ai);\
229 T r f(AI,AI);
230
231#define F ADvari<Double>&
232T2(F, operator+)
233T2(F, operator-)
234T2(F, operator*)
235T2(F, operator/)
236T2(F, atan2)
237T2(F, pow)
238T2(F, max)
239T2(F, min)
240T2(int, operator<)
241T2(int, operator<=)
242T2(int, operator==)
243T2(int, operator!=)
244T2(int, operator>=)
245T2(int, operator>)
246T1(operator+)
247T1(operator-)
248T1(abs)
249T1(acos)
250T1(acosh)
251T1(asin)
252T1(asinh)
253T1(atan)
254T1(atanh)
255T1(cos)
256T1(cosh)
257T1(exp)
258T1(log)
259T1(log10)
260T1(sin)
261T1(sinh)
262T1(sqrt)
263T1(tan)
264T1(tanh)
265T1(fabs)
266
267T F copy(AI);
268T F copy(Ai);
269
270#undef F
271#undef T2
272#undef T1
273#undef D
274#undef T
275#undef AI
276#undef Ai
277
278 template<typename Double>ADvari<Double>& ADf1(Double f, Double g, const IndepADvar<Double> &x);
279 template<typename Double>ADvari<Double>& ADf2(Double f, Double gx, Double gy,
280 const IndepADvar<Double> &x, const IndepADvar<Double> &y);
281 template<typename Double>ADvari<Double>& ADfn(Double f, int n,
282 const IndepADvar<Double> *x, const Double *g);
283
284 template<typename Double> IndepADvar<Double>& ADvar_operatoreq(IndepADvar<Double>*,
285 const ADvari<Double>&);
286 template<typename Double> ADvar<Double>& ADvar_operatoreq(ADvar<Double>*, const ADvari<Double>&);
287 template<typename Double> void AD_Const(const IndepADvar<Double>&);
288 template<typename Double> void AD_Const1(Double*, const IndepADvar<Double>&);
289 template<typename Double> ADvari<Double>& ADf1(Double, Double, const ADvari<Double>&);
290 template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
291 const ADvari<Double>&, const ADvari<Double>&);
292 template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
293 const IndepADvar<Double>&, const ADvari<Double>&);
294 template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
295 const ADvari<Double>&, const IndepADvar<Double>&);
296 template<typename Double> Double val(const ADvari<Double>&);
297
298 template<typename Double> class
299ADvari { // implementation of an ADvar
300 public:
301 typedef Double value_type;
304#ifdef RAD_AUTO_AD_Const
305 friend class IndepADvar<Double>;
306#ifdef RAD_Const_WARN
307 mutable const IndepADVar *padv;
308#else
309 protected:
310 mutable const IndepADVar *padv;
311#endif //RAD_Const_WARN
312 public:
313 inline void ADvari_padv(const IndepADVar *v) {
314 this->padv = v;
315 }
316#endif //RAD_AUTO_AD_Const
317#ifdef RAD_DEBUG
318 int gcgen;
319 int opno;
320 static int gcgen_cur, last_opno, zap_gcgen, zap_gcgen1, zap_opno;
321 static FILE *debug_file;
322#endif
323 Double Val; // result of this operation
324 mutable Double *aval; // adjoint -- partial of final result w.r.t. this Val
325 static ADvari *First_ADvari, **Last_ADvari;
327 private:
328 inline void Linkup() {
329 *Last_ADvari = this;
330 Last_ADvari = &this->Next;
331 }
332 public:
333 void *operator new(size_t len) {
334#ifdef RAD_DEBUG
335 ADvari *rv = (ADvari*)ADvari::adc.Memalloc(len);
336 rv->gcgen = gcgen_cur;
337 rv->opno = ++last_opno;
338 if (last_opno == zap_opno && gcgen_cur == zap_gcgen)
339 printf("Got to opno %d\n", last_opno);
340 return rv;
341#else
342 return ADvari::adc.Memalloc(len);
343#endif
344 }
345 void operator delete(void*) {} /*Should never be called.*/
346 inline ADvari(Double t): Val(t), aval(0) Padvinit { Linkup(); }
347 inline ADvari(): Val(0.), aval(0) Padvinit { Linkup(); }
348#ifdef RAD_AUTO_AD_Const
349 friend class ADcontext<Double>;
350 friend class ADvar<Double>;
351 friend class ADvar1<Double>;
352 friend class ADvar1s<Double>;
353 friend class ADvar2<Double>;
354 friend class ADvar2q<Double>;
355 friend class ConstADvar<Double>;
356 ADvari(const IndepADVar *, Double);
357#endif
358#define F friend
359#define R ADvari&
360#define Ai const ADvari&
361#define T1(r,f) F r f <>(Ai);
362#define T2(r,f) \
363F r f <>(Ai,Ai); \
364F r f <>(Ttype,Ai); \
365F r f <>(Ai,Ttype); \
366F r f <>(double,Ai); \
367F r f <>(Ai,double); \
368F r f <>(long,Ai); \
369F r f <>(Ai,long); \
370F r f <>(int,Ai); \
371F r f <>(Ai,int);
372 T1(R,operator+)
373 T2(R,operator+)
374 T1(R,operator-)
375 T2(R,operator-)
376 T2(R,operator*)
377 T2(R,operator/)
378 T1(R,abs)
379 T1(R,acos)
380 T1(R,acosh)
381 T1(R,asin)
382 T1(R,asinh)
383 T1(R,atan)
384 T1(R,atanh)
385 T2(R,atan2)
386 T2(R,max)
387 T2(R,min)
388 T1(R,cos)
389 T1(R,cosh)
390 T1(R,exp)
391 T1(R,log)
392 T1(R,log10)
393 T2(R,pow)
394 T1(R,sin)
395 T1(R,sinh)
396 T1(R,sqrt)
397 T1(R,tan)
398 T1(R,tanh)
399 T1(R,fabs)
400 T1(R,copy)
401 T2(int,operator<)
402 T2(int,operator<=)
403 T2(int,operator==)
404 T2(int,operator!=)
405 T2(int,operator>=)
406 T2(int,operator>)
407#undef T2
408#undef T1
409#undef Ai
410#undef R
411#undef F
412
413 friend ADvari& ADf1<>(Double f, Double g, const ADvari &x);
414 friend ADvari& ADf2<>(Double f, Double gx, Double gy, const ADvari &x, const ADvari &y);
415 friend ADvari& ADfn<>(Double f, int n, const IndepADVar *x, const Double *g);
416
417 inline operator Double() { return this->Val; }
418 inline operator Double() const { return this->Val; }
419
421 };
422
423 template<typename Double> class
424ADvar1: public ADvari<Double> { // simplest unary ops
425 public:
427 ADvar1(Double val1): ADVari(val1) {}
429 ADvar1(Double val1, const ADVari *c1): ADVari(val1), d(c1) {}
430 ADvar1(Double val1, const Double *a1, const ADVari *c1):
431 ADVari(val1), d(a1,this,c1) {}
432#ifdef RAD_AUTO_AD_Const
433 typedef typename ADVari::IndepADVar IndepADVar;
434 typedef ADvar<Double> ADVar;
435 ADvar1(const IndepADVar*, const IndepADVar&);
436 ADvar1(const IndepADVar*, const ADVari&);
437 ADvar1(const Double val1, const Double *a1, const ADVari *c1, const ADVar *v):
438 ADVari(val1), d(a1,this,c1) {
439 c1->padv = 0;
440 this->ADvari_padv(v);
441 }
442#endif
443 };
444
445
446 template<typename Double> class
448 private:
450 ConstADvari() {}; // prevent construction without value (?)
451 static ConstADvari *lastcad;
452 public:
453 friend class ADcontext<Double>;
457 inline void *operator new(size_t len) { return ConstADvari::cadc.Memalloc(len); }
458 inline ConstADvari(Double t): ADVari(t) { prevcad = lastcad; lastcad = this; }
459 static void aval_reset(void);
460 };
461
462 template<typename Double> class
463IndepADvar { // an independent ADvar
464 protected:
465 static void AD_Const(const IndepADvar&);
467 private:
469 /* private to prevent assignment */
470#ifdef RAD_AUTO_AD_Const
471 if (cv)
472 cv->padv = 0;
473 cv = new ADvar1<Double>(this,x);
474 return *this;
475#elif defined(RAD_EQ_ALIAS)
476 this->cv = x.cv;
477 return *this;
478#else
479 return ADvar_operatoreq(this,*x.cv);
480#endif //RAD_AUTO_AD_Const
481 }
482 public:
483 int Wantderiv(int);
484 typedef Double value_type;
485 friend class ADvar<Double>;
486 friend class ADcontext<Double>;
487 friend class ADvar1<Double>;
488 friend class ADvarn<Double>;
492 IndepADvar(double);
493 IndepADvar(int);
494 IndepADvar(long);
495 IndepADvar& operator= (Double);
496 inline int Wantderiv() { return 1; }
497#ifdef RAD_AUTO_AD_Const //{
498 inline IndepADvar(const IndepADvar &x) { cv = x.cv ? new ADvar1<Double>(this, x) : 0; };
499 inline IndepADvar() { cv = 0; }
500 inline ~IndepADvar() {
501 if (cv)
502 cv->padv = 0;
503 }
504#else
505 inline IndepADvar() {
506#ifndef RAD_EQ_ALIAS
507 cv = 0;
508#endif
509 }
510 inline ~IndepADvar() {}
511 friend IndepADvar& ADvar_operatoreq<>(IndepADvar*, const ADVari&);
512#endif //}
513
514#ifdef RAD_Const_WARN
515 inline operator ADVari&() const {
516 ADVari *tcv = this->cv;
517 if (tcv->opno < 0)
518 RAD_Const_Warn(tcv);
519 return *tcv;
520 }
521 inline operator ADVari*() const {
522 ADVari *tcv = this->cv;
523 if (tcv->opno < 0)
524 RAD_Const_Warn(tcv);
525 return tcv;
526 }
527#else //RAD_Const_WARN
528 inline operator ADVari&() const { return *this->cv; }
529 inline operator ADVari*() const { return this->cv; }
530#endif //RAD_Const_WARN
531
532 inline Double val() const {
533 return cv->Val;
534 }
535 inline Double adj() const {
536 return *cv->aval;
537 }
538 inline Double adj(int n) const {
539 return cv->aval[n];
540 }
541
542 friend void AD_Const1<>(Double*, const IndepADvar&);
543
544 friend ADVari& ADf1<>(Double, Double, const IndepADvar&);
545 friend ADVari& ADf2<>(Double, Double, Double, const IndepADvar&, const IndepADvar&);
546 friend ADVari& ADf2<>(Double, Double, Double, const ADVari&, const IndepADvar&);
547 friend ADVari& ADf2<>(Double, Double, Double, const IndepADvar&, const ADVari&);
548
549 static inline void Gradcomp() { ADcontext<Double>::Gradcomp(); }
550 static inline void aval_reset() { ConstADvari<Double>::aval_reset(); }
551 static inline void Weighted_Gradcomp(size_t n, ADVar **v, Double *w)
553 static inline void Weighted_GradcompVec(size_t n, size_t *np, ADVar ***v, Double **w)
555 static inline void Outvar_Gradcomp(ADVar &v)
557
558 /* We use #define to deal with bizarre templating rules that apparently */
559 /* require us to spell the some conversion explicitly */
560
561
562#define Ai const ADVari&
563#define AI const IndepADvar&
564#define D Double
565#define T2(r,f) \
566 r f <>(AI,AI);\
567 r f <>(Ai,AI);\
568 r f <>(AI,Ai);\
569 r f <>(Ttype,AI);\
570 r f <>(double,AI);\
571 r f <>(long,AI);\
572 r f <>(int,AI);\
573 r f <>(AI,Ttype);\
574 r f <>(AI,double);\
575 r f <>(AI,long);\
576 r f <>(AI,int);
577#define T1(f) friend ADVari& f<> (AI);
578
579#define F friend ADVari&
580T2(F, operator+)
581T2(F, operator-)
582T2(F, operator*)
583T2(F, operator/)
584T2(F, atan2)
585T2(F, max)
586T2(F, min)
587T2(F, pow)
588#undef F
589#define F friend int
590T2(F, operator<)
591T2(F, operator<=)
592T2(F, operator==)
593T2(F, operator!=)
594T2(F, operator>=)
595T2(F, operator>)
596
597T1(operator+)
598T1(operator-)
599T1(abs)
600T1(acos)
601T1(acosh)
602T1(asin)
603T1(asinh)
604T1(atan)
605T1(atanh)
606T1(cos)
607T1(cosh)
608T1(exp)
609T1(log)
610T1(log10)
611T1(sin)
612T1(sinh)
613T1(sqrt)
614T1(tan)
615T1(tanh)
616T1(fabs)
617T1(copy)
618
619#undef F
620#undef T1
621#undef T2
622#undef D
623#undef AI
624#undef Ai
625
626 };
627
628 template<typename Double> class
629ADvar: public IndepADvar<Double> { // an "active" variable
630 public:
632 template <typename U> struct apply { typedef ADvar<U> type; };
633
635 typedef typename IndepADVar::ADVari ADVari;
637 private:
638 void ADvar_ctr(Double d) {
639 ADVari *x;
640 if (ConstADVari::cadc.fpval_implies_const)
641 x = new ConstADVari(d);
642 else {
643#ifdef RAD_AUTO_AD_Const //{
644 x = new ADVari((IndepADVar*)this, d);
645 x->ADvari_padv(this);
646#else
647 x = new ADVari(d);
648#endif //}
649 }
650 this->cv = x;
651 }
652 public:
653 friend class ADvar1<Double>;
655 inline ADvar() { /* cv = 0; */ }
656 inline ADvar(Ttype d) { ADvar_ctr(d); }
657 inline ADvar(double i) { ADvar_ctr(Double(i)); }
658 inline ADvar(int i) { ADvar_ctr(Double(i)); }
659 inline ADvar(long i) { ADvar_ctr(Double(i)); }
660 inline ~ADvar() {}
661#ifdef RAD_AUTO_AD_Const //{{
662 inline ADvar(IndepADVar &x) {
663 this->cv = x.cv ? new ADVar1(this, x) : 0;
664 }
665 inline ADvar(ADVari &x) { this->cv = &x; x.ADvari_padv(this); }
666 inline ADvar& operator=(IndepADVar &x) {
667 if (this->cv)
668 this->cv->padv = 0;
669 this->cv = new ADVar1(this,x);
670 return *this;
671 }
672 inline ADvar& operator=(ADVari &x) {
673 if (this->cv)
674 this->cv->padv = 0;
675 this->cv = new ADVar1(this, x);
676 return *this;
677 }
678#else //}{
679 friend ADvar& ADvar_operatoreq<>(ADvar*, const ADVari&);
680#ifdef RAD_EQ_ALIAS
681 /* allow aliasing v and w after "v = w;" */
682 inline ADvar(const IndepADVar &x) {
683 this->cv = (ADVari*)x.cv;
684 }
685 inline ADvar(const ADVari &x) { this->cv = (ADVari*)&x; }
686 inline ADvar& operator=(IndepADVar &x) {
687 this->cv = (ADVari*)x.cv; return *this;
688 }
689 inline ADvar& operator=(const ADVari &x) { this->cv = (ADVari*)&x; return *this; }
690#else
692 if (x.cv) {
693 this->cv = new ADVar1(x.cv->Val, &this->cv->adc.One, x.cv);
694 }
695 else
696 this->cv = 0;
697 }
698 ADvar(const ADvar&x) {
699 if (x.cv) {
700 this->cv = new ADVar1(x.cv->Val, &this->cv->adc.One, (ADVari*)x.cv);
701 }
702 else
703 this->cv = 0;
704 }
705 ADvar(const ADVari &x) {
706 this->cv = new ADVar1(x.Val, &this->cv->adc.One, &x);
707 }
708 inline ADvar& operator=(const ADVari &x) { return ADvar_operatoreq(this,x); };
709#endif /* RAD_EQ_ALIAS */
710#endif /* RAD_AUTO_AD_Const */ //}}
711 ADvar& operator=(Double);
712 ADvar& operator+=(const ADVari&);
713 ADvar& operator+=(Double);
714 ADvar& operator-=(const ADVari&);
715 ADvar& operator-=(Double);
716 ADvar& operator*=(const ADVari&);
717 ADvar& operator*=(Double);
718 ADvar& operator/=(const ADVari&);
719 ADvar& operator/=(Double);
720 inline static bool get_fpval_implies_const(void)
721 { return ConstADVari::cadc.fpval_implies_const; }
722 inline static void set_fpval_implies_const(bool newval)
723 { ConstADVari::cadc.fpval_implies_const = newval; }
724 inline static bool setget_fpval_implies_const(bool newval) {
725 bool oldval = ConstADVari::cadc.fpval_implies_const;
726 ConstADVari::cadc.fpval_implies_const = newval;
727 return oldval;
728 }
729 static inline void Gradcomp() { ADcontext<Double>::Gradcomp(); }
730 static inline void aval_reset() { ConstADVari::aval_reset(); }
731 static inline void Weighted_Gradcomp(size_t n, ADvar **v, Double *w)
733 static inline void Weighted_GradcompVec(size_t n, size_t *np, ADvar ***v, Double **w)
735 static inline void Outvar_Gradcomp(ADvar &v)
737 };
738
739template<typename Double>
740 inline void AD_Const1(Double *notused, const IndepADvar<Double>&v)
742
743template<typename Double>
744 inline void AD_Const(const IndepADvar<Double>&v) { AD_Const1((Double*)0, v); }
745
746 template<typename Double> class
748 public:
750 typedef typename ADVar::ADVar1 ADVar1;
751 typedef typename ADVar::ADVari ADVari;
755 private: // disable op=
764 void ConstADvar_ctr(Double);
765 public:
766 ConstADvar(Ttype d) { ConstADvar_ctr(d); }
767 ConstADvar(double i) { ConstADvar_ctr(Double(i)); }
768 ConstADvar(int i) { ConstADvar_ctr(Double(i)); }
769 ConstADvar(long i) { ConstADvar_ctr(Double(i)); }
770 ConstADvar(const IndepADVar&);
771 ConstADvar(const ConstADvar&);
772 ConstADvar(const ADVari&);
773 inline ~ConstADvar() {}
774#ifdef RAD_NO_CONST_UPDATE
775 private:
776#endif
777 ConstADvar();
778 inline ConstADvar& operator=(Double d) {
779 this->cv->Val = d;
780 return *this;
781 }
783 this->cv->Val = d.Val;
784 return *this;
785 }
786 };
787
788#define ADvar_nd ADvar
789
790 template<typename Double> class
791ADvar1s: public ADvar1<Double> { // unary ops with partial "a"
792 public:
794 typedef typename ADVar1::ADVari ADVari;
795 Double a;
796 ADvar1s(Double val1, Double a1, const ADVari *c1): ADVar1(val1,&a,c1), a(a1) {}
797#ifdef RAD_AUTO_AD_Const
798 typedef typename ADVar1::ADVar ADVar;
799 ADvar1s(Double val1, Double a1, const ADVari *c1, const ADVar *v):
800 ADVar1(val1,&a,c1,v), a(a1) {}
801#endif
802 };
803
804 template<typename Double> class
805ADvar2: public ADvari<Double> { // basic binary ops
806 public:
809 ADvar2(Double val1): ADVari(val1) {}
810 DErp dL, dR;
811 ADvar2(Double val1, const ADVari *Lcv, const Double *Lc,
812 const ADVari *Rcv, const Double *Rc):
813 ADVari(val1) {
814 dR.next = DErp::LastDerp;
815 dL.next = &dR;
816 DErp::LastDerp = &dL;
817 dL.a = Lc;
818 dL.c = Lcv;
819 dR.a = Rc;
820 dR.c = Rcv;
821 dL.b = dR.b = this;
822 }
823#ifdef RAD_AUTO_AD_Const
824 typedef ADvar<Double> ADVar;
825 ADvar2(Double val1, const ADVari *Lcv, const Double *Lc,
826 const ADVari *Rcv, const Double *Rc, ADVar *v):
827 ADVari(val1) {
828 dR.next = DErp::LastDerp;
829 dL.next = &dR;
830 DErp::LastDerp = &dL;
831 dL.a = Lc;
832 dL.c = Lcv;
833 dR.a = Rc;
834 dR.c = Rcv;
835 dL.b = dR.b = this;
836 Lcv->padv = 0;
837 this->ADvari_padv(v);
838 }
839#endif
840 };
841
842 template<typename Double> class
843ADvar2q: public ADvar2<Double> { // binary ops with partials "a", "b"
844 public:
846 typedef typename ADVar2::ADVari ADVari;
847 typedef typename ADVar2::DErp DErp;
848 Double a, b;
849 ADvar2q(Double val1, Double Lp, Double Rp, const ADVari *Lcv, const ADVari *Rcv):
850 ADVar2(val1), a(Lp), b(Rp) {
851 this->dR.next = DErp::LastDerp;
852 this->dL.next = &this->dR;
853 DErp::LastDerp = &this->dL;
854 this->dL.a = &a;
855 this->dL.c = Lcv;
856 this->dR.a = &b;
857 this->dR.c = Rcv;
858 this->dL.b = this->dR.b = this;
859 }
860#ifdef RAD_AUTO_AD_Const
861 typedef typename ADVar2::ADVar ADVar;
862 ADvar2q(Double val1, Double Lp, Double Rp, const ADVari *Lcv,
863 const ADVari *Rcv, const ADVar *v):
864 ADVar2(val1), a(Lp), b(Rp) {
865 this->dR.next = DErp::LastDerp;
866 this->dL.next = &this->dR;
867 DErp::LastDerp = &this->dL;
868 this->dL.a = &a;
869 this->dL.c = Lcv;
870 this->dR.a = &b;
871 this->dR.c = Rcv;
872 this->dL.b = this->dR.b = this;
873 Lcv->padv = 0;
874 this->ADvari_padv(v);
875 }
876#endif
877 };
878
879 template<typename Double> class
880ADvarn: public ADvari<Double> { // n-ary ops with partials "a"
881 public:
885 int n;
886 Double *a;
888 ADvarn(Double val1, int n1, const IndepADVar *x, const Double *g):
889 ADVari(val1), n(n1) {
890 DErp *d1, *dlast;
891 Double *a1;
892 int i;
893
894 a1 = a = (Double*)ADVari::adc.Memalloc(n*sizeof(*a));
895 d1 = Da = (DErp*)ADVari::adc.Memalloc(n*sizeof(DErp));
896 dlast = DErp::LastDerp;
897 for(i = 0; i < n1; i++, d1++) {
898 d1->next = dlast;
899 dlast = d1;
900 a1[i] = g[i];
901 d1->a = &a1[i];
902 d1->b = this;
903 d1->c = x[i].cv;
904 }
905 DErp::LastDerp = dlast;
906 }
907 };
908
909template<typename Double>
910 inline ADvari<Double>& operator+(const ADvari<Double> &T) { return *(ADvari<Double>*)&T; }
911
912template<typename Double>
913 inline int operator<(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val < R.Val; }
914template<typename Double>
915 inline int operator<(const ADvari<Double> &L, Double R) { return L.Val < R; }
916template<typename Double>
917 inline int operator<(Double L, const ADvari<Double> &R) { return L < R.Val; }
918
919template<typename Double>
920 inline int operator<=(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val <= R.Val; }
921template<typename Double>
922 inline int operator<=(const ADvari<Double> &L, Double R) { return L.Val <= R; }
923template<typename Double>
924 inline int operator<=(Double L, const ADvari<Double> &R) { return L <= R.Val; }
925
926template<typename Double>
927 inline int operator==(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val == R.Val; }
928template<typename Double>
929 inline int operator==(const ADvari<Double> &L, Double R) { return L.Val == R; }
930template<typename Double>
931 inline int operator==(Double L, const ADvari<Double> &R) { return L == R.Val; }
932
933template<typename Double>
934 inline int operator!=(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val != R.Val; }
935template<typename Double>
936 inline int operator!=(const ADvari<Double> &L, Double R) { return L.Val != R; }
937template<typename Double>
938 inline int operator!=(Double L, const ADvari<Double> &R) { return L != R.Val; }
939
940template<typename Double>
941 inline int operator>=(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val >= R.Val; }
942template<typename Double>
943 inline int operator>=(const ADvari<Double> &L, Double R) { return L.Val >= R; }
944template<typename Double>
945 inline int operator>=(Double L, const ADvari<Double> &R) { return L >= R.Val; }
946
947template<typename Double>
948 inline int operator>(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val > R.Val; }
949template<typename Double>
950 inline int operator>(const ADvari<Double> &L, Double R) { return L.Val > R; }
951template<typename Double>
952 inline int operator>(Double L, const ADvari<Double> &R) { return L > R.Val; }
953
954template<typename Double>
955 inline void *ADcontext<Double>::Memalloc(size_t len) {
956 if (Mleft >= len)
957 return Mbase + (Mleft -= len);
958 return new_ADmemblock(len);
959 }
960
961template<typename Double>
962 inline Derp<Double>::Derp(const ADVari *c1): c(c1) {
963 next = LastDerp;
964 LastDerp = this;
965 }
966
967template<typename Double>
968 inline Derp<Double>::Derp(const Double *a1, const ADVari *c1): a(a1), c(c1) {
969 next = LastDerp;
970 LastDerp = this;
971 }
972
973
974template<typename Double>
975 inline Derp<Double>::Derp(const Double *a1, const ADVari *b1, const ADVari *c1): a(a1), b(b1), c(c1) {
976 next = LastDerp;
977 LastDerp = this;
978 }
979
980/**** radops ****/
981
982template<typename Double> Derp<Double> *Derp<Double>::LastDerp = 0;
983template<typename Double> ADcontext<Double> ADvari<Double>::adc;
984template<typename Double> const Double ADcontext<Double>::One = 1.;
985template<typename Double> const Double ADcontext<Double>::negOne = -1.;
986template<typename Double> CADcontext<Double> ConstADvari<Double>::cadc;
987template<typename Double> ConstADvari<Double> *ConstADvari<Double>::lastcad;
988
989template<typename Double> ADvari<Double>* ADvari<Double>::First_ADvari;
991
992#ifdef RAD_DEBUG
993#ifndef RAD_DEBUG_gcgen1
994#define RAD_DEBUG_gcgen1 -1
995#endif
996template<typename Double> int ADvari<Double>::gcgen_cur;
997template<typename Double> int ADvari<Double>::last_opno;
998template<typename Double> int ADvari<Double>::zap_gcgen;
999template<typename Double> int ADvari<Double>::zap_gcgen1 = RAD_DEBUG_gcgen1;
1000template<typename Double> int ADvari<Double>::zap_opno;
1001template<typename Double> FILE *ADvari<Double>::debug_file;
1002#endif
1003
1004
1005template<typename Double> ADcontext<Double>::ADcontext()
1006{
1007 First = new ADMemblock;
1008 First->next = 0;
1009 Busy = First;
1010 Free = Oldbusy = 0;
1011 Mbase = (char*)First->memblk;
1012 Mleft = sizeof(First->memblk);
1013 ncur = nmax = rad_need_reinit = 0;
1014#ifdef RAD_DEBUG_BLOCKKEEP
1015 rad_busy_blocks = 0;
1016 rad_mleft_save = 0;
1017 rad_Oldcurmb = 0;
1018#endif
1019 }
1020
1021template<typename Double> void*
1023{
1024 ADMemblock *mb, *mb0, *mb1, *mbf, *x;
1025#ifdef RAD_AUTO_AD_Const
1026 ADVari *a, *anext;
1028#ifdef RAD_Const_WARN
1029 ADVari *cv;
1030 int i, j;
1031#endif
1032#endif /*RAD_AUTO_AD_Const*/
1033
1034 if (rad_need_reinit && this == &ADVari::adc) {
1035 rad_need_reinit = 0;
1036 nmax = 0;
1037 DErp::LastDerp = 0;
1038 ADVari::Last_ADvari = &ADVari::First_ADvari;
1039#ifdef RAD_DEBUG_BLOCKKEEP
1040 Mleft = rad_mleft_save;
1041 if (Mleft < sizeof(First->memblk))
1042 _uninit_f2c(Mbase + Mleft,
1043 UninitType<Double>::utype,
1044 (sizeof(First->memblk) - Mleft)
1045 /sizeof(typename Sacado::ValueType<Double>::type));
1046 if ((mb = Busy->next)) {
1047 mb0 = rad_Oldcurmb;
1048 for(;; mb = mb->next) {
1049 _uninit_f2c(mb->memblk,
1050 UninitType<Double>::utype,
1051 sizeof(First->memblk)
1052 /sizeof(typename Sacado::ValueType<Double>::type));
1053 if (mb == mb0)
1054 break;
1055 }
1056 }
1057 rad_Oldcurmb = Busy;
1058 if (rad_busy_blocks >= RAD_DEBUG_BLOCKKEEP) {
1059 rad_busy_blocks = 0;
1060 rad_Oldcurmb = 0;
1061 mb0 = 0;
1062 mbf = Free;
1063 for(mb = Busy; mb != mb0; mb = mb1) {
1064 mb1 = mb->next;
1065 mb->next = mbf;
1066 mbf = mb;
1067 }
1068 Free = mbf;
1069 Busy = mb;
1070 Mbase = (char*)First->memblk;
1071 Mleft = sizeof(First->memblk);
1072 }
1073
1074#else /* !RAD_DEBUG_BLOCKKEEP */
1075
1076 mb0 = First;
1077 mbf = Free;
1078 for(mb = Busy; mb != mb0; mb = mb1) {
1079 mb1 = mb->next;
1080 mb->next = mbf;
1081 mbf = mb;
1082 }
1083 Free = mbf;
1084 Busy = mb;
1085 Mbase = (char*)First->memblk;
1086 Mleft = sizeof(First->memblk);
1087#endif /*RAD_DEBUG_BLOCKKEEP*/
1088#ifdef RAD_AUTO_AD_Const // {
1089 if ((a = ADVari::First_ADvari)) {
1090 do {
1091 anext = a->Next;
1092 if ((v = (IndepADvar<Double> *)a->padv)) {
1093#ifdef RAD_Const_WARN
1094 if ((i = a->opno) > 0)
1095 i = -i;
1096 j = a->gcgen;
1097 v->cv = cv = new ADVari(v, a->Val);
1098 cv->opno = i;
1099 cv->gcgen = j;
1100#else
1101 v->cv = new ADVari(v, a->Val);
1102#endif
1103 }
1104 }
1105 while((a = anext));
1106 }
1107#endif // } RAD_AUTO_AD_Const
1108 if (Mleft >= len)
1109 return Mbase + (Mleft -= len);
1110 }
1111
1112 if ((x = Free))
1113 Free = x->next;
1114 else
1115 x = new ADMemblock;
1116#ifdef RAD_DEBUG_BLOCKKEEP
1117 rad_busy_blocks++;
1118#endif
1119 x->next = Busy;
1120 Busy = x;
1121 return (Mbase = (char*)x->memblk) +
1122 (Mleft = sizeof(First->memblk) - len);
1123 }
1124
1125 template<typename Double> void
1127{
1128 ADMemblock *mb, *mb0, *mb1, *mbf;
1129 ADVari *a;
1131 Double *d, *de;
1132 size_t len;
1133
1134 if (!rad_need_reinit) {
1135 rad_mleft_save = Mleft;
1136 Oldbusy = Busy;
1137 }
1138 len = n * sizeof(Double);
1139 ncur = n;
1140 if (!nmax) {
1141 if (n > 0)
1142 goto aval_alloc;
1143 }
1144 else if (n > nmax) {
1145 mb0 = Oldbusy;
1146 mb = Busy;
1147 if (mb != mb0) {
1148 mbf = Free;
1149 do {
1150 mb1 = mb->next;
1151 mb->next = mbf;
1152 mbf = mb;
1153 mb = mb1;
1154 }
1155 while(mb != mb0);
1156 Free = mbf;
1157 }
1158 Mleft = rad_mleft_save;
1159 Busy = Oldbusy;
1160 aval_alloc:
1161 rad_need_reinit = 0;
1162 nmax = n;
1163 *ADVari::Last_ADvari = 0;
1164 for(a = ADVari::First_ADvari; a; a = a->Next) {
1165 d = a->aval = (Double*)Memalloc(len);
1166 de = d + n;
1167 do *d = 0.; while(++d < de);
1168 }
1170 d = c->aval = (Double*)Memalloc(len);
1171 de = d + n;
1172 do *d = 0.; while(++d < de);
1173 }
1174 }
1175 else if (!n) {
1176 nmax = 0;
1177 for(a = ADVari::First_ADvari; a; a = a->Next)
1178 a->aval = 0;
1179 }
1180 else for(a = ADVari::First_ADvari; a; a = a->Next) {
1181 d = a->aval;
1182 de = d + n;
1183 do *d = 0.; while(++d < de);
1184 }
1185 Mleft = 0;
1186 rad_need_reinit = 1;
1187#ifdef RAD_DEBUG
1188 if (ADVari::gcgen_cur == ADVari::zap_gcgen1) {
1189 const char *fname;
1190 if (!(fname = getenv("RAD_DEBUG_FILE")))
1191 fname = "rad_debug.out";
1192 else if (!*fname)
1193 fname = 0;
1194 if (fname)
1195 ADVari::debug_file = fopen(fname, "w");
1196 ADVari::zap_gcgen1 = -1;
1197 }
1198#endif
1199 }
1200
1201 template<typename Double> void
1203{
1204 DErp *d;
1205#ifdef RAD_AUTO_AD_Const
1206 ADVari *a, *anext;
1208#ifdef RAD_Const_WARN
1209 ADVari *cv;
1210 int i, j;
1211#endif
1212#endif /*RAD_AUTO_AD_Const*/
1213
1214 ADVari::adc.derp_init(1);
1215 if ((d = DErp::LastDerp) != 0) {
1216 *d->b->aval = 1;
1217#ifdef RAD_DEBUG
1218 if (ADVari::debug_file)
1219 do {
1220 fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
1221 d->c->opno, d->b->opno, &d->c->aval, *d->a, *d->b->aval);
1222 d->c->aval += *d->a * d->b->aval;
1223 fprintf(ADVari::debug_file, " = %g\n", *d->c->aval);
1224 } while((d = d->next));
1225 else
1226#endif
1227 do *d->c->aval += *d->a * *d->b->aval;
1228 while((d = d->next));
1229 }
1230#ifdef RAD_DEBUG
1231 if (ADVari::debug_file) {
1232 fclose(ADVari::debug_file);
1233 ADVari::debug_file = 0;
1234 }
1235 ADVari::gcgen_cur++;
1236 ADVari::last_opno = 0;
1237#endif
1238 }
1239
1240 template<typename Double> void
1242{
1243 size_t i;
1244 DErp *d;
1245#ifdef RAD_AUTO_AD_Const
1246 ADVari *a, *anext;
1248#ifdef RAD_Const_WARN
1249 ADVari *cv;
1250 int j, k;
1251#endif
1252#endif /*RAD_AUTO_AD_Const*/
1253
1254 ADVari::adc.derp_init(1);
1255
1256 if ((d = DErp::LastDerp) != 0) {
1257 for(i = 0; i < n; i++)
1258 *V[i]->cv->aval = w[i];
1259#ifdef RAD_DEBUG
1260 if (ADVari::debug_file)
1261 do {
1262 fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
1263 d->c->opno, d->b->opno, *d->c->aval, *d->a, *d->b->aval);
1264 *d->c->aval += *d->a * *d->b->aval;
1265 fprintf(ADVari::debug_file, " = %g\n", *d->c->aval);
1266 } while((d = d->next));
1267 else
1268#endif
1269 do *d->c->aval += *d->a * *d->b->aval;
1270 while((d = d->next));
1271 }
1272#ifdef RAD_DEBUG
1273 if (ADVari::debug_file) {
1274 fclose(ADVari::debug_file);
1275 ADVari::debug_file = 0;
1276 }
1277 if (ADVari::debug_file)
1278 fflush(ADVari::debug_file);
1279 ADVari::gcgen_cur++;
1280 ADVari::last_opno = 0;
1281#endif
1282 }
1283
1284 template<typename Double> void
1285ADcontext<Double>::Weighted_GradcompVec(size_t n, size_t *np, ADVar ***V, Double **w)
1286{
1287 size_t i, ii, ni;
1288 ADVar **Vi;
1289 DErp *d;
1290 Double A, *D, *D1, *De, *wi;
1291#ifdef RAD_AUTO_AD_Const
1292 ADVari *a, *anext;
1294#ifdef RAD_Const_WARN
1295 ADVari *cv;
1296 int j, k;
1297#endif
1298#endif /*RAD_AUTO_AD_Const*/
1299
1300 ADVari::adc.derp_init(n);
1301 if (!n)
1302 return;
1303
1304 if ((d = DErp::LastDerp) != 0) {
1305 for(i = 0; i < n; i++) {
1306 ni = np[i];
1307 wi = w[i];
1308 Vi = V[i];
1309 for(ii = 0; ii < ni; ++ii)
1310 Vi[ii]->cv->aval[i] = wi[ii];
1311 }
1312#ifdef RAD_DEBUG
1313 if (ADVari::debug_file)
1314 do {
1315 D = d->c->aval;
1316 De = D + n;
1317 D1 = d->b->aval;
1318 A = *d->a;
1319 for(i = 0; i < n; ++i, ++D, ++D1) {
1320 fprintf(ADVari::debug_file, "%d:\t%d\t%d\t%g + %g * %g",
1321 i, d->c->opno, d->b->opno, *D, A, *D1);
1322 *D += A * *D1;
1323 fprintf(ADVari::debug_file, " = %g\n", *D);
1324 }
1325 } while((d = d->next));
1326 else
1327#endif
1328 do {
1329 D = d->c->aval;
1330 De = D + n;
1331 D1 = d->b->aval;
1332 A = *d->a;
1333 do *D++ += A * *D1++; while(D < De);
1334 }
1335 while((d = d->next));
1336 }
1337#ifdef RAD_DEBUG
1338 if (ADVari::debug_file) {
1339 fclose(ADVari::debug_file);
1340 ADVari::debug_file = 0;
1341 }
1342 if (ADVari::debug_file)
1343 fflush(ADVari::debug_file);
1344 ADVari::gcgen_cur++;
1345 ADVari::last_opno = 0;
1346#endif
1347 }
1348
1349 template<typename Double> void
1351{
1352 Double w = 1;
1353 ADVar *v = &V;
1355 }
1356
1357 template<typename Double>
1359{
1360 cv = new ADVari(d);
1361 }
1362
1363 template<typename Double>
1365{
1366 cv = new ADVari(Double(i));
1367 }
1368
1369 template<typename Double>
1371{
1372 cv = new ADVari(Double(i));
1373 }
1374
1375 template<typename Double>
1377{
1378 cv = new ADVari(Double(i));
1379 }
1380
1381 template<typename Double>
1383{
1384 this->cv = new ConstADVari(0.);
1385 }
1386
1387 template<typename Double> void
1389{
1390 this->cv = new ConstADVari(d);
1391 }
1392
1393 template<typename Double>
1395{
1396 ConstADVari *y = new ConstADVari(x.cv->Val);
1397 DErp *d = new DErp(&x.adc.One, y, x.cv);
1398 this->cv = y;
1399 }
1400
1401 template<typename Double>
1403{
1404 ConstADVari *y = new ConstADVari(x.cv->Val);
1405 DErp *d = new DErp(&x.cv->adc.One, y, (ADVari*)x.cv);
1406 this->cv = y;
1407 }
1408
1409 template<typename Double>
1411{
1412 ConstADVari *y = new ConstADVari(x.Val);
1413 DErp *d = new DErp(&x.adc.One, y, &x);
1414 this->cv = y;
1415 }
1416
1417 template<typename Double>
1418 void
1420{
1421 typedef ConstADvari<Double> ConstADVari;
1422
1423 ConstADVari *ncv = new ConstADVari(v.val());
1424#ifdef RAD_AUTO_AD_Const
1425 v.cv->padv = 0;
1426#endif
1427 v.cv = ncv;
1428 }
1429
1430 template<typename Double>
1431 int
1433{
1434 return 1;
1435 }
1436
1437 template<typename Double>
1438 void
1440{
1442 while(x) {
1443 x->aval = 0;
1444 x = x->prevcad;
1445 }
1446 }
1447
1448#ifdef RAD_AUTO_AD_Const
1449
1450 template<typename Double>
1451ADvari<Double>::ADvari(const IndepADVar *x, Double d): Val(d), aval(0)
1452{
1453 this->ADvari_padv(x);
1454 Linkup();
1455 }
1456
1457 template<typename Double>
1458ADvar1<Double>::ADvar1(const IndepADVar *x, const IndepADVar &y):
1459 ADVari(y.cv->Val), d((const Double*)&ADcontext<Double>::One, (ADVari*)this, y.cv)
1460{
1461 this->ADvari_padv(x);
1462 }
1463
1464 template<typename Double>
1465ADvar1<Double>::ADvar1(const IndepADVar *x, const ADVari &y):
1466 ADVari(y.Val), d((const Double*)&ADcontext<Double>::One, this, &y)
1467{
1468 this->ADvari_padv(x);
1469 }
1470
1471#else /* !RAD_AUTO_AD_Const */
1472
1473 template<typename Double>
1474 IndepADvar<Double>&
1476{
1477 This->cv = new ADvar1<Double>(x.Val, &x.adc.One, &x);
1478 return *(IndepADvar<Double>*) This;
1479 }
1480
1481 template<typename Double>
1484{
1485 This->cv = new ADvar1<Double>(x.Val, &x.adc.One, &x);
1486 return *(ADvar<Double>*) This;
1487 }
1488
1489#endif /* RAD_AUTO_AD_Const */
1490
1491
1492 template<typename Double>
1495{
1496#ifdef RAD_AUTO_AD_Const
1497 if (this->cv)
1498 this->cv->padv = 0;
1499 this->cv = new ADVari(this,d);
1500#else
1501 this->cv = new ADVari(d);
1502#endif
1503 return *this;
1504 }
1505
1506 template<typename Double>
1509{
1510#ifdef RAD_AUTO_AD_Const
1511 if (this->cv)
1512 this->cv->padv = 0;
1513 this->cv = new ADVari(this,d);
1514#else
1515 this->cv = ConstADVari::cadc.fpval_implies_const
1516 ? new ConstADVari(d)
1517 : new ADVari(d);
1518#endif
1519 return *this;
1520 }
1521
1522 template<typename Double>
1525 return *(new ADvar1<Double>(-T.Val, &T.adc.negOne, &T));
1526 }
1527
1528 template<typename Double>
1529 ADvari<Double>&
1531 return *(new ADvar2<Double>(L.Val + R.Val, &L, &L.adc.One, &R, &L.adc.One));
1532 }
1533
1534#ifdef RAD_AUTO_AD_Const
1535#define RAD_ACA ,this
1536#else
1537#define RAD_ACA /*nothing*/
1538#endif
1539
1540 template<typename Double>
1541 ADvar<Double>&
1543 ADVari *Lcv = this->cv;
1544 this->cv = new ADvar2<Double>(Lcv->Val + R.Val, Lcv, &R.adc.One, &R, &R.adc.One RAD_ACA);
1545 return *this;
1546 }
1547
1548 template<typename Double>
1550operator+(const ADvari<Double> &L, Double R) {
1551 return *(new ADvar1<Double>(L.Val + R, &L.adc.One, &L));
1552 }
1553
1554 template<typename Double>
1555 ADvar<Double>&
1557 ADVari *tcv = this->cv;
1558 this->cv = new ADVar1(tcv->Val + R, &tcv->adc.One, tcv RAD_ACA);
1559 return *this;
1560 }
1561
1562 template<typename Double>
1564operator+(Double L, const ADvari<Double> &R) {
1565 return *(new ADvar1<Double>(L + R.Val, &R.adc.One, &R));
1566 }
1567
1568 template<typename Double>
1569 ADvari<Double>&
1571 return *(new ADvar2<Double>(L.Val - R.Val, &L, &L.adc.One, &R, &L.adc.negOne));
1572 }
1573
1574 template<typename Double>
1575 ADvar<Double>&
1577 ADVari *Lcv = this->cv;
1578 this->cv = new ADvar2<Double>(Lcv->Val - R.Val, Lcv, &R.adc.One, &R, &R.adc.negOne RAD_ACA);
1579 return *this;
1580 }
1581
1582 template<typename Double>
1584operator-(const ADvari<Double> &L, Double R) {
1585 return *(new ADvar1<Double>(L.Val - R, &L.adc.One, &L));
1586 }
1587
1588 template<typename Double>
1589 ADvar<Double>&
1591 ADVari *tcv = this->cv;
1592 this->cv = new ADVar1(tcv->Val - R, &tcv->adc.One, tcv RAD_ACA);
1593 return *this;
1594 }
1595
1596 template<typename Double>
1598operator-(Double L, const ADvari<Double> &R) {
1599 return *(new ADvar1<Double>(L - R.Val, &R.adc.negOne, &R));
1600 }
1601
1602 template<typename Double>
1603 ADvari<Double>&
1605 return *(new ADvar2<Double>(L.Val * R.Val, &L, &R.Val, &R, &L.Val));
1606 }
1607
1608 template<typename Double>
1609 ADvar<Double>&
1611 ADVari *Lcv = this->cv;
1612 this->cv = new ADvar2<Double>(Lcv->Val * R.Val, Lcv, &R.Val, &R, &Lcv->Val RAD_ACA);
1613 return *this;
1614 }
1615
1616 template<typename Double>
1618operator*(const ADvari<Double> &L, Double R) {
1619 return *(new ADvar1s<Double>(L.Val * R, R, &L));
1620 }
1621
1622 template<typename Double>
1623 ADvar<Double>&
1625 ADVari *Lcv = this->cv;
1626 this->cv = new ADvar1s<Double>(Lcv->Val * R, R, Lcv RAD_ACA);
1627 return *this;
1628 }
1629
1630 template<typename Double>
1632operator*(Double L, const ADvari<Double> &R) {
1633 return *(new ADvar1s<Double>(L * R.Val, L, &R));
1634 }
1635
1636 template<typename Double>
1637 ADvari<Double>&
1639 Double Lv = L.Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv;
1640 return *(new ADvar2q<Double>(q, pL, -q*pL, &L, &R));
1641 }
1642
1643 template<typename Double>
1644 ADvar<Double>&
1646 ADVari *Lcv = this->cv;
1647 Double Lv = Lcv->Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv;
1648 this->cv = new ADvar2q<Double>(q, pL, -q*pL, Lcv, &R RAD_ACA);
1649 return *this;
1650 }
1651
1652 template<typename Double>
1654operator/(const ADvari<Double> &L, Double R) {
1655 return *(new ADvar1s<Double>(L.Val / R, 1./R, &L));
1656 }
1657
1658 template<typename Double>
1659 ADvari<Double>&
1660operator/(Double L, const ADvari<Double> &R) {
1661 Double recip = 1. / R.Val;
1662 Double q = L * recip;
1663 return *(new ADvar1s<Double>(q, -q*recip, &R));
1664 }
1665
1666 template<typename Double>
1667 ADvar<Double>&
1669 ADVari *Lcv = this->cv;
1670 this->cv = new ADvar1s<Double>(Lcv->Val / R, 1./R, Lcv RAD_ACA);
1671 return *this;
1672 }
1673
1674 template<typename Double>
1677 Double t = v.Val;
1678 return *(new ADvar1s<Double>(std::acos(t), -1./std::sqrt(1. - t*t), &v));
1679 }
1680
1681 template<typename Double>
1682 ADvari<Double>&
1684 Double t = v.Val, t1 = std::sqrt(t*t - 1.);
1685 return *(new ADvar1s<Double>(std::log(t + t1), 1./t1, &v));
1686 }
1687
1688 template<typename Double>
1689 ADvari<Double>&
1691 Double t = v.Val;
1692 return *(new ADvar1s<Double>(std::asin(t), 1./std::sqrt(1. - t*t), &v));
1693 }
1694
1695 template<typename Double>
1696 ADvari<Double>&
1698 Double t = v.Val, td = 1., t1 = std::sqrt(t*t + 1.);
1699 if (t < 0.) {
1700 t = -t;
1701 td = -1.;
1702 }
1703 return *(new ADvar1s<Double>(td*std::log(t + t1), 1./t1, &v));
1704 }
1705
1706 template<typename Double>
1707 ADvari<Double>&
1709 Double t = v.Val;
1710 return *(new ADvar1s<Double>(std::atan(t), 1./(1. + t*t), &v));
1711 }
1712
1713 template<typename Double>
1714 ADvari<Double>&
1716 Double t = v.Val;
1717 return *(new ADvar1s<Double>(0.5*std::log((1.+t)/(1.-t)), 1./(1. - t*t), &v));
1718 }
1719
1720 template<typename Double>
1721 ADvari<Double>&
1723 Double x = L.Val, y = R.Val, t = x*x + y*y;
1724 return *(new ADvar2q<Double>(std::atan2(x,y), y/t, -x/t, &L, &R));
1725 }
1726
1727 template<typename Double>
1728 ADvari<Double>&
1729atan2(Double x, const ADvari<Double> &R) {
1730 Double y = R.Val, t = x*x + y*y;
1731 return *(new ADvar1s<Double>(std::atan2(x,y), -x/t, &R));
1732 }
1733
1734 template<typename Double>
1735 ADvari<Double>&
1736atan2(const ADvari<Double> &L, Double y) {
1737 Double x = L.Val, t = x*x + y*y;
1738 return *(new ADvar1s<Double>(std::atan2(x,y), y/t, &L));
1739 }
1740
1741 template<typename Double>
1742 ADvari<Double>&
1744 const ADvari<Double> &x = L.Val >= R.Val ? L : R;
1745 return *(new ADvar1<Double>(x.Val, &x.adc.One, &x));
1746 }
1747
1748 template<typename Double>
1749 ADvari<Double>&
1750max(Double L, const ADvari<Double> &R) {
1751 if (L >= R.Val)
1752 return *(new ADvari<Double>(L));
1753 return *(new ADvar1<Double>(R.Val, &R.adc.One, &R));
1754 }
1755
1756 template<typename Double>
1757 ADvari<Double>&
1758max(const ADvari<Double> &L, Double R) {
1759 if (L.Val >= R)
1760 return *(new ADvar1<Double>(L.Val, &L.adc.One, &L));
1761 return *(new ADvari<Double>(R));
1762 }
1763
1764 template<typename Double>
1765 ADvari<Double>&
1767 const ADvari<Double> &x = L.Val <= R.Val ? L : R;
1768 return *(new ADvar1<Double>(x.Val, &x.adc.One, &x));
1769 }
1770
1771 template<typename Double>
1772 ADvari<Double>&
1773min(Double L, const ADvari<Double> &R) {
1774 if (L <= R.Val)
1775 return *(new ADvari<Double>(L));
1776 return *(new ADvar1<Double>(R.Val, &R.adc.One, &R));
1777 }
1778
1779 template<typename Double>
1780 ADvari<Double>&
1781min(const ADvari<Double> &L, Double R) {
1782 if (L.Val <= R)
1783 return *(new ADvar1<Double>(L.Val, &L.adc.One, &L));
1784 return *(new ADvari<Double>(R));
1785 }
1786
1787 template<typename Double>
1788 ADvari<Double>&
1790 return *(new ADvar1s<Double>(std::cos(v.Val), -std::sin(v.Val), &v));
1791 }
1792
1793 template<typename Double>
1794 ADvari<Double>&
1796 return *(new ADvar1s<Double>(std::cosh(v.Val), std::sinh(v.Val), &v));
1797 }
1798
1799 template<typename Double>
1800 ADvari<Double>&
1802 ADvar1<Double>* rcv = new ADvar1<Double>(std::exp(v.Val), &v);
1803 rcv->d.a = &rcv->Val;
1804 rcv->d.b = rcv;
1805 return *rcv;
1806 }
1807
1808 template<typename Double>
1809 ADvari<Double>&
1811 Double x = v.Val;
1812 return *(new ADvar1s<Double>(std::log(x), 1. / x, &v));
1813 }
1814
1815 template<typename Double>
1816 ADvari<Double>&
1818 static double num = 1. / std::log(10.);
1819 Double x = v.Val;
1820 return *(new ADvar1s<Double>(std::log10(x), num / x, &v));
1821 }
1822
1823 template<typename Double>
1824 ADvari<Double>&
1826 Double x = L.Val, y = R.Val, t = std::pow(x,y);
1827 return *(new ADvar2q<Double>(t, y*t/x, t*std::log(x), &L, &R));
1828 }
1829
1830 template<typename Double>
1831 ADvari<Double>&
1832pow(Double x, const ADvari<Double> &R) {
1833 Double t = std::pow(x,R.Val);
1834 return *(new ADvar1s<Double>(t, t*std::log(x), &R));
1835 }
1836
1837 template<typename Double>
1838 ADvari<Double>&
1839pow(const ADvari<Double> &L, Double y) {
1840 Double x = L.Val, t = std::pow(x,y);
1841 return *(new ADvar1s<Double>(t, y*t/x, &L));
1842 }
1843
1844 template<typename Double>
1845 ADvari<Double>&
1847 return *(new ADvar1s<Double>(std::sin(v.Val), std::cos(v.Val), &v));
1848 }
1849
1850 template<typename Double>
1851 ADvari<Double>&
1853 return *(new ADvar1s<Double>(std::sinh(v.Val), std::cosh(v.Val), &v));
1854 }
1855
1856 template<typename Double>
1857 ADvari<Double>&
1859 Double t = std::sqrt(v.Val);
1860 return *(new ADvar1s<Double>(t, 0.5/t, &v));
1861 }
1862
1863 template<typename Double>
1864 ADvari<Double>&
1866 Double t = std::cos(v.Val);
1867 return *(new ADvar1s<Double>(std::tan(v.Val), 1./(t*t), &v));
1868 }
1869
1870 template<typename Double>
1871 ADvari<Double>&
1873 Double t = 1. / std::cosh(v.Val);
1874 return *(new ADvar1s<Double>(std::tanh(v.Val), t*t, &v));
1875 }
1876
1877 template<typename Double>
1878 ADvari<Double>&
1880 Double t, p;
1881 p = 1;
1882 if ((t = v.Val) < 0) {
1883 t = -t;
1884 p = -p;
1885 }
1886 return *(new ADvar1s<Double>(t, p, &v));
1887 }
1888
1889 template<typename Double>
1890 ADvari<Double>&
1891fabs(const ADvari<Double> &v) { // Synonym for "abs"
1892 // "fabs" is not the best choice of name,
1893 // but this name is used at Sandia.
1894 Double t, p;
1895 p = 1;
1896 if ((t = v.Val) < 0) {
1897 t = -t;
1898 p = -p;
1899 }
1900 return *(new ADvar1s<Double>(t, p, &v));
1901 }
1902
1903 template<typename Double>
1904 ADvari<Double>&
1905ADf1(Double f, Double g, const ADvari<Double> &x) {
1906 return *(new ADvar1s<Double>(f, g, &x));
1907 }
1908
1909 template<typename Double>
1910 inline ADvari<Double>&
1911ADf1(Double f, Double g, const IndepADvar<Double> &x) {
1912 return *(new ADvar1s<Double>(f, g, x.cv));
1913 }
1914
1915 template<typename Double>
1917ADf2(Double f, Double gx, Double gy, const ADvari<Double> &x, const ADvari<Double> &y) {
1918 return *(new ADvar2q<Double>(f, gx, gy, &x, &y));
1919 }
1920
1921 template<typename Double>
1923ADf2(Double f, Double gx, Double gy, const ADvari<Double> &x, const IndepADvar<Double> &y) {
1924 return *(new ADvar2q<Double>(f, gx, gy, &x, y.cv));
1925 }
1926
1927 template<typename Double>
1929ADf2(Double f, Double gx, Double gy, const IndepADvar<Double> &x, const ADvari<Double> &y) {
1930 return *(new ADvar2q<Double>(f, gx, gy, x.cv, &y));
1931 }
1932
1933 template<typename Double>
1935ADf2(Double f, Double gx, Double gy, const IndepADvar<Double> &x, const IndepADvar<Double> &y) {
1936 return *(new ADvar2q<Double>(f, gx, gy, x.cv, y.cv));
1937 }
1938
1939 template<typename Double>
1941ADfn(Double f, int n, const IndepADvar<Double> *x, const Double *g) {
1942 return *(new ADvarn<Double>(f, n, x, g));
1943 }
1944
1945 template<typename Double>
1946 inline ADvari<Double>&
1947ADfn(Double f, int n, const ADvar<Double> *x, const Double *g) {
1948 return ADfn<Double>(f, n, (IndepADvar<Double>*)x, g);
1949 }
1950
1951 template<typename Double>
1952 inline Double
1954 return x.Val;
1955 }
1956
1957#undef RAD_ACA
1958#define A (ADvari<Double>*)
1959#ifdef RAD_Const_WARN
1960#define C(x) (((x)->opno < 0) ? RAD_Const_Warn(x) : 0, *A x)
1961#else
1962#define C(x) *A x
1963#endif
1964#define T template<typename Double> inline
1965#define F ADvari<Double>&
1966#define Ai const ADvari<Double>&
1967#define AI const IndepADvar<Double>&
1968#define D Double
1969#define T2(r,f) \
1970 T r f(Ai L, AI R) { return f(L, C(R.cv)); }\
1971 T r f(AI L, Ai R) { return f(C(L.cv), R); }\
1972 T r f(AI L, AI R) { return f(C(L.cv), C(R.cv)); }\
1973 T r f(AI L, D R) { return f(C(L.cv), R); }\
1974 T r f(Ai L, Dtype R) { return f(L, (D)R); }\
1975 T r f(AI L, Dtype R) { return f(C(L.cv), (D)R); }\
1976 T r f(Ai L, long R) { return f(L, (D)R); }\
1977 T r f(AI L, long R) { return f(C(L.cv), (D)R); }\
1978 T r f(Ai L, int R) { return f(L, (D)R); }\
1979 T r f(AI L, int R) { return f(C(L.cv), (D)R); }\
1980 T r f(D L, AI R) { return f(L, C(R.cv)); }\
1981 T r f(Dtype L, Ai R) { return f((D)L, R); }\
1982 T r f(Dtype L, AI R) { return f((D)L, C(R.cv)); }\
1983 T r f(long L, Ai R) { return f((D)L, R); }\
1984 T r f(long L, AI R) { return f((D)L, C(R.cv)); }\
1985 T r f(int L, Ai R) { return f((D)L, R); }\
1986 T r f(int L, AI R) { return f((D)L, C(R.cv)); }
1987
1988T2(F, operator+)
1989T2(F, operator-)
1990T2(F, operator*)
1991T2(F, operator/)
1992T2(F, atan2)
1993T2(F, pow)
1994T2(F, max)
1995T2(F, min)
1996T2(int, operator<)
1997T2(int, operator<=)
1998T2(int, operator==)
1999T2(int, operator!=)
2000T2(int, operator>=)
2001T2(int, operator>)
2002
2003#undef T2
2004#undef D
2005
2006#define T1(f)\
2007 T F f(AI x) { return f(C(x.cv)); }
2008
2009T1(operator+)
2010T1(operator-)
2011T1(abs)
2012T1(acos)
2013T1(acosh)
2014T1(asin)
2015T1(asinh)
2016T1(atan)
2017T1(atanh)
2018T1(cos)
2019T1(cosh)
2020T1(exp)
2021T1(log)
2022T1(log10)
2023T1(sin)
2024T1(sinh)
2025T1(sqrt)
2026T1(tan)
2027T1(tanh)
2028T1(fabs)
2029
2031{
2032 return *(new ADvar1<Double>(x.cv->Val, &ADcontext<Double>::One, (ADvari<Double>*)x.cv));
2033 }
2034
2036{ return *(new ADvar1<Double>(x.Val, &ADcontext<Double>::One, (ADvari<Double>*)&x)); }
2037
2038#undef T1
2039#undef AI
2040#undef Ai
2041#undef F
2042#undef T
2043#undef A
2044#undef C
2045#undef Ttype
2046#undef Dtype
2047#undef Padvinit
2048
2049} /* namespace RadVec */
2050} /* namespace Sacado */
2051
2052#undef SNS
2053
2054#endif /* SACADO_TRADVEC_H */
int RAD_Const_Warn(void *v)
expr val()
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
#define T1(r, f)
#define R
#define A
#define T2(r, f)
#define D
#define Ttype
#define RAD_ACA
#define Ttype
#define T1(f)
#define R
#define T2(r, f)
#define Padvinit
static void Weighted_Gradcomp(size_t, ADVar **, Double *)
static void Outvar_Gradcomp(ADVar &)
static const Double negOne
static void aval_reset(void)
static void Weighted_GradcompVec(size_t, size_t *, ADVar ***, Double **)
ADmemblock< Double > ADMemblock
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)
ADvar & operator+=(const ADVari &)
static void Weighted_Gradcomp(size_t n, ADvar **v, Double *w)
static void set_fpval_implies_const(bool newval)
static void Outvar_Gradcomp(ADvar &v)
ADvar & operator=(const ADVari &x)
ADvar(const IndepADVar &x)
ADvar & operator-=(const ADVari &)
ConstADvari< Double > ConstADVari
ADvar1< Double > ADVar1
IndepADvar< Double > IndepADVar
IndepADVar::ADVari ADVari
static bool setget_fpval_implies_const(bool newval)
ADvar & operator*=(const ADVari &)
ADvar & operator/=(const ADVari &)
static bool get_fpval_implies_const(void)
static void Weighted_GradcompVec(size_t n, size_t *np, ADvar ***v, Double **w)
IndepADvar< Double > IndepADVar
static ADvari * First_ADvari
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*=(ADVari &)
ConstADvar & operator+=(Double)
ConstADvar & operator-=(Double)
ADVar::ConstADVari ConstADVari
ConstADvar & operator=(Double d)
ConstADvar & operator/=(Double)
ConstADvar & operator=(ADVari &d)
ConstADvar & operator-=(ADVari &)
ConstADvar & operator*=(Double)
ConstADvar & operator/=(ADVari &)
ConstADvar & operator+=(ADVari &)
static CADcontext< Double > cadc
static ConstADvari * lastcad
ADvari< Double > ADVari
static void AD_Const(const IndepADvar &)
static void Weighted_Gradcomp(size_t n, ADVar **v, Double *w)
static void Outvar_Gradcomp(ADVar &v)
IndepADvar & operator=(IndepADvar &x)
static void Weighted_GradcompVec(size_t n, size_t *np, ADVar ***v, Double **w)
const double y
const char * p
ADvari< Double > & acos(const ADvari< Double > &v)
ADvari< Double > & min(const ADvari< Double > &L, const ADvari< Double > &R)
ADvari< Double > & log(const ADvari< Double > &v)
IndepADvar< Double > & ADvar_operatoreq(IndepADvar< Double > *, const ADvari< Double > &)
ADvari< Double > & ADf1(Double f, Double g, const IndepADvar< Double > &x)
ADvari< Double > & asinh(const ADvari< Double > &v)
void AD_Const1(Double *, const IndepADvar< Double > &)
ADvari< Double > & sinh(const ADvari< Double > &v)
ADvari< Double > & sin(const ADvari< Double > &v)
ADvari< Double > & cos(const ADvari< Double > &v)
ADvari< Double > & atanh(const ADvari< Double > &v)
ADvari< Double > & exp(const ADvari< Double > &v)
ADvari< Double > & tan(const ADvari< Double > &v)
ADvari< Double > & operator-(const ADvari< Double > &T)
void AD_Const(const IndepADvar< Double > &)
ADvari< Double > & operator*(const ADvari< Double > &L, const ADvari< Double > &R)
ADvari< Double > & tanh(const ADvari< Double > &v)
ADvari< Double > & abs(const ADvari< Double > &v)
ADvari< Double > & pow(const ADvari< Double > &L, const ADvari< Double > &R)
ADvari< Double > & ADf2(Double f, Double gx, Double gy, const IndepADvar< Double > &x, const IndepADvar< Double > &y)
ADvari< Double > & atan(const ADvari< Double > &v)
ADvari< Double > & ADfn(Double f, int n, const IndepADvar< Double > *x, const Double *g)
ADvari< Double > & acosh(const ADvari< Double > &v)
ADvari< Double > & operator/(const ADvari< Double > &L, const ADvari< Double > &R)
ADvari< Double > & operator+(const ADvari< Double > &T)
ADvari< Double > & max(const ADvari< Double > &L, const ADvari< Double > &R)
ADvari< Double > & fabs(const ADvari< Double > &v)
ADvari< Double > & log10(const ADvari< Double > &v)
ADvari< Double > & asin(const ADvari< Double > &v)
ADvari< Double > & cosh(const ADvari< Double > &v)
ADvari< Double > & atan2(const ADvari< Double > &L, const ADvari< Double > &R)
ADvari< Double > & sqrt(const ADvari< Double > &v)
bool operator==(const Handle< T > &h1, const Handle< T > &h2)
Compare two handles.
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)
Definition uninit.c:94