SUBROUTINE FIXPQS(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,
     &     A,NFE,ARCLEN,YP,YOLD,YPOLD,QR,LENQR,PIVOT,PP,RHOVEC,Z0,DZ,
     &     T,WORK,SSPAR,PAR,IPAR)
C
C SUBROUTINE  FIXPQS  FINDS A FIXED POINT OR ZERO OF THE 
C N-DIMENSIONAL VECTOR FUNCTION  F(X), OR TRACKS A ZERO CURVE OF A 
C GENERAL HOMOTOPY MAP  RHO(A,X,LAMBDA).  FOR THE FIXED POINT PROBLEM
C F(X) IS ASSUMED TO BE A C2 MAP OF SOME BALL INTO ITSELF.  THE 
C EQUATION  X=F(X)  IS SOLVED BY FOLLOWING THE ZERO CURVE OF THE 
C HOMOTOPY MAP
C
C  LAMBDA*(X - F(X)) + (1 - LAMBDA)*(X - A),
C
C STARTING FROM  LAMBDA = 0, X = A.   THE CURVE IS PARAMETERIZED
C BY ARC LENGTH  S, AND IS FOLLOWED BY SOLVING THE ORDINARY 
C DIFFERENTIAL EQUATION  D(HOMOTOPY MAP)/DS = 0  FOR  
C Y(S) = (X(S),LAMBDA(S)).  THIS IS DONE BY USING A HERMITE CUBIC 
C PREDICTOR AND A CORRECTOR WHICH RETURNS TO THE ZERO CURVE IN A 
C HYPERPLANE PERPENDICULAR TO THE TANGENT TO THE ZERO CURVE AT THE 
C MOST RECENT POINT.
C
C FOR THE ZERO FINDING PROBLEM  F(X)  IS ASSUMED TO BE A C2 MAP SUCH
C THAT FOR SOME  R > 0,  X*F(X) >= 0  WHENEVER  NORM(X) = R.
C THE EQUATION  F(X) = 0  IS SOLVED BY FOLLOWING THE ZERO CURVE OF
C THE HOMOTOPY MAP
C
C  LAMBDA*F(X) + (1 - LAMBDA)*(X - A)
C
C EMANATING FROM  LAMBDA = 0, X = A.
C
C A  MUST BE AN INTERIOR POINT OF THE ABOVE MENTIONED BALLS.
C
C FOR THE CURVE TRACKING PROBLEM RHO(A,X,LAMBDA) IS ASSUMED TO 
C BE A C2 MAP FROM  E**M X [0,1) X E**N  INTO  E**N, WHICH FOR 
C ALMOST ALL PARAMETER VECTORS  A  IN SOME NONEMPTY OPEN SUBSET
C OF E**M SATISFIES
C
C  RANK [D RHO(A,X,LAMBDA)/D LAMBDA, D RHO(A,X,LAMBDA)/DX] = N
C
C FOR ALL POINTS  (X,LAMBDA)  SUCH THAT  RHO(A,X,LAMBDA) = 0.  IT IS
C FURTHER ASSUMED THAT
C
C         RANK [ D RHO(A,X0,0)/DX ] = N.
C
C WITH  A  FIXED, THE ZERO CURVE OF  RHO(A,X,LAMBDA)  EMANATING FROM
C LAMBDA = 0, X = X0  IS TRACKED UNTIL  LAMBDA = 1  BY SOLVING THE 
C ORDINARY DIFFERENTIAL EQUATION  D RHO(A,X(S),LAMBDA(S))/DS = 0
C FOR  Y(S) = (X(S),LAMBDA(S)), WHERE  S  IS ARC LENGTH ALONG THE
C ZERO CURVE.  ALSO THE HOMOTOPY MAP  RHO(A,X,LAMBDA)  IS ASSUMED TO
C BE CONSTRUCTED SUCH THAT
C
C         D LAMBDA(0)/DS > 0.
C
C FOR THE FIXED POINT AND ZERO FINDING PROBLEMS, THE USER MUST SUPPLY
C A SUBROUTINE  F(X,V)  WHICH EVALUATES F(X) AT X AND RETURNS THE
C VECTOR F(X) IN  V,  AND A SUBROUTINE  FJACS(X,QR,LENQR,PIVOT)  WHICH 
C EVALUATES THE (SYMMETRIC) JACOBIAN MATRIX OF F(X) AT X, AND RETURNS 
C THE SYMMETRIC JACOBIAN MATRIX IN PACKED SKYLINE STORAGE FORMAT IN  QR.
C LENQR  AND  PIVOT  DESCRIBE THE DATA STRUCTURE IN  QR.  FOR THE CURVE
C TRACKING PROBLEM, THE USER MUST SUPPLY A SUBROUTINE
C  RHO(A,LAMBDA,X,V,PAR,IPAR)  WHICH EVALUATES THE HOMOTOPY MAP RHO
C AT (A,X,LAMBDA) AND RETURNS THE VECTOR RHO(A,X,LAMBDA) IN  V, 
C AND A SUBROUTINE  RHOJS(A,LAMBDA,X,QR,LENQR,PIVOT,PP,PAR,IPAR)  WHICH
C RETURNS IN  QR  THE SYMMETRIC N X N JACOBIAN MATRIX [D RHO/DX]
C EVALUATED AT (A,X,LAMBDA) AND STORED IN PACKED SKYLINE FORMAT,
C AND RETURNS IN  PP  THE VECTOR  -(D RHO/D LAMBDA)  EVALUATED AT
C (A,X,LAMBDA).  LENQR  AND  PIVOT  DESCRIBE THE DATA STRUCTURE IN
C QR.
C *** NOTE THE MINUS SIGN IN THE DEFINITION OF  PP. ***
C
C
C FIXPQS DIRECTLY OR INDIRECTLY USES THE SUBROUTINES  D1MACH, F 
C (OR RHO), FJACS (OR RHOJS), GMFADS, MULTDS, PCGQS, ROOTQS, STEPQS,
C SOLVDS,  AND THE BLAS ROUTINES  DAXPY, DCOPY, DDOT, DNRM2, AND  DSCAL.  
C ONLY  D1MACH  CONTAINS MACHINE DEPENDENT CONSTANTS.  NO OTHER 
C MODIFICATIONS BY THE USER ARE REQUIRED.
C
C
C ON INPUT:
C
C N  IS THE DIMENSION OF X, F(X), AND RHO(A,X,LAMBDA).
C
C Y(1:N+1)  CONTAINS THE STARTING POINT FOR TRACKING THE HOMOTOPY MAP.
C    (Y(1),...,Y(N)) = A  FOR THE FIXED POINT AND ZERO FINDING 
C    PROBLEMS.  (Y(1),...,Y(N)) = X0  FOR THE CURVE TRACKING PROBLEM.
C    Y(N+1)  NEED NOT BE DEFINED BY THE USER.
C
C IFLAG  CAN BE -2, -1, 0, 2, OR 3.  IFLAG  SHOULD BE 0 ON THE FIRST
C    CALL TO  FIXPQS  FOR THE PROBLEM  X=F(X), -1 FOR THE PROBLEM
C    F(X)=0, AND -2 FOR THE PROBLEM  RHO(A,X,LAMBDA)=0.   IN CERTAIN
C    SITUATIONS  IFLAG  IS SET TO 2 OR 3 BY  FIXPQS, AND  FIXPQS  CAN
C    BE CALLED AGAIN WITHOUT CHANGING  IFLAG.
C
C ARCRE, ARCAE  ARE THE RELATIVE AND ABSOLUTE ERRORS, RESPECTIVELY,
C    ALLOWED THE ITERATION ALONG THE ZERO CURVE.  IF
C    ARC?E .LE. 0.0  ON INPUT, IT IS RESET TO  .5*SQRT(ANS?E).
C    NORMALLY  ARC?E  SHOULD BE CONSIDERABLY LARGER THAN  ANS?E.
C
C ANSRE, ANSAE  ARE THE RELATIVE AND ABSOLUTE ERROR VALUES USED FOR 
C    THE ANSWER AT  LAMBDA = 1.  THE ACCEPTED ANSWER  Y = (X,LAMBDA)
C    SATISFIES
C
C      |Y(1) - 1| .LE. ANSRE + ANSAE      .AND.
C  
C      ||DZ|| .LE. ANSRE*||Y|| + ANSAE      WHERE
C
C      DZ IS THE NEWTON STEP TO Y.
C
C TRACE  IS AN INTEGER SPECIFYING THE LOGICAL I/O UNIT FOR
C    INTERMEDIATE OUTPUT.  IF  TRACE .GT. 0  THE POINTS COMPUTED ON
C    THE ZERO CURVE ARE WRITTEN TO I/O UNIT  TRACE .
C
C A(1:*)  CONTAINS THE PARAMETER VECTOR  A.  FOR THE FIXED POINT
C    AND ZERO FINDING PROBLEMS,  A  NEED NOT BE INITIALIZED BY THE 
C    USER, AND IS ASSUMED TO HAVE LENGTH  N.  FOR THE CURVE
C    TRACKING PROBLEM,  A  MUST BE INITIALIZED BY THE USER.
C
C YP(1:N+1)  IS A WORK ARRAY CONTAINING THE TANGENT VECTOR TO THE
C    ZERO CURVE AT THE CURRENT POINT  Y.
C
C YOLD(1:N+1) IS A WORK ARRAY CONTAINING THE PREVIOUS POINT FOUND
C    ON THE ZERO CURVE.
C
C YPOLD(1:N+1) IS A WORK ARRAY CONTAINING THE TANGENT VECTOR TO
C    THE ZERO CURVE AT  YOLD.
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 LENQR  IS THE LENGTH OF THE N-DIMENSIONAL ARRAY  QR.  I.E.
C    IT IS THE NUMBER OF NON-ZERO ENTRIES IN THE JACOBIAN 
C    MATRIX [DF/DX] (OR [D RHO/DX]).
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), Z0(1:N+1), DZ(1:N+1), T(1:N+1)  ARE ALL WORK ARRAYS
C    USED BY  STEPQS, TANGQS, AND ROOTQS  TO CALCULATE THE TANGENT
C    VECTORS AND NEWTON STEPS.
C
C WORK(1:8*(N+1)+LENQR)  IS A WORK ARRAY USED BY THE CONJUGATE GRADIENT
C    ALGORITHM TO SOLVE LINEAR SYSTEMS.
C
C SSPAR(1:4) =  (HMIN, HMAX, BMIN, BMAX)  IS A VECTOR OF PARAMETERS 
C    USED FOR THE OPTIMAL STEP SIZE ESTIMATION.  A DEFAULT VALUE
C    CAN BE SPECIFIED FOR ANY OF THESE FOUR PARAMETERS BY SETTING IT
C    .LE. 0.0  ON INPUT.  SEE THE COMMENTS IN  STEPQS  FOR MORE
C    INFORMATION ABOUT THESE PARAMETERS.
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 , TRACE , A , LENQR  ARE UNCHANGED.
C
C Y(N+1) = LAMBDA, (Y(1),...,Y(N)) = X, AND  Y  IS AN APPROXIMATE
C    ZERO OF THE HOMOTOPY MAP.  NORMALLY  LAMBDA = 1  AND  X  IS A
C    FIXED POINT OR ZERO OF  F(X).   IN ABNORMAL SITUATIONS,  LAMBDA
C    MAY ONLY BE NEAR 1 AND  X  NEAR A FIXED POINT OR ZERO.
C
C IFLAG =
C
C   1   NORMAL RETURN.
C
C   2   SPECIFIED ERROR TOLERANCE CANNOT BE MET.  SOME OR ALL OF
C       ARCRE, ARCAE, ANSRE, ANSAE  HAVE BEEN INCREASED TO 
C       SUITABLE VALUES.  TO CONTINUE, JUST CALL  FIXPQS  AGAIN
C       WITHOUT CHANGING ANY PARAMETERS.
C
C   3   STEPQS  HAS BEEN CALLED 1000 TIMES.  TO CONTINUE, CALL
C       FIXPQS  AGAIN WITHOUT CHANGING ANY PARAMETERS.
C
C   4   JACOBIAN MATRIX DOES NOT HAVE FULL RANK.  THE ALGORITHM
C       HAS FAILED (THE ZERO CURVE OF THE HOMOTOPY MAP CANNOT BE
C       FOLLOWED ANY FURTHER).
C
C   5   THE TRACKING ALGORITHM HAS LOST THE ZERO CURVE OF THE 
C       HOMOTOPY MAP AND IS NOT MAKING PROGRESS.  THE ERROR 
C       TOLERANCES  ARC?E  AND  ANS?E  WERE TOO LENIENT.  THE PROBLEM 
C       SHOULD BE RESTRARTED BY CALLING  FIXPQS  WITH SMALLER ERROR 
C       TOLERANCES AND  IFLAG = 0 (-1, -2).
C
C   6   THE NEWTON ITERATION IN  STEPQS  OR  ROOTQS  FAILED TO
C       CONVERGE.  THE ERROR TOLERANCES  ANS?E  MAY BE TOO STRINGENT.
C
C   7   ILLEGAL INPUT PARAMETERS, A FATAL ERROR.
C
C ARCRE, ARCAE, ANSRE, ANSAE  ARE UNCHANGED AFTER A NORMAL RETURN
C    (IFLAG = 1).  THEY ARE INCREASED TO APPROPRIATE VALUES ON THE
C    RETURN  IFLAG = 2.
C
C NFE  IS THE NUMBER OF JACOBIAN EVALUATIONS.
C
C ARCLEN  IS THE APPROXIMATE LENGTH OF THE ZERO CURVE.  
C
C ***** DECLARATIONS *****
C
C     FUNCTION DECLARATIONS 
C 
        DOUBLE PRECISION D1MACH, DNRM2
C
C     LOCAL VARIABLES 
C
        DOUBLE PRECISION ABSERR, H, HOLD, RELERR, S, WK 
        INTEGER IFLAGC, ITER, JW, LIMITD, LIMIT, NP1, PCGWK
        LOGICAL CRASH, START       
C
C     SCALAR ARGUMENTS 
C
        DOUBLE PRECISION ARCRE, ARCAE, ANSRE, ANSAE, ARCLEN
        INTEGER N, IFLAG, TRACE, NFE, LENQR
C
C     ARRAY DECLARATIONS 
C
        DOUBLE PRECISION A(N), Y(N+1), YP(N+1), YOLD(N+1), YPOLD(N+1),
     &    QR(LENQR), PP(N), RHOVEC(N+1), Z0(N+1), DZ(N+1), T(N+1),
     &    WORK(8*(N+1)+LENQR), SSPAR(4), PAR(1)
        INTEGER PIVOT(N+2), IPAR(1)
C 
        SAVE
C
C ***** END OF DECLARATIONS *****
C
C LIMITD  IS AN UPPER BOUND ON THE NUMBER OF STEPS.  IT MAY BE 
C CHANGED BY CHANGING THE FOLLOWING PARAMETER STATEMENT:
        PARAMETER (LIMITD =1000)
C
C ***** FIRST EXECUTABLE STATEMENT *****
C
C CHECK IFLAG
C
        IF (N .LE. 0 .OR. ANSRE .LE. 0.0 .OR. ANSAE .LT. 0.0)
     &    IFLAG = 7
        IF (IFLAG .GE. -2 .AND. IFLAG .LE. 0) GO TO 10
        IF (IFLAG .EQ. 2) GO TO 50
        IF (IFLAG .EQ. 3) GO TO 40
C
C ONLY VALID INPUT FOR IFLAG IS -2, -1, 0, 2, 3.
C
        IFLAG = 7
        RETURN
C
C ***** INITIALIZATION BLOCK  *****
C
  10    ARCLEN = 0.0
        IF (ARCRE .LE. 0.0) ARCRE = .5*SQRT(ANSRE)
        IF (ARCAE .LE. 0.0) ARCAE = .5*SQRT(ANSAE)
        NFE=0
        IFLAGC = IFLAG
        NP1=N+1
        PCGWK = 2*N+3
C 
C SET INITIAL CONDITIONS FOR FIRST CALL TO STEPQS.
C
        START=.TRUE.
        CRASH=.FALSE.
        RELERR = ARCRE
        ABSERR = ARCAE
        HOLD=1.0
        H=0.1
        S=0.0
        YPOLD(NP1) = 1.0
        Y(NP1) = 0.0
        DO 20 JW=1,N
          YPOLD(JW)=0.0
  20    CONTINUE
C
C SET OPTIMAL STEP SIZE ESTIMATION PARAMETERS.
C
C     MINIMUM STEP SIZE HMIN
        IF (SSPAR(1) .LE. 0.0) SSPAR(1)= (SQRT(N+1.0)+4.0)*D1MACH(4)
C     MAXIMUM STEP SIZE HMAX
        IF (SSPAR(2) .LE. 0.0) SSPAR(2)= 1.0
C     MINIMUM STEP REDUCTION FACTOR BMIN
        IF (SSPAR(3) .LE. 0.0) SSPAR(3)= 0.1
C     MAXIMUM STEP EXPANSION FACTOR BMAX
        IF (SSPAR(4) .LE. 0.0) SSPAR(4)= 7.0
C
C LOAD  A  FOR THE FIXED POINT AND ZERO FINDING PROBLEMS.
C
        IF (IFLAGC .GE. -1) THEN
          CALL DCOPY(N,Y,1,A,1)
        ENDIF
C
  40    LIMIT=LIMITD
C
C ***** END OF INITIALIZATION BLOCK. *****
C
C ***** MAIN LOOP. *****
C
  50    DO 400 ITER=1,LIMIT
        IF (Y(NP1) .LT. 0.0) THEN
          ARCLEN = S
          IFLAG = 5
          RETURN
        END IF
C
C TAKE A STEP ALONG THE CURVE.
C
        CALL STEPQS(N,NFE,IFLAGC,LENQR,START,CRASH,HOLD,H,WK,RELERR,
     &    ABSERR,S,Y,YP,YOLD,YPOLD,A,QR,PIVOT,PP,RHOVEC,Z0,DZ,T,
     &    WORK,SSPAR,PAR,IPAR)
C
C PRINT LATEST POINT ON CURVE IF REQUESTED.
C
      IF (TRACE .GT. 0) THEN
	WRITE (TRACE,217) ITER,NFE,S,Y(NP1),(Y(JW),JW=1,N)
217     FORMAT(/' STEP',I5,3X,'NFE =',I5,3X,'ARC LENGTH =',F9.4,3X,
     &  'LAMBDA =',F7.4,5X,'X vector:'/1P,(1X,6E12.4))
      ENDIF
C
C CHECK IF THE STEP WAS SUCCESSFUL.
C
        IF (IFLAGC .GT. 0) THEN
          ARCLEN=S
          IFLAG=IFLAGC
          RETURN
        END IF
C
        IF (CRASH) THEN
C
C         RETURN CODE FOR ERROR TOLERANCE TOO SMALL.
C      
          IFLAG=2
C
C         CHANGE ERROR TOLERANCES.
C
          IF (ARCRE .LT. RELERR) THEN
            ARCRE=RELERR
            ANSRE=RELERR
          ENDIF
          IF (ARCAE .LT. ABSERR) ARCAE=ABSERR
C
C         CHANGE LIMIT ON NUMBER OF ITERATIONS.
C
          LIMIT = LIMIT - ITER
          RETURN
        END IF
C
C IF  LAMBDA >= 1.0,  USE  ROOTQS  TO FIND SOLUTION.
C
        IF (Y(NP1) .GE. 1.0) GOTO 500
C
  400   CONTINUE
C
C ***** END OF MAIN LOOP *****
C
C DID NOT CONVERGE IN  LIMIT  ITERATIONS, SET  IFLAG  AND RETURN.
C
        ARCLEN = S
        IFLAG = 3
        RETURN
C
C ***** FINAL STEP -- FIND SOLUTION AT LAMBDA=1 *****
C
C SAVE  YOLD  FOR ARC LENGTH CALCULATION LATER.
C
  500   CALL DCOPY(NP1,YOLD,1,T,1)
C
C FIND SOLUTION.
C
        CALL ROOTQS(N,NFE,IFLAGC,LENQR,ANSRE,ANSAE,Y,YP,YOLD,YPOLD,
     &    A,QR,PIVOT,PP,RHOVEC,Z0,DZ,WORK(PCGWK),PAR,IPAR)
C
C CHECK IF SOLUTION WAS FOUND AND SET  IFLAG  ACCORDINGLY.
C
        IFLAG=1
C
C     SET ERROR FLAG IF ROOTQS COULD NOT GET THE POINT ON THE ZERO
C     CURVE AT  LAMBDA = 1.0 .
C
        IF (IFLAGC .GT. 0) IFLAG=IFLAGC
C
C CALCULATE FINAL ARC LENGTH.
C
        CALL DCOPY(NP1,Y,1,DZ,1)
        WK=-1.0
        CALL DAXPY(NP1,WK,T,1,DZ,1)
        ARCLEN=S - HOLD + DNRM2(NP1,DZ,1)
C
C ***** END OF FINAL STEP *****
C
        RETURN
C
C ***** END OF SUBROUTINE FIXPQS *****
        END