source: trunk/ADOL-C/examples/odexam.cpp @ 160

Last change on this file since 160 was 42, checked in by awalther, 10 years ago

set svn keywords property

  • Property svn:keywords set to Author Date Id Revision
File size: 3.5 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     odexam.cpp
4 Revision: $Id: odexam.cpp 42 2009-07-15 18:37:17Z kulshres $
5 Contents: Nonlinear ordinary differential equation based on the
6           Robertson test problem, described in the manual
7
8 Copyright (c) Andrea Walther, Andreas Griewank, Andreas Kowarz,
9               Hristo Mitev, Sebastian Schlenkrich, Jean Utke, Olaf Vogel
10 
11 This file is part of ADOL-C. This software is provided as open source.
12 Any use, reproduction, or distribution of the software constitutes
13 recipient's acceptance of the terms of the accompanying license file.
14 
15---------------------------------------------------------------------------*/
16
17/****************************************************************************/
18/*                                                                 INCLUDES */
19#include <adouble.h>            // use of active double
20#include <drivers/odedrivers.h> // use of "Easy To Use" ODE drivers
21#include <adalloc.h>            // use of ADOL-C allocation utilities
22#include <taping.h>             // use of taping
23
24// NOTICE: If one wants to include all ADOL-C interfaces without
25//         getting in trouble to find out the right header files
26//
27//         #include "adolc.h"
28//
29//         will do the right work.
30
31#include <iostream>
32using namespace std;
33
34/****************************************************************************/
35/*                                                          ADOUBLE ROUTINE */
36void tracerhs(short int tag, double* py, double* pyprime) {
37    adouble y[3];                // we left the parameters passive
38    adouble yprime[3];           
39
40    trace_on(tag);
41    y[0] <<= py[0]; y[1] <<= py[1];   y[2] <<= py[2];                              // initialize and mark independents
42    yprime[0] = -sin(y[2]) + 1.0e8*y[2]*(1.0-1.0/y[0]);
43    yprime[1] = -10.0*y[0] + 3.0e7*y[2]*(1-y[1]);
44    yprime[2] = -yprime[0] - yprime[1];
45    yprime[0] >>= pyprime[0]; yprime[1] >>= pyprime[1];  yprime[2] >>= pyprime[2]; // mark and pass dependents
46    trace_off(1);                // write tape array onto a file
47} // end tracerhs
48
49/****************************************************************************/
50/*                                                             MAIN PROGRAM */
51int main() {
52    int i,j,deg;
53    int n = 3;
54    double py[3];
55    double pyp[3];
56
57    cout << "MODIFIED ROBERTSON TEST PROBLEM (ADOL-C Documented Example)\n\n";
58    cout << "degree of Taylor series =?\n";
59    cin >> deg;
60
61    short** nz = new short*[n];
62    double  **X;
63    double ***Z;
64    double ***B;
65
66    X = myalloc2(n,deg+1);
67    Z = myalloc3(n,n,deg);
68    B = myalloc3(n,n,deg);
69
70    for(i=0; i<n; i++) {
71        py[i]   = (i==0) ? 1.0 : 0.0;       // Initialize the base point
72        X[i][0] = py[i];                    // and the Taylor coefficient;
73        nz[i]   = new short[n];             // set up sparsity array
74    } // end for
75
76    tracerhs(1,py,pyp);                   // trace RHS with tag = 1
77
78    forode(1,n,deg,X);                    // compute deg coefficients
79    reverse(1,n,n,deg-1,Z,nz);            // U defaults to the identity
80    accode(n,deg-1,Z,B,nz);
81
82    cout << "nonzero pattern:\n";
83    for(i=0; i<n; i++) {
84        for(j=0; j<n; j++)
85            cout << nz[i][j]<<"\t";
86        cout <<"\n";
87    } // end for
88    cout << "\n";
89    cout << " 4 = transcend , 3 = rational , 2 = polynomial ,"
90    << " 1 = linear , 0 = zero \n";
91    cout << " negative number k indicate that entries of all"
92    << " B_j with j < -k vanish \n";
93
94    return 1;
95} // end main
Note: See TracBrowser for help on using the repository browser.