NGSolve  5.3
autodiffdiff.hpp
1 #ifndef FILE_AUTODIFFDIFF
2 #define FILE_AUTODIFFDIFF
3 
4 /**************************************************************************/
5 /* File: autodiffdiff.hpp */
6 /* Author: Joachim Schoeberl */
7 /* Date: 13. June. 05 */
8 /**************************************************************************/
9 
10 namespace ngstd
11 {
12  using ngcore::IfPos;
13 
14 // Automatic second differentiation datatype
15 
16 
22 template <int D, typename SCAL = double>
24 {
25  SCAL val;
26  SCAL dval[D?D:1];
27  SCAL ddval[D?D*D:1];
28 public:
29 
31 
32 
34  AutoDiffDiff () throw() { ; }
35 
37  AutoDiffDiff (const AutoDiffDiff & ad2) throw()
38  {
39  val = ad2.val;
40  for (int i = 0; i < D; i++)
41  dval[i] = ad2.dval[i];
42  for (int i = 0; i < D*D; i++)
43  ddval[i] = ad2.ddval[i];
44  }
45 
47  AutoDiffDiff (SCAL aval) throw()
48  {
49  val = aval;
50  for (int i = 0; i < D; i++)
51  dval[i] = 0;
52  for (int i = 0; i < D*D; i++)
53  ddval[i] = 0;
54  }
55 
57  AutoDiffDiff (const AutoDiff<D, SCAL> & ad2) throw()
58  {
59  val = ad2.Value();
60  for (int i = 0; i < D; i++)
61  dval[i] = ad2.DValue(i);
62  for (int i = 0; i < D*D; i++)
63  ddval[i] = 0;
64  }
65 
67  AutoDiffDiff (SCAL aval, int diffindex) throw()
68  {
69  val = aval;
70  for (int i = 0; i < D; i++)
71  dval[i] = 0;
72  for (int i = 0; i < D*D; i++)
73  ddval[i] = 0;
74  dval[diffindex] = 1;
75  }
76 
77  INLINE AutoDiffDiff (SCAL aval, const SCAL * grad)
78  {
79  val = aval;
80  LoadGradient (grad);
81  for (int i = 0; i < D*D; i++)
82  ddval[i] = 0;
83  }
84 
85  INLINE AutoDiffDiff (SCAL aval, const SCAL * grad, const SCAL * hesse)
86  {
87  val = aval;
88  LoadGradient (grad);
89  LoadHessian (hesse);
90  }
91 
93  AutoDiffDiff & operator= (SCAL aval) throw()
94  {
95  val = aval;
96  for (int i = 0; i < D; i++)
97  dval[i] = 0;
98  for (int i = 0; i < D*D; i++)
99  ddval[i] = 0;
100  return *this;
101  }
102 
103  INLINE void StoreGradient (SCAL * p) const
104  {
105  for (int i = 0; i < D; i++)
106  p[i] = dval[i];
107  }
108 
109  INLINE void LoadGradient (const SCAL * p)
110  {
111  for (int i = 0; i < D; i++)
112  dval[i] = p[i];
113  }
114 
115  INLINE void StoreHessian (SCAL * p) const
116  {
117  for (int i = 0; i < D*D; i++)
118  p[i] = ddval[i];
119  }
120 
121  INLINE void LoadHessian (const SCAL * p)
122  {
123  for (int i = 0; i < D*D; i++)
124  ddval[i] = p[i];
125  }
126 
128  SCAL Value() const throw() { return val; }
129 
131  SCAL DValue (int i) const throw() { return dval[i]; }
132 
133  AutoDiff<D,SCAL> DValueAD (int i) const
134  {
135  AutoDiff<D,SCAL> r(dval[i]);
136  for (int j = 0; j < D; j++)
137  r.DValue(j) = ddval[i*D+j];
138  return r;
139  }
140 
142  SCAL DDValue (int i) const throw() { return ddval[i]; }
143 
145  SCAL DDValue (int i, int j) const throw() { return ddval[i*D+j]; }
146 
148  SCAL & Value() throw() { return val; }
149 
151  SCAL & DValue (int i) throw() { return dval[i]; }
152 
154  SCAL & DDValue (int i) throw() { return ddval[i]; }
155 
157  SCAL & DDValue (int i, int j) throw() { return ddval[i*D+j]; }
158 
159  explicit operator AutoDiff<D,SCAL> () const
160  { return AutoDiff<D,SCAL> (val, &dval[0]); }
161 
164  {
165  val += y.val;
166  for (int i = 0; i < D; i++)
167  dval[i] += y.dval[i];
168  for (int i = 0; i < D*D; i++)
169  ddval[i] += y.ddval[i];
170  return *this;
171  }
172 
175  {
176  val -= y.val;
177  for (int i = 0; i < D; i++)
178  dval[i] -= y.dval[i];
179  for (int i = 0; i < D*D; i++)
180  ddval[i] -= y.ddval[i];
181  return *this;
182  }
183 
186  {
187  for (int i = 0; i < D*D; i++)
188  ddval[i] = val * y.ddval[i] + y.val * ddval[i];
189 
190  for (int i = 0; i < D; i++)
191  for (int j = 0; j < D; j++)
192  ddval[i*D+j] += dval[i] * y.dval[j] + dval[j] * y.dval[i];
193 
194  for (int i = 0; i < D; i++)
195  {
196  dval[i] *= y.val;
197  dval[i] += val * y.dval[i];
198  }
199  val *= y.val;
200  return *this;
201  }
202 
204  AutoDiffDiff<D, SCAL> & operator*= (const SCAL & y) throw()
205  {
206  for ( int i = 0; i < D*D; i++ )
207  ddval[i] *= y;
208  for (int i = 0; i < D; i++)
209  dval[i] *= y;
210  val *= y;
211  return *this;
212  }
213 
215  AutoDiffDiff<D, SCAL> & operator/= (const SCAL & y) throw()
216  {
217  SCAL iy = 1.0 / y;
218  for ( int i = 0; i < D*D; i++ )
219  ddval[i] *= iy;
220  for (int i = 0; i < D; i++)
221  dval[i] *= iy;
222  val *= iy;
223  return *this;
224  }
225 
227  bool operator== (SCAL val2) throw()
228  {
229  return val == val2;
230  }
231 
233  bool operator!= (SCAL val2) throw()
234  {
235  return val != val2;
236  }
237 
239  bool operator< (SCAL val2) throw()
240  {
241  return val < val2;
242  }
243 
245  bool operator> (SCAL val2) throw()
246  {
247  return val > val2;
248  }
249 };
250 
251 
253 
255 template<int D, typename SCAL>
256 inline ostream & operator<< (ostream & ost, const AutoDiffDiff<D, SCAL> & x)
257 {
258  ost << x.Value() << ", D = ";
259  for (int i = 0; i < D; i++)
260  ost << x.DValue(i) << " ";
261  ost << ", DD = ";
262  for (int i = 0; i < D*D; i++)
263  ost << x.DDValue(i) << " ";
264  return ost;
265 }
266 
268 template<int D, typename SCAL>
269 inline AutoDiffDiff<D, SCAL> operator+ (const AutoDiffDiff<D, SCAL> & x, const AutoDiffDiff<D, SCAL> & y) throw()
270 {
271  AutoDiffDiff<D, SCAL> res;
272  res.Value () = x.Value()+y.Value();
273  for (int i = 0; i < D; i++)
274  res.DValue(i) = x.DValue(i) + y.DValue(i);
275  for (int i = 0; i < D*D; i++)
276  res.DDValue(i) = x.DDValue(i) + y.DDValue(i);
277  return res;
278 }
279 
280 
282 template<int D, typename SCAL>
283 inline AutoDiffDiff<D, SCAL> operator- (const AutoDiffDiff<D, SCAL> & x, const AutoDiffDiff<D, SCAL> & y) throw()
284 {
285  AutoDiffDiff<D, SCAL> res;
286  res.Value() = x.Value()-y.Value();
287  for (int i = 0; i < D; i++)
288  res.DValue(i) = x.DValue(i) - y.DValue(i);
289  for (int i = 0; i < D*D; i++)
290  res.DDValue(i) = x.DDValue(i) - y.DDValue(i);
291  return res;
292 }
293 
294 
296 template<int D, typename SCAL, typename SCAL2,
297  typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
298 inline AutoDiffDiff<D, SCAL> operator+ (SCAL2 x, const AutoDiffDiff<D, SCAL> & y) throw()
299 {
300  AutoDiffDiff<D, SCAL> res;
301  res.Value() = x+y.Value();
302  for (int i = 0; i < D; i++)
303  res.DValue(i) = y.DValue(i);
304  for (int i = 0; i < D*D; i++)
305  res.DDValue(i) = y.DDValue(i);
306  return res;
307 }
308 
310 template<int D, typename SCAL, typename SCAL2,
311  typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
312 inline AutoDiffDiff<D, SCAL> operator+ (const AutoDiffDiff<D, SCAL> & y, SCAL2 x) throw()
313 {
314  AutoDiffDiff<D, SCAL> res;
315  res.Value() = x+y.Value();
316  for (int i = 0; i < D; i++)
317  res.DValue(i) = y.DValue(i);
318  for (int i = 0; i < D*D; i++)
319  res.DDValue(i) = y.DDValue(i);
320  return res;
321 }
322 
323 
325 template<int D, typename SCAL>
326 inline AutoDiffDiff<D, SCAL> operator- (const AutoDiffDiff<D, SCAL> & x) throw()
327 {
328  AutoDiffDiff<D, SCAL> res;
329  res.Value() = -x.Value();
330  for (int i = 0; i < D; i++)
331  res.DValue(i) = -x.DValue(i);
332  for (int i = 0; i < D*D; i++)
333  res.DDValue(i) = -x.DDValue(i);
334  return res;
335 }
336 
338 template<int D, typename SCAL, typename SCAL2,
339  typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
340 inline AutoDiffDiff<D, SCAL> operator- (const AutoDiffDiff<D, SCAL> & x, SCAL2 y) throw()
341 {
342  AutoDiffDiff<D, SCAL> res;
343  res.Value() = x.Value()-y;
344  for (int i = 0; i < D; i++)
345  res.DValue(i) = x.DValue(i);
346  for (int i = 0; i < D*D; i++)
347  res.DDValue(i) = x.DDValue(i);
348  return res;
349 }
350 
352 template<int D, typename SCAL, typename SCAL2,
353  typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
354 inline AutoDiffDiff<D, SCAL> operator- (SCAL2 x, const AutoDiffDiff<D, SCAL> & y) throw()
355 {
356  AutoDiffDiff<D, SCAL> res;
357  res.Value() = x-y.Value();
358  for (int i = 0; i < D; i++)
359  res.DValue(i) = -y.DValue(i);
360  for (int i = 0; i < D*D; i++)
361  res.DDValue(i) = -y.DDValue(i);
362  return res;
363 }
364 
365 
367 template<int D, typename SCAL, typename SCAL2,
368  typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
369 inline AutoDiffDiff<D, SCAL> operator* (SCAL2 x, const AutoDiffDiff<D, SCAL> & y) throw()
370 {
371  AutoDiffDiff<D, SCAL> res;
372  res.Value() = x*y.Value();
373  for (int i = 0; i < D; i++)
374  res.DValue(i) = x*y.DValue(i);
375  for (int i = 0; i < D*D; i++)
376  res.DDValue(i) = x*y.DDValue(i);
377  return res;
378 }
379 
381 template<int D, typename SCAL, typename SCAL2,
382  typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
383 inline AutoDiffDiff<D, SCAL> operator* (const AutoDiffDiff<D, SCAL> & y, SCAL2 x) throw()
384 {
385  AutoDiffDiff<D, SCAL> res;
386  res.Value() = x*y.Value();
387  for (int i = 0; i < D; i++)
388  res.DValue(i) = x*y.DValue(i);
389  for (int i = 0; i < D*D; i++)
390  res.DDValue(i) = x*y.DDValue(i);
391  return res;
392 }
393 
395 template<int D, typename SCAL>
396 inline AutoDiffDiff<D, SCAL> operator* (const AutoDiffDiff<D, SCAL> & x, const AutoDiffDiff<D, SCAL> & y) throw()
397 {
398  AutoDiffDiff<D, SCAL> res;
399  SCAL hx = x.Value();
400  SCAL hy = y.Value();
401 
402  res.Value() = hx*hy;
403  for (int i = 0; i < D; i++)
404  res.DValue(i) = hx*y.DValue(i) + hy*x.DValue(i);
405 
406  for (int i = 0; i < D; i++)
407  for (int j = 0; j < D; j++)
408  res.DDValue(i,j) = hx * y.DDValue(i,j) + hy * x.DDValue(i,j)
409  + x.DValue(i) * y.DValue(j) + x.DValue(j) * y.DValue(i);
410 
411  return res;
412 }
413 
414 
415 
416 template<int D, typename SCAL>
417 inline AutoDiffDiff<D, SCAL> Inv (const AutoDiffDiff<D, SCAL> & x)
418 {
419  AutoDiffDiff<D, SCAL> res(1.0 / x.Value());
420  for (int i = 0; i < D; i++)
421  res.DValue(i) = -x.DValue(i) / (x.Value() * x.Value());
422 
423  SCAL fac1 = 2/(x.Value()*x.Value()*x.Value());
424  SCAL fac2 = 1/sqr(x.Value());
425  for (int i = 0; i < D; i++)
426  for (int j = 0; j < D; j++)
427  res.DDValue(i,j) = fac1*x.DValue(i)*x.DValue(j) - fac2*x.DDValue(i,j);
428  return res;
429 }
430 
431 
432 template<int D, typename SCAL>
433 inline AutoDiffDiff<D, SCAL> operator/ (const AutoDiffDiff<D, SCAL> & x, const AutoDiffDiff<D, SCAL> & y)
434 {
435  return x * Inv (y);
436 }
437 
438 template<int D, typename SCAL, typename SCAL2,
439  typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
440 inline AutoDiffDiff<D, SCAL> operator/ (const AutoDiffDiff<D, SCAL> & x, SCAL2 y)
441 {
442  return (1/y) * x;
443 }
444 
445 template<int D, typename SCAL, typename SCAL2,
446  typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
447 inline AutoDiffDiff<D, SCAL> operator/ (SCAL2 x, const AutoDiffDiff<D, SCAL> & y)
448 {
449  return x * Inv(y);
450 }
451 
452 
453 template<int D, typename SCAL>
454 inline AutoDiffDiff<D, SCAL> sqrt (const AutoDiffDiff<D, SCAL> & x)
455 {
456  AutoDiffDiff<D, SCAL> res;
457  res.Value() = sqrt(x.Value());
458  for (int j = 0; j < D; j++)
459  res.DValue(j) = IfZero(x.DValue(j),SCAL{0.},0.5 / res.Value() * x.DValue(j));
460 
461 
462  for (int i = 0; i < D; i++)
463  for (int j = 0; j < D; j++)
464  res.DDValue(i,j) = IfZero(x.DDValue(i,j)+x.DValue(i) * x.DValue(j),SCAL{0.},0.5/res.Value() * x.DDValue(i,j) - 0.25 / (x.Value()*res.Value()) * x.DValue(i) * x.DValue(j));
465 
466  return res;
467 }
468 
469 // df(u)/dx = exp(x) * du/dx
470 // d^2 f(u) / dx^2 = exp(x) * (du/dx)^2 + exp(x) * d^2u /dx^2
471 template <int D, typename SCAL>
472 INLINE AutoDiffDiff<D, SCAL> exp (AutoDiffDiff<D, SCAL> x)
473 {
474  AutoDiffDiff<D, SCAL> res;
475  res.Value() = exp(x.Value());
476  for (int k = 0; k < D; k++)
477  res.DValue(k) = x.DValue(k) * res.Value();
478  for (int k = 0; k < D; k++)
479  for (int l = 0; l < D; l++)
480  res.DDValue(k,l) = (x.DValue(k) * x.DValue(l)+x.DDValue(k,l)) * res.Value();
481  return res;
482 }
483 
484 using std::pow;
485 template <int D, typename SCAL>
486 INLINE AutoDiffDiff<D,SCAL> pow (AutoDiffDiff<D,SCAL> x, AutoDiffDiff<D,SCAL> y )
487 {
488  return exp(log(x)*y);
489 }
490 
491 template <int D, typename SCAL>
492 INLINE AutoDiffDiff<D, SCAL> log (AutoDiffDiff<D, SCAL> x)
493 {
494  AutoDiffDiff<D, SCAL> res;
495  res.Value() = log(x.Value());
496  SCAL xinv = 1.0/x.Value();
497  for (int k = 0; k < D; k++)
498  res.DValue(k) = x.DValue(k) * xinv;
499  for (int k = 0; k < D; k++)
500  for (int l = 0; l < D; l++)
501  res.DDValue(k,l) = -xinv*xinv*x.DValue(k) * x.DValue(l) + xinv * x.DDValue(k,l);
502  return res;
503 }
504 
505 
506 
507 template <int D, typename SCAL>
508 INLINE AutoDiffDiff<D, SCAL> sin (AutoDiffDiff<D, SCAL> x)
509 {
510  AutoDiffDiff<D, SCAL> res;
511  SCAL s = sin(x.Value());
512  SCAL c = cos(x.Value());
513 
514  res.Value() = s;
515  for (int k = 0; k < D; k++)
516  res.DValue(k) = x.DValue(k) * c;
517  for (int k = 0; k < D; k++)
518  for (int l = 0; l < D; l++)
519  res.DDValue(k,l) = -s * x.DValue(k) * x.DValue(l) + c * x.DDValue(k,l);
520  return res;
521 }
522 
523 
524 template <int D, typename SCAL>
525 INLINE AutoDiffDiff<D, SCAL> cos (AutoDiffDiff<D, SCAL> x)
526 {
527  AutoDiffDiff<D, SCAL> res;
528  SCAL s = sin(x.Value());
529  SCAL c = cos(x.Value());
530 
531  res.Value() = c;
532  for (int k = 0; k < D; k++)
533  res.DValue(k) = -s * x.DValue(k);
534  for (int k = 0; k < D; k++)
535  for (int l = 0; l < D; l++)
536  res.DDValue(k,l) = -c * x.DValue(k) * x.DValue(l) - s * x.DDValue(k,l);
537  return res;
538 }
539 
540 template <int D, typename SCAL>
541 INLINE AutoDiffDiff<D, SCAL> tan (AutoDiffDiff<D, SCAL> x)
542 { return sin(x) / cos(x); }
543 
544 
545 template <int D, typename SCAL>
546 INLINE AutoDiffDiff<D, SCAL> atan (AutoDiffDiff<D, SCAL> x)
547 {
548  AutoDiffDiff<D, SCAL> res;
549  SCAL a = atan(x.Value());
550  res.Value() = a;
551  for (int k = 0; k < D; k++)
552  res.DValue(k) = x.DValue(k)/(1+x.Value()*x.Value()) ;
553  for (int k = 0; k < D; k++)
554  for (int l = 0; l < D; l++)
555  res.DDValue(k,l) = -2*x.Value()/((1+x.Value()*x.Value())*(1+x.Value()*x.Value())) * x.DValue(k) * x.DValue(l) + x.DDValue(k,l)/(1+x.Value()*x.Value());
556  return res;
557 }
558 
559 template <int D, typename SCAL>
560 INLINE AutoDiffDiff<D, SCAL> atan2 (AutoDiffDiff<D, SCAL> x,AutoDiffDiff<D, SCAL> y)
561 {
562  AutoDiffDiff<D, SCAL> res;
563  SCAL a = atan2(x.Value(), y.Value());
564  res.Value() = a;
565  for (int k = 0; k < D; k++)
566  res.DValue(k) = (x.Value()*y.DValue(k)-y.Value()*x.DValue(k))/(y.Value()*y.Value()+x.Value()*x.Value());
567 
568  for (int k = 0; k < D; k++)
569  for (int l = 0; l < D; l++)
570  res.DDValue(k,l) = (x.DValue(k)*y.DValue(l)+x.Value()*y.DDValue(l,k) - y.DValue(k)*x.DValue(l) - y.Value()*x.DDValue(l,k))/(y.Value()*y.Value()+x.Value()*x.Value()) - 2 * (x.Value()*y.DValue(k)-y.Value()*x.DValue(k)) * (x.Value()*x.DValue(k) + y.Value()*y.DValue(k))/( (y.Value()*y.Value()+x.Value()*x.Value()) * (y.Value()*y.Value()+x.Value()*x.Value()) );
571  return res;
572 }
573 
574 
575 
576 using std::acos;
577 template <int D, typename SCAL>
578 INLINE AutoDiffDiff<D,SCAL> acos (AutoDiffDiff<D,SCAL> x)
579 {
580  AutoDiffDiff<D,SCAL> res;
581  SCAL a = acos(x.Value());
582  res.Value() = a;
583  auto omaa = 1-x.Value()*x.Value();
584  auto s = sqrt(omaa);
585  SCAL da = -1 / s;
586  SCAL dda = -x.Value() / (s*omaa);
587  for (int k = 0; k < D; k++)
588  res.DValue(k) = x.DValue(k)*da;
589  for (int k = 0; k < D; k++)
590  for (int l = 0; l < D; l++)
591  res.DDValue(k,l) = dda * x.DValue(k) * x.DValue(l) + da * x.DDValue(k,l);
592 
593  return res;
594 }
595 
596 
597 using std::acos;
598 template <int D, typename SCAL>
599 INLINE AutoDiffDiff<D,SCAL> asin (AutoDiffDiff<D,SCAL> x)
600 {
601  AutoDiffDiff<D,SCAL> res;
602  SCAL a = asin(x.Value());
603  res.Value() = a;
604  auto omaa = 1-x.Value()*x.Value();
605  auto s = sqrt(omaa);
606  SCAL da = 1 / s;
607  SCAL dda = x.Value() / (s*omaa);
608  for (int k = 0; k < D; k++)
609  res.DValue(k) = x.DValue(k)*da;
610  for (int k = 0; k < D; k++)
611  for (int l = 0; l < D; l++)
612  res.DDValue(k,l) = dda * x.DValue(k) * x.DValue(l) + da * x.DDValue(k,l);
613 
614  return res;
615 }
616 
617 
618 template <int D, typename SCAL>
619 INLINE AutoDiffDiff<D, SCAL> sinh (AutoDiffDiff<D, SCAL> x)
620 {
621  AutoDiffDiff<D, SCAL> res;
622  SCAL sh = sinh(x.Value());
623  SCAL ch = cosh(x.Value());
624 
625  res.Value() = sh;
626  for (int k = 0; k < D; k++)
627  res.DValue(k) = x.DValue(k) * ch;
628  for (int k = 0; k < D; k++)
629  for (int l = 0; l < D; l++)
630  res.DDValue(k,l) = sh * x.DValue(k) * x.DValue(l) + ch * x.DDValue(k,l);
631  return res;
632 }
633 
634 
635 template <int D, typename SCAL>
636 INLINE AutoDiffDiff<D, SCAL> cosh (AutoDiffDiff<D, SCAL> x)
637 {
638  AutoDiffDiff<D, SCAL> res;
639  SCAL sh = sinh(x.Value());
640  SCAL ch = cosh(x.Value());
641 
642  res.Value() = ch;
643  for (int k = 0; k < D; k++)
644  res.DValue(k) = sh * x.DValue(k);
645  for (int k = 0; k < D; k++)
646  for (int l = 0; l < D; l++)
647  res.DDValue(k,l) = ch * x.DValue(k) * x.DValue(l) + sh * x.DDValue(k,l);
648  return res;
649 }
650 
651 using std::floor;
652 template<int D, typename SCAL>
653 INLINE AutoDiffDiff<D,SCAL> floor (const AutoDiffDiff<D,SCAL> & x)
654 {
655  return floor(x.Value());
656 }
657 
658 using std::ceil;
659 template<int D, typename SCAL>
660 INLINE AutoDiffDiff<D,SCAL> ceil (const AutoDiffDiff<D,SCAL> & x)
661 {
662  return ceil(x.Value());
663 }
664 
665 
666 template <int D, typename SCAL, typename TB, typename TC>
667 auto IfPos (AutoDiffDiff<D,SCAL> a, TB b, TC c) -> decltype(IfPos (a.Value(), b, c))
668 {
669  return IfPos (a.Value(), b, c);
670 }
671 
672 template <int D, typename SCAL>
673 INLINE AutoDiffDiff<D,SCAL> IfPos (SCAL /* SIMD<double> */ a, AutoDiffDiff<D,SCAL> b, AutoDiffDiff<D,SCAL> c)
674 {
675  AutoDiffDiff<D,SCAL> res;
676  res.Value() = IfPos (a, b.Value(), c.Value());
677  for (int j = 0; j < D; j++)
678  {
679  res.DValue(j) = IfPos (a, b.DValue(j), c.DValue(j));
680  res.DDValue(j) = IfPos (a, b.DDValue(j), c.DDValue(j));
681  }
682  return res;
683 }
684 
685 template <int D, typename SCAL, typename TC>
686 INLINE AutoDiffDiff<D,SCAL> IfPos (SCAL /* SIMD<double> */ a, AutoDiffDiff<D,SCAL> b, TC c)
687 {
688  return IfPos (a, b, AutoDiffDiff<D,SCAL> (c));
689 }
690 
691 
692 
693 
695 
696 }
697 
698 
699 #endif
ngstd::AutoDiffDiff::AutoDiffDiff
AutoDiffDiff(const AutoDiffDiff &ad2)
copy constructor
Definition: autodiffdiff.hpp:37
ngstd::AutoDiffDiff::DDValue
SCAL DDValue(int i, int j) const
returns partial derivative
Definition: autodiffdiff.hpp:145
ngstd::AutoDiffDiff::DDValue
SCAL & DDValue(int i)
accesses partial derivative
Definition: autodiffdiff.hpp:154
ngstd::AutoDiffDiff::operator!=
bool operator!=(SCAL val2)
different values
Definition: autodiffdiff.hpp:233
ngstd::AutoDiffDiff::AutoDiffDiff
AutoDiffDiff(const AutoDiff< D, SCAL > &ad2)
initial object with value and derivative
Definition: autodiffdiff.hpp:57
ngstd::AutoDiffDiff::operator>
bool operator>(SCAL val2)
greater
Definition: autodiffdiff.hpp:245
ngstd::AutoDiffDiff::operator+=
AutoDiffDiff< D, SCAL > & operator+=(const AutoDiffDiff< D, SCAL > &y)
add autodiffdiff object
Definition: autodiffdiff.hpp:163
ngstd::AutoDiffRec
Definition: autodiff.hpp:16
ngstd::AutoDiffDiff::Value
SCAL Value() const
returns value
Definition: autodiffdiff.hpp:128
ngstd::AutoDiffDiff::AutoDiffDiff
AutoDiffDiff(SCAL aval)
initial object with constant value
Definition: autodiffdiff.hpp:47
ngstd::AutoDiffDiff::AutoDiffDiff
AutoDiffDiff(SCAL aval, int diffindex)
init object with (val, e_diffindex)
Definition: autodiffdiff.hpp:67
ngstd::AutoDiffDiff::operator-=
AutoDiffDiff< D, SCAL > & operator-=(const AutoDiffDiff< D, SCAL > &y)
subtract autodiffdiff object
Definition: autodiffdiff.hpp:174
ngstd::AutoDiffDiff::DDValue
SCAL DDValue(int i) const
returns partial derivative
Definition: autodiffdiff.hpp:142
ngstd::AutoDiffDiff::operator==
bool operator==(SCAL val2)
same value
Definition: autodiffdiff.hpp:227
ngstd
namespace for standard data types and algorithms.
Definition: ngstd.hpp:55
ngstd::AutoDiffDiff::DDValue
SCAL & DDValue(int i, int j)
accesses partial derivative
Definition: autodiffdiff.hpp:157
ngstd::AutoDiffDiff::AutoDiffDiff
AutoDiffDiff()
elements are undefined
Definition: autodiffdiff.hpp:34
ngstd::operator<<
ostream & operator<<(ostream &ost, const AutoDiffDiff< D, SCAL > &x)
Prints AudoDiffDiff.
Definition: autodiffdiff.hpp:256
ngstd::AutoDiffDiff
Datatype for automatic differentiation.
Definition: autodiffdiff.hpp:23
ngstd::AutoDiffDiff::operator=
AutoDiffDiff & operator=(SCAL aval)
assign constant value
Definition: autodiffdiff.hpp:93
ngstd::AutoDiffDiff::DValue
SCAL & DValue(int i)
accesses partial derivative
Definition: autodiffdiff.hpp:151
ngstd::AutoDiffDiff::DValue
SCAL DValue(int i) const
returns partial derivative
Definition: autodiffdiff.hpp:131
ngstd::AutoDiffDiff::operator<
bool operator<(SCAL val2)
less
Definition: autodiffdiff.hpp:239
ngstd::sqr
INLINE AutoDiffVec< D, SCAL > sqr(const AutoDiffVec< D, SCAL > &x)
AutoDiffVec times AutoDiffVec.
Definition: autodiff.hpp:251
ngstd::AutoDiffDiff::operator*=
AutoDiffDiff< D, SCAL > & operator*=(const AutoDiffDiff< D, SCAL > &y)
multiply with autodiffdiff object
Definition: autodiffdiff.hpp:185
ngstd::AutoDiffDiff::Value
SCAL & Value()
access value
Definition: autodiffdiff.hpp:148
ngstd::AutoDiffDiff::operator/=
AutoDiffDiff< D, SCAL > & operator/=(const SCAL &y)
divide by scalar
Definition: autodiffdiff.hpp:215