SUBROUTINE FIXPDF(N,Y,IFLAG,ARCTOL,EPS,TRACE,A,NDIMA,
& NFE,ARCLEN)
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) 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)
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 Whichever of the routines F, FJAC, RHOA, RHOJAC are required
C should be supplied as external subroutines, conforming with the
C interfaces in the module HOMOTOPY.
C FIXPDF directly or indirectly uses the subroutines
C F (or RHOA ), FJAC (or RHOJAC ), FODE, ROOT,
C SINTRP, STEPS, the LAPACK routine DGEQPF, auxiliary LAPACK
C routines, and the BLAS functions DCOPY, DDOT, DGEMV, DGER,
C DNRM2, DSCAL, DSWAP, IDAMAX.
C The module REAL_PRECISION specifies 64-bit real arithmetic, which
C the user may want to change.
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 array 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*SQRT(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, used for the curve tracking problem,
C and must be N for the fixed point and zero finding problems.
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 8 Memory allocation error, fatal.
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 matrix evaluations).
C
C ARCLEN is the length of the path followed.
C
C
C Automatic work arrays:
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:3*N+3), 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
USE HOMOTOPY
USE REAL_PRECISION
C
INTEGER, INTENT(IN)::N,NDIMA,TRACE
REAL (KIND=R8), DIMENSION(:), INTENT(IN OUT)::A,Y
INTEGER, INTENT(IN OUT)::IFLAG
REAL (KIND=R8), INTENT(IN OUT)::ARCTOL,EPS
INTEGER, INTENT(OUT)::NFE
REAL (KIND=R8), INTENT(OUT)::ARCLEN
C
C LOCAL VARIABLES.
REAL (KIND=R8), SAVE:: CURSW,CURTOL,EPSSTP,EPST,H,HOLD,
& S,S99,SA,SB,SOUT,SQNP1,XOLD,Y1SOUT
INTEGER, SAVE:: IFLAGC,ITER,IVC,JW,K,KGI,KOLD,
& KSTEPS,LCODE,LIMIT,NFEC,NP1
LOGICAL, SAVE:: CRASH,START,ST99
C
C ***** ARRAY DECLARATIONS. *****
C
C ARRAYS NEEDED BY THE ODE SUBROUTINE STEPS .
REAL (KIND=R8), ALLOCATABLE, SAVE:: P(:),PHI(:,:),WT(:),YP(:)
REAL (KIND=R8), SAVE:: ALPHAS(12),G(13),GI(11),W(12)
INTEGER, SAVE:: IV(10)
C
C ARRAYS NEEDED BY FIXPDF , FODE , AND LAPACK ROUTINES.
REAL (KIND=R8), DIMENSION(:), ALLOCATABLE, SAVE:: YPOLD
REAL (KIND=R8):: ALPHA(3*N+3),AOLD(NDIMA),QR(N,N+1),TZ(N+1)
INTEGER:: PIVOT(N+1)
C
C ***** END OF DIMENSIONAL INFORMATION. *****
C
C LIMITD IS AN UPPER BOUND ON THE NUMBER OF STEPS. IT MAY BE
C CHANGED BY CHANGING THE FOLLOWING PARAMETER STATEMENT:
INTEGER, PARAMETER:: LIMITD=1000
C
INTERFACE
SUBROUTINE FODE(S,Y,YP,YPOLD,A,QR,ALPHA,TZ,PIVOT,NFE,N,IFLAG)
USE REAL_PRECISION
REAL (KIND=R8):: S
INTEGER:: IFLAG,N,NFE
REAL (KIND=R8):: A(:),Y(:),YP(N+1),YPOLD(N+1)
REAL (KIND=R8):: ALPHA(3*N+3),QR(N,N+1),TZ(N+1)
INTEGER, DIMENSION(N+1):: PIVOT
END SUBROUTINE FODE
SUBROUTINE STEPS(F,NEQN,Y,X,H,EPS,WT,START,HOLD,K,KOLD,CRASH,
& PHI,P,YP,ALPHA,W,G,KSTEPS,XOLD,IVC,IV,KGI,GI, FPWA1,FPWA2,
& FPWA3,FPWA4,FPWA5,IFPWA1,IFPC1,IFPC2)
USE REAL_PRECISION
EXTERNAL F
REAL (KIND=R8):: ALPHA,EPS,FPWA1,FPWA2,FPWA3,FPWA4,FPWA5,
& G,GI,H,HOLD,P,PHI,W,WT,X,XOLD,Y,YP
INTEGER:: IFPC1,IFPC2,IFPWA1,IV,IVC,K,KGI,KOLD,KSTEPS,NEQN
LOGICAL:: CRASH,START
DIMENSION Y(:),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),
& ALPHA(12),W(12),G(13),GI(11),IV(10),
& FPWA1(NEQN),FPWA2(:),FPWA3(NEQN-1,NEQN),FPWA4(3*NEQN),
& FPWA5(NEQN),IFPWA1(NEQN)
END SUBROUTINE STEPS
END INTERFACE
C
C
C : : : : : : : : : : : : : : : : : : : : :
IF (N .LE. 0 .OR. EPS .LE. 0.0 .OR. N+1 .NE. SIZE(Y)
& .OR. NDIMA .NE. SIZE(A) .OR.
& ((IFLAG .EQ. -1 .OR. IFLAG .EQ. 0) .AND. N .NE. SIZE(A)))
& 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))
IF (ALLOCATED(P)) DEALLOCATE(P)
IF (ALLOCATED(PHI)) DEALLOCATE(PHI)
IF (ALLOCATED(WT)) DEALLOCATE(WT)
IF (ALLOCATED(YP)) DEALLOCATE(YP)
IF (ALLOCATED(YPOLD)) DEALLOCATE(YPOLD)
ALLOCATE(P(NP1),PHI(NP1,16),WT(NP1),YP(NP1),YPOLD(NP1),
& STAT=JW)
IF (JW /= 0) THEN
IFLAG=8
RETURN
END IF
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
YPOLD(2:NP1)=0.0
YP(2:NP1)=0.0
C LOAD A FOR THE FIXED POINT AND ZERO FINDING PROBLEMS.
IF (IFLAGC .GE. -1) THEN
A=Y(2:NP1)
ENDIF
30 LIMIT=LIMITD
C
C ***** END OF INITIALIZATION BLOCK. *****
C
35 MAIN_LOOP: DO ITER=1,LIMIT ! ***** MAIN LOOP. *****
IF (Y(1) .LT. 0.0) THEN
ARCLEN=ARCLEN+S
IFLAG=5
CALL CLEANUP ; RETURN
ENDIF
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
AOLD=A
CALL RHOA(A,Y(1),Y(2:NP1))
C IF NEW AND OLD A DIFFER BY TOO MUCH, TRACKING SHOULD NOT CONTINUE.
IF (ANY(ABS(A-AOLD) .GT. 1.0+ABS(AOLD))) THEN
ARCLEN=ARCLEN+S
IFLAG=5
CALL CLEANUP ; RETURN
ENDIF
ELSE
CALL F(Y(2:NP1),YP(1:N))
AOLD=A
IF (IFLAGC .EQ. -1) THEN
A=Y(1)*YP(1:N)/(1.0 - Y(1)) + Y(2:NP1)
ELSE
A=(Y(2:NP1) - Y(1)*YP(1:N))/(1.0 - Y(1))
ENDIF
C IF NEW AND OLD A DIFFER BY TOO MUCH, TRACKING SHOULD NOT CONTINUE.
IF (ANY(ABS(A-AOLD) .GT. 1.0+ABS(AOLD))) THEN
ARCLEN=ARCLEN+S
IFLAG=5
CALL CLEANUP ; RETURN
ENDIF
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
WHERE (ABS(YP-YPOLD) .LE. CURTOL)
WT=(ABS(Y)+1.0)
ELSEWHERE
WT=(ABS(Y)+1.0)*EPST
END WHERE
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)
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:'/(1X,6ES12.4))
ENDIF
NFE=NFEC
C CHECK IF THE STEP WAS SUCCESSFUL.
IF (IFLAGC .EQ. 4) THEN
ARCLEN=ARCLEN+S
IFLAG=4
CALL CLEANUP ; RETURN
ENDIF
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
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 Y=WT
YPOLD=YP
S=S99
GO TO 90
ENDIF
C
END DO MAIN_LOOP ! ***** 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
Y=WT
CALL CLEANUP ; RETURN
C
CONTAINS
SUBROUTINE CLEANUP
IF (ALLOCATED(P)) DEALLOCATE(P)
IF (ALLOCATED(PHI)) DEALLOCATE(PHI)
IF (ALLOCATED(WT)) DEALLOCATE(WT)
IF (ALLOCATED(YP)) DEALLOCATE(YP)
IF (ALLOCATED(YPOLD)) DEALLOCATE(YPOLD)
END SUBROUTINE CLEANUP
END SUBROUTINE FIXPDF