share/hompack/fortran/stepnf.f

   1       SUBROUTINE STEPNF(N,NFE,IFLAG,START,CRASH,HOLD,H,RELERR,
   2      &   ABSERR,S,Y,YP,YOLD,YPOLD,A,QR,ALPHA,TZ,PIVOT,W,WP,
   3      &   Z0,Z1,SSPAR,PAR,IPAR)
   4 C
   5 C  STEPNF  TAKES ONE STEP ALONG THE ZERO CURVE OF THE HOMOTOPY MAP
   6 C USING A PREDICTOR-CORRECTOR ALGORITHM.  THE PREDICTOR USES A HERMITE
   7 C CUBIC INTERPOLANT, AND THE CORRECTOR RETURNS TO THE ZERO CURVE ALONG
   8 C THE FLOW NORMAL TO THE DAVIDENKO FLOW.  STEPNF  ALSO ESTIMATES A
   9 C STEP SIZE H FOR THE NEXT STEP ALONG THE ZERO CURVE.  NORMALLY
  10 C  STEPNF  IS USED INDIRECTLY THROUGH  FIXPNF , AND SHOULD BE CALLED
  11 C DIRECTLY ONLY IF IT IS NECESSARY TO MODIFY THE STEPPING ALGORITHM'S
  12 C PARAMETERS.
  13 C
  14 C ON INPUT:
  15 C
  16 C N = DIMENSION OF X AND THE HOMOTOPY MAP.
  17 C
  18 C NFE = NUMBER OF JACOBIAN MATRIX EVALUATIONS.
  19 C
  20 C IFLAG = -2, -1, OR 0, INDICATING THE PROBLEM TYPE.
  21 C
  22 C START = .TRUE. ON FIRST CALL TO  STEPNF , .FALSE. OTHERWISE.
  23 C
  24 C HOLD = ||Y - YOLD||; SHOULD NOT BE MODIFIED BY THE USER.
  25 C
  26 C H = UPPER LIMIT ON LENGTH OF STEP THAT WILL BE ATTEMPTED.  H  MUST BE
  27 C    SET TO A POSITIVE NUMBER ON THE FIRST CALL TO  STEPNF .
  28 C    THEREAFTER  STEPNF  CALCULATES AN OPTIMAL VALUE FOR  H , AND  H
  29 C    SHOULD NOT BE MODIFIED BY THE USER.
  30 C
  31 C RELERR, ABSERR = RELATIVE AND ABSOLUTE ERROR VALUES.  THE ITERATION IS
  32 C    CONSIDERED TO HAVE CONVERGED WHEN A POINT W=(LAMBDA,X) IS FOUND
  33 C    SUCH THAT
  34 C
  35 C    ||Z|| <= RELERR*||W|| + ABSERR  ,          WHERE
  36 C
  37 C    Z IS THE NEWTON STEP TO W=(LAMBDA,X).
  38 C
  39 C S = (APPROXIMATE) ARC LENGTH ALONG THE HOMOTOPY ZERO CURVE UP TO
  40 C    Y(S) = (LAMBDA(S), X(S)).
  41 C
  42 C Y(1:N+1) = PREVIOUS POINT (LAMBDA(S), X(S)) FOUND ON THE ZERO CURVE OF
  43 C    THE HOMOTOPY MAP.
  44 C
  45 C YP(1:N+1) = UNIT TANGENT VECTOR TO THE ZERO CURVE OF THE HOMOTOPY MAP
  46 C    AT  Y .
  47 C
  48 C YOLD(1:N+1) = A POINT BEFORE  Y  ON THE ZERO CURVE OF THE HOMOTOPY MAP.
  49 C
  50 C YPOLD(1:N+1) = UNIT TANGENT VECTOR TO THE ZERO CURVE OF THE HOMOTOPY
  51 C    MAP AT  YOLD .
  52 C
  53 C A(1:*) = PARAMETER VECTOR IN THE HOMOTOPY MAP.
  54 C
  55 C QR(1:N,1:N+2), ALPHA(1:N), TZ(1:N+1), PIVOT(1:N+1), W(1:N+1),
  56 C    WP(1:N+1)  ARE WORK ARRAYS USED FOR THE QR FACTORIZATION (IN THE
  57 C    NEWTON STEP CALCULATION) AND THE INTERPOLATION.
  58 C
  59 C Z0(1:N+1), Z1(1:N+1)  ARE WORK ARRAYS USED FOR THE ESTIMATION OF THE
  60 C    NEXT STEP SIZE  H .
  61 C
  62 C SSPAR(1:8) = (LIDEAL, RIDEAL, DIDEAL, HMIN, HMAX, BMIN, BMAX, P)  IS
  63 C    A VECTOR OF PARAMETERS USED FOR THE OPTIMAL STEP SIZE ESTIMATION.
  64 C
  65 C PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS,
  66 C    WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES
  67 C    RHO, RHOJAC.
  68 C
  69 C ON OUTPUT:
  70 C
  71 C N , A , SSPAR  ARE UNCHANGED.
  72 C
  73 C NFE  HAS BEEN UPDATED.
  74 C
  75 C IFLAG
  76 C    = -2, -1, OR 0 (UNCHANGED) ON A NORMAL RETURN.
  77 C
  78 C    = 4 IF A JACOBIAN MATRIX WITH RANK < N HAS OCCURRED.  THE
  79 C        ITERATION WAS NOT COMPLETED.
  80 C
  81 C    = 6 IF THE ITERATION FAILED TO CONVERGE.  W  CONTAINS THE LAST
  82 C        NEWTON ITERATE.
  83 C
  84 C START = .FALSE. ON A NORMAL RETURN.
  85 C
  86 C CRASH
  87 C    = .FALSE. ON A NORMAL RETURN.
  88 C
  89 C    = .TRUE. IF THE STEP SIZE  H  WAS TOO SMALL.  H  HAS BEEN
  90 C      INCREASED TO AN ACCEPTABLE VALUE, WITH WHICH  STEPNF  MAY BE
  91 C      CALLED AGAIN.
  92 C
  93 C    = .TRUE. IF  RELERR  AND/OR  ABSERR  WERE TOO SMALL.  THEY HAVE
  94 C      BEEN INCREASED TO ACCEPTABLE VALUES, WITH WHICH  STEPNF  MAY
  95 C      BE CALLED AGAIN.
  96 C
  97 C HOLD = ||Y - YOLD||.
  98 C
  99 C H = OPTIMAL VALUE FOR NEXT STEP TO BE ATTEMPTED.  NORMALLY  H  SHOULD
 100 C    NOT BE MODIFIED BY THE USER.
 101 C
 102 C RELERR, ABSERR  ARE UNCHANGED ON A NORMAL RETURN.
 103 C
 104 C S = (APPROXIMATE) ARC LENGTH ALONG THE ZERO CURVE OF THE HOMOTOPY MAP
 105 C    UP TO THE LATEST POINT FOUND, WHICH IS RETURNED IN  Y .
 106 C
 107 C Y, YP, YOLD, YPOLD  CONTAIN THE TWO MOST RECENT POINTS AND TANGENT
 108 C    VECTORS FOUND ON THE ZERO CURVE OF THE HOMOTOPY MAP.
 109 C
 110 C
 111 C CALLS  D1MACH , DNRM2 , TANGNF .
 112 C
 113       DOUBLE PRECISION ABSERR,D1MACH,DCALC,DD001,DD0011,DD01,
 114      &   DD011,DELS,DNRM2,F0,F1,FOURU,FP0,FP1,H,HFAIL,HOLD,HT,
 115      &   LCALC,QOFS,RCALC,RELERR,RHOLEN,S,TEMP,TWOU
 116       INTEGER IFLAG,ITNUM,J,JUDY,LITFH,N,NFE,NP1
 117       LOGICAL CRASH,FAIL,START
 118 C
 119 C ***** ARRAY DECLARATIONS. *****
 120 C
 121       DOUBLE PRECISION Y(N+1),YP(N+1),YOLD(N+1),YPOLD(N+1),A(N),
 122      &   QR(N,N+2),ALPHA(N),TZ(N+1),W(N+1),WP(N+1),Z0(N+1),
 123      &   Z1(N+1),SSPAR(8),PAR(1)
 124       INTEGER PIVOT(N+1),IPAR(1)
 125 C
 126 C ***** END OF DIMENSIONAL INFORMATION. *****
 127 C
 128 C THE LIMIT ON THE NUMBER OF NEWTON ITERATIONS ALLOWED BEFORE REDUCING
 129 C THE STEP SIZE  H  MAY BE CHANGED BY CHANGING THE FOLLOWING PARAMETER
 130 C STATEMENT:
 131       PARAMETER (LITFH=4)
 132 C
 133 C DEFINITION OF HERMITE CUBIC INTERPOLANT VIA DIVIDED DIFFERENCES.
 134 C
 135       DD01(F0,F1,DELS)=(F1-F0)/DELS
 136       DD001(F0,FP0,F1,DELS)=(DD01(F0,F1,DELS)-FP0)/DELS
 137       DD011(F0,F1,FP1,DELS)=(FP1-DD01(F0,F1,DELS))/DELS
 138       DD0011(F0,FP0,F1,FP1,DELS)=(DD011(F0,F1,FP1,DELS) -
 139      &                            DD001(F0,FP0,F1,DELS))/DELS
 140       QOFS(F0,FP0,F1,FP1,DELS,S)=((DD0011(F0,FP0,F1,FP1,DELS)*(S-DELS) +
 141      &   DD001(F0,FP0,F1,DELS))*S + FP0)*S + F0
 142 C
 143 C
 144       TWOU=2.0*D1MACH(4)
 145       FOURU=TWOU+TWOU
 146       NP1=N+1
 147       CRASH=.TRUE.
 148 C THE ARCLENGTH  S  MUST BE NONNEGATIVE.
 149       IF (S .LT. 0.0) RETURN
 150 C IF STEP SIZE IS TOO SMALL, DETERMINE AN ACCEPTABLE ONE.
 151       IF (H .LT. FOURU*(1.0+S)) THEN
 152         H=FOURU*(1.0+S)
 153         RETURN
 154       ENDIF
 155 C IF ERROR TOLERANCES ARE TOO SMALL, INCREASE THEM TO ACCEPTABLE VALUES.
 156       TEMP=DNRM2(NP1,Y,1)
 157       IF (.5*(RELERR*TEMP+ABSERR) .GE. TWOU*TEMP) GO TO 40
 158       IF (RELERR .NE. 0.0) THEN
 159         RELERR=FOURU*(1.0+FOURU)
 160         ABSERR=MAX(ABSERR,0.0D0)
 161       ELSE
 162         ABSERR=FOURU*TEMP
 163       ENDIF
 164       RETURN
 165  40   CRASH=.FALSE.
 166       IF (.NOT. START) GO TO 300
 167 C
 168 C *****  STARTUP SECTION(FIRST STEP ALONG ZERO CURVE.  *****
 169 C
 170       FAIL=.FALSE.
 171       START=.FALSE.
 172 C DETERMINE SUITABLE INITIAL STEP SIZE.
 173       H=MIN(H, .10D0, SQRT(SQRT(RELERR*TEMP+ABSERR)))
 174 C USE LINEAR PREDICTOR ALONG TANGENT DIRECTION TO START NEWTON ITERATION.
 175       YPOLD(1)=1.0
 176       DO 50 J=2,NP1
 177         YPOLD(J)=0.0
 178  50   CONTINUE
 179       CALL TANGNF(S,Y,YP,YPOLD,A,QR,ALPHA,TZ,PIVOT,NFE,N,IFLAG,
 180      &            PAR,IPAR)
 181       IF (IFLAG .GT. 0) RETURN
 182  70   DO 80 J=1,NP1
 183         TEMP=Y(J) + H * YP(J)
 184         W(J)=TEMP
 185         Z0(J)=TEMP
 186  80   CONTINUE
 187       DO 200 JUDY=1,LITFH
 188         RHOLEN=-1.0
 189 C CALCULATE THE NEWTON STEP  TZ  AT THE CURRENT POINT  W .
 190         CALL TANGNF(RHOLEN,W,WP,YPOLD,A,QR,ALPHA,TZ,PIVOT,NFE,N,IFLAG,
 191      &              PAR,IPAR)
 192         IF (IFLAG .GT. 0) RETURN
 193 C
 194 C TAKE NEWTON STEP AND CHECK CONVERGENCE.
 195         DO 90 J=1,NP1
 196           W(J)=W(J) + TZ(J)
 197  90     CONTINUE
 198         ITNUM=JUDY
 199 C COMPUTE QUANTITIES USED FOR OPTIMAL STEP SIZE ESTIMATION.
 200         IF (JUDY .EQ. 1) THEN
 201           LCALC=DNRM2(NP1,TZ,1)
 202           RCALC=RHOLEN
 203           DO 110 J=1,NP1
 204             Z1(J)=W(J)
 205 110       CONTINUE
 206         ELSE IF (JUDY .EQ. 2) THEN
 207           LCALC=DNRM2(NP1,TZ,1)/LCALC
 208           RCALC=RHOLEN/RCALC
 209         ENDIF
 210 C GO TO MOP-UP SECTION AFTER CONVERGENCE.
 211         IF (DNRM2(NP1,TZ,1) .LE. RELERR*DNRM2(NP1,W,1)+ABSERR)
 212      &                                                 GO TO 600
 213 C
 214 200   CONTINUE
 215 C
 216 C NO CONVERGENCE IN  LITFH  ITERATIONS.  REDUCE  H  AND TRY AGAIN.
 217       IF (H .LE. FOURU*(1.0 + S)) THEN
 218         IFLAG=6
 219         RETURN
 220       ENDIF
 221       H=.5 * H
 222       GO TO 70
 223 C
 224 C ***** END OF STARTUP SECTION. *****
 225 C
 226 C ***** PREDICTOR SECTION. *****
 227 C
 228 300   FAIL=.FALSE.
 229 C COMPUTE POINT PREDICTED BY HERMITE INTERPOLANT.  USE STEP SIZE  H
 230 C COMPUTED ON LAST CALL TO  STEPNF .
 231 320   DO 330 J=1,NP1
 232         TEMP=QOFS(YOLD(J),YPOLD(J),Y(J),YP(J),HOLD,HOLD+H)
 233         W(J)=TEMP
 234         Z0(J)=TEMP
 235 330   CONTINUE
 236 C
 237 C ***** END OF PREDICTOR SECTION. *****
 238 C
 239 C ***** CORRECTOR SECTION. *****
 240 C
 241       DO 500 JUDY=1,LITFH
 242         RHOLEN=-1.0
 243 C CALCULATE THE NEWTON STEP  TZ  AT THE CURRENT POINT  W .
 244         CALL TANGNF(RHOLEN,W,WP,YP,A,QR,ALPHA,TZ,PIVOT,NFE,N,IFLAG,
 245      &              PAR,IPAR)
 246         IF (IFLAG .GT. 0) RETURN
 247 C
 248 C TAKE NEWTON STEP AND CHECK CONVERGENCE.
 249         DO 420 J=1,NP1
 250            W(J)=W(J) + TZ(J)
 251 420     CONTINUE
 252         ITNUM=JUDY
 253 C COMPUTE QUANTITIES USED FOR OPTIMAL STEP SIZE ESTIMATION.
 254         IF (JUDY .EQ. 1) THEN
 255           LCALC=DNRM2(NP1,TZ,1)
 256           RCALC=RHOLEN
 257           DO 440 J=1,NP1
 258             Z1(J)=W(J)
 259 440       CONTINUE
 260         ELSE IF (JUDY .EQ. 2) THEN
 261           LCALC=DNRM2(NP1,TZ,1)/LCALC
 262           RCALC=RHOLEN/RCALC
 263         ENDIF
 264 C GO TO MOP-UP SECTION AFTER CONVERGENCE.
 265         IF (DNRM2(NP1,TZ,1) .LE. RELERR*DNRM2(NP1,W,1)+ABSERR)
 266      &                                                 GO TO 600
 267 C
 268 500   CONTINUE
 269 C
 270 C NO CONVERGENCE IN  LITFH  ITERATIONS.  RECORD FAILURE AT CALCULATED  H ,
 271 C SAVE THIS STEP SIZE, REDUCE  H  AND TRY AGAIN.
 272       FAIL=.TRUE.
 273       HFAIL=H
 274       IF (H .LE. FOURU*(1.0 + S)) THEN
 275         IFLAG=6
 276         RETURN
 277       ENDIF
 278       H=.5 * H
 279       GO TO 320
 280 C
 281 C ***** END OF CORRECTOR SECTION. *****
 282 C
 283 C ***** MOP-UP SECTION. *****
 284 C
 285 C YOLD  AND  Y  ALWAYS CONTAIN THE LAST TWO POINTS FOUND ON THE ZERO
 286 C CURVE OF THE HOMOTOPY MAP.  YPOLD  AND  YP  CONTAIN THE TANGENT
 287 C VECTORS TO THE ZERO CURVE AT  YOLD  AND  Y , RESPECTIVELY.
 288 C
 289 600   DO 620 J=1,NP1
 290         YOLD(J)=Y(J)
 291         YPOLD(J)=YP(J)
 292         Y(J)=W(J)
 293         YP(J)=WP(J)
 294         W(J)=Y(J) - YOLD(J)
 295 620   CONTINUE
 296 C UPDATE ARC LENGTH.
 297       HOLD=DNRM2(NP1,W,1)
 298       S=S+HOLD
 299 C
 300 C ***** END OF MOP-UP SECTION. *****
 301 C
 302 C ***** OPTIMAL STEP SIZE ESTIMATION SECTION. *****
 303 C
 304 C CALCULATE THE DISTANCE FACTOR  DCALC .
 305 700   DO 710 J=1,NP1
 306         TZ(J)=Z0(J) - Y(J)
 307         W(J)=Z1(J) - Y(J)
 308 710   CONTINUE
 309       DCALC=DNRM2(NP1,TZ,1)
 310       IF (DCALC .NE. 0.0) DCALC=DNRM2(NP1,W,1)/DCALC
 311 C
 312 C THE OPTIMAL STEP SIZE HBAR IS DEFINED BY
 313 C
 314 C   HT=HOLD * [MIN(LIDEAL/LCALC, RIDEAL/RCALC, DIDEAL/DCALC)]**(1/P)
 315 C
 316 C     HBAR = MIN [ MAX(HT, BMIN*HOLD, HMIN), BMAX*HOLD, HMAX ]
 317 C
 318 C IF CONVERGENCE HAD OCCURRED AFTER 1 ITERATION, SET THE CONTRACTION
 319 C FACTOR  LCALC  TO ZERO.
 320       IF (ITNUM .EQ. 1) LCALC = 0.0
 321 C FORMULA FOR OPTIMAL STEP SIZE.
 322       IF (LCALC+RCALC+DCALC .EQ. 0.0) THEN
 323         HT = SSPAR(7) * HOLD
 324       ELSE
 325         HT = (1.0/MAX(LCALC/SSPAR(1), RCALC/SSPAR(2), DCALC/SSPAR(3)))
 326      &       **(1.0/SSPAR(8)) * HOLD
 327       ENDIF
 328 C  HT  CONTAINS THE ESTIMATED OPTIMAL STEP SIZE.  NOW PUT IT WITHIN
 329 C REASONABLE BOUNDS.
 330       H=MIN(MAX(HT,SSPAR(6)*HOLD,SSPAR(4)), SSPAR(7)*HOLD, SSPAR(5))
 331       IF (ITNUM .EQ. 1) THEN
 332 C IF CONVERGENCE HAD OCCURRED AFTER 1 ITERATION, DON'T DECREASE  H .
 333         H=MAX(H,HOLD)
 334       ELSE IF (ITNUM .EQ. LITFH) THEN
 335 C IF CONVERGENCE REQUIRED THE MAXIMUM  LITFH  ITERATIONS, DON'T
 336 C INCREASE  H .
 337         H=MIN(H,HOLD)
 338       ENDIF
 339 C IF CONVERGENCE DID NOT OCCUR IN  LITFH  ITERATIONS FOR A PARTICULAR
 340 C H = HFAIL , DON'T CHOOSE THE NEW STEP SIZE LARGER THAN  HFAIL .
 341       IF (FAIL) H=MIN(H,HFAIL)
 342 C
 343 C
 344       RETURN
 345       END