      SUBROUTINE FIXPDF(N,Y,IFLAG,ARCTOL,EPS,TRACE,A,NDIMA,NFE,
     $     ARCLEN,YP,YPOLD,QR,ALPHA,TZ,PIVOT,WT,PHI,P,PAR,IPAR)
C
C SUBROUTINE  FIXPDF  FINDS A FIXED POINT OR ZERO OF THE
C N-DIMENSIONAL VECTOR FUNCTION F(X), OR TRACKS A ZERO CURVE
C OF A GENERAL HOMOTOPY MAP RHO(A,LAMBDA,X).  FOR THE FIXED
C POINT PROBLEM F(X) IS ASSUMED TO BE A C2 MAP OF SOME BALL
C INTO ITSELF.  THE EQUATION  X = F(X)  IS SOLVED BY
C FOLLOWING THE ZERO CURVE OF THE 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)).
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
C OF 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
C FROM  LAMBDA = 0, X = X0  IS TRACKED UNTIL  LAMBDA = 1  BY
C SOLVING THE ORDINARY DIFFERENTIAL EQUATION
C D RHO(A,LAMBDA(S),X(S))/DS = 0  FOR  Y(S) = (LAMBDA(S), X(S)),
C WHERE S IS ARC LENGTH ALONG THE ZERO CURVE.  ALSO THE HOMOTOPY
C MAP RHO(A,LAMBDA,X) IS ASSUMED TO BE CONSTRUCTED SUCH THAT
C
C              D LAMBDA(0)/DS > 0  .
C
C THIS CODE IS BASED ON THE ALGORITHM IN L. T. WATSON, A
C GLOBALLY CONVERGENT ALGORITHM FOR COMPUTING FIXED POINTS OF
C C2 MAPS, APPL. MATH. COMPUT., 5 (1979) 297-311.
C
C
C FOR THE FIXED POINT AND ZERO FINDING PROBLEMS, THE USER
C MUST SUPPLY A SUBROUTINE  F(X,V)  WHICH EVALUATES F(X) AT X
C AND RETURNS THE VECTOR F(X) IN V, AND A SUBROUTINE  FJAC(X,V,K)
C WHICH RETURNS IN V THE KTH COLUMN OF THE JACOBIAN MATRIX OF
C F(X) EVALUATED AT X.  FOR THE CURVE TRACKING PROBLEM, THE USER MUST
C SUPPLY A SUBROUTINE  RHOA(V,LAMBDA,X,PAR,IPAR)  WHICH GIVEN
C (LAMBDA,X) RETURNS A PARAMETER VECTOR A IN V SUCH THAT
C RHO(A,LAMBDA,X)=0, AND A SUBROUTINE  RHOJAC(A,LAMBDA,X,V,K,PAR,IPAR)
C WHICH RETURNS IN V THE KTH COLUMN OF THE N X (N+1) JACOBIAN
C MATRIX [D RHO/D LAMBDA, D RHO/DX] EVALUATED AT (A,LAMBDA,X).
C  FIXPDF  DIRECTLY OR INDIRECTLY USES THE SUBROUTINES
C  STEPS , SINTRP , ROOT , FODE , F (OR  RHOA ),
C FJAC (OR  RHOJAC ), DCPOSE , D1MACH , AND THE BLAS FUNCTIONS
C DDOT  AND  DNRM2 .  ONLY  D1MACH  CONTAINS MACHINE
C DEPENDENT CONSTANTS.  NO OTHER MODIFICATIONS BY THE USER ARE
C REQUIRED.
C
C ***WARNING:  THIS SUBROUTINE IS GENERALLY MORE ROBUST THAN  FIXPNF
C AND  FIXPQF , BUT MAY BE SLOWER THAN THOSE SUBROUTINES BY A
C FACTOR OF TWO.
C
C
C ON INPUT:
C
C N  IS THE DIMENSION OF X, F(X), AND RHO(A,LAMBDA,X).
C
C Y  IS AN ARRRAY OF LENGTH  N + 1.  (Y(2),...,Y(N+1)) = A  IS THE
C    STARTING POINT FOR THE ZERO CURVE FOR THE FIXED POINT AND
C    ZERO FINDING PROBLEMS.  (Y(2),...,Y(N+1)) = X0  FOR THE CURVE
C    TRACKING PROBLEM.
C
C IFLAG  CAN BE -2, -1, 0, 2, OR 3.  IFLAG  SHOULD BE 0 ON THE
C    FIRST CALL TO  FIXPDF  FOR THE PROBLEM  X=F(X), -1 FOR THE
C    PROBLEM  F(X)=0, AND -2 FOR THE PROBLEM  RHO(A,LAMBDA,X)=0.
C    IN CERTAIN SITUATIONS  IFLAG  IS SET TO 2 OR 3 BY  FIXPDF,
C    AND  FIXPDF  CAN BE CALLED AGAIN WITHOUT CHANGING  IFLAG.
C
C ARCTOL  IS THE LOCAL ERROR ALLOWED THE ODE SOLVER WHEN
C    FOLLOWING THE ZERO CURVE.  IF  ARCTOL .LE. 0.0  ON INPUT
C    IT IS RESET TO  .5*DSQRT(EPS).  NORMALLY  ARCTOL  SHOULD
C    BE CONSIDERABLY LARGER THAN  EPS.
C
C EPS  IS THE LOCAL ERROR ALLOWED THE ODE SOLVER WHEN VERY
C    NEAR THE FIXED POINT(ZERO).  EPS  IS APPROXIMATELY THE
C    MIXED ABSOLUTE AND RELATIVE ERROR IN THE COMPUTED FIXED
C    POINT(ZERO).
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:NDIMA) 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  HAS LENGTH  NDIMA  AND MUST BE INITIALIZED
C    BY THE USER.
C
C NDIMA  IS THE DIMENSION OF  A, AND IS USED ONLY FOR THE CURVE
C    TRACKING PROBLEM.
C
C YP(1:N+1) IS A WORK ARRAY CONTAINING THE CURRENT TANGENT
C    VECTOR TO THE ZERO CURVE.
C
C YPOLD(1:N+1) IS A WORK ARRAY CONTAINING THE PREVIOUS TANGENT
C    VECTOR TO THE ZERO CURVE.
C
C QR(1:N,1:N+1), ALPHA(1:N), TZ(1:N+1), AND PIVOT(1:N+1) ARE
C    ALL WORK ARRAYS USED BY  FODE  TO CALCULATE THE TANGENT
C    VECTOR YP.
C
C WT(1:N+1), PHI(1:N+1,1:16), AND P(1:N+1) ARE ALL WORK ARRAYS
C    USED BY THE ODE SUBROUTINE  STEPS  .
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    RHOA, RHOJAC.
C
C Y, ARCTOL, EPS, ARCLEN, NFE, AND IFLAG SHOULD ALL BE
C VARIABLES IN THE CALLING PROGRAM.
C
C
C ON OUTPUT:
C
C N  AND  TRACE  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(ZERO) OF F(X).  IN ABNORMAL SITUATIONS LAMBDA
C    MAY ONLY BE NEAR 1 AND X IS NEAR A FIXED POINT(ZERO).
C
C IFLAG =
C  -2   CAUSES  FIXPDF  TO INITIALIZE EVERYTHING FOR THE PROBLEM
C       RHO(A,LAMBDA,X) = 0 (USE ON FIRST CALL).
C
C  -1   CAUSES  FIXPDF  TO INITIALIZE EVERYTHING FOR THE PROBLEM
C       F(X) = 0 (USE ON FIRST CALL).
C
C   0   CAUSES  FIXPDF  TO INITIALIZE EVERYTHING FOR THE PROBLEM
C       X = F(X) (USE ON FIRST CALL).
C
C   1   NORMAL RETURN.
C
C   2   SPECIFIED ERROR TOLERANCE CANNOT BE MET.  EPS HAS BEEN
C       INCREASED TO A SUITABLE VALUE.  TO CONTINUE, JUST CALL
C       FIXPDF  AGAIN WITHOUT CHANGING ANY PARAMETERS.
C
C   3   STEPS  HAS BEEN CALLED 1000 TIMES.  TO CONTINUE, CALL
C       FIXPDF  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   EPS  (OR  ARCTOL ) IS TOO LARGE.  THE PROBLEM SHOULD BE
C       RESTARTED BY CALLING  FIXPDF  WITH A SMALLER  EPS  (OR
C       ARCTOL ) AND  IFLAG = 0 (-1, -2).
C
C   6   I - DF(X)  IS NEARLY SINGULAR AT THE FIXED POINT (DF(X) IS
C       NEARLY SINGULAR AT THE ZERO, OR  D RHO(A,LAMBDA,X)/DX  IS
C       NEARLY SINGULAR AT  LAMBDA = 1 ).  ANSWER MAY NOT BE
C       ACCURATE.
C
C   7   ILLEGAL INPUT PARAMETERS, A FATAL ERROR.
C
C ARCTOL = EPS AFTER A NORMAL RETURN (IFLAG = 1).
C
C EPS  IS UNCHANGED AFTER A NORMAL RETURN (IFLAG = 1).  IT IS
C    INCREASED TO AN APPROPRIATE VALUE ON THE RETURN IFLAG = 2.
C
C A  WILL (NORMALLY) HAVE BEEN MODIFIED.
C
C NFE  IS THE NUMBER OF FUNCTION EVALUATIONS (= NUMBER OF
C    JACOBIAN EVALUATIONS).
C
C ARCLEN  IS THE LENGTH OF THE PATH FOLLOWED.
C
C
      DOUBLE PRECISION AOLD,ARCLEN,ARCTOL,CURSW,CURTOL,EPS,
     1  EPSSTP,EPST,H,HOLD,S,S99,SA,SB,SOUT,SQNP1,XOLD,Y1SOUT
      INTEGER IFLAG,IFLAGC,ITER,IVC,J,JUDY,JW,K,KGI,KOLD,
     1  KSTEPS,LCODE,LIMIT,LIMITD,N,NDIMA,NFE,NFEC,NP1,TRACE
      LOGICAL START,CRASH,ST99
C
C *****  ARRAY DECLARATIONS.  *****
C
C ARRAYS NEEDED BY THE ODE SUBROUTINE  STEPS .
      DOUBLE PRECISION Y(N+1),WT(N+1),PHI(N+1,16),P(N+1),YP(N+1),
     1     ALPHAS(12),W(12),G(13),GI(11)
      INTEGER IV(10)
C
C ARRAYS NEEDED BY  FIXPDF , FODE , AND  DCPOSE .
      DOUBLE PRECISION YPOLD(N+1),A(N),QR(N,N+1),ALPHA(N),TZ(N+1),
     $     PAR(1)
      INTEGER PIVOT(N+1),IPAR(1)
C
C *****  END OF DIMENSIONAL INFORMATION.  *****
C
      SAVE
      EXTERNAL FODE
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
C :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :
      IF (N .LE. 0  .OR.  EPS .LE. 0.0 ) IFLAG=7
      IF (IFLAG .GE. -2  .AND.  IFLAG .LE. 0) GO TO 10
      IF (IFLAG .EQ. 2) GO TO 35
      IF (IFLAG .EQ. 3) GO TO 30
C ONLY VALID INPUT FOR  IFLAG  IS -2, -1, 0, 2, 3.
      IFLAG=7
      RETURN
C
C *****  INITIALIZATION BLOCK.  *****
C
10    ARCLEN=0.0
      S=0.0
      IF (ARCTOL .LE. 0.0) ARCTOL=.5*SQRT(EPS)
      NFEC=0
      IFLAGC=IFLAG
      NP1=N+1
      SQNP1=SQRT(DBLE(NP1))
C
C SWITCH FROM THE TOLERANCE  ARCTOL  TO THE (FINER) TOLERANCE  EPS  IF
C THE CURVATURE OF ANY COMPONENT OF  Y  EXCEEDS  CURSW.
C
      CURSW=10.0
C
      ST99=.FALSE.
      START=.TRUE.
      CRASH=.FALSE.
      HOLD=1.0
      H=.1
      EPSSTP=ARCTOL
      KSTEPS=0
C SET INITIAL CONDITIONS FOR ORDINARY DIFFERENTIAL EQUATION.
      YPOLD(1)=1.0
      YP(1)=1.0
      Y(1)=0.0
      DO 20 J=2,NP1
        YPOLD(J)=0.0
        YP(J)=0.0
20    CONTINUE
C LOAD  A  FOR THE FIXED POINT AND ZERO FINDING PROBLEMS.
      IF (IFLAGC .GE. -1) THEN
        DO 23 J=2,NP1
          A(J-1)=Y(J)
23      CONTINUE
      ENDIF
30    LIMIT=LIMITD
C
C *****  END OF INITIALIZATION BLOCK.  *****
C
C
C *****  MAIN LOOP.  *****
C
35    DO 150 ITER=1,LIMIT
      IF (Y(1) .LT. 0.0) THEN
40      ARCLEN=ARCLEN+S
        IFLAG=5
        RETURN
      ENDIF
50    IF (S .LE. 7.0*SQNP1) GO TO 80
C ARC LENGTH IS GETTING TOO LONG, THE PROBLEM WILL BE
C RESTARTED WITH A DIFFERENT  A  VECTOR.
      ARCLEN=ARCLEN+S
      S=0.0
60    START=.TRUE.
      CRASH=.FALSE.
C COMPUTE A NEW  A  VECTOR.
      IF (IFLAGC .EQ. -2) THEN
        DO 63 JW=1,NDIMA
          QR(JW,1)=A(JW)
63      CONTINUE
        CALL RHOA(A,Y(1),Y(2),PAR,IPAR)
        DO 65 JW=1,NDIMA
          AOLD=QR(JW,1)
C IF NEW AND OLD  A  DIFFER BY TOO MUCH, TRACKING SHOULD NOT CONTINUE.
          IF (ABS(A(JW)-AOLD) .GT. 1.0+ABS(AOLD)) THEN
            ARCLEN=ARCLEN+S
            IFLAG=5
            RETURN
          ENDIF
65      CONTINUE
      ELSE
        CALL F(Y(2),YP)
        DO 70 JW=1,N
          AOLD=A(JW)
          IF (IFLAGC .EQ. -1) THEN
            A(JW)=Y(1)*YP(JW)/(1.0 - Y(1)) + Y(JW+1)
          ELSE
            A(JW)=(Y(JW+1) - Y(1)*YP(JW))/(1.0 - Y(1))
          ENDIF
C IF NEW AND OLD  A  DIFFER BY TOO MUCH, TRACKING SHOULD NOT CONTINUE.
          IF (ABS(A(JW)-AOLD) .GT. 1.0+ABS(AOLD)) THEN
            ARCLEN=ARCLEN+S
            IFLAG=5
            RETURN
          ENDIF
70      CONTINUE
      ENDIF
      GO TO 100
80    IF (Y(1) .LE. .99  .OR. ST99) GO TO 100
C WHEN LAMBDA REACHES .99, THE PROBLEM WILL BE RESTARTED WITH
C A NEW  A  VECTOR.
90    ST99=.TRUE.
      EPSSTP=EPS
      ARCTOL=EPS
      GO TO 60
C
C SET DIFFERENT ERROR TOLERANCE FOR HIGH CURVATURE COMPONENTS OF THE
C TRAJECTORY Y(S).
100   CURTOL=CURSW*HOLD
      EPST=EPS/EPSSTP
      DO 110 JW=1,NP1
        IF (ABS(YP(JW)-YPOLD(JW)) .LE. CURTOL) THEN
          WT(JW)=(ABS(Y(JW))+1.0)
        ELSE
          WT(JW)=(ABS(Y(JW))+1.0)*EPST
        ENDIF
110   CONTINUE
C
C TAKE A STEP ALONG THE CURVE.
      CALL STEPS(FODE,NP1,Y,S,H,EPSSTP,WT,START,HOLD,K,KOLD,CRASH,
     +     PHI,P,YP,ALPHAS,W,G,KSTEPS,XOLD,IVC,IV,KGI,GI,
     +     YPOLD,A,QR,ALPHA,TZ,PIVOT,NFEC,IFLAGC,PAR,IPAR)
C PRINT LATEST POINT ON CURVE IF REQUESTED.
      IF (TRACE .GT. 0) THEN
        WRITE (TRACE,117) ITER,NFEC,S,Y(1),(Y(JW),JW=2,NP1)
117     FORMAT(/' STEP',I5,3X,'NFE =',I5,3X,'ARC LENGTH =',F9.4,3X,
     $  'LAMBDA =',F7.4,5X,'X vector:'/1P,(1X,6E12.4))
      ENDIF
      NFE=NFEC
C CHECK IF THE STEP WAS SUCCESSFUL.
      IF (IFLAGC .EQ. 4) THEN
        ARCLEN=ARCLEN+S
        IFLAG=4
        RETURN
      ENDIF
120   IF (CRASH) THEN
C RETURN CODE FOR ERROR TOLERANCE TOO SMALL.
        IFLAG=2
C CHANGE ERROR TOLERANCES.
        EPS=EPSSTP
        IF (ARCTOL .LT. EPS) ARCTOL=EPS
C CHANGE LIMIT ON NUMBER OF ITERATIONS.
        LIMIT=LIMIT-ITER
        RETURN
      ENDIF
C
130   IF (Y(1) .GE. 1.0) THEN
        IF (ST99) GO TO 160
C
C IF LAMBDA .GE. 1.0 BUT THE PROBLEM HAS NOT BEEN RESTARTED
C WITH A NEW  A  VECTOR, BACK UP AND RESTART.
C
        S99=S-.5*HOLD
C GET AN APPROXIMATE ZERO Y(S) WITH  Y(1)=LAMBDA .LT. 1.0  .
135     CALL SINTRP(S,Y,S99,WT,YP,NP1,KOLD,PHI,IVC,IV,KGI,GI,
     $     ALPHAS,G,W,XOLD,P)
        IF (WT(1) .LT. 1.0) GO TO 140
        S99=.5*(S-HOLD+S99)
        GO TO 135
C
140     DO 144 JUDY=1,NP1
          Y(JUDY)=WT(JUDY)
          YPOLD(JUDY)=YP(JUDY)
144     CONTINUE
        S=S99
        GO TO 90
      ENDIF
C
150   CONTINUE
C
C *****  END OF MAIN LOOP.  *****
C
C LAMBDA HAS NOT REACHED 1 IN 1000 STEPS.
      IFLAG=3
      RETURN
C
C
C USE INVERSE INTERPOLATION TO GET THE ANSWER AT LAMBDA = 1.0 .
C
160   SA=S-HOLD
      SB=S
      LCODE=1
170   CALL ROOT(SOUT,Y1SOUT,SA,SB,EPS,EPS,LCODE)
C ROOT  FINDS S SUCH THAT Y(1)(S) = LAMBDA = 1 .
      IF (LCODE .GT. 0) GO TO 190
      CALL SINTRP(S,Y,SOUT,WT,YP,NP1,KOLD,PHI,IVC,IV,KGI,GI,
     $     ALPHAS,G,W,XOLD,P)
      Y1SOUT=WT(1)-1.0
      GO TO 170
190   IFLAG=1
C SET IFLAG = 6 IF  ROOT  COULD NOT GET  LAMBDA = 1.0  .
      IF (LCODE .GT. 2) IFLAG=6
      ARCLEN=ARCLEN+SA
C LAMBDA(SA) = 1.0 .
      CALL SINTRP(S,Y,SA,WT,YP,NP1,KOLD,PHI,IVC,IV,KGI,GI,
     $     ALPHAS,G,W,XOLD,P)
C
      DO 210 J=1,NP1
210   Y(J)=WT(J)
      RETURN
      END
