C ALGORITHM 799, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 26, NO. 1, March, 2000, P. 19-- 45. #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # C/ # C/Drivers/ # C/Drivers/data # C/Drivers/driver.c # C/Drivers/res # C/Src/ # C/Src/revolve.h # C/Src/src.c # Doc/ # Doc/makefile_c # Doc/makefile_f # Fortran77/ # Fortran77/Drivers/ # Fortran77/Drivers/data # Fortran77/Drivers/driver.f # Fortran77/Drivers/res # Fortran77/Src/ # Fortran77/Src/src.f # This archive created: Mon Sep 4 12:16:34 2000 export PATH; PATH=/bin:$PATH if test ! -d 'C' then mkdir 'C' fi cd 'C' if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test -f 'data' then echo shar: will not over-write existing file "'data'" else cat << SHAR_EOF > 'data' 14 3 3 SHAR_EOF fi # end of overwriting check if test -f 'driver.c' then echo shar: will not over-write existing file "'driver.c'" else cat << SHAR_EOF > 'driver.c' #include #include #include #include #include "revolve.h" #include #include int main() { int check, capo, fine, steps, snaps; int info; enum action whatodo; capo = 0; snaps = 0; printf(" ENTER: STEPS, SNAPS, INFO \n"); scanf("%i",&steps); scanf("%i",&snaps); scanf("%i",&info); printf(" %d %d %d \n",steps,snaps,info); fine = steps + capo; check = -1; /* Neccessary for first call */ do { whatodo = revolve(&check, &capo, &fine, snaps, &info); if ((whatodo == takeshot) && (info > 1)) printf(" takeshot at %6d \n",capo); if ((whatodo == advance) && (info > 2)) printf(" advance to %7d \n",capo); if ((whatodo == firsturn) && (info > 2)) printf(" firsturn at %6d \n",capo); if ((whatodo == youturn) && (info > 2)) printf(" youturn at %7d \n",capo); if ((whatodo == restore) && (info > 2)) printf(" restore at %7d \n",capo); if (whatodo == error) { printf(" irregular termination of revolve \n"); switch(info) { case 10: printf(" number of checkpoints stored exceeds checkup, \n"); printf(" increase constant 'checkup' and recompile \n"); break; case 11: printf(" number of checkpoints stored = %d exceeds snaps = %d, \n" ,check+1,snaps); printf(" ensure 'snaps' > 0 and increase initial 'fine' \n"); break; case 12: printf(" error occurs in numforw \n"); break; case 13: printf(" enhancement of 'fine', 'snaps' checkpoints stored, \n"); printf(" increase 'snaps'\n"); break; case 14: printf(" number of snaps exceeds snapsup, "); printf(" increase constant 'snapsup' and recompile \n"); break; case 15: printf(" number of reps exceeds repsup, "); printf(" increase constant 'repsup' and recompile \n"); } } } while((whatodo != terminate) && (whatodo != error)); return info; } SHAR_EOF fi # end of overwriting check if test -f 'res' then echo shar: will not over-write existing file "'res'" else cat << SHAR_EOF > 'res' ENTER: STEPS, SNAPS, INFO 14 3 3 prediction of needed forward steps: 27 => slowdown factor: 1.9286 takeshot at 0 advance to 5 takeshot at 5 advance to 10 takeshot at 10 advance to 13 firsturn at 13 restore at 10 advance to 12 youturn at 12 restore at 10 advance to 11 youturn at 11 restore at 10 youturn at 10 restore at 5 advance to 7 takeshot at 7 advance to 9 youturn at 9 restore at 7 advance to 8 youturn at 8 restore at 7 youturn at 7 restore at 5 advance to 6 youturn at 6 restore at 5 youturn at 5 restore at 0 advance to 1 takeshot at 1 advance to 2 takeshot at 2 advance to 4 youturn at 4 restore at 2 advance to 3 youturn at 3 restore at 2 youturn at 2 restore at 1 youturn at 1 restore at 0 youturn at 0 advances: 27 takeshots: 6 commands: 47 SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test -f 'revolve.h' then echo shar: will not over-write existing file "'revolve.h'" else cat << SHAR_EOF > 'revolve.h' #ifndef _REVOLVE_H_ #define _REVOLVE_H_ enum action { advance, takeshot, restore, firsturn, youturn, terminate, error}; int maxrange(int ss, int tt); int adjustsize(int* steps, int* snaps, int* reps); enum action revolve(int* check,int* capo,int* fine,int snaps,int* info); #endif SHAR_EOF fi # end of overwriting check if test -f 'src.c' then echo shar: will not over-write existing file "'src.c'" else cat << SHAR_EOF > 'src.c' /* ----- * The function REVOLVE coded below is meant to be used as a * * "controller" for running a time-dependent applications program * * in the reverse mode with checkpointing described in the paper * * "Achieving logarithmic Growth in temporal and spatial complexity * * in reverse automatic differentiation", Optimization Methods and * * Software, Vol.1 pp. 35-54. * * A postscript source of that paper can be found in the ftp sites * * info.mcs.anl.gov and nbtf02.math.tu-dresden.de. * * Apart from REVOLVE this file contains five auxiliary routines * * NUMFORW, EXPENSE, MAXRANGE, and ADJUST. * * * *--------------------------------------------------------------------* * * * To utilize REVOLVE the user must have procedures for * * - Advancing the state of the modeled system to a certain time. * * - Saving the current state onto a stack of snapshots. * * - Restoring the the most recently saved snapshot and * * restarting the forward simulation from there. * * - Initializing the adjoints at the end of forward sweep. * * - Performing one combined forward and adjoint step. * * Through an encoding of its return value REVOLVE asks the * * calling program to perform one of these 'actions', which we will * * refer to as * * * * 'advance', 'takeshot', 'restore', 'firsturn' and 'youturn' .* * There are two other return values, namely * * 'terminate' and 'error' * * which indicate a regular or faulty termination of the calls * * to REVOLVE. * * * * The action 'firsturn' includes a 'youturn', in that it requires * * -advancing through the last time-step with recording * * of intermediates * * -initializing the adjoint values (possibly after * * performing some IO) * * -reversing the last time step using the record just written * * The action 'firsturn' is obtained when the difference FINE-CAPO * * has been reduced to 1 for the first time. * * * *--------------------------------------------------------------------* * * * The calling sequence is * * * * REVOLVE(CHECK,CAPO,FINE,SNAPS,INFO) * * * * with the return value being one of the actions to be taken. The * * calling parameters are all integers with the following meaning * * * * CHECK number of checkpoint being written or retrieved * * CAPO beginning of subrange currently being processed * * FINE end of subrange currently being processed * * SNAPS upper bound on number of checkpoints taken * * INFO determines how much information will be printed * * and contains information about an error occured * * * * Since REVOLVE involves only a few integer operations its * * run-time is truly negligible within any nontrivial application. * * * * The parameter SNAPS is selected by the user (possibly with the * * help of the routines EXPENSE and ADJUST described below ) and * * remains unchanged throughout. * * * * The pair (CAPO,FINE) always represents the initial and final * * state of the subsequence of time steps currently being traversed * * backwards. * * * * The conditions * * CHECK >= -1 and CAPO <= FINE * * are necessary and sufficient for a regular response of REVOLVE. * * If either condition is violated the value 'error' is returned. * * * * The first call to REVOLVE must be with CHECK=-1 so that * * appropriate initializations can be performed internally. * * * * When CHECK =-1 and CAPO = FINE then 'terminate' is returned as * * action value. This combination necessarily arises after a * * sufficiently large number of calls to REVOLVE, which depends * * only on the initial difference FINE-CAPO. * * * * The last parameter INFO determines how much information about * * the actions performed will be printed. When INFO =0 no * * information is sent to standard output. When INFO > 0 REVOLVE * * produces an output that contains a prediction of the number of * * forward steps and of the factor by which the execution will slow * * down. When an error occurs, the return value of INFO contains * * information about the reason: * * * * INFO = 10: number of checkpoints stored exceeds CHECKUP, * * increase constant CHECKUP and recompile * * INFO = 11: number of checkpoints stored exceeds SNAPS, ensure * * SNAPS greater than 0 and increase initial FINE * * INFO = 12: error occurs in NUMFORW * * INFO = 13: enhancement of FINE, SNAPS checkpoints stored, * * SNAPS must be increased * * INFO = 14: number of SNAPS exceeds CHECKUP, increase constant * * CHECKUP and recompile * * INFO = 15: number of REPS exceeds REPSUP, increase constant * * REPSUP and recompile * * * *--------------------------------------------------------------------* * * * Some further explanations and motivations: * * * * There is an implicit bound on CHECK through the dimensioning of * * the integer array CH[CHEKUP] with CHECKUP = 64 being the default.* * If anybody wants to have that even larger he must change the * * source. Also for the variable REPS an upper bound REPSUP is * * defined. The default value equals 64. If during a call to * * TREEVERSE a (CHECKUP+1)-st checkpoint would normally be called * * for then control is returned after an appropriate error message. * * When the calculated REPS exceeds REPSUP also an error message * * occurs. * * During the forward sweep the user is free to change the last * * three parameters from call to call, except that FINE may never * * be less than the current value of CAPO. This may be useful when * * the total number of time STEPS to be taken is not a priori * * known. The choice FINE=CAPO+1 initiates the reverse sweep, which * * happens automatically if is left constant as CAPO is eventually * * moved up to it. Once the first reverse or restore action has * * been taken only the last two parameters should be changed. * * * *--------------------------------------------------------------------* * * * The necessary number of forward steps without recording is * * calculated by the function * * * * NUMFORW(STEPS,SNAPS) * * * * STEPS denotes the total number of time steps, i.e. FINE-CAPO * * during the first call of REVOLVE. When SNAPS is less than 1 an * * error message will be given and -1 is returned as value. * * * *--------------------------------------------------------------------* * * * To choose an appropriated value of SNAPS the function * * * * EXPENSE(STEPS,SNAPS) * * * * estimates the run-time factor incurred by REVOLVE for a * * particular value of SNAPS. The ratio NUMFORW(STEPS,SNAPS)/STEPS * * is returned. This ratio corresponds to the run-time factor of * * the execution relative to the run-time of one forward time step. * * * *--------------------------------------------------------------------* * * * The auxiliary function * * * * MAXRANGE(SNAPS,REPS) * * * * returns the integer (SNAPS+REPS)!/(SNAPS!REPS!) provided * * SNAPS >=0, REPS >= 0. Otherwise there will be appropriate error * * messages and the value -1 will be returned. If the binomial * * expression is not representable as a signed 4 byte integer, * * greater than 2^31-1, this maximal value is returned and a * * warning message printed. * * * *--------------------------------------------------------------------* * * * Furthermore, the function * * * * ADJUST(STEPS) * * * * is provided. It can be used to determine a value of SNAPS so * * that the increase in spatial complexity equals approximately the * * increase in temporal complexity. For that ADJUST computes a * * return value satisfying SNAPS ~= log_4 (STEPS) because of the * * theory developed in the paper mentioned above. * * * *--------------------------------------------------------------------*/ #include "revolve.h" #include #include #define checkup 64 #define repsup 64 #define MAXINT 2147483647 struct { int advances; int takeshots; int commands; } numbers; /* ************************************************************************* */ int numforw(int steps, int snaps) { int reps, range, num; if (snaps < 1) { printf(" error occurs in numforw: snaps < 1\n"); return -1; } if (snaps > checkup) { printf(" number of snaps=%d exceeds checkup \n",snaps); printf(" redefine 'checkup' \n"); return -1; } reps = 0; range = 1; while(range < steps) { reps += 1; range = range*(reps + snaps)/reps; } if (reps > repsup) { printf(" number of reps=%d exceeds repsup \n",reps); printf(" redefine 'repsup' \n"); return -1; } num = reps * steps - range*reps/(snaps+1); return num; } /* ************************************************************************* */ double expense(int steps, int snaps) { double ratio; if (snaps < 1) { printf(" error occurs in expense: snaps < 0\n"); return -1; } if (steps < 1) { printf(" error occurs in expense: steps < 0\n"); return -1; } ratio = ((double) numforw(steps,snaps)); if (ratio == -1) return -1; ratio = ratio/steps; return ratio; } /* ************************************************************************* */ int maxrange(int ss, int tt) { int i, ires; double res = 1.0; if((tt<0) || (ss<0)) { printf("error in MAXRANGE: negative parameter"); return -1; } for(i=1; i<= tt; i++) { res *= (ss + i); res /= i; if (res > MAXINT) { ires=MAXINT; printf("warning from MAXRANGE: returned maximal integer %d\n",ires); return ires; } } ires = res; return ires; } /* ************************************************************************* */ int adjust(int steps) { int snaps, s, reps; snaps = 1; reps = 1; s = 0; while( maxrange(snaps+s, reps+s) > steps ) s--; while( maxrange(snaps+s, reps+s) < steps ) s++; snaps += s; reps += s ; s = -1; while( maxrange(snaps,reps) >= steps ) { if (snaps > reps) { snaps -= 1; s = 0; } else { reps -= 1; s = 1; } } if ( s == 0 ) snaps += 1 ; if ( s == 1 ) reps += 1; return snaps; } /* ************************************************************************* */ enum action revolve(int* check,int* capo,int* fine,int snaps,int* info) { static int turn, reps, range, ch[checkup], oldsnaps, oldfine; int ds, oldcapo, num, bino1, bino2, bino3, bino4, bino5; /* (*capo,*fine) is the time range currently under consideration */ /* ch[j] is the number of the state that is stored in checkpoint j */ numbers.commands += 1; if ((*check < -1) || (*capo > *fine)) return error; if ((*check == -1) && (*capo < *fine)) { if (*check == -1) turn = 0; /* initialization of turn counter */ *ch = *capo-1; } switch(*fine-*capo) { case 0: /* reduce capo to previous checkpoint, unless done */ if(*check == -1 || *capo==*ch ) { *check -= 1; if (*info > 0) { printf(" \n advances: %5d",numbers.advances); printf(" \n takeshots: %4d",numbers.takeshots); printf(" \n commands: %5d \n",numbers.commands); } return terminate; } else { *capo = ch[*check]; oldfine = *fine; return restore; } case 1: /* (possibly first) combined forward/reverse step */ *fine -= 1; if(*check >= 0 && ch[*check] == *capo) *check -= 1; if(turn == 0) { turn = 1; oldfine = *fine; return firsturn; } else { oldfine = *fine; return youturn; } default: if(*check == -1 || ch[*check] != *capo) { *check += 1 ; if(*check >= checkup) { *info = 10; return error; } if(*check+1 > snaps) { *info = 11; return error; } ch[*check] = *capo; if (*check == 0) { numbers.advances = 0; numbers.takeshots = 0; numbers.commands = 1; oldsnaps = snaps; if (snaps > checkup) { *info = 14; return error; } if (*info > 0) { num = numforw(*fine-*capo,snaps); if (num == -1) { *info = 12; return error; } printf(" prediction of needed forward steps: %8d => \n",num); printf(" slowdown factor: %8.4f \n\n",((double) num)/(*fine-*capo)); } } numbers.takeshots += 1; oldfine = *fine; return takeshot; } else { if ((oldfine < *fine) && (snaps == *check+1)) { *info = 13; return error; } oldcapo = *capo; ds = snaps - *check; if (ds < 1) { *info = 11; return error; } reps = 0; range = 1; while(range < *fine - *capo) { reps += 1; range = range*(reps + ds)/reps; } if (reps > repsup) { *info = 15; return error; } if (snaps != oldsnaps) { if (snaps > checkup) { *info = 14; return error; } } bino1 = range*reps/(ds+reps); bino2 = (ds > 1) ? bino1*ds/(ds+reps-1) : 1; if (ds == 1) bino3 = 0; else bino3 = (ds > 2) ? bino2*(ds-1)/(ds+reps-2) : 1; bino4 = bino2*(reps-1)/ds; if (ds < 3) bino5 = 0; else bino5 = (ds > 3) ? bino3*(ds-2)/reps : 1; if (*fine-*capo <= bino1 + bino3) *capo = *capo+bino4; else { if (*fine-*capo >= range - bino5) *capo = *capo + bino1; else *capo = *fine-bino2-bino3; } if (*capo == oldcapo) *capo = oldcapo+1; numbers.advances = numbers.advances + *capo - oldcapo; oldfine = *fine; return advance; } } free(ch); } SHAR_EOF fi # end of overwriting check cd .. cd .. if test ! -d 'Doc' then mkdir 'Doc' fi cd 'Doc' if test -f 'makefile_c' then echo shar: will not over-write existing file "'makefile_c'" else cat << SHAR_EOF > 'makefile_c' CC = gcc MCC = gcc LIBS = -lm all: driver src.o : src.c $(CC) -c -O3 src.c driver : src.o driver.o $(CC) -o driver src.o driver.o $(LIBS) driver.o : driver.c $(CC) -c -O3 driver.c clean: rm *.o; rm driver; SHAR_EOF fi # end of overwriting check if test -f 'makefile_f' then echo shar: will not over-write existing file "'makefile_f'" else cat << SHAR_EOF > 'makefile_f' CC = f77 MCC = f77 LIBS = -lm all: driver src.o : src.f $(CC) -c src.f driver : src.o driver.o $(CC) -o driver src.o driver.o $(LIBS) driver.o : driver.f $(CC) -c driver.f clean: rm *.o; rm driver; SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Fortran77' then mkdir 'Fortran77' fi cd 'Fortran77' if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test -f 'data' then echo shar: will not over-write existing file "'data'" else cat << SHAR_EOF > 'data' 14 3 3 SHAR_EOF fi # end of overwriting check if test -f 'driver.f' then echo shar: will not over-write existing file "'driver.f'" else cat << SHAR_EOF > 'driver.f' PROGRAM MAIN C C .. Parameters .. REAL TAKSHT,ADVAN,FSTURN,YUTURN PARAMETER (TAKSHT=1,ADVAN=2,FSTURN=3,YUTURN=4) REAL RESTRE,TRMATE,ERROR PARAMETER (RESTRE=5,TRMATE=6,ERROR=7) C .. C .. Local Scalars .. INTEGER CAPO,CHECK,FINE,INFO,SNAPS,STEPS,WHATDO C .. C .. External Functions .. INTEGER REVOLV EXTERNAL REVOLV C .. WRITE (*,FMT=*) 'ENTER: STEPS, SNAPS, INFO' READ (*,FMT=*) STEPS,SNAPS,INFO WRITE (*,FMT=*) STEPS,SNAPS,INFO CAPO = 0 FINE = STEPS + CAPO CHECK = -1 10 CONTINUE WHATDO = REVOLV(CHECK,CAPO,FINE,SNAPS,INFO) IF ((WHATDO.EQ.TAKSHT) .AND. (INFO.GT.1)) THEN WRITE (*,FMT=9000) CAPO END IF IF ((WHATDO.EQ.ADVAN) .AND. (INFO.GT.2)) THEN WRITE (*,FMT=9010) CAPO END IF IF ((WHATDO.EQ.FSTURN) .AND. (INFO.GT.2)) THEN WRITE (*,FMT=9020) CAPO END IF IF ((WHATDO.EQ.YUTURN) .AND. (INFO.GT.2)) THEN WRITE (*,FMT=9030) CAPO END IF IF ((WHATDO.EQ.RESTRE) .AND. (INFO.GT.2)) THEN WRITE (*,FMT=9040) CAPO END IF IF (WHATDO.EQ.ERROR) THEN WRITE (*,FMT=*) ' irregular termination of treeverse' IF (INFO.EQ.10) THEN WRITE (*,FMT=*) + ' number of checkpoints stored exceeds CHEKUP,' WRITE (*,FMT=*) ' increase constant CHEKUP and recompile' END IF IF (INFO.EQ.11) THEN WRITE (*,FMT=*) ' number of checkpoints stored = ', + CHECK + 1,' exceeds SNAPS,' WRITE (*,FMT=*) ' ensure SNAPS > 0 and ', + 'increase initial FINE' END IF IF (INFO.EQ.12) WRITE (*,FMT=*) ' error occurs in NUMFRW' IF (INFO.EQ.13) THEN WRITE (*,FMT=*) ' enhancement of FINE, SNAPS = ', + CHECK + 1,'checkpoints stored, increase SNAPS' END IF IF (INFO.EQ.14) THEN WRITE (*,FMT=*) ' number of SNAPS = ',SNAPS, + ' exceeds CHEKUP,' WRITE (*,FMT=*) ' increase constant CHEKUP and recompile' END IF IF (INFO.EQ.15) THEN WRITE (*,FMT=*) ' number of reps exceeds REPSUP, ' WRITE (*,FMT=*) ' increase constant REPSUP and recompile' END IF END IF IF ((WHATDO.EQ.TRMATE) .OR. (WHATDO.EQ.ERROR)) THEN GO TO 20 ELSE GO TO 10 END IF 20 CONTINUE 9000 FORMAT (' takeshot at',I6) 9010 FORMAT (' advance to',I7) 9020 FORMAT (' firsturn at',I6) 9030 FORMAT (' youturn at',I7) 9040 FORMAT (' restore at',I7) END SHAR_EOF fi # end of overwriting check if test -f 'res' then echo shar: will not over-write existing file "'res'" else cat << SHAR_EOF > 'res' ENTER: STEPS, SNAPS, INFO 14 3 3 prediction of needed forward steps: 27 => slowdown factor: 1.9286 takeshot at 0 advance to 5 takeshot at 5 advance to 10 takeshot at 10 advance to 13 firsturn at 13 restore at 10 advance to 12 youturn at 12 restore at 10 advance to 11 youturn at 11 restore at 10 youturn at 10 restore at 5 advance to 7 takeshot at 7 advance to 9 youturn at 9 restore at 7 advance to 8 youturn at 8 restore at 7 youturn at 7 restore at 5 advance to 6 youturn at 6 restore at 5 youturn at 5 restore at 0 advance to 1 takeshot at 1 advance to 2 takeshot at 2 advance to 4 youturn at 4 restore at 2 advance to 3 youturn at 3 restore at 2 youturn at 2 restore at 1 youturn at 1 restore at 0 youturn at 0 advances: 27 takeshots: 6 commands: 47 SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test -f 'src.f' then echo shar: will not over-write existing file "'src.f'" else cat << SHAR_EOF > 'src.f' C The function REVOLV coded below is meant to be used as a * C "controller" for running a time-dependent applications program * C in the reverse mode with checkpointing described in the paper * C "Achieving logarithmic Growth in temporal and spatial complexity * C in reverse automatic differentiation", Optimization Methods and * C Software, Vol.1 pp. 35-54. * C A postscript source of that paper can be found in the ftp sites * C info.mcs.anl.gov and nbtf02.math.tu-dresden.de. * C Apart from REVOLV this file contains five auxiliary routines * C NUMFRW, EXPNS, MAXRGE, and ADJUST. * C * C--------------------------------------------------------------------* C * C To utilize REVOLV the user must have procedures for * C - Advancing the state of the modeled system to a certain time. * C - Saving the current state onto a stack of snapshots. * C - Restoring the the most recently saved snapshot and * C restarting the forward simulation from there. * C - Initializing the adjoints at the end of forward sweep. * C - Performing one combined forward and adjoint step. * C Through an encoding of its return value REVOLV asks the calling * C program to perform one of these 'actions', which we will * C refer to as * C * C 'advance', 'takeshot', 'restore', 'firsturn' and 'youturn' .* C There are two other return values, namely * C 'terminate' and 'error' * C which indicate a regular or faulty termination of the calls * C to REVOLV. * C * C The action 'firsturn' includes a 'youturn', in that it requires * C -advancing through the last time-step with recording * C of intermediates * C -initializing the adjoint values (possibly after * C performing some IO) * C -reversing the last time step using the record just written * C The action 'firsturn' is obtained when the difference FINE-CAPO * C has been reduced to 1 for the first time. * C * C--------------------------------------------------------------------* C * C The calling sequence is * C * C REVOLV(CHECK,CAPO,FINE,SNAPS,INFO) * C * C with the return value being one of the actions to be taken. The * C calling parameters are all integers with the following meaning * C * C CHECK number of checkpoint being written or retrieved * C CAPO beginning of subrange currently being processed * C FINE end of subrange currently being processed * C SNAPS upper bound on number of checkpoints taken * C INFO determines how much information will be printed * C and contains information about an error occured * C * C Since REVOLV involves only a few integer operations its run-time * C is truly negligible within any nontrivial application. * C * C The parameter SNAPS is selected by the user (possibly with the * C help of the routines EXPNS and ADJUST described below ) and * C remains unchanged throughout. * C * C The pair (CAPO,FINE) always represents the initial and final * C state of the subsequence of time steps currently being traversed * C backwards. * C * C The conditions * C CHECK >= -1 and CAPO <= FINE * C are necessary and sufficient for a regular response of REVOLV. * C If either condition is violated the value 'error' is returned. * C * C The first call to REVOLV must be with CHECK=-1 so that * C appropriate initializations can be performed internally. * C * C When CHECK =-1 and CAPO = FINE then 'terminate' is returned as * C action value. This combination necessarily arises after a * C sufficiently large number of calls to REVOLV, which depends only * C on the initial difference FINE-CAPO. * C * C The last parameter INFO determines how much information about * C the actions performed will be printed. When INFO =0 no * C information is sent to standard output. When INFO > 0 REVOLV * C produces an output that contains a prediction of the number of * C forward steps and of the factor by which the execution will slow * C down. When an error occurs, the return value of INFO contains * C information about the reason: * C INFO = 10: number of checkpoints stored exceeds CHECKUP, * C increase constant CHECKUP and recompile * C INFO = 11: number of checkpoints stored exceeds SNAPS, ensure * C SNAPS greater than 0 and increase initial FINE * C INFO = 12: error occurs in NUMFORW * C INFO = 13: enhancement of FINE, SNAPS checkpoints stored, * C SNAPS must be increased * C INFO = 14: number of SNAPS exceeds CHECKUP, increase constant * C CHECKUP and recompile * C INFO = 15: number of REPS exceeds REPSUP, increase constant * C REPSUP and recompile * C * C * C--------------------------------------------------------------------* C * C Some further explanations and motivations: * C * C There is an implicit bound on CHECK through the dimensioning of * C the integer array CH[CHEKUP] with CHECKUP = 64 being the default.* C If anybody wants to have that even larger he must change the * C source. Also for the variable REPS an upper bound REPSUP is * C defined. The default value equals 64. If during a call to * C TREEVERSE a (CHECKUP+1)-st checkpoint would normally be called * C for then control is returned after an appropriate error message. * C When the calculated REPS exceeds REPSUP also an error message * C occurs. * C During the forward sweep the user is free to change the last * C three parameters from call to call, except that FINE may never * C be less than the current value of CAPO. This may be useful when * C the total number of time STEPS to be taken is not a priori * C known. The choice FINE=CAPO+1 initiates the reverse sweep, which * C happens automatically if is left constant as CAPO is eventually * C moved up to it. Once the first reverse or restore action has * C been taken only the last two parameters should be changed. * C * C--------------------------------------------------------------------* C * C The necessary number of forward steps without recording is * C calculated by the function * C * C NUMFRW(STEPS,SNAPS) * C * C STEPS denotes the total number of time steps, i.e. FINE-CAPO * C during the first call of REVOLV. When SNAPS is less than 1 an * C error message will be given and -1 is returned as value. * C * C--------------------------------------------------------------------* C * C To choose an appropriated value of SNAPS the function * C * C EXPNS(STEPS,SNAPS) * C * C estimates the run-time factor incurred by REVOLV for a * C particular value of SNAPS. The ratio NUMFRW(STEPS,SNAPS)/STEPS * C is returned. This ratio corresponds to the run-time factor of * C the execution relative to the run-time of one forward time step. * C * C--------------------------------------------------------------------* C * C The auxiliary function * C * C MAXRGE(SNAPS,REPS) * C * C returns the integer (SNAPS+REPS)!/(SNAPS!REPS!) provided * C SNAPS >=0, REPS >= 0. Otherwise there will be appropriate error * C messages and the value -1 will be returned. If the binomial * C expression is not representable as a signed 4 byte integer, * C greater than 2^31-1, this maximal value is returned and a * C warning message printed. * C * C--------------------------------------------------------------------* C * C Furthermore, the function * C * C ADJUST(STEPS) * C * C is provided. It can be used to determine a value of SNAPS so * C that the increase in spatial complexity equals approximately the * C increase in temporal complexity. For that ADJUST computes a * C return value satisfying SNAPS ~= log_4 (STEPS) because of the * C theory developed in the paper mentioned above. * C * C--------------------------------------------------------------------* INTEGER FUNCTION NUMFRW(STEPS,SNAPS) C .. Parameters .. INTEGER CHEKUP,REPSUP PARAMETER (CHEKUP=64,REPSUP=64) C .. C .. Scalar Arguments .. INTEGER SNAPS,STEPS C .. C .. Local Scalars .. INTEGER RANGE,REPS C .. IF (SNAPS.LT.1) THEN NUMFRW = -1 ELSE REPS = 0 RANGE = 1 10 CONTINUE IF (RANGE.LT.STEPS) THEN REPS = REPS + 1 RANGE = RANGE* (REPS+SNAPS)/REPS GO TO 10 END IF IF (REPS.GT.REPSUP) THEN WRITE (*,FMT=*) ' number of reps exceeds REPSUP, ' WRITE (*,FMT=*) ' increase constant REPSUP and recompile' NUMFRW = -1 RETURN END IF IF (SNAPS.LE.CHEKUP) THEN NUMFRW = REPS*STEPS - RANGE*REPS/(SNAPS+1) END IF END IF END C--------------------------------------------------------------------* DOUBLE PRECISION FUNCTION EXPNS(STEPS,SNAPS) C .. Scalar Arguments .. INTEGER SNAPS,STEPS C .. C .. External Functions .. INTEGER NUMFRW EXTERNAL NUMFRW C .. C .. Intrinsic Functions .. INTRINSIC DBLE C .. IF (SNAPS.LT.1) THEN WRITE (*,FMT=*) 'error occurs in EXPNS: SNAPS < 1' EXPNS = -1 ELSE IF (SNAPS.LT.1) THEN WRITE (*,FMT=*) 'error occurs in EXPNS: SNAPS < 1' EXPNS = -1 ELSE EXPNS = DBLE(NUMFRW(STEPS,SNAPS)) IF (EXPNS.NE.-1) THEN EXPNS = EXPNS/STEPS END IF END IF END IF END C--------------------------------------------------------------------* INTEGER FUNCTION MAXRGE(SS,TT) C .. Scalar Arguments .. INTEGER SS,TT C .. C .. Local Scalars .. DOUBLE PRECISION RES INTEGER I C .. RES = 1. IF (TT.LT.0 .OR. SS.LT.0) THEN WRITE (*,FMT=*) 'error in MAXRGE: negative parameter ' MAXRGE = -1 GO TO 30 END IF DO 10 I = 1,TT RES = RES* (SS+I) RES = RES/I IF (RES.GE.2.0**31) GO TO 20 10 CONTINUE 20 CONTINUE IF (RES.LT.2.0**31-2) THEN MAXRGE = RES ELSE MAXRGE = 2.0**31 - 3 WRITE (*,FMT=*) + 'warning from MAXRGE: returned maximal integer' WRITE (*,FMT=*) MAXRGE END IF 30 CONTINUE RETURN END C--------------------------------------------------------------------* INTEGER FUNCTION ADJUST(STEPS) C .. Scalar Arguments .. INTEGER STEPS C .. C .. Local Scalars .. INTEGER REPS,S,SNAPS C .. C .. External Functions .. INTEGER MAXRGE EXTERNAL MAXRGE C .. SNAPS = 1 REPS = 1 S = 0 10 IF (MAXRGE(SNAPS+S,REPS+S).GT.STEPS) THEN S = S - 1 GO TO 10 END IF 20 IF (MAXRGE(SNAPS+S,REPS+S).LT.STEPS) THEN S = S + 1 GO TO 20 END IF SNAPS = SNAPS + S REPS = REPS + S S = -1 30 IF (MAXRGE(SNAPS,REPS).GE.STEPS) THEN IF (SNAPS.GT.REPS) THEN SNAPS = SNAPS - 1 S = 0 ELSE REPS = REPS - 1 S = 1 END IF GO TO 30 END IF IF (S.EQ.0) SNAPS = SNAPS + 1 IF (S.EQ.1) REPS = REPS + 1 ADJUST = SNAPS RETURN END C--------------------------------------------------------------------* INTEGER FUNCTION REVOLV(CHECK,CAPO,FINE,SNAPS,INFO) C (CAPO ,FINE) is the time range currently under consideration C .. Parameters .. INTEGER CHEKUP,REPSUP PARAMETER (CHEKUP=64,REPSUP=64) INTEGER TAKSHT,ADVAN,FSTURN,YUTURN PARAMETER (TAKSHT=1,ADVAN=2,FSTURN=3,YUTURN=4) INTEGER RESTRE,TRMATE,ERROR PARAMETER (RESTRE=5,TRMATE=6,ERROR=7) C .. C .. Scalar Arguments .. INTEGER CAPO,CHECK,FINE,INFO,REPS,SNAPS C .. C .. Scalars in Common .. INTEGER NUMADV,NUMCMD,NUMTKS,OLDFIN,OLDSNP,TURN C .. C .. Arrays in Common .. INTEGER CH(0:CHEKUP) C .. C .. Local Scalars .. INTEGER DS,I,NUM,OLDCPO,RANGE REAL BINO1, BINO2, BINO3, BINO4, BINO5 C .. C .. External Functions .. INTEGER NUMFRW EXTERNAL NUMFRW C .. C .. Intrinsic Functions .. INTRINSIC DBLE C .. C .. Common blocks .. COMMON /INFNUM/NUMADV,NUMTKS,NUMCMD COMMON /SNREP/CH,TURN,OLDSNP,OLDFIN C .. IF (INFO.GT.3) INFO = 3 NUMCMD = NUMCMD + 1 IF ((CHECK.LT.-1) .AND. (CAPO.GT.FINE)) THEN REVOLV = ERROR GO TO 30 END IF IF ((CHECK.EQ.-1) .AND. (CAPO.LT.FINE)) THEN IF (CHECK.EQ.-1) TURN = 0 DO 10 I = 0,CHEKUP CH(I) = 0 10 CONTINUE CH(0) = CAPO - 1 END IF C IF ((FINE-CAPO).EQ.0) THEN IF ((CHECK.EQ. (-1)) .OR. (CAPO.EQ.CH(0))) THEN CHECK = CHECK - 1 IF (INFO.GT.0) THEN WRITE (*,FMT=*) WRITE (*,FMT=9000) NUMADV WRITE (*,FMT=9010) NUMTKS WRITE (*,FMT=9020) NUMCMD END IF REVOLV = TRMATE GO TO 30 ELSE CAPO = CH(CHECK) OLDFIN = FINE REVOLV = RESTRE GO TO 30 END IF END IF IF ((FINE-CAPO).EQ.1) THEN FINE = FINE - 1 IF ((CHECK.GE.0) .AND. (CH(CHECK).EQ.CAPO)) THEN CHECK = CHECK - 1 END IF IF (TURN.EQ.0) THEN OLDFIN = FINE REVOLV = FSTURN TURN = 1 ELSE OLDFIN = FINE REVOLV = YUTURN END IF GO TO 30 END IF IF ((CHECK.EQ. (-1)) .OR. (CH(CHECK).NE.CAPO)) THEN CHECK = CHECK + 1 IF (CHECK.GE.CHEKUP) THEN INFO = 10 REVOLV = ERROR GO TO 30 ELSE IF (CHECK+1.GT.SNAPS) THEN INFO = 11 REVOLV = ERROR GO TO 30 END IF CH(CHECK) = CAPO IF (CHECK.EQ.0) THEN NUMADV = 0 NUMTKS = 0 NUMCMD = 1 OLDSNP = SNAPS IF (SNAPS.GT.CHEKUP) THEN INFO = 14 REVOLV = ERROR GO TO 30 END IF IF (INFO.GT.0) THEN NUM = NUMFRW(FINE-CAPO,SNAPS) IF (NUM.EQ.-1) THEN INFO = 12 REVOLV = ERROR GO TO 30 END IF WRITE (*,FMT=9030) NUM WRITE (*,FMT=9040) DBLE(NUM)/ (FINE-CAPO) WRITE (*,FMT=*) END IF END IF NUMTKS = NUMTKS + 1 OLDFIN = FINE REVOLV = TAKSHT GO TO 30 END IF ELSE IF ((OLDFIN.LT.FINE) .AND. (SNAPS.EQ.CHECK+1)) THEN INFO = 13 REVOLV = ERROR GO TO 30 END IF OLDCPO = CAPO DS = SNAPS - CHECK IF (DS.LT.1) THEN INFO = 11 REVOLV = ERROR GO TO 30 END IF REPS = 0 RANGE = 1 20 CONTINUE IF (RANGE.LT.FINE-CAPO) THEN REPS = REPS + 1 RANGE = RANGE* (REPS+DS)/REPS GO TO 20 END IF IF (REPS.GT.REPSUP) THEN INFO = 15 REVOLV = ERROR GO TO 30 END IF IF (SNAPS.NE.OLDSNP) THEN IF (SNAPS.GT.CHEKUP) THEN INFO = 14 REVOLV = ERROR GO TO 30 END IF OLDSNP = SNAPS END IF BINO1 = RANGE*REPS/(DS+REPS) IF (DS .GT. 1) THEN BINO2 = BINO1*DS/(DS+REPS-1) ELSE BINO2 = 1 ENDIF IF (DS .EQ. 1) THEN BINO3 = 0 ELSE IF (DS .GT. 2) THEN BINO3 = BINO2*(DS-1)/(DS+REPS-2) ELSE BINO3 = 1 ENDIF ENDIF BINO4 = BINO2*(REPS-1)/DS IF (DS .LT. 3) THEN BINO5 = 0 ELSE IF (DS .GT. 3) THEN BINO5 = BINO3*(DS-1)/REPS ELSE BINO5 = 1 ENDIF ENDIF IF (FINE-CAPO.LE.BINO1+BINO3) THEN CAPO = CAPO + BINO4 ELSE IF (FINE-CAPO.GE.RANGE-BINO5) THEN CAPO = CAPO + BINO1 ELSE CAPO = FINE - BINO2 - BINO3 END IF END IF IF (CAPO.EQ.OLDCPO) CAPO = OLDCPO + 1 NUMADV = NUMADV + CAPO - OLDCPO REVOLV = ADVAN END IF 30 CONTINUE RETURN 9000 FORMAT (' advances:',I8) 9010 FORMAT (' takeshots:',I7) 9020 FORMAT (' commands:',I8) 9030 FORMAT (' prediction of needed forward steps: ',I7,' => ') 9040 FORMAT (' slowdown factor: ',F8.4) END SHAR_EOF fi # end of overwriting check cd .. cd .. # End of shell archive exit 0