1 /*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File: sgenmain.cpp
4 Revision: $Id: sgenmain.cpp 511 2014-05-14 13:01:01Z kulshres $
5 Contents: Scalar Generic Main File:
6 for use with function modules containing several scalar
7 examples
8 (e.g. the determinant example in sfunc_determinant.cpp)
9
10 Each << function module >> contains:
11
12 (1) const char* const controlFileName
13 (2) int indepDim;
14 (3) void initProblemParameters( void )
15 (4) void initIndependents( double* indeps )
16 (5) double originalScalarFunction( double* indeps )
17 (6) double tapingScalarFunction( int tag, double* indeps )
18
19 Copyright (c) Andrea Walther, Andreas Griewank, Andreas Kowarz,
20 Hristo Mitev, Sebastian Schlenkrich, Jean Utke, Olaf Vogel
21
22 This file is part of ADOL-C. This software is provided as open source.
23 Any use, reproduction, or distribution of the software constitutes
24 recipient's acceptance of the terms of the accompanying license file.
25
26 ---------------------------------------------------------------------------*/
27 #define _SGENMAIN_C_
28
29 /****************************************************************************/
30 /* INCLUDES */
31 #include <adolc/adolc.h>
32 #include "../clock/myclock.h"
33
34 #include <cstdlib>
35 #include <time.h>
36
37
38 /****************************************************************************/
39 /* MACROS */
40 #define TIMEFORMAT " %12.6E units, %12.6E seconds\n"
41
42
43 /****************************************************************************/
44 /* EXTERNAL STUFF FROM FUNCTION MODULES*/
45
46 /*--------------------------------------------------------------------------*/
47 /* Control file name */
48 const extern char* controlFileName;
49
50 /*--------------------------------------------------------------------------*/
51 /* Dimensions */
52 extern int indepDim;
53
54 /*--------------------------------------------------------------------------*/
55 /* Init Problem Parameters */
56 extern void initProblemParameters( void );
57
58 /*--------------------------------------------------------------------------*/
59 /* Initialize indeps */
60 extern void initIndependents( double* indeps );
61
62 /*--------------------------------------------------------------------------*/
63 /* Original scalar function */
64 extern double originalScalarFunction( double* indeps );
65
66 /*--------------------------------------------------------------------------*/
67 /* Taping scalar function */
68 extern double tapingScalarFunction( int tag, double* indeps );
69
70
71 /****************************************************************************/
72 /* CONTROL STUFF */
73 enum controlParameter {
74 cpDimension,
75 cpAverageCount,
76 cpDegree,
77 cpVecCountFW,
78 cpVecCountRV,
79 cpVecCountTR,
80 cpZosFW,
81 cpFosFW,
82 cpHosFW,
83 cpFovFW,
84 cpHovFW,
85 cpFosRV,
86 cpHosRV,
87 cpFovRV,
88 cpHovRV,
89 cpFunction,
90 cpJacobian,
91 cpVecJac,
92 cpJacVec,
93 cpHessian,
94 cpHessVec,
95 cpLagHessVec,
96 cpTensor,
97 cpInvTensor,
98 cpCount
99 };
100
101
102 /****************************************************************************/
103 /* PROVIDE RANDOM INITs */
104 //unsigned short int dx[3]; /* variable needed by erand48(.) */
105
initRand(void)106 void initRand ( void ) /* a function to initialize dx using actual time */
107 { struct tm s;
108 time_t t;
109 time(&t);
110 s=*localtime(&t);
111 srand(s.tm_sec*s.tm_min);
112 /* dx[0]=rand();
113 dx[1]=rand();
114 dx[2]=rand();*/
115 }
116
117
118 /****************************************************************************/
119 /* MAIN PROGRAM */
main()120 int main() {
121 int i, j, k;
122 int tag = 1; /* tape tag */
123 int taskCount = 0;
124
125 int pFW, pRV, pTR, degree, keep; /* forward/reverse parameters */
126 int evalCount; /* # of evaluations */
127
128
129 /****************************************************************************/
130 /* READ CONTROL PARAMETERS FROM FILE */
131 int controlParameters[cpCount];
132 FILE* controlFile;
133
134 /*------------------------------------------------------------------------*/
135 /* open file to read */
136 if ((controlFile = fopen(controlFileName,"r")) == NULL) {
137 fprintf(stdout,"ERROR: Could not open control file %s\n",
138 controlFileName);
139 exit(-1);
140 }
141
142 /*------------------------------------------------------------------------*/
143 /* read all values */
144 for (i=0; i<cpCount; i++)
145 fscanf(controlFile,"%d%*[^\n]",&controlParameters[i]);
146
147 indepDim = controlParameters[cpDimension];
148 pFW = controlParameters[cpVecCountFW];
149 pRV = controlParameters[cpVecCountRV];
150 pTR = controlParameters[cpVecCountTR];
151 degree = controlParameters[cpDegree];
152 evalCount = controlParameters[cpAverageCount];
153
154 /*------------------------------------------------------------------------*/
155 /* close control file */
156 fclose(controlFile);
157
158
159 /****************************************************************************/
160 /* VARIABLES & INITIALIZATION */
161
162 /*------------------------------------------------------------------------*/
163 /* Initialize all problem parameters (including dimension) */
164 initProblemParameters();
165
166 /*------------------------------------------------------------------------*/
167 /* Initialize the independent variables */
168 double* indeps = new double[indepDim];
169 initIndependents(indeps);
170
171 /*------------------------------------------------------------------------*/
172 /* Check main parameters */
173 if (evalCount <= 0) {
174 fprintf(stdout," # of evaluations to average over = ? ");
175 fscanf(stdin,"%d",&evalCount);
176 fprintf(stdout,"\n");
177 }
178
179 if ((degree <= 1) &&
180 (controlParameters[cpHosFW] || controlParameters[cpHovFW] ||
181 controlParameters[cpHosRV] || controlParameters[cpHovRV] ||
182 controlParameters[cpTensor])) {
183 fprintf(stdout," degree = ? ");
184 fscanf(stdin,"%d",°ree);
185 fprintf(stdout,"\n");
186 }
187 keep = degree + 1;
188
189 if ((pFW < 1) &&
190 (controlParameters[cpFovFW] || controlParameters[cpHovFW])) {
191 fprintf(stdout," # of vectors in vector forward mode = ? ");
192 fscanf(stdin,"%d",&pFW);
193 fprintf(stdout,"\n");
194 }
195
196 if ((pRV < 1) &&
197 (controlParameters[cpFovRV] || controlParameters[cpHovRV])) {
198 fprintf(stdout," # of vectors in vector reverse mode = ? ");
199 fscanf(stdin,"%d",&pRV);
200 fprintf(stdout,"\n");
201 }
202
203 if ((pTR < 1) &&
204 (controlParameters[cpTensor])) {
205 fprintf(stdout," # of vectors in tensor mode = ? ");
206 fscanf(stdin,"%d",&pTR);
207 fprintf(stdout,"\n");
208 }
209
210 /*------------------------------------------------------------------------*/
211 /* Necessary variable */
212 double depOrig=0.0, depTape; /* function value */
213 double ***XPPP, **XPP;
214 double ***YPPP, **YPP, *YP;
215 double ***ZPPP, **ZPP, *ZP;
216 double *UP, u;
217 double *VP;
218 double *WP;
219 double *JP;
220 short **nzPP;
221 int retVal=0; /* return value */
222 double t00, t01, t02, t03; /* time values */
223 double **TPP;
224 double **SPP;
225 double **HPP;
226 int dim;
227
228
229 /****************************************************************************/
230 /* NORMALIZE TIMER */
231
232
233
234 /****************************************************************************/
235 /* 0. ORIGINAL FUNCTION EVALUATION */
236 /* ---> always */
237 fprintf(stdout,"\nTASK %d: Original function evaluation\n",
238 taskCount++);
239
240 t00 = myclock();
241 for (i=0; i<evalCount; i++)
242 depOrig = originalScalarFunction(indeps);
243 t01 = myclock();
244
245 double timeUnit;
246 if (t01-t00) {
247 timeUnit = 1.0/(t01-t00);
248 fprintf(stdout," ");
249 fprintf(stdout,TIMEFORMAT,1.0,
250 (t01-t00)/evalCount);
251 } else {
252 fprintf(stdout," !!! zero timing !!!\n");
253 fprintf(stdout," set time unit to 1.0\n");
254 timeUnit = 1;
255 }
256
257
258 /****************************************************************************/
259 /* 1. TAPING THE FUNCTION */
260 /* ---> always */
261 fprintf(stdout,"--------------------------------------------------------");
262 fprintf(stdout,"\nTASK %d: Taping the function\n",
263 taskCount++);
264
265 t00 = myclock();
266 /* NOTE: taping will be performed ONCE only */
267 depTape = tapingScalarFunction(tag,indeps);
268 t01 = myclock();
269
270 size_t tape_stats[STAT_SIZE];
271 tapestats(tag,tape_stats);
272
273 fprintf(stdout,"\n independents %zu\n",tape_stats[NUM_INDEPENDENTS]);
274 fprintf(stdout," dependents %zu\n",tape_stats[NUM_DEPENDENTS]);
275 fprintf(stdout," operations %zu\n",tape_stats[NUM_OPERATIONS]);
276 fprintf(stdout," operations buffer size %zu\n",tape_stats[OP_BUFFER_SIZE]);
277 fprintf(stdout," locations buffer size %zu\n",tape_stats[LOC_BUFFER_SIZE]);
278 fprintf(stdout," constants buffer size %zu\n",tape_stats[VAL_BUFFER_SIZE]);
279 fprintf(stdout," maxlive %zu\n",tape_stats[NUM_MAX_LIVES]);
280 fprintf(stdout," valstack size %zu\n\n",tape_stats[TAY_STACK_SIZE]);
281
282 fprintf(stdout," ");
283 fprintf(stdout,TIMEFORMAT,(t01-t00)*timeUnit*evalCount,
284 (t01-t00));
285
286 /****************************************************************************/
287 /* 2. ZOS_FORWARD */
288 if (controlParameters[cpZosFW]) {
289 fprintf(stdout,"--------------------------------------------------------");
290 fprintf(stdout,"\nTASK %d: forward(tag, m=1, n=%d, keep, X[n], Y[m])\n",
291 taskCount++,indepDim);
292 fprintf(stdout," ---> zos_forward\n");
293
294 /*----------------------------------------------------------------------*/
295 /* NO KEEP */
296 t00 = myclock();
297 for (i=0; i<evalCount; i++)
298 retVal = forward(tag,1,indepDim,0,indeps,&depTape);
299 t01 = myclock();
300
301 fprintf(stdout," NO KEEP");
302 fprintf(stdout,TIMEFORMAT,(t01-t00)*timeUnit,
303 (t01-t00)/evalCount);
304
305 /*----------------------------------------------------------------------*/
306 /* KEEP */
307 t02 = myclock();
308 for (i=0; i<evalCount; i++)
309 retVal = forward(tag,1,indepDim,1,indeps,&depTape);
310 t03 = myclock();
311
312 fprintf(stdout," KEEP ");
313 fprintf(stdout,TIMEFORMAT,(t03-t02)*timeUnit,
314 (t03-t02)/evalCount);
315
316 /*----------------------------------------------------------------------*/
317 /* Debug infos */
318 if (controlParameters[cpZosFW] > 1) {
319 fprintf(stdout,"\n Return value: %d\n",retVal);
320 fprintf(stdout," Should be the same values:\n");
321 fprintf(stdout," (original) %12.8E =? %12.8E (forward from tape)\n",
322 depOrig,depTape);
323 }
324 }
325
326
327 /****************************************************************************/
328 /* 3. FOS_FORWARD */
329 if (controlParameters[cpFosFW]) {
330 fprintf(stdout,"--------------------------------------------------------");
331 fprintf(stdout,"\nTASK %d: forward(tag, m=1, n=%d, d=1, keep, X[n][d+1], Y[d+1])\n",
332 taskCount++,indepDim);
333 fprintf(stdout," ---> fos_forward\n");
334
335 /*----------------------------------------------------------------------*/
336 /* Allocation & initialisation of tensors */
337 XPP = new double*[indepDim];
338 for (i=0; i<indepDim; i++) {
339 XPP[i] = new double[2];
340 XPP[i][0] = indeps[i];
341 XPP[i][1] = (double)rand();
342 }
343 YP = new double[2];
344
345 /*----------------------------------------------------------------------*/
346 /* NO KEEP */
347 t00 = myclock();
348 for (i=0; i<evalCount; i++)
349 retVal = forward(tag,1,indepDim,1,0,XPP,YP);
350 t01 = myclock();
351
352 fprintf(stdout," NO KEEP");
353 fprintf(stdout,TIMEFORMAT,(t01-t00)*timeUnit,
354 (t01-t00)/evalCount);
355
356 /*----------------------------------------------------------------------*/
357 /* KEEP */
358 t02 = myclock();
359 for (i=0; i<evalCount; i++)
360 retVal = forward(tag,1,indepDim,1,2,XPP,YP);
361 t03 = myclock();
362
363 fprintf(stdout," KEEP ");
364 fprintf(stdout,TIMEFORMAT,(t03-t02)*timeUnit,
365 (t03-t02)/evalCount);
366
367 /*----------------------------------------------------------------------*/
368 /* Debug infos */
369 if (controlParameters[cpFosFW] > 1) {
370 fprintf(stdout,"\n Return value: %d\n",retVal);
371 fprintf(stdout," Should be the same values:\n");
372 fprintf(stdout," (original) %12.8E =? %12.8E (forward from tape)\n",
373 depOrig,YP[0]);
374 }
375
376 /*----------------------------------------------------------------------*/
377 /* Free tensors */
378 for (i=0; i<indepDim; i++)
379 delete[] XPP[i];
380 delete[] XPP;
381 delete[] YP;
382 }
383
384
385 /****************************************************************************/
386 /* 4. HOS_FORWARD */
387 if (controlParameters[cpHosFW]) {
388 fprintf(stdout,"--------------------------------------------------------");
389 fprintf(stdout,"\nTASK %d: forward(tag, m=1, n=%d, d=%d, keep, X[n][d+1], Y[d+1])\n",
390 taskCount++,indepDim,degree);
391 fprintf(stdout," ---> hos_forward\n");
392
393 /*----------------------------------------------------------------------*/
394 /* Allocation & initialisation of tensors */
395 XPP = new double*[indepDim];
396 for (i=0; i<indepDim; i++) {
397 XPP[i] = new double[1+degree];
398 XPP[i][0] = indeps[i];
399 for (j=1; j<=degree; j++)
400 XPP[i][j] = (double)rand();
401 }
402 YP = new double[1+degree];
403
404 /*----------------------------------------------------------------------*/
405 /* NO KEEP */
406 t00 = myclock();
407 for (i=0; i<evalCount; i++)
408 retVal = forward(tag,1,indepDim,degree,0,XPP,YP);
409 t01 = myclock();
410
411 fprintf(stdout," NO KEEP");
412 fprintf(stdout,TIMEFORMAT,(t01-t00)*timeUnit,
413 (t01-t00)/evalCount);
414
415 /*----------------------------------------------------------------------*/
416 /* KEEP */
417 t02 = myclock();
418 for (i=0; i<evalCount; i++)
419 retVal = forward(tag,1,indepDim,degree,keep,XPP,YP);
420 t03 = myclock();
421
422 fprintf(stdout," KEEP ");
423 fprintf(stdout,TIMEFORMAT,(t03-t02)*timeUnit,
424 (t03-t02)/evalCount);
425
426 /*----------------------------------------------------------------------*/
427 /* Debug infos */
428 if (controlParameters[cpHosFW] > 1) {
429 fprintf(stdout,"\n Return value: %d\n",retVal);
430 fprintf(stdout," Should be the same values:\n");
431 fprintf(stdout," (original) %12.8E =? %12.8E (forward from tape)\n",
432 depOrig,YP[0]);
433 }
434
435 /*----------------------------------------------------------------------*/
436 /* Free tensors */
437 for (i=0; i<indepDim; i++)
438 delete[] XPP[i];
439 delete[] XPP;
440 delete[] YP;
441 }
442
443
444 /****************************************************************************/
445 /* 5. FOV_FORWARD */
446 if (controlParameters[cpFovFW]) {
447 fprintf(stdout,"--------------------------------------------------------");
448 fprintf(stdout,"\nTASK %d: forward(tag, m=1, n=%d, p=%d, x[n], X[n][p], y[m], Y[m][p])\n",
449 taskCount++,indepDim,pFW);
450 fprintf(stdout," ---> fov_forward\n");
451
452 /*----------------------------------------------------------------------*/
453 /* Allocation & initialisation of tensors */
454 XPP = new double*[indepDim];
455 for (i=0; i<indepDim; i++) {
456 XPP[i] = new double[pFW];
457 for (j=0; j<pFW; j++)
458 XPP[i][j] = (double)rand();
459 }
460 YP = new double[1];
461 YPP = new double*[1];
462 YPP[0] = new double[pFW];
463
464 /*----------------------------------------------------------------------*/
465 /* always NO KEEP */
466 t00 = myclock();
467 for (i=0; i<evalCount; i++)
468 retVal = forward(tag,1,indepDim,pFW,indeps,XPP,YP,YPP);
469 t01 = myclock();
470
471 fprintf(stdout," (NO KEEP)");
472 fprintf(stdout,TIMEFORMAT,(t01-t00)*timeUnit,
473 (t01-t00)/evalCount);
474
475 /*----------------------------------------------------------------------*/
476 /* Debug infos */
477 if (controlParameters[cpFovFW] > 1) {
478 fprintf(stdout,"\n Return value: %d\n",retVal);
479 }
480
481 /*----------------------------------------------------------------------*/
482 /* Free tensors */
483 for (i=0; i<indepDim; i++)
484 delete[] XPP[i];
485 delete[] XPP;
486 delete[] YP;
487 delete[] YPP[0];
488 delete[] YPP;
489 }
490
491
492 /****************************************************************************/
493 /* 6. HOV_FORWARD */
494 if (controlParameters[cpHovFW]) {
495 fprintf(stdout,"--------------------------------------------------------");
496 fprintf(stdout,"\nTASK %d: forward(tag, m=1, n=%d, d=%d, p=%d, x[n], X[n][p][d], y[m], Y[m][p][d])\n",
497 taskCount++,indepDim,degree,pFW);
498 fprintf(stdout," ---> hov_forward\n");
499
500 /*----------------------------------------------------------------------*/
501 /* Allocation & initialisation of tensors */
502 XPPP = new double**[indepDim];
503 for (i=0; i<indepDim; i++) {
504 XPPP[i] = new double*[pFW];
505 for (j=0; j<pFW; j++) {
506 XPPP[i][j] = new double[degree];
507 for (k=0; k<degree; k++)
508 XPPP[i][j][k] = (double)rand();
509 }
510 }
511 YP = new double[1];
512 YPPP = new double**[1];
513 YPPP[0] = new double*[pFW];
514 for (j=0; j<pFW; j++)
515 YPPP[0][j] = new double[degree];
516
517 /*----------------------------------------------------------------------*/
518 /* always NO KEEP */
519 t00 = myclock();
520 for (i=0; i<evalCount; i++)
521 retVal = forward(tag,1,indepDim,degree,pFW,indeps,XPPP,YP,YPPP);
522 t01 = myclock();
523
524 fprintf(stdout," (NO KEEP)");
525 fprintf(stdout,TIMEFORMAT,(t01-t00)*timeUnit,
526 (t01-t00)/evalCount);
527
528 /*----------------------------------------------------------------------*/
529 /* Debug infos */
530 if (controlParameters[cpHovFW] > 1) {
531 fprintf(stdout,"\n Return value: %d\n",retVal);
532 }
533
534 /*----------------------------------------------------------------------*/
535 /* Free tensors */
536 for (i=0; i<indepDim; i++) {
537 for (j=0; j<pFW; j++)
538 delete[] XPPP[i][j];
539 delete[] XPPP[i];
540 }
541 delete[] XPPP;
542 delete[] YP;
543 for (j=0; j<pFW; j++)
544 delete[] YPPP[0][j];
545 delete[] YPPP[0];
546 delete[] YPPP;
547 }
548
549
550 /****************************************************************************/
551 /* 7. FOS_REVERSE */
552 if (controlParameters[cpFosRV]) {
553 fprintf(stdout,"--------------------------------------------------------");
554 fprintf(stdout,"\nTASK %d: reverse(tag, m=1, n=%d, d=0, u, Z[n])\n",
555 taskCount++,indepDim);
556 fprintf(stdout," ---> fos_reverse\n");
557
558 /*----------------------------------------------------------------------*/
559 /* Allocation & initialisation of tensors */
560 ZP = new double[indepDim];
561 u = (double)rand();
562
563 /*----------------------------------------------------------------------*/
564 /* Forward with keep*/
565 forward(tag,1,indepDim,1,indeps,&depTape);
566
567 /*----------------------------------------------------------------------*/
568 /* Reverse */
569 t00 = myclock();
570 for (i=0; i<evalCount; i++)
571 retVal = reverse(tag,1,indepDim,0,u,ZP);
572 t01 = myclock();
573
574 fprintf(stdout," ");
575 fprintf(stdout,TIMEFORMAT,(t01-t00)*timeUnit,
576 (t01-t00)/evalCount);
577
578 /*----------------------------------------------------------------------*/
579 /* Debug infos */
580 if (controlParameters[cpFosRV] > 1) {
581 fprintf(stdout,"\n Return value: %d\n",retVal);
582 }
583
584 /*----------------------------------------------------------------------*/
585 /* Free tensors */
586 delete[] ZP;
587 }
588
589
590 /****************************************************************************/
591 /* 8. HOS_REVERSE */
592 if (controlParameters[cpHosRV]) {
593 fprintf(stdout,"--------------------------------------------------------");
594 fprintf(stdout,"\nTASK %d: reverse(tag, m=1, n=%d, d=%d, u, Z[n][d+1])\n",
595 taskCount++,indepDim,degree);
596 fprintf(stdout," ---> hos_reverse\n");
597
598 /*----------------------------------------------------------------------*/
599 /* Allocation & initialisation of tensors */
600 ZPP = new double*[indepDim];
601 for (i=0; i<indepDim; i++)
602 ZPP[i] = new double[degree+1];
603 u = (double)rand();
604 XPP = new double*[indepDim];
605 for (i=0; i<indepDim; i++) {
606 XPP[i] = new double[1+degree];
607 XPP[i][0] = indeps[i];
608 for (j=1; j<=degree; j++)
609 XPP[i][j] = (double)rand();
610 }
611 YP = new double[1+degree];
612
613 /*----------------------------------------------------------------------*/
614 /* Forward with keep*/
615 forward(tag,1,indepDim,degree,keep,XPP,YP);
616
617 /*----------------------------------------------------------------------*/
618 /* Reverse */
619 t00 = myclock();
620 for (i=0; i<evalCount; i++)
621 retVal = reverse(tag,1,indepDim,degree,u,ZPP);
622 t01 = myclock();
623
624 fprintf(stdout," ");
625 fprintf(stdout,TIMEFORMAT,(t01-t00)*timeUnit,
626 (t01-t00)/evalCount);
627
628 /*----------------------------------------------------------------------*/
629 /* Debug infos */
630 if (controlParameters[cpHosRV] > 1) {
631 fprintf(stdout,"\n Return value: %d\n",retVal);
632 }
633
634 /*----------------------------------------------------------------------*/
635 /* Free tensors */
636 for (i=0; i<indepDim; i++)
637 delete[] ZPP[i];
638 delete[] ZPP;
639 for (i=0; i<indepDim; i++)
640 delete[] XPP[i];
641 delete[] XPP;
642 delete[] YP;
643 }
644
645
646 /****************************************************************************/
647 /* 9. FOV_REVERSE */
648 if (controlParameters[cpFovRV]) {
649 fprintf(stdout,"--------------------------------------------------------");
650 fprintf(stdout,"\nTASK %d: reverse(tag, m=1, n=%d, d=0, p=%d, U[p], Z[p][n])\n",
651 taskCount++,indepDim,pRV);
652 fprintf(stdout," ---> fov_reverse\n");
653
654 /*----------------------------------------------------------------------*/
655 /* Allocation & initialisation of tensors */
656 ZPP = new double*[pRV];
657 for (i=0; i<pRV; i++)
658 ZPP[i] = new double[indepDim];
659 UP = new double[pRV];
660 for (i=0; i<pRV; i++)
661 UP[i] = (double)rand();
662
663 /*----------------------------------------------------------------------*/
664 /* Forward with keep*/
665 forward(tag,1,indepDim,1,indeps,&depTape);
666
667 /*----------------------------------------------------------------------*/
668 /* Reverse */
669 t00 = myclock();
670 for (i=0; i<evalCount; i++)
671 retVal = reverse(tag,1,indepDim,0,pRV,UP,ZPP);
672 t01 = myclock();
673
674 fprintf(stdout," ");
675 fprintf(stdout,TIMEFORMAT,(t01-t00)*timeUnit,
676 (t01-t00)/evalCount);
677
678 /*----------------------------------------------------------------------*/
679 /* Debug infos */
680 if (controlParameters[cpFovRV] > 1) {
681 fprintf(stdout,"\n Return value: %d\n",retVal);
682 }
683
684 /*----------------------------------------------------------------------*/
685 /* Free tensors */
686 for (i=0; i<pRV; i++)
687 delete[] ZPP[i];
688 delete[] ZPP;
689 delete[] UP;
690 }
691
692
693 /****************************************************************************/
694 /* 10. HOV_REVERSE */
695 if (controlParameters[cpHovRV]) {
696 fprintf(stdout,"--------------------------------------------------------");
697 fprintf(stdout,"\nTASK %d: reverse(tag, m=1, n=%d, d=%d, p=%d, U[p], Z[p][n][d+1], nz[p][n])\n",
698 taskCount++,indepDim,degree,pRV);
699 fprintf(stdout," ---> hov_reverse\n");
700
701 /*----------------------------------------------------------------------*/
702 /* Allocation & initialisation of tensors */
703 ZPPP = new double**[pRV];
704 for (i=0; i<pRV; i++) {
705 ZPPP[i] = new double*[indepDim];
706 for (j=0; j<indepDim; j++)
707 ZPPP[i][j] = new double[degree+1];
708 }
709 UP = new double[pRV];
710 for (i=0; i<pRV; i++)
711 UP[i] = (double)rand();
712 XPP = new double*[indepDim];
713 for (i=0; i<indepDim; i++) {
714 XPP[i] = new double[1+degree];
715 XPP[i][0] = indeps[i];
716 for (j=1; j<=degree; j++)
717 XPP[i][j] = (double)rand();
718 }
719 YP = new double[1+degree];
720 nzPP = new short*[pRV];
721 for (i=0; i<pRV; i++)
722 nzPP[i] = new short[indepDim];
723
724 /*----------------------------------------------------------------------*/
725 /* Forward with keep*/
726 forward(tag,1,indepDim,degree,keep,XPP,YP);
727
728 /*----------------------------------------------------------------------*/
729 /* Reverse without nonzero pattern*/
730 t00 = myclock();
731 for (i=0; i<evalCount; i++)
732 retVal = reverse(tag,1,indepDim,degree,pRV,UP,ZPPP);
733 t01 = myclock();
734
735 fprintf(stdout," ");
736 fprintf(stdout,TIMEFORMAT,(t01-t00)*timeUnit,
737 (t01-t00)/evalCount);
738
739 /*----------------------------------------------------------------------*/
740 /* Reverse with nonzero pattern*/
741 t00 = myclock();
742 for (i=0; i<evalCount; i++)
743 retVal = reverse(tag,1,indepDim,degree,pRV,UP,ZPPP,nzPP);
744 t01 = myclock();
745
746 fprintf(stdout," (NZ)");
747 fprintf(stdout,TIMEFORMAT,(t01-t00)*timeUnit,
748 (t01-t00)/evalCount);
749
750 /*----------------------------------------------------------------------*/
751 /* Debug infos */
752 if (controlParameters[cpHovRV] > 1) {
753 fprintf(stdout,"\n Return value: %d\n",retVal);
754 }
755
756 /*----------------------------------------------------------------------*/
757 /* Free tensors */
758 for (i=0; i<pRV; i++) {
759 for (j=0; j<indepDim; j++)
760 delete[] ZPPP[i][j];
761 delete[] ZPPP[i];
762 delete[] nzPP[i];
763 }
764 delete[] ZPPP;
765 delete[] nzPP;
766 delete[] UP;
767 for (i=0; i<indepDim; i++)
768 delete[] XPP[i];
769 delete[] XPP;
770 delete[] YP;
771 }
772
773
774 /****************************************************************************/
775 /* 11. FUNCTION */
776 if (controlParameters[cpFunction]) {
777 fprintf(stdout,"--------------------------------------------------------");
778 fprintf(stdout,"\nTASK %d: function(tag, m=1, n=%d, X[n], Y[m])\n",
779 taskCount++,indepDim);
780
781 /*----------------------------------------------------------------------*/
782 /* Function evaluation */
783 t00 = myclock();
784 for (i=0; i<evalCount; i++)
785 retVal = function(tag,1,indepDim,indeps,&depTape);
786 t01 = myclock();
787
788 fprintf(stdout," ");
789 fprintf(stdout,TIMEFORMAT,(t01-t00)*timeUnit,
790 (t01-t00)/evalCount);
791
792 /*----------------------------------------------------------------------*/
793 /* Debug infos */
794 if (controlParameters[cpFunction] > 1) {
795 fprintf(stdout,"\n Return value: %d\n",retVal);
796 fprintf(stdout," Should be the same values:\n");
797 fprintf(stdout," (original) %12.8E =? %12.8E (forward from tape)\n",
798 depOrig,depTape);
799 }
800 }
801
802
803 /****************************************************************************/
804 /* 12. JACOBIAN */
805 if (controlParameters[cpJacobian]) {
806 fprintf(stdout,"--------------------------------------------------------");
807 fprintf(stdout,"\nTASK %d: gradient(tag, n=%d, X[n], G[n])\n",
808 taskCount++,indepDim);
809
810 /*----------------------------------------------------------------------*/
811 /* Allocation & initialisation of tensors */
812 JP = new double[indepDim];
813
814 /*----------------------------------------------------------------------*/
815 /* Gradient evaluation */
816 t00 = myclock();
817 for (i=0; i<evalCount; i++)
818 retVal = gradient(tag,indepDim,indeps,JP);
819 t01 = myclock();
820
821 fprintf(stdout," ");
822 fprintf(stdout,TIMEFORMAT,(t01-t00)*timeUnit,
823 (t01-t00)/evalCount);
824
825 /*----------------------------------------------------------------------*/
826 /* Debug infos */
827 if (controlParameters[cpJacobian] > 1) {
828 fprintf(stdout,"\n Return value: %d\n",retVal);
829 }
830
831 /*----------------------------------------------------------------------*/
832 /* Free tensors */
833 delete[] JP;
834 }
835
836
837 /****************************************************************************/
838 /* 13. VECJAC */
839 if (controlParameters[cpVecJac]) {
840 fprintf(stdout,"--------------------------------------------------------");
841 fprintf(stdout,"\nTASK %d: vec_jac(tag, m=1, n=%d, repeat, X[n], U[m], V[n])\n",
842 taskCount++,indepDim);
843
844 /*----------------------------------------------------------------------*/
845 /* Allocation & initialisation of tensors */
846 UP = new double[1];
847 UP[0] = (double)rand();
848 VP = new double[indepDim];
849
850 /*----------------------------------------------------------------------*/
851 /* Evaluation without repeat */
852 t00 = myclock();
853 for (i=0; i<evalCount; i++)
854 retVal = vec_jac(tag,1,indepDim,0,indeps,UP,VP);
855 t01 = myclock();
856
857 fprintf(stdout,"(no repeat)");
858 fprintf(stdout,TIMEFORMAT,(t01-t00)*timeUnit,
859 (t01-t00)/evalCount);
860
861 /*----------------------------------------------------------------------*/
862 /* Evaluation with repeat */
863 t00 = myclock();
864 for (i=0; i<evalCount; i++)
865 retVal = vec_jac(tag,1,indepDim,1,indeps,UP,VP);
866 t01 = myclock();
867
868 fprintf(stdout," (repeat)");
869 fprintf(stdout,TIMEFORMAT,(t01-t00)*timeUnit,
870 (t01-t00)/evalCount);
871
872 /*----------------------------------------------------------------------*/
873 /* Debug infos */
874 if (controlParameters[cpVecJac] > 1) {
875 fprintf(stdout,"\n Return value: %d\n",retVal);
876 }
877
878 /*----------------------------------------------------------------------*/
879 /* Free tensors */
880 delete[] UP;
881 delete[] VP;
882 }
883
884
885 /****************************************************************************/
886 /* 14. JACVEC */
887 if (controlParameters[cpJacVec]) {
888 fprintf(stdout,"--------------------------------------------------------");
889 fprintf(stdout,"\nTASK %d: jac_vec(tag, m=1, n=%d, X[n], V[n], U[m])\n",
890 taskCount++,indepDim);
891
892 /*----------------------------------------------------------------------*/
893 /* Allocation & initialisation of tensors */
894 UP = new double[1];
895 VP = new double[indepDim];
896 for (i=0; i<indepDim; i++)
897 VP[i] = (double)rand();
898
899 /*----------------------------------------------------------------------*/
900 /* Evaluation */
901 t00 = myclock();
902 for (i=0; i<evalCount; i++)
903 retVal = jac_vec(tag,1,indepDim,indeps,VP,UP);
904 t01 = myclock();
905
906 fprintf(stdout," ");
907 fprintf(stdout,TIMEFORMAT,(t01-t00)*timeUnit,
908 (t01-t00)/evalCount);
909
910 /*----------------------------------------------------------------------*/
911 /* Debug infos */
912 if (controlParameters[cpJacVec] > 1) {
913 fprintf(stdout,"\n Return value: %d\n",retVal);
914 }
915
916 /*----------------------------------------------------------------------*/
917 /* Free tensors */
918 delete[] UP;
919 delete[] VP;
920 }
921
922
923 /****************************************************************************/
924 /* 15. HESSIAN */
925 if (controlParameters[cpHessian]) {
926 fprintf(stdout,"--------------------------------------------------------");
927 fprintf(stdout,"\nTASK %d: hessian(tag, n=%d, X[n], lower triangle of H[n][n])\n",
928 taskCount++,indepDim);
929
930 /*----------------------------------------------------------------------*/
931 /* Allocation & initialisation of tensors */
932 HPP = new double*[indepDim];
933 for (i=0; i<indepDim; i++)
934 HPP[i] = new double[indepDim];
935
936 /*----------------------------------------------------------------------*/
937 /* Evaluation */
938 t00 = myclock();
939 for (i=0; i<evalCount; i++)
940 retVal = hessian(tag,indepDim,indeps,HPP);
941 t01 = myclock();
942
943 fprintf(stdout," ");
944 fprintf(stdout,TIMEFORMAT,(t01-t00)*timeUnit,
945 (t01-t00)/evalCount);
946
947 /*----------------------------------------------------------------------*/
948 /* Debug infos */
949 if (controlParameters[cpHessian] > 1) {
950 fprintf(stdout,"\n Return value: %d\n",retVal);
951 }
952
953 /*----------------------------------------------------------------------*/
954 /* Free tensors */
955 for (i=0; i<indepDim; i++)
956 delete[] HPP[i];
957 delete[] HPP;
958 }
959
960
961 /****************************************************************************/
962 /* 16. HESSVEC */
963 if (controlParameters[cpHessVec]) {
964 fprintf(stdout,"--------------------------------------------------------");
965 fprintf(stdout,"\nTASK %d: hess_vec(tag, n=%d, X[n], V[n], W[n])\n",
966 taskCount++,indepDim);
967
968 /*----------------------------------------------------------------------*/
969 /* Allocation & initialisation of tensors */
970 VP = new double[indepDim];
971 for (i=0; i<indepDim; i++)
972 VP[i] = (double)rand();
973 WP = new double[indepDim];
974
975 /*----------------------------------------------------------------------*/
976 /* Evaluation */
977 t00 = myclock();
978 for (i=0; i<evalCount; i++)
979 retVal = hess_vec(tag,indepDim,indeps,VP,WP);
980 t01 = myclock();
981
982 fprintf(stdout," ");
983 fprintf(stdout,TIMEFORMAT,(t01-t00)*timeUnit,
984 (t01-t00)/evalCount);
985
986 /*----------------------------------------------------------------------*/
987 /* Debug infos */
988 if (controlParameters[cpHessVec] > 1) {
989 fprintf(stdout,"\n Return value: %d\n",retVal);
990 }
991
992 /*----------------------------------------------------------------------*/
993 /* Free tensors */
994 delete[] VP;
995 delete[] WP;
996 }
997
998
999 /****************************************************************************/
1000 /* 17. LAGHESSVEC */
1001 if (controlParameters[cpLagHessVec]) {
1002 fprintf(stdout,"--------------------------------------------------------");
1003 fprintf(stdout,"\nTASK %d: lagra_hess_vec(tag, m=1, n=%d, X[n], U[m], V[n], W[n])\n",
1004 taskCount++,indepDim);
1005
1006 /*----------------------------------------------------------------------*/
1007 /* Allocation & initialisation of tensors */
1008 UP = new double[1];
1009 UP[0] = (double)rand();
1010 VP = new double[indepDim];
1011 for (i=0; i<indepDim; i++)
1012 VP[i] = (double)rand();
1013 WP = new double[indepDim];
1014
1015 /*----------------------------------------------------------------------*/
1016 /* Evaluation */
1017 t00 = myclock();
1018 for (i=0; i<evalCount; i++)
1019 retVal = lagra_hess_vec(tag,1,indepDim,indeps,UP,VP,WP);
1020 t01 = myclock();
1021
1022 fprintf(stdout," ");
1023 fprintf(stdout,TIMEFORMAT,(t01-t00)*timeUnit,
1024 (t01-t00)/evalCount);
1025
1026 /*----------------------------------------------------------------------*/
1027 /* Debug infos */
1028 if (controlParameters[cpLagHessVec] > 1) {
1029 fprintf(stdout,"\n Return value: %d\n",retVal);
1030 }
1031
1032 /*----------------------------------------------------------------------*/
1033 /* Free tensors */
1034 delete[] VP;
1035 delete[] WP;
1036 delete[] UP;
1037 }
1038
1039
1040 /****************************************************************************/
1041 /* 18. TENSOR */
1042 if (controlParameters[cpTensor]) {
1043 fprintf(stdout,"--------------------------------------------------------");
1044 fprintf(stdout,"\nTASK %d: tensor_eval(tag, m =1, n=%d, d=%d, p=%d, X[n], tensor[m][dim], S[n][p])\n",
1045 taskCount++,indepDim,degree, pTR);
1046 fprintf(stdout,"\n dim = ((p+d) over d)\n");
1047
1048 /*----------------------------------------------------------------------*/
1049 /* Allocation & initialisation of tensors */
1050 dim = binomi(pTR+degree,degree);
1051 TPP = new double*[1];
1052 TPP[0] = new double[dim];
1053 SPP = new double*[indepDim];
1054 for (i=0; i<indepDim; i++) {
1055 SPP[i] = new double[pTR];
1056 for (j=0; j<pTR; j++)
1057 SPP[i][j]=(i==j)?1.0:0.0;
1058 }
1059
1060 /*----------------------------------------------------------------------*/
1061 /* tensor evaluation */
1062 t00 = myclock();
1063 for (i=0; i<evalCount; i++)
1064 tensor_eval(tag,1,indepDim,degree,pTR,indeps,TPP,SPP);
1065 t01 = myclock();
1066
1067 fprintf(stdout," ");
1068 fprintf(stdout,TIMEFORMAT,(t01-t00)*timeUnit,
1069 (t01-t00)/evalCount);
1070
1071 /*----------------------------------------------------------------------*/
1072 /* Debug infos */
1073 if (controlParameters[cpTensor] > 1) {}
1074
1075 /*----------------------------------------------------------------------*/
1076 /* Free tensors */
1077 delete[] TPP[0];
1078 delete[] TPP;
1079 for (i=0; i<indepDim; i++)
1080 delete[] SPP[i];
1081 delete[] SPP;
1082 }
1083
1084
1085 /****************************************************************************/
1086 /* 19. INVERSE TENSOR */
1087 if (controlParameters[cpInvTensor] && (1==indepDim)) {
1088 fprintf(stdout,"--------------------------------------------------------");
1089 fprintf(stdout,"\nTASK %d: inverse_tensor_eval(tag, m=n=1, d=%d, p=%d, X[n], tensor[m][dim], S[n][p])\n",
1090 taskCount++,degree, pTR);
1091 fprintf(stdout,"\n dim = ((p+d) over d)\n");
1092
1093 /*----------------------------------------------------------------------*/
1094 /* Allocation & initialisation of tensors */
1095 dim = binomi(pTR+degree,degree);
1096 TPP = new double*[1];
1097 TPP[0] = new double[dim];
1098 SPP = new double*[1];
1099 SPP[0] = new double[pTR];
1100 for (j=0; j<pTR; j++)
1101 SPP[0][j]=(0==j)?1.0:0.0;
1102
1103 /*----------------------------------------------------------------------*/
1104 /* tensor evaluation */
1105 t00 = myclock();
1106 for (i=0; i<evalCount; i++)
1107 inverse_tensor_eval(tag,1,degree,pTR,indeps,TPP,SPP);
1108 t01 = myclock();
1109
1110 fprintf(stdout," ");
1111 fprintf(stdout,TIMEFORMAT,(t01-t00)*timeUnit,
1112 (t01-t00)/evalCount);
1113
1114 /*----------------------------------------------------------------------*/
1115 /* Debug infos */
1116 if (controlParameters[cpInvTensor] > 1) {}
1117
1118 /*----------------------------------------------------------------------*/
1119 /* Free tensors */
1120 delete[] TPP[0];
1121 delete[] TPP;
1122 delete[] SPP[0];
1123 delete[] SPP;
1124 }
1125
1126 return 1;
1127 }
1128
1129 #undef _SGENMAIN_C_
1130
1131
1132