Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Sacado_trad2.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Sacado Package
5// Copyright (2007) 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// Extension of the RAD package (Reverse Automatic Differentiation) --
31// a package specialized for function and gradient evaluations -- to
32// Hessian-vector products.
33// Written in 2007 by David M. Gay at Sandia National Labs, Albuquerque, NM.
34
35#ifndef SACADO_TRAD2_H
36#define SACADO_TRAD2_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 be here when SACADO_NAMESPACE is #defined
67#endif
68#endif
69
70namespace Sacado {
71namespace Rad2 {
72
73// -DRAD_NO_USING_STDCC is needed, e.g., with Sun CC 5.7
74#ifndef RAD_NO_USING_STDCC
75 // Bring math functions into scope
76 using std::exp;
77 using std::log;
78 using std::log10;
79 using std::sqrt;
80 using std::cos;
81 using std::sin;
82 using std::tan;
83 using std::acos;
84 using std::asin;
85 using std::atan;
86 using std::cosh;
87 using std::sinh;
88 using std::tanh;
89 using std::abs;
90 using std::fabs;
91 using std::atan2;
92 using std::pow;
93#endif
94
95#ifdef RAD_AUTO_AD_Const
96#undef RAD_DEBUG_BLOCKKEEP
97#else
98#ifdef RAD_DEBUG_BLOCKKEEP
99#if !(RAD_DEBUG_BLOCKKEEP > 0)
100#undef RAD_DEBUG_BLOCKKEEP
101#else
102extern "C" void _uninit_f2c(void *x, int type, long len);
103
104template <typename T>
105struct UninitType {};
106
107template <>
108struct UninitType<float> {
109 static const int utype = 4;
110};
111
112template <>
113struct UninitType<double> {
114 static const int utype = 5;
115};
116
117template <typename T>
118struct UninitType< std::complex<T> > {
119 static const int utype = UninitType<T>::utype + 2;
120};
121
122#endif /*RAD_DEBUG_BLOCKKEEP > 0*/
123#endif /*RAD_DEBUG_BLOCKKEEP*/
124#endif /*RAD_AUTO_AD_Const*/
125
127
128 template<typename T> class
130 public:
131 typedef double dtype;
132 typedef T ttype;
133 };
134 template<> class
136 public:
139 };
140
141#define Dtype typename DoubleAvoid<Double>::dtype
142#define Ttype typename DoubleAvoid<Double>::ttype
143
144 template<typename Double> class IndepADvar;
145 template<typename Double> class ConstADvar;
146 template<typename Double> class ConstADvari;
147 template<typename Double> class ADvar;
148 template<typename Double> class ADvar1;
149 template<typename Double> class ADvar1g;
150 template<typename Double> class ADvar1s;
151 template<typename Double> class ADvar2;
152 template<typename Double> class ADvar2g;
153 template<typename Double> class ADvar2q;
154 template<typename Double> class ADvari;
155 template<typename Double> class ADvari_block;
156 template<typename Double> class ADvarn;
157 template<typename Double> class Derp;
158
159 template<typename Double> struct
160ADmemblock { // We get memory in ADmemblock chunks and never give it back,
161 // but reuse it once computations start anew after call(s) on
162 // ADcontext::Gradcomp() or ADcontext::Weighted_Gradcomp().
164 Double memblk[2000];
165 };
166
167 template<typename Double> class
169 public:
171 enum { Gulp = 1021 };
174 ADVari *pADvari[Gulp];
175 };
176
177 template<typename Double> class
178ADcontext { // A singleton class: one instance in radops.c
180 typedef ADvari <Double> ADVari;
181 typedef ADvar1 <Double> ADVar1;
182 typedef ADvar2 <Double> ADVar2;
183 typedef ADvarn <Double> ADVarn;
190 char *Mbase;
191 size_t Mleft;
192 ADVari **Ailimit, **Ainext;
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
205 int rad_busy_blocks;
206 ADMemblock *rad_Oldcurmb;
207#endif
208 public:
209 static const Double One, negOne;
210 ADcontext();
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*);
216 static inline void Gradcomp() { Gradcomp(1); }
217 inline void ADvari_record(ADVari *x) {
218 if (Ainext >= Ailimit)
219 new_ADvari_block();
220 *Ainext++ = x;
221 }
222 };
223
224 template<typename Double> class
225CADcontext: public ADcontext<Double> { // for possibly constant ADvar values
226 protected:
227 bool fpval_implies_const;
228 public:
229 friend class ADvar<Double>;
230 CADcontext(): ADcontext<Double>() { fpval_implies_const = false; }
231 };
232
233 template<typename Double> class
234Derp { // one derivative-propagation operation
235 public:
236 friend class ADvarn<Double>;
238 static Derp *LastDerp;
240 const Double *a;
241 const ADVari *b;
242 mutable ADVari *c;
243 Derp(){};
244 Derp(const ADVari *);
245 Derp(const Double *, const ADVari *);
246 Derp(const Double *, const ADVari *, const ADVari *);
247 /* c->aval += a * b->aval; */
248 };
249
250
251// Now we use #define to overcome bad design in the C++ templating system
252
253#define Ai const ADvari<Double>&
254#define AI const IndepADvar<Double>&
255#define T template<typename Double>
256#define D Double
257#define T1(f) \
258T F f (AI); \
259T F f (Ai);
260#define T2(r,f) \
261 T r f(Ai,Ai); \
262 T r f(Ai,D); \
263 T r f(Ai,Dtype); \
264 T r f(Ai,long); \
265 T r f(Ai,int); \
266 T r f(D,Ai); \
267 T r f(Dtype,Ai); \
268 T r f(long,Ai); \
269 T r f(int,Ai); \
270 T r f(AI,D); \
271 T r f(AI,Dtype); \
272 T r f(AI,long); \
273 T r f(AI,int); \
274 T r f(D,AI); \
275 T r f(Dtype,AI); \
276 T r f(long,AI); \
277 T r f(int,AI); \
278 T r f(Ai,AI);\
279 T r f(AI,Ai);\
280 T r f(AI,AI);
281
282#define F ADvari<Double>&
283T2(F, operator+)
284T2(F, operator-)
285T2(F, operator*)
286T2(F, operator/)
287T2(F, atan2)
288T2(F, pow)
289T2(F, max)
290T2(F, min)
291T2(int, operator<)
292T2(int, operator<=)
293T2(int, operator==)
294T2(int, operator!=)
295T2(int, operator>=)
296T2(int, operator>)
297T1(operator+)
298T1(operator-)
299T1(abs)
300T1(acos)
301T1(acosh)
302T1(asin)
303T1(asinh)
304T1(atan)
305T1(atanh)
306T1(cos)
307T1(cosh)
308T1(exp)
309T1(log)
310T1(log10)
311T1(sin)
312T1(sinh)
313T1(sqrt)
314T1(tan)
315T1(tanh)
316T1(fabs)
317
318T F copy(AI);
319T F copy(Ai);
320
321#undef F
322#undef T2
323#undef T1
324#undef D
325#undef T
326#undef AI
327#undef Ai
328
329} /* namespace Rad2 */
330} /* namespace Sacado */
331#define SNS Sacado::Rad2
332namespace std { // Moved here from bottom for use in testing nesting of Rad with itself.
333 using SNS::exp;
334 using SNS::log;
335 using SNS::log10;
336 using SNS::sqrt;
337 using SNS::cos;
338 using SNS::sin;
339 using SNS::tan;
340 using SNS::acos;
341 using SNS::asin;
342 using SNS::atan;
343 using SNS::cosh;
344 using SNS::sinh;
345 using SNS::tanh;
346 using SNS::abs;
347 using SNS::fabs;
348 using SNS::atan2;
349 using SNS::pow;
350}
351#undef SNS
352namespace Sacado {
353namespace Rad2 {
354
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,
358 template<typename Double>ADvari<Double>& ADfn(Double f, int n,
359 const IndepADvar<Double> *x, const Double *g);
360
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>&);
374
388
389 template<typename Double> ADvari<Double>&
390ADf1(Double f, Double g, Double h, const ADvari<Double> &x);
391
392 template<typename Double> ADvari<Double>&
393ADf2(Double f, Double gx, Double gy, Double hxx,
394 Double hxy, Double hyy, const ADvari<Double> &x, const ADvari<Double> &y);
395
396 template<typename Double> ADvari<Double>&
397ADfn(Double f, int n, const IndepADvar<Double> *x, const Double *g, const Double *h);
398
399 template<typename Double> class
400ADvari { // implementation of an ADvar
401 public:
402 typedef Double value_type;
405#ifdef RAD_AUTO_AD_Const
406 friend class IndepADvar<Double>;
407#ifdef RAD_Const_WARN
408 mutable const IndepADVar *padv;
409#else
410 protected:
411 mutable const IndepADVar *padv;
412#endif //RAD_Const_WARN
413 private:
414 ADvari *Next;
415 static ADvari *First_ADvari, **Last_ADvari;
416 public:
417 inline void ADvari_padv(const IndepADVar *v) {
418 *Last_ADvari = this;
419 Last_ADvari = &this->Next;
420 this->padv = v;
421 }
422#endif //RAD_AUTO_AD_Const
423#ifdef RAD_DEBUG
424 int gcgen;
425 int opno;
426 static int gcgen_cur, last_opno, zap_gcgen, zap_gcgen1, zap_opno;
427 static FILE *debug_file;
428#endif
431 Double Val; // result of this operation
432 mutable Double aval; // adjoint -- partial of final result w.r.t. this Val
433 mutable Double dO; // deriv of op w.r.t. t in x + t*p
434 mutable Double aO; // adjoint (in Hv computation) of op
435 mutable Double adO; // adjoint (in Hv computation) of dO
436 void *operator new(size_t len) {
437#ifdef RAD_DEBUG
438 ADvari *rv = (ADvari*)ADvari::adc.Memalloc(len);
439 rv->gcgen = gcgen_cur;
440 rv->opno = ++last_opno;
441 if (last_opno == zap_opno && gcgen_cur == zap_gcgen)
442 printf("");
443 return rv;
444#else
445 return ADvari::adc.Memalloc(len);
446#endif
447 }
448 void operator delete(void*) {} /*Should never be called.*/
449 inline ADvari(Advari_Opclass oc, Double t):
450 opclass(oc), Val(t), aval(0.), dO(0.)
451 { if (oc != Hv_const) ADvari::adc.ADvari_record(this); }
452 inline ADvari(Advari_Opclass oc, Double t, Double ta):
453 opclass(oc), Val(t), aval(ta), dO(0.)
454 { if (oc != Hv_const) ADvari::adc.ADvari_record(this); }
455 private:
456 inline ADvari(): Val(0.), aval(0.), dO(0.) {} // prevent construction without value (?)
457 public:
458 friend class ConstADvari<Double>;
459#ifdef RAD_AUTO_AD_Const
460 friend class ADcontext<Double>;
461 friend class ADvar<Double>;
462 friend class ADvar1<Double>;
463 friend class ADvar1s<Double>;
464 friend class ADvar2<Double>;
465 friend class ADvar2q<Double>;
466 friend class ConstADvar<Double>;
467 ADvari(const IndepADVar *, Double);
468 ADvari(const IndepADVar *, Double, Double);
469 ADvari(const IndepADVar *, Double, Double, int);
470#endif
471#define F friend
472#define R ADvari&
473#define Ai const ADvari&
474#define T1(r,f) F r f <>(Ai);
475#define T2(r,f) \
476F r f <>(Ai,Ai); \
477F r f <>(Ttype,Ai); \
478F r f <>(Ai,Ttype); \
479F r f <>(double,Ai); \
480F r f <>(Ai,double); \
481F r f <>(long,Ai); \
482F r f <>(Ai,long); \
483F r f <>(int,Ai); \
484F r f <>(Ai,int);
485 T1(R,operator+)
486 T2(R,operator+)
487 T1(R,operator-)
488 T2(R,operator-)
489 T2(R,operator*)
490 T2(R,operator/)
491 T1(R,abs)
492 T1(R,acos)
493 T1(R,acosh)
494 T1(R,asin)
495 T1(R,asinh)
496 T1(R,atan)
497 T1(R,atanh)
498 T2(R,atan2)
499 T2(R,max)
500 T2(R,min)
501 T1(R,cos)
502 T1(R,cosh)
503 T1(R,exp)
504 T1(R,log)
505 T1(R,log10)
506 T2(R,pow)
507 T1(R,sin)
508 T1(R,sinh)
509 T1(R,sqrt)
510 T1(R,tan)
511 T1(R,tanh)
512 T1(R,fabs)
513 T1(R,copy)
514 T2(int,operator<)
515 T2(int,operator<=)
516 T2(int,operator==)
517 T2(int,operator!=)
518 T2(int,operator>=)
519 T2(int,operator>)
520#undef T2
521#undef T1
522#undef Ai
523#undef R
524#undef F
525
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);
531
532 inline operator Double() { return this->Val; }
533 inline operator Double() const { return this->Val; }
534 };
535
536 template<typename Double> class
537ADvar1: public ADvari<Double> { // simplest unary ops
538 public:
541 ADvar1(Advari_Opclass oc, Double val1, const Double *a1, const ADVari *c1):
542 ADVari(oc,val1), d(a1,this,c1) {}
543#ifdef RAD_AUTO_AD_Const
544 typedef typename ADVari::IndepADVar IndepADVar;
545 typedef ADvar<Double> ADVar;
546 ADvar1(const IndepADVar*, const IndepADVar&);
547 ADvar1(const IndepADVar*, const ADVari&);
548 ADvar1(Advari_Opclass oc, const Double val1, const Double *a1,
549 const ADVari *c1, const ADVar *v):
550 ADVari(oc,val1), d(a1,this,c1) {
551 c1->padv = 0;
552 this->ADvari_padv(v);
553 }
554#endif
555 };
556
557
558 template<typename Double> class
560 private:
562 ConstADvari() {}; // prevent construction without value (?)
564 public:
568 inline void *operator new(size_t len) { return ConstADvari::cadc.Memalloc(len); }
569 inline ConstADvari(Double t): ADVari(Hv_copy, t) { prevcad = lastcad; lastcad = this; }
570 static void aval_reset(void);
571 };
572
573
574 template<typename Double> class
575IndepADvar { // an independent ADvar
576 private:
578 /* private to prevent assignment */
579#ifdef RAD_AUTO_AD_Const
580 if (cv)
581 cv->padv = 0;
582 return new ADvar1<Double>(this,x);
583#elif defined(RAD_EQ_ALIAS)
584 this->cv = x.cv;
585 return *this;
586#else
587 return ADvar_operatoreq(this,*x.cv);
588#endif //RAD_AUTO_AD_Const
589 }
590 protected:
591 static void AD_Const(const IndepADvar&);
593 public:
594 typedef Double value_type;
595 friend class ADvar<Double>;
596 friend class ADcontext<Double>;
597 friend class ADvar1<Double>;
598 friend class ADvarn<Double>;
602 IndepADvar(double);
603 IndepADvar(int);
604 IndepADvar(long);
605 IndepADvar& operator= (Double);
606#ifdef RAD_AUTO_AD_Const
607 friend IndepADvar& ADvar_operatoreq<>(IndepADvar*, const ADVari&);
608 inline IndepADvar(const IndepADvar &x) { cv = x.cv ? new ADvar1<Double>(this, x) : 0; };
609 inline IndepADvar(ADVari *ncv) { cv = ncv; }
610 inline IndepADvar() { cv = 0; }
611 inline ~IndepADvar() {
612 if (cv)
613 cv->padv = 0;
614 }
615#else
616 inline IndepADvar() {
617#ifndef RAD_EQ_ALIAS
618 cv = 0;
619#endif
620 }
621 inline ~IndepADvar() {}
622 friend IndepADvar& ADvar_operatoreq<>(IndepADvar*, const ADVari&);
623#endif
624
625#ifdef RAD_Const_WARN
626 inline operator ADVari&() const {
627 ADVari *tcv = this->cv;
628 if (tcv->opno < 0)
629 RAD_Const_Warn(tcv);
630 return *tcv;
631 }
632 inline operator ADVari*() const {
633 ADVari *tcv = this->cv;
634 if (tcv->opno < 0)
635 RAD_Const_Warn(tcv);
636 return tcv;
637 }
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
642
643 Double val() const { return cv->Val; }
644 Double adj() const { return cv->aval; }
645
646 friend void AD_Const1<>(Double*, const IndepADvar&);
647
648 friend ADVari& ADf1<>(Double, Double, const IndepADvar&);
649 friend ADVari& ADf2<>(Double, Double, Double, const IndepADvar&, const IndepADvar&);
650 friend ADVari& ADf2<>(Double, Double, Double, const ADVari&, const IndepADvar&);
651 friend ADVari& ADf2<>(Double, Double, Double, const IndepADvar&, const ADVari&);
652
653 static inline void Gradcomp(int wantgrad)
654 { ADcontext<Double>::Gradcomp(wantgrad); }
655 static inline void Gradcomp()
657 static inline void Hvprod(int n, ADVar **vp, Double *v, Double *hv)
658 { ADcontext<Double>::Hvprod(n, vp, v, hv); }
659 static inline void aval_reset() { ConstADvari<Double>::aval_reset(); }
660 static inline void Weighted_Gradcomp(int n, ADVar **v, Double *w)
662
663 /* We use #define to deal with bizarre templating rules that apparently */
664 /* require us to spell the some conversion explicitly */
665
666
667#define Ai const ADVari&
668#define AI const IndepADvar&
669#define D Double
670#define T2(r,f) \
671 r f <>(AI,AI);\
672 r f <>(Ai,AI);\
673 r f <>(AI,Ai);\
674 r f <>(Ttype,AI);\
675 r f <>(double,AI);\
676 r f <>(long,AI);\
677 r f <>(int,AI);\
678 r f <>(AI,Ttype);\
679 r f <>(AI,double);\
680 r f <>(AI,long);\
681 r f <>(AI,int);
682#define T1(f) friend ADVari& f<> (AI);
683
684#define F friend ADVari&
685T2(F, operator+)
686T2(F, operator-)
687T2(F, operator*)
688T2(F, operator/)
689T2(F, atan2)
690T2(F, max)
691T2(F, min)
692T2(F, pow)
693#undef F
694#define F friend int
695T2(F, operator<)
696T2(F, operator<=)
697T2(F, operator==)
698T2(F, operator!=)
699T2(F, operator>=)
700T2(F, operator>)
701
702T1(operator+)
703T1(operator-)
704T1(abs)
705T1(acos)
706T1(acosh)
707T1(asin)
708T1(asinh)
709T1(atan)
710T1(atanh)
711T1(cos)
712T1(cosh)
713T1(exp)
714T1(log)
715T1(log10)
716T1(sin)
717T1(sinh)
718T1(sqrt)
719T1(tan)
720T1(tanh)
721T1(fabs)
722T1(copy)
723
724#undef F
725#undef T1
726#undef T2
727#undef D
728#undef AI
729#undef Ai
730
731 };
732
733 template<typename Double> class
734ADvar: public IndepADvar<Double> { // an "active" variable
735 public:
737 template <typename U> struct apply { typedef ADvar<U> type; };
738
742 private:
743 void ADvar_ctr(Double d) {
744 ADVari *x;
745 if (ConstADVari::cadc.fpval_implies_const)
746 x = new ConstADVari(d);
747 else {
748#ifdef RAD_AUTO_AD_Const
749 x = new ADVari((IndepADVar*)this, d);
750 x->ADvari_padv(this);
751#else
752 x = new ADVari(Hv_const, d);
753#endif
754 }
755 this->cv = x;
756 }
757 public:
758 friend class ADvar1<Double>;
760 ADvar() { /* cv = 0; */ }
761 ADvar(Ttype d) { ADvar_ctr(d); }
762 ADvar(double i) { ADvar_ctr(Double(i)); }
763 ADvar(int i) { ADvar_ctr(Double(i)); }
764 ADvar(long i) { ADvar_ctr(Double(i)); }
765 inline ~ADvar() {}
766#ifdef RAD_AUTO_AD_Const
767 ADvar(IndepADVar &x) {
768 this->cv = x.cv ? new ADVar1(this, x) : 0;
769 }
770 ADvar(ADVari &x) :IndepADVar(&x) { x.ADvari_padv((IndepADVar*)this);}
771 inline ADvar& operator=(IndepADVar &x) {
772 if (this->cv)
773 this->cv->padv = 0;
774 this->cv = new ADVar1(this,x);
775 return *this;
776 }
777 inline ADvar& operator=(ADVari &x) {
778 if (this->cv)
779 this->cv->padv = 0;
780 this->cv = new ADVar1(this, x);
781 return *this;
782 }
783#else
784 friend ADvar& ADvar_operatoreq<>(ADvar*, const ADVari&);
785#ifdef RAD_EQ_ALIAS
786 /* allow aliasing v and w after "v = w;" */
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; }
791#else
792 friend ADvar& ADvar_operatoreq<>(ADvar*, const ADVari&);
793 ADvar(const IndepADVar &x) { this->cv = x.cv ? new ADVar1(x.cv->Val, &this->cv->adc.One, x.cv) : 0; }
794 ADvar(const ADvar&x) { this->cv = x.cv ?
795 new ADVar1(Hv_copy, x.cv->Val, &this->cv->adc.One, (ADVari*)x.cv) : 0; }
796 ADvar(const ADVari &x) { this->cv = new ADVar1(Hv_copy, x.Val, &this->cv->adc.One, &x); }
797 inline ADvar& operator=(IndepADVar &x) { return ADvar_operatoreq(this,*x.cv); };
798 inline ADvar& operator=(const ADVari &x) { return ADvar_operatoreq(this,x); };
799#endif /* RAD_EQ_ALIAS */
800#endif /* RAD_AUTO_AD_Const */
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);
810 inline static bool get_fpval_implies_const(void)
811 { return ConstADVari::cadc.fpval_implies_const; }
812 inline static void set_fpval_implies_const(bool newval)
813 { ConstADVari::cadc.fpval_implies_const = newval; }
814 inline static bool setget_fpval_implies_const(bool newval) {
815 bool oldval = ConstADVari::cadc.fpval_implies_const;
816 ConstADVari::cadc.fpval_implies_const = newval;
817 return oldval;
818 }
819 static inline void Gradcomp(int wantgrad)
820 { ADcontext<Double>::Gradcomp(wantgrad); }
821 static inline void Gradcomp()
823 static inline void aval_reset() { ConstADVari::aval_reset(); }
824 static inline void Weighted_Gradcomp(int n, ADvar **v, Double *w)
826 };
827
828template<typename Double>
829 inline void AD_Const1(Double *notused, const IndepADvar<Double>&v)
831
832template<typename Double>
833 inline void AD_Const(const IndepADvar<Double>&v) { AD_Const1((Double*)0, v); }
834
835 template<typename Double> class
837 public:
839 typedef typename ADVar::ADVari ADVari;
843 private: // disable op=
852 void ConstADvar_ctr(Double);
853 public:
854 ConstADvar(Ttype d) { ConstADvar_ctr(d); }
855 ConstADvar(double i) { ConstADvar_ctr(Double(i)); }
856 ConstADvar(int i) { ConstADvar_ctr(Double(i)); }
857 ConstADvar(long i) { ConstADvar_ctr(Double(i)); }
858 ConstADvar(const IndepADVar&);
859 ConstADvar(const ConstADvar&);
860 ConstADvar(const ADVari&);
861 inline ~ConstADvar() {}
862#ifdef RAD_NO_CONST_UPDATE
863 private:
864#endif
865 ConstADvar();
866 inline ConstADvar& operator=(Double d) { this->cv->Val = d; return *this; }
867 inline ConstADvar& operator=(ADVari& d) { this->cv->Val = d.Val; return *this; }
868 };
869
870 template<typename Double> class
871ADvar1s: public ADvar1<Double> { // unary ops with partials
872 public:
874 typedef typename ADVar1::ADVari ADVari;
875 Double pL; // deriv of op w.r.t. left operand L
876 ADvar1s(Double val1, Double a1, const ADVari *c1):
877 ADVar1(Hv_timesL,val1,&pL,c1), pL(a1) {}
878#ifdef RAD_AUTO_AD_Const
879 Double a;
880 typedef typename ADVar1::ADVar ADVar;
881 ADvar1s(Double val1, Double a1, const ADVari *c1, const ADVar *v):
882 ADVar1(Hv_timesL,val1,&a,c1,v), a(a1) {}
883#endif
884 };
885
886 template<typename Double> class
887ADvar1g: public ADvar1<Double> { // unary ops with partials
888 public:
890 typedef typename ADVar1::ADVari ADVari;
891 Double pL; // deriv of op w.r.t. left operand L
892 Double pL2; // partial of op w.r.t. L,L
893 ADvar1g(Double val1, Double d1, Double d2, const ADVari *c1):
894 ADVar1(Hv_unary,val1,&pL,c1), pL(d1), pL2(d2) {}
895 };
896
897 template<typename Double> class
898ADvar2: public ADvari<Double> { // basic binary ops
899 public:
902 DErp dL, dR;
903 ADvar2(Advari_Opclass oc, Double val1, const ADVari *Lcv, const Double *Lc,
904 const ADVari *Rcv, const Double *Rc):
905 ADVari(oc,val1) {
906 dR.next = DErp::LastDerp;
907 dL.next = &dR;
908 DErp::LastDerp = &dL;
909 dL.a = Lc;
910 dL.c = (ADVari*)Lcv;
911 dR.a = Rc;
912 dR.c = (ADVari*)Rcv;
913 dL.b = dR.b = this;
914 }
915#ifdef RAD_AUTO_AD_Const
916 typedef ADvar<Double> ADVar;
917 ADvar2(Advari_Opclass oc, Double val1, const ADVari *Lcv, const Double *Lc,
918 const ADVari *Rcv, const Double *Rc, ADVar *v):
919 ADVari(oc,val1) {
920 dR.next = DErp::LastDerp;
921 dL.next = &dR;
922 DErp::LastDerp = &dL;
923 dL.a = Lc;
924 dL.c = Lcv;
925 dR.a = Rc;
926 dR.c = Rcv;
927 dL.b = dR.b = this;
928 Lcv->padv = 0;
929 this->ADvari_padv(v);
930 }
931#endif
932 };
933
934 template<typename Double> class
935ADvar2q: public ADvar2<Double> { // binary ops with partials
936 public:
938 typedef typename ADVar2::ADVari ADVari;
939 typedef typename ADVar2::DErp DErp;
940 Double pL; // deriv of op w.r.t. left operand L
941 Double pR; // deriv of op w.r.t. right operand R
942 Double pLR; // second partial w.r.t. L,R
943 Double pR2; // second partial w.r.t. R,R
944 ADvar2q(Double val1, Double Lp, Double Rp, Double LR, Double R2,
945 const ADVari *Lcv, const ADVari *Rcv):
946 ADVar2(Hv_quotLR,val1,Lcv,&pL,Rcv,&pR),
947 pL(Lp), pR(Rp), pLR(LR), pR2(R2) {}
948#ifdef RAD_AUTO_AD_Const
949 typedef typename ADVar2::ADVar ADVar;
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) {}
954#endif
955 };
956
957 template<typename Double> class
958ADvar2g: public ADvar2<Double> { // general binary ops with partials
959 public:
961 typedef typename ADVar2::ADVari ADVari;
962 Double pL; // deriv of op w.r.t. left operand L
963 Double pR; // deriv of op w.r.t. right operand R
964 Double pL2; // second partial w.r.t. L,L
965 Double pLR; // second partial w.r.t. L,R
966 Double pR2; // second partial w.r.t. R,R
967 ADvar2g(Double val1, Double Lp, Double Rp, Double L2, Double LR, Double R2,
968 const ADVari *Lcv, const ADVari *Rcv):
969 ADVar2(Hv_binary,val1,Lcv,&pL,Rcv,&pR),
970 pL(Lp), pR(Rp), pL2(L2), pLR(LR), pR2(R2) { }
971 };
972
973 template<typename Double> class
974ADvarn: public ADvari<Double> { // n-ary ops with partials g and
975 // 2nd partials h (lower triangle, rowwise)
976 public:
981 int n;
982 Double *G, *H;
984 ADvarn(Double val1, int n1, const IndepADVar *x, const Double *g, const Double *h):
985 ADVari(Hv_nary,val1), n(n1) {
986 DErp *d1, *dlast;
987 Double *a1;
988 int i, nh;
989
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++) {
994 d1->next = dlast;
995 dlast = d1;
996 a1[i] = g[i];
997 d1->a = &a1[i];
998 d1->b = this;
999 d1->c = x[i].cv;
1000 }
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++)
1005 a1[i] = h[i];
1006 }
1007 };
1008
1009template<typename Double>
1010 inline ADvari<Double>& operator+(const ADvari<Double> &T) { return *(ADvari<Double>*)&T; }
1011
1012template<typename Double>
1013 inline int operator<(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val < R.Val; }
1014template<typename Double>
1015 inline int operator<(const ADvari<Double> &L, Double R) { return L.Val < R; }
1016template<typename Double>
1017 inline int operator<(Double L, const ADvari<Double> &R) { return L < R.Val; }
1018
1019template<typename Double>
1020 inline int operator<=(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val <= R.Val; }
1021template<typename Double>
1022 inline int operator<=(const ADvari<Double> &L, Double R) { return L.Val <= R; }
1023template<typename Double>
1024 inline int operator<=(Double L, const ADvari<Double> &R) { return L <= R.Val; }
1025
1026template<typename Double>
1027 inline int operator==(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val == R.Val; }
1028template<typename Double>
1029 inline int operator==(const ADvari<Double> &L, Double R) { return L.Val == R; }
1030template<typename Double>
1031 inline int operator==(Double L, const ADvari<Double> &R) { return L == R.Val; }
1032
1033template<typename Double>
1034 inline int operator!=(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val != R.Val; }
1035template<typename Double>
1036 inline int operator!=(const ADvari<Double> &L, Double R) { return L.Val != R; }
1037template<typename Double>
1038 inline int operator!=(Double L, const ADvari<Double> &R) { return L != R.Val; }
1039
1040template<typename Double>
1041 inline int operator>=(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val >= R.Val; }
1042template<typename Double>
1043 inline int operator>=(const ADvari<Double> &L, Double R) { return L.Val >= R; }
1044template<typename Double>
1045 inline int operator>=(Double L, const ADvari<Double> &R) { return L >= R.Val; }
1046
1047template<typename Double>
1048 inline int operator>(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val > R.Val; }
1049template<typename Double>
1050 inline int operator>(const ADvari<Double> &L, Double R) { return L.Val > R; }
1051template<typename Double>
1052 inline int operator>(Double L, const ADvari<Double> &R) { return L > R.Val; }
1053
1054template<typename Double>
1055 inline void *ADcontext<Double>::Memalloc(size_t len) {
1056 if (Mleft >= len)
1057 return Mbase + (Mleft -= len);
1058 return new_ADmemblock(len);
1059 }
1060
1061template<typename Double>
1062 inline Derp<Double>::Derp(const ADVari *c1): c((ADVari*)c1) {
1063 next = LastDerp;
1064 LastDerp = this;
1065 }
1066
1067template<typename Double>
1068 inline Derp<Double>::Derp(const Double *a1, const ADVari *c1): a(a1), c((ADVari*)c1) {
1069 next = LastDerp;
1070 LastDerp = this;
1071 }
1072
1073template<typename Double>
1074 inline Derp<Double>::Derp(const Double *a1, const ADVari *b1, const ADVari *c1):
1075 a(a1), b(b1), c((ADVari*)c1) {
1076 next = LastDerp;
1077 LastDerp = this;
1078 }
1079
1080/**** radops ****/
1081
1082template<typename Double> Derp<Double> *Derp<Double>::LastDerp = 0;
1083template<typename Double> ADcontext<Double> ADvari<Double>::adc;
1084template<typename Double> const Double ADcontext<Double>::One = 1.;
1085template<typename Double> const Double ADcontext<Double>::negOne = -1.;
1086template<typename Double> CADcontext<Double> ConstADvari<Double>::cadc;
1087template<typename Double> ConstADvari<Double> *ConstADvari<Double>::lastcad;
1088
1089#ifdef RAD_AUTO_AD_Const
1090template<typename Double> ADvari<Double>* ADvari<Double>::First_ADvari;
1092#endif
1093
1094#ifdef RAD_DEBUG
1095#ifndef RAD_DEBUG_gcgen1
1096#define RAD_DEBUG_gcgen1 -1
1097#endif
1098template<typename Double> int ADvari<Double>::gcgen_cur;
1099template<typename Double> int ADvari<Double>::last_opno;
1100template<typename Double> int ADvari<Double>::zap_gcgen;
1101template<typename Double> int ADvari<Double>::zap_gcgen1 = RAD_DEBUG_gcgen1;
1102template<typename Double> int ADvari<Double>::zap_opno;
1103template<typename Double> FILE *ADvari<Double>::debug_file;
1104#endif
1105
1106
1107 template<typename Double>
1109{
1110 ADVari_block *fb;
1111
1112 First = (ADMemblock*)First0;
1113 First->next = 0;
1114 Busy = First;
1115 Free = 0;
1116 Mbase = (char*)First->memblk;
1117 Mleft = sizeof(First->memblk);
1118 AiFirst = Aibusy = fb = (ADVari_block*)First1;
1119 Aifree = 0;
1120 Ainext = fb->pADvari;
1121 fb->next = fb->prev = 0;
1122 fb->limit = Ailimit = fb->pADvari + ADVari_block::Gulp;
1123 rad_need_reinit = 0;
1124#ifdef RAD_DEBUG_BLOCKKEEP
1125 rad_busy_blocks = 0;
1126 rad_mleft_save = 0;
1127 rad_Oldcurmb = 0;
1128#endif
1129 }
1130
1131template<typename Double> void*
1133{
1134 ADMemblock *mb, *mb0, *mb1, *mbf, *x;
1135 ADVari_block *b;
1136#ifdef RAD_AUTO_AD_Const
1137 ADVari *a, *anext;
1139#ifdef RAD_Const_WARN
1140 ADVari *cv;
1141 int i, j;
1142#endif
1143#endif /*RAD_AUTO_AD_Const*/
1144
1145 if ((rad_need_reinit & 1) && this == &ADVari::adc) {
1146 rad_need_reinit &= ~1;
1147 DErp::LastDerp = 0;
1148 Aibusy = b = AiFirst;
1149 Aifree = b->next;
1150 b->next = b->prev = 0;
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))
1155 _uninit_f2c(Mbase + Mleft,
1156 UninitType<Double>::utype,
1157 (sizeof(First->memblk) - Mleft)
1158 /sizeof(typename Sacado::ValueType<Double>::type));
1159 if ((mb = Busy->next)) {
1160 if (!(mb0 = rad_Oldcurmb))
1161 mb0 = (ADMemblock*)First0;
1162 for(;; mb = mb->next) {
1163 _uninit_f2c(mb->memblk,
1164 UninitType<Double>::utype,
1165 sizeof(First->memblk)
1166 /sizeof(typename Sacado::ValueType<Double>::type));
1167 if (mb == mb0)
1168 break;
1169 }
1170 }
1171 rad_Oldcurmb = Busy;
1172 if (rad_busy_blocks >= RAD_DEBUG_BLOCKKEEP) {
1173 rad_busy_blocks = 0;
1174 rad_Oldcurmb = 0;
1175 mb0 = (ADMemblock*)First0;
1176 mbf = Free;
1177 for(mb = Busy; mb != mb0; mb = mb1) {
1178 mb1 = mb->next;
1179 mb->next = mbf;
1180 mbf = mb;
1181 }
1182 Free = mbf;
1183 Busy = mb;
1184 Mbase = (char*)First->memblk;
1185 Mleft = sizeof(First->memblk);
1186 }
1187
1188#else /* !RAD_DEBUG_BLOCKKEEP */
1189
1190 mb0 = First;
1191 mbf = Free;
1192 for(mb = Busy; mb != mb0; mb = mb1) {
1193 mb1 = mb->next;
1194 mb->next = mbf;
1195 mbf = mb;
1196 }
1197 Free = mbf;
1198 Busy = mb;
1199 Mbase = (char*)First->memblk;
1200 Mleft = sizeof(First->memblk);
1201#ifdef RAD_AUTO_AD_Const
1202 if (ADVari::adc.rad_need_reinit & 2) {
1203 ADVari::adc.rad_need_reinit &= ~2;
1204 *ADVari::Last_ADvari = 0;
1205 ADVari::Last_ADvari = &ADVari::First_ADvari;
1206 anext = ADVari::First_ADvari;
1207 if (anext) {
1208 while((a = anext)) {
1209 anext = a->Next;
1210 if ((v = (IndepADvar<Double> *)a->padv)) {
1211#ifdef RAD_Const_WARN
1212 if ((i = a->opno) > 0)
1213 i = -i;
1214 j = a->gcgen;
1215 v->cv = cv = new ADVari(v, a->Val, a->aval);
1216 cv->opno = i;
1217 cv->gcgen = j;
1218#else
1219 v->cv = new ADVari(v, a->Val, a->aval);
1220#endif
1221 }
1222 }
1223 }
1224 }
1225#endif /*RAD_AUTO_AD_Const*/
1226#endif /*RAD_DEBUG_BLOCKKEEP*/
1227 if (Mleft >= len)
1228 return Mbase + (Mleft -= len);
1229 }
1230
1231 if ((x = Free))
1232 Free = x->next;
1233 else
1234 x = new ADMemblock;
1235#ifdef RAD_DEBUG_BLOCKKEEP
1236 rad_busy_blocks++;
1237#endif
1238 x->next = Busy;
1239 Busy = x;
1240 return (Mbase = (char*)x->memblk) +
1241 (Mleft = sizeof(First->memblk) - len);
1242 }
1243
1244template<typename Double> void
1246{
1247 ADVari_block *ob, *nb;
1248 ob = Aibusy;
1249 ob->limit = Ailimit; // should be unnecessary, but harmless
1250 if ((nb = Aifree))
1251 Aifree = nb->next;
1252 else
1253 nb = new ADVari_block;
1254 Aibusy = Aibusy->next = nb;
1255 nb->limit = Ailimit = (Ainext = nb->pADvari) + ADVari_block::Gulp;
1256 ob->next = nb;
1257 nb->prev = ob;
1258 nb->next = 0;
1259 }
1260
1261template<typename Double> void
1263{
1264 DErp *d;
1265
1266 if (ADVari::adc.rad_need_reinit) {
1267 for(d = DErp::LastDerp; d; d = d->next)
1268 d->c->aval = 0;
1269 }
1270 if (!(ADVari::adc.rad_need_reinit & 1)) {
1271 ADVari::adc.rad_need_reinit = 1;
1272 ADVari::adc.rad_mleft_save = ADVari::adc.Mleft;
1273 ADVari::adc.Mleft = 0;
1274 }
1275#ifdef RAD_DEBUG
1276 if (ADVari::gcgen_cur == ADVari::zap_gcgen1) {
1277 const char *fname;
1278 if (!(fname = getenv("RAD_DEBUG_FILE")))
1279 fname = "rad_debug.out";
1280 else if (!*fname)
1281 fname = 0;
1282 if (fname)
1283 ADVari::debug_file = fopen(fname, "w");
1284 ADVari::zap_gcgen1 = -1;
1285 }
1286#endif
1287 if ((d = DErp::LastDerp) != 0) {
1288#ifdef RAD_AUTO_AD_Const
1289 ADVari::adc.rad_need_reinit |= 2;
1290#endif /*RAD_AUTO_AD_Const*/
1291 if (wantgrad) {
1292 d->b->aval = 1;
1293#ifdef RAD_DEBUG
1294 if (ADVari::debug_file)
1295 do {
1296 fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
1297 d->c->opno, d->b->opno, d->c->aval,
1298 *d->a, d->b->aval);
1299 d->c->aval += *d->a * d->b->aval;
1300 fprintf(ADVari::debug_file, " = %g\n", d->c->aval);
1301 fflush(ADVari::debug_file);
1302 } while((d = d->next));
1303 else
1304#endif
1305 do d->c->aval += *d->a * d->b->aval;
1306 while((d = d->next));
1307 DErp::LastDerp = 0;
1308 }
1309 }
1310#ifdef RAD_DEBUG
1311 if (ADVari::debug_file) {
1312 fclose(ADVari::debug_file);
1313 ADVari::debug_file = 0;
1314 }
1315#endif //RAD_DEBUG
1316#ifdef RAD_DEBUG
1317 ADVari::gcgen_cur++;
1318 ADVari::last_opno = 0;
1319#endif
1320 }
1321
1322template<typename Double> void
1324{
1325 DErp *d;
1326 int i;
1327#ifdef RAD_Const_WARN
1328 ADVari *cv;
1329 int j;
1330#endif
1331#ifdef RAD_AUTO_AD_Const
1332 ADVari *a, *anext;
1334#endif /*RAD_AUTO_AD_Const*/
1335
1336 if (ADVari::adc.rad_need_reinit) {
1337 for(d = DErp::LastDerp; d; d = d->next)
1338 d->c->aval = 0;
1339 }
1340 if (!(ADVari::adc.rad_need_reinit & 1)) {
1341 ADVari::adc.rad_need_reinit = 1;
1342 ADVari::adc.rad_mleft_save = ADVari::adc.Mleft;
1343 ADVari::adc.Mleft = 0;
1344 }
1345#ifdef RAD_DEBUG
1346 if (ADVari::gcgen_cur == ADVari::zap_gcgen1) {
1347 const char *fname;
1348 if (!(fname = getenv("RAD_DEBUG_FILE")))
1349 fname = "rad_debug.out";
1350 else if (!*fname)
1351 fname = 0;
1352 if (fname)
1353 ADVari::debug_file = fopen(fname, "w");
1354 ADVari::zap_gcgen1 = -1;
1355 }
1356#endif
1357 if ((d = DErp::LastDerp) != 0) {
1358 for(i = 0; i < n; i++)
1359 V[i]->cv->aval = w[i];
1360#ifdef RAD_DEBUG
1361 if (ADVari::debug_file)
1362 do {
1363 fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
1364 d->c->opno, d->b->opno, d->c->aval, *d->a, d->b->aval);
1365 d->c->aval += *d->a * d->b->aval;
1366 fprintf(ADVari::debug_file, " = %g\n", d->c->aval);
1367 fflush(ADVari::debug_file);
1368 } while((d = d->next));
1369 else
1370#endif
1371 do d->c->aval += *d->a * d->b->aval;
1372 while((d = d->next));
1373 }
1374#ifdef RAD_DEBUG
1375 if (ADVari::debug_file) {
1376 fclose(ADVari::debug_file);
1377 ADVari::debug_file = 0;
1378 }
1379#endif //RAD_DEBUG
1380#ifdef RAD_AUTO_AD_Const
1381 *ADVari::Last_ADvari = 0;
1382 ADVari::Last_ADvari = &ADVari::First_ADvari;
1383 if ((anext = ADVari::First_ADvari) && !(ADVari::adc.rad_need_reinit & 2)) {
1384 ADVari::adc.rad_need_reinit = 3;
1385 while((a = anext)) {
1386 anext = a->Next;
1387 if ((v = (IndepADvar<Double> *)a->padv)) {
1388#ifdef RAD_Const_WARN
1389 if ((i = a->opno) > 0)
1390 i = -i;
1391 j = a->gcgen;
1392 v->cv = cv = new ADVari(v, a->Val, a->aval);
1393 cv->opno = i;
1394 cv->gcgen = j;
1395#else
1396 v->cv = new ADVari(v, a->Val, a->aval);
1397#endif
1398 }
1399 }
1400 DErp::LastDerp = 0;
1401 }
1402#endif /*RAD_AUTO_AD_Const*/
1403#ifdef RAD_DEBUG
1404 ADVari::gcgen_cur++;
1405 ADVari::last_opno = 0;
1406#endif
1407 }
1408
1409 template<typename Double>
1411{
1412
1413 ADVari *x = new ADVari(Hv_const, d);
1414 cv = x;
1415 }
1416
1417 template<typename Double>
1419{
1420
1421 ADVari *x = new ADVari(Hv_const, Double(i));
1422 cv = x;
1423 }
1424
1425 template<typename Double>
1427{
1428
1429 ADVari *x = new ADVari(Hv_const, Double(i));
1430 cv = x;
1431 }
1432
1433 template<typename Double>
1435{
1436
1437 ADVari *x = new ADVari(Hv_const, Double(i));
1438 cv = x;
1439 }
1440
1441 template<typename Double>
1443{
1444 ConstADVari *x = new ConstADVari(0.);
1445 this->cv = x;
1446 }
1447
1448 template<typename Double> void
1450{
1451 ConstADVari *x = new ConstADVari(d);
1452 this->cv = x;
1453 }
1454
1455 template<typename Double>
1457{
1458 ConstADVari *y = new ConstADVari(x.cv->Val);
1459 DErp *d = new DErp(&x.adc.One, y, x.cv);
1460 this->cv = y;
1461 }
1462
1463 template<typename Double>
1465{
1466 ConstADVari *y = new ConstADVari(x.cv->Val);
1467 DErp *d = new DErp(&x.cv->adc.One, y, (ADVari*)x.cv);
1468 this->cv = y;
1469 }
1470
1471 template<typename Double>
1473{
1474 ConstADVari *y = new ConstADVari(x.Val);
1475 DErp *d = new DErp(&x.adc.One, y, &x);
1476 this->cv = y;
1477 }
1478
1479 template<typename Double>
1480 void
1482{
1483 typedef ConstADvari<Double> ConstADVari;
1484
1485 ConstADVari *ncv = new ConstADVari(v.val());
1486#ifdef RAD_AUTO_AD_Const
1487 if (v.cv)
1488 v.cv->padv = 0;
1489#endif
1490 v.cv = ncv;
1491 }
1492
1493 template<typename Double>
1494 void
1496{
1498 while(x) {
1499 x->aval = 0;
1500 x = x->prevcad;
1501 }
1502 }
1503
1504#ifdef RAD_AUTO_AD_Const
1505
1506 template<typename Double>
1507ADvari<Double>::ADvari(const IndepADVar *x, Double d): Val(d), aval(0.)
1508{
1509 opclass = Hv_const;
1510 this->ADvari_padv(x);
1511 }
1512
1513 template<typename Double>
1514ADvari<Double>::ADvari(const IndepADVar *x, Double d, Double g): Val(d), aval(g)
1515{
1516 opclass = Hv_const;
1517 this->ADvari_padv(x);
1518 }
1519
1520 template<typename Double>
1521ADvar1<Double>::ADvar1(const IndepADVar *x, const IndepADVar &y):
1522 ADVari(Hv_copy, y.cv->Val), d((const Double*)&ADcontext<Double>::One, (ADVari*)this, y.cv)
1523{
1524 this->ADvari_padv(x);
1525 }
1526
1527 template<typename Double>
1528ADvar1<Double>::ADvar1(const IndepADVar *x, const ADVari &y):
1529 ADVari(Hv_copy, y.Val), d((const Double*)&ADcontext<Double>::One, this, &y)
1530{
1531 this->ADvari_padv(x);
1532 }
1533
1534#endif /* RAD_AUTO_AD_Const */
1535
1536 template<typename Double>
1537 IndepADvar<Double>&
1539{ This->cv = new ADvar1<Double>(Hv_copy, x.Val, &x.adc.One, &x);
1540 return *(IndepADvar<Double>*) This; }
1541
1542 template<typename Double>
1545{ This->cv = new ADvar1<Double>(Hv_copy, x.Val, &x.adc.One, &x); return *(ADvar<Double>*) This; }
1546
1547 template<typename Double>
1550{
1551#ifdef RAD_AUTO_AD_Const
1552 ADVari *ncv = new ADVari(this, d);
1553 if (this->cv)
1554 this->cv->padv = 0;
1555 this->cv = ncv;
1556#else
1557 this->cv = new ADVari(Hv_const, d);
1558#endif
1559 return *this;
1560 }
1561
1562 template<typename Double>
1565{
1566#ifdef RAD_AUTO_AD_Const
1567 ADVari *nv = new ADVari(this, d);
1568 if (this->cv)
1569 this->cv->padv = 0;
1570 this->cv = nv;
1571#else
1572 this->cv = ConstADVari::cadc.fpval_implies_const
1573 ? new ConstADVari(d)
1574 : new ADVari(Hv_const, d);
1575#endif
1576 return *this;
1577 }
1578
1579 template<typename Double>
1582 return *(new ADvar1<Double>(Hv_negate, -T.Val, &T.adc.negOne, &T));
1583 }
1584
1585 template<typename Double>
1586 ADvari<Double>&
1588 return *(new ADvar2<Double>(Hv_plusLR, L.Val + R.Val, &L, &L.adc.One, &R, &L.adc.One));
1589 }
1590
1591#ifdef RAD_AUTO_AD_Const
1592#define RAD_ACA ,this
1593#else
1594#define RAD_ACA /*nothing*/
1595#endif
1596
1597 template<typename Double>
1598 ADvar<Double>&
1600 ADVari *Lcv = this->cv;
1601 this->cv = new ADvar2<Double>(Hv_plusLR, Lcv->Val + R.Val, Lcv,
1602 &R.adc.One, &R, &R.adc.One RAD_ACA);
1603 return *this;
1604 }
1605
1606 template<typename Double>
1608operator+(const ADvari<Double> &L, Double R) {
1609 return *(new ADvar1<Double>(Hv_copy, L.Val + R, &L.adc.One, &L));
1610 }
1611
1612 template<typename Double>
1613 ADvar<Double>&
1615 ADVari *tcv = this->cv;
1616 this->cv = new ADVar1(Hv_copy, tcv->Val + R, &tcv->adc.One, tcv RAD_ACA);
1617 return *this;
1618 }
1619
1620 template<typename Double>
1622operator+(Double L, const ADvari<Double> &R) {
1623 return *(new ADvar1<Double>(Hv_copy, L + R.Val, &R.adc.One, &R));
1624 }
1625
1626 template<typename Double>
1627 ADvari<Double>&
1629 return *(new ADvar2<Double>(Hv_minusLR, L.Val - R.Val, &L, &L.adc.One, &R, &L.adc.negOne));
1630 }
1631
1632 template<typename Double>
1633 ADvar<Double>&
1635 ADVari *Lcv = this->cv;
1636 this->cv = new ADvar2<Double>(Hv_minusLR,Lcv->Val - R.Val, Lcv,
1637 &R.adc.One, &R, &R.adc.negOne RAD_ACA);
1638 return *this;
1639 }
1640
1641 template<typename Double>
1643operator-(const ADvari<Double> &L, Double R) {
1644 return *(new ADvar1<Double>(Hv_copy, L.Val - R, &L.adc.One, &L));
1645 }
1646
1647 template<typename Double>
1648 ADvar<Double>&
1650 ADVari *tcv = this->cv;
1651 this->cv = new ADVar1(Hv_copy, tcv->Val - R, &tcv->adc.One, tcv RAD_ACA);
1652 return *this;
1653 }
1654
1655 template<typename Double>
1657operator-(Double L, const ADvari<Double> &R) {
1658 return *(new ADvar1<Double>(Hv_negate, L - R.Val, &R.adc.negOne, &R));
1659 }
1660
1661 template<typename Double>
1662 ADvari<Double>&
1664 return *(new ADvar2<Double>(Hv_timesLR, L.Val * R.Val, &L, &R.Val, &R, &L.Val));
1665 }
1666
1667 template<typename Double>
1668 ADvar<Double>&
1670 ADVari *Lcv = this->cv;
1671 this->cv = new ADvar2<Double>(Hv_timesLR, Lcv->Val * R.Val, Lcv,
1672 &R.Val, &R, &Lcv->Val RAD_ACA);
1673 return *this;
1674 }
1675
1676 template<typename Double>
1678operator*(const ADvari<Double> &L, Double R) {
1679 return *(new ADvar1s<Double>(L.Val * R, R, &L));
1680 }
1681
1682 template<typename Double>
1683 ADvar<Double>&
1685 ADVari *Lcv = this->cv;
1686 this->cv = new ADvar1s<Double>(Lcv->Val * R, R, Lcv RAD_ACA);
1687 return *this;
1688 }
1689
1690 template<typename Double>
1692operator*(Double L, const ADvari<Double> &R) {
1693 return *(new ADvar1s<Double>(L * R.Val, L, &R));
1694 }
1695
1696 template<typename Double>
1697 ADvari<Double>&
1699 Double Lv = L.Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv, qpL = q*pL;
1700 return *(new ADvar2q<Double>(q, pL, -qpL, -pL*pL, 2.*pL*qpL, &L, &R));
1701 }
1702
1703 template<typename Double>
1704 ADvar<Double>&
1706 ADVari *Lcv = this->cv;
1707 Double Lv = Lcv->Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv, qpL = q*pL;
1708 this->cv = new ADvar2q<Double>(q, pL, -qpL, -pL*pL, 2.*pL*qpL, Lcv, &R RAD_ACA);
1709 return *this;
1710 }
1711
1712 template<typename Double>
1714operator/(const ADvari<Double> &L, Double R) {
1715 return *(new ADvar1s<Double>(L.Val / R, 1./R, &L));
1716 }
1717
1718 template<typename Double>
1719 ADvari<Double>&
1720operator/(Double L, const ADvari<Double> &R) {
1721 Double recip = 1. / R.Val;
1722 Double q = L * recip;
1723 Double d1 = -q*recip;
1724 return *(new ADvar1g<Double>(q, d1, -q*d1, &R));
1725 }
1726
1727 template<typename Double>
1728 ADvar<Double>&
1730 ADVari *Lcv = this->cv;
1731 this->cv = new ADvar1s<Double>(Lcv->Val / R, 1./R, Lcv RAD_ACA);
1732 return *this;
1733 }
1734
1735 template<typename Double>
1738 Double t = v.Val, t1 = 1. - t*t, d1 = -1./std::sqrt(t1);
1739 return *(new ADvar1g<Double>(std::acos(t), d1, t*d1/t1, &v));
1740 }
1741
1742 template<typename Double>
1743 ADvari<Double>&
1745 Double d1, t, t1, t2;
1746 t = v.Val;
1747 t1 = std::sqrt(t2 = t*t - 1.);
1748 d1 = 1. / t1;
1749 return *(new ADvar1g<Double>(std::log(t + t1), d1, -t*d1/t2, &v));
1750 }
1751
1752 template<typename Double>
1753 ADvari<Double>&
1755 Double d1, t, t1;
1756 t = v.Val;
1757 d1 = 1. / std::sqrt(t1 = 1. - t*t);
1758 return *(new ADvar1g<Double>(std::asin(t), d1, t*d1/t1, &v));
1759 }
1760
1761 template<typename Double>
1762 ADvari<Double>&
1764 Double d1, t, t1, t2, td;
1765 t = v.Val;
1766 t1 = std::sqrt(t2 = t*t + 1.);
1767 d1 = 1. / t1;
1768 td = 1.;
1769 if (t < 0.)
1770 td = -1.;
1771 return *(new ADvar1g<Double>(td*std::log(t*td + t1), d1, -(t/t2)*d1, &v));
1772 }
1773
1774 template<typename Double>
1775 ADvari<Double>&
1777 Double t = v.Val, d1 = 1./(1. + t*t);
1778 return *(new ADvar1g<Double>(std::atan(t), d1, -(t+t)*d1*d1, &v));
1779 }
1780
1781 template<typename Double>
1782 ADvari<Double>&
1784 Double t = v.Val, d1 = 1./(1. - t*t);
1785 return *(new ADvar1g<Double>(0.5*std::log((1.+t)/(1.-t)), d1, (t+t)*d1*d1, &v));
1786 }
1787
1788 template<typename Double>
1789 ADvari<Double>&
1791 Double R2, t, t2, x, x2, y, y2;
1792 x = L.Val;
1793 y = R.Val;
1794 x2 = x*x;
1795 y2 = y*y;
1796 t = 1./(x2 + y2);
1797 t2 = t*t;
1798 R2 = 2.*t2*x*y;
1799 return *(new ADvar2g<Double>(std::atan2(x,y), y*t, -x*t, -R2, t2*(x2 - y2), R2, &L, &R));
1800 }
1801
1802 template<typename Double>
1803 ADvari<Double>&
1804atan2(Double x, const ADvari<Double> &R) {
1805 Double t, x2, y, y2;
1806 y = R.Val;
1807 x2 = x*x;
1808 y2 = y*y;
1809 t = 1./(x2 + y2);
1810 return *(new ADvar1g<Double>(std::atan2(x,y), -x*t, 2.*t*t*x*y, &R));
1811 }
1812
1813 template<typename Double>
1814 ADvari<Double>&
1815atan2(const ADvari<Double> &L, Double y) {
1816 Double t, x, x2, y2;
1817 x = L.Val;
1818 x2 = x*x;
1819 y2 = y*y;
1820 t = 1./(x2 + y2);
1821 return *(new ADvar1g<Double>(std::atan2(x,y), y*t, -2.*t*t*x*y, &L));
1822 }
1823
1824 template<typename Double>
1825 ADvari<Double>&
1827 const ADvari<Double> &x = L.Val >= R.Val ? L : R;
1828 return *(new ADvar1<Double>(Hv_copy, x.Val, &x.adc.One, &x));
1829 }
1830
1831 template<typename Double>
1832 ADvari<Double>&
1833max(Double L, const ADvari<Double> &R) {
1834 if (L >= R.Val)
1835 return *(new ADvari<Double>(Hv_const, L));
1836 return *(new ADvar1<Double>(Hv_copy, R.Val, &R.adc.One, &R));
1837 }
1838
1839 template<typename Double>
1840 ADvari<Double>&
1841max(const ADvari<Double> &L, Double R) {
1842 if (L.Val >= R)
1843 return *(new ADvar1<Double>(Hv_copy, L.Val, &L.adc.One, &L));
1844 return *(new ADvari<Double>(Hv_const, R));
1845 }
1846
1847 template<typename Double>
1848 ADvari<Double>&
1850 const ADvari<Double> &x = L.Val <= R.Val ? L : R;
1851 return *(new ADvar1<Double>(Hv_copy, x.Val, &x.adc.One, &x));
1852 }
1853
1854 template<typename Double>
1855 ADvari<Double>&
1856min(Double L, const ADvari<Double> &R) {
1857 if (L <= R.Val)
1858 return *(new ADvari<Double>(Hv_const, L));
1859 return *(new ADvar1<Double>(Hv_copy, R.Val, &R.adc.One, &R));
1860 }
1861
1862 template<typename Double>
1863 ADvari<Double>&
1864min(const ADvari<Double> &L, Double R) {
1865 if (L.Val <= R)
1866 return *(new ADvar1<Double>(Hv_copy, L.Val, &L.adc.One, &L));
1867 return *(new ADvari<Double>(Hv_const, R));
1868 }
1869
1870 template<typename Double>
1871 ADvari<Double>&
1873 Double t = std::cos(v.Val);
1874 return *(new ADvar1g<Double>(t, -std::sin(v.Val), -t, &v));
1875 }
1876
1877 template<typename Double>
1878 ADvari<Double>&
1880 Double t = std::cosh(v.Val);
1881 return *(new ADvar1g<Double>(t, std::sinh(v.Val), t, &v));
1882 }
1883
1884 template<typename Double>
1885 ADvari<Double>&
1887 Double t = std::exp(v.Val);
1888 return *(new ADvar1g<Double>(t, t, t, &v));
1889 }
1890
1891 template<typename Double>
1892 ADvari<Double>&
1894 Double x = v.Val, d1 = 1. / x;
1895 return *(new ADvar1g<Double>(std::log(x), d1, -d1*d1, &v));
1896 }
1897
1898 template<typename Double>
1899 ADvari<Double>&
1901 static double num = 1. / std::log(10.);
1902 Double d1, t, x;
1903 x = v.Val;
1904 t = 1. / x;
1905 d1 = num * t;
1906 return *(new ADvar1g<Double>(std::log10(x), d1, -d1*t, &v));
1907 }
1908
1909 template<typename Double>
1910 ADvari<Double>&
1912 Double dx, dy, t, x, xlog, xym1, y;
1913 x = L.Val;
1914 y = R.Val;
1915 t = std::pow(x,y);
1916 xym1 = t / x;
1917 xlog = std::log(x);
1918 dx = y*xym1;
1919 dy = t * xlog;
1920 return *(new ADvar2g<Double>(t, dx, dy, (y-1.)*dx/x, xym1*(1. + y*xlog), dy*xlog, &L, &R));
1921 }
1922
1923 template<typename Double>
1924 ADvari<Double>&
1925pow(Double x, const ADvari<Double> &R) {
1926 Double dy, t, xlog, y;
1927 y = R.Val;
1928 t = std::pow(x,y);
1929 xlog = std::log(x);
1930 dy = t * xlog;
1931 return *(new ADvar1g<Double>(t, dy, dy*xlog, &R));
1932 }
1933
1934 template<typename Double>
1935 ADvari<Double>&
1936pow(const ADvari<Double> &L, Double y) {
1937 Double dx, t, x;
1938 x = L.Val;
1939 t = std::pow(x,y);
1940 dx = y*t/x;
1941 return *(new ADvar1g<Double>(t, dx, (y-1.)*dx/x, &L));
1942 }
1943
1944 template<typename Double>
1945 ADvari<Double>&
1947 Double t = std::sin(v.Val);
1948 return *(new ADvar1g<Double>(t, std::cos(v.Val), -t, &v));
1949 }
1950
1951 template<typename Double>
1952 ADvari<Double>&
1954 Double t = std::sinh(v.Val);
1955 return *(new ADvar1g<Double>(t, std::cosh(v.Val), t, &v));
1956 }
1957
1958 template<typename Double>
1959 ADvari<Double>&
1961 Double t = std::sqrt(v.Val);
1962 Double d1 = 0.5 / t;
1963 return *(new ADvar1g<Double>(t, d1, -0.5*d1/v.Val, &v));
1964 }
1965
1966 template<typename Double>
1967 ADvari<Double>&
1969 Double d1, rv, t;
1970 rv = std::tan(v.Val);
1971 t = 1. / std::cos(v.Val);
1972 d1 = t*t;
1973 return *(new ADvar1g<Double>(rv, d1, (rv+rv)*d1, &v));
1974 }
1975
1976 template<typename Double>
1977 ADvari<Double>&
1979 Double d1, rv, t;
1980 rv = std::tanh(v.Val);
1981 t = 1. / std::cosh(v.Val);
1982 d1 = t*t;
1983 return *(new ADvar1g<Double>(rv, d1, -(rv+rv)*d1, &v));
1984 }
1985
1986 template<typename Double>
1987 ADvari<Double>&
1989 Double t, p;
1990 p = 1.;
1991 if ((t = v.Val) < 0) {
1992 t = -t;
1993 p = -p;
1994 }
1995 return *(new ADvar1g<Double>(t, p, 0., &v));
1996 }
1997
1998 template<typename Double>
1999 ADvari<Double>&
2000fabs(const ADvari<Double> &v) { // Synonym for "abs"
2001 // "fabs" is not the best choice of name,
2002 // but this name is used at Sandia.
2003 Double t, p;
2004 p = 1.;
2005 if ((t = v.Val) < 0) {
2006 t = -t;
2007 p = -p;
2008 }
2009 return *(new ADvar1g<Double>(t, p, 0., &v));
2010 }
2011
2012 template<typename Double>
2013 ADvari<Double>&
2014ADf1(Double f, Double g, Double h, const ADvari<Double> &x) {
2015 return *(new ADvar1g<Double>(f, g, h, &x));
2016 }
2017
2018 template<typename Double>
2019 inline ADvari<Double>&
2020ADf1(Double f, Double g, Double h, const IndepADvar<Double> &x) {
2021 return *(new ADvar1g<Double>(f, g, h, x.cv));
2022 }
2023
2024 template<typename Double>
2025 ADvari<Double>&
2026ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
2027 const ADvari<Double> &x, const ADvari<Double> &y) {
2028 return *(new ADvar2g<Double>(f, gx, gy, hxx, hxy, hyy, &x, &y));
2029 }
2030
2031 template<typename Double>
2033ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
2034 const ADvari<Double> &x, const IndepADvar<Double> &y) {
2035 return *(new ADvar2g<Double>(f, gx, gy, hxx, hxy, hyy, &x, y.cv));
2036 }
2037
2038 template<typename Double>
2039 ADvari<Double>&
2040ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
2041 const IndepADvar<Double> &x, const ADvari<Double> &y) {
2042 return *(new ADvar2g<Double>(f, gx, gy, hxx, hxy, hyy, x.cv, &y));
2043 }
2044
2045 template<typename Double>
2046 ADvari<Double>&
2047ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
2048 const IndepADvar<Double> &x, const IndepADvar<Double> &y) {
2049 return *(new ADvar2g<Double>(f, gx, gy, hxx, hxy, hyy, x.cv, y.cv));
2050 }
2051
2052 template<typename Double>
2053 ADvari<Double>&
2054ADfn(Double f, int n, const IndepADvar<Double> *x, const Double *g, const Double *h) {
2055 return *(new ADvarn<Double>(f, n, x, g, h));
2056 }
2057
2058 template<typename Double>
2059 inline ADvari<Double>&
2060ADfn(Double f, int n, const ADvar<Double> *x, const Double *g, const Double *h) {
2061 return ADfn<Double>(f, n, (IndepADvar<Double>*)x, g, h);
2062 }
2063
2064 template<typename Double>
2065 void
2066ADcontext<Double>::Hvprod(int n, ADvar<Double> **x, Double *v, Double *hv)
2067{
2068 ADVari *a, *aL, *aR, **ap, **ape;
2069 ADVari_block *b, *b0;
2070 DErp *d;
2071 Double aO, adO, *g, *h, *h0, t, tL, tR;
2072 int i, j, k, m;
2073 for(i = 0; i < n; i++) {
2074 a = x[i]->cv;
2075 a->dO = v[i];
2076 a->aO = a->adO = 0.;
2077 }
2078 ADVari::adc.Aibusy->limit = ADVari::adc.Ainext;
2079 a = 0;
2080 for(b0 = 0, b = ADVari::adc.AiFirst; b; b0 = b, b = b->next) {
2081 ap = b->pADvari;
2082 ape = b->limit;
2083 while(ap < ape) {
2084 a = *ap++;
2085 a->aO = a->adO = 0.;
2086 switch(a->opclass) {
2087 case Hv_copy:
2088 a->dO = ((ADVar1*)a)->d.c->dO;
2089 break;
2090 case Hv_binary:
2091 a->dO = ((ADVar2g*)a)->pL * ((ADVar2g*)a)->dL.c->dO
2092 + ((ADVar2g*)a)->pR * ((ADVar2g*)a)->dR.c->dO;
2093 break;
2094 case Hv_unary:
2095 a->dO = ((ADVar1g*)a)->pL * ((ADVar1g*)a)->d.c->dO;
2096 break;
2097 case Hv_negate:
2098 a->dO = -((ADVar1*)a)->d.c->dO;
2099 break;
2100 case Hv_plusLR:
2101 a->dO = ((ADVar2*)a)->dL.c->dO + ((ADVar2*)a)->dR.c->dO;
2102 break;
2103 case Hv_minusLR:
2104 a->dO = ((ADVar2*)a)->dL.c->dO - ((ADVar2*)a)->dR.c->dO;
2105 break;
2106 case Hv_timesL:
2107 a->dO = ((ADVar1s*)a)->pL * ((ADVar1s*)a)->d.c->dO;
2108 break;
2109 case Hv_timesLR:
2110 a->dO = ((ADVar2*)a)->dR.c->Val * ((ADVar2*)a)->dL.c->dO
2111 + ((ADVar2*)a)->dL.c->Val * ((ADVar2*)a)->dR.c->dO;
2112 break;
2113 case Hv_quotLR:
2114 a->dO = ((ADVar2q*)a)->pL * ((ADVar2q*)a)->dL.c->dO
2115 + ((ADVar2q*)a)->pR * ((ADVar2q*)a)->dR.c->dO;
2116 break;
2117 case Hv_nary:
2118 d = ((ADVarn*)a)->D;
2119 m = ((ADVarn*)a)->n;
2120 g = ((ADVarn*)a)->G;
2121 t = 0.;
2122 for(i = 0; i < m; i++)
2123 t += g[i] * d[i].c->dO;
2124 a->dO = t;
2125 break;
2126 case Hv_const:
2127 ;
2128 }
2129 }
2130 }
2131 if (a)
2132 a->adO = 1.;
2133 for(b = b0; b; b = b->prev) {
2134 ape = b->pADvari;
2135 ap = b->limit;
2136 while(ap > ape) {
2137 a = *--ap;
2138 aO = a->aO;
2139 adO = a->adO;
2140 switch(a->opclass) {
2141 case Hv_copy:
2142 aL = ((ADVar1*)a)->d.c;
2143 aL->aO += aO;
2144 aL->adO += adO;
2145 break;
2146 case Hv_binary:
2147 aL = ((ADVar2g*)a)->dL.c;
2148 aR = ((ADVar2g*)a)->dR.c;
2149 tL = adO*aL->dO;
2150 tR = adO*aR->dO;
2151 aL->aO += aO*((ADVar2g*)a)->pL
2152 + tL*((ADVar2g*)a)->pL2
2153 + tR*((ADVar2g*)a)->pLR;
2154 aR->aO += aO*((ADVar2g*)a)->pR
2155 + tL*((ADVar2g*)a)->pLR
2156 + tR*((ADVar2g*)a)->pR2;
2157 aL->adO += adO * ((ADVar2g*)a)->pL;
2158 aR->adO += adO * ((ADVar2g*)a)->pR;
2159 break;
2160 case Hv_unary:
2161 aL = ((ADVar1g*)a)->d.c;
2162 aL->aO += aO * ((ADVar1g*)a)->pL
2163 + adO * aL->dO * ((ADVar1g*)a)->pL2;
2164 aL->adO += adO * ((ADVar1g*)a)->pL;
2165 break;
2166 case Hv_negate:
2167 aL = ((ADVar1*)a)->d.c;
2168 aL->aO -= aO;
2169 aL->adO -= adO;
2170 break;
2171 case Hv_plusLR:
2172 aL = ((ADVar2*)a)->dL.c;
2173 aR = ((ADVar2*)a)->dR.c;
2174 aL->aO += aO;
2175 aL->adO += adO;
2176 aR->aO += aO;
2177 aR->adO += adO;
2178 break;
2179 case Hv_minusLR:
2180 aL = ((ADVar2*)a)->dL.c;
2181 aR = ((ADVar2*)a)->dR.c;
2182 aL->aO += aO;
2183 aL->adO += adO;
2184 aR->aO -= aO;
2185 aR->adO -= adO;
2186 break;
2187 case Hv_timesL:
2188 aL = ((ADVar1s*)a)->d.c;
2189 aL->aO += aO * (tL = ((ADVar1s*)a)->pL);
2190 aL->adO += adO * tL;
2191 break;
2192 case Hv_timesLR:
2193 aL = ((ADVar2*)a)->dL.c;
2194 aR = ((ADVar2*)a)->dR.c;
2195 aL->aO += aO * (tL = aR->Val) + adO*aR->dO;
2196 aR->aO += aO * (tR = aL->Val) + adO*aL->dO;
2197 aL->adO += adO * tL;
2198 aR->adO += adO * tR;
2199 break;
2200 case Hv_quotLR:
2201 aL = ((ADVar2q*)a)->dL.c;
2202 aR = ((ADVar2q*)a)->dR.c;
2203 tL = adO*aL->dO;
2204 tR = adO*aR->dO;
2205 aL->aO += aO*((ADVar2q*)a)->pL
2206 + tR*((ADVar2q*)a)->pLR;
2207 aR->aO += aO*((ADVar2q*)a)->pR
2208 + tL*((ADVar2q*)a)->pLR
2209 + tR*((ADVar2q*)a)->pR2;
2210 aL->adO += adO * ((ADVar2q*)a)->pL;
2211 aR->adO += adO * ((ADVar2q*)a)->pR;
2212 break;
2213 case Hv_nary:
2214 d = ((ADVarn*)a)->D;
2215 m = ((ADVarn*)a)->n;
2216 g = ((ADVarn*)a)->G;
2217 h0 = ((ADVarn*)a)->H;
2218 for(i = 0; i < m; i++) {
2219 aL = d[i].c;
2220 aL->adO += adO * (t = g[i]);
2221 aL->aO += t*aO;
2222 t = adO * aL->dO;
2223 for(h = h0, j = 0; j <= i; j++)
2224 d[j].c->aO += t * *h++;
2225 h0 = h--;
2226 for(k = j; j < m; j++)
2227 d[j].c->aO += t * *(h += k++);
2228 }
2229 case Hv_const:
2230 ;
2231 }
2232 }
2233 }
2234 for(i = 0; i < n; i++) {
2235 a = x[i]->cv;
2236 a->dO = 0.;
2237 hv[i] = a->aO;
2238 }
2239 }
2240
2241 template<typename Double>
2242 inline Double
2244 return x.Val;
2245 }
2246
2247#undef RAD_ACA
2248#define A (ADvari<Double>*)
2249#ifdef RAD_Const_WARN
2250#define C(x) (((x)->opno < 0) ? RAD_Const_Warn(x) : 0, *A x)
2251#else
2252#define C(x) *A x
2253#endif
2254#define T template<typename Double> inline
2255#define F ADvari<Double>&
2256#define Ai const ADvari<Double>&
2257#define AI const IndepADvar<Double>&
2258#define D Double
2259#define T2(r,f) \
2260 T r f(Ai L, AI R) { return f(L, C(R.cv)); }\
2261 T r f(AI L, Ai R) { return f(C(L.cv), R); }\
2262 T r f(AI L, AI R) { return f(C(L.cv), C(R.cv)); }\
2263 T r f(AI L, D R) { return f(C(L.cv), R); }\
2264 T r f(Ai L, Dtype R) { return f(L, (D)R); }\
2265 T r f(AI L, Dtype R) { return f(C(L.cv), (D)R); }\
2266 T r f(Ai L, long R) { return f(L, (D)R); }\
2267 T r f(AI L, long R) { return f(C(L.cv), (D)R); }\
2268 T r f(Ai L, int R) { return f(L, (D)R); }\
2269 T r f(AI L, int R) { return f(C(L.cv), (D)R); }\
2270 T r f(D L, AI R) { return f(L, C(R.cv)); }\
2271 T r f(Dtype L, Ai R) { return f((D)L, R); }\
2272 T r f(Dtype L, AI R) { return f((D)L, C(R.cv)); }\
2273 T r f(long L, Ai R) { return f((D)L, R); }\
2274 T r f(long L, AI R) { return f((D)L, C(R.cv)); }\
2275 T r f(int L, Ai R) { return f((D)L, R); }\
2276 T r f(int L, AI R) { return f((D)L, C(R.cv)); }
2277
2278T2(F, operator+)
2279T2(F, operator-)
2280T2(F, operator*)
2281T2(F, operator/)
2282T2(F, atan2)
2283T2(F, pow)
2284T2(F, max)
2285T2(F, min)
2286T2(int, operator<)
2287T2(int, operator<=)
2288T2(int, operator==)
2289T2(int, operator!=)
2290T2(int, operator>=)
2291T2(int, operator>)
2292
2293#undef T2
2294#undef D
2295
2296#define T1(f)\
2297 T F f(AI x) { return f(C(x.cv)); }
2298
2299T1(operator+)
2300T1(operator-)
2301T1(abs)
2302T1(acos)
2303T1(acosh)
2304T1(asin)
2305T1(asinh)
2306T1(atan)
2307T1(atanh)
2308T1(cos)
2309T1(cosh)
2310T1(exp)
2311T1(log)
2312T1(log10)
2313T1(sin)
2314T1(sinh)
2315T1(sqrt)
2316T1(tan)
2317T1(tanh)
2318T1(fabs)
2319
2321{ return *(new ADvar1<Double>(Hv_copy, x.cv->Val, &ADcontext<Double>::One, (ADvari<Double>*)x.cv)); }
2322
2324{ return *(new ADvar1<Double>(Hv_copy, x.Val, &ADcontext<Double>::One, (ADvari<Double>*)&x)); }
2325
2326#undef T1
2327#undef AI
2328#undef Ai
2329#undef F
2330#undef T
2331#undef A
2332#undef C
2333#undef Ttype
2334#undef Dtype
2335
2336} /* namespace Rad2 */
2337} /* namespace Sacado */
2338
2339#endif /* SACADO_TRAD2_H */
int RAD_Const_Warn(void *v)
expr expr dx(i)
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 T2(r, f)
#define D
#define Ttype
#define RAD_ACA
ADvari< Double > ADVari
static void aval_reset(void)
ADvari_block< Double > ADVari_block
ADvar1g< Double > ADVar1g
ADvar2g< Double > ADVar2g
ADvar2q< Double > ADVar2q
static void Weighted_Gradcomp(int, ADVar **, Double *)
static void Hvprod(int, ADVar **, Double *, Double *)
ADvar1< Double > ADVar1
ADvar1s< Double > ADVar1s
void * Memalloc(size_t len)
ADmemblock< Double > ADMemblock
ADvarn< Double > ADVarn
ADvar2< Double > ADVar2
static const Double negOne
void * new_ADmemblock(size_t)
void ADvari_record(ADVari *x)
ADvar1(Advari_Opclass oc, Double val1, const Double *a1, const ADVari *c1)
ADvari< Double > ADVari
ADvar1< Double > ADVar1
ADvar1g(Double val1, Double d1, Double d2, const ADVari *c1)
ADvar1s(Double val1, Double a1, const ADVari *c1)
ADvar1< Double > ADVar1
ADvar2(Advari_Opclass oc, Double val1, const ADVari *Lcv, const Double *Lc, const ADVari *Rcv, const Double *Rc)
ADvari< Double > ADVari
ADvar2< Double > ADVar2
ADvar2g(Double val1, Double Lp, Double Rp, Double L2, Double LR, Double R2, const ADVari *Lcv, const ADVari *Rcv)
ADvar2< Double > ADVar2
ADvar2q(Double val1, Double Lp, Double Rp, Double LR, Double R2, const ADVari *Lcv, const ADVari *Rcv)
ConstADvari< Double > ConstADVari
static void Weighted_Gradcomp(int n, ADvar **v, Double *w)
ADvar & operator/=(const ADVari &)
ADvar & operator=(const ADVari &x)
IndepADvar< Double > IndepADVar
static void Gradcomp()
ADvar & operator*=(const ADVari &)
static void Gradcomp(int wantgrad)
ADvar(const IndepADVar &x)
static void set_fpval_implies_const(bool newval)
ADvar & operator-=(const ADVari &)
ADvar & operator+=(const ADVari &)
ADvari< Double > ADVari
ADvar(const ADVari &x)
ADvar1< Double > ADVar1
static bool get_fpval_implies_const(void)
static bool setget_fpval_implies_const(bool newval)
ADvar & operator=(IndepADVar &x)
ADvar(const ADvar &x)
static void aval_reset()
ScalarType< value_type >::type scalar_type
Advari_Opclass opclass
IndepADvar< Double > IndepADVar
ADvari(Advari_Opclass oc, Double t)
static ADcontext< Double > adc
ADvari(Advari_Opclass oc, Double t, Double ta)
ADVari::IndepADVar IndepADVar
ADvari< Double > ADVari
ADvar< Double > ADVar
ADvarn(Double val1, int n1, const IndepADVar *x, const Double *g, const Double *h)
ConstADvar & operator=(ADVari &d)
ADVar::IndepADVar IndepADVar
ConstADvar & operator*=(Double)
ConstADvar & operator+=(ADVari &)
ConstADvar & operator-=(Double)
ConstADvar & operator*=(ADVari &)
ConstADvar & operator/=(ADVari &)
ADVar::ConstADVari ConstADVari
ConstADvar & operator-=(ADVari &)
ConstADvar & operator+=(Double)
ConstADvar & operator/=(Double)
ConstADvar & operator=(Double d)
static ConstADvari * lastcad
static void aval_reset(void)
static CADcontext< Double > cadc
static Derp * LastDerp
ADvari< Double > ADVari
ADvari< Double > * cv
IndepADvar & operator=(IndepADvar &x)
static void AD_Const(const IndepADvar &)
static void Hvprod(int n, ADVar **vp, Double *v, Double *hv)
static void Weighted_Gradcomp(int n, ADVar **v, Double *w)
static void Gradcomp(int wantgrad)
const double y
const char * p
ADvari< Double > & operator+(const ADvari< Double > &T)
ADvari< Double > & atan(const ADvari< Double > &v)
ADvari< Double > & abs(const ADvari< Double > &v)
void AD_Const1(Double *, const IndepADvar< Double > &)
ADvari< Double > & exp(const ADvari< Double > &v)
void AD_Const(const IndepADvar< Double > &)
ADvari< Double > & atanh(const ADvari< Double > &v)
ADvari< Double > & atan2(const ADvari< Double > &L, const ADvari< Double > &R)
ADvari< Double > & asinh(const ADvari< Double > &v)
ADvari< Double > & max(const ADvari< Double > &L, const ADvari< Double > &R)
ADvari< Double > & acosh(const ADvari< Double > &v)
ADvari< Double > & pow(const ADvari< Double > &L, const ADvari< Double > &R)
int operator<=(const ADvari< Double > &L, const ADvari< Double > &R)
ADvari< Double > & operator*(const ADvari< Double > &L, const ADvari< Double > &R)
ADvari< Double > & tan(const ADvari< Double > &v)
ADvari< Double > & tanh(const ADvari< Double > &v)
ADvari< Double > & log10(const ADvari< Double > &v)
ADvari< Double > & ADf2(Double f, Double gx, Double gy, const IndepADvar< Double > &x, const IndepADvar< Double > &y)
ADvari< Double > & min(const ADvari< Double > &L, const ADvari< Double > &R)
int operator>(const ADvari< Double > &L, const ADvari< Double > &R)
ADvari< Double > & log(const ADvari< Double > &v)
ADvari< Double > & sqrt(const ADvari< Double > &v)
int operator==(const ADvari< Double > &L, const ADvari< Double > &R)
ADvari< Double > & fabs(const ADvari< Double > &v)
int operator>=(const ADvari< Double > &L, const ADvari< Double > &R)
ADvari< Double > & operator/(const ADvari< Double > &L, const ADvari< Double > &R)
ADvari< Double > & sinh(const ADvari< Double > &v)
ADvari< Double > & operator-(const ADvari< Double > &T)
ADvari< Double > & ADfn(Double f, int n, const IndepADvar< Double > *x, const Double *g)
ADvari< Double > & ADf1(Double f, Double g, const IndepADvar< Double > &x)
int operator<(const ADvari< Double > &L, const ADvari< Double > &R)
ADvari< Double > & acos(const ADvari< Double > &v)
ADvari< Double > & sin(const ADvari< Double > &v)
ADvari< Double > & cosh(const ADvari< Double > &v)
IndepADvar< Double > & ADvar_operatoreq(IndepADvar< Double > *, const ADvari< Double > &)
ADvari< Double > & cos(const ADvari< Double > &v)
ADvari< Double > & asin(const ADvari< Double > &v)
int operator!=(const ADvari< Double > &L, const ADvari< Double > &R)
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