Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Sacado_radops.cpp
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// Support routines for the RAD package (Reverse Automatic Differentiation) --
31// a package specialized for function and gradient evaluations.
32// Written in 2004 by David M. Gay at Sandia National Labs, Albuquerque, NM.
33
34#include "Sacado_rad.hpp"
35
36namespace Sacado {
37namespace Radnt {
38
39#ifdef RAD_AUTO_AD_Const
40ADvari *ADvari::First_ADvari, **ADvari::Last_ADvari = &ADvari::First_ADvari;
41#undef RAD_DEBUG_BLOCKKEEP
42#else
43#ifdef RAD_DEBUG_BLOCKKEEP
44#if !(RAD_DEBUG_BLOCKKEEP > 0)
45#undef RAD_DEBUG_BLOCKKEEP
46#else
47extern "C" void _uninit_f2c(void *x, int type, long len);
48#define TYDREAL 5
49static ADmemblock *rad_Oldcurmb;
50static int rad_busy_blocks;
51#endif /*RAD_DEBUG_BLOCKKEEP > 0*/
52#endif /*RAD_DEBUG_BLOCKKEEP*/
53#endif /*RAD_AUTO_AD_Const*/
54
55Derp *Derp::LastDerp = 0;
56ADcontext ADvari::adc;
57CADcontext ConstADvari::cadc;
58ConstADvari *ConstADvari::lastcad;
59static int rad_need_reinit;
60#ifdef RAD_DEBUG_BLOCKKEEP
61static size_t rad_mleft_save;
62#endif
63
64const double CADcontext::One = 1., CADcontext::negOne = -1.;
65
67{
68 First.next = 0;
69 Busy = &First;
70 Free = 0;
71 Mbase = (char*)First.memblk;
72 Mleft = sizeof(First.memblk);
73 }
74
75 void*
77{
78 ADmemblock *mb, *mb0, *mb1, *mbf, *x;
79#ifdef RAD_AUTO_AD_Const
80 ADvari *a, *anext;
81 IndepADvar *v;
82#endif /* RAD_AUTO_AD_Const */
83
84 if (rad_need_reinit && this == &ADvari::adc) {
87#ifdef RAD_DEBUG_BLOCKKEEP
88 Mleft = rad_mleft_save;
89 if (Mleft < sizeof(First.memblk))
91 (sizeof(First.memblk) - Mleft)/sizeof(double));
92 if ((mb = Busy->next)) {
93 if (!(mb0 = rad_Oldcurmb))
94 mb0 = &First;
95 for(;; mb = mb->next) {
97 sizeof(First.memblk)/sizeof(double));
98 if (mb == mb0)
99 break;
100 }
101 }
102 rad_Oldcurmb = Busy;
103 if (rad_busy_blocks >= RAD_DEBUG_BLOCKKEEP) {
104 rad_busy_blocks = 0;
105 rad_Oldcurmb = 0;
106 mb0 = &First;
107 mbf = Free;
108 for(mb = Busy; mb != mb0; mb = mb1) {
109 mb1 = mb->next;
110 mb->next = mbf;
111 mbf = mb;
112 }
113 Free = mbf;
114 Busy = mb;
115 Mbase = (char*)First.memblk;
116 Mleft = sizeof(First.memblk);
117 }
118
119#else /* !RAD_DEBUG_BLOCKKEEP */
120
121 mb0 = &ADvari::adc.First;
122 mbf = ADvari::adc.Free;
123 for(mb = ADvari::adc.Busy; mb != mb0; mb = mb1) {
124 mb1 = mb->next;
125 mb->next = mbf;
126 mbf = mb;
127 }
128 ADvari::adc.Free = mbf;
129 ADvari::adc.Busy = mb;
132#ifdef RAD_AUTO_AD_Const
133 *ADvari::Last_ADvari = 0;
134 ADvari::Last_ADvari = &ADvari::First_ADvari;
135 if ((anext = ADvari::First_ADvari)) {
136 while((a = anext)) {
137 anext = a->Next;
138 if ((v = a->padv))
139 v->cv = new ADvari(v, a->Val);
140 }
141 }
142#endif /*RAD_AUTO_AD_Const*/
143#endif /* RAD_DEBUG_BLOCKKEEP */
144 if (Mleft >= len)
145 return Mbase + (Mleft -= len);
146 }
147
148 if ((x = Free))
149 Free = x->next;
150 else
151 x = new ADmemblock;
152#ifdef RAD_DEBUG_BLOCKKEEP
153 rad_busy_blocks++;
154#endif
155 x->next = Busy;
156 Busy = x;
157 return (Mbase = (char*)x->memblk) +
158 (Mleft = sizeof(First.memblk) - len);
159 }
160
161 void
163{
164 Derp *d;
165
166 if (rad_need_reinit && wantgrad) {
167 for(d = Derp::LastDerp; d; d = d->next)
168 d->c->aval = 0;
169 }
170 else {
171 rad_need_reinit = 1;
172#ifdef RAD_DEBUG_BLOCKKEEP
173 rad_mleft_save = ADvari::adc.Mleft;
174#endif
175 ADvari::adc.Mleft = 0;
176 }
177
178 if ((d = Derp::LastDerp) && wantgrad) {
179 d->b->aval = 1;
180 do d->c->aval += *d->a * d->b->aval;
181 while((d = d->next));
182 }
183 }
184
185 void
187{
188 Derp *d;
189 int i;
190
191 if (rad_need_reinit) {
192 for(d = Derp::LastDerp; d; d = d->next)
193 d->c->aval = 0;
194 }
195 else {
196 rad_need_reinit = 1;
197#ifdef RAD_DEBUG_BLOCKKEEP
198 rad_mleft_save = ADvari::adc.Mleft;
199#endif
200 ADvari::adc.Mleft = 0;
201 }
202 if ((d = Derp::LastDerp)) {
203 for(i = 0; i < n; i++)
204 v[i]->cv->aval = w[i];
205 do d->c->aval += *d->a * d->b->aval;
206 while((d = d->next));
207 }
208 }
209
210#ifdef RAD_AUTO_AD_Const
211
213{
214 cv = new ADvari(this,d);
215 }
216
218{
219 cv = new ADvari(this,d);
220 }
221
223{
224 cv = new ADvari(this,d);
225 }
226
227#else
230{
231 ADvari *x = new ADvari(d);
232 cv = x;
233 }
234
236{
237 ADvari *x = new ADvari((double)d);
238 cv = x;
239 }
240
242{
243 ADvari *x = new ADvari((double)d);
244 cv = x;
245 }
246
247#endif /*RAD_AUTO_AD_Const*/
248
249 void
251{
252#ifdef RAD_AUTO_AD_Const
253 ADvari *x = new ADvari(this,d);
254#else
256 ? new ConstADvari(d)
257 : new ADvari(d);
258#endif
259 cv = x;
260 }
261
263{
264 ConstADvari *x = new ConstADvari(0.);
265 cv = x;
266 }
267
268 void
270{
271 ConstADvari *x = new ConstADvari(d);
272 cv = x;
273 }
275{
276 ConstADvari *y = new ConstADvari(x.Val);
277 new Derp(&CADcontext::One, y, &x); /*for side effect; value not used */
278 cv = y;
279 }
280
281 void
283{
284 ConstADvari *ncv = new ConstADvari(v.val());
285#ifdef RAD_AUTO_AD_Const
286 v.cv->padv = 0;
287#endif
288 ((IndepADvar*)&v)->cv = ncv;
289 }
290
291 void
293{
295 while(x) {
296 x->aval = 0;
297 x = x->prevcad;
298 }
299 }
300
301#ifdef RAD_AUTO_AD_Const
302
303ADvari::ADvari(const IndepADvar *x, double d): Val(d), aval(0.)
304{
305 *Last_ADvari = this;
306 Last_ADvari = &Next;
307 padv = (IndepADvar*)x;
308 }
309
310ADvar1::ADvar1(const IndepADvar *x, const IndepADvar &y):
311 ADvari(y.cv->Val), d(&CADcontext::One, this, y.cv)
312{
313 *ADvari::Last_ADvari = this;
314 ADvari::Last_ADvari = &Next;
315 padv = (IndepADvar*)x;
316 }
317
318ADvar1::ADvar1(const IndepADvar *x, const ADvari &y):
319 ADvari(y.Val), d(&CADcontext::One, this, &y)
320{
321 *ADvari::Last_ADvari = this;
322 ADvari::Last_ADvari = &Next;
323 padv = (IndepADvar*)x;
324 }
325
326#else
327
328 ADvar&
330{ This->cv = new ADvar1(x.Val, &CADcontext::One, (ADvari*)&x); return *This; }
331
334{ This->cv = new ADvar1(x.Val, &CADcontext::One, (ADvari*)&x); return *This; }
335
336#endif /*RAD_AUTO_AD_Const*/
337
340{
341#ifdef RAD_AUTO_AD_Const
342 if (cv)
343 cv->padv = 0;
344 cv = new ADvari(this,d);
345#else
346 cv = new ADvari(d);
347#endif
348 return *this;
349 }
350
351 ADvar&
353{
354#ifdef RAD_AUTO_AD_Const
355 if (cv)
356 cv->padv = 0;
357 cv = new ADvari(this,d);
358#else
360 ? new ConstADvari(d)
361 : new ADvari(d);
362#endif
363 return *this;
364 }
365
366 ADvari&
367operator-(const ADvari &T) {
368 return *(new ADvar1(-T.Val, &CADcontext::negOne, &T));
369 }
370
371 ADvari&
372operator+(const ADvari &L, const ADvari &R) {
373 return *(new ADvar2(L.Val + R.Val, &L, &CADcontext::One, &R, &CADcontext::One));
374 }
375
376 ADvar&
378 ADvari *Lcv = cv;
379#ifdef RAD_AUTO_AD_Const
380 Lcv->padv = 0;
381#endif
382 cv = new ADvar2(Lcv->Val + R.Val, Lcv, &CADcontext::One, &R, &CADcontext::One);
383 return *this;
384 }
385
386 ADvari&
387operator+(const ADvari &L, double R) {
388 return *(new ADvar1(L.Val + R, &CADcontext::One, &L));
389 }
390
391 ADvar&
393 ADvari *tcv = cv;
394#ifdef RAD_AUTO_AD_Const
395 tcv->padv = 0;
396#endif
397 cv = new ADvar1(tcv->Val + R, &CADcontext::One, tcv);
398 return *this;
399 }
400
401 ADvari&
402operator+(double L, const ADvari &R) {
403 return *(new ADvar1(L + R.Val, &CADcontext::One, &R));
404 }
405
406 ADvari&
407operator-(const ADvari &L, const ADvari &R) {
408 return *(new ADvar2(L.Val - R.Val, &L, &CADcontext::One, &R, &CADcontext::negOne));
409 }
410
411 ADvar&
413 ADvari *Lcv = cv;
414#ifdef RAD_AUTO_AD_Const
415 Lcv->padv = 0;
416#endif
417 cv = new ADvar2(Lcv->Val - R.Val, Lcv, &CADcontext::One, &R, &CADcontext::negOne);
418 return *this;
419 }
420
421 ADvari&
422operator-(const ADvari &L, double R) {
423 return *(new ADvar1(L.Val - R, &CADcontext::One, &L));
424 }
425
426 ADvar&
428 ADvari *tcv = cv;
429#ifdef RAD_AUTO_AD_Const
430 tcv->padv = 0;
431#endif
432 cv = new ADvar1(tcv->Val - R, &CADcontext::One, tcv);
433 return *this;
434 }
435
436 ADvari&
437operator-(double L, const ADvari &R) {
438 return *(new ADvar1(L - R.Val, &CADcontext::negOne, &R));
439 }
440
441 ADvari&
442operator*(const ADvari &L, const ADvari &R) {
443 return *(new ADvar2(L.Val * R.Val, &L, &R.Val, &R, &L.Val));
444 }
445
446 ADvar&
448 ADvari *Lcv = cv;
449#ifdef RAD_AUTO_AD_Const
450 Lcv->padv = 0;
451#endif
452 cv = new ADvar2(Lcv->Val * R.Val, Lcv, &R.Val, &R, &Lcv->Val);
453 return *this;
454 }
455
456 ADvari&
457operator*(const ADvari &L, double R) {
458 return *(new ADvar1s(L.Val * R, R, &L));
459 }
460
461 ADvar&
463 ADvari *Lcv = cv;
464#ifdef RAD_AUTO_AD_Const
465 Lcv->padv = 0;
466#endif
467 cv = new ADvar1s(Lcv->Val * R, R, Lcv);
468 return *this;
469 }
470
471 ADvari&
472operator*(double L, const ADvari &R) {
473 return *(new ADvar1s(L * R.Val, L, &R));
474 }
475
476 ADvari&
477operator/(const ADvari &L, const ADvari &R) {
478 double Lv = L.Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv;
479 return *(new ADvar2q(q, pL, -q*pL, &L, &R));
480 }
481
482 ADvar&
484 ADvari *Lcv = cv;
485#ifdef RAD_AUTO_AD_Const
486 Lcv->padv = 0;
487#endif
488 double Lv = Lcv->Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv;
489 cv = new ADvar2q(q, pL, -q*pL, Lcv, &R);
490 return *this;
491 }
492
493 ADvari&
494operator/(const ADvari &L, double R) {
495 return *(new ADvar1s(L.Val / R, 1./R, &L));
496 }
497
498 ADvari&
499operator/(double L, const ADvari &R) {
500 double recip = 1. / R.Val;
501 double q = L * recip;
502 return *(new ADvar1s(q, -q*recip, &R));
503 }
504
505 ADvar&
507 ADvari *Lcv = cv;
508#ifdef RAD_AUTO_AD_Const
509 Lcv->padv = 0;
510#endif
511 cv = new ADvar1s(Lcv->Val / R, 1./R, Lcv);
512 return *this;
513 }
514
515 ADvari&
516acos(const ADvari &v) {
517 double t = v.Val;
518 return *(new ADvar1s(acos(t), -1./sqrt(1. - t*t), &v));
519 }
520
521 ADvari&
522acosh(const ADvari &v) {
523 double t = v.Val, t1 = sqrt(t*t - 1.);
524 return *(new ADvar1s(log(t + t1), 1./t1, &v));
525 }
526
527 ADvari&
528asin(const ADvari &v) {
529 double t = v.Val;
530 return *(new ADvar1s(asin(t), 1./sqrt(1. - t*t), &v));
531 }
532
533 ADvari&
534asinh(const ADvari &v) {
535 double t = v.Val, td = 1., t1 = sqrt(t*t + 1.);
536 if (t < 0.) {
537 t = -t;
538 td = -1.;
539 }
540 return *(new ADvar1s(td*log(t + t1), 1./t1, &v));
541 }
542
543 ADvari&
544atan(const ADvari &v) {
545 double t = v.Val;
546 return *(new ADvar1s(atan(t), 1./(1. + t*t), &v));
547 }
548
549 ADvari&
550atanh(const ADvari &v) {
551 double t = v.Val;
552 return *(new ADvar1s(0.5*log((1.+t)/(1.-t)), 1./(1. - t*t), &v));
553 }
554
555 ADvari&
556max(const ADvari &L, const ADvari &R) {
557 const ADvari &x = L.Val >= R.Val ? L : R;
558 return *(new ADvar1(x.Val, &CADcontext::One, &x));
559 }
560
561 ADvari&
562max(double L, const ADvari &R) {
563 if (L >= R.Val)
564 return *(new ADvari(L));
565 return *(new ADvar1(R.Val, &CADcontext::One, &R));
566 }
567
568 ADvari&
569max(const ADvari &L, double R) {
570 if (L.Val >= R)
571 return *(new ADvar1(L.Val, &CADcontext::One, &L));
572 return *(new ADvari(R));
573 }
574
575 ADvari&
576min(const ADvari &L, const ADvari &R) {
577 const ADvari &x = L.Val <= R.Val ? L : R;
578 return *(new ADvar1(x.Val, &CADcontext::One, &x));
579 }
580
581 ADvari&
582min(double L, const ADvari &R) {
583 if (L <= R.Val)
584 return *(new ADvari(L));
585 return *(new ADvar1(R.Val, &CADcontext::One, &R));
586 }
587
588 ADvari&
589min(const ADvari &L, double R) {
590 if (L.Val <= R)
591 return *(new ADvar1(L.Val, &CADcontext::One, &L));
592 return *(new ADvari(R));
593 }
594
595 ADvari&
596atan2(const ADvari &L, const ADvari &R) {
597 double x = L.Val, y = R.Val, t = x*x + y*y;
598 return *(new ADvar2q(atan2(x,y), y/t, -x/t, &L, &R));
599 }
600
601 ADvari&
602atan2(double x, const ADvari &R) {
603 double y = R.Val, t = x*x + y*y;
604 return *(new ADvar1s(atan2(x,y), -x/t, &R));
605 }
606
607 ADvari&
608atan2(const ADvari &L, double y) {
609 double x = L.Val, t = x*x + y*y;
610 return *(new ADvar1s(atan2(x,y), y/t, &L));
611 }
612
613 ADvari&
614cos(const ADvari &v) {
615 return *(new ADvar1s(cos(v.Val), -sin(v.Val), &v));
616 }
617
618 ADvari&
619cosh(const ADvari &v) {
620 return *(new ADvar1s(cosh(v.Val), sinh(v.Val), &v));
621 }
622
623 ADvari&
624exp(const ADvari &v) {
625 ADvar1* rcv = new ADvar1(exp(v.Val), &v);
626 rcv->d.a = &rcv->Val;
627 rcv->d.b = rcv;
628 return *rcv;
629 }
630
631 ADvari&
632log(const ADvari &v) {
633 double x = v.Val;
634 return *(new ADvar1s(log(x), 1. / x, &v));
635 }
636
637 ADvari&
638log10(const ADvari &v) {
639 static double num = 1. / log(10.);
640 double x = v.Val;
641 return *(new ADvar1s(log10(x), num / x, &v));
642 }
643
644 ADvari&
645pow(const ADvari &L, const ADvari &R) {
646 double x = L.Val, y = R.Val, t = pow(x,y);
647 return *(new ADvar2q(t, y*t/x, t*log(x), &L, &R));
648 }
649
650 ADvari&
651pow(double x, const ADvari &R) {
652 double t = pow(x,R.Val);
653 return *(new ADvar1s(t, t*log(x), &R));
654 }
655
656 ADvari&
657pow(const ADvari &L, double y) {
658 double x = L.Val, t = pow(x,y);
659 return *(new ADvar1s(t, y*t/x, &L));
660 }
661
662 ADvari&
663sin(const ADvari &v) {
664 return *(new ADvar1s(sin(v.Val), cos(v.Val), &v));
665 }
666
667 ADvari&
668sinh(const ADvari &v) {
669 return *(new ADvar1s(sinh(v.Val), cosh(v.Val), &v));
670 }
671
672 ADvari&
673sqrt(const ADvari &v) {
674 double t = sqrt(v.Val);
675 return *(new ADvar1s(t, 0.5/t, &v));
676 }
677
678 ADvari&
679tan(const ADvari &v) {
680 double t = cos(v.Val);
681 return *(new ADvar1s(tan(v.Val), 1./(t*t), &v));
682 }
683
684 ADvari&
685tanh(const ADvari &v) {
686 double t = 1. / cosh(v.Val);
687 return *(new ADvar1s(tanh(v.Val), t*t, &v));
688 }
689
690 ADvari&
691fabs(const ADvari &v) { // "fabs" is not the best choice of name,
692 // but this name is used at Sandia.
693 double t, p;
694 p = 1;
695 if ((t = v.Val) < 0) {
696 t = -t;
697 p = -p;
698 }
699 return *(new ADvar1s(t, p, &v));
700 }
701
702 ADvari&
703ADf1(double f, double g, const ADvari &x) {
704 return *(new ADvar1s(f, g, &x));
705 }
706
707 ADvari&
708ADf2(double f, double gx, double gy, const ADvari &x, const ADvari &y) {
709 return *(new ADvar2q(f, gx, gy, &x, &y));
710 }
711
712ADvarn::ADvarn(double val1, int n1, const ADvar *x, const double *g): ADvari(val1), n(n1)
713{
714 Derp *d1, *dlast;
715 double *a1;
716 int i;
717
718 a1 = a = (double*)ADvari::adc.Memalloc(n*sizeof(*a));
719 d1 = Da = (Derp*)ADvari::adc.Memalloc(n*sizeof(Derp));
720 dlast = Derp::LastDerp;
721 for(i = 0; i < n1; i++, d1++) {
722 d1->next = dlast;
723 dlast = d1;
724 a1[i] = g[i];
725 d1->a = &a1[i];
726 d1->b = this;
727 d1->c = x[i].cv;
728 }
729 Derp::LastDerp = dlast;
730 }
731
732 ADvari&
733ADfn(double f, int n, const ADvar *x, const double *g) {
734 return *(new ADvarn(f, n, x, g));
735 }
736
737} // namespace Radnt
738} // namespace Sacado
#define R
void * new_ADmemblock(size_t)
static void Weighted_Gradcomp(int, ADvar **, double *)
ADvar & operator+=(const ADvari &)
ADvar & operator/=(const ADvari &)
ADvar & operator*=(const ADvari &)
ADvar & operator-=(const ADvari &)
void ADvar_ctr(double d)
ADvar & operator=(const ADvari &x)
static ADcontext adc
ADvarn(double val1, int n1, const ADvar *x, const double *g)
static const double One
static const double negOne
static ConstADvari * lastcad
const double * a
static Derp * LastDerp
const ADvari * c
const ADvari * b
IndepADvar & operator=(const IndepADvar &x)
friend void AD_Const(const IndepADvar &)
const double y
const char * p
ADvari & exp(const ADvari &v)
ADvari & atan(const ADvari &v)
ADvari & sinh(const ADvari &v)
ADvari & ADfn(double f, int n, const ADvar *x, const double *g)
ADvari & sqrt(const ADvari &v)
ADvari & tanh(const ADvari &v)
ADvari & ADf2(double f, double gx, double gy, const ADvari &x, const ADvari &y)
ADvari & asin(const ADvari &v)
ADvari & asinh(const ADvari &v)
ADvari & sin(const ADvari &v)
ADvari & min(const ADvari &L, const ADvari &R)
static int rad_need_reinit
ADvari & ADf1(double f, double g, const ADvari &x)
ADvari & tan(const ADvari &v)
ADvari & operator-(const ADvari &T)
ADvari & operator/(const ADvari &L, const ADvari &R)
ADvari & max(const ADvari &L, const ADvari &R)
ADvari & operator+(ADvari &T)
ADvari & acosh(const ADvari &v)
ADvari & acos(const ADvari &v)
ADvari & fabs(const ADvari &v)
ADvari & atanh(const ADvari &v)
ADvari & cos(const ADvari &v)
ADvar & ADvar_operatoreq(ADvar *This, const ADvari &x)
ADvari & log10(const ADvari &v)
ADvari & cosh(const ADvari &v)
ADvari & log(const ADvari &v)
ADvari & operator*(const ADvari &L, const ADvari &R)
ADvari & pow(const ADvari &L, const ADvari &R)
ADvari & atan2(const ADvari &L, const ADvari &R)
#define TYDREAL
Definition uninit.c:38
void _uninit_f2c(void *x, int type, long len)
Definition uninit.c:94