SUBROUTINE ROOTNX(N,NFE,IFLAG,RELERR,ABSERR,Y,YP,YOLD,
& YPOLD,A,GOFW,TZ,W,WP)
C
C ROOTNX is an expert user version of ROOTN(F|S), written using the
C reverse call protocol. All matrix data structures and numerical linear
C algebra are the responsibility of the calling program. ROOTNX
C indicates to the calling program, via flags, at which points
C RHO(A,LAMBDA,X) and [ D RHO(A,LAMBDA,X)/D LAMBDA, D RHO(A,LAMBDA,X)/DX ]
C must be evaluated, and what linear algebra must be done with these
C functions. ROOTNX solves the following problem: given
C YOLD = (LAMBDA(S1),X(S1)), YPOLD = (LAMBDA'(S1),X'(S1)),
C Y = (LAMBDA(S2),X(S2)), YP = (LAMBDA'(S2),X'(S2)), S1 < S2, and a
C continuous function G(Y) = G(LAMBDA,X) such that G(YOLD) G(Y) <= 0,
C find the point Y(S) = (LAMBDA(S),X(S)), S1 <= S <= S2, such that
C G(Y(S)) = 0. ROOTNX alternates between secant estimates of Y(S)
C and Newton iteration until convergence.
C
C The following interface block should be inserted in the calling
C program:
C
C INTERFACE
C SUBROUTINE ROOTNX(N,NFE,IFLAG,RELERR,ABSERR,Y,YP,YOLD,
C & YPOLD,A,GOFW,TZ,W,WP)
C USE HOMOTOPY
C USE REAL_PRECISION
C INTEGER, INTENT(IN):: N
C INTEGER, INTENT(IN OUT):: NFE,IFLAG
C REAL (KIND=R8), INTENT(IN):: RELERR,ABSERR
C REAL (KIND=R8), DIMENSION(:), INTENT(IN):: A
C REAL (KIND=R8), INTENT(IN OUT):: GOFW
C REAL (KIND=R8), DIMENSION(:), INTENT(IN OUT):: Y,YP,YOLD,YPOLD,
C & TZ,W,WP
C END SUBROUTINE ROOTNX
C END INTERFACE
C
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 ROOTNX . ROOTNX 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,...,5, indicate to ROOTNX
C where to resume after a reverse call. The calling program
C must not modify IFLAG after a reverse call.
C
C RELERR, ABSERR = relative and absolute error values. The iteration is
C considered to have converged when a point Y=(LAMBDA,X) is found
C such that
C
C | G(Y(S)) | <= RELERR*||Y|| + ABSERR and
C
C ||Z|| <= RELERR*||Y|| + ABSERR , where
C
C Z is the Newton step to Y=(LAMBDA,X).
C
C Y(1:N+1) = point (LAMBDA(S), X(S)) on zero curve of 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 different from Y on the zero curve.
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 GOFW = G(W), the value requested by some reverse calls.
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 YPOLD > 0, ||WP|| = 1,
C and TZ is the minimum norm solution of
C D RHO(A,W)/DW TZ = - RHO(A,W).
C
C
C ON OUTPUT:
C
C N , RELERR , ABSERR , A are unchanged.
C
C NFE has been updated.
C
C IFLAG
C = -52, -51, or -50 requests the calling program to
C return the unit tangent vector in WP and the normal flow Newton
C step in TZ , all evaluated at the point W .
C
C = 0-10*R, -1-10*R, or -2-10*R, R = 1,...,4, requests the calling
C program to return in GOFW the scalar function G(Y) evaluated
C at W . The calling program must not modify IFLAG after a
C reverse call.
C
C = -2, -1, or 0 (unchanged) on a normal return.
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. Y and YOLD contain
C the last two points found on the zero curve.
C
C = 7 if input arguments or array sizes are invalid, or IFLAG was
C changed during a reverse call.
C
C Y is the point on the zero curve of the homotopy map such that
C G(Y) = 0.
C
C
C Calls DNRM2 , ROOT .
C
USE HOMOTOPY
USE REAL_PRECISION
INTEGER, INTENT(IN):: N
INTEGER, INTENT(IN OUT):: NFE,IFLAG
REAL (KIND=R8), INTENT(IN):: RELERR,ABSERR
REAL (KIND=R8), DIMENSION(:), INTENT(IN):: A
REAL (KIND=R8), INTENT(IN OUT):: GOFW
REAL (KIND=R8), DIMENSION(:), INTENT(IN OUT):: Y,YP,YOLD,YPOLD,
& TZ,W,WP
C
C ***** LOCAL VARIABLES. *****
C
REAL (KIND=R8), SAVE:: AERR,DELS,F0,F1,FP0,FP1,GOFY,
& GOFYOLD,GOFYP,RERR,S,SA,SB,SOUT,U
REAL (KIND=R8):: DD001,DD0011,DD01,DD011,DNRM2,QOFS
INTEGER, SAVE:: IFLAGC,JUDY,JW,LCODE,LIMIT,NP1
LOGICAL, SAVE:: BRACK
C
C ***** END OF SPECIFICATION INFORMATION. *****
C
C DEFINITION OF HERMITE CUBIC INTERPOLANT VIA DIVIDED DIFFERENCES.
C
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 ((MOD(-IFLAG,10) > 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.
& (IFLAG < -2 .AND. IFLAG /= IFLAGC)) 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 (100,110,130,210,200), ABS(IFLAG)/10
ENDIF
U=EPSILON(1.0_R8)
RERR=MAX(RELERR,U)
AERR=MAX(ABSERR,0.0_R8)
C
C THE LIMIT ON THE NUMBER OF ITERATIONS ALLOWED MAY BE CHANGED BY
C CHANGING THE FOLLOWING STATEMENT:
LIMIT=2*(INT(ABS(LOG10(AERR+RERR)))+1)
C
TZ=Y - YOLD
DELS=DNRM2(NP1,TZ,1)
C EVALUATE G AT YOLD AND Y .
W = YOLD
IFLAG = IFLAGC - 10
IFLAGC = IFLAG
RETURN
100 GOFYOLD = GOFW
W = Y
IFLAG = IFLAGC - 20
IFLAGC = IFLAG
RETURN
110 GOFY = GOFW
C
C USING TWO POINTS AND TANGENTS ON THE HOMOTOPY ZERO CURVE, CONSTRUCT
C THE HERMITE CUBIC INTERPOLANT Q(S). THEN USE ROOT TO FIND THE S
C CORRESPONDING TO G(Q(S)) = 0 . THE TWO POINTS ON THE ZERO CURVE ARE
C ALWAYS CHOSEN TO BRACKET G(Y(S)) = 0, WITH THE BRACKETING INTERVAL
C ALWAYS BEING [0, DELS].
C
SA=0.0
SB=DELS
LCODE=1
130 DO
CALL ROOT(SOUT,GOFW,SA,SB,RERR,AERR,LCODE)
IF (LCODE .GT. 0) EXIT
DO JW=1,NP1
W(JW)=QOFS(YOLD(JW),YPOLD(JW),Y(JW),YP(JW),DELS,SOUT)
END DO
C REQUEST VALUE G(Q(SOUT)) BY REVERSE CALL.
IFLAG = IFLAGC - 30
IFLAGC = IFLAG
RETURN
END DO
C IF G(Y(S)) = 0 WERE BRACKETED, ROOT CANNOT FAIL.
IF (LCODE .GT. 2) THEN
IFLAG=6
RETURN
ENDIF
C
C CALCULATE Q(SA) AS THE INITIAL POINT FOR A NEWTON ITERATION.
DO JW=1,NP1
W(JW)=QOFS(YOLD(JW),YPOLD(JW),Y(JW),YP(JW),DELS,SA)
END DO
C
C ***** END OF CALCULATION OF CUBIC INTERPOLANT. *****
C
C TANGENT INFORMATION YP IS NO LONGER NEEDED. HEREAFTER, YP
C REPRESENTS THE MOST RECENT POINT WHICH IS ON THE OPPOSITE SIDE OF
C THE SOLUTION TO G(Y(S)) = 0 FROM Y.
C
C PREPARE FOR MAIN LOOP.
C
YP=YOLD
GOFYP=GOFYOLD
C
C INITIALIZE BRACK TO INDICATE THAT THE POINTS Y AND YOLD BRACKET
C G(Y(S)) = 0, THUS YOLD = YP .
C
BRACK = .TRUE.
C
C ***** MAIN LOOP. *****
JUDY=1 ! DO JUDY = 1,LIMIT
170 IF (JUDY > LIMIT) GO TO 600
C CALCULATE NEWTON STEP TZ AT CURRENT ESTIMATE W .
IFLAG = IFLAGC - 50
IFLAGC = IFLAG
NFE = NFE+1
RETURN
C
C NEXT POINT = CURRENT POINT + NEWTON STEP.
C
200 W = W + TZ
C
C GET FUNCTION VALUE G(W) .
C
IFLAG = IFLAGC - 40
IFLAGC = IFLAG
RETURN
C
C CHECK FOR CONVERGENCE.
C
210 SA = RERR*DNRM2(NP1,W,1)+AERR
IF ((ABS(GOFW) .LE. SA) .AND. (DNRM2(NP1,TZ,1) .LE. SA)) THEN
Y = W
IFLAG = IFLAGC
RETURN
ENDIF
C
C PREPARE FOR NEXT ITERATION.
C
IF (ABS(GOFW) .LE. SA) THEN
YPOLD=WP
GO TO 590 ! CYCLE
ENDIF
C
C UPDATE Y AND YOLD .
C
YOLD=Y
Y=W
GOFYOLD=GOFY
GOFY=GOFW
C
C UPDATE YP SUCH THAT YP IS THE MOST RECENT POINT
C OPPOSITE OF G(Y(S)) = 0 FROM Y . SET BRACK = .TRUE. IFF
C Y AND YOLD BRACKET G(Y(S)) = 0 SO THAT YP = YOLD .
C
IF ( GOFY * GOFYOLD .GT. 0) THEN
BRACK = .FALSE.
ELSE
BRACK = .TRUE.
YP=YOLD
GOFYP=GOFYOLD
END IF
C
C COMPUTE DELS = ||Y-YP||.
C
TZ=Y - YP
DELS=DNRM2(NP1,TZ,1)
C
C COMPUTE TZ FOR THE LINEAR PREDICTOR W = Y + TZ,
C WHERE TZ = SA*(YOLD-Y).
C
S = ABS(GOFY - GOFYOLD)
IF (S .GE. 1.0) THEN
SA = GOFY/(GOFY - GOFYOLD)
TZ = SA*(YOLD - Y)
ELSE IF (ANY(ABS(GOFY*(YOLD-Y)) .GE. S*HUGE(1.0_R8))) THEN
TZ = DELS
ELSE
SA = GOFY/(GOFY - GOFYOLD)
TZ = SA*(YOLD - Y)
END IF
C
C TO INSURE STABILITY, THE LINEAR PREDICTION MUST BE NO FARTHER
C FROM Y THAN YP IS. THIS IS GUARANTEED IF BRACK = .TRUE.
C IF LINEAR PREDICTION IS TOO FAR AWAY, USE BRACKETING POINTS
C TO COMPUTE LINEAR PREDICTION.
C
IF (.NOT. BRACK) THEN
IF (DNRM2(NP1,TZ,1) .GT. DELS) THEN
C
C COMPUTE TZ = SA*(YP-Y).
C
SA = GOFY/(GOFY - GOFYP)
TZ = SA*(YP - Y)
END IF
END IF
C
C COMPUTE ESTIMATE W = Y + TZ AND SAVE OLD TANGENT VECTOR.
C
W = W + TZ
YPOLD = WP
590 JUDY=JUDY+1
GO TO 170 ! END DO
C
C ***** END OF MAIN LOOP. *****
C
C THE ALTERNATING SECANT ESTIMATION AND NEWTON ITERATION
C HAS NOT CONVERGED IN LIMIT STEPS. ERROR RETURN.
600 IFLAG=6
RETURN
END SUBROUTINE ROOTNX