SUBROUTINE ROOTQS(N,NFE,IFLAG,LENQR,RELERR,ABSERR,Y,YP,YOLD,
     &     YPOLD,A,QR,PIVOT,PP,RHOVEC,Z,DZ,WORK,PAR,IPAR)
C
C ROOTQS  FINDS THE POINT  YBAR = (XBAR, 1)  ON THE ZERO CURVE OF THE
C HOMOTOPY MAP.  IT STARTS WITH TWO POINTS  YOLD=(XOLD,LAMBDAOLD)  AND
C Y=(X,LAMBDA)  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 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 LENQR = THE LENGTH OF THE ONE-DIMENSIONAL ARRAY  QR.
C
C RELERR, ABSERR = RELATIVE AND ABSOLUTE ERROR VALUES.  THE ITERATION IS
C    CONSIDERED TO HAVE CONVERGED WHEN A POINT Y=(X,LAMBDA) IS FOUND 
C    SUCH THAT
C
C    |Y(N+1) - 1| <= RELERR + ABSERR              AND
C
C    ||DZ|| <= RELERR*||Y|| + ABSERR,           WHERE
C
C    DZ  IS THE NEWTON STEP TO Y.
C
C Y(1:N+1) = POINT (X(S),LAMBDA(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 QR(1:LENQR)  IS A WORK ARRAY CONTAINING THE N X N SYMMETRIC
C    JACOBIAN MATRIX WITH RESPECT TO X STORED IN PACKED SKYLINE
C    STORAGE FORMAT.  LENQR  AND  PIVOT  DESCRIBE THE DATA
C    STRUCTURE IN  QR.  (SEE SUBROUTINE  PCGQS  FOR A DESCRIPTION
C    OF THIS DATA STRUCTURE).
C
C PIVOT(1:N+2)  IS A WORK ARRAY WHOSE FIRST N+1 COMPONENTS CONTAIN
C    THE INDICES OF THE DIAGONAL ELEMENTS OF THE N X N SYMMETRIC 
C    JACOBIAN MATRIX (WITH RESPECT TO X) WITHIN  QR.
C
C PP(1:N)  IS A WORK ARRAY CONTAINING THE NEGATIVE OF THE LAST COLUMN
C    OF THE JACOBIAN MATRIX  -[D RHO/D LAMBDA].
C
C RHOVEC(1:N+1), Z(1:N+1), DZ(1:N+1)  ARE ALL WORK ARRAYS
C    USED TO CALCULATE THE NEWTON STEPS.
C
C WORK(1:6*(N+1)+LENQR)  IS A WORK ARRAY USED BY THE CONJUGATE GRADIENT
C    ALGORITHM TO SOLVE LINEAR SYSTEMS.
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, RHOJS.
C
C
C ON OUTPUT:
C
C N, LENQR, RELERR, ABSERR, 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 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 OBTAINED BY NEWTON STEPS, AND  YP
C        CONTAINS A POINT OPPOSITE OF THE HYPERPLANE  LAMBDA=1  FROM  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, DNRM2, F (OR RHO), 
C        FJACS (OR RHOJS), PCGQS, ROOT
C
C ***** DECLARATIONS ***** 
C
C     FUNCTION DECLARATIONS 
C
        DOUBLE PRECISION D1MACH, DNRM2, QOFS
C
C     LOCAL VARIABLES 
C
        DOUBLE PRECISION AERR, DD001, DD0011, DD01, DD011, DELS,
     &     LAMBDA, ONE, P0, P1, PP0, PP1, QSOUT, RERR, S, SA, SB,
     &     SIGMA, SOUT, U, ZERO
        INTEGER ISTEP, I, J, LCODE, LIMIT, NP1, ZU
        LOGICAL BRACK
C
C     SCALAR ARGUMENTS 
C
        DOUBLE PRECISION RELERR, ABSERR
        INTEGER N, NFE, IFLAG, LENQR
C
C     ARRAY DECLARATIONS 
C
        DOUBLE PRECISION Y(N+1), YP(N+1), YOLD(N+1), YPOLD(N+1), A(N),
     &     QR(LENQR), PP(N), RHOVEC(N+1), Z(N+1), DZ(N+1), 
     &     WORK(6*(N+1)+LENQR), PAR(1)      
        INTEGER PIVOT(N+2), 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 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
        LIMIT = 2*(INT(-LOG10(AERR+RERR*DNRM2(NP1,Y,1)))+1)
        ZU=N+2
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(NP1),YPOLD(NP1),Y(NP1),YP(NP1),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 ***** 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 SET STARTING POINTS FOR CONJUGATE GRADIENT ALGORITHM TO ZERO.
C 
          DO 70 J=ZU,ZU+2*N+1
            WORK(J) = 0.0
  70      CONTINUE
C
C COMPUTE NEWTON STEP.
C
C       COMPUTE QR = [D RHO/DX], RHOVEC = RHO, -PP = (D RHO/D LAMBDA).
C
          LAMBDA = Z(NP1)
          IF (IFLAG .EQ. -2) THEN
C 
C           CURVE TRACKING PROBLEM.
C
            CALL RHOJS(A,LAMBDA,Z,QR,LENQR,PIVOT,PP,PAR,IPAR)
            CALL RHO(A,LAMBDA,Z,RHOVEC,PAR,IPAR)
          ELSE IF (IFLAG .EQ. -1) THEN
C
C           ZERO FINDING PROBLEM.
C
            CALL FJACS(Z,QR,LENQR,PIVOT)
            CALL DSCAL(LENQR,LAMBDA,QR,1)
            SIGMA = 1.0-LAMBDA
            DO 80 J=1,N
              QR(PIVOT(J))=QR(PIVOT(J))+SIGMA              
  80        CONTINUE
            CALL DCOPY(N,Z,1,RHOVEC,1)
            CALL DAXPY(N,-ONE,A,1,RHOVEC,1)
            CALL F(Z,PP)
            CALL DSCAL(N,-ONE,PP,1)
            CALL DAXPY(N,ONE,RHOVEC,1,PP,1)
            CALL DAXPY(N,-LAMBDA,PP,1,RHOVEC,1)
          ELSE
C
C           FIXED POINT PROBLEM.
C
            CALL FJACS(Z,QR,LENQR,PIVOT)
            CALL DSCAL(LENQR,-LAMBDA,QR,1)
            DO 90 J=1,N
              QR(PIVOT(J))=QR(PIVOT(J))+1.0
  90        CONTINUE
            CALL DCOPY(N,Z,1,RHOVEC,1)
            CALL DAXPY(N,-ONE,A,1,RHOVEC,1)
            CALL F(Z,PP)
            CALL DAXPY(N,-ONE,A,1,PP,1)
            CALL DAXPY(N,-LAMBDA,PP,1,RHOVEC,1)
          END IF
          RHOVEC(NP1) = 0.0       
          NFE = NFE+1
C
C       SOLVE SYSTEM TO FIND NEWTON STEP.
C
          CALL PCGQS(N,QR,LENQR,PIVOT,PP,YP,RHOVEC,DZ,WORK,IFLAG)
          IF (IFLAG .GT. 0) RETURN
C
C       TAKE NEWTON STEP.
C
          CALL DAXPY(NP1,ONE,DZ,1,Z,1)
C
C       CHECK FOR CONVERGENCE. 
C
          IF ((ABS(Z(NP1)-1.0) .LE. RERR+AERR) .AND. 
     &        (DNRM2(NP1,DZ,1) .LE. RERR*DNRM2(N,Z,1)+AERR)) THEN
             RETURN
          END IF
C
C PREPARE FOR NEXT ITERATION.
C
C       IF LAMBDA COMPONENT OF  Z=1, THEN DO NOT COMPUTE A 
C       NEW PREDICTOR, BUT RATHER CONTINUE WITH ANOTHER NEWTON
C       ITERATION.
C
          IF (ABS(Z(NP1)-1.0) .LT. RERR+AERR) GOTO 300
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 OPPOSITE 
C       OF  LAMBDA=1  FROM  Y.  SET  BRACK = .TRUE.  IFF  Y & YOLD  
C       BRACKET  LAMBDA=1  SO THAT  YPOLD=YOLD. 
C
          IF ((YOLD(NP1)-1.0)*(Y(NP1)-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 = DZ + Y,
C           WHERE  DZ = SA*(YOLD-Y).
C
          SA = (1.0-Y(NP1))/(YOLD(NP1)-Y(NP1))
          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(NP1))/(YPOLD(NP1)-Y(NP1))
              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 = DZ+Y.
C
          CALL DCOPY(NP1,Y,1,Z,1)
          CALL DAXPY(NP1,ONE,DZ,1,Z,1)
  300   CONTINUE
C
C ***** END OF MAIN LOOP. *****
C
C THE ALTERNATING OSCULATORY LINEAR PREDICTION AND NEWTON 
C CORRECTION HAS NOT CONVERGED IN  LIMIT  STEPS.  ERROR RETURN.
        IFLAG=6
        RETURN
C
C ***** END OF SUBROUTINE ROOTQS *****
        END