source: trunk/ADOL-C/examples/additional_examples/hessmat/hessmat.cpp @ 171

Last change on this file since 171 was 171, checked in by kulshres, 9 years ago

Squashed merge branch 'master' of 'gitclone' into svn

  • 'master' of 'gitclone': (84 commits) adjust example makefiles and include paths get rid of the symlink in the src subdirectory

details of the commits:
commit c9e4bc332d2363f737fc2e8a8fcfc2e43ddb9d15
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Mon Oct 4 15:43:47 2010 +0200

adjust example makefiles and include paths

include paths in example sources were wrong for some time now
simplify makefile rules too, there is really no need for checking SPARSE
adjust include paths in makefiles.

Signed-off-by: Kshitij Kulshreshtha <kshitij@…>

commit e6e1963e41e097fd5b4a79cd1611c12f6868dc94
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Mon Oct 4 15:41:25 2010 +0200

get rid of the symlink in the src subdirectory

windows doesn't like symlinks and make infinite depth directories
we now create a symlink for build in the directory parallel to src
adjust all makefiles.am accordingly for build

Signed-off-by: Kshitij Kulshreshtha <kshitij@…>

  • Property svn:keywords set to Author Date Id Revision
File size: 7.9 KB
RevLine 
[40]1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     hessmat.cpp
[42]4 Revision: $Id: hessmat.cpp 171 2010-10-04 13:57:19Z kulshres $
[40]5 Contents: example for testing the routines:
6           hov_wk_forward    ( = Higher Order Vector forward With Keep )
7           hos_ov_reverse    ( = Higher Order Scalar reverse over vectors)
8
9 Copyright (c) Andrea Walther, Andreas Kowarz, 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 */
[171]19#include <adolc/adolc.h>
[40]20
21#include <stdlib.h>
22#include <iostream>
23using namespace std;
24
25/****************************************************************************/
26/*                                                                     MAIN */
27int main() {
28    int i,j,l,m,n,d,q,bd, keep;
29
30    /*--------------------------------------------------------------------------*/
31    /* inputs */
32    cout << "vector x Hessian x matrix for the function \n\n";
33    cout << " y[0] = cos(x[0])* ...*cos(x[n]) \n";
34    cout << " y[1] = x[0]^n \n";
35    cout << " y[2] = condassign(y[i],y[0]>y[1],y[1],y[0]) \n";
36    cout << " y[3] = sin(x[0])+ ...+sin(x[n]) \n";
37    cout << " y[4] = exp(x[0])- ...-exp(x[n]) \n";
38    cout << " y[5] = pow(y[1],3) \n";
39    cout << " y[6] += y[5]*y[4] \n";
40    cout << " y[7] -= y[6]*y[5] \n";
41    cout << " y[j] = 1/x[0]/ .../x[n], j > 3 \n\n";
42
43    cout << " Number of independents = ?\n ";
44    cin >> n;
45    cout << " Number of dependents =  ?\n ";
46    cin >> m;
47    cout << " Degree d (for forward) =  ?\n";
48    cin >> d;
49    cout << " keep (degree of corresponding reverse = keep-1) =  ?\n";
50    cout << "       keep <= d+1 must be valid \n";
51    cin >> keep;
52    cout << " Number of directions =  ?\n ";
53    cin >> q;
54
55    /*--------------------------------------------------------------------------*/
56    /* allocations and inits */
57
58    double* xp = new double[n];                      /* passive indeps        */
59    double* yp = new double[m];                      /* passive depends       */
60
61    /* vector x Hessian x matrix = Upp x H x XPPP */
62
63    double* Up = myalloc(m);        /* vector on left-hand side                */
64    double** Upp = myalloc(m,d+1);   /* vector on left-hand side                */
65    double*** Xppp = myalloc(n,q,d); /* matrix on right-hand side             */
66    double*** Zppp = myalloc(q,n,d+1); /* result of Up x H x XPPP             */
67
68    double*** Yppp = myalloc(m,q,d); /* results of needed hos_wk_forward      */
69
70    /* check results with usual lagra-Hess-vec  */
71
72    double** Xpp = myalloc(n,d);
73    double** V   = myalloc(n,q);
74    double** W   = myalloc(q,n);
75    double** H   = myalloc(n,n);
76    double** Ypp = myalloc(m,d);
77    double** Zpp = myalloc(n,d+1);
78    /* inits */
79
80    for (l=0; l<d; l++)                    /* first everything is set to zero */
81        for (i=0; i<n; i++)
82            for (j=0;j<q;j++)
83                Xppp[i][j][l] = 0;
84
85    /* now carthesian directions as choosen as    */
86    /* matrix on right-hand side of Up x H x XPPP */
87    bd = (n<q)?n:q;
88    for (j=0;j<bd;j++)
89        Xppp[j][j][0] = 1;
90
91    for (i=0; i<m; i++)         /* vector on left-hand side of Up x H x XPPP  */
92    {Up[i] = 1;                /* is initialised with 1's                    */
93        Upp[i][0] = 1;
94        for (j=1;j<=d;j++)
95            Upp[i][j] = 0;
96    }
97
98    for (i=0; i<n; i++)                    /* first everything is set to zero */
99        for (j=0;j<d;j++)
100            Xpp[i][j] = 0;
101    Xpp[0][0] = 1;                 /* now one carthesian direction as choosen */
102    /* as vector for lagra-Hess-vec            */
103
104    for (i=0; i<n; i++)                            /* inits of passive indeps */
105        xp[i] = (i+1.0)/(2.0+i);
106
107    for (i=0; i<n; i++) {
108        for (j=0;j<q;j++)
109            V[i][j] = 0;
110        if (i < q)
111            V[i][i] = 1;
112    }
113
114    /*--------------------------------------------------------------------------*/
115    trace_on(1);                                      /* tracing the function */
116
117    adouble* x = new adouble[n];                     /* active indeps         */
118    adouble* y = new adouble[m];                     /* active depends        */
119    for(i=0;i<m;i++)
120        y[i] = 1;
121    for (i=0; i<n; i++) {
122        x[i] <<= xp[i];
123        y[0] *= cos(x[i]);
124    }
125
126    for(i=1;i<m;i++)
127        for(j=0;j<n;j++) {
128            switch (i) {
129                case 1 :
130                    y[i] *= x[0];
131                    break;
132                case 2 :
133                    condassign(y[i],y[0]>y[1],y[1],y[0]);
134                    break;
135                case 3 :
136                    y[i] -= sin(x[j]);
137                    break;
138                case 4 :
139                    y[i] -= exp(x[j]);
140                    break;
141                case 5 :
142                    y[5] = pow(y[1],3);
143                case 6 :
144                    y[6] += y[5]*y[4];
145                case 7 :
146                    y[7] -= y[6]*y[5];
147                default :
148                    y[i] /= x[j];
149            }
150        }
151    for (i=0; i<m; i++)
152        y[i] >>= yp[i] ;
153    trace_off();
154
155
156    /*--------------------------------------------------------------------------*/
157    /* work on the tape */
158
159
160    /* compute results of lagra_hess_vec */
161    /* the following is equal to calls inside of lagra_hess_vec(..) */
162    /* direct calls to the basic routines hos_forward and hos_reverse */
163    /* seem to be faster than call of lagra_hess_vec(..) */
164    /* at least in some of our test cases */
165
166    hos_forward(1,m,n,d,keep,xp,Xpp,yp,Ypp);
167    hos_reverse(1,m,n,keep-1,Up,Zpp);
168
169    printf("\n Results of hos_reverse:\n\n");
170
171    for (i=0; i<=d; i++) {
172        printf(" d = %d \n",i);
173        for (j=0;j<n;j++)
174            printf(" %6.3f ",Zpp[j][i]);
175        printf("\n");
176    }
177
178    /* The new drivers. First, hov_wk_forward(..) is called.
179       So far, it was impossible to store the results of
180       a higher-order-vector (=hov) forward in order to perform
181       a corresponding reverse sweep (for no particular reason.
182       Now we have hov with keep (=wk) and the results needed on
183       the way back are stored in a specific tape */
184
185    hov_wk_forward(1,m,n,d,keep,q,xp,Xppp,yp,Yppp);
186
187    /* The corresponding reverse sweep
188       So far we had only a higher-order-scalar (=hos, scalar because
189       only one vector on the left-hand-side) for a scalar forward
190       call.
191       Now, we use the stored vector information (= hos vector)
192       to compute multiple lagra_hess_vec at once */
193
194    hos_ov_reverse(1,m,n,keep-1,q,Upp,Zppp);
195
196    printf("\n Results of hosv_reverse:\n");
197
198    for (l=0; l<q; l++) {
199        for (i=0; i<=d; i++) {
200            printf(" d = %d \n",i);
201            for (j=0;j<n;j++)
202                printf(" %6.3f ",Zppp[l][j][i]);
203            printf("\n");
204        }
205        printf("\n\n");
206    }
207
208    if (m==1) {
209        printf("hess_mat:\n");
210        hess_mat(1,n,q,xp,V,W);
211        for (i=0; i<q; i++) {
212            for (j=0;j<n;j++)
213                printf(" %6.3f ",W[i][j]);
214            printf("\n");
215        }
216        printf("hessian2:\n");
217        hessian2(1,n,xp,H);
218        for (i=0; i<n; i++) {
219            for (j=0;j<n;j++)
220                printf(" %6.3f ",H[i][j]);
221            printf("\n");
222        }
223    }
224
225    myfree(Zpp);
226    myfree(Ypp);
227    myfree(H);
228    myfree(W);
229    myfree(V);
230    myfree(Xpp);
231    myfree(Zppp);
232    myfree(Yppp);
233    myfree(Xppp);
234    myfree(yp);
235    myfree(xp);
236    myfree(Up);
237    return 1;
238}
239
240
241/****************************************************************************/
242/*                                                               THAT'S ALL */
243
244
245
246
Note: See TracBrowser for help on using the repository browser.