source: trunk/ADOL-C/examples/additional_examples/lighthouse/cubic.cpp @ 42

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

set svn keywords property

  • Property svn:keywords set to Author Date Id Revision
File size: 5.1 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     cubic.cpp
4 Revision: $Id: cubic.cpp 42 2009-07-15 18:37:17Z awalther $
5 Contents: example for cubic lighthouse example of Griewank's Book
6           using Cardan's formula with two tapes
7 
8
9 Copyright (c) Andrea Walther, Andreas Griewank, Andreas Kowarz,
10               Hristo Mitev, Sebastian Schlenkrich, Jean Utke, Olaf Vogel
11 
12 This file is part of ADOL-C. This software is provided as open source.
13 Any use, reproduction, or distribution of the software constitutes
14 recipient's acceptance of the terms of the accompanying license file.
15 
16---------------------------------------------------------------------------*/
17
18/****************************************************************************/
19/*                                                                 INCLUDES */
20#include <adolc.h>
21
22#include <math.h>
23#define PI 3.1415926536
24
25
26/****************************************************************************/
27/*                                                          ADOUBLE ROUTINE */
28adouble activeCubicLighthouse1( adouble t ) {
29    adouble p, q, d, r, u, u1,u2, v, c;
30    /*---------------------*/
31    p = tan(t);
32    q = p - 0.2;
33    p /= 3.0;
34    d = q*q;
35    d -= p*p*p;
36    /* 1. branch ----------*/
37    r = sqrt(d);
38    u = q + r;
39    u1 = pow(fabs(u),1.0/3.0);
40    u2 = -u1;
41    condassign(u,u,u1,u2);
42    v = q - r;
43    u1 = pow(fabs(v),1.0/3.0);
44    u2 = -u1;
45    condassign(v,v,u1,u2);
46    c = u + v;
47    /*---------------------*/
48    c += 2.0;
49    return c;
50}
51
52/****************************************************************************/
53/*                                                          ADOUBLE ROUTINE */
54adouble activeCubicLighthouse2( adouble t ) {
55    adouble p, q, d, r, u, v, c, a, z, b;
56    /*---------------------*/
57    p = tan(t);
58    q = p - 0.2;
59    p /= 3.0;
60    d = q*q;
61    d -= p*p*p;
62    /* 2. branch ----------*/
63    p = fabs(p);
64    p = sqrt(p);
65    q /= p*p*p;
66    a = acos(q);
67    a /= 3.0;
68    z = cos(a);
69    b = a + PI/3.0;
70    b = -cos(b);
71    z = fmin(z,b);
72    b = a - PI/3.0;
73    b = -cos(b);
74    z = fmin(z,b);
75    z = 2.0*z*p;
76    /*---------------------*/
77    z += 2.0;
78    return z;
79}
80
81/****************************************************************************/
82/*                                                             MAIN PROGRAM */
83int main() {
84    int i, vc;
85    int tag1 = 1, tag2 = 2;
86    double z, z1, z2, t, tmin, tmax, tdist, dz;
87
88    /*--------------------------------------------------------------------------*/
89    /* Preparation */
90    fprintf(stdout,"CUBIC LIGHTHOUSE Using CARDAN (ADOL-C Example)\n\n");
91    tmin = 0.15;
92    tmax = 0.24;
93    fprintf(stdout,"How many values = ? \n");
94    scanf("%d",&vc);
95
96    /*--------------------------------------------------------------------------*/
97    t = 0.1;
98    adouble az,at;
99    trace_on(tag1);
100    at <<= t;
101    az = activeCubicLighthouse1(at);
102    az >>= z;
103    trace_off();
104    trace_on(tag2);
105    at <<= t;
106    az = activeCubicLighthouse2(at);
107    az >>= z;
108    trace_off();
109
110    /*--------------------------------------------------------------------------*/
111    int tape_stats[STAT_SIZE];
112
113    tapestats(tag1,tape_stats);
114
115    fprintf(stdout,"\n    independents            %d\n",tape_stats[NUM_INDEPENDENTS]);
116    fprintf(stdout,"    dependents              %d\n",tape_stats[NUM_DEPENDENTS]);
117    fprintf(stdout,"    operations              %d\n",tape_stats[NUM_OPERATIONS]);
118    fprintf(stdout,"    operations buffer size  %d\n",tape_stats[OP_BUFFER_SIZE]);
119    fprintf(stdout,"    locations buffer size   %d\n",tape_stats[LOC_BUFFER_SIZE]);
120    fprintf(stdout,"    constants buffer size   %d\n",tape_stats[VAL_BUFFER_SIZE]);
121    fprintf(stdout,"    maxlive                 %d\n",tape_stats[NUM_MAX_LIVES]);
122    fprintf(stdout,"    valstack size           %d\n\n",tape_stats[TAY_STACK_SIZE]);
123
124    tapestats(tag2,tape_stats);
125
126    fprintf(stdout,"\n    independents            %d\n",tape_stats[NUM_INDEPENDENTS]);
127    fprintf(stdout,"    dependents              %d\n",tape_stats[NUM_DEPENDENTS]);
128    fprintf(stdout,"    operations              %d\n",tape_stats[NUM_OPERATIONS]);
129    fprintf(stdout,"    operations buffer size  %d\n",tape_stats[OP_BUFFER_SIZE]);
130    fprintf(stdout,"    locations buffer size   %d\n",tape_stats[LOC_BUFFER_SIZE]);
131    fprintf(stdout,"    constants buffer size   %d\n",tape_stats[VAL_BUFFER_SIZE]);
132    fprintf(stdout,"    maxlive                 %d\n",tape_stats[NUM_MAX_LIVES]);
133    fprintf(stdout,"    valstack size           %d\n\n",tape_stats[TAY_STACK_SIZE]);
134
135    /*--------------------------------------------------------------------------*/
136    tdist = (tmax-tmin)/((double) (vc-1));
137    t = tmin;
138    for (i=0; i<vc; i++) {
139        function(tag1,1,1,&t,&z1);
140        function(tag2,1,1,&t,&z2);
141        if (!(z1==z1)) // check for NaN
142        { gradient(tag2,1,&t,&dz);
143            fprintf(stdout,"%e %e %e\n",t,z2,dz);
144        } else {
145            gradient(tag1,1,&t,&dz);
146            fprintf(stdout,"%e %e %e\n",t,z1,dz);
147        }
148        t += tdist;
149    }
150
151    /*--------------------------------------------------------------------------*/
152    return 1;
153}
154
155
156
157
Note: See TracBrowser for help on using the repository browser.