share/hompack/fortran/stepns.f

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