SUBROUTINE STEPNX(N,NFE,IFLAG,START,CRASH,HOLD,H,RELERR,
& ABSERR,S,Y,YP,YOLD,YPOLD,A,TZ,W,WP,RHOLEN,SSPAR)
C
C STEPNX takes one step along the zero curve of the homotopy map
C using a predictor-corrector algorithm. The predictor uses a Hermite
C cubic interpolant, and the corrector returns to the zero curve along
C the flow normal to the Davidenko flow. STEPNX also estimates a
C step size H for the next step along the zero curve. STEPNX is an
C expert user version of STEPN(F|S), written using the reverse call
C protocol. All matrix data structures and numerical linear algebra
C are the responsibility of the calling program. STEPNX indicates to
C the calling program, via flags, at which points RHO(A,LAMBDA,X) and
C [ D RHO(A,LAMBDA,X)/D LAMBDA, D RHO(A,LAMBDA,X)/DX ] must be
C evaluated, and what linear algebra must be done with these functions.
C Out of range arguments can also be signaled to STEPNX , which will
C attempt to modify its steplength algorithm to reflect this
C information.
C
C The following interface block should be inserted in the calling
C program:
C
C INTERFACE
C SUBROUTINE STEPNX(N,NFE,IFLAG,START,CRASH,HOLD,H,RELERR,
C & ABSERR,S,Y,YP,YOLD,YPOLD,A,TZ,W,WP,RHOLEN,SSPAR)
C USE HOMOTOPY
C USE REAL_PRECISION
C INTEGER, INTENT(IN):: N
C INTEGER, INTENT(IN OUT):: NFE,IFLAG
C LOGICAL, INTENT(IN OUT):: START,CRASH
C REAL (KIND=R8), INTENT(IN OUT):: HOLD,H,RELERR,ABSERR,S,RHOLEN,
C & SSPAR(8)
C REAL (KIND=R8), DIMENSION(:), INTENT(IN):: A
C REAL (KIND=R8), DIMENSION(:), INTENT(IN OUT):: Y,YP,YOLD,YPOLD,
C & TZ,W,WP
C REAL (KIND=R8), DIMENSION(:), ALLOCATABLE, SAVE:: Z0,Z1
C END SUBROUTINE STEPNX
C END INTERFACE
C
C ON INPUT:
C
C N = dimension of X and the homotopy map.
C
C NFE = number of Jacobian matrix evaluations.
C
C IFLAG = -2, -1, or 0, indicating the problem type, on the first
C call to STEPNX . STEPNX does not distinguish between
C these values, but they are permitted for consistency with
C the rest of HOMPACK.
C
C = 0-10*R, -1-10*R, or -2-10*R, R = 1,2,3, indicate to STEPNX
C where to resume after a reverse call. The calling program
C must not modify IFLAG after a reverse call, except as
C noted next.
C
C = -40, -41, or -42, used for a final call to deallocate working
C storage, after all path tracking is finished. START and
C IFLAG are reset on return.
C
C = -100-10*R, -101-10*R, -102-10*R, R = 1,2,3, indicate to
C STEPNX where to resume after a reverse call, and that the
C requested evaluation point was out of range. STEPNX will
C reduce H and try again.
C
C START = .TRUE. on first call to STEPNX , .FALSE. otherwise.
C
C HOLD = ||Y - YOLD||; should not be modified by the user.
C
C H = upper limit on length of step that will be attempted. H must be
C set to a positive number on the first call to STEPNX .
C Thereafter STEPNX calculates an optimal value for H , and H
C should not be modified by the user.
C
C RELERR, ABSERR = relative and absolute error values. The iteration is
C considered to have converged when a point W=(LAMBDA,X) is found
C such that
C
C ||Z|| <= RELERR*||W|| + ABSERR , where
C
C Z is the Newton step to W=(LAMBDA,X).
C
C S = (approximate) arc length along the homotopy zero curve up to
C Y(S) = (LAMBDA(S), X(S)).
C
C Y(1:N+1) = previous point (LAMBDA(S), X(S)) found on the zero curve of
C the homotopy map.
C
C YP(1:N+1) = unit tangent vector to the zero curve of the homotopy map
C at Y .
C
C YOLD(1:N+1) = a point before Y on the zero curve of the homotopy map.
C
C YPOLD(1:N+1) = unit tangent vector to the zero curve of the homotopy
C map at YOLD .
C
C A(:) = parameter vector in the homotopy map.
C
C TZ(1:N+1), W(1:N+1), and WP(1:N+1) are work arrays used for the
C Newton step calculation and the interpolation. On reentry after
C a reverse call, WP and TZ contain the tangent vector and
C Newton step, respectively, at the point W . Precisely,
C D RHO(A,W)/DW WP = 0, WP^T YP > 0, ||WP|| = 1,
C and TZ is the minimum norm solution of
C D RHO(A,W)/DW TZ = - RHO(A,W).
C
C RHOLEN = ||RHO(A,W)||_2 is required by some reverse calls.
C
C SSPAR(1:8) = (LIDEAL, RIDEAL, DIDEAL, HMIN, HMAX, BMIN, BMAX, P) is
C a vector of parameters used for the optimal step size estimation.
C If SSPAR(J) .LE. 0.0 on input, it is reset to a default value
C by STEPNX . Otherwise the input value of SSPAR(J) is used.
C See the comments below in STEPNX for more information about
C these constants.
C
C
C ON OUTPUT:
C
C N and A are unchanged.
C
C NFE has been updated.
C
C IFLAG
C = -22, -21, -20, -32, -31, or -30 requests the calling program to
C return the unit tangent vector in WP , the normal flow Newton
C step in TZ , and the 2-norm of the homotopy map in RHOLEN ,
C all evaluated at the point W .
C
C = -12, -11, or -10 requests the calling program to return in WP
C the unit tangent vector at W .
C
C = -2, -1, or 0 (unchanged) on a normal return after a successful
C step.
C
C = 4 if a Jacobian matrix with rank < N has occurred. The
C iteration was not completed.
C
C = 6 if the iteration failed to converge. W contains the last
C Newton iterate.
C
C = 7 if input arguments or array sizes are invalid, or IFLAG was
C changed during a reverse call.
C
C START = .FALSE. on a normal return.
C
C CRASH
C = .FALSE. on a normal return.
C
C = .TRUE. if the step size H was too small. H has been
C increased to an acceptable value, with which STEPNX may be
C called again.
C
C = .TRUE. if RELERR and/or ABSERR were too small. They have
C been increased to acceptable values, with which STEPNX may
C be called again.
C
C HOLD = ||Y - YOLD||.
C
C H = optimal value for next step to be attempted. Normally H should
C not be modified by the user.
C
C RELERR, ABSERR are unchanged on a normal return.
C
C S = (approximate) arc length along the zero curve of the homotopy map
C up to the latest point found, which is returned in Y .
C
C Y, YP, YOLD, YPOLD contain the two most recent points and tangent
C vectors found on the zero curve of the homotopy map.
C
C SSPAR may have been changed to default values.
C
C
C Z0(1:N+1), Z1(1:N+1) are allocatable work arrays used for the
C estimation of the next step size H .
C
C Calls DNRM2 .
C
USE HOMOTOPY
USE REAL_PRECISION
INTEGER, INTENT(IN):: N
INTEGER, INTENT(IN OUT):: NFE,IFLAG
LOGICAL, INTENT(IN OUT):: START,CRASH
REAL (KIND=R8), INTENT(IN OUT):: HOLD,H,RELERR,ABSERR,S,RHOLEN,
& SSPAR(8)
REAL (KIND=R8), DIMENSION(:), INTENT(IN):: A
REAL (KIND=R8), DIMENSION(:), INTENT(IN OUT):: Y,YP,YOLD,YPOLD,
& TZ,W,WP
REAL (KIND=R8), DIMENSION(:), ALLOCATABLE, SAVE:: Z0,Z1
C
C ***** LOCAL VARIABLES. *****
C
REAL (KIND=R8), SAVE:: DCALC,DELS,F0,F1,FOURU,FP0,FP1,
& HFAIL,HT,LCALC,RCALC,TEMP,TWOU
INTEGER, SAVE:: IFLAGC,ITNUM,J,JUDY,NP1
LOGICAL, SAVE:: FAIL
C
C ***** END OF SPECIFICATION INFORMATION. *****
C
C THE LIMIT ON THE NUMBER OF NEWTON ITERATIONS ALLOWED BEFORE REDUCING
C THE STEP SIZE H MAY BE CHANGED BY CHANGING THE FOLLOWING PARAMETER
C STATEMENT:
INTEGER, PARAMETER:: LITFH=4
C
C DEFINITION OF HERMITE CUBIC INTERPOLANT VIA DIVIDED DIFFERENCES.
C
REAL (KIND=R8):: DD001,DD0011,DD01,DD011,DNRM2,QOFS
DD01(F0,F1,DELS)=(F1-F0)/DELS
DD001(F0,FP0,F1,DELS)=(DD01(F0,F1,DELS)-FP0)/DELS
DD011(F0,F1,FP1,DELS)=(FP1-DD01(F0,F1,DELS))/DELS
DD0011(F0,FP0,F1,FP1,DELS)=(DD011(F0,F1,FP1,DELS) -
& DD001(F0,FP0,F1,DELS))/DELS
QOFS(F0,FP0,F1,FP1,DELS,S)=((DD0011(F0,FP0,F1,FP1,DELS)*(S-DELS) +
& DD001(F0,FP0,F1,DELS))*S + FP0)*S + F0
C
C
NP1=N+1
IF (IFLAG > 0) RETURN
IF ((START .AND. IFLAG < -2) .OR. SIZE(Y) /= NP1 .OR.
& SIZE(YP) /= NP1 .OR. SIZE(YOLD) /= NP1 .OR.
& SIZE(YPOLD) /= NP1 .OR. SIZE(TZ) /= NP1 .OR.
& SIZE(W) /= NP1 .OR. SIZE(WP) /= NP1 .OR.
& (.NOT. START .AND. -MOD(-IFLAG,100) /= IFLAGC .AND.
& ABS(IFLAG)/10 /= 4)) THEN
IFLAG=7
RETURN
ENDIF
IFLAGC=-MOD(-IFLAG,10)
C
C PICK UP EXECUTION WEHRE IT LEFT OFF AFTER A REVERSE CALL.
C
IF (IFLAG < -2) THEN
GO TO (50,100,400,700), MOD(ABS(IFLAG),100)/10
ENDIF
TWOU=2.0*EPSILON(1.0_R8)
FOURU=TWOU+TWOU
CRASH=.TRUE.
C THE ARCLENGTH S MUST BE NONNEGATIVE.
IF (S .LT. 0.0) RETURN
C IF STEP SIZE IS TOO SMALL, DETERMINE AN ACCEPTABLE ONE.
IF (H .LT. FOURU*(1.0+S)) THEN
H=FOURU*(1.0+S)
RETURN
ENDIF
C IF ERROR TOLERANCES ARE TOO SMALL, INCREASE THEM TO ACCEPTABLE VALUES.
TEMP=DNRM2(NP1,Y,1)+1.0
IF (.5*(RELERR*TEMP+ABSERR) .LT. TWOU*TEMP) THEN
IF (RELERR .NE. 0.0) THEN
RELERR=FOURU*(1.0+FOURU)
ABSERR=MAX(ABSERR,0.0_R8)
ELSE
ABSERR=FOURU*TEMP
ENDIF
RETURN
ENDIF
CRASH=.FALSE.
IF (.NOT. START) GO TO 300
C
C ***** STARTUP SECTION (FIRST STEP ALONG ZERO CURVE). *****
C
FAIL=.FALSE.
START=.FALSE.
IF (ALLOCATED(Z0)) DEALLOCATE(Z0)
IF (ALLOCATED(Z1)) DEALLOCATE(Z1)
ALLOCATE(Z0(NP1),Z1(NP1))
C
C SET OPTIMAL STEP SIZE ESTIMATION PARAMETERS.
C LET Z[K] DENOTE THE NEWTON ITERATES ALONG THE FLOW NORMAL TO THE
C DAVIDENKO FLOW AND Y THEIR LIMIT.
C IDEAL CONTRACTION FACTOR: ||Z[2] - Z[1]|| / ||Z[1] - Z[0]||
IF (SSPAR(1) .LE. 0.0) SSPAR(1)= .5
C IDEAL RESIDUAL FACTOR: ||RHO(A, Z[1])|| / ||RHO(A, Z[0])||
IF (SSPAR(2) .LE. 0.0) SSPAR(2)= .01
C IDEAL DISTANCE FACTOR: ||Z[1] - Y|| / ||Z[0] - Y||
IF (SSPAR(3) .LE. 0.0) SSPAR(3)= .5
C MINIMUM STEP SIZE HMIN .
IF (SSPAR(4) .LE. 0.0) SSPAR(4)=(SQRT(N+1.0)+4.0)*EPSILON(1.0_R8)
C MAXIMUM STEP SIZE HMAX .
IF (SSPAR(5) .LE. 0.0) SSPAR(5)= 1.0
C MINIMUM STEP SIZE REDUCTION FACTOR BMIN .
IF (SSPAR(6) .LE. 0.0) SSPAR(6)= .1_R8
C MAXIMUM STEP SIZE EXPANSION FACTOR BMAX .
IF (SSPAR(7) .LE. 0.0) SSPAR(7)= 3.0
C ASSUMED OPERATING ORDER P .
IF (SSPAR(8) .LE. 0.0) SSPAR(8)= 2.0
C
C DETERMINE SUITABLE INITIAL STEP SIZE.
H=MIN(H, .10_R8, SQRT(SQRT(RELERR*TEMP+ABSERR)))
C USE LINEAR PREDICTOR ALONG TANGENT DIRECTION TO START NEWTON ITERATION.
YPOLD(1)=1.0
YPOLD(2:NP1)=0.0
C REQUEST TANGENT VECTOR AT Y VIA REVERSE CALL.
W=Y
YP=YPOLD
IFLAG=IFLAGC-10
IFLAGC=IFLAG
NFE=NFE+1
RETURN
50 YP=WP
C IF THE STARTING POINT IS OUT OF RANGE, GIVE UP.
IF (IFLAG .LE. -100) THEN
IFLAG=6
RETURN
ENDIF
70 W=Y + H*YP
Z0=W
JUDY=1 ! DO JUDY=1,LITFH
80 IF (JUDY > LITFH) GO TO 200
C REQUEST THE CALCULATION OF THE NEWTON STEP TZ AT THE CURRENT
C POINT W VIA REVERSE CALL.
IFLAG=IFLAGC-20
IFLAGC=IFLAG
NFE=NFE+1
RETURN
100 IF (IFLAG .LE. -100) GO TO 200
C
C TAKE NEWTON STEP AND CHECK CONVERGENCE.
W=W + TZ
ITNUM=JUDY
C COMPUTE QUANTITIES USED FOR OPTIMAL STEP SIZE ESTIMATION.
IF (JUDY .EQ. 1) THEN
LCALC=DNRM2(NP1,TZ,1)
RCALC=RHOLEN
Z1=W
ELSE IF (JUDY .EQ. 2) THEN
LCALC=DNRM2(NP1,TZ,1)/LCALC
RCALC=RHOLEN/RCALC
ENDIF
C GO TO MOP-UP SECTION AFTER CONVERGENCE.
IF (DNRM2(NP1,TZ,1) .LE. RELERR*DNRM2(NP1,W,1)+ABSERR)
& GO TO 600
C
JUDY=JUDY+1
GO TO 80 ! END DO
C
C NO CONVERGENCE IN LITFH ITERATIONS. REDUCE H AND TRY AGAIN.
200 IF (H .LE. FOURU*(1.0 + S)) THEN
IFLAG=6
RETURN
ENDIF
H=.5 * H
GO TO 70
C
C ***** END OF STARTUP SECTION. *****
C
C ***** PREDICTOR SECTION. *****
C
300 FAIL=.FALSE.
C COMPUTE POINT PREDICTED BY HERMITE INTERPOLANT. USE STEP SIZE H
C COMPUTED ON LAST CALL TO STEPNX .
320 DO J=1,NP1
W(J)=QOFS(YOLD(J),YPOLD(J),Y(J),YP(J),HOLD,HOLD+H)
END DO
Z0=W
C
C ***** END OF PREDICTOR SECTION. *****
C
C ***** CORRECTOR SECTION. *****
C
JUDY=1 ! CORRECTOR: DO JUDY=1,LITFH
350 IF (JUDY > LITFH) GO TO 500
C REQUEST THE CALCULATION OF THE NEWTON STEP TZ AT THE CURRENT
C POINT W VIA REVERSE CALL.
IFLAG=IFLAGC-30
IFLAGC=IFLAG
NFE=NFE+1
RETURN
400 IF (IFLAG .LE. -100) GO TO 500
C
C TAKE NEWTON STEP AND CHECK CONVERGENCE.
W=W + TZ
ITNUM=JUDY
C COMPUTE QUANTITIES USED FOR OPTIMAL STEP SIZE ESTIMATION.
IF (JUDY .EQ. 1) THEN
LCALC=DNRM2(NP1,TZ,1)
RCALC=RHOLEN
Z1=W
ELSE IF (JUDY .EQ. 2) THEN
LCALC=DNRM2(NP1,TZ,1)/LCALC
RCALC=RHOLEN/RCALC
ENDIF
C GO TO MOP-UP SECTION AFTER CONVERGENCE.
IF (DNRM2(NP1,TZ,1) .LE. RELERR*DNRM2(NP1,W,1)+ABSERR)
& GO TO 600
C
JUDY=JUDY+1
GO TO 350 ! END DO CORRECTOR
C
C NO CONVERGENCE IN LITFH ITERATIONS. RECORD FAILURE AT CALCULATED H ,
C SAVE THIS STEP SIZE, REDUCE H AND TRY AGAIN.
500 FAIL=.TRUE.
HFAIL=H
IF (H .LE. FOURU*(1.0 + S)) THEN
IFLAG=6
RETURN
ENDIF
H=.5 * H
GO TO 320
C
C ***** END OF CORRECTOR SECTION. *****
C
C ***** MOP-UP SECTION. *****
C
C YOLD AND Y ALWAYS CONTAIN THE LAST TWO POINTS FOUND ON THE ZERO
C CURVE OF THE HOMOTOPY MAP. YPOLD AND YP CONTAIN THE TANGENT
C VECTORS TO THE ZERO CURVE AT YOLD AND Y , RESPECTIVELY.
C
600 YPOLD=YP
YOLD=Y
Y=W
YP=WP
W=Y - YOLD
C UPDATE ARC LENGTH.
HOLD=DNRM2(NP1,W,1)
S=S+HOLD
C
C ***** END OF MOP-UP SECTION. *****
C
C ***** OPTIMAL STEP SIZE ESTIMATION SECTION. *****
C
C CALCULATE THE DISTANCE FACTOR DCALC .
TZ=Z0 - Y
W=Z1 - Y
DCALC=DNRM2(NP1,TZ,1)
IF (DCALC .NE. 0.0) DCALC=DNRM2(NP1,W,1)/DCALC
C
C THE OPTIMAL STEP SIZE HBAR IS DEFINED BY
C
C HT=HOLD * [MIN(LIDEAL/LCALC, RIDEAL/RCALC, DIDEAL/DCALC)]**(1/P)
C
C HBAR = MIN [ MAX(HT, BMIN*HOLD, HMIN), BMAX*HOLD, HMAX ]
C
C IF CONVERGENCE HAD OCCURRED AFTER 1 ITERATION, SET THE CONTRACTION
C FACTOR LCALC TO ZERO.
IF (ITNUM .EQ. 1) LCALC = 0.0
C FORMULA FOR OPTIMAL STEP SIZE.
IF (LCALC+RCALC+DCALC .EQ. 0.0) THEN
HT = SSPAR(7) * HOLD
ELSE
HT = (1.0/MAX(LCALC/SSPAR(1), RCALC/SSPAR(2), DCALC/SSPAR(3)))
& **(1.0/SSPAR(8)) * HOLD
ENDIF
C HT CONTAINS THE ESTIMATED OPTIMAL STEP SIZE. NOW PUT IT WITHIN
C REASONABLE BOUNDS.
H=MIN(MAX(HT,SSPAR(6)*HOLD,SSPAR(4)), SSPAR(7)*HOLD, SSPAR(5))
IF (ITNUM .EQ. 1) THEN
C IF CONVERGENCE HAD OCCURRED AFTER 1 ITERATION, DON'T DECREASE H .
H=MAX(H,HOLD)
ELSE IF (ITNUM .EQ. LITFH) THEN
C IF CONVERGENCE REQUIRED THE MAXIMUM LITFH ITERATIONS, DON'T
C INCREASE H .
H=MIN(H,HOLD)
ENDIF
C IF CONVERGENCE DID NOT OCCUR IN LITFH ITERATIONS FOR A PARTICULAR
C H = HFAIL , DON'T CHOOSE THE NEW STEP SIZE LARGER THAN HFAIL .
IF (FAIL) H=MIN(H,HFAIL)
C
C
IFLAG=IFLAGC
RETURN
C CLEAN UP ALLOCATED WORKING STORAGE.
700 START=.TRUE.
IFLAG=IFLAGC
IF (ALLOCATED(Z0)) DEALLOCATE(Z0)
IF (ALLOCATED(Z1)) DEALLOCATE(Z1)
RETURN
END SUBROUTINE STEPNX