TAP_support/ADFirstAidKit/adBinomial.c

b4daa24319 Shre*0001 /*
                0002  * TAPENADE Automatic Differentiation Engine
                0003  * Copyright (C) 1999-2021 Inria
                0004  * See the LICENSE.md file in the project root for more information.
                0005  *
                0006  */
                0007
                0008 #include <stdio.h>
                0009 #include "adBinomial.h"
                0010
                0011 static int stack1[15] ;
                0012 static int stack2[297] ;
                0013 static int is1 = -1 ;
                0014 static int is2 = -1 ;
                0015
                0016 #define BOTTOM stack1[3*is1]
                0017 #define MAXSNP stack1[3*is1+1]
                0018 #define OFFSET stack1[3*is1+2]
                0019 #define LENGTH stack2[3*is2]
                0020 #define CURPOS stack2[3*is2+1]
                0021 #define CKPCUT stack2[3*is2+2]
                0022
                0023 void adBinomial_init(int length, int nbSnap, int firstStep) {
                0024   if (length<=0) {
                0025     printf("Error: Cannot reverse a sequence of length %i\n",length) ;
                0026   } else if (nbSnap==0 && length>=2) {
                0027     printf("Error: Cannot reverse a sequence of length %i with no snapshot\n",length) ;
                0028   } else if (is1>4 || is2+nbSnap+1>98) {
                0029     printf("Error: Binomial-Checkpointing memory exceeded !\n") ;
                0030   } else {
                0031     is1 = is1+1 ;
                0032     is2 = is2+1 ;
                0033     BOTTOM = is2 ;
                0034     MAXSNP = nbSnap ;
                0035     OFFSET = firstStep ;
                0036     LENGTH= length ;
                0037     CURPOS = 0 ;
                0038     CKPCUT = 0 ;
                0039   }
                0040 }
                0041
                0042 void adBinomial_setCut() {
                0043   int length = LENGTH ;
                0044   int nbSnap = MAXSNP-(is2-BOTTOM)+1 ;
                0045   if (length<=1) {
                0046     CKPCUT = 0 ;
                0047   } else if (nbSnap==1) {
                0048     CKPCUT = length-1 ;
                0049   } else {
                0050     int minRecomp = 1 ;
                0051     int eta = nbSnap+1 ;
                0052     while (eta<length) {
                0053       minRecomp++ ;
                0054       eta = (eta*(nbSnap+minRecomp))/minRecomp ;
                0055     }
                0056     CKPCUT = (length*minRecomp)/(minRecomp+nbSnap) ;
                0057     if (CKPCUT==0)
                0058       CKPCUT=1 ;
                0059     else if (CKPCUT>=length)
                0060       CKPCUT=length-1 ;
                0061   }
                0062 }
                0063
                0064 int adBinomial_next(int *action, int *step) {
                0065   int i ;
                0066
                0067 //  Begin "only for debug" part:
                0068 //  printf("\n") ;
                0069 //  printf("STACK1:\n") ;
                0070 //  for (i=is1 ; i>=0 ; i--)
                0071 //    printf("%i snapshots, stack2 bottom:%i (offset:%i)\n",
                0072 //           stack1[3*i+1],stack1[3*i],stack1[3*i+2]) ;
                0073 //  printf("-------------------\n") ;
                0074 //  printf("STACK2:\n") ;
                0075 //  for (i=is2 ; i>=0 ; i--)
                0076 //    printf("%i: R( ,%i) %i/%i\n",
                0077 //           i,stack2[3*i],stack2[3*i+1],stack2[3*i+2]) ;
                0078 //  //    printf("%i: at %i, nextckp %i, end %i\n",
                0079 //  //           i,stack2[3*i+1],stack2[3*i+2],stack2[3*i]) ;
                0080 //  printf("-------------------\n") ;
                0081 //  End "only for debug" part.
                0082
                0083   if (LENGTH<=0 && is2==BOTTOM) {
                0084     *step = -1 ;
                0085     *action = -1 ;
                0086     is1-- ;
                0087     is2-- ;
                0088     return 0 ;
                0089   } else {
                0090     *step = 1 ;
                0091     for (i=BOTTOM+1 ; i<=is2 ; i++) *step += stack2[3*i+1] ;
                0092     if (CURPOS==-1) {
                0093       *action = (LENGTH==1)?POPSNAP:LOOKSNAP ;
                0094       CURPOS = 0 ;
                0095     } else {
                0096       if (CURPOS==LENGTH-1) {
                0097         *action = (*step==stack2[3*BOTTOM])?FIRSTTURN:TURN ;
                0098         if (CURPOS==0 && is2>BOTTOM) {
                0099           is2-- ;
                0100           LENGTH = CKPCUT ;
                0101         } else {
                0102           LENGTH = CURPOS ;
                0103         }
                0104         CURPOS = -1 ;
                0105         adBinomial_setCut() ;
                0106       } else if (CURPOS==CKPCUT) {
                0107         int remainingLength = LENGTH-CURPOS ;
                0108         *action = PUSHSNAP ;
                0109         is2++ ;
                0110         LENGTH = remainingLength ;
                0111         CURPOS = 0 ;
                0112         adBinomial_setCut() ;
                0113         (*step)-- ;
                0114       } else {
                0115         *action = ADVANCE ;
                0116         CURPOS++ ;
                0117       }
                0118     }
                0119     return 1 ;
                0120   }
                0121 }
                0122
                0123 void adBinomial_resize() {
                0124   int step = 1, i ;
                0125   for (i=BOTTOM+1 ; i<=is2 ; ++i) step += stack2[3*i+1] ;
                0126   printf("Binomial iteration exits on step %i before expected %i\n",
                0127          step-1, stack2[3*BOTTOM]) ;
                0128   stack2[3*BOTTOM] = step-1 ;
                0129   LENGTH = CURPOS ;
                0130   --(CURPOS) ;
                0131 }
                0132
                0133 /****************** INTERFACE CALLED FROM FORTRAN *******************/
                0134
                0135 void adbinomial_init_(int *length, int *nbSnap, int *firstStep) {
                0136   adBinomial_init(*length, *nbSnap, *firstStep) ;
                0137 }
                0138
                0139 int adbinomial_next_(int *action, int *step) {
                0140   return adBinomial_next(action, step) ;
                0141 }
                0142
                0143 void adbinomial_resize_() {
                0144   adBinomial_resize() ;
                0145 }