[40] | 1 | /*---------------------------------------------------------------------------- |
---|
| 2 | ADOL-C -- Automatic Differentiation by Overloading in C++ |
---|
| 3 | File: powexam.cpp |
---|
[42] | 4 | Revision: $Id: powexam.cpp 171 2010-10-04 13:57:19Z kulshres $ |
---|
[40] | 5 | Contents: computation of n-th power, described in the manual |
---|
| 6 | |
---|
| 7 | Copyright (c) Andrea Walther, Andreas Griewank, Andreas Kowarz, |
---|
| 8 | Hristo Mitev, Sebastian Schlenkrich, Jean Utke, Olaf Vogel |
---|
| 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 | /****************************************************************************/ |
---|
| 17 | /* INCLUDES */ |
---|
[171] | 18 | #include <adolc/adolc.h> /* use of ALL ADOL-C interfaces */ |
---|
[40] | 19 | |
---|
| 20 | #include <iostream> |
---|
| 21 | using namespace std; |
---|
| 22 | |
---|
| 23 | /****************************************************************************/ |
---|
| 24 | /* ADOUBLE ROUTINE */ |
---|
| 25 | adouble power(adouble x, int n) { |
---|
| 26 | adouble z = 1; |
---|
| 27 | |
---|
| 28 | if (n>0) /* Recursion and branches */ |
---|
| 29 | { int nh = n/2; /* that do not depend on */ |
---|
| 30 | z = power(x,nh); /* adoubles are fine !!!! */ |
---|
| 31 | z *= z; |
---|
| 32 | if (2*nh != n) |
---|
| 33 | z *= x; |
---|
| 34 | return z; |
---|
| 35 | } /* end if */ |
---|
| 36 | else { |
---|
| 37 | if (n==0) /* The local adouble z dies */ |
---|
| 38 | return z; /* as it goes out of scope. */ |
---|
| 39 | else |
---|
| 40 | return 1/power(x,-n); |
---|
| 41 | } /* end else */ |
---|
| 42 | } /* end power */ |
---|
| 43 | |
---|
| 44 | /****************************************************************************/ |
---|
| 45 | /* MAIN PROGRAM */ |
---|
| 46 | int main() { |
---|
| 47 | int i,tag = 1; |
---|
| 48 | int n; |
---|
| 49 | |
---|
| 50 | cout << "COMPUTATION OF N-TH POWER (ADOL-C Documented Example)\n\n"; |
---|
| 51 | cout << "monomial degree=? \n"; /* input the desired degree */ |
---|
| 52 | cin >> n; |
---|
| 53 | /* allocations and initializations */ |
---|
| 54 | double** X; |
---|
| 55 | double** Y; |
---|
| 56 | X = myalloc2(1,n+4); |
---|
| 57 | Y = myalloc2(1,n+4); |
---|
| 58 | X[0][0] = 0.5; /* function value = 0. coefficient */ |
---|
| 59 | X[0][1] = 1.0; /* first derivative = 1. coefficient */ |
---|
| 60 | for(i=0; i<n+2; i++) |
---|
| 61 | X[0][i+2] = 0; /* further coefficients */ |
---|
| 62 | double** Z; /* used for checking consistency */ |
---|
| 63 | Z = myalloc2(1,n+2); /* between forward and reverse */ |
---|
| 64 | |
---|
| 65 | adouble y,x; /* declare active variables */ |
---|
| 66 | /* beginning of active section */ |
---|
| 67 | trace_on(tag); /* tag = 1 and keep = 0 */ |
---|
| 68 | x <<= X[0][0]; /* only one independent var */ |
---|
| 69 | y = power(x,n); /* actual function call */ |
---|
| 70 | y >>= Y[0][0]; /* only one dependent adouble */ |
---|
| 71 | trace_off(); /* no global adouble has died */ |
---|
| 72 | /* end of active section */ |
---|
| 73 | double u[1]; /* weighting vector */ |
---|
| 74 | u[0]=1; /* for reverse call */ |
---|
| 75 | for(i=0; i<n+2; i++) /* note that keep = i+1 in call */ |
---|
| 76 | { forward(tag,1,1,i,i+1,X,Y); /* evaluate the i-the derivative */ |
---|
| 77 | if (i==0) |
---|
| 78 | cout << Y[0][i] << " - " << y.value() << " = " << Y[0][i]-y.value() |
---|
| 79 | << " (should be 0)\n"; |
---|
| 80 | else { |
---|
| 81 | Z[0][i] = Z[0][i-1]/i; /* scale derivative to Taylorcoeff. */ |
---|
| 82 | cout << Y[0][i] << " - " << Z[0][i] << " = " << Y[0][i]-Z[0][i] |
---|
| 83 | << " (should be 0)\n"; |
---|
| 84 | } |
---|
| 85 | reverse(tag,1,1,i,u,Z); /* evaluate the (i+1)-st deriv. */ |
---|
| 86 | } /* end for */ |
---|
| 87 | |
---|
| 88 | return 1; |
---|
| 89 | } /* end main */ |
---|
| 90 | |
---|