source: trunk/ADOL-C/examples/additional_examples/helm/helm-diff-exam.cpp @ 624

Last change on this file since 624 was 88, checked in by utke, 10 years ago

compiler message: ISO C++ does not support the `%le' printf format

  • 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:     helm-diff-exam.cpp
4 Revision: $Id: helm-diff-exam.cpp 88 2010-02-17 18:33:52Z kulshres $
5 Contents: example for  Helmholtz energy example
6           Computes gradient using divide differences
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 <sys/types.h>
20#include <stdio.h>
21
22#include <math.h>
23#include <cstdlib>
24
25
26/****************************************************************************/
27/*                                                    CONSTANTS & VARIABLES */
28#define delta 0.000001
29#define TE    0.01
30#define R     sqrt(2.0)
31
32
33/****************************************************************************/
34/*                                                         HELMHOLTZ ENERGY */
35double energy( int n, double x[], double bv[] ) {
36    double he, xax, bx, tem;
37    int i,j;
38    xax = 0;
39    bx  = 0;
40    he  = 0;
41    for (i=0; i<n; i++) {
42        he += x[i]*log(x[i]);
43        bx +=  bv[i]*x[i];
44        tem = (2.0/(1.0+i+i))*x[i];
45        for (j=0; j<i; j++)
46            tem += (1.0/(1.0+i+j))*x[j];
47        xax += x[i]*tem;
48    }
49    xax *= 0.5;
50    he   = 1.3625E-3*(he-TE*log(1.0-bx));
51    he   = he - log((1+bx*(1+R))/(1+bx*(1-R)))*xax/bx;
52    return he;
53}
54
55
56/****************************************************************************/
57/*                                                                     MAIN */
58/*
59   This program computes first order directional derivatives
60   for the helmholtz energy function */
61int main(int argc, char *argv[]) {
62    int nf, n, j, l;
63    double result1, result2;
64    double q, jd, r;
65    double *x, *bv;
66
67    fprintf(stdout,"HELM-DIFF-EXAM (ADOL-C Example)\n\n");
68    fprintf(stdout," # of independents/10 =? \n ");
69    scanf("%d",&nf);
70
71    /*--------------------------------------------------------------------------*/
72    n = 10 * nf;                                            /* Initilizations */
73    x   = (double*) malloc(n*sizeof(double));
74    bv  = (double*) malloc(n*sizeof(double));
75
76    r = 1.0/n;
77    for (j=0; j<n; j++) {
78        jd     = j;
79        bv[j]  = 0.02*(1.0+fabs(sin(jd)));
80        x[j]   = r*sqrt(1.0+jd);
81    }
82
83    /*--------------------------------------------------------------------------*/
84    result2 = energy(n,x,bv);                                    /* basepoint */
85    fprintf(stdout,"%14.6E -- energy\n",result2);
86
87    /*--------------------------------------------------------------------------*/
88    for (l=0; l<n; l++)                            /* directional derivatives */
89    { x[l]    = x[l]+delta;
90        result1 = energy(n,x,bv);
91        x[l]    = x[l]-delta;
92        q       = (result1-result2)/delta;
93        fprintf(stdout,"%3d: %14.6E,  \n",l,q);
94    }
95    fprintf(stdout,"%14.6E -- energy\n",result2);
96
97    free((char*) bv);
98    free((char*) x);
99
100    return 0;
101}
102
103
104/****************************************************************************/
105/*                                                               THAT'S ALL */
106
Note: See TracBrowser for help on using the repository browser.