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