1 | /*---------------------------------------------------------------------------- |
---|
2 | ADOL-C -- Automatic Differentiation by Overloading in C++ |
---|
3 | File: speelpenning.cpp |
---|
4 | Revision: $Id: speelpenning.cpp 171 2010-10-04 13:57:19Z kulshres $ |
---|
5 | Contents: speelpennings example, 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 */ |
---|
18 | #include <adolc/adouble.h> // use of active doubles |
---|
19 | #include <adolc/drivers/drivers.h> // use of "Easy to Use" drivers |
---|
20 | // gradient(.) and hessian(.) |
---|
21 | #include <adolc/taping.h> // use of taping |
---|
22 | |
---|
23 | #include <iostream> |
---|
24 | using namespace std; |
---|
25 | |
---|
26 | #include <cstdlib> |
---|
27 | #include <math.h> |
---|
28 | |
---|
29 | /****************************************************************************/ |
---|
30 | /* MAIN PROGRAM */ |
---|
31 | int main() { |
---|
32 | int n,i,j; |
---|
33 | int tape_stats[STAT_SIZE]; |
---|
34 | |
---|
35 | cout << "SPEELPENNINGS PRODUCT (ADOL-C Documented Example)\n\n"; |
---|
36 | cout << "number of independent variables = ? \n"; |
---|
37 | cin >> n; |
---|
38 | |
---|
39 | double *xp = new double[n]; |
---|
40 | double yp = 0.0; |
---|
41 | adouble *x = new adouble[n]; |
---|
42 | adouble y = 1; |
---|
43 | |
---|
44 | for(i=0; i<n; i++) |
---|
45 | xp[i] = (i+1.0)/(2.0+i); // some initialization |
---|
46 | |
---|
47 | trace_on(1); // tag = 1, keep = 0 by default |
---|
48 | for(i=0; i<n; i++) { |
---|
49 | x[i] <<= xp[i]; // or x <<= xp outside the loop |
---|
50 | y *= x[i]; |
---|
51 | } // end for |
---|
52 | y >>= yp; |
---|
53 | delete[] x; |
---|
54 | trace_off(1); |
---|
55 | |
---|
56 | tapestats(1,tape_stats); // reading of tape statistics |
---|
57 | cout<<"maxlive "<<tape_stats[NUM_MAX_LIVES]<<"\n"; |
---|
58 | // ..... print other tape stats |
---|
59 | |
---|
60 | double* g = new double[n]; |
---|
61 | gradient(1,n,xp,g); // gradient evaluation |
---|
62 | |
---|
63 | double** H = (double**) malloc(n*sizeof(double*)); |
---|
64 | for(i=0; i<n; i++) |
---|
65 | H[i] = (double*)malloc((i+1)*sizeof(double)); |
---|
66 | hessian(1,n,xp,H); // H equals (n-1)g since g is |
---|
67 | double errg = 0; // homogeneous of degree n-1. |
---|
68 | double errh = 0; |
---|
69 | for(i=0; i<n; i++) |
---|
70 | errg += fabs(g[i]-yp/xp[i]); // vanishes analytically. |
---|
71 | for(i=0; i<n; i++) { |
---|
72 | for(j=0; j<n; j++) { |
---|
73 | if (i>j) // lower half of hessian |
---|
74 | errh += fabs(H[i][j]-g[i]/xp[j]); |
---|
75 | } // end for |
---|
76 | } // end for |
---|
77 | cout << yp-1/(1.0+n) << " error in function \n"; |
---|
78 | cout << errg <<" error in gradient \n"; |
---|
79 | cout << errh <<" consistency check \n"; |
---|
80 | |
---|
81 | return 0; |
---|
82 | } // end main |
---|
83 | |
---|