source: trunk/ADOL-C/examples/additional_examples/lighthouse/cubic-2.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: 3.8 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     cubic-2.cpp
4 Revision: $Id: cubic-2.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 conditionals
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 <adolc.h>
20
21#include <math.h>
22#define PI 3.1415926536
23
24
25/****************************************************************************/
26/*                                                          ADOUBLE ROUTINE */
27adouble activeCubicLighthouse( adouble t ) {
28    adouble p, q, d, r, u, u1, u2, v, a, b, c, z;
29    /*---------------------*/
30    p = tan(t);
31    q = p - 0.2;
32    p /= 3.0;
33    d = q*q;
34    d -= p*p*p;
35    /* 1. branch ----------*/
36    r = sqrt(d);
37    u = q + r;
38    u1 = pow(fabs(u),1.0/3.0);
39    u2 = -u1;
40    condassign(u,u,u1,u2);
41    v = q - r;
42    u1 = pow(fabs(v),1.0/3.0);
43    u2 = -u1;
44    condassign(v,v,u1,u2);
45    c = u + v;
46    /* 2. branch ----------*/
47    p = fabs(p);
48    p = sqrt(p);
49    q /= p*p*p;
50    a = acos(q);
51    a /= 3.0;
52    z = cos(a);
53    b = a + PI/3.0;
54    b = -cos(b);
55    z = fmin(z,b);
56    b = a - PI/3.0;
57    b = -cos(b);
58    z = fmin(z,b);
59    z = 2.0*z*p;
60    /*---------------------*/
61    condassign(z,d,c);
62    z += 2.0;
63    return z;
64}
65
66/****************************************************************************/
67/*                                                             MAIN PROGRAM */
68int main() {
69    int i, vc;
70    int tag = 1;
71    double z, t, tmin, tmax, tdist, dz;
72
73    /*--------------------------------------------------------------------------*/
74    /* Preparation */
75    fprintf(stdout,"CUBIC LIGHTHOUSE Using CARDAN (ADOL-C Example)\n\n");
76    tmin = 0.15;
77    tmax = 0.24;
78    fprintf(stdout,"How many values = ? \n");
79    scanf("%d",&vc);
80
81    /*--------------------------------------------------------------------------*/
82    t = 0.1;
83    adouble az,at;
84    trace_on(tag);
85    at <<= t;
86    az = activeCubicLighthouse(at);
87    az >>= z;
88    trace_off();
89
90    /*--------------------------------------------------------------------------*/
91    int tape_stats[STAT_SIZE];
92    tapestats(tag,tape_stats);
93
94    fprintf(stdout,"\n    independents            %d\n",tape_stats[NUM_INDEPENDENTS]);
95    fprintf(stdout,"    dependents              %d\n",tape_stats[NUM_DEPENDENTS]);
96    fprintf(stdout,"    operations              %d\n",tape_stats[NUM_OPERATIONS]);
97    fprintf(stdout,"    operations buffer size  %d\n",tape_stats[OP_BUFFER_SIZE]);
98    fprintf(stdout,"    locations buffer size   %d\n",tape_stats[LOC_BUFFER_SIZE]);
99    fprintf(stdout,"    constants buffer size   %d\n",tape_stats[VAL_BUFFER_SIZE]);
100    fprintf(stdout,"    maxlive                 %d\n",tape_stats[NUM_MAX_LIVES]);
101    fprintf(stdout,"    valstack size           %d\n\n",tape_stats[TAY_STACK_SIZE]);
102
103
104    /*--------------------------------------------------------------------------*/
105    tdist = (tmax-tmin)/((double) (vc-1));
106    t = tmin;
107    for (i=0; i<vc; i++) {
108        function(tag,1,1,&t,&z);
109        gradient(tag,1,&t,&dz);
110        fprintf(stdout,"%e %e %e\n",t,z,dz);
111        t += tdist;
112    }
113
114    /*--------------------------------------------------------------------------*/
115    return 1;
116}
117
118
119
120
Note: See TracBrowser for help on using the repository browser.