1 /*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File: revolve.c
4 Revision: $Id: revolve.c 704 2016-05-17 11:49:55Z kulshres $
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 occurred *
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 <adolc/revolve.h>
187 #include "taping_p.h"
188
189 #define MAXINT 2147483647
190
191 #ifndef _OPENMP
192 revolve_nums revolve_numbers;
193 #else
194 revolve_nums *revolve_numbers = NULL;
195 #endif
196
197 /* ************************************************************************* */
198
numforw(int steps,int snaps)199 int 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
expense(int steps,int snaps)229 double 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
maxrange(int ss,int tt)249 int 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
adjust(int steps)273 int 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
revolve(int * check,int * capo,int * fine,int snaps,int * info)304 enum 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
498