source: trunk/ADOL-C/examples/speelpenning.cpp @ 103

Last change on this file since 103 was 103, checked in by awalther, 9 years ago

correction in taylor.c, 1 -> tag

  • Property svn:keywords set to Author Date Id Revision
File size: 2.9 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     speelpenning.cpp
4 Revision: $Id: speelpenning.cpp 103 2010-05-12 09:31:41Z awalther $
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 <adouble.h>            // use of active doubles
19#include <drivers/drivers.h>    // use of "Easy to Use" drivers
20// gradient(.) and hessian(.)
21#include <taping.h>             // use of taping
22
23#include <iostream>
24using namespace std;
25
26#include <cstdlib>
27#include <math.h>
28
29/****************************************************************************/
30/*                                                             MAIN PROGRAM */
31int 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
Note: See TracBrowser for help on using the repository browser.