source: trunk/ADOL-C/src/drivers/drivers.c @ 61

Last change on this file since 61 was 61, checked in by awalther, 10 years ago

inclusion of error function for gcc compiler

  • Property svn:keywords set to Author Date Id Revision
File size: 9.7 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     drivers/drivers.c
4 Revision: $Id: drivers.c 61 2009-12-07 14:49:34Z awalther $
5 Contents: Easy to use drivers for optimization and nonlinear equations
6           (Implementation of the C/C++ callable interfaces).
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#include <drivers/drivers.h>
17#include <interfaces.h>
18#include <adalloc.h>
19
20#include <math.h>
21
22BEGIN_C_DECLS
23
24/****************************************************************************/
25/*                         DRIVERS FOR OPTIMIZATION AND NONLINEAR EQUATIONS */
26
27/*--------------------------------------------------------------------------*/
28/*                                                                 function */
29/* function(tag, m, n, x[n], y[m])                                          */
30int function(short tag,
31             int m,
32             int n,
33             double* argument,
34             double* result) {
35    int rc= -1;
36
37    rc= zos_forward(tag,m,n,0,argument,result);
38
39    return rc;
40}
41
42/*--------------------------------------------------------------------------*/
43/*                                                                 gradient */
44/* gradient(tag, n, x[n], g[n])                                             */
45int gradient(short tag,
46             int n,
47             const double* argument,
48             double* result) {
49    int rc= -1;
50    double one = 1.0;
51
52    rc = zos_forward(tag,1,n,1,argument,result);
53    if(rc < 0)
54        return rc;
55    MINDEC(rc, fos_reverse(tag,1,n,&one,result));
56    return rc;
57}
58
59/*--------------------------------------------------------------------------*/
60/*                                                                          */
61/* vec_jac(tag, m, n, repeat, x[n], u[m], v[n])                             */
62int vec_jac(short tag,
63            int m,
64            int n,
65            int repeat,
66            double* argument,
67            double* lagrange,
68            double* row) {
69    int rc= -1;
70    double *y = NULL;
71
72    if(!repeat) {
73        y = myalloc1(m);
74        rc = zos_forward(tag,m,n,1, argument, y);
75        if(rc < 0) return rc;
76    }
77    MINDEC(rc, fos_reverse(tag,m,n,lagrange,row));
78    if (!repeat) myfree1(y);
79    return rc;
80}
81
82/*--------------------------------------------------------------------------*/
83/*                                                                 jacobian */
84/* jacobian(tag, m, n, x[n], J[m][n])                                       */
85
86int jacobian(short tag,
87             int depen,
88             int indep,
89             double *argument,
90             double **jacobian) {
91    int rc;
92    double *result, **I;
93
94    result = myalloc1(depen);
95
96    if (indep/2 < depen) {
97        I = myallocI2(indep);
98        rc = fov_forward(tag,depen,indep,indep,argument,I,result,jacobian);
99        myfreeI2(indep, I);
100    } else {
101        I = myallocI2(depen);
102        rc = zos_forward(tag,depen,indep,1,argument,result);
103        if (rc < 0) return rc;
104        MINDEC(rc,fov_reverse(tag,depen,indep,depen,I,jacobian));
105        myfreeI2(depen, I);
106    }
107    return rc;
108}
109
110/*--------------------------------------------------------------------------*/
111/*                                                           large_jacobian */
112/* large_jacobian(tag, m, n, k, x[n], y[m], J[m][n])                        */
113
114int large_jacobian(short tag,
115                   int depen,
116                   int indep,
117                   int runns,
118                   double *argument,
119                   double *result,
120                   double **jacobian)
121{
122    int rc, dirs, i;
123    double **I;
124
125        I = myallocI2(indep);
126    if (runns > indep) runns = indep;
127    if (runns < 1)     runns = 1;
128    dirs = indep / runns;
129    if (indep % runns) ++dirs;
130    for (i=0; i<runns-1; ++i) {
131        rc = fov_offset_forward(tag, depen, indep, dirs, i * dirs, argument,
132                I, result, jacobian);
133    }
134    dirs = indep - (runns-1) * dirs;
135    rc = fov_offset_forward(tag, depen, indep, dirs, indep-dirs, argument,
136                     I, result, jacobian);
137    myfreeI2(indep, I);
138    return rc;
139}
140
141/*--------------------------------------------------------------------------*/
142/*                                                                  jac_vec */
143/* jac_vec(tag, m, n, x[n], v[n], u[m]);                                    */
144int jac_vec(short tag,
145            int m,
146            int n,
147            double* argument,
148            double* tangent,
149            double* column) {
150    int rc= -1;
151    double *y;
152
153    y = myalloc1(m);
154
155    rc = fos_forward(tag, m, n, 0, argument, tangent, y, column);
156    myfree1(y);
157
158    return rc;
159}
160
161/*--------------------------------------------------------------------------*/
162/*                                                                 hess_vec */
163/* hess_vec(tag, n, x[n], v[n], w[n])                                       */
164int hess_vec(short tag,
165             int n,
166             double *argument,
167             double *tangent,
168             double *result) {
169    double one = 1.0;
170    return lagra_hess_vec(tag,1,n,argument,tangent,&one,result);
171}
172
173/*--------------------------------------------------------------------------*/
174/*                                                                 hess_mat */
175/* hess_mat(tag, n, q, x[n], V[n][q], W[n][q])                              */
176int hess_mat(short tag,
177             int n,
178             int q,
179             double *argument,
180             double **tangent,
181             double **result) {
182    int rc;
183    int i,j;
184    double y;
185    double*** Xppp;
186    double*** Yppp;
187    double*** Zppp;
188    double**  Upp;
189
190    Xppp = myalloc3(n,q,1);   /* matrix on right-hand side  */
191    Yppp = myalloc3(1,q,1);   /* results of hos_wk_forward  */
192    Zppp = myalloc3(q,n,2);   /* result of Up x H x XPPP */
193    Upp  = myalloc2(1,2);     /* vector on left-hand side */
194
195    for (i = 0; i < n; ++i)
196        for (j = 0; j < q; ++j)
197            Xppp[i][j][0] = tangent[i][j];
198
199    Upp[0][0] = 1;
200    Upp[0][1] = 0;
201
202    rc = hov_wk_forward(tag, 1, n, 1, 2, q, argument, Xppp, &y, Yppp);
203    MINDEC(rc, hos_ov_reverse(tag, 1, n, 1, q, Upp, Zppp));
204
205    for (i = 0; i < q; ++i)
206        for (j = 0; j < n; ++j)
207            result[j][i] = Zppp[i][j][1];
208
209    myfree2(Upp);
210    myfree3(Zppp);
211    myfree3(Yppp);
212    myfree3(Xppp);
213
214    return rc;
215}
216
217/*--------------------------------------------------------------------------*/
218/*                                                                  hessian */
219/* hessian(tag, n, x[n], lower triangle of H[n][n])                         */
220/* uses Hessian-vector product                                              */
221int hessian(short tag,
222            int n,
223            double* argument,
224            double** hess) {
225    int rc= 3;
226    int i,j;
227    double *v = myalloc1(n);
228    double *w = myalloc1(n);
229    for(i=0;i<n;i++) v[i] = 0;
230    for(i=0;i<n;i++) {
231        v[i] = 1;
232        MINDEC(rc, hess_vec(tag, n, argument, v, w));
233        if( rc < 0) {
234            free((char *)v);
235            free((char *) w);
236            return rc;
237        }
238        for(j=0;j<=i;j++)
239            hess[i][j] = w[j];
240        v[i] = 0;
241    }
242
243    free((char *)v);
244    free((char *) w);
245    return rc;
246    /* Note that only the lower triangle of hess is filled */
247}
248
249/*--------------------------------------------------------------------------*/
250/*                                                                 hessian2 */
251/* hessian2(tag, n, x[n], lower triangle of H[n][n])                        */
252/* uses Hessian-matrix product                                              */
253int hessian2(short tag,
254             int n,
255             double* argument,
256             double** hess) {
257    int rc;
258    int i,j;
259
260    double*** Xppp = myalloc3(n,n,1);   /* matrix on right-hand side  */
261    double*   y    = myalloc1(1);       /* results of function evaluation */
262    double*** Yppp = myalloc3(1,n,1);   /* results of hos_wk_forward  */
263    double*** Zppp = myalloc3(n,n,2);   /* result of Up x H x XPPP */
264    double**  Upp  = myalloc2(1,2);     /* vector on left-hand side */
265
266    for (i=0; i<n; i++) {
267        for (j=0;j<n;j++)
268            Xppp[i][j][0] = 0;
269        Xppp[i][i][0] = 1;
270    }
271
272    Upp[0][0] = 1;
273    Upp[0][1] = 0;
274
275    rc = hov_wk_forward(tag,1,n,1,2,n,argument,Xppp,y,Yppp);
276    MINDEC(rc,hos_ov_reverse(tag,1,n,1,n,Upp,Zppp));
277
278    for (i=0; i<n; i++)
279        for (j=0;j<=i;j++)
280            hess[i][j] = Zppp[i][j][1];
281
282    myfree2(Upp);
283    myfree3(Zppp);
284    myfree3(Yppp);
285    myfree1(y);
286    myfree3(Xppp);
287    return rc;
288    /* Note that only the lower triangle of hess is filled */
289}
290
291/*--------------------------------------------------------------------------*/
292/*                                                           lagra_hess_vec */
293/* lagra_hess_vec(tag, m, n, x[n], v[n], u[m], w[n])                        */
294int lagra_hess_vec(short tag,
295                   int m,
296                   int n,
297                   double *argument,
298                   double *tangent,
299                   double *lagrange,
300                   double *result) {
301    int rc=-1;
302    int i;
303    int degree = 1;
304    int keep = degree+1;
305    double **X, *y, *y_tangent;
306
307    X = myalloc2(n,2);
308    y = myalloc1(m);
309    y_tangent = myalloc1(m);
310
311    rc = fos_forward(tag, m, n, keep, argument, tangent, y, y_tangent);
312
313    if(rc < 0) return rc;
314
315    MINDEC(rc, hos_reverse(tag, m, n, degree, lagrange, X));
316
317    for(i = 0; i < n; ++i)
318        result[i] = X[i][1];
319
320    myfree1(y_tangent);
321    myfree1(y);
322    myfree2(X);
323
324    return rc;
325}
326
327END_C_DECLS
Note: See TracBrowser for help on using the repository browser.