source: trunk/ADOL-C/examples/additional_examples/timing/sfunc_helmholtz.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.3 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     sfunc_helmholtz.cpp
4 Revision: $Id: sfunc_helmholtz.cpp 42 2009-07-15 18:37:17Z awalther $
5 Contents: function module containing Helmholtz energy function
6
7   Each << function module >> contains:
8         
9     (1) const char* const controlFileName
10     (2) int indepDim;
11     (3) void initProblemParameters( void )
12     (4) void initIndependents( double* indeps )
13     (5) double originalScalarFunction( double* indeps )
14     (6) double tapingScalarFunction( int tag, double* indeps )   
15 
16 Copyright (c) Andrea Walther, Andreas Griewank, Andreas Kowarz,
17               Hristo Mitev, Sebastian Schlenkrich, Jean Utke, Olaf Vogel
18 
19 This file is part of ADOL-C. This software is provided as open source.
20 Any use, reproduction, or distribution of the software constitutes
21 recipient's acceptance of the terms of the accompanying license file.
22 
23---------------------------------------------------------------------------*/
24#define _SFUNC_HELMHOLTZ_C_
25
26
27/****************************************************************************/
28/*                                                                 INCLUDES */
29#include <adolc.h>
30#include <cmath>
31
32
33/****************************************************************************/
34/*                                                         GLOBAL VARIABLES */
35
36/*--------------------------------------------------------------------------*/
37/*                                                        Control file name */
38const char* controlFileName = "helmholtzexam.ctrl";
39
40/*--------------------------------------------------------------------------*/
41/*                                                               Dimensions */
42int indepDim;
43
44/*--------------------------------------------------------------------------*/
45/*                                       Other problem dependent parameters */
46double *bv = NULL;
47const double R = sqrt(2.0);
48const double TE= 0.01; /* originally 0.0 */
49
50
51/****************************************************************************/
52/*                                                  INIT PROBLEM PARAMETERS */
53void initProblemParameters( void ) {
54    fprintf(stdout,"HELMHOLTZ ENERGY (ADOL-C Example)\n\n");
55    if (indepDim <= 0) {
56        fprintf(stdout,"    number of independent variables = ? ");
57        fscanf(stdin,"%d",&indepDim);
58        fprintf(stdout,"\n");
59    }
60}
61
62
63/****************************************************************************/
64/*                                                        INITIALIZE INDEPs */
65void initIndependents( double* indeps ) {
66    int i;
67    double r = 1.0/indepDim;
68    if (bv)
69        delete[] bv;
70    bv = new double[indepDim];
71    for (i=0; i<indepDim; i++) {
72        indeps[i] = r*sqrt(1.0+i);
73        bv[i]     = 0.02*(1.0+fabs(sin(double(i))));
74    }
75}
76
77
78/****************************************************************************/
79/*                                                 ORIGINAL SCALAR FUNCTION */
80
81/*--------------------------------------------------------------------------*/
82/*                                                         Helmholtz energy */
83double helmholtz( int dim, double* indeps, double* bv ) {
84    int i,j;
85    double he;
86    double xax, bx, tem;
87    xax = 0;
88    bx  = 0;
89    he  = 0;
90    for (i=0; i<dim; i++) {
91        he += indeps[i]*log(indeps[i]);
92        bx += bv[i]*indeps[i];
93        tem = (2.0/(1.0+i+i))*indeps[i];
94        for (j=0; j<i; j++)
95            tem += (1.0/(1.0+i+j))*indeps[j];
96        xax += indeps[i]*tem;
97    }
98    xax *= 0.5;
99    he   = 1.3625E-3*(he-TE*log(1.0-bx));
100    he   = he - log((1+bx*(1+R))/(1+bx*(1-R)))*xax/bx;
101    return he;
102}
103
104/*--------------------------------------------------------------------------*/
105/*                                                   The interface function */
106double originalScalarFunction( double* indeps ) {
107    return helmholtz(indepDim, indeps, bv);
108}
109
110
111/****************************************************************************/
112/*                                                   TAPING SCALAR FUNCTION */
113
114/*--------------------------------------------------------------------------*/
115/*                                                  active Helmholtz energy */
116adouble activeHelmholtz( int dim, adouble* indeps, double* bv ) {
117    int i,j;
118    adouble he;
119    adouble xax, bx, tem;
120    xax = 0;
121    bx  = 0;
122    he  = 0;
123    for (i=0; i<dim; i++) {
124        he += indeps[i]*log(indeps[i]);
125        bx += bv[i]*indeps[i];
126        tem = (2.0/(1.0+i+i))*indeps[i];
127        for (j=0; j<i; j++)
128            tem += (1.0/(1.0+i+j))*indeps[j];
129        xax += indeps[i]*tem;
130    }
131    xax *= 0.5;
132    he   = 1.3625E-3*(he-TE*log(1.0-bx));
133    he   = he - log((1+bx*(1+R))/(1+bx*(1-R)))*xax/bx;
134    return he;
135}
136
137/*--------------------------------------------------------------------------*/
138/*                                                   The interface function */
139double tapingScalarFunction( int tag, double* indeps ) {
140    int i;
141    trace_on(tag);
142    adouble* activeIndeps = new adouble[indepDim];
143    adouble* aIP = activeIndeps;
144    double*  iP  = indeps;
145    for (i=0; i<indepDim; i++)
146        *aIP++ <<= *iP++;
147    adouble ares = activeHelmholtz(indepDim, activeIndeps, bv);
148    double res = 0;
149    ares >>= res;
150    trace_off();
151    return res;
152}
153
154#undef _SFUNC_HELMHOLTZ_C_
155
156
157
158
159
Note: See TracBrowser for help on using the repository browser.