SUBROUTINE FIXPQF(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,A,
     $     NFE,ARCLEN,YP,YOLD,YPOLD,QT,R,F0,F1,Z0,DZ,W,T,YSAV,
     $     SSPAR,PAR,IPAR)
C
C SUBROUTINE  FIXPQF  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,LAMBDA,X).  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) = (LAMBDA(S), X(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
C SUCH 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,LAMBDA,X) 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,LAMBDA,X)/D LAMBDA, D RHO(A,LAMBDA,X)/DX] = N
C
C FOR ALL POINTS  (LAMBDA,X)  SUCH THAT  RHO(A,LAMBDA,X) = 0.  IT IS
C FURTHER ASSUMED THAT
C
C         RANK [ D RHO(A,0,X0)/DX ] = N.
C
C WITH  A  FIXED, THE ZERO CURVE OF  RHO(A,LAMBDA,X)  EMANATING FROM
C LAMBDA = 0, X = X0  IS TRACKED UNTIL  LAMBDA = 1  BY SOLVING THE
C ORDINARY DIFFERENTIAL EQUATION    D RHO(A,LAMBDA(S),X(S))/DS = 0
C FOR  Y(S) = (LAMBDA(S), X(S)), WHERE  S  IS ARC LENGTH ALONG THE
C ZERO CURVE.  ALSO THE HOMOTOPY MAP  RHO(A,LAMBDA,X)  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  FJAC(X,V,K)  WHICH RETURNS IN  V
C THE KTH COLUMN OF THE JACOBIAN MATRIX OF F(X) EVALUATED AT X.  FOR
C THE CURVE TRACKING PROBLEM, THE USER MUST SUPPLY A SUBROUTINE
C RHO(A,LAMBDA,X,V,PAR,IPAR)  WHICH EVALUATES THE HOMOTOPY MAP  RHO AT
C (A,LAMBDA,X)  AND RETURNS THE VECTOR  RHO(A,LAMBDA,X)  IN  V, AND
C A SUBROUTINE  RHOJAC(A,LAMBDA,X,V,K,PAR,IPAR)  WHICH RETURNS IN  V
C THE KTH COLUMN OF THE  N X (N+1)  JACOBIAN MATRIX
C [D RHO/D LAMBDA, D RHO/DX]  EVALUATED AT  (A,LAMBDA,X).  FIXPQF
C DIRECTLY OR INDIRECTLY USES THE SUBROUTINES  D1MACH, F (OR RHO),
C FJAC (OR RHOJAC), QRFAQF, QRSLQF, ROOT, ROOTQF, STEPQF, TANGQF,
C UPQRQF  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,LAMBDA,X).
C
C Y(1:N+1)  CONTAINS THE STARTING POINT FOR TRACKING THE HOMOTOPY MAP.
C    (Y(2),...,Y(N+1)) = A  FOR THE FIXED POINT AND ZERO FINDING
C    PROBLEMS.  (Y(2),...,Y(N+1)) = X0  FOR THE CURVE TRACKING PROBLEM.
C    Y(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  FIXPQF  FOR THE PROBLEM  X=F(X), -1 FOR THE PROBLEM
C    F(X)=0, AND -2 FOR THE PROBLEM  RHO(A,LAMBDA,X)=0.   IN CERTAIN
C    SITUATIONS  IFLAG  IS SET TO 2 OR 3 BY  FIXPQF, AND  FIXPQF  CAN
C    BE CALLED AGAIN WITHOUT CHANGING  IFLAG.
C
C ARCRE, ARCAE  ARE THE RELATIVE AND ABSOLUTE ERRORS, RESPECTIVELY,
C    ALLOWED THE QUASI-NEWTON 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 = (LAMBDA, X)
C    SATISFIES
C
C      |Y(1) - 1| .LE. ANSRE + ANSAE      .AND.
C
C      ||DZ|| .LE. ANSRE*||Y|| + ANSAE      WHERE
C
C      DZ IS THE QUASI-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 QT(1:N+1,1:N+1), R((N+1)*(N+2)/2), F0(1:N+1), F1(1:N+1), Z0(1:N+1),
C    DZ(1:N+1), W(1:N+1), T(1:N+1), YSAV(1:N+1)  ARE ALL WORK ARRAYS
C    USED BY  STEPQF, TANGQF AND ROOTQF TO CALCULATE THE TANGENT
C    VECTORS AND QUASI-NEWTON STEPS.
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  STEPQF  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, RHOJAC.
C
C
C ON OUTPUT:
C
C N , TRACE , A  ARE UNCHANGED.
C
C Y(1) = LAMBDA, (Y(2),...,Y(N+1)) = 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  FIXPQF  AGAIN
C       WITHOUT CHANGING ANY PARAMETERS.
C
C   3   STEPQF  HAS BEEN CALLED 1000 TIMES.  TO CONTINUE, CALL
C       FIXPQF  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  FIXPQF  WITH SMALLER ERROR
C       TOLERANCES AND  IFLAG = 0 (-1, -2).
C
C   6   THE QUASI-NEWTON ITERATION IN  STEPQF  OR  ROOTQF  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
        LOGICAL CRASH, START
C
C     SCALAR ARGUMENTS
C
        DOUBLE PRECISION ARCRE, ARCAE, ANSRE, ANSAE, ARCLEN
        INTEGER N,IFLAG,TRACE,NFE
C
C     ARRAY DECLARATIONS
C
        DOUBLE PRECISION A(N), Y(N+1), YP(N+1), YOLD(N+1), YPOLD(N+1),
     $    QT(N+1,N+1), R((N+1)*(N+2)/2), F0(N+1), F1(N+1), Z0(N+1),
     $    DZ(N+1), W(N+1), T(N+1), YSAV(N+1), SSPAR(4), PAR(1)
        INTEGER 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
C
C SET INITIAL CONDITIONS FOR FIRST CALL TO STEPQF.
C
        START=.TRUE.
        CRASH=.FALSE.
        RELERR = ARCRE
        ABSERR = ARCAE
        HOLD=1.0
        H=0.1
        S=0.0
        YPOLD(1) = 1.0
        Y(1) = 0.0
        DO 20 JW=2,NP1
          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(2),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(1) .LT. 0.0) THEN
          ARCLEN = S
          IFLAG = 5
          RETURN
        END IF
C
C TAKE A STEP ALONG THE CURVE.
C
        CALL STEPQF(N,NFE,IFLAGC,START,CRASH,HOLD,H,WK,
     $    RELERR,ABSERR,S,Y,YP,YOLD,YPOLD,A,QT,R,F0,F1,Z0,DZ,
     $    W,T,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(1),(Y(JW),JW=2,NP1)
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
          END IF
          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 ROOTQF TO FIND SOLUTION.
C
        IF (Y(1) .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,YSAV,1)
C
C FIND SOLUTION.
C
        CALL ROOTQF(N,NFE,IFLAGC,ANSRE,ANSAE,Y,YP,YOLD,
     $    YPOLD,A,QT,R,DZ,Z0,W,T,F0,F1,PAR,IPAR)
C
C CHECK IF SOLUTION WAS FOUND AND SET  IFLAG  ACCORDINGLY.
C
        IFLAG=1
C
C     SET ERROR FLAG IF ROOTQF 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,YSAV,1,DZ,1)
        ARCLEN = S - HOLD + DNRM2(NP1,DZ,1)
C
C ***** END OF FINAL STEP *****
C
        RETURN
C
C ***** END OF SUBROUTINE FIXPQF *****
        END