source: trunk/ADOL-C/examples/additional_examples/timing/sgenmain.cpp @ 171

Last change on this file since 171 was 171, checked in by kulshres, 9 years ago

Squashed merge branch 'master' of 'gitclone' into svn

  • 'master' of 'gitclone': (84 commits) adjust example makefiles and include paths get rid of the symlink in the src subdirectory

details of the commits:
commit c9e4bc332d2363f737fc2e8a8fcfc2e43ddb9d15
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Mon Oct 4 15:43:47 2010 +0200

adjust example makefiles and include paths

include paths in example sources were wrong for some time now
simplify makefile rules too, there is really no need for checking SPARSE
adjust include paths in makefiles.

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

commit e6e1963e41e097fd5b4a79cd1611c12f6868dc94
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Mon Oct 4 15:41:25 2010 +0200

get rid of the symlink in the src subdirectory

windows doesn't like symlinks and make infinite depth directories
we now create a symlink for build in the directory parallel to src
adjust all makefiles.am accordingly for build

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

  • Property svn:keywords set to Author Date Id Revision
File size: 42.2 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     sgenmain.cpp
4 Revision: $Id: sgenmain.cpp 171 2010-10-04 13:57:19Z 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 */
48const extern char* controlFileName;
49
50/*--------------------------------------------------------------------------*/
51/*                                                               Dimensions */
52extern int indepDim;
53
54/*--------------------------------------------------------------------------*/
55/*                                                  Init Problem Parameters */
56extern void initProblemParameters( void );
57
58/*--------------------------------------------------------------------------*/
59/*                                                        Initialize indeps */
60extern void initIndependents( double* indeps );
61
62/*--------------------------------------------------------------------------*/
63/*                                                 Original scalar function */
64extern double originalScalarFunction( double* indeps );
65
66/*--------------------------------------------------------------------------*/
67/*                                                   Taping scalar function */
68extern double tapingScalarFunction( int tag, double* indeps );
69
70
71/****************************************************************************/
72/*                                                            CONTROL STUFF */
73enum 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
106void 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 */
120int 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",&degree);
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    int tape_stats[STAT_SIZE];
271    tapestats(tag,tape_stats);
272
273    fprintf(stdout,"\n    independents            %d\n",tape_stats[NUM_INDEPENDENTS]);
274    fprintf(stdout,"    dependents              %d\n",tape_stats[NUM_DEPENDENTS]);
275    fprintf(stdout,"    operations              %d\n",tape_stats[NUM_OPERATIONS]);
276    fprintf(stdout,"    operations buffer size  %d\n",tape_stats[OP_BUFFER_SIZE]);
277    fprintf(stdout,"    locations buffer size   %d\n",tape_stats[LOC_BUFFER_SIZE]);
278    fprintf(stdout,"    constants buffer size   %d\n",tape_stats[VAL_BUFFER_SIZE]);
279    fprintf(stdout,"    maxlive                 %d\n",tape_stats[NUM_MAX_LIVES]);
280    fprintf(stdout,"    valstack size           %d\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
Note: See TracBrowser for help on using the repository browser.