source: trunk/ADOL-C/include/adolc/adoublecuda.h @ 408

Last change on this file since 408 was 408, checked in by kulshres, 7 years ago

update revision id stamps

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