source: trunk/ADOL-C/src/interfaces.cpp @ 708

Last change on this file since 708 was 542, checked in by kulshres, 5 years ago

make dvlparms.h totally internal

This file does not contain anything that the external
user interface requires. This is a totally internal
matter, and need not even be installed into ${includedir}

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

  • Property svn:keywords set to Author Date Id Revision
File size: 15.2 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     interfaces.cpp
4 Revision: $Id: interfaces.cpp 542 2014-08-15 14:11:25Z kulshres $
5 Contents: Genuine C++ Interfaces to ADOL-C forward & reverse calls.
6 
7 Copyright (c) Andrea Walther, Andreas Griewank, Andreas Kowarz,
8               Hristo Mitev, Sebastian Schlenkrich, Jean Utke, Olaf Vogel
9 
10 This file is part of ADOL-C. This software is provided as open source.
11 Any use, reproduction, or distribution of the software constitutes
12 recipient's acceptance of the terms of the accompanying license file.
13 
14----------------------------------------------------------------------------*/
15
16#include <adolc/interfaces.h>
17#include <adolc/adalloc.h>
18#include "dvlparms.h"
19
20/****************************************************************************/
21/*                                                                   MACROS */
22#define fabs(x) ((x) > 0 ? (x) : -(x))
23#define ceil(x) ((int)((x)+1) - (int)((x) == (int)(x)))
24
25extern "C" void adolc_exit(int errorcode, const char *what, const char* function, const char *file, int line);
26
27/****************************************************************************/
28/*                                           FORWARD MODE, overloaded calls */
29
30/****************************************************************************/
31/*                                                             general call */
32/*                                                                          */
33int forward( short  tag,
34             int    m,
35             int    n,
36             int    d,
37             int    keep,
38             double **X,
39             double **Y)
40/* forward(tag, m, n, d, keep, X[n][d+1], Y[m][d+1])                        */
41{ /* olvo 980729 general ec */
42    static double *x, *y, *xp, *yp;
43    static int maxn, maxm;
44    int rc = -1, i, k;
45    if (n > maxn) {
46        if (x)
47            myfree1(x);
48        if (xp)
49            myfree1(xp);
50        x  = myalloc1(maxn = n);
51        xp = myalloc1(maxn);
52    }
53    if (m > maxm) {
54        if (y)
55            myfree1(y);
56        if (yp)
57            myfree1(yp);
58        y  = myalloc1(maxm = m);
59        yp = myalloc1(maxm);
60    }
61
62    /*------------------------------------------------------------------------*/
63    /* prepare input */
64    for (i=0; i<n; i++) {
65        x[i] = X[i][0];
66        if (d == 1)
67            xp[i] = X[i][1];
68        else
69            for (k=0; k<d; k++)
70                X[i][k] = X[i][k+1];
71    }
72
73    /*------------------------------------------------------------------------*/
74    /* function calls */
75    if (d == 0)
76        rc = zos_forward(tag,m,n,keep,x,y);
77    else
78        if (d == 1)
79            rc = fos_forward(tag,m,n,keep,x,xp,y,yp);
80        else
81            rc = hos_forward(tag,m,n,d,keep,x,X,y,Y);
82
83    /*------------------------------------------------------------------------*/
84    /* prepare output */
85    for (i=0; i<n; i++)
86        if (d > 1) {
87            for (k=d; k>0; k--)
88                X[i][k] = X[i][k-1];
89            X[i][0] = x[i];
90        }
91
92    for (i=0; i<m; i++) {
93        if (d == 1)
94            Y[i][1] = yp[i];
95        else
96            for (k=d; k>0; k--)
97                Y[i][k] = Y[i][k-1];
98        Y[i][0] = y[i];
99    }
100
101    return rc;
102}
103
104
105/****************************************************************************/
106/*         Y can be one dimensional if m=1                                  */
107/*                                                                          */
108int forward( short  tag,
109             int    m,
110             int    n,
111             int    d,
112             int    keep,
113             double **X,
114             double *Y)
115/* forward(tag, 1, n, d, keep, X[n][d+1], Y[d+1]), m=1                      */
116{ /* olvo 980729 general ec */
117    static double *x, *xp;
118    static int maxn;
119    double y;
120    int rc= -1, i, k;
121
122    if (m == 1) {
123        if (n > maxn) {
124            if (x)
125                myfree1(x);
126            if (xp)
127                myfree1(xp);
128            x  = myalloc1(maxn = n);
129            xp = myalloc1(maxn);
130        }
131
132        /*----------------------------------------------------------------------*/
133        /* prepare input */
134        for (i=0; i<n; i++) {
135            x[i] = X[i][0];
136            if (d == 1)
137                xp[i] = X[i][1];
138            else
139                for (k=0; k<d; k++)
140                    X[i][k] = X[i][k+1];
141        }
142
143        /*----------------------------------------------------------------------*/
144        /* function calls */
145        if (d == 0)
146            rc = zos_forward(tag,m,n,keep,x,&y);
147        else
148            if (d == 1)
149                rc = fos_forward(tag,m,n,keep,x,xp,&y,Y);
150            else
151                rc = hos_forward(tag,m,n,d,keep,x,X,&y,&Y);
152
153        /*----------------------------------------------------------------------*/
154        /* prepare output */
155        for (i=0; i<n; i++)
156            if (d > 1) {
157                for (k=d; k>0; k--)
158                    X[i][k] = X[i][k-1];
159                X[i][0] = x[i];
160            }
161
162        for (k=d; k>0; k--)
163            Y[k] = Y[k-1];
164        Y[0] = y;
165    } else {
166        fprintf(DIAG_OUT,"ADOL-C error: wrong Y dimension in forward \n");
167        adolc_exit(-1,"",__func__,__FILE__,__LINE__);
168    }
169
170    return rc;
171}
172
173
174/****************************************************************************/
175/*         X and Y can be one dimensional if d = 0                          */
176/*                                                                          */
177int forward( short  tag,
178             int    m,
179             int    n,
180             int    d,
181             int    keep,
182             double *X,
183             double *Y)
184/* forward(tag, m, n, 0, keep, X[n], Y[m]), d=0                             */
185{ int rc = -1;
186
187    if (d != 0) {
188        fprintf(DIAG_OUT,"ADOL-C error:  wrong X and Y dimensions in forward \n");
189        adolc_exit(-1,"",__func__,__FILE__,__LINE__);
190    } else
191        rc = zos_forward(tag,m,n,keep,X,Y);
192
193    return rc;
194}
195
196
197/****************************************************************************/
198/*         X and Y can be one dimensional if d omitted                      */
199/*                                                                          */
200int forward(short  tag,
201            int    m,
202            int    n,
203            int    keep,
204            double *X,
205            double *Y)
206/* forward(tag, m, n, keep, X[n], Y[m])                                     */
207{ return zos_forward(tag,m,n,keep,X,Y);
208}
209
210
211/****************************************************************************/
212/*                                                             general call */
213/*                                                                          */
214int forward( short  tag,
215             int    m,
216             int    n,
217             int    d,
218             int    p,
219             double *x,
220             double ***X,
221             double *y,
222             double ***Y)
223/* forward(tag, m, n, d, p, x[n], X[n][p][d], y[m], Y[m][p][d])             */
224{ return hov_forward(tag,m,n,d,p,x,X,y,Y);
225}
226
227
228/****************************************************************************/
229/*                                                             general call */
230/*                                                                          */
231int forward( short  tag,
232             int    m,
233             int    n,
234             int    p,
235             double *x,
236             double **X,
237             double *y,
238             double **Y)
239/* forward(tag, m, n, p, x[n], X[n][p], y[m], Y[m][p])                      */
240{ return fov_forward(tag,m,n,p,x,X,y,Y);
241}
242
243
244/****************************************************************************/
245/*                                           REVERSE MODE, overloaded calls */
246
247/****************************************************************************/
248/*                                                             general call */
249/*                                                                          */
250int reverse( short  tag,
251             int    m,
252             int    n,
253             int    d,
254             double *u,
255             double **Z)
256/* reverse(tag, m, n, d, u[m], Z[n][d+1])                                   */
257{ return hos_reverse(tag,m,n,d,u,Z);
258}
259
260
261/****************************************************************************/
262/*         u can be a scalar if m=1                                         */
263/*                                                                          */
264int reverse( short  tag,
265             int    m,
266             int    n,
267             int    d,
268             double u,
269             double **Z)
270/* reverse(tag, 1, n, 0, u, Z[n][d+1]), m=1 => u scalar                     */
271{ int rc=-1;
272
273    if (m != 1) {
274        fprintf(DIAG_OUT,"ADOL-C error:  wrong u dimension in scalar-reverse \n");
275        adolc_exit(-1,"",__func__,__FILE__,__LINE__);
276    } else
277        rc = hos_reverse(tag,m,n,d,&u,Z);
278
279    return rc;
280}
281
282
283/****************************************************************************/
284/*         Z can be vector if d = 0; Done by specialized code               */
285/*                                                                          */
286int reverse( short  tag,
287             int    m,
288             int    n,
289             int    d,
290             double *u,
291             double *Z)
292/* reverse(tag, m, n, 0, u[m], Z[n]), d=0                                   */
293{ if (d != 0) {
294        fprintf(DIAG_OUT,"ADOL-C error:  wrong Z dimension in scalar-reverse \n");
295        adolc_exit(-1,"",__func__,__FILE__,__LINE__);
296    }
297
298    return fos_reverse(tag,m,n,u,Z);
299}
300
301
302/****************************************************************************/
303/*         u and Z can be scalars if m=1 and d=0;                           */
304/*                                                                          */
305int reverse( short  tag,
306             int    m,
307             int    n,
308             int    d,
309             double u,
310             double *Z)
311/* reverse(tag, 1, n, 0, u, Z[n]), m=1 and d=0 => u and Z scalars           */
312{ int rc=-1;
313
314    if (m != 1 || d != 0 ) {
315        fprintf(DIAG_OUT,"ADOL-C error:  wrong u or Z dimension in scalar-reverse \n");
316        adolc_exit(-1,"",__func__,__FILE__,__LINE__);
317    } else
318        rc = fos_reverse(tag,m,n,&u,Z);
319    \
320    return rc;
321}
322
323
324/****************************************************************************/
325/*                                                             general call */
326/*                                                                          */
327int reverse( short  tag,
328             int    m,
329             int    n,
330             int    d,
331             int    q,
332             double **U,
333             double ***Z,
334             short  **nz)
335/* reverse(tag, m, n, d, q, U[q][m], Z[q][n][d+1], nz[q][n])                */
336{ return hov_reverse(tag,m,n,d,q,U,Z,nz);
337}
338
339
340/****************************************************************************/
341/*         U can be a vector if m=1                                         */
342/*                                                                          */
343int reverse( short  tag,
344             int    m,
345             int    n,
346             int    d,
347             int    q,
348             double *U,
349             double ***Z,
350             short  **nz)
351/* reverse(tag, 1, n, d, q, U[q], Z[q][n][d+1], nz[q][n]), m=1 => u vector  */
352{ int rc=-1;
353
354    if (m != 1) {
355        fprintf(DIAG_OUT,"ADOL-C error:  wrong U dimension in vector-reverse \n");
356        adolc_exit(-1,"",__func__,__FILE__,__LINE__);
357    } else { /* olvo 980727 ??? */
358        /* double** upp = new double*[nrows]; */
359        double **upp = (double**) malloc(q*sizeof(double*));
360        for (int i=0; i<q; i++)
361            upp[i] = &U[i];
362        rc=hov_reverse(tag,m,n,d,q,upp,Z,nz);
363        /* delete[] upp; */
364        free((char*)upp);
365    }
366    return rc;
367}
368
369
370/****************************************************************************/
371/*                                                                          */
372/*         If d=0 then Z may be matrix; Done by specialized code            */
373/*                                                                          */
374int reverse( short  tag,
375             int    m,
376             int    n,
377             int    d,
378             int    q,
379             double **U,
380             double **Z)
381/* reverse(tag, m, n, 0, q, U[q][m], Z[q][n]), d=0 => Z matrix              */
382{ int rc=-1;
383
384    if (d != 0) {
385        fprintf(DIAG_OUT,"ADOL-C error:  wrong degree in vector-reverse \n");
386        adolc_exit(-1,"",__func__,__FILE__,__LINE__);
387    } else
388        rc = fov_reverse(tag,m,n,q,U,Z);
389
390    return rc;
391}
392
393
394/****************************************************************************/
395/*                                                                          */
396/*         d=0 may be omitted, then Z may be a matrix; specialized code     */
397/*                                                                          */
398int reverse( short  tag,
399             int    m,
400             int    n,
401             int    q,
402             double **U,
403             double **Z)
404/* reverse(tag, m, n, q, U[q][m], Z[q][n]), d=0 => Z matrix                 */
405{ int rc=-1;
406
407    rc = fov_reverse(tag,m,n,q,U,Z);
408
409    return rc;
410}
411
412
413/****************************************************************************/
414/*                                                                          */
415/*         If m=1 and d=0 then U can be vector and Z a matrix but no nz.    */
416/*                                                                          */
417int reverse( short  tag,
418             int    m,
419             int    n,
420             int    d,
421             int    q,
422             double *U,
423             double **Z)
424/* reverse(tag, 1, n, 0, q, U[q], Z[q][n]),
425                            m=1 and d=0 => U vector and Z matrix but no nz  */
426{ int rc=-1;
427
428    /* olvo 981126 ??? what's that: */
429    /* (++d)--; *//* degre is reserved for the future use. Ingore this line */
430
431    if ((m != 1) || (d != 0)) {
432        fprintf(DIAG_OUT,"ADOL-C error:  wrong U dimension in vector-reverse \n");
433        adolc_exit(-1,"",__func__,__FILE__,__LINE__);
434    } else { /* olvo 980727 ??? */
435        /* double ** upp = new double*[nrows]; */
436        double **upp = (double**) malloc(q*sizeof(double*));
437        for (int i=0; i<q; i++)
438            upp[i] = &U[i];
439        rc = fov_reverse(tag,m,n,q,upp,Z);
440        /* delete[] upp; */
441        free((char*) upp);
442    }
443
444    return rc;
445}
446
447
448/****************************************************************************/
449/*                                                                          */
450/*         If p and U are omitted they default to m and I so that as above  */
451/*                                                                          */
452int reverse( short  tag,
453             int    m,
454             int    n,
455             int    d,
456             double ***Z,
457             short  **nz)
458/* reverse(tag, m, n, d, Z[p][n][d+1], nz[p][n]),
459           If p and U are omitted they default to m and I                   */
460{ static int depax;
461    static double** I;
462    if (m adolc_compsize depax) {
463        if (depax)
464            myfreeI2(depax,I);
465        I = myallocI2(depax = m);
466    }
467    return hov_reverse(tag,m,n,d,m,I,Z,nz);
468}
469
470/****************************************************************************/
471/*                                                               THAT'S ALL */
Note: See TracBrowser for help on using the repository browser.