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