1 /* ---------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3
4 Revision: $Id: adoublecuda.h 709 2016-08-09 12:51:09Z kulshres $
5 Contents: adoublecuda.h contains the class of adouble specifically
6 suited to be used within CUDA environment
7
8 Copyright (c) Alina Koniaeva, Andrea Walther
9
10 This file is part of ADOL-C. This software is provided as open source.
11 Any use, reproduction, or distribution of the software constitutes
12 recipient's acceptance of the terms of the accompanying license file.
13
14 ---------------------------------------------------------------------------*/
15
16 #if !defined(ADOLC_ADOUBLECUDA_H)
17 #define ADOLC_ADOUBLECUDA_H 1
18
19 #include <cstdlib>
20 #include <iostream>
21 #include <cmath>
22 using std::cout;
23 using std::ostream;
24 using std::istream;
25
26 #include <cuda_runtime.h>
27 #include <math_constants.h>
28
29 namespace adtlc {
30
31
32 #if defined(NUMBER_DIRECTIONS)
33 __managed__ size_t ADOLC_numDir = NUMBER_DIRECTIONS;
34 # if defined(DYNAMIC_DIRECTIONS)
35 # define ADVAL_DECL *adval
36 # define ADVAL_TYPE_ADV const double* adv
37 # else
38 # define ADVAL_DECL adval[NUMBER_DIRECTIONS]
39 # define ADVAL_TYPE_ADV const double adv[NUMBER_DIRECTIONS]
40 # endif
41 # define ADVAL_TYPE double*
42 # define FOR_I_EQ_0_LT_NUMDIR for (size_t i = 0; i < ADOLC_numDir; i++)
43 # define ADVAL_I adval[i]
44 # define ADV_I adv[i]
45 # define V_I v[i]
46 #else
47 # define ADVAL_DECL adval
48 # define ADVAL adval
49 # define ADVAL_TYPE_ADV double adv
50 # define ADVAL_TYPE double
51 # define FOR_I_EQ_0_LT_NUMDIR
52 # define ADVAL_I adval
53 # define ADV_I adv
54 # define V_I v
55 #endif
56
57
58 #define ADOLC_MATH_NSP std
59
makeNaN()60 inline __device__ double makeNaN() {
61 return CUDART_NAN;
62 }
63
makeInf()64 inline __device__ double makeInf() {
65 return CUDART_INF;
66 }
67
68
69 #define CUDADEV __device__ inline
70 #define CUDAHOST __host__ inline
71 #define CUDAHOSTDEV __host__ __device__ inline
72
73 class adouble {
74 public:
75 // ctors
76 CUDADEV adouble();
77 CUDADEV adouble(const double v);
78 CUDADEV adouble(const double v, ADVAL_TYPE_ADV);
79 CUDADEV adouble(const adouble& a);
80 #if defined(DYNAMIC_DIRECTIONS)
81 CUDADEV ~adouble();
82 #endif
83 /******************* temporary results ******************************/
84 // sign
85 CUDADEV adouble operator - () const;
86 CUDADEV adouble operator + () const;
87
88 // addition
89 CUDADEV adouble operator + (const double v) const;
90 CUDADEV adouble operator + (const adouble& a) const;
91 CUDADEV friend
92 adouble operator + (const double v, const adouble& a);
93
94 // substraction
95 CUDADEV adouble operator - (const double v) const;
96 CUDADEV adouble operator - (const adouble& a) const;
97 CUDADEV friend
98 adouble operator - (const double v, const adouble& a);
99
100 // multiplication
101 CUDADEV adouble operator * (const double v) const;
102 CUDADEV adouble operator * (const adouble& a) const;
103 CUDADEV friend
104 adouble operator * (const double v, const adouble& a);
105
106 // division
107 CUDADEV adouble operator / (const double v) const;
108 CUDADEV adouble operator / (const adouble& a) const;
109 CUDADEV friend
110 adouble operator / (const double v, const adouble& a);
111
112 // inc/dec
113 CUDADEV adouble operator ++ ();
114 CUDADEV adouble operator ++ (int);
115 CUDADEV adouble operator -- ();
116 CUDADEV adouble operator -- (int);
117
118 // functions
119 CUDADEV friend adouble tan(const adouble &a);
120 CUDADEV friend adouble exp(const adouble &a);
121 CUDADEV friend adouble log(const adouble &a);
122 CUDADEV friend adouble sqrt(const adouble &a);
123 CUDADEV friend adouble sin(const adouble &a);
124 CUDADEV friend adouble cos(const adouble &a);
125 CUDADEV friend adouble asin(const adouble &a);
126 CUDADEV friend adouble acos(const adouble &a);
127 CUDADEV friend adouble atan(const adouble &a);
128
129 CUDADEV friend adouble atan2(const adouble &a, const adouble &b);
130 CUDADEV friend adouble pow(const adouble &a, double v);
131 CUDADEV friend adouble pow(const adouble &a, const adouble &b);
132 CUDADEV friend adouble pow(double v, const adouble &a);
133 CUDADEV friend adouble log10(const adouble &a);
134
135 CUDADEV friend adouble sinh (const adouble &a);
136 CUDADEV friend adouble cosh (const adouble &a);
137 CUDADEV friend adouble tanh (const adouble &a);
138 #if defined(ATRIG_ERF)
139 CUDADEV friend adouble asinh (const adouble &a);
140 CUDADEV friend adouble acosh (const adouble &a);
141 CUDADEV friend adouble atanh (const adouble &a);
142 #endif
143 CUDADEV friend adouble fabs (const adouble &a);
144 CUDADEV friend adouble ceil (const adouble &a);
145 CUDADEV friend adouble floor (const adouble &a);
146 CUDADEV friend adouble fmax (const adouble &a, const adouble &b);
147 CUDADEV friend adouble fmax (double v, const adouble &a);
148 CUDADEV friend adouble fmax (const adouble &a, double v);
149 CUDADEV friend adouble fmin (const adouble &a, const adouble &b);
150 CUDADEV friend adouble fmin (double v, const adouble &a);
151 CUDADEV friend adouble fmin (const adouble &a, double v);
152 CUDADEV friend adouble ldexp (const adouble &a, const adouble &b);
153 CUDADEV friend adouble ldexp (const adouble &a, const double v);
154 CUDADEV friend adouble ldexp (const double v, const adouble &a);
155 CUDADEV friend double frexp (const adouble &a, int* v);
156 #if defined(ATRIG_ERF)
157 CUDADEV friend adouble erf (const adouble &a);
158 #endif
159
160
161 /******************* nontemporary results ***************************/
162 // assignment
163 CUDADEV void operator = (const double v);
164 CUDADEV void operator = (const adouble& a);
165
166 // addition
167 CUDADEV void operator += (const double v);
168 CUDADEV void operator += (const adouble& a);
169
170 // substraction
171 CUDADEV void operator -= (const double v);
172 CUDADEV void operator -= (const adouble& a);
173
174 // multiplication
175 CUDADEV void operator *= (const double v);
176 CUDADEV void operator *= (const adouble& a);
177
178 // division
179 CUDADEV void operator /= (const double v);
180 CUDADEV void operator /= (const adouble& a);
181
182 // not
183 CUDADEV int operator ! () const;
184
185 // comparision
186 CUDADEV int operator != (const adouble&) const;
187 CUDADEV int operator != (const double) const;
188 CUDADEV friend int operator != (const double, const adouble&);
189
190 CUDADEV int operator == (const adouble&) const;
191 CUDADEV int operator == (const double) const;
192 CUDADEV friend int operator == (const double, const adouble&);
193
194 CUDADEV int operator <= (const adouble&) const;
195 CUDADEV int operator <= (const double) const;
196 CUDADEV friend int operator <= (const double, const adouble&);
197
198 CUDADEV int operator >= (const adouble&) const;
199 CUDADEV int operator >= (const double) const;
200 CUDADEV friend int operator >= (const double, const adouble&);
201
202 CUDADEV int operator > (const adouble&) const;
203 CUDADEV int operator > (const double) const;
204 CUDADEV friend int operator > (const double, const adouble&);
205
206 CUDADEV int operator < (const adouble&) const;
207 CUDADEV int operator < (const double) const;
208 CUDADEV friend int operator < (const double, const adouble&);
209
210 /******************* getter / setter ********************************/
211 CUDAHOSTDEV double getValue() const;
212 CUDAHOSTDEV void setValue(const double v);
213 CUDAHOSTDEV ADVAL_TYPE getADValue() const;
214 CUDAHOSTDEV void setADValue(ADVAL_TYPE v);
215 #if defined(NUMBER_DIRECTIONS)
216 CUDAHOSTDEV double getADValue(const unsigned int p) const;
217 CUDAHOSTDEV void setADValue(const unsigned int p, const double v);
218 #endif
219
220 /******************* i/o operations *********************************/
221 CUDAHOST friend ostream& operator << ( ostream&, const adouble& );
222 CUDAHOST friend istream& operator >> ( istream&, adouble& );
223
224 private:
225 // internal variables
226 double val;
227 double ADVAL_DECL;
228 };
229
230 /******************************* ctors ************************************/
adouble()231 CUDADEV adouble::adouble() {
232 #if defined(DYNAMIC_DIRECTIONS)
233 adval = new double[ADOLC_numDir];
234 #endif
235 }
236
adouble(const double v)237 CUDADEV adouble::adouble(const double v) : val(v) {
238 #if defined(DYNAMIC_DIRECTIONS)
239 adval = new double[ADOLC_numDir];
240 #endif
241 FOR_I_EQ_0_LT_NUMDIR
242 ADVAL_I = 0.0;
243 }
244
adouble(const double v,ADVAL_TYPE_ADV)245 CUDADEV adouble::adouble(const double v, ADVAL_TYPE_ADV) : val(v) {
246 #if defined(DYNAMIC_DIRECTIONS)
247 adval = new double[ADOLC_numDir];
248 #endif
249 FOR_I_EQ_0_LT_NUMDIR
250 ADVAL_I=ADV_I;
251 }
252
adouble(const adouble & a)253 CUDADEV adouble::adouble(const adouble& a) : val(a.val) {
254 #if defined(DYNAMIC_DIRECTIONS)
255 adval = new double[ADOLC_numDir];
256 #endif
257 FOR_I_EQ_0_LT_NUMDIR
258 ADVAL_I=a.ADVAL_I;
259 }
260
261 /******************************* dtors ************************************/
262 #if defined(DYNAMIC_DIRECTIONS)
~adouble()263 CUDADEV adouble::~adouble() {
264 delete[] adval;
265 }
266 #endif
267
268 /************************* temporary results ******************************/
269 // sign
270 CUDADEV adouble adouble::operator - () const {
271 adouble tmp;
272 tmp.val=-val;
273 FOR_I_EQ_0_LT_NUMDIR
274 tmp.ADVAL_I=-ADVAL_I;
275 return tmp;
276 }
277
278 CUDADEV adouble adouble::operator + () const {
279 return *this;
280 }
281
282 // addition
283 CUDADEV adouble adouble::operator + (const double v) const {
284 return adouble(val+v, adval);
285 }
286
287 CUDADEV adouble adouble::operator + (const adouble& a) const {
288 adouble tmp;
289 tmp.val=val+a.val;
290 FOR_I_EQ_0_LT_NUMDIR
291 tmp.ADVAL_I=ADVAL_I+a.ADVAL_I;
292 return tmp;
293 }
294
295 CUDADEV adouble operator + (const double v, const adouble& a) {
296 return adouble(v+a.val, a.adval);
297 }
298
299 // subtraction
300 CUDADEV adouble adouble::operator - (const double v) const {
301 return adouble(val-v, adval);
302 }
303
304 CUDADEV adouble adouble::operator - (const adouble& a) const {
305 adouble tmp;
306 tmp.val=val-a.val;
307 FOR_I_EQ_0_LT_NUMDIR
308 tmp.ADVAL_I=ADVAL_I-a.ADVAL_I;
309 return tmp;
310 }
311
312 CUDADEV adouble operator - (const double v, const adouble& a) {
313 adouble tmp;
314 tmp.val=v-a.val;
315 FOR_I_EQ_0_LT_NUMDIR
316 tmp.ADVAL_I=-a.ADVAL_I;
317 return tmp;
318 }
319
320 // multiplication
321 CUDADEV adouble adouble::operator * (const double v) const {
322 adouble tmp;
323 tmp.val=val*v;
324 FOR_I_EQ_0_LT_NUMDIR
325 tmp.ADVAL_I=ADVAL_I*v;
326 return tmp;
327 }
328
329 CUDADEV adouble adouble::operator * (const adouble& a) const {
330 adouble tmp;
331 tmp.val=val*a.val;
332 FOR_I_EQ_0_LT_NUMDIR
333 tmp.ADVAL_I=ADVAL_I*a.val+val*a.ADVAL_I;
334 return tmp;
335 }
336
337 CUDADEV adouble operator * (const double v, const adouble& a) {
338 adouble tmp;
339 tmp.val=v*a.val;
340 FOR_I_EQ_0_LT_NUMDIR
341 tmp.ADVAL_I=v*a.ADVAL_I;
342 return tmp;
343 }
344
345 // division
346 CUDADEV adouble adouble::operator / (const double v) const {
347 adouble tmp;
348 tmp.val=val/v;
349 FOR_I_EQ_0_LT_NUMDIR
350 tmp.ADVAL_I=ADVAL_I/v;
351 return tmp;
352 }
353
354 CUDADEV adouble adouble::operator / (const adouble& a) const {
355 adouble tmp;
356 tmp.val=val/a.val;
357 FOR_I_EQ_0_LT_NUMDIR
358 tmp.ADVAL_I=(ADVAL_I*a.val-val*a.ADVAL_I)/(a.val*a.val);
359 return tmp;
360 }
361
362 CUDADEV adouble operator / (const double v, const adouble& a) {
363 adouble tmp;
364 tmp.val=v/a.val;
365 FOR_I_EQ_0_LT_NUMDIR
366 tmp.ADVAL_I=(-v*a.ADVAL_I)/(a.val*a.val);
367 return tmp;
368 }
369
370 // inc/dec
371 CUDADEV adouble adouble::operator ++ () {
372 ++val;
373 return *this;
374 }
375
376 CUDADEV adouble adouble::operator ++ (int) {
377 adouble tmp;
378 tmp.val=val++;
379 FOR_I_EQ_0_LT_NUMDIR
380 tmp.ADVAL_I=ADVAL_I;
381 return tmp;
382 }
383
384 CUDADEV adouble adouble::operator -- () {
385 --val;
386 return *this;
387 }
388
389 CUDADEV adouble adouble::operator -- (int) {
390 adouble tmp;
391 tmp.val=val--;
392 FOR_I_EQ_0_LT_NUMDIR
393 tmp.ADVAL_I=ADVAL_I;
394 return tmp;
395 }
396
397 // functions
tan(const adouble & a)398 CUDADEV adouble tan(const adouble& a) {
399 adouble tmp;
400 double tmp2;
401 tmp.val=ADOLC_MATH_NSP::tan(a.val);
402 tmp2=ADOLC_MATH_NSP::cos(a.val);
403 tmp2*=tmp2;
404 FOR_I_EQ_0_LT_NUMDIR
405 tmp.ADVAL_I=a.ADVAL_I/tmp2;
406 return tmp;
407 }
408
exp(const adouble & a)409 CUDADEV adouble exp(const adouble &a) {
410 adouble tmp;
411 tmp.val=ADOLC_MATH_NSP::exp(a.val);
412 FOR_I_EQ_0_LT_NUMDIR
413 tmp.ADVAL_I=tmp.val*a.ADVAL_I;
414 return tmp;
415 }
416
log(const adouble & a)417 CUDADEV adouble log(const adouble &a) {
418 adouble tmp;
419 tmp.val=ADOLC_MATH_NSP::log(a.val);
420 FOR_I_EQ_0_LT_NUMDIR
421 if (a.val>0) tmp.ADVAL_I=a.ADVAL_I/a.val;
422 else if (a.val==0 && a.ADVAL_I != 0.0) {
423 int sign = (a.ADVAL_I < 0) ? -1 : 1;
424 tmp.ADVAL_I=sign* makeInf();
425 } else tmp.ADVAL_I=makeNaN();
426 return tmp;
427 }
428
sqrt(const adouble & a)429 CUDADEV adouble sqrt(const adouble &a) {
430 adouble tmp;
431 tmp.val=ADOLC_MATH_NSP::sqrt(a.val);
432 FOR_I_EQ_0_LT_NUMDIR
433 if (a.val>0) tmp.ADVAL_I=a.ADVAL_I/(tmp.val*2);
434 else if (a.val==0.0 && a.ADVAL_I != 0.0) {
435 int sign = (a.ADVAL_I < 0) ? -1 : 1;
436 tmp.ADVAL_I=sign * makeInf();
437 }
438 else tmp.ADVAL_I=makeNaN();
439 return tmp;
440 }
441
sin(const adouble & a)442 CUDADEV adouble sin(const adouble &a) {
443 adouble tmp;
444 double tmp2;
445 tmp.val=ADOLC_MATH_NSP::sin(a.val);
446 tmp2=ADOLC_MATH_NSP::cos(a.val);
447 FOR_I_EQ_0_LT_NUMDIR
448 tmp.ADVAL_I=tmp2*a.ADVAL_I;
449 return tmp;
450 }
451
cos(const adouble & a)452 CUDADEV adouble cos(const adouble &a) {
453 adouble tmp;
454 double tmp2;
455 tmp.val=ADOLC_MATH_NSP::cos(a.val);
456 tmp2=-ADOLC_MATH_NSP::sin(a.val);
457 FOR_I_EQ_0_LT_NUMDIR
458 tmp.ADVAL_I=tmp2*a.ADVAL_I;
459 return tmp;
460 }
461
asin(const adouble & a)462 CUDADEV adouble asin(const adouble &a) {
463 adouble tmp;
464 tmp.val=ADOLC_MATH_NSP::asin(a.val);
465 double tmp2=ADOLC_MATH_NSP::sqrt(1-a.val*a.val);
466 FOR_I_EQ_0_LT_NUMDIR
467 tmp.ADVAL_I=a.ADVAL_I/tmp2;
468 return tmp;
469 }
470
acos(const adouble & a)471 CUDADEV adouble acos(const adouble &a) {
472 adouble tmp;
473 tmp.val=ADOLC_MATH_NSP::acos(a.val);
474 double tmp2=-ADOLC_MATH_NSP::sqrt(1-a.val*a.val);
475 FOR_I_EQ_0_LT_NUMDIR
476 tmp.ADVAL_I=a.ADVAL_I/tmp2;
477 return tmp;
478 }
479
atan(const adouble & a)480 CUDADEV adouble atan(const adouble &a) {
481 adouble tmp;
482 tmp.val=ADOLC_MATH_NSP::atan(a.val);
483 double tmp2=1+a.val*a.val;
484 tmp2=1/tmp2;
485 if (tmp2!=0)
486 FOR_I_EQ_0_LT_NUMDIR
487 tmp.ADVAL_I=a.ADVAL_I*tmp2;
488 else
489 FOR_I_EQ_0_LT_NUMDIR
490 tmp.ADVAL_I=0.0;
491 return tmp;
492 }
493
atan2(const adouble & a,const adouble & b)494 CUDADEV adouble atan2(const adouble &a, const adouble &b) {
495 adouble tmp;
496 tmp.val=ADOLC_MATH_NSP::atan2(a.val, b.val);
497 double tmp2=a.val*a.val;
498 double tmp3=b.val*b.val;
499 double tmp4=tmp3/(tmp2+tmp3);
500 if (tmp4!=0)
501 FOR_I_EQ_0_LT_NUMDIR
502 tmp.ADVAL_I=(a.ADVAL_I*b.val-a.val*b.ADVAL_I)/tmp3*tmp4;
503 else
504 FOR_I_EQ_0_LT_NUMDIR
505 tmp.ADVAL_I=0.0;
506 return tmp;
507 }
508
pow(const adouble & a,double v)509 CUDADEV adouble pow(const adouble &a, double v) {
510 adouble tmp;
511 tmp.val=ADOLC_MATH_NSP::pow(a.val, v);
512 double tmp2=v*ADOLC_MATH_NSP::pow(a.val, v-1);
513 FOR_I_EQ_0_LT_NUMDIR
514 tmp.ADVAL_I=tmp2*a.ADVAL_I;
515 return tmp;
516 }
517
pow(const adouble & a,const adouble & b)518 CUDADEV adouble pow(const adouble &a, const adouble &b) {
519 adouble tmp;
520 tmp.val=ADOLC_MATH_NSP::pow(a.val, b.val);
521 double tmp2=b.val*ADOLC_MATH_NSP::pow(a.val, b.val-1);
522 double tmp3=ADOLC_MATH_NSP::log(a.val)*tmp.val;
523 FOR_I_EQ_0_LT_NUMDIR
524 tmp.ADVAL_I=tmp2*a.ADVAL_I+tmp3*b.ADVAL_I;
525 return tmp;
526 }
527
pow(double v,const adouble & a)528 CUDADEV adouble pow(double v, const adouble &a) {
529 adouble tmp;
530 tmp.val=ADOLC_MATH_NSP::pow(v, a.val);
531 double tmp2=tmp.val*ADOLC_MATH_NSP::log(v);
532 FOR_I_EQ_0_LT_NUMDIR
533 tmp.ADVAL_I=tmp2*a.ADVAL_I;
534 return tmp;
535 }
536
log10(const adouble & a)537 CUDADEV adouble log10(const adouble &a) {
538 adouble tmp;
539 tmp.val=ADOLC_MATH_NSP::log10(a.val);
540 double tmp2=ADOLC_MATH_NSP::log((double)10)*a.val;
541 FOR_I_EQ_0_LT_NUMDIR
542 tmp.ADVAL_I=a.ADVAL_I/tmp2;
543 return tmp;
544 }
545
sinh(const adouble & a)546 CUDADEV adouble sinh (const adouble &a) {
547 adouble tmp;
548 tmp.val=ADOLC_MATH_NSP::sinh(a.val);
549 double tmp2=ADOLC_MATH_NSP::cosh(a.val);
550 FOR_I_EQ_0_LT_NUMDIR
551 tmp.ADVAL_I=a.ADVAL_I*tmp2;
552 return tmp;
553 }
554
cosh(const adouble & a)555 CUDADEV adouble cosh (const adouble &a) {
556 adouble tmp;
557 tmp.val=ADOLC_MATH_NSP::cosh(a.val);
558 double tmp2=ADOLC_MATH_NSP::sinh(a.val);
559 FOR_I_EQ_0_LT_NUMDIR
560 tmp.ADVAL_I=a.ADVAL_I*tmp2;
561 return tmp;
562 }
563
tanh(const adouble & a)564 CUDADEV adouble tanh (const adouble &a) {
565 adouble tmp;
566 tmp.val=ADOLC_MATH_NSP::tanh(a.val);
567 double tmp2=ADOLC_MATH_NSP::cosh(a.val);
568 tmp2*=tmp2;
569 FOR_I_EQ_0_LT_NUMDIR
570 tmp.ADVAL_I=a.ADVAL_I/tmp2;
571 return tmp;
572 }
573
574 #if defined(ATRIG_ERF)
asinh(const adouble & a)575 CUDADEV adouble asinh (const adouble &a) {
576 adouble tmp;
577 tmp.val=ADOLC_MATH_NSP_ERF::asinh(a.val);
578 double tmp2=ADOLC_MATH_NSP::sqrt(a.val*a.val+1);
579 FOR_I_EQ_0_LT_NUMDIR
580 tmp.ADVAL_I=a.ADVAL_I/tmp2;
581 return tmp;
582 }
583
acosh(const adouble & a)584 CUDADEV adouble acosh (const adouble &a) {
585 adouble tmp;
586 tmp.val=ADOLC_MATH_NSP_ERF::acosh(a.val);
587 double tmp2=ADOLC_MATH_NSP::sqrt(a.val*a.val-1);
588 FOR_I_EQ_0_LT_NUMDIR
589 tmp.ADVAL_I=a.ADVAL_I/tmp2;
590 return tmp;
591 }
592
atanh(const adouble & a)593 CUDADEV adouble atanh (const adouble &a) {
594 adouble tmp;
595 tmp.val=ADOLC_MATH_NSP_ERF::atanh(a.val);
596 double tmp2=1-a.val*a.val;
597 FOR_I_EQ_0_LT_NUMDIR
598 tmp.ADVAL_I=a.ADVAL_I/tmp2;
599 return tmp;
600 }
601 #endif
602
fabs(const adouble & a)603 CUDADEV adouble fabs (const adouble &a) {
604 adouble tmp;
605 tmp.val=ADOLC_MATH_NSP::fabs(a.val);
606 int as=0;
607 if (a.val>0) as=1;
608 if (a.val<0) as=-1;
609 if (as!=0)
610 FOR_I_EQ_0_LT_NUMDIR
611 tmp.ADVAL_I=a.ADVAL_I*as;
612 else
613 FOR_I_EQ_0_LT_NUMDIR {
614 as=0;
615 if (a.ADVAL_I>0) as=1;
616 if (a.ADVAL_I<0) as=-1;
617 tmp.ADVAL_I=a.ADVAL_I*as;
618 }
619 return tmp;
620 }
621
ceil(const adouble & a)622 CUDADEV adouble ceil (const adouble &a) {
623 adouble tmp;
624 tmp.val=ADOLC_MATH_NSP::ceil(a.val);
625 FOR_I_EQ_0_LT_NUMDIR
626 tmp.ADVAL_I=0.0;
627 return tmp;
628 }
629
floor(const adouble & a)630 CUDADEV adouble floor (const adouble &a) {
631 adouble tmp;
632 tmp.val=ADOLC_MATH_NSP::floor(a.val);
633 FOR_I_EQ_0_LT_NUMDIR
634 tmp.ADVAL_I=0.0;
635 return tmp;
636 }
637
fmax(const adouble & a,const adouble & b)638 CUDADEV adouble fmax (const adouble &a, const adouble &b) {
639 adouble tmp;
640 double tmp2=a.val-b.val;
641 if (tmp2<0) {
642 tmp.val=b.val;
643 FOR_I_EQ_0_LT_NUMDIR
644 tmp.ADVAL_I=b.ADVAL_I;
645 } else {
646 tmp.val=a.val;
647 if (tmp2>0) {
648 FOR_I_EQ_0_LT_NUMDIR
649 tmp.ADVAL_I=a.ADVAL_I;
650 } else {
651 FOR_I_EQ_0_LT_NUMDIR
652 {
653 if (a.ADVAL_I<b.ADVAL_I) tmp.ADVAL_I=b.ADVAL_I;
654 else tmp.ADVAL_I=a.ADVAL_I;
655 }
656 }
657 }
658 return tmp;
659 }
660
fmax(double v,const adouble & a)661 CUDADEV adouble fmax (double v, const adouble &a) {
662 adouble tmp;
663 double tmp2=v-a.val;
664 if (tmp2<0) {
665 tmp.val=a.val;
666 FOR_I_EQ_0_LT_NUMDIR
667 tmp.ADVAL_I=a.ADVAL_I;
668 } else {
669 tmp.val=v;
670 if (tmp2>0) {
671 FOR_I_EQ_0_LT_NUMDIR
672 tmp.ADVAL_I=0.0;
673 } else {
674 FOR_I_EQ_0_LT_NUMDIR
675 {
676 if (a.ADVAL_I>0) tmp.ADVAL_I=a.ADVAL_I;
677 else tmp.ADVAL_I=0.0;
678 }
679 }
680 }
681 return tmp;
682 }
683
fmax(const adouble & a,double v)684 CUDADEV adouble fmax (const adouble &a, double v) {
685 adouble tmp;
686 double tmp2=a.val-v;
687 if (tmp2<0) {
688 tmp.val=v;
689 FOR_I_EQ_0_LT_NUMDIR
690 tmp.ADVAL_I=0.0;
691 } else {
692 tmp.val=a.val;
693 if (tmp2>0) {
694 FOR_I_EQ_0_LT_NUMDIR
695 tmp.ADVAL_I=a.ADVAL_I;
696 } else {
697 FOR_I_EQ_0_LT_NUMDIR
698 {
699 if (a.ADVAL_I>0) tmp.ADVAL_I=a.ADVAL_I;
700 else tmp.ADVAL_I=0.0;
701 }
702 }
703 }
704 return tmp;
705 }
706
fmin(const adouble & a,const adouble & b)707 CUDADEV adouble fmin (const adouble &a, const adouble &b) {
708 adouble tmp;
709 double tmp2=a.val-b.val;
710 if (tmp2<0) {
711 tmp.val=a.val;
712 FOR_I_EQ_0_LT_NUMDIR
713 tmp.ADVAL_I=a.ADVAL_I;
714 } else {
715 tmp.val=b.val;
716 if (tmp2>0) {
717 FOR_I_EQ_0_LT_NUMDIR
718 tmp.ADVAL_I=b.ADVAL_I;
719 } else {
720 FOR_I_EQ_0_LT_NUMDIR
721 {
722 if (a.ADVAL_I<b.ADVAL_I) tmp.ADVAL_I=a.ADVAL_I;
723 else tmp.ADVAL_I=b.ADVAL_I;
724 }
725 }
726 }
727 return tmp;
728 }
729
fmin(double v,const adouble & a)730 CUDADEV adouble fmin (double v, const adouble &a) {
731 adouble tmp;
732 double tmp2=v-a.val;
733 if (tmp2<0) {
734 tmp.val=v;
735 FOR_I_EQ_0_LT_NUMDIR
736 tmp.ADVAL_I=0.0;
737 } else {
738 tmp.val=a.val;
739 if (tmp2>0) {
740 FOR_I_EQ_0_LT_NUMDIR
741 tmp.ADVAL_I=a.ADVAL_I;
742 } else {
743 FOR_I_EQ_0_LT_NUMDIR
744 {
745 if (a.ADVAL_I<0) tmp.ADVAL_I=a.ADVAL_I;
746 else tmp.ADVAL_I=0.0;
747 }
748 }
749 }
750 return tmp;
751 }
752
fmin(const adouble & a,double v)753 CUDADEV adouble fmin (const adouble &a, double v) {
754 adouble tmp;
755 double tmp2=a.val-v;
756 if (tmp2<0) {
757 tmp.val=a.val;
758 FOR_I_EQ_0_LT_NUMDIR
759 tmp.ADVAL_I=a.ADVAL_I;
760 } else {
761 tmp.val=v;
762 if (tmp2>0) {
763 FOR_I_EQ_0_LT_NUMDIR
764 tmp.ADVAL_I=0.0;
765 } else {
766 FOR_I_EQ_0_LT_NUMDIR
767 {
768 if (a.ADVAL_I<0) tmp.ADVAL_I=a.ADVAL_I;
769 else tmp.ADVAL_I=0.0;
770 }
771 }
772 }
773 return tmp;
774 }
775
ldexp(const adouble & a,const adouble & b)776 CUDADEV adouble ldexp (const adouble &a, const adouble &b) {
777 return a*pow(2.,b);
778 }
779
ldexp(const adouble & a,const double v)780 CUDADEV adouble ldexp (const adouble &a, const double v) {
781 return a*ADOLC_MATH_NSP::pow(2.,v);
782 }
783
ldexp(const double v,const adouble & a)784 CUDADEV adouble ldexp (const double v, const adouble &a) {
785 return v*pow(2.,a);
786 }
787
frexp(const adouble & a,int * v)788 CUDADEV double frexp (const adouble &a, int* v) {
789 return ADOLC_MATH_NSP::frexp(a.val, v);
790 }
791
792 #if defined(ATRIG_ERF)
erf(const adouble & a)793 CUDADEV adouble erf (const adouble &a) {
794 adouble tmp;
795 tmp.val=ADOLC_MATH_NSP_ERF::erf(a.val);
796 double tmp2 = 2.0 /
797 ADOLC_MATH_NSP_ERF::sqrt(ADOLC_MATH_NSP::acos(-1.0)) *
798 ADOLC_MATH_NSP_ERF::exp(-a.val*a.val);
799 FOR_I_EQ_0_LT_NUMDIR
800 tmp.ADVAL_I=tmp2*a.ADVAL_I;
801 return tmp;
802 }
803 #endif
804
805
806 /******************* nontemporary results *********************************/
807 CUDADEV void adouble::operator = (const double v) {
808 val=v;
809 FOR_I_EQ_0_LT_NUMDIR
810 ADVAL_I=0.0;
811 }
812
813 CUDADEV void adouble::operator = (const adouble& a) {
814 val=a.val;
815 FOR_I_EQ_0_LT_NUMDIR
816 ADVAL_I=a.ADVAL_I;
817 }
818
819 CUDADEV void adouble::operator += (const double v) {
820 val+=v;
821 }
822
823 CUDADEV void adouble::operator += (const adouble& a) {
824 val=val+a.val;
825 FOR_I_EQ_0_LT_NUMDIR
826 ADVAL_I+=a.ADVAL_I;
827 }
828
829 CUDADEV void adouble::operator -= (const double v) {
830 val-=v;
831 }
832
833 CUDADEV void adouble::operator -= (const adouble& a) {
834 val=val-a.val;
835 FOR_I_EQ_0_LT_NUMDIR
836 ADVAL_I-=a.ADVAL_I;
837 }
838
839 CUDADEV void adouble::operator *= (const double v) {
840 val=val*v;
841 FOR_I_EQ_0_LT_NUMDIR
842 ADVAL_I*=v;
843 }
844
845 CUDADEV void adouble::operator *= (const adouble& a) {
846 FOR_I_EQ_0_LT_NUMDIR
847 ADVAL_I=ADVAL_I*a.val+val*a.ADVAL_I;
848 val*=a.val;
849 }
850
851 CUDADEV void adouble::operator /= (const double v) {
852 val/=v;
853 FOR_I_EQ_0_LT_NUMDIR
854 ADVAL_I/=v;
855 }
856
857 CUDADEV void adouble::operator /= (const adouble& a) {
858 FOR_I_EQ_0_LT_NUMDIR
859 ADVAL_I=(ADVAL_I*a.val-val*a.ADVAL_I)/(a.val*a.val);
860 val=val/a.val;
861 }
862
863 // not
864 CUDADEV int adouble::operator ! () const {
865 return val==0.0;
866 }
867
868 // comparision
869 CUDADEV int adouble::operator != (const adouble &a) const {
870 return val!=a.val;
871 }
872
873 CUDADEV int adouble::operator != (const double v) const {
874 return val!=v;
875 }
876
877 CUDADEV int operator != (const double v, const adouble &a) {
878 return v!=a.val;
879 }
880
881 CUDADEV int adouble::operator == (const adouble &a) const {
882 return val==a.val;
883 }
884
885 CUDADEV int adouble::operator == (const double v) const {
886 return val==v;
887 }
888
889 CUDADEV int operator == (const double v, const adouble &a) {
890 return v==a.val;
891 }
892
893 CUDADEV int adouble::operator <= (const adouble &a) const {
894 return val<=a.val;
895 }
896
897 CUDADEV int adouble::operator <= (const double v) const {
898 return val<=v;
899 }
900
901 CUDADEV int operator <= (const double v, const adouble &a) {
902 return v<=a.val;
903 }
904
905 CUDADEV int adouble::operator >= (const adouble &a) const {
906 return val>=a.val;
907 }
908
909 CUDADEV int adouble::operator >= (const double v) const {
910 return val>=v;
911 }
912
913 CUDADEV int operator >= (const double v, const adouble &a) {
914 return v>=a.val;
915 }
916
917 CUDADEV int adouble::operator > (const adouble &a) const {
918 return val>a.val;
919 }
920
921 CUDADEV int adouble::operator > (const double v) const {
922 return val>v;
923 }
924
925 CUDADEV int operator > (const double v, const adouble &a) {
926 return v>a.val;
927 }
928
929 CUDADEV int adouble::operator < (const adouble &a) const {
930 return val<a.val;
931 }
932
933 CUDADEV int adouble::operator < (const double v) const {
934 return val<v;
935 }
936
937 CUDADEV int operator < (const double v, const adouble &a) {
938 return v<a.val;
939 }
940
941 /******************* getter / setter **************************************/
getValue()942 CUDAHOSTDEV double adouble::getValue() const {
943 return val;
944 }
945
setValue(const double v)946 CUDAHOSTDEV void adouble::setValue(const double v) {
947 val=v;
948 }
949
getADValue()950 CUDAHOSTDEV ADVAL_TYPE adouble::getADValue() const {
951 return (ADVAL_TYPE)adval;
952 }
953
setADValue(ADVAL_TYPE v)954 CUDAHOSTDEV void adouble::setADValue(ADVAL_TYPE v) {
955 FOR_I_EQ_0_LT_NUMDIR
956 ADVAL_I=V_I;
957 }
958
959 # if defined(NUMBER_DIRECTIONS)
getADValue(const unsigned int p)960 CUDAHOSTDEV double adouble::getADValue(const unsigned int p) const {
961 unsigned int locp = p;
962 if (locp>=ADOLC_numDir) {
963 locp = ADOLC_numDir -1;
964 }
965 return adval[locp];
966 }
967
setADValue(const unsigned int p,const double v)968 CUDAHOSTDEV void adouble::setADValue(const unsigned int p, const double v) {
969 unsigned int locp = p;
970 if (locp>=ADOLC_numDir) {
971 locp = ADOLC_numDir - 1;
972 }
973 adval[locp]=v;
974 }
975 # endif
976
977 #if defined(NUMBER_DIRECTIONS)
setNumDir(const unsigned int p)978 void setNumDir(const unsigned int p) {
979 #if !defined(DYNAMIC_DIRECTIONS)
980 if (p>NUMBER_DIRECTIONS) ADOLC_numDir=NUMBER_DIRECTIONS;
981 else ADOLC_numDir=p;
982 #else
983 ADOLC_numDir = p;
984 #endif
985 }
986 #endif
987
988 /******************* i/o operations ***************************************/
989 CUDAHOST ostream& operator << ( ostream& out, const adouble& a) {
990 out << "Value: " << a.val;
991 #if !defined(NUMBER_DIRECTIONS)
992 out << " ADValue: ";
993 #else
994 out << " ADValues (" << ADOLC_numDir << "): ";
995 #endif
996 FOR_I_EQ_0_LT_NUMDIR
997 out << a.ADVAL_I << " ";
998 out << "(a)";
999 return out;
1000 }
1001
1002 CUDAHOST istream& operator >> ( istream& in, adouble& a) {
1003 char c;
1004 do {
1005 in >> c;
1006 } while (c!=':' && !in.eof());
1007 in >> a.val;
1008 #if !defined(NUMBER_DIRECTIONS)
1009 do in >> c;
1010 while (c!=':' && !in.eof());
1011 #else
1012 unsigned int num;
1013 do in >> c;
1014 while (c!='(' && !in.eof());
1015 in >> num;
1016 if (num>NUMBER_DIRECTIONS) {
1017 cout << "ADOL-C error: to many directions in input\n";
1018 exit(-1);
1019 }
1020 do in >> c;
1021 while (c!=')' && !in.eof());
1022 #endif
1023 FOR_I_EQ_0_LT_NUMDIR
1024 in >> a.ADVAL_I;
1025 do in >> c;
1026 while (c!=')' && !in.eof());
1027 return in;
1028 }
1029 }
1030 /****************************************************************************/
1031 /* end traceless gpu implementation first order derivatives */
1032 /****************************************************************************/
1033 #endif
1034