SUBROUTINE FIXPDS(N,Y,IFLAG,ARCTOL,EPS,TRACE,A,NDIMA,NFE,
& ARCLEN,MODE,LENQR)
C
C Subroutine FIXPDS 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,X,LAMBDA). 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) = (X(S), LAMBDA(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,X,LAMBDA) is assumed to
C be a C2 map from E**M X E**N X [0,1) 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
C from LAMBDA = 0, X = X0 is tracked until LAMBDA = 1 by
C solving the ordinary differential equation
C D RHO(A,X(S),LAMBDA(S))/DS = 0 for Y(S) = (X(S), LAMBDA(S)),
C where S is arc length along the zero curve. Also the homotopy
C map RHO(A,X,LAMBDA) 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
C FJACS(X) which evaluates, if
C MODE = 1,
C the (symmetric) Jacobian matrix of F(X) at X, and returns the
C symmetric Jacobian matrix in packed skyline storage format in
C QR, or if
C MODE = 2,
C returns the (nonsymmetric) Jacobian matrix in sparse row format
C in QR. The MODE 1 format is defined by QR, LENQR, ROWPOS; the
C MODE 2 format is defined by QR, LENQR, ROWPOS, COLPOS.
C
C For the curve tracking problem, the user must supply a subroutine
C RHOA(V,LAMBDA,X) which given (X,LAMBDA) returns a
C parameter vector A in V such that RHO(A,X,LAMBDA)=0, and a
C subroutine RHOJS(A,LAMBDA,X) which, if
C MODE = 1,
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). This data structure is described by QR, LENQR,
C ROWPOS, PP. *** Note the minus sign in the definition of PP. *** If
C MODE = 2,
C the (nonsymmetric) N X (N+1) Jacobian matrix [D RHO/DX, D RHO/DLAMBDA]
C evaluated at (A,X,LAMBDA) is returned in a data structure described
C by QR, LENQR, ROWPOS, COLPOS.
C
C Whichever of the routines F, FJACS, RHOA, RHOJS are required
C should be supplied as external subroutines, conforming with the
C interfaces in the module HOMOTOPY.
C
C
C Functions and subroutines directly or indirectly called by FIXPDS:
C DLAIC1 and DLAMCH (LAPACK), F (or RHOA ), FJACS (or RHOJS ),
C FODEDS , GMFADS , GMRES , GMRILUDS , ILUFDS , ILUSOLVDS , MULTDS ,
C MULT2DS , PCGDS , ROOT , SINTRP , SOLVDS , STEPDS , and the BLAS
C functions DDOT , DNRM2. The module REAL_PRECISION specifies 64-bit
C real arithmetic, which the user may want to change.
C
C ***Warning: this subroutine is generally more robust than FIXPNS
C and FIXPQS, 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,X,LAMBDA).
C
C Y is an array of length N + 1. (Y(1),...,Y(N)) = A is the
C starting point for the zero curve for the fixed point and
C zero finding problems. (Y(1),...,Y(N)) = 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 FIXPDS for the problem X=F(X), -1 for the
C problem F(X)=0, and -2 for the problem RHO(A,X,LAMBDA)=0.
C In certain situations IFLAG is set to 2 or 3 by FIXPDS,
C and FIXPDS 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 MODE = 1 if the Jacobian matrix is symmetric and stored in a packed
C skyline format;
C = 2 if the Jacobian matrix is stored in a sparse row format.
C
C LENQR is the number of nonzero entries in the sparse Jacobian
C matrices, used to determine the sparse matrix data structures.
C
C A, 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),...,Y(N)) = X, Y(N+1) = LAMBDA, 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 FIXPDS to initialize everything for the problem
C RHO(A,X,LAMBDA) = 0 (use on first call).
C
C -1 causes FIXPDS to initialize everything for the problem
C F(X) = 0 (use on first call).
C
C 0 causes FIXPDS 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 FIXPDS again without changing any parameters.
C
C 3 STEPDS has been called 1000 times. To continue, call
C FIXPDS again without changing any parameters.
C
C 4 Jacobian matrix does not have full rank or has a zero on the
C diagonal, and/or the conjugate gradient iteration for the
C kernel of the Jacobian matrix failed to converge. The
C algorithm has failed (the zero curve of the homotopy map
C cannot be followed any further).
C
C 5 EPS (or ARCTOL ) is too large. The problem should be
C restarted by calling FIXPDS 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,X,LAMBDA)/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
C Allocatable and 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:LENQR), PP(1:N), ROWPOS(1:N+2), COLPOS(1:LENQR) are all work
C arrays used to define the sparse Jacobian matrices, allocated
C here, and distributed via the module HOMOTOPY .
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 STEPDS .
C
USE HOMOTOPY, QR => QRSPARSE
USE REAL_PRECISION
INTEGER, INTENT(IN)::LENQR,MODE,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. *****
C
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 ARRAYS NEEDED BY THE ODE SUBROUTINE STEPDS .
REAL (KIND=R8), SAVE:: ALPHAS(12),G(13),GI(11),W(12)
REAL (KIND=R8), ALLOCATABLE, SAVE:: P(:),PHI(:,:),WT(:),YP(:)
INTEGER, SAVE:: IV(10)
C
C ARRAYS NEEDED BY FIXPDS , FODEDS , AND PCGDS .
REAL (KIND=R8), ALLOCATABLE, DIMENSION(:), SAVE:: AOLD,YPOLD
C
C ***** END OF DIMENSIONAL INFORMATION. *****
C
INTERFACE
SUBROUTINE FODEDS(S,Y,YP,N,IFLAG,YPOLD,A,NDIMA,LENQR,MODE,NFE)
USE HOMOTOPY, QR => QRSPARSE
USE REAL_PRECISION
INTEGER:: IFLAG,LENQR,MODE,N,NDIMA,NFE
REAL (KIND=R8):: A(NDIMA),S,Y(N+1),YP(N+1),YPOLD(N+1)
END SUBROUTINE FODEDS
END INTERFACE
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
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))
& .OR. MODE .LE. 0 .OR. MODE .GE. 3)
& 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(REAL(NP1,KIND=R8))
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 ALLOCATE SAVED WORK ARRAYS.
IF (ALLOCATED(AOLD)) DEALLOCATE(AOLD)
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(AOLD(NDIMA),P(NP1),PHI(NP1,16),WT(NP1),YP(NP1),
& YPOLD(NP1))
C SET INITIAL CONDITIONS FOR ORDINARY DIFFERENTIAL EQUATION.
YPOLD(NP1)=1.0
YP(NP1)=1.0
Y(NP1)=0.0
YPOLD(1:N)=0.0
YP(1:N)=0.0
C LOAD A FOR THE FIXED POINT AND ZERO FINDING PROBLEMS.
IF (IFLAGC .GE. -1) THEN
A=Y(1:N)
ENDIF
30 LIMIT=LIMITD
C ALLOCATE ARRAYS FOR SPARSE JACOBIAN MATRIX DATA STRUCTURE.
35 SELECT CASE (MODE)
CASE (1)
IF (.NOT. ALLOCATED(QR)) ALLOCATE(QR(LENQR))
IF (.NOT. ALLOCATED(ROWPOS)) ALLOCATE(ROWPOS(N+2))
IF (.NOT. ALLOCATED(PP)) ALLOCATE(PP(N))
CASE (2)
IF (.NOT. ALLOCATED(QR)) ALLOCATE(QR(LENQR))
IF (.NOT. ALLOCATED(ROWPOS)) ALLOCATE(ROWPOS(N+2))
IF (.NOT. ALLOCATED(COLPOS)) ALLOCATE(COLPOS(LENQR))
IF ((.NOT. ALLOCATED(PP)) .AND. (IFLAGC .GE. -1))
& ALLOCATE(PP(N))
END SELECT
C
C ***** END OF INITIALIZATION BLOCK. *****
C
MAIN_LOOP: DO ITER=1,LIMIT ! ***** MAIN LOOP. *****
IF (Y(NP1) .LT. 0.0) THEN
ARCLEN=ARCLEN+S
IFLAG=5
CALL CLEANUPALL
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(NP1),Y(1:N))
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 CLEANUPALL
RETURN
ENDIF
ELSE
CALL F(Y(1:N),YP(1:N))
AOLD=A
IF (IFLAGC .EQ. -1) THEN
A=Y(NP1)*YP(1:N)/(1.0 - Y(NP1)) + Y(1:N)
ELSE
A=(Y(1:N) - Y(NP1)*YP(1:N))/(1.0 - Y(NP1))
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 CLEANUPALL
RETURN
ENDIF
ENDIF
GO TO 100
80 IF (Y(NP1) .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 STEPDS(FODEDS,NP1,Y,S,H,EPSSTP,WT,START,HOLD,K,KOLD,CRASH,
& PHI,P,YP,ALPHAS,W,G,KSTEPS,XOLD,IVC,IV,KGI,GI,
& IFLAGC,YPOLD,A,NDIMA,LENQR,MODE,NFEC)
C PRINT LATEST POINT ON CURVE IF REQUESTED.
IF (TRACE .GT. 0) THEN
WRITE (TRACE,117) ITER,NFEC,S,Y(NP1),(Y(JW),JW=1,N)
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 CLEANUPALL
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
CALL CLEANUP
RETURN
ENDIF
C
IF (Y(NP1) .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(NP1)=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(NP1) .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
CALL CLEANUP
RETURN
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(NP1)(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(NP1)-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 CLEANUPALL
RETURN
C
CONTAINS
SUBROUTINE CLEANUPALL
IF (ALLOCATED(AOLD)) DEALLOCATE(AOLD)
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)
CALL CLEANUP
RETURN
END SUBROUTINE CLEANUPALL
SUBROUTINE CLEANUP
IF (ALLOCATED(QR)) DEALLOCATE(QR)
IF (ALLOCATED(ROWPOS)) DEALLOCATE(ROWPOS)
IF (ALLOCATED(COLPOS)) DEALLOCATE(COLPOS)
IF (ALLOCATED(PP)) DEALLOCATE(PP)
IF (ALLOCATED(PAR)) DEALLOCATE(PAR)
IF (ALLOCATED(IPAR)) DEALLOCATE(IPAR)
RETURN
END SUBROUTINE CLEANUP
END SUBROUTINE FIXPDS