source: trunk/ADOL-C/examples/additional_examples/ext_diff_func/ext_diff_func.cpp @ 40

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

adapted directoy structure

File size: 3.8 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     ext_diff_func.cpp
4 Revision: $Id: ext_diff_func.cpp 37 2009-05-28 12:56:44Z awalther $
5 Contents: example for external differentiated functions
6
7 Copyright (c) Andrea Walther
8 
9 This file is part of ADOL-C. This software is provided as open source.
10 Any use, reproduction, or distribution of the software constitutes
11 recipient's acceptance of the terms of the accompanying license file.
12 
13---------------------------------------------------------------------------*/
14#include <math.h>
15#include <adolc.h>
16
17#define h 0.01
18#define steps 100
19
20// time step function
21// double version
22int euler_step(int n, double *yin, int m, double *yout);
23
24// adouble version
25void euler_step_act(int n, adouble *yin, int m, adouble *yout);
26
27
28// versions for usage as external differentiated function
29int zos_for_euler_step(int n, double *yin, int m, double *yout);
30int fos_rev_euler_step(int n, double *u, int m, double *z);
31
32ext_diff_fct *edf;
33int tag_full, tag_part, tag_ext_fct;
34
35
36int main()
37{
38  // time interval
39  double t0, tf;
40
41  // state, double and adouble version
42  adouble y[2];
43  adouble ynew[2];
44  int n, m;
45
46  // control, double and adouble version
47  adouble con[2];
48  double conp[2];
49
50  // target value;
51  double f;
52
53  //variables for derivative caluclation
54  double yp[2], ynewp[2];
55  double u[2], z[2];
56  double grad[2];
57
58
59  int i,j;
60 
61  // tape identifiers
62  tag_full = 1;
63  tag_part = 2;
64  tag_ext_fct = 3;
65
66  // two input variables for external differentiated function
67  n = 2;
68  // two output variables for external differentiated function
69  m = 2;
70
71  // time interval
72  t0 = 0.0;
73  tf = 1.0;
74
75  //control
76  conp[0] = 1.0;
77  conp[1] = 1.0;
78
79  trace_on(tag_full);
80    con[0] <<= conp[0];
81    con[1] <<= conp[1];
82    y[0] = con[0];
83    y[1] = con[1];
84 
85    for(i=0;i<steps;i++)
86      {
87        euler_step_act(n,y,m,ynew);
88        for(j=0;j<2;j++)
89          y[j] = ynew[j];
90      }
91    y[0] + y[1] >>= f;
92  trace_off(1);
93
94  gradient(tag_full,2,conp,grad);
95 
96  printf(" full taping:\n gradient=( %f, %f)\n\n",grad[0],grad[1]);
97
98  // Now using external function facilities
99
100  // tape external differentiated function
101
102  trace_on(tag_ext_fct);
103    y[0] <<= conp[0];
104    y[1] <<= conp[1];
105 
106    euler_step_act(2,y,2,ynew);
107    ynew[0] >>= f;
108    ynew[1] >>= f;
109  trace_off(1);
110
111
112  // register external function
113  edf = reg_ext_fct(euler_step);
114
115  // information for Zero-Order-Scalar (=zos) forward
116//   yp = new double[2];
117//   ynewp = new double[2];
118  edf->zos_forward = zos_for_euler_step;
119  edf->dp_x = yp;
120  edf->dp_y = ynewp;
121  // information for First-Order-Scalar (=fos) reverse
122  edf->fos_reverse = fos_rev_euler_step;
123  edf->dp_U = u;
124  edf->dp_Z = z;
125
126  trace_on(tag_part);
127    con[0] <<= conp[0];
128    con[1] <<= conp[1];
129    y[0] = con[0];
130    y[1] = con[1];
131 
132    for(i=0;i<steps;i++)
133      {
134        call_ext_fct(edf, 2, yp, y, 2, ynewp, ynew);
135        for(j=0;j<2;j++)
136          y[j] = ynew[j];
137      }
138    y[0] + y[1] >>= f;
139  trace_off(1);
140  gradient(tag_part,2,conp,grad);
141 
142  printf(" taping with external function facility:\n gradient=( %f, %f)\n\n",grad[0],grad[1]);
143
144  return 0;
145}
146
147void euler_step_act(int n, adouble *yin, int m, adouble *yout)
148{
149
150 // Euler step, adouble version
151 
152 yout[0] = yin[0]+h*yin[0];
153 yout[1] = yin[1]+h*2*yin[1];
154}
155
156int euler_step(int n, double *yin, int m, double *yout)
157{
158
159 // Euler step, double version
160 yout[0] = yin[0]+h*yin[0];
161 yout[1] = yin[1]+h*2*yin[1];
162
163 return 1;
164}
165
166int zos_for_euler_step(int n, double *yin, int m, double *yout)
167{
168  int rc;
169
170  rc = zos_forward(tag_ext_fct, 2, 2, 0, yin, yout);
171
172  return rc;
173}
174
175int fos_rev_euler_step(int n, double *u, int m, double *z)
176{
177  int rc;
178
179  zos_forward(tag_ext_fct, 2, 2, 1, edf->dp_x, edf->dp_y);
180  rc = fos_reverse(tag_ext_fct, 2, 2, u, z);
181
182  return rc;
183}
Note: See TracBrowser for help on using the repository browser.