TAP_support/ADFirstAidKit/adFixedPoint.c

feb7fa5d1e dngo*0001 #include "adFixedPoint.h"
                0002 #include <stdio.h>
                0003
                0004 // Up to 5 levels of nested fixed-point loops should be enough:
                0005 double refCumuls[5] = {0.0, 0.0, 0.0, 0.0, 0.0} ;
                0006 double prevCumuls[5] = {0.0, 0.0, 0.0, 0.0, 0.0} ;
                0007 int adjIters[5] = {0, 0, 0, 0, 0} ;
                0008 int fpDepth = -1 ;
                0009
                0010 // this function was written by L Hascoet, since the function name is fixed
                0011 // and the alg is not appropriate for the STREAMICE fixed-point problem, this
                0012 // is not used but kept for analysis
                0013 /*
                0014 int adFixedPoint_notReduced(double cumul, float reduction) {
                0015   // Dan's method: stop when cumul has reduced by a factor "reduction"
                0016   if (cumul<0.0) {
                0017     // Begin 1st iteration of a new adjoint FP loop.
                0018     // Given cumul is always -1.0.
                0019     // Just prepare for 2nd iteration, that will set the reference refCumuls[] :
                0020     if (fpDepth>=(5-1)) return 0 ; // protect from going out of bounds.
                0021     ++fpDepth ;
                0022     refCumuls[fpDepth] = -1.0 ;
                0023     prevCumuls[fpDepth] = -1.0 ;
                0024     adjIters[fpDepth] = 1 ;
                0025     return 1 ;
                0026   } else {
                0027     // Begin 2nd of any following iteration:
                0028     int iterate ;
                0029     int growth ;
                0030     if (refCumuls[fpDepth]<0.0) {
                0031       // Begin 2nd iteration of current adjoint FP loop.
                0032       // Set reference refCumuls[] to the given cumul:
                0033       refCumuls[fpDepth] = cumul ;
                0034       prevCumuls[fpDepth] = cumul ;
                0035       // Almost always iterate, except if cumul is really small (quadratic convergence?):
                0036       iterate = (cumul >  1.e-10) ;
                0037       growth = 0 ;
                0038     } else {
                0039       // Begin 3rd or any following iteration.
                0040       // Compare with reference refCumuls:
                0041       iterate = (cumul > reduction*refCumuls[fpDepth]) ;
                0042       growth = (adjIters[fpDepth]>5 && cumul>prevCumuls[fpDepth]) ;
                0043       prevCumuls[fpDepth]=cumul ;
                0044     }
                0045     if (iterate && !growth) {
                0046       ++(adjIters[fpDepth]) ;
                0047       printf("%i adjoint iterations (reduced %e -> %e)\n", adjIters[fpDepth], refCumuls[fpDepth], cumul) ;
                0048       return 1 ;
                0049     } else {
                0050       if (growth) {
                0051           printf("%i adjoint iterations (reduced %e -> %e, TERMINATED)\n", adjIters[fpDepth], refCumuls[fpDepth], cumul) ; }
                0052       else {
                0053           printf("%i adjoint iterations (reduced %e -> %e), CONVERGED\n", adjIters[fpDepth], refCumuls[fpDepth], cumul) ; }
                0054
                0055       if (fpDepth<0) return 0 ; // protect from going out of bounds.
                0056       --fpDepth ;
                0057       return 0 ;
                0058     }
                0059   }
                0060 }
                0061 */
                0062
                0063 int adFixedPoint_tooLarge(double cumul, float minCumul) {
                0064   // Naive method: stop when cumul becomes less than given minCumul
                0065   if (cumul<0.0) {
                0066     // Begin 1st iteration of a new adjoint FP loop.
                0067     // Given cumul is always -1.0.
                0068     if (fpDepth>=(5-1)) return 0 ; // protect from going out of bounds.
                0069     ++fpDepth ;
                0070     adjIters[fpDepth] = 1 ;
                0071     return 1 ;
                0072   } else {
                0073     // Begin 2nd of any following iteration:
                0074     int iterate = (cumul>minCumul) ;
                0075     if (iterate) {
                0076       ++(adjIters[fpDepth]) ;
                0077       return 1 ;
                0078     } else {
                0079 /*       printf("%i adjoint iterations (residual %e < %e)\n", adjIters[fpDepth], cumul, minCumul) ; */
                0080       if (fpDepth<0) return 0 ; // protect from going out of bounds.
                0081       --fpDepth ;
                0082       return 0 ;
                0083     }
                0084   }
                0085 }
                0086
                0087 /****************** INTERFACE CALLED FROM FORTRAN *******************/
                0088
                0089 /*
                0090 int adfixedpoint_notreduced_(double *cumul, float *reduction) {
                0091   return adFixedPoint_notReduced(*cumul, *reduction) ;
                0092 }
                0093 */
                0094
                0095 int adfixedpoint_toolarge_(double *cumul, float *minCumul) {
                0096   return adFixedPoint_tooLarge(*cumul, *minCumul) ;
                0097 }