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
!
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
!
SUBROUTINE FIXPNF(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,A,
& SSPAR,NFE,ARCLEN, POLY_SWITCH)
C
C Subroutine FIXPNF 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)) using a Hermite cubic predictor and a
C corrector which returns to the zero curve along the flow normal
C to the Davidenko flow (which consists of the integral curves of
C D(HOMOTOPY MAP)/DS ).
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
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 a
C subroutine RHOJAC(A,LAMBDA,X,V,K) which returns in V the Kth
C column of the N X (N+1) Jacobian matrix [D RHO/D LAMBDA, D RHO/DX]
C evaluated at (A,LAMBDA,X). FIXPNF directly or indirectly uses
C the subroutines F (or RHO ), FJAC (or RHOJAC ),
C ROOT, ROOTNF, STEPNF, the LAPACK routines DGEQPF, DORMQR,
C their auxiliary routines, and the BLAS routines DCOPY,
C DDOT, DGEMM, DGEMV, DGER, DNRM2, DSCAL, DSWAP, DTRMM, DTRMV,
C IDAMAX. The module REAL_PRECISION specifies 64-bit
C real arithmetic, 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(:) 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 FIXPNF 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 FIXPNF,
C and FIXPNF can be called again without changing IFLAG.
C
C ARCRE , ARCAE are the relative and absolute errors, respectively,
C allowed the normal flow 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 ||Z|| .LE. ANSRE*||X|| + ANSAE where
C
C (.,Z) is the 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:8) = (LIDEAL, RIDEAL, DIDEAL, HMIN, HMAX, BMIN, BMAX, P) is
C a vector of parameters used for the optimal step size estimation.
C If SSPAR(J) .LE. 0.0 on input, it is reset to a default value
C by FIXPNF . Otherwise the input value of SSPAR(J) is used.
C See the comments below and in STEPNF for more information about
C these constants.
C
C POLY_SWITCH is an optional logical variable used only by the driver
C POLSYS1H for polynomial systems.
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(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 FIXPNF to initialize everything for the problem
C RHO(A,LAMBDA,X) = 0 (use on first call).
C
C -1 causes FIXPNF to initialize everything for the problem
C F(X) = 0 (use on first call).
C
C 0 causes FIXPNF 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. Some or all of
C ARCRE , ARCAE , ANSRE , ANSAE have been increased to
C suitable values. To continue, just call FIXPNF again
C without changing any parameters.
C
C 3 STEPNF has been called 1000 times. To continue, call
C FIXPNF 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 tolerances
C ARC?E and ANS?E were too lenient. The problem should be
C restarted by calling FIXPNF with smaller error tolerances
C and IFLAG = 0 (-1, -2).
C
C 6 The normal flow Newton iteration in STEPNF or ROOTNF
C failed to converge. The error tolerances ANS?E may be too
C 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 function evaluations (= number of
C Jacobian matrix 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 tangent vector to
C the 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 QR(1:N,1:N+2), ALPHA(1:3*N+3), TZ(1:N+1), PIVOT(1:N+1), W(1:N+1),
C WP(1:N+1), Z0(1:N+1), Z1(1:N+1) are all work arrays used by
C STEPNF to calculate the tangent vectors and Newton steps.
C
C
USE REAL_PRECISION
INTEGER, INTENT(IN)::N,TRACE
REAL (KIND=R8), DIMENSION(:), INTENT(IN OUT)::A,Y
INTEGER, INTENT(IN OUT)::IFLAG
REAL (KIND=R8), INTENT(IN OUT)::ANSAE,ANSRE,ARCAE,ARCRE,
& SSPAR(8)
INTEGER, INTENT(OUT)::NFE
REAL (KIND=R8), INTENT(OUT)::ARCLEN
LOGICAL, INTENT(IN), OPTIONAL::POLY_SWITCH
C
C LOCAL VARIABLES.
REAL (KIND=R8), SAVE:: ABSERR,CURTOL,H,HOLD,RELERR,S
INTEGER, SAVE:: IFLAGC,ITER,JW,LIMIT,NC,NFEC,NP1
LOGICAL, SAVE:: CRASH,POLSYS,START
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 ALLOCATABLE AND AUTOMATIC ARRAYS.
REAL (KIND=R8), DIMENSION(:), ALLOCATABLE, SAVE:: YOLD,YP,YPOLD
REAL (KIND=R8):: ALPHA(3*N+3),QR(N,N+2),TZ(N+1),
& W(N+1),WP(N+1),Z0(N+1),Z1(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
C SWITCH FROM THE TOLERANCE ARC?E TO THE (FINER) TOLERANCE ANS?E IF
C THE CURVATURE OF ANY COMPONENT OF Y EXCEEDS CURSW.
REAL (KIND=R8), PARAMETER:: CURSW=10.0
C
INTERFACE
SUBROUTINE STEPNF(N,NFE,IFLAG,START,CRASH,HOLD,H,RELERR,
& ABSERR,S,Y,YP,YOLD,YPOLD,A,QR,ALPHA,TZ,PIVOT,W,WP,
& Z0,Z1,SSPAR)
USE REAL_PRECISION
REAL (KIND=R8):: ABSERR,H,HOLD,RELERR,S
INTEGER:: IFLAG,N,NFE
LOGICAL:: CRASH,START
REAL (KIND=R8):: A(:),ALPHA(3*N+3),QR(N,N+2),SSPAR(8),TZ(N+1),
& W(N+1),WP(N+1),Y(:),YOLD(N+1),YP(N+1),YPOLD(N+1),
& Z0(N+1),Z1(N+1)
INTEGER:: PIVOT(N+1)
END SUBROUTINE STEPNF
SUBROUTINE ROOTNF(N,NFE,IFLAG,RELERR,ABSERR,Y,YP,YOLD,
& YPOLD,A,QR,ALPHA,TZ,PIVOT,W,WP)
USE REAL_PRECISION
REAL (KIND=R8):: ABSERR,RELERR
INTEGER:: IFLAG,N,NFE
REAL (KIND=R8):: A(:),ALPHA(3*N+3),QR(N,N+2),TZ(N+1),W(N+1),
& WP(N+1),Y(:),YOLD(N+1),YP(N+1),YPOLD(N+1)
INTEGER:: PIVOT(N+1)
END SUBROUTINE ROOTNF
END INTERFACE
C
C
C : : : : : : : : : : : : : : : : : : : : : : : :
C TEST LOGICAL SWITCH TO REFLECT INTENDED USAGE OF FIXPNF.
IF (PRESENT(POLY_SWITCH)) THEN
POLSYS=.TRUE.
ELSE
POLSYS=.FALSE.
ENDIF
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 20
IF (IFLAG .EQ. 2) GO TO 120
IF (IFLAG .EQ. 3) GO TO 90
C ONLY VALID INPUT FOR IFLAG IS -2, -1, 0, 2, 3.
IFLAG=7
RETURN
C
C ***** INITIALIZATION BLOCK. *****
C
20 ARCLEN=0.0
IF (ARCRE .LE. 0.0) ARCRE=.5*SQRT(ANSRE)
IF (ARCAE .LE. 0.0) ARCAE=.5*SQRT(ANSAE)
NC=N
NFEC=0
IFLAGC=IFLAG
NP1=N+1
IF (ALLOCATED(YP)) DEALLOCATE(YP)
IF (ALLOCATED(YOLD)) DEALLOCATE(YOLD)
IF (ALLOCATED(YPOLD)) DEALLOCATE(YPOLD)
ALLOCATE(YP(NP1),YOLD(NP1),YPOLD(NP1),STAT=JW)
IF (JW /= 0) THEN
IFLAG=8
RETURN
END IF
C SET INITIAL CONDITIONS FOR FIRST CALL TO STEPNF .
START=.TRUE.
CRASH=.FALSE.
HOLD=1.0
H=.1
S=0.0
YPOLD(1)=1.0
YP(1)=1.0
Y(1)=0.0
YPOLD(2:NP1)=0.0
YP(2:NP1)=0.0
C SET OPTIMAL STEP SIZE ESTIMATION PARAMETERS.
C LET Z[K] DENOTE THE NEWTON ITERATES ALONG THE FLOW NORMAL TO THE
C DAVIDENKO FLOW AND Y THEIR LIMIT.
C IDEAL CONTRACTION FACTOR: ||Z[2] - Z[1]|| / ||Z[1] - Z[0]||
IF (SSPAR(1) .LE. 0.0) SSPAR(1)= .5
C IDEAL RESIDUAL FACTOR: ||RHO(A, Z[1])|| / ||RHO(A, Z[0])||
IF (SSPAR(2) .LE. 0.0) SSPAR(2)= .01
C IDEAL DISTANCE FACTOR: ||Z[1] - Y|| / ||Z[0] - Y||
IF (SSPAR(3) .LE. 0.0) SSPAR(3)= .5
C MINIMUM STEP SIZE HMIN .
IF (SSPAR(4) .LE. 0.0) SSPAR(4)=(SQRT(N+1.0)+4.0)*EPSILON(1.0_R8)
C MAXIMUM STEP SIZE HMAX .
IF (SSPAR(5) .LE. 0.0) SSPAR(5)= 1.0
C MINIMUM STEP SIZE REDUCTION FACTOR BMIN .
IF (SSPAR(6) .LE. 0.0) SSPAR(6)= .1_R8
C MAXIMUM STEP SIZE EXPANSION FACTOR BMAX .
IF (SSPAR(7) .LE. 0.0) SSPAR(7)= 3.0
C ASSUMED OPERATING ORDER P .
IF (SSPAR(8) .LE. 0.0) SSPAR(8)= 2.0
C
C LOAD A FOR THE FIXED POINT AND ZERO FINDING PROBLEMS.
IF (IFLAGC .GE. -1) THEN
A=Y(2:NP1)
ENDIF
90 LIMIT=LIMITD
C
C ***** END OF INITIALIZATION BLOCK. *****
C
120 MAIN_LOOP: DO ITER=1,LIMIT ! ***** MAIN LOOP. *****
IF (Y(1) .LT. 0.0) THEN
ARCLEN=S
IFLAG=5
CALL CLEANUP ; RETURN
ENDIF
C
C SET DIFFERENT ERROR TOLERANCE IF THE TRAJECTORY Y(S) HAS ANY HIGH
C CURVATURE COMPONENTS.
CURTOL=CURSW*HOLD
RELERR=ARCRE
ABSERR=ARCAE
IF (ANY(ABS(YP-YPOLD) .GT. CURTOL)) THEN
RELERR=ANSRE
ABSERR=ANSAE
ENDIF
C
C TAKE A STEP ALONG THE CURVE.
CALL STEPNF(NC,NFEC,IFLAGC,START,CRASH,HOLD,H,RELERR,ABSERR,
& S,Y,YP,YOLD,YPOLD,A,QR,ALPHA,TZ,PIVOT,W,WP,Z0,Z1,SSPAR)
C PRINT LATEST POINT ON CURVE IF REQUESTED.
IF (TRACE .GT. 0) THEN
WRITE (TRACE,217) ITER,NFEC,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
NFE=NFEC
C CHECK IF THE STEP WAS SUCCESSFUL.
IF (IFLAGC .GT. 0) THEN
ARCLEN=S
IFLAG=IFLAGC
CALL CLEANUP ; RETURN
ENDIF
IF (CRASH) THEN
C RETURN CODE FOR ERROR TOLERANCE TOO SMALL.
IFLAG=2
C CHANGE ERROR TOLERANCES.
IF (ARCRE .LT. RELERR) ARCRE=RELERR
IF (ANSRE .LT. RELERR) ANSRE=RELERR
IF (ARCAE .LT. ABSERR) ARCAE=ABSERR
IF (ANSAE .LT. ABSERR) ANSAE=ABSERR
C CHANGE LIMIT ON NUMBER OF ITERATIONS.
LIMIT=LIMIT-ITER
RETURN
ENDIF
C
IF (Y(1) .GE. 1.0) THEN
C
C USE HERMITE CUBIC INTERPOLATION AND NEWTON ITERATION TO GET THE
C ANSWER AT LAMBDA = 1.0 .
C
C SAVE YOLD FOR ARC LENGTH CALCULATION LATER.
Z0=YOLD
CALL ROOTNF(NC,NFEC,IFLAGC,ANSRE,ANSAE,Y,YP,YOLD,YPOLD,
& A,QR,ALPHA,TZ,PIVOT,W,WP)
C
NFE=NFEC
IFLAG=1
C SET ERROR FLAG IF ROOTNF COULD NOT GET THE POINT ON THE ZERO
C CURVE AT LAMBDA = 1.0 .
IF (IFLAGC .GT. 0) IFLAG=IFLAGC
C CALCULATE FINAL ARC LENGTH.
W=Y-Z0
ARCLEN=S - HOLD + DNRM2(NP1,W,1)
CALL CLEANUP ; RETURN
ENDIF
C
C FOR POLYNOMIAL SYSTEMS AND THE POLSYS1H HOMOTOPY MAP,
C D LAMBDA/DS .GE. 0 NECESSARILY. THIS CONDITION IS FORCED HERE IF
C THE POLY_SWITCH VARIABLE IS PRESENT.
C
IF (POLSYS) THEN
IF (YP(1) .LT. 0.0) THEN
C REVERSE TANGENT DIRECTION SO D LAMBDA/DS = YP(1) > 0 .
YP=-YP
YPOLD=YP
C FORCE STEPNF TO USE THE LINEAR PREDICTOR FOR THE NEXT STEP ONLY.
START=.TRUE.
ENDIF
ENDIF
C
END DO MAIN_LOOP ! ***** END OF MAIN LOOP. *****
C
C LAMBDA HAS NOT REACHED 1 IN 1000 STEPS.
IFLAG=3
ARCLEN=S
RETURN
C
CONTAINS
SUBROUTINE CLEANUP
IF (ALLOCATED(YP)) DEALLOCATE(YP)
IF (ALLOCATED(YOLD)) DEALLOCATE(YOLD)
IF (ALLOCATED(YPOLD)) DEALLOCATE(YPOLD)
END SUBROUTINE CLEANUP
END SUBROUTINE FIXPNF