source: trunk/ADOL-C/src/revolve.c @ 71

Last change on this file since 71 was 42, checked in by awalther, 10 years ago

set svn keywords property

  • Property svn:keywords set to Author Date Id Revision
File size: 22.5 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     revolve.c
4 Revision: $Id: revolve.c 42 2009-07-15 18:37:17Z awalther $
5 Contents: optimal binomial checkpointing adapted for ADOL-C
6
7 Copyright (c) Andrea Walther, Andreas Griewank, Andreas Kowarz
8 
9 This file is part of ADOL-C. This software is provided as open source.
10 Any use, reproduction, or distribution of the software constitutes
11 recipient's acceptance of the terms of the accompanying license file.
12 
13---------------------------------------------------------------------------*/
14
15/* -----
16*   The function REVOLVE coded below is meant to be used as a        *
17*   "controller" for running a time-dependent applications program   *
18*   in the reverse mode with checkpointing described in the paper    *
19*   "Achieving logarithmic Growth in temporal and spatial complexity *
20*   in reverse automatic differentiation", Optimization Methods and  *
21*   Software,  Vol.1 pp. 35-54.                                      *
22*   A postscript source of that paper can be found in the ftp sites  *
23*        info.mcs.anl.gov and nbtf02.math.tu-dresden.de.             *
24*   Apart from REVOLVE this file contains five auxiliary routines    *
25*   NUMFORW, EXPENSE, MAXRANGE, and ADJUST.                          *
26*                                                                    *
27*--------------------------------------------------------------------*
28*                                                                    *
29*   To utilize REVOLVE the user must have procedures for             *
30*     - Advancing the state of the modeled system to a certain time. *
31*     - Saving the current state onto a stack of snapshots.          *
32*     - Restoring the the most recently saved snapshot and           *
33*       restarting the forward simulation from there.                *
34*     - Initializing the adjoints at the end of forward sweep.       *
35*     - Performing one combined forward and adjoint step.            *
36*   Through an encoding of its return value REVOLVE asks the         *
37*   calling program to perform one of these 'actions', which we will *
38*   refer to as                                                      *
39*                                                                    *
40*       'advance', 'takeshot', 'restore', 'firsturn' and 'youturn'  .*
41*   There are two other return values, namely                        *
42*       'terminate'   and     'error'                                *
43*   which indicate a regular or faulty termination of the calls      *
44*   to REVOLVE.                                                      *
45*                                                                    *
46*   The action 'firsturn' includes a 'youturn', in that it requires  *
47*     -advancing through the last time-step with recording           *
48*      of intermediates                                              *
49*     -initializing the adjoint values (possibly after               *
50*      performing some IO)                                           *
51*     -reversing the last time step using the record just written    *
52*   The action 'firsturn' is obtained when the difference FINE-CAPO  *
53*   has been reduced to 1 for the first time.                        *
54*                                                                    *
55*--------------------------------------------------------------------*
56*                                                                    *
57*   The calling sequence is                                          *
58*                                                                    *
59*               REVOLVE(CHECK,CAPO,FINE,SNAPS,INFO)                  *
60*                                                                    *
61*   with the return value being one of the actions to be taken. The  *
62*   calling parameters are all integers with the following meaning   *
63*                                                                    *
64*         CHECK     number of checkpoint being written or retrieved  *
65*         CAPO      beginning of subrange currently being processed  *
66*         FINE      end of subrange currently being processed        *
67*         SNAPS     upper bound on number of checkpoints taken       *
68*         INFO      determines how much information will be printed  *
69*                   and contains information about an error occured  *
70*                                                                    *
71*   Since REVOLVE involves only a few integer operations its         *
72*   run-time is truly negligible within any nontrivial application.  *
73*                                                                    *
74*   The parameter SNAPS is selected by the user (possibly with the   *
75*   help of the routines EXPENSE and ADJUST described below ) and    *
76*   remains unchanged throughout.                                    *
77*                                                                    *
78*   The pair (CAPO,FINE) always represents the initial and final     *
79*   state of the subsequence of time steps currently being traversed *
80*   backwards.                                                       *
81*                                                                    *
82*   The conditions                                                   *
83*                    CHECK >= -1      and     CAPO <= FINE           *
84*   are necessary and sufficient for a regular response of REVOLVE.  *
85*   If either condition is violated the value 'error' is returned.   *
86*                                                                    *
87*   The first call to REVOLVE must be with CHECK=-1 so that          *
88*   appropriate initializations can be performed internally.         *
89*                                                                    *
90*   When CHECK =-1 and CAPO = FINE  then 'terminate' is returned as  *
91*   action value. This combination necessarily arises after a        *
92*   sufficiently large number of calls to REVOLVE, which depends     *
93*   only on the initial difference FINE-CAPO.                        *
94*                                                                    *
95*   The last parameter INFO determines how much information about    *
96*   the actions performed will be printed. When INFO =0 no           *
97*   information is sent to standard output. When INFO > 0 REVOLVE    *
98*   produces an output that contains a prediction of the number of   *   
99*   forward steps and of the factor by which the execution will slow *   
100*   down. When an error occurs, the return value of INFO contains    *
101*   information about the reason:                                    *
102*                                                                    *
103*     INFO = 10: number of checkpoints stored exceeds CHECKUP,       *
104*                increase constant CHECKUP and recompile             *
105*     INFO = 11: number of checkpoints stored exceeds SNAPS, ensure  *
106*                SNAPS greater than 0 and increase initial FINE      *
107*     INFO = 12: error occurs in NUMFORW                             *
108*     INFO = 13: enhancement of FINE, SNAPS checkpoints stored,      *
109*                SNAPS must be increased                             *
110*     INFO = 14: number of SNAPS exceeds CHECKUP, increase constant  *
111*                CHECKUP and recompile                               *
112*     INFO = 15: number of REPS exceeds REPSUP, increase constant    *
113*                REPSUP and recompile                                *
114*                                                                    *
115*--------------------------------------------------------------------*
116*                                                                    *
117*   Some further explanations and motivations:                       *
118*                                                                    *
119*   There is an implicit bound on CHECK through the dimensioning of  *
120*   the integer array CH[CHEKUP] with CHECKUP = 64 being the default.*
121*   If anybody wants to have that even larger he must change the     *
122*   source. Also for the variable REPS an upper bound REPSUP is      *
123*   defined. The default value equals 64. If during a call to        *
124*   TREEVERSE a (CHECKUP+1)-st checkpoint would normally be called   *
125*   for then control is returned after an appropriate error message. *
126*   When the calculated REPS exceeds REPSUP also an error message    *
127*   occurs.                                                          *
128*   During the forward sweep the user is free to change the last     *
129*   three parameters from call to call, except that FINE may never   *
130*   be less than the current value of CAPO. This may be useful when  *
131*   the total number of time STEPS to be taken is not a priori       *
132*   known. The choice FINE=CAPO+1 initiates the reverse sweep, which *
133*   happens automatically if is left constant as CAPO is eventually  *
134*   moved up to it. Once the first reverse or restore action has     *
135*   been taken only the last two parameters should be changed.       *
136*                                                                    *
137*--------------------------------------------------------------------*
138*                                                                    *
139*   The necessary number of forward steps without recording is       *
140*   calculated by the function                                       *
141*                                                                    *
142*                      NUMFORW(STEPS,SNAPS)                          *
143*                                                                    *
144*   STEPS denotes the total number of time steps, i.e. FINE-CAPO     *
145*   during the first call of REVOLVE. When SNAPS is less than 1 an   *
146*   error message will be given and -1 is returned as value.         *
147*                                                                    *
148*--------------------------------------------------------------------*
149*                                                                    *
150*   To choose an appropriated value of SNAPS the function            *
151*                                                                    *
152*                      EXPENSE(STEPS,SNAPS)                          *
153*                                                                    *
154*   estimates the run-time factor incurred by REVOLVE for a          *
155*   particular value of SNAPS. The ratio NUMFORW(STEPS,SNAPS)/STEPS  *
156*   is returned. This ratio corresponds to the run-time factor of    *
157*   the execution relative to the run-time of one forward time step. *
158*                                                                    *
159*--------------------------------------------------------------------*
160*                                                                    *
161*   The auxiliary function                                           *
162*                                                                    *
163*                      MAXRANGE(SNAPS,REPS)                          *
164*                                                                    *
165*   returns the integer (SNAPS+REPS)!/(SNAPS!REPS!) provided         *
166*   SNAPS >=0, REPS >= 0. Otherwise there will be appropriate error  *
167*   messages and the value -1 will be returned. If the binomial      *
168*   expression is not representable as a  signed 4 byte integer,     *
169*   greater than 2^31-1, this maximal value is returned and a        *
170*   warning message printed.                                         *
171*                                                                    *
172*--------------------------------------------------------------------*
173*                                                                    *
174*   Furthermore, the function                                        *
175*                                                                    *
176*                      ADJUST(STEPS)                                 *
177*                                                                    *
178*   is provided. It can be used to determine a value of SNAPS so     *
179*   that the increase in spatial complexity equals approximately the *
180*   increase in temporal complexity. For that ADJUST computes a      *
181*   return value satisfying SNAPS ~= log_4 (STEPS) because of the    *
182*   theory developed in the paper mentioned above.                   *
183*                                                                    *
184*--------------------------------------------------------------------*/
185
186#include <revolve.h>
187#include <taping_p.h>
188
189#define MAXINT 2147483647
190
191#ifndef _OPENMP
192revolve_nums revolve_numbers;
193#else
194revolve_nums *revolve_numbers = NULL;
195#endif
196
197/* ************************************************************************* */
198
199int numforw(int steps, int snaps) {
200    int reps, range, num;
201
202    if (snaps < 1) {
203        printf(" error occurs in numforw: snaps < 1\n");
204        return -1;
205    }
206    if (snaps > ADOLC_CHECKUP) {
207        printf(" number of snaps=%d exceeds ADOLC_CHECKUP \n",snaps);
208        printf(" redefine 'ADOLC_CHECKUP' \n");
209        return -1;
210    }
211    reps = 0;
212    range = 1;
213    while(range < steps) {
214        reps += 1;
215        range = range*(reps + snaps)/reps;
216    }
217    printf("range =  %d \n",range);
218    if (reps > ADOLC_REPSUP) {
219        printf(" number of reps=%d exceeds ADOLC_REPSUP \n",reps);
220        printf(" redefine 'ADOLC_REPSUP' \n");
221        return -1;
222    }
223    num = reps * steps - range*reps/(snaps+1);
224    return num;
225}
226
227/* ************************************************************************* */
228
229double expense(int steps, int snaps) {
230    double ratio;
231
232    if (snaps < 1) {
233        printf(" error occurs in expense: snaps < 0\n");
234        return -1;
235    }
236    if (steps < 1) {
237        printf(" error occurs in expense: steps < 0\n");
238        return -1;
239    }
240    ratio = ((double) numforw(steps,snaps));
241    if (ratio == -1)
242        return -1;
243    ratio = ratio/steps;
244    return ratio;
245}
246
247/* ************************************************************************* */
248
249int maxrange(int ss, int tt) {
250    int i, ires;
251    double res = 1.0;
252
253    if((tt<0) || (ss<0)) {
254        printf("error in MAXRANGE: negative parameter");
255        return -1;
256    }
257    for(i=1; i<= tt; i++) {
258        res *= (ss + i);
259        res /= i;
260        if (res > MAXINT) {
261            ires=MAXINT;
262            printf("warning from MAXRANGE: returned maximal integer %d\n",
263                    ires);
264            return ires;
265        }
266    }
267    ires = res;
268    return ires;
269}
270
271/* ************************************************************************* */
272
273int adjust(int steps) {
274    int snaps, s, reps;
275
276    snaps = 1;
277    reps = 1;
278    s = 0;
279    while( maxrange(snaps+s, reps+s) > steps )
280        s--;
281    while( maxrange(snaps+s, reps+s) < steps )
282        s++;
283    snaps += s;
284    reps += s ;
285    s = -1;
286    while( maxrange(snaps,reps) >= steps ) {
287        if (snaps > reps) {
288            snaps -= 1;
289            s = 0;
290        } else {
291            reps -= 1;
292            s = 1;
293        }
294    }
295    if ( s == 0 )
296        snaps += 1 ;
297    if ( s == 1 )
298        reps += 1;
299    return snaps;
300}
301
302/* ************************************************************************* */
303
304enum revolve_action revolve
305(int* check,int* capo,int* fine,int snaps,int* info) {
306    int ds, oldcapo, num, bino1, bino2, bino3, bino4, bino5, bino6;
307    /* (*capo,*fine) is the time range currently under consideration */
308    /* ch[j] is the number of the state that is stored in checkpoint j */
309    ADOLC_OPENMP_THREAD_NUMBER;
310
311    ADOLC_OPENMP_GET_THREAD_NUMBER;
312    REVOLVE_NUMBERS.commands += 1;
313    if ((*check < -1) || (*capo > *fine)) {
314        *info = 9;
315        return revolve_error;
316    }
317    if ((*check == -1) && (*capo < *fine)) {
318        if (*check == -1)
319            REVOLVE_NUMBERS.turn = 0;   /* initialization of turn counter */
320        *REVOLVE_NUMBERS.ch = *capo-1;
321    }
322    switch(*fine-*capo) {
323        case 0:   /* reduce capo to previous checkpoint, unless done  */
324            if(*check == -1 || *capo==*REVOLVE_NUMBERS.ch ) {
325                *check -= 1;
326                if (*info > 0) {
327                    printf(" \n advances: %5d",REVOLVE_NUMBERS.advances);
328                    printf(" \n takeshots: %4d",REVOLVE_NUMBERS.takeshots);
329                    printf(" \n commands: %5d \n",REVOLVE_NUMBERS.commands);
330                }
331                return revolve_terminate;
332            } else {
333                *capo = REVOLVE_NUMBERS.ch[*check];
334                REVOLVE_NUMBERS.oldfine = *fine;
335                return revolve_restore;
336            }
337        case 1/* (possibly first) combined forward/reverse step */
338            *fine -= 1;
339            if(*check >= 0 && REVOLVE_NUMBERS.ch[*check] == *capo)
340                *check -= 1;
341            if(REVOLVE_NUMBERS.turn == 0) {
342                REVOLVE_NUMBERS.turn = 1;
343                REVOLVE_NUMBERS.oldfine = *fine;
344                return revolve_firsturn;
345            } else {
346                REVOLVE_NUMBERS.oldfine = *fine;
347                return revolve_youturn;
348            }
349        default:
350            if(*check == -1 || REVOLVE_NUMBERS.ch[*check] != *capo) {
351                *check += 1 ;
352                if(*check >= ADOLC_CHECKUP) {
353                    *info = 10;
354                    return revolve_error;
355                }
356                if(*check+1 > snaps) {
357                    *info = 11;
358                    return revolve_error;
359                }
360                REVOLVE_NUMBERS.ch[*check] = *capo;
361                if (*check == 0) {
362                    REVOLVE_NUMBERS.advances = 0;
363                    REVOLVE_NUMBERS.takeshots = 0;
364                    REVOLVE_NUMBERS.commands = 1;
365                    REVOLVE_NUMBERS.oldsnaps = snaps;
366                    if (snaps > ADOLC_CHECKUP) {
367                        *info = 14;
368                        return revolve_error;
369                    }
370                    if (*info > 0) {
371                        num = numforw(*fine-*capo,snaps);
372                        if (num == -1) {
373                            *info = 12;
374                            return revolve_error;
375                        }
376                        printf(" prediction of needed forward steps: %8d => "
377                                "\n",num);
378                        printf(" slowdown factor: %8.4f \n\n",
379                                ((double) num)/(*fine-*capo));
380                    }
381                }
382                REVOLVE_NUMBERS.takeshots += 1;
383                REVOLVE_NUMBERS.oldfine = *fine;
384                return revolve_takeshot;
385            } else {
386                if ((REVOLVE_NUMBERS.oldfine < *fine) &&
387                        (snaps == *check+1))
388                {
389                    *info = 13;
390                    return revolve_error;
391                }
392                oldcapo = *capo;
393                ds = snaps - *check;
394                if (ds < 1) {
395                    *info = 11;
396                    return revolve_error;
397                }
398                REVOLVE_NUMBERS.reps = 0;
399                REVOLVE_NUMBERS.range = 1;
400                while(REVOLVE_NUMBERS.range < *fine - *capo) {
401                    REVOLVE_NUMBERS.reps += 1;
402                    REVOLVE_NUMBERS.range = REVOLVE_NUMBERS.range *
403                        (REVOLVE_NUMBERS.reps + ds) / REVOLVE_NUMBERS.reps;
404                }
405                if (REVOLVE_NUMBERS.reps > ADOLC_REPSUP) {
406                    *info = 15;
407                    return revolve_error;
408                }
409                if (snaps != REVOLVE_NUMBERS.oldsnaps) {
410                    if (snaps > ADOLC_CHECKUP) {
411                        *info = 14;
412                        return revolve_error;
413                    }
414                }
415
416                bino1 = REVOLVE_NUMBERS.range * REVOLVE_NUMBERS.reps /
417                    (ds+REVOLVE_NUMBERS.reps);
418                bino2 = (ds > 1) ? bino1*ds/(ds+REVOLVE_NUMBERS.reps-1) : 1;
419                if (ds == 1)
420                    bino3 = 0;
421                else
422                    bino3 = (ds > 2) ? bino2 * (ds - 1) /
423                        (ds + REVOLVE_NUMBERS.reps - 2) : 1;
424                bino4 = bino2*(REVOLVE_NUMBERS.reps-1)/ds;
425                if (ds < 3)
426                    bino5 = 0;
427                else
428                    bino5 = (ds > 3) ? bino3*(ds-2)/REVOLVE_NUMBERS.reps : 1;
429
430                bino6 = 0;
431
432                /* range = beta(c,r) >= l (r -> min)
433                 * bino1 = beta(c,r-1)
434                 * bino2 = beta(c-1,r-1)
435                 * bino3 = beta(c-2,r-1)
436                 * bino4 = beta(c,r-2)
437                 * bino5 = beta(c-3,r) */
438
439                /* new version by A. Kowarz
440                 * l^ as large as possible
441                 *         bino6 = beta(c-1,r-2)
442                 
443                        if (ds < 1)
444                           bino6 = 0;
445                        else
446                           bino6 = (ds > 1) ? bino2*(reps-1)/(ds+reps-2) : 1;
447                 
448                        if (*fine-*capo>=range-bino5)
449                           *capo += bino1;
450                        else
451                           if (*fine-*capo>bino1+bino2)
452                              *capo = *fine-bino2-bino3;
453                           else
454                              if (*fine-*capo>=bino1+bino6)
455                                 *capo += bino1-bino3;
456                              else
457                                 *capo = *fine-bino1+bino4; */
458
459                /* new version by A. Kowarz
460                 * l^ as small as possible
461                 *         bino6 = beta(c-1,r) */
462
463                bino6 = bino1*ds/REVOLVE_NUMBERS.reps;
464
465                if (*fine-*capo<=bino1+bino3)
466                    *capo += bino4;
467                else
468                    if (*fine-*capo<bino1+bino2)
469                        *capo = *fine-bino2-bino3;
470                    else
471                        if (*fine-*capo<=bino1+bino2+bino5)
472                            *capo += bino1-bino3;
473                        else
474                            *capo = *fine-bino6;
475
476                /* original by A. Walther
477                 
478                        if (*fine-*capo <= bino1 + bino3)
479                          *capo = *capo+bino4;
480                        else
481                         {
482                          if (*fine-*capo >= range - bino5)
483                            *capo = *capo + bino1;
484                          else
485                             *capo = *fine-bino2-bino3;
486                         } */
487
488                if (*capo == oldcapo)
489                    *capo = oldcapo+1;
490                REVOLVE_NUMBERS.advances = REVOLVE_NUMBERS.advances +
491                    *capo - oldcapo;
492                REVOLVE_NUMBERS.oldfine = *fine;
493                return revolve_advance;
494            }
495    }
496}
497
Note: See TracBrowser for help on using the repository browser.