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

Last change on this file since 348 was 348, checked in by kulshres, 7 years ago

correct compilation error due to condassign

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

  • Property svn:keywords set to Author Date Id Revision
File size: 7.9 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     hessmat.cpp
4 Revision: $Id: hessmat.cpp 348 2012-09-15 21:44:27Z kulshres $
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 */
19#include <adolc/adolc.h>
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],(adouble)(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.