SUBROUTINE FIXPQF(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,A,
& SSPAR,NFE,ARCLEN)
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) 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) 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 F (or RHO),
C FJAC (or RHOJAC), ROOT, ROOTQF, STEPQF, TANGQF, UPQRQF,
C the LAPACK routines DGEQRF, DORGQR, their auxiliary routines,
C and the BLAS routines DCOPY, DDOT, DGEMM, DGEMV,
C DGER, DNRM2, DSCAL, DTPMV, DTPSV, DTRMM, and DTRMV.
C The module REAL_PRECISION specifies 64-bit real arithmetic,
C which the user may want to change.
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(:) 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 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
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 8 Memory allocation error, fatal.
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
C Allocatable and automatic work arrays:
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 Q(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
C ***** DECLARATIONS *****
USE REAL_PRECISION
C
C FUNCTION DECLARATIONS
C
INTERFACE
FUNCTION DNRM2(N,X,STRIDE)
USE REAL_PRECISION
INTEGER:: N,STRIDE
REAL (KIND=R8):: DNRM2,X(N)
END FUNCTION DNRM2
END INTERFACE
C
C LOCAL VARIABLES
C
REAL (KIND=R8), SAVE:: ABSERR, H, HOLD, RELERR, S, WK
INTEGER, SAVE:: IFLAGC, ITER, JW, LIMIT, NP1
LOGICAL, SAVE:: CRASH, START
C
C SCALAR ARGUMENTS
C
REAL (KIND=R8):: ARCRE, ARCAE, ANSRE, ANSAE, ARCLEN
INTEGER:: N,IFLAG,TRACE,NFE
C
C ARRAY DECLARATIONS
C
REAL (KIND=R8), DIMENSION(:), ALLOCATABLE, SAVE::
& R,YOLD,YP,YPOLD
REAL (KIND=R8), DIMENSION(:,:), ALLOCATABLE, SAVE:: Q
REAL (KIND=R8):: A(:), DZ(N+1), F0(N+1), F1(N+1),
& SSPAR(4), T(N+1), W(N+1), Y(:), YSAV(N+1), Z0(N+1)
C
C ***** END OF DECLARATIONS *****
INTERFACE
SUBROUTINE STEPQF(N,NFE,IFLAG,START,CRASH,HOLD,H,
& WK,RELERR,ABSERR,S,Y,YP,YOLD,YPOLD,A,Q,R,
& F0,F1,Z0,DZ,W,T,SSPAR)
USE REAL_PRECISION
INTEGER:: N, NFE, IFLAG
LOGICAL:: START, CRASH
REAL (KIND=R8):: HOLD, H, WK, RELERR, ABSERR, S
REAL (KIND=R8):: A(:), DZ(N+1), F0(N+1), F1(N+1),
& Q(N+1,N+1), R((N+1)*(N+2)/2), SSPAR(4), T(N+1), W(N+1),
& Y(:), YOLD(N+1), YP(N+1), YPOLD(N+1), Z0(N+1)
END SUBROUTINE STEPQF
SUBROUTINE ROOTQF(N,NFE,IFLAG,RELERR,ABSERR,Y,YP,YOLD,
& YPOLD,A,Q,R,DZ,Z,W,T,F0,F1)
USE REAL_PRECISION
REAL (KIND=R8):: RELERR, ABSERR
INTEGER:: N, NFE, IFLAG
REAL (KIND=R8):: A(:), DZ(N+1), F0(N+1), F1(N+1),
& Q(N+1,N+1), R((N+1)*(N+2)/2), T(N+1), W(N+1),
& Y(:), YOLD(N+1), YP(N+1), YPOLD(N+1), Z(N+1)
END SUBROUTINE ROOTQF
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 ***** FIRST EXECUTABLE STATEMENT *****
C
C CHECK IFLAG
C
IF (N .LE. 0 .OR. ANSRE .LE. 0.0 .OR. ANSAE .LT. 0.0
& .OR. N+1 .NE. SIZE(Y) .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 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
IF (ALLOCATED(Q)) DEALLOCATE(Q)
IF (ALLOCATED(R)) DEALLOCATE(R)
IF (ALLOCATED(YOLD)) DEALLOCATE(YOLD)
IF (ALLOCATED(YP)) DEALLOCATE(YP)
IF (ALLOCATED(YPOLD)) DEALLOCATE(YPOLD)
ALLOCATE(Q(NP1,NP1),R(NP1*(N+2)/2),YOLD(NP1),YP(NP1),YPOLD(NP1),
& STAT=JW)
IF (JW /= 0) THEN
IFLAG=8
RETURN
END IF
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
YPOLD(2:NP1)=0.0
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)*EPSILON(1.0_R8)
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_R8
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
A=Y(2:NP1)
ENDIF
C
40 LIMIT=LIMITD
C
C ***** END OF INITIALIZATION BLOCK. *****
C
50 DO ITER=1,LIMIT ! ***** MAIN LOOP. *****
IF (Y(1) .LT. 0.0) THEN
ARCLEN = S
IFLAG = 5
CALL CLEANUP ; 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,Q,R,F0,F1,Z0,DZ,
& W,T,SSPAR)
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:'/(1X,6ES12.4))
ENDIF
C
C CHECK IF THE STEP WAS SUCCESSFUL.
C
IF (IFLAGC .GT. 0) THEN
ARCLEN=S
IFLAG=IFLAGC
CALL CLEANUP ; 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
END DO ! ***** 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 YSAV=YOLD
C
C FIND SOLUTION.
C
CALL ROOTQF(N,NFE,IFLAGC,ANSRE,ANSAE,Y,YP,YOLD,
& YPOLD,A,Q,R,DZ,Z0,W,T,F0,F1)
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
DZ = Y - YSAV
ARCLEN = S - HOLD + DNRM2(NP1,DZ,1)
C
C ***** END OF FINAL STEP *****
C
CALL CLEANUP ; RETURN
C
CONTAINS
SUBROUTINE CLEANUP
IF (ALLOCATED(Q)) DEALLOCATE(Q)
IF (ALLOCATED(R)) DEALLOCATE(R)
IF (ALLOCATED(YOLD)) DEALLOCATE(YOLD)
IF (ALLOCATED(YP)) DEALLOCATE(YP)
IF (ALLOCATED(YPOLD)) DEALLOCATE(YPOLD)
END SUBROUTINE CLEANUP
END SUBROUTINE FIXPQF