Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Sacado_radops2.cpp
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// Support routines for Hessian-vector products via the RAD package
31// (Reverse Automatic Differentiation) -- a package specialized for
32// function and gradient evaluations, and herewith extended for Hessian-
33// vector products. This specialization is for several Hessian-vector
34// products at one point, where the vectors are determined on the fly,
35// e.g., by the conjugate-gradient algorithm. Where the vectors are known
36// in advance, nestings of RAD and FAD might be more efficient.
37// RAD Hessian-vector product extension written in 2007 by David M. Gay at
38// Sandia National Labs, Albuquerque, NM.
39
40#include "Sacado_rad2.hpp"
41
42#ifndef SACADO_NO_NAMESPACE
43namespace Sacado {
44namespace Rad2d { // "2" for 2nd derivatives, "d" for "double"
45#endif
46
47#ifdef RAD_AUTO_AD_Const
48ADvari *ADvari::First_ADvari, **ADvari::Last_ADvari = &ADvari::First_ADvari;
49#undef RAD_DEBUG_BLOCKKEEP
50#else
51#ifdef RAD_DEBUG_BLOCKKEEP
52#if !(RAD_DEBUG_BLOCKKEEP > 0)
53#undef RAD_DEBUG_BLOCKKEEP
54#else
55extern "C" void _uninit_f2c(void *x, int type, long len);
56#define TYDREAL 5
57static ADmemblock *rad_Oldcurmb;
58static int rad_busy_blocks;
59#endif /*RAD_DEBUG_BLOCKKEEP > 0*/
60#endif /*RAD_DEBUG_BLOCKKEEP*/
61#endif /*RAD_AUTO_AD_Const*/
62
63Derp *Derp::LastDerp = 0;
64ADcontext ADvari::adc;
65CADcontext ConstADvari::cadc;
66ConstADvari *ConstADvari::lastcad;
67static int rad_need_reinit;
68#ifdef RAD_DEBUG_BLOCKKEEP
69static size_t rad_mleft_save;
70#endif
71
72const double CADcontext::One = 1., CADcontext::negOne = -1.;
73
75{
76 First.next = 0;
77 Busy = &First;
78 Free = 0;
79 Mbase = (char*)First.memblk;
80 Mleft = sizeof(First.memblk);
81 Aibusy = &AiFirst;
82 Aifree = 0;
86 }
87
88 void*
90{
91 ADmemblock *mb, *mb0, *mb1, *mbf, *x;
92 ADvari_block *b;
93#ifdef RAD_AUTO_AD_Const
94 ADvari *a, *anext;
95 IndepADvar *v;
96#endif /* RAD_AUTO_AD_Const */
97
98 if (rad_need_reinit && this == &ADvari::adc) {
100 Derp::LastDerp = 0;
101 Aibusy = b = &AiFirst;
102 Aifree = b->next;
103 b->next = b->prev = 0;
105#ifdef RAD_DEBUG_BLOCKKEEP
106 Mleft = rad_mleft_save;
107 if (Mleft < sizeof(First.memblk))
109 (sizeof(First.memblk) - Mleft)/sizeof(double));
110 if ((mb = Busy->next)) {
111 if (!(mb0 = rad_Oldcurmb))
112 mb0 = &First;
113 for(;; mb = mb->next) {
115 sizeof(First.memblk)/sizeof(double));
116 if (mb == mb0)
117 break;
118 }
119 }
120 rad_Oldcurmb = Busy;
121 if (rad_busy_blocks >= RAD_DEBUG_BLOCKKEEP) {
122 rad_busy_blocks = 0;
123 rad_Oldcurmb = 0;
124 mb0 = &First;
125 mbf = Free;
126 for(mb = Busy; mb != mb0; mb = mb1) {
127 mb1 = mb->next;
128 mb->next = mbf;
129 mbf = mb;
130 }
131 Free = mbf;
132 Busy = mb;
133 Mbase = (char*)First.memblk;
134 Mleft = sizeof(First.memblk);
135 }
136
137#else /* !RAD_DEBUG_BLOCKKEEP */
138
139 mb0 = &ADvari::adc.First;
140 mbf = ADvari::adc.Free;
141 for(mb = ADvari::adc.Busy; mb != mb0; mb = mb1) {
142 mb1 = mb->next;
143 mb->next = mbf;
144 mbf = mb;
145 }
146 ADvari::adc.Free = mbf;
147 ADvari::adc.Busy = mb;
150#ifdef RAD_AUTO_AD_Const
151 *ADvari::Last_ADvari = 0;
152 ADvari::Last_ADvari = &ADvari::First_ADvari;
153 if ((anext = ADvari::First_ADvari)) {
154 while((a = anext)) {
155 anext = a->Next;
156 if ((v = a->padv))
157 v->cv = new ADvari(v, a->Val);
158 }
159 }
160#endif /*RAD_AUTO_AD_Const*/
161#endif /* RAD_DEBUG_BLOCKKEEP */
162 if (Mleft >= len)
163 return Mbase + (Mleft -= len);
164 }
165
166 if ((x = Free))
167 Free = x->next;
168 else
169 x = new ADmemblock;
170#ifdef RAD_DEBUG_BLOCKKEEP
171 rad_busy_blocks++;
172#endif
173 x->next = Busy;
174 Busy = x;
175 return (Mbase = (char*)x->memblk) +
176 (Mleft = sizeof(First.memblk) - len);
177 }
178
179 void
181{
182 ADvari_block *ob, *nb;
183 ob = Aibusy;
184 ob->limit = Ailimit; // should be unnecessary, but harmless
185 if ((nb = Aifree))
186 Aifree = nb->next;
187 else
188 nb = new ADvari_block;
189 Aibusy = Aibusy->next = nb;
191 ob->next = nb;
192 nb->prev = ob;
193 nb->next = 0;
194 }
195
196 void
198{
199 Derp *d;
200
201 if (rad_need_reinit && wantgrad) {
202 for(d = Derp::LastDerp; d; d = d->next)
203 d->c->aval = 0;
204 }
205 else {
206 rad_need_reinit = 1;
207#ifdef RAD_DEBUG_BLOCKKEEP
208 rad_mleft_save = ADvari::adc.Mleft;
209#endif
210 ADvari::adc.Mleft = 0;
211 }
212
213 if ((d = Derp::LastDerp) && wantgrad) {
214 d->b->aval = 1;
215 do d->c->aval += *d->a * d->b->aval;
216 while((d = d->next));
217 }
218 }
219
220 void
222{
223 Derp *d;
224 int i;
225
226 if (rad_need_reinit) {
227 for(d = Derp::LastDerp; d; d = d->next)
228 d->c->aval = 0;
229 }
230 else {
231 rad_need_reinit = 1;
232#ifdef RAD_DEBUG_BLOCKKEEP
233 rad_mleft_save = ADvari::adc.Mleft;
234#endif
235 ADvari::adc.Mleft = 0;
236 }
237 if ((d = Derp::LastDerp)) {
238 for(i = 0; i < n; i++)
239 v[i]->cv->aval = w[i];
240 do d->c->aval += *d->a * d->b->aval;
241 while((d = d->next));
242 }
243 }
244
245#ifdef RAD_AUTO_AD_Const
246
248{
249 cv = new ADvari(this, d);
250 }
251
253{
254 cv = new ADvari(this, (double)d);
255 }
256
258{
259 cv = new ADvari(this, (double)d);
260 }
261
262#else
265{
266 ADvari *x = new ADvari(Hv_const, d);
267 cv = x;
268 }
269
271{
272 ADvari *x = new ADvari(Hv_const, (double)d);
273 cv = x;
274 }
275
277{
278 ADvari *x = new ADvari(Hv_const, (double)d);
279 cv = x;
280 }
281
282#endif /*RAD_AUTO_AD_Const*/
283
284 void
286{
287#ifdef RAD_AUTO_AD_Const
288 ADvari *x = new ADvari(this, d);
289#else
291 ? new ConstADvari(d)
292 : new ADvari(Hv_const, d);
293#endif
294 cv = x;
295 }
296
298{
299 ConstADvari *x = new ConstADvari(0.);
300 cv = x;
301 }
302
303 void
305{
306 ConstADvari *x = new ConstADvari(d);
307 cv = x;
308 }
310{
311 ConstADvari *y = new ConstADvari(x.Val);
312 new Derp(&CADcontext::One, y, &x); /* for side effect; value not used */
313 cv = y;
314 }
315
316 void
318{
319 ConstADvari *ncv = new ConstADvari(v.val());
320 ((IndepADvar*)&v)->cv = ncv;
321 }
322
323 void
325{
327 while(x) {
328 x->aval = 0;
329 x = x->prevcad;
330 }
331 }
332
333#ifdef RAD_AUTO_AD_Const
334
335ADvari::ADvari(const IndepADvar *x, double d): Val(d), aval(0.)
336{
338 *Last_ADvari = this;
339 Last_ADvari = &Next;
340 padv = (IndepADvar*)x;
341 }
342
343ADvar1::ADvar1(const IndepADvar *x, const IndepADvar &y):
344 ADvari(Hv_copy, y.cv->Val), d(&CADcontext::One, this, y.cv)
345{
346 *ADvari::Last_ADvari = this;
347 ADvari::Last_ADvari = &Next;
348 padv = (IndepADvar*)x;
349 }
350
351ADvar1::ADvar1(const IndepADvar *x, const ADvari &y):
352 ADvari(Hv_copy, y.Val), d(&CADcontext::One, this, &y)
353{
354 *ADvari::Last_ADvari = this;
355 ADvari::Last_ADvari = &Next;
356 padv = (IndepADvar*)x;
357 }
358
359#else
360
361 ADvar&
363{ This->cv = new ADvar1(Hv_copy, x.Val, &CADcontext::One, (ADvari*)&x); return *This; }
364
367{ This->cv = new ADvar1(Hv_copy, x.Val, &CADcontext::One, (ADvari*)&x); return *This; }
368
369#endif /*RAD_AUTO_AD_Const*/
370
371ADvar2q::ADvar2q(double val1, double Lp, double Rp,
372 double LR, double R2, const ADvari *Lcv, const ADvari *Rcv):
373 ADvar2(Hv_quotLR,val1,Lcv,&pL,Rcv,&pR),
374 pL(Lp), pR(Rp), pLR(LR), pR2(R2) { }
375
376ADvar2g::ADvar2g(double val1, double Lp, double Rp,
377 double L2, double LR, double R2, const ADvari *Lcv, const ADvari *Rcv):
378 ADvar2(Hv_binary,val1,Lcv,&pL,Rcv,&pR),
379 pL(Lp), pR(Rp), pL2(L2), pLR(LR), pR2(R2) { }
380
383{
384#ifdef RAD_AUTO_AD_Const
385 if (cv)
386 cv->padv = 0;
387 cv = new ADvari(this,d);
388#else
389 cv = new ADvari(Hv_const, d);
390#endif
391 return *this;
392 }
393
394 ADvar&
396{
397#ifdef RAD_AUTO_AD_Const
398 if (cv)
399 cv->padv = 0;
400 cv = new ADvari(this, d);
401#else
403 ? new ConstADvari(d)
404 : new ADvari(Hv_const, d);
405#endif
406 return *this;
407 }
408
409 ADvari&
411 return *(new ADvar1(Hv_negate, -T.Val, &CADcontext::negOne, &T));
412 }
413
414 ADvari&
415operator+(const ADvari &L, const ADvari &R) {
416 return *(new ADvar2(Hv_plusLR,L.Val + R.Val, &L, &CADcontext::One, &R, &CADcontext::One));
417 }
418
419 ADvar&
421 ADvari *Lcv = cv;
422#ifdef RAD_AUTO_AD_Const
423 Lcv->padv = 0;
424#endif
425 cv = new ADvar2(Hv_plusLR, Lcv->Val + R.Val, Lcv, &CADcontext::One, &R, &CADcontext::One);
426 return *this;
427 }
428
429 ADvari&
430operator+(const ADvari &L, double R) {
431 return *(new ADvar1(Hv_copy, L.Val + R, &CADcontext::One, &L));
432 }
433
434 ADvar&
436 ADvari *tcv = cv;
437#ifdef RAD_AUTO_AD_Const
438 tcv->padv = 0;
439#endif
440 cv = new ADvar1(Hv_copy, tcv->Val + R, &CADcontext::One, tcv);
441 return *this;
442 }
443
444 ADvari&
445operator+(double L, const ADvari &R) {
446 return *(new ADvar1(Hv_copy, L + R.Val, &CADcontext::One, &R));
447 }
448
449 ADvari&
450operator-(const ADvari &L, const ADvari &R) {
451 return *(new ADvar2(Hv_minusLR,L.Val - R.Val, &L, &CADcontext::One, &R, &CADcontext::negOne));
452 }
453
454 ADvar&
456 ADvari *Lcv = cv;
457#ifdef RAD_AUTO_AD_Const
458 Lcv->padv = 0;
459#endif
460 cv = new ADvar2(Hv_minusLR, Lcv->Val - R.Val, Lcv, &CADcontext::One, &R, &CADcontext::negOne);
461 return *this;
462 }
463
464 ADvari&
465operator-(const ADvari &L, double R) {
466 return *(new ADvar1(Hv_copy, L.Val - R, &CADcontext::One, &L));
467 }
468
469 ADvar&
471 ADvari *tcv = cv;
472#ifdef RAD_AUTO_AD_Const
473 tcv->padv = 0;
474#endif
475 cv = new ADvar1(Hv_copy, tcv->Val - R, &CADcontext::One, tcv);
476 return *this;
477 }
478
479 ADvari&
480operator-(double L, const ADvari &R) {
481 return *(new ADvar1(Hv_negate, L - R.Val, &CADcontext::negOne, &R));
482 }
483
484 ADvari&
485operator*(const ADvari &L, const ADvari &R) {
486 return *(new ADvar2(Hv_timesLR, L.Val * R.Val, &L, &R.Val, &R, &L.Val));
487 }
488
489 ADvar&
491 ADvari *Lcv = cv;
492#ifdef RAD_AUTO_AD_Const
493 Lcv->padv = 0;
494#endif
495 cv = new ADvar2(Hv_timesLR, Lcv->Val * R.Val, Lcv, &R.Val, &R, &Lcv->Val);
496 return *this;
497 }
498
499 ADvari&
500operator*(const ADvari &L, double R) {
501 return *(new ADvar1s(L.Val * R, R, &L));
502 }
503
504 ADvar&
506 ADvari *Lcv = cv;
507#ifdef RAD_AUTO_AD_Const
508 Lcv->padv = 0;
509#endif
510 cv = new ADvar1s(Lcv->Val * R, R, Lcv);
511 return *this;
512 }
513
514 ADvari&
515operator*(double L, const ADvari &R) {
516 return *(new ADvar1s(L * R.Val, L, &R));
517 }
518
519 ADvari&
520operator/(const ADvari &L, const ADvari &R) {
521 double Lv = L.Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv, qpL = q*pL;
522 return *(new ADvar2q(q, pL, -qpL, -pL*pL, 2.*pL*qpL, &L, &R));
523 }
524
525 ADvar&
527 ADvari *Lcv = cv;
528#ifdef RAD_AUTO_AD_Const
529 Lcv->padv = 0;
530#endif
531 double Lv = Lcv->Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv, qpL = q*pL;
532 cv = new ADvar2q(q, pL, -qpL, -pL*pL, 2.*pL*qpL, Lcv, &R);
533 return *this;
534 }
535
536 ADvari&
537operator/(const ADvari &L, double R) {
538 return *(new ADvar1s(L.Val / R, 1./R, &L));
539 }
540
541 ADvari&
542operator/(double L, const ADvari &R) {
543 double recip = 1. / R.Val;
544 double q = L * recip;
545 double d1 = -q*recip;
546 return *(new ADvar1g(q, d1, -q*d1, &R));
547 }
548
549 ADvar&
551 ADvari *Lcv = cv;
552#ifdef RAD_AUTO_AD_Const
553 Lcv->padv = 0;
554#endif
555 cv = new ADvar1s(Lcv->Val / R, 1./R, Lcv);
556 return *this;
557 }
558
559 ADvari&
560acos(const ADvari &v) {
561 double t = v.Val;
562 double t1 = 1. - t*t;
563 double d1 = -1. / sqrt(t1);
564 return *(new ADvar1g(acos(t), d1, t*d1/t1, &v));
565 }
566
567 ADvari&
568acosh(const ADvari &v) {
569 double d1, t, t1, t2;
570 t = v.Val;
571 t1 = sqrt(t2 = t*t - 1.);
572 d1 = 1./t1;
573 return *(new ADvar1g(log(t + t1), d1, -t*d1/t2, &v));
574 }
575
576 ADvari&
577asin(const ADvari &v) {
578 double d1, t, t1;
579 t = v.Val;
580 d1 = 1. / sqrt(t1 = 1. - t*t);
581 return *(new ADvar1g(asin(t), d1, t*d1/t1, &v));
582 }
583
584 ADvari&
585asinh(const ADvari &v) {
586 double d1, t, t1, t2, td;
587 t = v.Val;
588 t1 = sqrt(t2 = t*t + 1.);
589 d1 = 1. / t1;
590 td = 1.;
591 if (t < 0.)
592 td = -1.;
593 return *(new ADvar1g(td*log(t*td + t1), d1, -(t/t2)*d1, &v));
594 }
595
596 ADvari&
597atan(const ADvari &v) {
598 double d1, t;
599 t = v.Val;
600 d1 = 1./(1. + t*t);
601 return *(new ADvar1g(atan(t), d1, -(t+t)*d1*d1, &v));
602 }
603
604 ADvari&
605atanh(const ADvari &v) {
606 double t = v.Val;
607 double d1 = 1./(1. - t*t);
608 return *(new ADvar1g(0.5*log((1.+t)/(1.-t)), d1, (t+t)*d1*d1, &v));
609 }
610
611 ADvari&
612max(const ADvari &L, const ADvari &R) {
613 const ADvari &x = L.Val >= R.Val ? L : R;
614 return *(new ADvar1(Hv_copy, x.Val, &CADcontext::One, &x));
615 }
616
617 ADvari&
618max(double L, const ADvari &R) {
619 if (L >= R.Val)
620 return *(new ADvari(Hv_const, L));
621 return *(new ADvar1(Hv_copy, R.Val, &CADcontext::One, &R));
622 }
623
624 ADvari&
625max(const ADvari &L, double R) {
626 if (L.Val >= R)
627 return *(new ADvar1(Hv_copy, L.Val, &CADcontext::One, &L));
628 return *(new ADvari(Hv_const, R));
629 }
630
631 ADvari&
632min(const ADvari &L, const ADvari &R) {
633 const ADvari &x = L.Val <= R.Val ? L : R;
634 return *(new ADvar1(Hv_copy, x.Val, &CADcontext::One, &x));
635 }
636
637 ADvari&
638min(double L, const ADvari &R) {
639 if (L <= R.Val)
640 return *(new ADvari(Hv_const, L));
641 return *(new ADvar1(Hv_copy, R.Val, &CADcontext::One, &R));
642 }
643
644 ADvari&
645min(const ADvari &L, double R) {
646 if (L.Val <= R)
647 return *(new ADvar1(Hv_copy, L.Val, &CADcontext::One, &L));
648 return *(new ADvari(Hv_const, R));
649 }
650
651 ADvari&
652atan2(const ADvari &L, const ADvari &R) {
653 double x = L.Val, y = R.Val;
654 double x2 = x*x, y2 = y*y;
655 double t = 1./(x2 + y2);
656 double t2 = t*t;
657 double R2 = 2.*t2*x*y;
658 return *(new ADvar2g(atan2(x,y), y*t, -x*t, -R2, t2*(x2 - y2), R2, &L, &R));
659 }
660
661 ADvari&
662atan2(double x, const ADvari &R) {
663 double y = R.Val;
664 double x2 = x*x, y2 = y*y;
665 double t = 1./(x2 + y2);
666 return *(new ADvar1g(atan2(x,y), -x*t, 2.*t*t*x*y, &R));
667 }
668
669 ADvari&
670atan2(const ADvari &L, double y) {
671 double x = L.Val;
672 double x2 = x*x, y2 = y*y;
673 double t = 1./(x2 + y2);
674 return *(new ADvar1g(atan2(x,y), y*t, -2.*t*t*x*y, &L));
675 }
676
677 ADvari&
678cos(const ADvari &v) {
679 double t = cos(v.Val);
680 return *(new ADvar1g(t, -sin(v.Val), -t, &v));
681 }
682
683 ADvari&
684cosh(const ADvari &v) {
685 double t = cosh(v.Val);
686 return *(new ADvar1g(t, sinh(v.Val), t, &v));
687 }
688
689 ADvari&
690exp(const ADvari &v) {
691 double t = exp(v.Val);
692 return *(new ADvar1g(t, t, t, &v));
693 }
694
695 ADvari&
696log(const ADvari &v) {
697 double x = v.Val;
698 double d1 = 1. / x;
699 return *(new ADvar1g(log(x), d1, -d1*d1, &v));
700 }
701
702 ADvari&
703log10(const ADvari &v) {
704 static double num = 1. / log(10.);
705 double x = v.Val, t = 1. / x;
706 double d1 = num * t;
707 return *(new ADvar1g(log10(x), d1, -d1*t, &v));
708 }
709
710 ADvari&
711pow(const ADvari &L, const ADvari &R) {
712 double x = L.Val, y = R.Val, t = pow(x,y);
713 double xym1 = t/x;
714 double xlog = log(x);
715 double dx = y*xym1;
716 double dy = t * xlog;
717 return *(new ADvar2g(t, dx, dy, (y-1.)*dx/x, xym1*(1. + y*xlog), dy*xlog, &L, &R));
718 }
719
720 ADvari&
721pow(double x, const ADvari &R) {
722 double y = R.Val, t = pow(x,y);
723 double xlog = log(x);
724 double dy = t * xlog;
725 return *(new ADvar1g(t, dy, dy*xlog, &R));
726 }
727
728 ADvari&
729pow(const ADvari &L, double y) {
730 double x = L.Val, t = pow(x,y);
731 double dx = y*t/x;
732 return *(new ADvar1g(t, dx, (y-1.)*dx/x, &L));
733 }
734
735 ADvari&
736sin(const ADvari &v) {
737 double t = sin(v.Val);
738 return *(new ADvar1g(t, cos(v.Val), -t, &v));
739 }
740
741 ADvari&
742sinh(const ADvari &v) {
743 double t = sinh(v.Val);
744 return *(new ADvar1g(t, cosh(v.Val), t, &v));
745 }
746
747 ADvari&
748sqrt(const ADvari &v) {
749 double t = sqrt(v.Val);
750 double d1 = 0.5/t;
751 return *(new ADvar1g(t, d1, -0.5*d1/v.Val, &v));
752 }
753
754 ADvari&
755tan(const ADvari &v) {
756 double d1, rv, t;
757 rv = tan(v.Val);
758 t = 1. / cos(v.Val);
759 d1 = t*t;
760 return *(new ADvar1g(rv, d1, (rv+rv)*d1, &v));
761 }
762
763 ADvari&
764tanh(const ADvari &v) {
765 double d1, rv, t;
766 rv = tanh(v.Val);
767 t = 1. / cosh(v.Val);
768 d1 = t*t;
769 return *(new ADvar1g(rv, d1, -(rv+rv)*d1, &v));
770 }
771
772 ADvari&
773fabs(const ADvari &v) { // "fabs" is not the best choice of name,
774 // but this name is used at Sandia.
775 double t, p;
776 p = 1;
777 if ((t = v.Val) < 0) {
778 t = -t;
779 p = -p;
780 }
781 return *(new ADvar1g(t, p, 0., &v));
782 }
783
784 ADvari&
785ADf1(double f, double g, double h, const ADvari &x) {
786 return *(new ADvar1g(f, g, h, &x));
787 }
788
789 ADvari&
790ADf2(double f, double gx, double gy, double hxx, double hxy, double hyy,
791 const ADvari &x, const ADvari &y) {
792 return *(new ADvar2g(f, gx, gy, hxx, hxy, hyy, &x, &y));
793 }
794
795 ADvarn::ADvarn(double val1, int n1, const ADvar *x, const double *g, const double *h):
796 ADvari(Hv_nary,val1), n(n1)
797{
798 Derp *d1, *dlast;
799 double *a1;
800 int i, nh;
801
802 a1 = G = (double*)ADvari::adc.Memalloc(n1*sizeof(*g));
803 d1 = D = (Derp*)ADvari::adc.Memalloc(n1*sizeof(Derp));
804 dlast = Derp::LastDerp;
805 for(i = 0; i < n1; i++, d1++) {
806 d1->next = dlast;
807 dlast = d1;
808 a1[i] = g[i];
809 d1->a = &a1[i];
810 d1->b = this;
811 d1->c = x[i].cv;
812 }
813 Derp::LastDerp = dlast;
814 nh = (n1*(n1+1)) >> 1;
815 a1 = H = (double*)ADvari::adc.Memalloc(nh * sizeof(*h));
816 for(i = 0; i < nh; i++)
817 a1[i] = h[i];
818 }
819
820 ADvari&
821ADfn(double f, int n, const ADvar *x, const double *g, const double *h) {
822 return *(new ADvarn(f, n, x, g, h));
823 }
824
825 void
826ADcontext::Hvprod(int n, ADvar **x, double *v, double *hv)
827{
828 ADvari *a, *aL, *aR, **ap, **ape;
829 ADvari_block *b, *b0;
830 Derp *d;
831 double aO, adO, *g, *h, *h0, t, tL, tR;
832 int i, j, k, m;
833 for(i = 0; i < n; i++) {
834 a = x[i]->cv;
835 a->dO = v[i];
836 a->aO = a->adO = 0.;
837 }
839 a = 0;
840 for(b0 = 0, b = &ADvari::adc.AiFirst; b; b0 = b, b = b->next) {
841 ap = b->pADvari;
842 ape = b->limit;
843 while(ap < ape) {
844 a = *ap++;
845 a->aO = a->adO = 0.;
846 switch(a->opclass) {
847 case Hv_copy:
848 a->dO = ((ADvar1*)a)->d.c->dO;
849 break;
850 case Hv_binary:
851 a->dO = ((ADvar2g*)a)->pL * ((ADvar2g*)a)->dL.c->dO
852 + ((ADvar2g*)a)->pR * ((ADvar2g*)a)->dR.c->dO;
853 break;
854 case Hv_unary:
855 a->dO = ((ADvar1g*)a)->pL * ((ADvar1g*)a)->d.c->dO;
856 break;
857 case Hv_negate:
858 a->dO = -((ADvar1*)a)->d.c->dO;
859 break;
860 case Hv_plusLR:
861 a->dO = ((ADvar2*)a)->dL.c->dO + ((ADvar2*)a)->dR.c->dO;
862 break;
863 case Hv_minusLR:
864 a->dO = ((ADvar2*)a)->dL.c->dO - ((ADvar2*)a)->dR.c->dO;
865 break;
866 case Hv_timesL:
867 a->dO = ((ADvar1s*)a)->pL * ((ADvar1s*)a)->d.c->dO;
868 break;
869 case Hv_timesLR:
870 a->dO = ((ADvar2*)a)->dR.c->Val * ((ADvar2*)a)->dL.c->dO
871 + ((ADvar2*)a)->dL.c->Val * ((ADvar2*)a)->dR.c->dO;
872 break;
873 case Hv_quotLR:
874 a->dO = ((ADvar2q*)a)->pL * ((ADvar2q*)a)->dL.c->dO
875 + ((ADvar2q*)a)->pR * ((ADvar2q*)a)->dR.c->dO;
876 case Hv_const: // This case does not arise; the intent here
877 // is to eliminate a red-herring compiler warning.
878 break;
879 case Hv_nary:
880 d = ((ADvarn*)a)->D;
881 m = ((ADvarn*)a)->n;
882 g = ((ADvarn*)a)->G;
883 t = 0.;
884 for(i = 0; i < m; i++)
885 t += g[i] * d[i].c->dO;
886 a->dO = t;
887 }
888 }
889 }
890 if (a)
891 a->adO = 1.;
892 for(b = b0; b; b = b->prev) {
893 ape = b->pADvari;
894 ap = b->limit;
895 while(ap > ape) {
896 a = *--ap;
897 aO = a->aO;
898 adO = a->adO;
899 switch(a->opclass) {
900 case Hv_copy:
901 aL = ((ADvar1*)a)->d.c;
902 aL->aO += aO;
903 aL->adO += adO;
904 break;
905 case Hv_binary:
906 aL = ((ADvar2g*)a)->dL.c;
907 aR = ((ADvar2g*)a)->dR.c;
908 tL = adO*aL->dO;
909 tR = adO*aR->dO;
910 aL->aO += aO*((ADvar2g*)a)->pL
911 + tL*((ADvar2g*)a)->pL2
912 + tR*((ADvar2g*)a)->pLR;
913 aR->aO += aO*((ADvar2g*)a)->pR
914 + tL*((ADvar2g*)a)->pLR
915 + tR*((ADvar2g*)a)->pR2;
916 aL->adO += adO * ((ADvar2g*)a)->pL;
917 aR->adO += adO * ((ADvar2g*)a)->pR;
918 break;
919 case Hv_unary:
920 aL = ((ADvar1g*)a)->d.c;
921 aL->aO += aO * ((ADvar1g*)a)->pL
922 + adO * aL->dO * ((ADvar1g*)a)->pL2;
923 aL->adO += adO * ((ADvar1g*)a)->pL;
924 break;
925 case Hv_negate:
926 aL = ((ADvar1*)a)->d.c;
927 aL->aO -= aO;
928 aL->adO -= adO;
929 break;
930 case Hv_plusLR:
931 aL = ((ADvar2*)a)->dL.c;
932 aR = ((ADvar2*)a)->dR.c;
933 aL->aO += aO;
934 aL->adO += adO;
935 aR->aO += aO;
936 aR->adO += adO;
937 break;
938 case Hv_minusLR:
939 aL = ((ADvar2*)a)->dL.c;
940 aR = ((ADvar2*)a)->dR.c;
941 aL->aO += aO;
942 aL->adO += adO;
943 aR->aO -= aO;
944 aR->adO -= adO;
945 break;
946 case Hv_timesL:
947 aL = ((ADvar1s*)a)->d.c;
948 aL->aO += aO * (tL = ((ADvar1s*)a)->pL);
949 aL->adO += adO * tL;
950 break;
951 case Hv_timesLR:
952 aL = ((ADvar2*)a)->dL.c;
953 aR = ((ADvar2*)a)->dR.c;
954 aL->aO += aO * (tL = aR->Val) + adO*aR->dO;
955 aR->aO += aO * (tR = aL->Val) + adO*aL->dO;
956 aL->adO += adO * tL;
957 aR->adO += adO * tR;
958 break;
959 case Hv_quotLR:
960 aL = ((ADvar2q*)a)->dL.c;
961 aR = ((ADvar2q*)a)->dR.c;
962 tL = adO*aL->dO;
963 tR = adO*aR->dO;
964 aL->aO += aO*((ADvar2q*)a)->pL
965 + tR*((ADvar2q*)a)->pLR;
966 aR->aO += aO*((ADvar2q*)a)->pR
967 + tL*((ADvar2q*)a)->pLR
968 + tR*((ADvar2q*)a)->pR2;
969 aL->adO += adO * ((ADvar2q*)a)->pL;
970 aR->adO += adO * ((ADvar2q*)a)->pR;
971 case Hv_const: // This case does not arise; the intent here
972 // is to eliminate a red-herring compiler warning.
973 break;
974 case Hv_nary:
975 d = ((ADvarn*)a)->D;
976 m = ((ADvarn*)a)->n;
977 g = ((ADvarn*)a)->G;
978 h0 = ((ADvarn*)a)->H;
979 for(i = 0; i < m; i++) {
980 aL = d[i].c;
981 aL->adO += adO * (t = g[i]);
982 aL->aO += t*aO;
983 t = adO * aL->dO;
984 for(h = h0, j = 0; j <= i; j++)
985 d[j].c->aO += t * *h++;
986 h0 = h--;
987 for(k = j; j < m; j++)
988 d[j].c->aO += t * *(h += k++);
989 }
990 }
991 }
992 }
993 for(i = 0; i < n; i++) {
994 a = x[i]->cv;
995 a->dO = 0.;
996 hv[i] = a->aO;
997 }
998 }
999
1000#ifndef SACADO_NO_NAMESPACE
1001} // namespace Rad2d
1002} // namespace Sacado
1003#endif
fabs(expr.val())
asinh(expr.val())
log(expr.val())
tan(expr.val())
cos(expr.val())
expr expr dx(i)
cosh(expr.val())
acos(expr.val())
sin(expr.val())
sinh(expr.val())
log10(expr.val())
exp(expr.val())
atan(expr.val())
acosh(expr.val())
atan2(expr1.val(), expr2.val())
sqrt(expr.val())
tanh(expr.val())
atanh(expr.val())
asin(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 T
Definition: Sacado_rad.hpp:573
#define R
Definition: Sacado_rad.hpp:181
static void Weighted_Gradcomp(int, ADvar **, double *)
void * new_ADmemblock(size_t)
static void Hvprod(int, ADvar **, double *, double *)
ADvar1(Advari_Opclass oc, double val1, const double *a1, const ADvari *c1)
ADvar2g(double val1, double Lp, double Rp, double L2, double LR, double R2, const ADvari *Lcv, const ADvari *Rcv)
ADvar2q(double val1, double Lp, double Rp, double LR, double R2, const ADvari *Lcv, const ADvari *Rcv)
ADvar & operator-=(const ADvari &)
ADvar & operator=(const ADvari &x)
ADvar & operator/=(const ADvari &)
ADvar & operator*=(const ADvari &)
ADvar & operator+=(const ADvari &)
void ADvar_ctr(double d)
Advari_Opclass opclass
static ADcontext adc
ADvarn(double val1, int n1, const ADvar *x, const double *g, const double *h)
static const double negOne
static const double One
static ConstADvari * lastcad
static CADcontext cadc
static void aval_reset(void)
const ADvari * b
static Derp * LastDerp
const double * a
IndepADvar & operator=(const IndepADvar &x)
friend void AD_Const(const IndepADvar &)
const double y
const char * p
ADvari & ADf2(double f, double gx, double gy, double hxx, double hxy, double hyy, const ADvari &x, const ADvari &y)
ADvari & operator+(ADvari &T)
static int rad_need_reinit
ADvar & ADvar_operatoreq(ADvar *This, const ADvari &x)
ADvari & operator-(const ADvari &T)
ADvari & ADf1(double f, double g, double h, const ADvari &x)
ADvari & operator*(const ADvari &L, const ADvari &R)
ADvari & ADfn(double f, int n, const ADvar *x, const double *g, const double *h)
ADvari & operator/(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