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
25 extern "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 /* */
forward(short tag,int m,int n,int d,int keep,double ** X,double ** Y)33 int 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 /* */
forward(short tag,int m,int n,int d,int keep,double ** X,double * Y)108 int 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 /* */
forward(short tag,int m,int n,int d,int keep,double * X,double * Y)177 int 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 /* */
forward(short tag,int m,int n,int keep,double * X,double * Y)200 int 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 /* */
forward(short tag,int m,int n,int d,int p,double * x,double *** X,double * y,double *** Y)214 int 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 /* */
forward(short tag,int m,int n,int p,double * x,double ** X,double * y,double ** Y)231 int 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 /* */
reverse(short tag,int m,int n,int d,double * u,double ** Z)250 int 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 /* */
reverse(short tag,int m,int n,int d,double u,double ** Z)264 int 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 /* */
reverse(short tag,int m,int n,int d,double * u,double * Z)286 int 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 /* */
reverse(short tag,int m,int n,int d,double u,double * Z)305 int 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 /* */
reverse(short tag,int m,int n,int d,int q,double ** U,double *** Z,short ** nz)327 int 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 /* */
reverse(short tag,int m,int n,int d,int q,double * U,double *** Z,short ** nz)343 int 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 /* */
reverse(short tag,int m,int n,int d,int q,double ** U,double ** Z)374 int 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 /* */
reverse(short tag,int m,int n,int q,double ** U,double ** Z)398 int 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 /* */
reverse(short tag,int m,int n,int d,int q,double * U,double ** Z)417 int 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 /* */
reverse(short tag,int m,int n,int d,double *** Z,short ** nz)452 int 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 */
472