SUBROUTINE ROOTQF(N,NFE,IFLAG,RELERR,ABSERR,Y,YP,YOLD,
     $     YPOLD,A,QT,R,DZ,Z,W,T,F0,F1,PAR,IPAR)
C
C ROOTQF  FINDS THE POINT  YBAR = (1, XBAR)  ON THE ZERO CURVE OF THE
C HOMOTOPY MAP.  IT STARTS WITH TWO POINTS  YOLD=(LAMBDAOLD,XOLD)  AND
C Y=(LAMBDA,X)  SUCH THAT  LAMBDAOLD < 1 <= LAMBDA, AND ALTERNATES
C BETWEEN USING A SECANT METHOD TO FIND A PREDICTED POINT ON THE
C HYPERPLANE  LAMBDA=1, AND TAKING A QUASI-NEWTON STEP TO RETURN TO THE
C ZERO CURVE OF THE HOMOTOPY MAP.
C
C
C ON INPUT:
C
C N = DIMENSION OF X.
C
C NFE = NUMBER OF JACOBIAN MATRIX EVALUATIONS.
C
C IFLAG = -2, -1, OR 0, INDICATING THE PROBLEM TYPE.
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    |Y(1) - 1| <= RELERR + ABSERR              AND
C
C    ||DZ|| <= RELERR*||Y|| + ABSERR,           WHERE
C
C    DZ  IS THE QUASI-NEWTON STEP TO Y.
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(1:*) = PARAMETER VECTOR IN THE HOMOTOPY MAP.
C
C QT(1:N+1,1:N+1) CONTAINS  Q  TRANSPOSE OF THE QR FACTORIZATION OF
C    THE AUGMENTED JACOBIAN MATRIX EVALUATED AT THE POINT Y.
C
C R((N+1)*(N+2)/2) CONTAINS THE UPPER TRIANGLE OF THE R PART OF
C    OF THE QR FACTORIZATION, STORED BY ROWS.
C
C DZ(1:N+1), Z(1:N+1), W(1:N+1), T(1:N+1), F0(1:N+1), F1(1:N+1)
C    ARE WORK ARRAYS USED FOR THE QUASI-NEWTON STEP AND THE SECANT
C    STEP.
C
C PAR(1:*)  AND  IPAR(1:*)  ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS,
C    WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES
C    RHO, RHOJAC.
C
C
C ON OUTPUT:
C
C N, RELERR, ABSERR, AND A  ARE UNCHANGED.
C
C NFE  HAS BEEN UPDATED.
C
C IFLAG
C    = -2, -1, OR 0 (UNCHANGED) ON A NORMAL RETURN.
C
C    = 4 IF A SINGULAR JACOBIAN MATRIX 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 OBTAINED BY QUASI-NEWTON STEPS, AND  YP
C        CONTAINS A POINT OPPOSITE OF THE HYPERPLANE  LAMBDA=1  FROM
C        Y.
C
C Y  IS THE POINT ON THE ZERO CURVE OF THE HOMOTOPY MAP AT  LAMBDA = 1.
C
C YP  AND  YOLD  CONTAIN POINTS NEAR THE SOLUTION.
C
C CALLS  D1MACH, DAXPY, DCOPY, DDOT, DNRM2, F (OR RHO),
C        QRSLQF, ROOT, UPQRQF.
C
C ***** DECLARATIONS *****
C
C     FUNCTION DECLARATIONS
C
        DOUBLE PRECISION D1MACH, DDOT, DNRM2, QOFS
C
C     LOCAL VARIABLES
C
        DOUBLE PRECISION AERR, DD001, DD0011, DD01, DD011, DELS, ETA,
     $     ONE, P0, P1, PP0, PP1, QSOUT, RERR, S, SA, SB, SOUT,
     $     U, ZERO
        INTEGER ISTEP, I, LCODE, LIMIT,NP1
        LOGICAL BRACK
C
C     SCALAR ARGUMENTS
C
        DOUBLE PRECISION RELERR, ABSERR
        INTEGER N, NFE, IFLAG
C
C     ARRAY DECLARATIONS
C
        DOUBLE PRECISION Y(N+1), YP(N+1), YOLD(N+1), YPOLD(N+1), A(N),
     $     QT(N+1:N+1), R((N+1)*(N+2)/2), DZ(N+1), Z(N+1), W(N+1),
     $     T(N+1), F0(N+1), F1(N+1), PAR(1)
        INTEGER IPAR(1)
C
C ***** END OF DECLARATIONS *****
C
C
C DEFINITION OF HERMITE CUBIC INTERPOLANT VIA DIVIDED DIFFERENCES.
C
        DD01(P0,P1,DELS)=(P1-P0)/DELS
        DD001(P0,PP0,P1,DELS)=(DD01(P0,P1,DELS)-PP0)/DELS
        DD011(P0,P1,PP1,DELS)=(PP1-DD01(P0,P1,DELS))/DELS
        DD0011(P0,PP0,P1,PP1,DELS)=(DD011(P0,P1,PP1,DELS) -
     $                              DD001(P0,PP0,P1,DELS))/DELS
        QOFS(P0,PP0,P1,PP1,DELS,S)=((DD0011(P0,PP0,P1,PP1,DELS)*
     $     (S-DELS) + DD001(P0,PP0,P1,DELS))*S + PP0)*S + P0
C
C ***** FIRST EXECUTABLE STATEMENT *****
C
C ***** INITIALIZATION *****
C
C ETA = PARAMETER FOR BROYDEN'S UPDATE.
C LIMIT = MAXIMUM NUMBER OF ITERATIONS ALLOWED.
C
        ONE=1.0
        ZERO=0.0
        U=D1MACH(4)
        RERR=MAX(RELERR,U)
        AERR=MAX(ABSERR,ZERO)
        NP1=N+1
        ETA = 100.0*U
        LIMIT = 2*(INT(-LOG10(AERR+RERR*DNRM2(NP1,Y,1)))+1)
C
C F0 = (RHO(Y), YP*Y) TRANSPOSE.
C
        IF (IFLAG .EQ. -2) THEN
C
C         CURVE TRACKING PROBLEM.
C
          CALL RHO(A,Y(1),Y(2),F0,PAR,IPAR)
        ELSE IF (IFLAG .EQ. -1) THEN
C
C         ZERO FINDING PROBLEM.
C
          CALL F(Y(2),F0)
          DO 10 I=1,N
            F0(I) = Y(1)*F0(I) + (1.0-Y(1))*(Y(I+1)-A(I))
10        CONTINUE
        ELSE
C
C         FIXED POINT PROBLEM.
C
          CALL F(Y(2),F0)
          DO 20 I=1,N
            F0(I) = Y(1)*(A(I)-F0(I))+Y(I+1)-A(I)
20        CONTINUE
        END IF
        F0(NP1) = DDOT(NP1,YP,1,Y,1)
C
C ***** END OF INITIALIZATION BLOCK *****
C
C ***** COMPUTE FIRST INTERPOLANT WITH A HERMITE CUBIC *****
C
C FIND DISTANCE BETWEEN Y AND YOLD.  DZ=||Y-YOLD||.
C
        CALL DCOPY(NP1,Y,1,DZ,1)
        CALL DAXPY(NP1,-ONE,YOLD,1,DZ,1)
        DELS=DNRM2(NP1,DZ,1)
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  LAMBDA = 1.  THE TWO POINTS ON THE ZERO CURVE ARE
C ALWAYS CHOSEN TO BRACKET  LAMBDA=1, WITH THE BRACKETING INTERVAL
C ALWAYS BEING [0, DELS].
C
        SA=0.0
        SB=DELS
        LCODE=1
40      CALL ROOT(SOUT,QSOUT,SA,SB,RERR,AERR,LCODE)
          IF (LCODE .GT. 0) GO TO 50
          QSOUT=QOFS(YOLD(1),YPOLD(1),Y(1),YP(1),DELS,SOUT) - 1.0
          GO TO 40
C
C     IF  LAMBDA = 1  WERE BRACKETED,  ROOT  CANNOT FAIL.
C
50      IF (LCODE .GT. 2) THEN
          IFLAG=6
          RETURN
        ENDIF
C
C CALCULATE  Q(SA)  AS THE INITIAL POINT FOR A NEWTON ITERATION.
C
        DO 60 I=1,NP1
          Z(I)=QOFS(YOLD(I),YPOLD(I),Y(I),YP(I),DELS,SA)
60      CONTINUE
C
C CALCULATE DZ = Z-Y.
C
        CALL DCOPY(NP1,Z,1,DZ,1)
        CALL DAXPY(NP1,-ONE,Y,1,DZ,1)
C
C ***** END OF CALCULATION OF CUBIC INTERPOLANT *****
C
C TANGENT INFORMATION  YPOLD  IS NO LONGER NEEDED.  HEREAFTER,  YPOLD
C REPRESENTS THE MOST RECENT POINT WHICH IS ON THE OPPOSITE SIDE OF
C LAMBDA=1  FROM  Y.
C
C ***** PREPARE FOR MAIN LOOP *****
C
        CALL DCOPY(NP1,YOLD,1,YPOLD,1)
C
C INITIALIZE  BRACK  TO INDICATE THAT THE POINTS  Y  AND  YOLD  BRACKET
C LAMBDA=1,  THUS YOLD = YPOLD.
C
        BRACK = .TRUE.
C
C ***** MAIN LOOP *****
C
        DO 300 ISTEP=1,LIMIT
C
C UPDATE JACOBIAN MATRIX.
C
C       F1=(RHO(Z), YP*Z) TRANSPOSE.
C
          IF (IFLAG .EQ. -2) THEN
            CALL RHO(A,Z(1),Z(2),F1,PAR,IPAR)
          ELSE IF (IFLAG .EQ. -1) THEN
            CALL F(Z(2),F1)
            DO 80 I=1,N
              F1(I) = Z(1)*F1(I) + (1-Z(1))*(Z(I+1)-A(I))
  80        CONTINUE
          ELSE
            CALL F(Z(2),F1)
            DO 90 I=1,N
              F1(I) = Z(1)*(A(I)-F1(I))+Z(I+1)-A(I)
  90        CONTINUE
          END IF
          F1(NP1) = DDOT(NP1,YP,1,Z,1)
C
C
C PERFORM BROYDEN UPDATE.
C
          CALL UPQRQF(NP1,ETA,DZ,F0,F1,QT,R,W,T)
C
C QUASI-NEWTON STEP.
C
C       COMPUTE NEWTON STEP.
C
          CALL DCOPY(N,F1,1,DZ,1)
          CALL DSCAL(N,-ONE,DZ,1)
          DZ(NP1) = 0.0
          CALL QRSLQF(QT,R,DZ,W,NP1)
C
C       TAKE NEWTON STEP.
C
          CALL DCOPY(NP1,Z,1,W,1)
          CALL DAXPY(NP1,ONE,DZ,1,Z,1)
C
C       CHECK FOR CONVERGENCE.
C
          IF ((ABS(Z(1)-1.0) .LE. RERR+AERR) .AND.
     $        (DNRM2(NP1,DZ,1) .LE. RERR*DNRM2(N,Z(2),1)+AERR)) THEN
             CALL DCOPY(NP1,Z,1,Y,1)
             RETURN
          END IF
C
C PREPARE FOR NEXT ITERATION.
C
C       F0 = F1.
C
          CALL DCOPY(NP1,F1,1,F0,1)
C
C       IF  Z(1) = 1.0  THEN PERFORM QUASI-NEWTON ITERATION AGAIN
C       WITHOUT COMPUTING A NEW PREDICTOR.
C
          IF (ABS(Z(1)-1.0) .LE. RERR+AERR) THEN
             CALL DCOPY(NP1,Z,1,DZ,1)
             CALL DAXPY(NP1,-ONE,W,1,DZ,1)
             GOTO 300
          END IF
C
C       UPDATE  Y  AND  YOLD.
C
          CALL DCOPY(NP1,Y,1,YOLD,1)
          CALL DCOPY(NP1,Z,1,Y,1)
C
C       UPDATE  YPOLD  SUCH THAT  YPOLD  IS THE MOST RECENT POINT
C       OPPOSITE OF  LAMBDA=1  FROM  Y.  SET  BRACK = .TRUE.  IFF
C       Y & YOLD  BRACKET  LAMBDA=1  SO THAT  YPOLD=YOLD.
C
          IF ((Y(1)-1.0)*(YOLD(1)-1.0) .GT. 0) THEN
            BRACK = .FALSE.
          ELSE
            BRACK = .TRUE.
            CALL DCOPY(NP1,YOLD,1,YPOLD,1)
          END IF
C
C       COMPUTE DELS = ||Y-YPOLD||.
C
          CALL DCOPY(NP1,Y,1,DZ,1)
          CALL DAXPY(NP1,-ONE,YPOLD,1,DZ,1)
          DELS=DNRM2(NP1,DZ,1)
C
C       COMPUTE  DZ  FOR THE LINEAR PREDICTOR   Z = Y + DZ,
C           WHERE  DZ = SA*(YOLD-Y).
C
          SA = (1.0-Y(1))/(YOLD(1)-Y(1))
          CALL DCOPY(NP1,YOLD,1,DZ,1)
          CALL DAXPY(NP1,-ONE,Y,1,DZ,1)
          CALL DSCAL(NP1,SA,DZ,1)
C
C       TO INSURE STABILITY, THE LINEAR PREDICTION MUST BE NO FARTHER
C       FROM  Y  THAN  YPOLD  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,DZ,1) .GT. DELS) THEN
C
C             COMPUTE  DZ = SA*(YPOLD-Y).
C
              SA = (1.0-Y(1))/(YPOLD(1)-Y(1))
              CALL DCOPY(NP1,YPOLD,1,DZ,1)
              CALL DAXPY(NP1,-ONE,Y,1,DZ,1)
              CALL DSCAL(NP1,SA,DZ,1)
            END IF
          END IF
C
C       COMPUTE PREDICTOR Z = Y+DZ, AND DZ = NEW Z  - OLD Z (USED FOR
C         QUASI-NEWTON UPDATE).
C
          CALL DAXPY(NP1,ONE,DZ,1,Z,1)
          CALL DCOPY(NP1,Z,1,DZ,1)
          CALL DAXPY(NP1,-ONE,W,1,DZ,1)
  300   CONTINUE
C
C ***** END OF MAIN LOOP. *****
C
C THE ALTERNATING OSCULATORY LINEAR PREDICTION AND QUASI-NEWTON
C CORRECTION HAS NOT CONVERGED IN  LIMIT  STEPS.  ERROR RETURN.
        IFLAG=6
        RETURN
C
C ***** END OF SUBROUTINE ROOTQF *****
        END