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