source: trunk/ADOL-C/examples/additional_examples/scal/scalexam.cpp @ 171

Last change on this file since 171 was 171, checked in by kulshres, 9 years ago

Squashed merge branch 'master' of 'gitclone' into svn

  • 'master' of 'gitclone': (84 commits) adjust example makefiles and include paths get rid of the symlink in the src subdirectory

details of the commits:
commit c9e4bc332d2363f737fc2e8a8fcfc2e43ddb9d15
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Mon Oct 4 15:43:47 2010 +0200

adjust example makefiles and include paths

include paths in example sources were wrong for some time now
simplify makefile rules too, there is really no need for checking SPARSE
adjust include paths in makefiles.

Signed-off-by: Kshitij Kulshreshtha <kshitij@…>

commit e6e1963e41e097fd5b4a79cd1611c12f6868dc94
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Mon Oct 4 15:41:25 2010 +0200

get rid of the symlink in the src subdirectory

windows doesn't like symlinks and make infinite depth directories
we now create a symlink for build in the directory parallel to src
adjust all makefiles.am accordingly for build

Signed-off-by: Kshitij Kulshreshtha <kshitij@…>

  • Property svn:keywords set to Author Date Id Revision
File size: 7.2 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     scalexam.cpp
4 Revision: $Id: scalexam.cpp 171 2010-10-04 13:57:19Z kulshres $
5 Contents:
6          This program can be used to verify the consistency and
7          correctness of derivatives computed by ADOL-C in its forward
8          and reverse mode. 
9          Ther use is required to select one integer input id.
10          For positive n = id the monomial x^n is evaluated recursively
11          at x=0.5 and all its nonzero Taylor coeffcients at this point
12          are evaluated in the forward and reverse mode.
13          A negative choice of id >= -9 leads to one of nine
14          identities, whose derivatives should be trivial. These identities
15          may be used to check the correctness of particular code segments
16          in the ADOL-C sources uni5_.c and *o_rev.c. No timings are
17          performed in this example program.
18
19 Copyright (c) Andrea Walther, Andreas Griewank, Andreas Kowarz,
20               Hristo Mitev, Sebastian Schlenkrich, Jean Utke, Olaf Vogel
21 
22 This file is part of ADOL-C. This software is provided as open source.
23 Any use, reproduction, or distribution of the software constitutes
24 recipient's acceptance of the terms of the accompanying license file.
25 
26---------------------------------------------------------------------------*/
27
28/****************************************************************************/
29/*                                                                 INCLUDES */
30#include <adolc/adolc.h>
31
32#include <math.h>
33#include <iostream>
34using namespace std;
35
36/****************************************************************************/
37/*                                                                    POWER */
38/* The monomial evaluation routine which has been obtained from
39   the original version by retyping all `doubles' as `adoubles' */
40adouble power( adouble x, int n ) {
41    adouble z = 1;
42    if (n > 0) {
43        int nh =n/2;
44        z = power(x,nh);
45        z *= z;
46        if (2*nh != n)
47            z *= x;
48        return z;
49    } else
50        if (n == 0)
51            return z;
52        else
53            return 1.0/power(x,-n);
54}
55
56
57/****************************************************************************/
58/*                                                                     MAIN */
59int main() {
60    int n, i, id;
61    int tag = 0;
62    /*--------------------------------------------------------------------------*/
63    fprintf(stdout,"SCALEXAM (ADOL-C Example)\n\n");
64    fprintf(stdout,"problem number(-1 .. -10) / degree of monomial =? \n");
65    scanf("%d",&id);
66    n = id >0 ? id : 3;
67
68    double *xp,*yp;
69    xp = new double[n+4];
70    yp = new double[n+4];
71    yp[0] = 0;
72    xp[0] = 0.5;
73    xp[1] = 1.0;
74
75    /*--------------------------------------------------------------------------*/
76    int dum = 1;
77    trace_on(tag,dum);   // Begin taping all calculations with 'adoubles'
78    adouble y,x;
79    x <<= xp[0];
80    if (id >= 0) {
81        fprintf(stdout,"Evaluate and differentiate recursive power routine \n");
82        y = power(x,n);
83    } else {
84        fprintf(stdout,
85                "Check Operations and Functions by Algebraic Identities \n");
86        switch (id) {
87            case -1 :
88                fprintf(stdout,
89                        "Addition/Subtraction: y = x + x - (2.0/3)*x - x/3 \n");
90                y =  x + x - (2.0/3)*x - x/3 ;
91                break;
92            case -2 :
93                fprintf(stdout,"Multiplication/divison:  y = x*x/x \n");
94                y = x*x/x;
95                break;
96            case -3 :
97                fprintf(stdout,"Square root and power: y = sqrt(pow(x,2)) \n");
98                y = sqrt(pow(x,2));
99                break;
100            case -4 :
101                fprintf(stdout,"Exponential and log: y = exp(log(log(exp(x)))) \n");
102                y = exp(log(log(exp(x))));
103                break;
104            case -5 :
105                fprintf(stdout,"Trig identity: y = x + sin(2*x)-2*cos(x)*sin(x) \n");
106                y =  x + sin(2.0*x)-2.0*cos(x)*sin(x);
107                break;
108            case -6 :
109                fprintf(stdout,"Check out quadrature macro \n");
110                y = exp(myquad(myquad(exp(x))));
111                break;
112            case -7 :
113                fprintf(stdout,"Arcsin: y = sin(asin(acos(cos(x)))) \n");
114                y = sin(asin(acos(cos(x))));
115                break;
116            case -8 :
117                fprintf(stdout,
118                        "Hyperbolic tangent: y = x + tanh(x)-sinh(x)/cosh(x) \n");
119                y = x + tanh(x)-sinh(x)/cosh(x) ;
120                break;
121            case -9 :
122                fprintf(stdout,"Absolute value: y = x + fabs(x) - fabs(-x) \n");
123                y = x + fabs(-x) - fabs(x);
124                break;
125            case -10 :
126                fprintf(stdout,"atan2: y = atan2(sin(x-0.5+pi),cos(x-0.5+pi)) \n");
127                y = atan2(sin(x),cos(x));
128                break;
129            default :
130                fprintf(stdout," Please select problem number >= -10 \n");
131                exit(-1);
132        }
133    }
134    y >>= yp[0];
135    trace_off();  // The (partial) execution trace is completed.
136
137    /*--------------------------------------------------------------------------*/
138    if( id < 0 )
139        fprintf(stdout,"Round-off error: %14.6E\n",(y-x).value());
140
141    /*--------------------------------------------------------------------------*/
142    int tape_stats[STAT_SIZE];
143    tapestats(tag,tape_stats);
144
145    fprintf(stdout,"\n    independents            %d\n",tape_stats[NUM_INDEPENDENTS]);
146    fprintf(stdout,"    dependents              %d\n",tape_stats[NUM_DEPENDENTS]);
147    fprintf(stdout,"    operations              %d\n",tape_stats[NUM_OPERATIONS]);
148    fprintf(stdout,"    operations buffer size  %d\n",tape_stats[OP_BUFFER_SIZE]);
149    fprintf(stdout,"    locations buffer size   %d\n",tape_stats[LOC_BUFFER_SIZE]);
150    fprintf(stdout,"    constants buffer size   %d\n",tape_stats[VAL_BUFFER_SIZE]);
151    fprintf(stdout,"    maxlive                 %d\n",tape_stats[NUM_MAX_LIVES]);
152    fprintf(stdout,"    valstack size           %d\n\n",tape_stats[TAY_STACK_SIZE]);
153
154    /*--------------------------------------------------------------------------*/
155    double *res;
156    res = new double[n+2];
157    double u[1];
158    u[0] = 1;
159    fprintf(stdout,
160            "\nThe two Taylor coefficients in each row should agree\n\n");
161
162    double ***V = (double***)new double**[1];
163    V[0] = new double*[1];
164    V[0][0] = new double[n+2];
165    double **U = new double*[1];
166    U[0] = new double[1];
167    U[0][0] = 1;
168    double** xpoint = &xp;
169    double** ypoint = &yp;
170    double** respoint = &res;
171
172    // tape_doc(tag,depen,indep,*xpoint,*respoint);
173
174    fprintf(stdout," \n \t   forward  \t    reverse  \n");
175    for (i=0; i < n+2; i++) {
176        xp[i+2]=0;
177        forward(tag,1,1,i,i+1,xpoint,respoint);
178        fprintf(stdout,"%d\t%14.6E\t\t%14.6E\n",i,res[i],yp[i]);
179        reverse(tag,1,1,i,u,ypoint); // call higher order scalar reverse
180        reverse(tag,1,1,i,1,U,V);
181        yp[i+1] = yp[i]/(i+1);
182        if (V[0][0][i] != yp[i])
183            fprintf(stdout,"%d-th component in error %14.6E\n",i,V[0][0][i]-yp[i]);
184    }
185    cout << "\nWhen n<0 all rows except the first two should vanish \n";
186
187    return 1;
188}
Note: See TracBrowser for help on using the repository browser.