C coefficients are modified by subroutine SCLGNP . The problem
C is solved with the scaled coefficients and scaled variables.
C The coefficients are returned scaled.
C
C If the projective transformation is specified, then
C essentially the system is reformulated in homogeneous
C coordinates, Z(1), ..., Z(N+1), and solved in complex
C projective space. The resulting solutions are
C untransformed via
C
C X(J) = Z(J)/Z(N+1) J=1, ..., N.
C
C On return,
C
C ROOTS(1,J,M) = real part of X(J) for the Mth path,
C
C ROOTS(2,J,M) = imaginary part of X(J) for the Mth path,
C
C for J=1, ..., N, and
C
C ROOTS(1,N+1,M) = real part of Z(N+1) for the Mth path,
C
C ROOTS(2,N+1,M) = imaginary part of Z(N+1) for the Mth path.
C
C If ROOTS(*,N+1,M) is small, then the associated solution
C should be regarded as being "near infinity". Note that,
C when the projective transformation has been specified, the
C ROOTS values have been untransformed -- that is, divided
C through by Z(N+1) -- unless such division would have caused
C overflow. In this latter case, the affected components of
C ROOTS are set to the largest floating point number (machine
C infinity).
C
C The code can be modified easily to solve systems with complex
C coefficients, COEF . Only the subroutines INITP and FFUNP
C need be changed.
C
C The FORTRAN COMPLEX declaration is not used in POLSYS1H.
C Complex variables are represented by real arrays with first
C index dimensioned 2 and complex operations are evoked by
C subroutine calls.
C
C The total number of paths that will be tracked (if
C IFLG2(M)=-2 for all M) is equal to the "total degree" of the
C system, TOTDG. TOTDG is equal to the products of the
C degrees of all the equations in the system. The degree of
C an equation is the maximum of the degrees of its terms. The
C degree of a term is the sum of the degrees of the variables.
C Thus, TOTDG = IDEG(1) * ... * IDEG(N) where IDEG(J) =
C MAX {JDEG(J,K) | K=1,...,NUMT(J)} where JDEG(J,K) = KDEG(J,1,K) +
C ... + KDEG(J,N,K).
C
C IFLG1 determines whether the system is to be automatically
C scaled by POLSYS1H and whether the projective transformation
C of the system is to be automatically evoked by POLSYS1H. See
c "ON INPUT" below.
C
C IFLG2, EPSBIG, EPSSML, and SSPAR tell the path tracker
C FIXPNF which paths to track and set parameters for the path
C tracker.
C
C NUMRR tells POLSYS1H how many multiples of 1000 steps to try
C before abandoning a path.
C
C The output consists of IFLG1, and of LAMBDA, ROOTS, ARCLEN, and
C NFE for each path. IFLG1 returns input data error information.
C ROOTS gives the solutions themselves, while LAMBDA, ARCLEN,
C and NFE give information about the associated paths.
C
C
C The following subroutines are used directly or indirectly by
C POLSYS1H:
C Special for POLSYS1H:
C INITP , STRPTP , OTPUTP , RHO , RHOJAC ,
C HFUNP , HFUN1P , GFUNP , FFUNP ,
C MULP , POWP , DIVP , SCLGNP .
C From the general HOMPACK routines:
C FIXPNF , ROOT , ROOTNF , STEPNF , TANGNF .
C From LAPACK routines:
C DGEQPF , DGEQRF , DORMQR .
C From BLAS routines:
C DCOPY , DDOT , DGEMM , DGEMV , DGER ,
C DNRM2 , DSCAL , DSWAP , DTRMM , DTRMV , DTRSV ,
C IDAMAX , LSAME , XERBLA .
C
C ON INPUT:
C
C N is the number of equations and variables.
C
C NUMT(1:N) is an integer array. NUMT(J) is the number of terms
C in the Jth equation for J=1 to N.
C
C COEF(1:N,1:) is a real array. COEF(J,K) is
C the Kth coefficient of the Jth equation for J=1 to N,
C K=1 to NUMT(J). The second dimension must be greater than or equal
C to the maximum number of terms in each equation. In other words,
C SIZE(COEF,DIM=2) .GE. MAXT = MAX {NUMT(J) | J=1, ..., N} .
C
C KDEG(1:N,1:N+1,1:) is an integer array.
C KDEG(J,L,K) is the degree of the Lth variable in the Kth
C term of the Jth equation for J=1 to N, L=1 to N, K=1 to NUMT(J).
C SIZE(KDEG,DIM=3) .GE. MAXT = MAX {NUMT(J) | J=1, ..., N} .
C
C IFLG1 =
C 00 if the problem is to be solved without
C calling POLSYS1H' scaling routine, SCLGNP, and
C without using the projective transformtion.
C
C 01 if scaling but no projective transformation is to be used.
C
C 10 if no scaling but projective transformation is to be used.
C
C 11 if both scaling and projective transformation are to be used.
C
C IFLG2(1:TOTDG) is an integer array. If IFLG2(M) = -2, then the
C Mth path is tracked. Otherwise the Mth path is skipped.
C Thus, to find all solutions set IFLG2(M) = -2 for M=1,...,TOTDG.
C Selected paths can be rerun by setting IFLG2(M)=-2 for
C the paths to be rerun and IFLG2(M).NE.-2 for the others.
C
C EPSBIG is the local error tolerance allowed the path tracker along
C the path. ARCRE and ARCAE (in FIXPNF ) are set to EPSBIG.
C
C EPSSML is the accuracy desired for the final solution. ANSRE and
C ANSAE (in FIXPNF ) are set to EPSSML.
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 in FIXPNF and in STEPNF for more information
C about these constants.
C
C NUMRR is the number of multiples of 1000 steps that will be tried
C before abandoning a path.
C
C
C ON OUTPUT:
C
C N, NUMT, COEF, KDEG, EPSBIG, EPSSML, and NUMRR are unchanged.
C
C IFLG1=
C -1 if NUMT is incorrectly dimensioned or invalid.
C -2 if COEF is incorrectly dimensioned.
C -3 if KDEG is incorrectly dimensioned or invalid.
C -4 if any of IFLG2, LAMBDA, ROOTS, ARCLEN, or NFE are
C incorrectly dimensioned.
C -5 if the global work arrays IPAR and PAR could not be
C allocated.
C -6 if IFLG1 on input is not 00 or 01 or 10 or 11.
C Unchanged otherwise.
C
C IFLG2(1:TOTDG) gives information about how the Mth path terminated:
C IFLG2(M) =
C 1 Normal return.
C
C 2 Specified error tolerance cannot be met. Increase EPSBIG
C and EPSSML and rerun.
C
C 3 Maximum number of steps exceeded. To track the path further,
C increase NUMRR and rerun the path. However, the path may
C be diverging, if the LAMBDA value is near 1 and the ROOTS
C values are large.
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 EPSBIG and EPSSML were too lenient. The problem should be
C restarted with smaller error tolerances.
C
C 6 The normal flow Newton iteration in STEPNF or ROOTNF
C failed to converge. The error tolerances EPSBIG or EPSSML
C may be too stringent.
C
C 7 Illegal input parameters, a fatal error.
C
C LAMBDA(M) is the final LAMBDA value for the Mth path, M = 1, ...,
C TOTDG, where LAMBDA is the continuation parameter.
C
C ROOTS(1,J,M), ROOTS(2,J,M) are the real and imaginary parts
C of the Jth variable respectively, for J = 1,...,N, for
C the Mth path, for M = 1,...,TOTDG. If IFLG1 = 10 or 11, then
C ROOTS(1,N+1,M) and ROOTS(2,N+1,M) are the real and
C imaginary parts respectively of the projective
C coordinate of the solution.
C
C ARCLEN(M) is the arc length of the Mth path for M = 1, ..., TOTDG.
C
C NFE(M) is the number of Jacobian matrix evaluations required to
C track the Mth path for M =1, ..., TOTDG.
C
C ----------------------------------------------------------------------
USE HOMOTOPY
USE REAL_PRECISION
INTERFACE
SUBROUTINE INITP(IFLG1,N,NUMT,KDEG,COEF,
& IDEG,FACV,CL,PDG,QDG,R)
USE HOMOTOPY
USE REAL_PRECISION
INTEGER, INTENT(IN):: IFLG1,N,NUMT(:)
INTEGER, INTENT(IN OUT):: KDEG(:,:,:)
REAL (KIND=R8), INTENT(IN OUT):: COEF(:,:)
INTEGER, INTENT(OUT):: IDEG(N)
REAL (KIND=R8), INTENT(OUT)::
& FACV(N),CL(2,N+1),PDG(2,N),QDG(2*N),R(2,N)
END SUBROUTINE INITP
C
SUBROUTINE STRPTP(N,ICOUNT,IDEG,R,X)
USE REAL_PRECISION
INTEGER:: N,ICOUNT(N),IDEG(N)
REAL (KIND=R8):: R(2,N),X(2*N)
END SUBROUTINE STRPTP
C
! SUBROUTINE FIXPNF(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,A,
! & SSPAR,NFE,ARCLEN,POLY_SWITCH)
! 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
! END SUBROUTINE FIXPNF
C
SUBROUTINE OTPUTP(N,NUMPAT,CL,FACV,CLX,X,XNP1)
USE REAL_PRECISION
INTEGER, INTENT(IN):: N,NUMPAT
REAL (KIND=R8), INTENT(IN):: CL(2,N+1),FACV(N)
REAL (KIND=R8), INTENT(IN OUT):: CLX(2,N),X(2,N),XNP1(2)
END SUBROUTINE OTPUTP
END INTERFACE
C
C TYPE DECLARATIONS FOR INPUT AND OUTPUT
C
INTEGER, INTENT(IN):: N,NUMT(:),NUMRR
REAL (KIND=R8), INTENT(IN OUT):: COEF(:,:),SSPAR(8)
INTEGER, INTENT(IN OUT):: KDEG(:,:,:),IFLG1,IFLG2(:)
REAL (KIND=R8), INTENT(IN):: EPSBIG,EPSSML
REAL (KIND=R8), INTENT(OUT):: LAMBDA(:),ROOTS(:,:,:),ARCLEN(:)
INTEGER, INTENT(OUT):: NFE(:)
C
C TYPE DECLARATIONS FOR LOCAL VARIABLES
C
INTEGER:: I,ICOUNT(N),IDEG(N),IDUMMY,IFLAG,IJ,
& IPROFF(15),J,LIPAR(15),LPAR(25),MAXT,N2,N2P1,
& NNFE,NP1,NUMPAT,PROFF(25),TOTDG,TRACE
REAL (KIND=R8):: AARCLN,ANSAE,ANSRE,ARCAE,ARCRE,CL(2,N+1),
& FACV(N),PDG(2,N),QDG(2*N),R(2,N),XNP1(2),Y(2*N+1)
C
C ----------------------------------------------------------------------
N2=2*N
NP1=N+1
N2P1=N2+1
C
C CHECK THAT DIMENSIONS ARE VALID.
C
IF ((SIZE(NUMT) /= N) .OR. ANY(NUMT .LE. 0)) THEN
IFLG1=-1
RETURN
END IF
MAXT = MAXVAL(NUMT)
IF ((SIZE(COEF,DIM=1) /= N) .OR. (SIZE(COEF,DIM=2) < MAXT)) THEN
IFLG1=-2
RETURN
END IF
KDEG = ABS(KDEG)
IF ((SIZE(KDEG,DIM=1) /= N) .OR. (SIZE(KDEG,DIM=2) /= NP1) .OR.
& (SIZE(KDEG,DIM=3) < MAXT) ) THEN
IFLG1=-3
RETURN
END IF
DO J=1,N
IDEG(J)=MAXVAL(SUM(KDEG(J,1:N,1:NUMT(J)),DIM=1))
END DO
TOTDG = PRODUCT(IDEG)
IF ((SIZE(IFLG2) < TOTDG) .OR. (SIZE(LAMBDA) < TOTDG) .OR.
& (SIZE(ROOTS,DIM=3) < TOTDG) .OR. (SIZE(ARCLEN) < TOTDG) .OR.
& (SIZE(NFE) < TOTDG) .OR.
& (IFLG1 <= 1 .AND. SIZE(ROOTS,DIM=2) /= N) .OR.
& (IFLG1 >= 10 .AND. SIZE(ROOTS,DIM=2) /= NP1)) THEN
IFLG1=-4
RETURN
END IF
IF (IFLG1 /= 0 .AND. IFLG1 /= 1 .AND.
& IFLG1 /= 10 .AND. IFLG1 /= 11) THEN
IFLG1=-6
RETURN
END IF
C
C ALLOCATE THE GLOBAL WORK ARRAYS IPAR AND PAR, USED TO COMMUNICATE
C DATA BETWEEN SUBROUTINES VIA THE MODULE HOMOTOPY.
C
ALLOCATE(IPAR(42 + 2*N + N*(N+1)*MAXT),
& PAR(2 + 28*N + 6*N**2 + 7*N*MAXT + 4*N**2*MAXT),STAT=IJ)
IF (IJ .NE. 0) THEN
IFLG1=-5
RETURN
END IF
C
C INITIALIZATION
C
CALL INITP(IFLG1,N,NUMT,KDEG,COEF,
& IDEG,FACV,CL,PDG,QDG,R)
C
C INTEGER VARIABLES AND ARRAYS TO BE PASSED IN IPAR:
C
C IPAR INDEX VARIABLE NAME LENGTH
C ---------- ------------- -----------------
C 1 N 1
C 2 MAXT 1
C 3 PROFF 25
C 4 IPROFF 15
C 5 IDEG N
C 6 NUMT N
C 7 KDEG N*(N+1)*MAXT
C
C
C DOUBLE PRECISION VARIABLES AND ARRAYS TO BE PASSED IN PAR:
C
C PAR INDEX VARIABLE NAME LENGTH
C ---------- ------------- -----------------
C 1 PDG 2*N
C 2 CL 2*(N+1)
C 3 COEF N*MAXT
C 4 H N2
C 5 DHX N2*N2
C 6 DHT N2
C 7 XDGM1 2*N
C 8 XDG 2*N
C 9 G 2*N
C 10 DG 2*N
C 11 PXDGM1 2*N
C 12 PXDG 2*N
C 13 F 2*N
C 14 DF 2*N*(N+1)
C 15 XX 2*N*(N+1)*MAXT
C 16 TRM 2*N*MAXT
C 17 DTRM 2*N*(N+1)*MAXT
C 18 CLX 2*N
C 19 DXNP1 2*N
C
C SET LENGTHS OF VARIABLES
LIPAR(1)=1
LIPAR(2)=1
LIPAR(3)=25
LIPAR(4)=15
LIPAR(5)=N
LIPAR(6)=N
LIPAR(7)=N*(N+1)*MAXT
LPAR( 1)=2*N
LPAR( 2)=2*NP1
LPAR( 3)=N*MAXT
LPAR( 4)=N2
LPAR( 5)=N2*N2
LPAR( 6)=N2
LPAR( 7)=2*N
LPAR( 8)=2*N
LPAR( 9)=2*N
LPAR(10)=2*N
LPAR(11)=2*N
LPAR(12)=2*N
LPAR(13)=2*N
LPAR(14)=2*N*NP1
LPAR(15)=2*N*NP1*MAXT
LPAR(16)=2*N*MAXT
LPAR(17)=2*N*NP1*MAXT
LPAR(18)=2*N
LPAR(19)=2*N
C
C PROFF AND IPROFF ARE OFFSETS THAT DEFINE THE VARIABLES LISTED ABOVE
PROFF(1)=1
DO I=2,19
PROFF(I)=PROFF(I-1)+LPAR(I-1)
END DO
IPROFF(1)=1
DO I=2,7
IPROFF(I)=IPROFF(I-1)+LIPAR(I-1)
END DO
C
C DEFINE VARIABLES
IPAR(1)=N
IPAR(2)=MAXT
IPAR(IPROFF(3):IPROFF(3)+18) = PROFF(1:19)
IPAR(IPROFF(4):IPROFF(4)+ 6) = IPROFF(1:7)
IPAR(IPROFF(5):IPROFF(5)+N-1) = IDEG(1:N)
IPAR(IPROFF(6):IPROFF(6)+N-1) = NUMT(1:N)
IPAR(IPROFF(7):IPROFF(7)+LIPAR(7)-1) =
& PACK(KDEG(:,:,1:MAXT),.TRUE.)
PAR(PROFF(1):PROFF(1)+LPAR(1)-1) = PACK(PDG,.TRUE.)
PAR(PROFF(2):PROFF(2)+LPAR(2)-1) = PACK(CL,.TRUE.)
PAR(PROFF(3):PROFF(3)+LPAR(3)-1) = PACK(COEF(:,1:MAXT),.TRUE.)
C
C ICOUNT IS A COUNTER USED BY "STRPTP"
ICOUNT(1)=0
ICOUNT(2:N)=1
C
C PATHS LOOP -- ITERATE THROUGH PATHS
C
PATHS: DO NUMPAT = 1,TOTDG
C GET A START POINT, Y, FOR THE PATH.
Y(1) = 0.0
CALL STRPTP(N,ICOUNT,IDEG,R,Y(2:N2P1))
C CHECK WHETHER PATH IS TO BE FOLLOWED.
IFLAG = IFLG2(NUMPAT)
IF (IFLAG .NE. -2) CYCLE PATHS
ARCRE = EPSBIG
ARCAE = ARCRE
ANSRE = EPSSML
ANSAE = ANSRE
TRACE = 0
C TRACK A HOMOTOPY PATH.
DO IDUMMY=1,MAX(NUMRR,1)
CALL FIXPNF(N2,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,
& QDG,SSPAR,NNFE,AARCLN, POLY_SWITCH=.TRUE.)
IF (IFLAG .NE. 2 .AND. IFLAG .NE. 3) EXIT
END DO
C UNSCALE AND UNTRANSFORM COMPUTED SOLUTION.
CALL OTPUTP(N,NUMPAT,CL,FACV,
& PAR(PROFF(18):PROFF(18)+LPAR(18)-1),Y(2:N2P1),XNP1)
LAMBDA(NUMPAT) = Y(1)
ROOTS(1,1:N,NUMPAT) = Y(2:N2P1:2)
ROOTS(2,1:N,NUMPAT) = Y(3:N2P1:2)
ROOTS(1:2,NP1,NUMPAT) = XNP1
C
ARCLEN(NUMPAT)= AARCLN
NFE(NUMPAT) = NNFE
IFLG2(NUMPAT) = IFLAG
END DO PATHS
C CLEAN UP WORK SPACE.
IF (ALLOCATED(IPAR)) DEALLOCATE(IPAR)
IF (ALLOCATED(PAR)) DEALLOCATE(PAR)
RETURN
END SUBROUTINE POLSYS1H
END MODULE HOMPACK90
SUBROUTINE DIVP(XXXX,YYYY,ZZZZ,IERR)
C
C THIS SUBROUTINE PERFORMS DIVISION OF COMPLEX NUMBERS:
C ZZZZ = XXXX/YYYY
C
C ON INPUT:
C
C XXXX IS AN ARRAY OF LENGTH TWO REPRESENTING THE FIRST COMPLEX
C NUMBER, WHERE XXXX(1) = REAL PART OF XXXX AND XXXX(2) =
C IMAGINARY PART OF XXXX.
C
C YYYY IS AN ARRAY OF LENGTH TWO REPRESENTING THE SECOND COMPLEX
C NUMBER, WHERE YYYY(1) = REAL PART OF YYYY AND YYYY(2) =
C IMAGINARY PART OF YYYY.
C
C ON OUTPUT:
C
C ZZZZ IS AN ARRAY OF LENGTH TWO REPRESENTING THE RESULT OF
C THE DIVISION, ZZZZ = XXXX/YYYY, WHERE ZZZZ(1) =
C REAL PART OF ZZZZ AND ZZZZ(2) = IMAGINARY PART OF ZZZZ.
C
C IERR =
C 1 IF DIVISION WOULD HAVE CAUSED OVERFLOW. IN THIS CASE, THE
C APPROPRIATE PARTS OF ZZZZ ARE SET EQUAL TO THE LARGEST
C FLOATING POINT NUMBER, AS GIVEN BY FUNCTION HUGE .
C
C 0 IF DIVISION DOES NOT CAUSE OVERFLOW.
C
C DECLARATION OF INPUT
USE REAL_PRECISION
REAL (KIND=R8), DIMENSION(2), INTENT(IN):: XXXX,YYYY
C
C DECLARATION OF OUTPUT
INTEGER, INTENT(OUT):: IERR
REAL (KIND=R8), DIMENSION(2), INTENT(OUT):: ZZZZ
C
C DECLARATION OF VARIABLES
REAL (KIND=R8):: DENOM,XNUM
C
IERR = 0
DENOM = YYYY(1)*YYYY(1) + YYYY(2)*YYYY(2)
XNUM = XXXX(1)*YYYY(1) + XXXX(2)*YYYY(2)
IF (ABS(DENOM) .GE. 1.0 .OR. ( ABS(DENOM) .LT. 1.0 .AND.
& ABS(XNUM)/HUGE(1.0_R8) .LT. ABS(DENOM) ) ) THEN
ZZZZ(1) = XNUM/DENOM
ELSE
ZZZZ(1) = HUGE(1.0_R8)
IERR =1
END IF
XNUM = XXXX(2)*YYYY(1) - XXXX(1)*YYYY(2)
IF (ABS(DENOM) .GE. 1.0 .OR. ( ABS(DENOM) .LT. 1.0 .AND.
& ABS(XNUM)/HUGE(1.0_R8) .LT. ABS(DENOM) ) ) THEN
ZZZZ(2) = XNUM/DENOM
ELSE
ZZZZ(2) = HUGE(1.0_R8)
IERR =1
END IF
RETURN
END SUBROUTINE DIVP
SUBROUTINE FFUNP(N,NUMT,MAXT,KDEG,COEF,CL,X,
& XX,TRM,DTRM,CLX,DXNP1,F,DF)
C
C FFUNP EVALUATES THE SYSTEM "F(X)=0" AND ITS PARTIAL
C DERIVATIVES, USING THE "TABLEAU" INPUT: N,NUMT,KDEG,COEF.
C
C FFUNP CAN BE MADE MORE EFFICIENT BY CUSTOMIZING IT TO
C PARTICULAR SYSTEM TYPES. FOR EXAMPLE,
C IF X(1)**2 AND X(1)**3 ARE USED IN SEVERAL
C EQUATIONS, THE CURRENT FFUNP RECOMPUTES BOTH OF THESE FOR
C EACH EQUATION. BUT (OF COURSE) WE CAN COMPUTE
C X1SQ=X(1)**2 AND X1CU=XSQ(1)*X(1), AND
C USE THESE IN EACH OF THE EQUATIONS.
C
C THE PART OF THE CODE BELOW LABELED "BLOCK A" CAN BE
C CUSTOMIZED IN THIS WAY. (THE CODE OUTSIDE OF
C BLOCK A CONCERNS THE PROJECTIVE TRANSFORMATION AND NEED NOT
C BE CHANGED.) HOWEVER, BLOCK A REQUIRES THE HOMOGENEOUS FORM
C OF THE POLYNOMIALS RATHER THAN THE STANDARD FORM. FURTHER,
C THE PARTIAL DERIVATIVES WITH RESPECT TO ALL N+1 PROJECTIVE
C VARIABLES MUST BE COMPUTED. MORE EXPLICITLY,
C THE ORIGINAL SYSTEM, F(X)=0, IS GIVEN IN "NON-HOMOGENEOUS FORM" AS
C DESCRIBED IN SUBROUTINE POLSYS. F(X) IS
C REPRESENTED IN "HOMOGENEOUS FORM" AS FOLLOWS:
C
C NUMT(J)
C
C F(J) = SUM TRM(J,K)
C
C K=1
C
C WHERE TRM(J,K)=COEF(J,K) * XX(J,1,K)*XX(J,2,K)* ... *XX(J,N+1,K)
C
C WITH XX(J,L,K) = X(L)**KDEG(J,L,K) FOR J=1 TO N, L=1 TO N, AND
C K=1 TO NUMT(J), AND WITH XX(J,N+1,K) = XNP1**KDEG(J,N+1,K) FOR J=1 TO
C N AND K=1 TO NUMT(J), WHERE XNP1 IS THE "HOMOGENEOUS COORDINATE,"
C KDEG(J,N+1,K)=IDEG(J)-(KDEG(J,1,K)+ ... + KDEG(J,N,K)),
C AND IDEG(J) THE DEGREE OF THE J-TH EQUATION. XNP1 IS GENERATED
C FROM X AND CL BEFORE BLOCK A.
C
C IN THIS DISCUSSION WE HAVE OMITTED, FOR SIMPLICITY OF
C EXPOSITION, THE LEADING INDEX, WHICH DIFFERENTIATES THE
C REAL AND IMAGINARY PARTS. HOWEVER, THIS INDEX MUST NOT BE
C OMITTED IN THE CODE.
C
C WE COMPLETE THE EXPOSITION OF "REPLACING BLOCK A WITH MORE EFFICIENT
C CODE" WITH AN EXPLICIT EXAMPLE. FIRST, THE SYSTEM IS DESCRIBED.
C THEN THE CODE THAT SHOULD BE USED IS GIVEN (COMMENTED OUT).
C IN TESTS POLSYS WITH THE MORE EFFICIENT FFUNP RAN ABOUT TWICE AS
C FAST AS WITH THE GENERIC FFUNP .
C
C HERE IS THE SYSTEM TO BE SOLVED:
C
C F(1) = COEF(1,1) * X(1)**4
C & + COEF(1,2) * X(1)**3 * X(2)
C & + COEF(1,3) * X(1)**3
C & + COEF(1,4) * X(1)
C & + COEF(1,5)
C F(2) = COEF(2,1) * X(1) * X(2)**2
C & + COEF(2,2) X(2)**2
C & + COEF(2,3)
C
C THE REPLACEMENT CODE REQUIRES THE FOLLOWING DECLARATIONS:
C DOUBLE PRECISION X1SQ,X1CU,X2SQ,X3SQ,X3CU,
C & TEMPA,TEMPB,TEMPC,TEMPD,TEMPE,TEMPF
C DIMENSION X1SQ(2),X1CU(2),X2SQ(2),X3SQ(2),X3CU(2),
C & TEMPA(2),TEMPB(2),TEMPC(2),TEMPD(2),TEMPE(2),TEMPF(2)
C
C HERE IS CODE TO REPLACE BLOCK A:
C
C****************** BEGIN BLOCK A *******************
C
C CALL MULP(X(1,1),X(1,1),X1SQ)
C CALL MULP(X1SQ ,X(1,1),X1CU)
C CALL MULP(X(1,2),X(1,2),X2SQ)
C CALL MULP(XNP1, XNP1, X3SQ)
C CALL MULP(X3SQ ,XNP1, X3CU)
C
C DO 1 I=1,2
C TEMPA(I)= COEF(1,1) * X(I,1)
C & + COEF(1,2) * X(I,2)
C & + COEF(1,3) * XNP1(I)
C TEMPB(I)= COEF(1,4) * X(I,1)
C & + COEF(1,5) * XNP1(I)
C 1 CONTINUE
C
C CALL MULP(X1SQ, TEMPA,TEMPC)
C CALL MULP(X(1,1),TEMPC,TEMPD)
C CALL MULP(X3SQ, TEMPB,TEMPE)
C CALL MULP(XNP1, TEMPE,TEMPF)
C
C DO 2 I=1,2
C F(I,1)=TEMPD(I) + TEMPF(I)
C DF(I,1,1)= 3. *TEMPC(I) + COEF(1,1)*X1CU(I) + COEF(1,4)*X3CU(I)
C DF(I,1,2)= COEF(1,2) * X1CU(I)
C DF(I,1,3)= COEF(1,3)*X1CU(I) + 3. *TEMPE(I) + COEF(1,5)*X3CU(I)
C
C TEMPA(I) = COEF(2,1) * X(I,1) + COEF(2,2) * XNP1(I)
C 2 CONTINUE
C
C CALL MULP(TEMPA,X(1,2),TEMPB)
C CALL MULP(TEMPB,X(1,2),TEMPC)
C
C DO 3 I=1,2
C F(I,2) = TEMPC(I) + COEF(2,3) * X3CU(I)
C DF(I,2,1) = COEF(2,1) * X2SQ(I)
C DF(I,2,2) = 2. * TEMPB(I)
C DF(I,2,3) = COEF(2,2) * X2SQ(I) + COEF(2,3) * 3. * X3SQ(I)
C 3 CONTINUE
C****************** END OF BLOCK A *******************
C
C ON INPUT:
C
C N IS THE NUMBER OF EQUATIONS AND VARIABLES.
C
C NUMT(J) IS THE NUMBER OF TERMS IN THE JTH EQUATION.
C
C MAXT IS AN UPPER BOUND ON NUMT(J) FOR J=1 TO N.
C
C KDEG(J,L,K) IS THE DEGREE OF THE L-TH VARIABLE IN THE K-TH TERM
C OF THE J-TH EQUATION.
C
C COEF(J,K) IS THE K-TH COEFFICIENT OF THE J-TH EQUATION.
C
C CL IS USED TO DEFINE THE PROJECTIVE TRANSFORMATION. IF
C THE PROJECTIVE TRANSFORMATION IS NOT SPECIFIED, THEN CL
C CONTAINS DUMMY VALUES.
C
C X(1,J), X(2,J) ARE THE REAL AND IMAGINARY PARTS RESPECTIVELY OF
C THE J-TH INDEPENDENT VARIABLE.
C
C XX, TRM, DTRM, CLX, DXNP1 ARE WORKSPACE VARIABLES.
C
C ON OUTPUT:
C
C F(1,J), F(2,J) ARE THE REAL AND IMAGINARY PARTS RESPECTIVELY OF
C THE J-TH EQUATION.
C
C DF(1,J,K), DF(2,J,K) ARE THE REAL AND IMAGINARY PARTS RESPECTIVELY
C OF THE K-TH PARTIAL DERIVATIVE OF THE J-TH EQUATION.
C
C
C VARIABLES: XNP1,TEMP1,TEMP2.
C
C NOTE: XNP1(1), XNP1(2) ARE THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE PROJECTIVE VARIABLE. XNP1 IS UNITY
C IF THE PROJECTIVE TRANSFORMATION IS NOT SPECIFIED.
C
C SUBROUTINES: MULP,POWP,DIVP.
C
USE REAL_PRECISION
C
INTERFACE
SUBROUTINE DIVP(XXXX,YYYY,ZZZZ,IERR)
USE REAL_PRECISION
REAL (KIND=R8), DIMENSION(2), INTENT(IN):: XXXX,YYYY
INTEGER, INTENT(OUT):: IERR
REAL (KIND=R8), DIMENSION(2), INTENT(OUT):: ZZZZ
END SUBROUTINE DIVP
SUBROUTINE MULP(XXXX,YYYY,ZZZZ)
USE REAL_PRECISION
REAL (KIND=R8), DIMENSION(2), INTENT(IN):: XXXX,YYYY
REAL (KIND=R8), DIMENSION(2), INTENT(OUT):: ZZZZ
END SUBROUTINE MULP
SUBROUTINE POWP(NNNN,XXXX,YYYY)
USE REAL_PRECISION
INTEGER, INTENT(IN):: NNNN
REAL (KIND=R8), DIMENSION(2), INTENT(IN):: XXXX
REAL (KIND=R8), DIMENSION(2), INTENT(IN OUT):: YYYY
END SUBROUTINE POWP
END INTERFACE
C
C DECLARATION OF INPUT AND OUTPUT:
INTEGER, INTENT(IN):: N,NUMT(N),MAXT,KDEG(N,N+1,MAXT)
REAL (KIND=R8), INTENT(IN):: COEF(N,MAXT),CL(2,N+1),X(2,N)
REAL (KIND=R8), INTENT(IN OUT)::
& XX(2,N,N+1,MAXT),TRM(2,N,MAXT),DTRM(2,N,N+1,MAXT)
REAL (KIND=R8), INTENT(OUT)::
& CLX(2,N),DXNP1(2,N),F(2,N),DF(2,N,N+1)
C
C DECLARATION OF LOCAL VARIABLES:
INTEGER:: IERR,J,K,L,M,NNNN,NP1
REAL (KIND=R8), DIMENSION(2):: TEMP1,TEMP2,XNP1
C
NP1=N+1
C
C GENERATE XNP1, THE PROJECTIVE COORDINATE, AND ITS DERIVATIVES.
DO J=1,N
CALL MULP(CL(1,J),X(1,J),CLX(1,J))
END DO
C
XNP1(1:2)=CL(1:2,NP1) + SUM(CLX(1:2,1:N),DIM=2)
DXNP1(1:2,1:N)=CL(1:2,1:N)
C
C****************** BEGIN BLOCK A *******************
C
C "BLOCK A" TAKES X AND XNP1 AS INPUT AND RETURNS F
C AND DF AS OUTPUT. F IS THE HOMOGENEOUS FORM OF THE
C ORIGINAL F, AND DF CONSISTS OF THE PARTIAL
C DERIVATIVES OF THE HOMOGENEOUS FORM OF F WITH RESPECT
C TO THE N+1 VARIABLES X(1), ... ,X(N), XNP1.
C
C BEGIN "COMPUTE F"
C
DO J=1,N
DO K=1,NUMT(J)
CALL POWP(KDEG(J,NP1,K),XNP1, XX(1,J,NP1,K))
DO L=1,N
CALL POWP(KDEG(J, L,K),X(1,L),XX(1,J, L,K))
END DO
END DO
END DO
TRM = 0.0
DO J=1,N
DO K=1,NUMT(J)
TRM(1,J,K)=COEF(J,K)
DO L=1,NP1
CALL MULP(XX(1,J,L,K), TRM(1,J,K),TEMP1)
TRM(1:2,J,K ) = TEMP1(1:2)
END DO
END DO
END DO
F(1:2,1:N) = SUM(TRM(1:2,1:N,:),DIM=3)
C
C END OF "COMPUTE F"
C
C BEGIN "COMPUTE DF"
C
J_LOOP: DO J=1,N
K_LOOP: DO K=1,NUMT(J)
M_LOOP: DO M=1,NP1
C
C IF TERM DOES NOT INCLUDE X(M), SET PARTIAL DERIVATIVE OF TERM
C EQUAL TO ZERO.
IF(KDEG(J,M,K) .EQ. 0) THEN
DTRM(1:2,J,M,K) = 0.0
ELSE
C
C IF TERM DOES INCLUDE X(M), TRY COMPUTING THE PARTIAL BY DIVIDING
C THE TERM BY X(M).
IF(M.LE.N) CALL DIVP(TRM(1,J,K),X(1,M),DTRM(1,J,M,K),IERR)
IF(M.EQ.NP1) CALL DIVP(TRM(1,J,K),XNP1,DTRM(1,J,M,K),IERR)
IF (IERR .EQ. 0) THEN
DTRM(1:2,J,M,K)=KDEG(J,M,K)*DTRM(1:2,J,M,K)
ELSE
C
C IF DIVISION WOULD CAUSE OVERFLOW, GENERATE THE PARTIAL BY
C THE POLYNOMIAL FORMULA.
DTRM(1,J,M,K)=COEF(J,K)
DTRM(2,J,M,K)=0.0
DO L=1,NP1
IF (L .EQ. M) CYCLE
CALL MULP(XX(1,J,L,K),DTRM(1,J,M,K),TEMP1)
DTRM(1:2,J,M,K)=TEMP1(1:2)
END DO
NNNN=KDEG(J,M,K)-1
IF (M .LE. N) CALL POWP(NNNN,X(1,M),TEMP2)
IF (M .EQ. NP1) CALL POWP(NNNN,XNP1 ,TEMP2)
CALL MULP(TEMP2,TEMP1,DTRM(1,J,M,K))
DTRM(1:2,J,M,K)=KDEG(J,M,K)*DTRM(1:2,J,M,K)
END IF
END IF
END DO M_LOOP
END DO K_LOOP
END DO J_LOOP
DO J=1,N
DF(1:2,J,1:NP1) = SUM(DTRM(1:2,J,1:NP1,1:NUMT(J)), DIM=3)
END DO
C
C END OF "COMPUTE DF"
C******************* END BLOCK A ********************
C
C CONVERT DF TO BE PARTIALS WITH RESPECT TO X(1), ... ,X(N),
C BY APPLYING THE CHAIN RULE WITH XNP1 CONSIDERED A FUNCTION OF
C OF X(1), ... ,X(N).
C
DO J=1,N
DO K=1,N
CALL MULP(DF(1,J,NP1),DXNP1(1,K),TEMP1)
DF(1:2,J,K)=DF(1:2,J,K)+TEMP1(1:2)
END DO
END DO
RETURN
END SUBROUTINE FFUNP
SUBROUTINE FODE(S,Y,YP,YPOLD,A,QR,ALPHA,TZ,PIVOT,NFE,N,IFLAG)
C
C SUBROUTINE FODE IS USED BY SUBROUTINE STEPS TO SPECIFY THE
C ORDINARY DIFFERENTIAL EQUATION DY/DS = G(S,Y) , WHOSE SOLUTION
C IS THE ZERO CURVE OF THE HOMOTOPY MAP. S = ARC LENGTH,
C YP = DY/DS, AND Y(S) = (LAMBDA(S), X(S)) .
C
C CALLS DGEQPF , DNRM2 .
C
USE HOMOTOPY
USE REAL_PRECISION
REAL (KIND=R8):: DNRM2,S,YPNORM
INTEGER:: I,IFLAG,IK,K,KP1,LW,N,NFE,NP1
C
C ***** ARRAY DECLARATIONS. *****
C
REAL (KIND=R8):: A(:),Y(:),YP(N+1),YPOLD(N+1)
C
C ARRAYS FOR COMPUTING THE JACOBIAN MATRIX AND ITS KERNEL.
REAL (KIND=R8):: ALPHA(3*N+3),QR(N,N+1),TZ(N+1)
INTEGER, DIMENSION(N+1):: PIVOT
C
C ***** END OF DIMENSIONAL INFORMATION. *****
C
C
NP1=N+1
NFE=NFE+1
C NFE CONTAINS THE NUMBER OF JACOBIAN EVALUATIONS.
C * * * * * * * * * * * * * * * * *
C
C COMPUTE THE JACOBIAN MATRIX, STORE IT IN QR.
C
IF (IFLAG .EQ. -2) THEN
C
C QR = ( D RHO(A,LAMBDA,X)/D LAMBDA , D RHO(A,LAMBDA,X)/DX ) .
C
DO K=1,NP1
CALL RHOJAC(A,Y(1),Y(2:NP1),QR(:,K),K)
END DO
ELSE
CALL F(Y(2:NP1),TZ(1:N))
QR(:,1)=A
IF (IFLAG .EQ. 0) THEN
C
C QR = ( A - F(X), I - LAMBDA*DF(X) ) .
C
QR(:,1)=QR(:,1)-TZ(1:N)
DO K=1,N
CALL FJAC(Y(2:NP1),TZ(1:N),K)
KP1=K+1
QR(:,KP1)=-Y(1)*TZ(1:N)
QR(K,KP1)=1.0+QR(K,KP1)
END DO
ELSE
C
C QR = ( F(X) - X + A, LAMBDA*DF(X) + (1 - LAMBDA)*I ) .
C
QR(:,1)=QR(:,1)+TZ(1:N)-Y(2:NP1)
DO K=1,N
CALL FJAC(Y(2:NP1),TZ(1:N),K)
KP1=K+1
QR(:,KP1)=Y(1)*TZ(1:N)
QR(K,KP1)=1.0-Y(1)+QR(K,KP1)
END DO
ENDIF
ENDIF
C
C * * * * * * * * * * * * * * * * *
C REDUCE THE JACOBIAN MATRIX TO UPPER TRIANGULAR FORM.
C
PIVOT = 0
C
CALL DGEQPF(N,NP1,QR,N,PIVOT,YP,ALPHA,K)
C
IF (ABS(QR(N,N)) .LE. ABS(QR(1,1))*EPSILON(1.0_R8)) THEN
IFLAG=4
RETURN
ENDIF
C COMPUTE KERNEL OF JACOBIAN, WHICH SPECIFIES YP=DY/DS.
TZ(NP1)=1.0
DO LW=1,N
I=NP1-LW
IK=I+1
TZ(I)=-DOT_PRODUCT(QR(I,IK:NP1),TZ(IK:NP1))/QR(I,I)
END DO
YPNORM=DNRM2(NP1,TZ,1)
YP(PIVOT)=TZ/YPNORM
IF (DOT_PRODUCT(YP,YPOLD) .LT. 0.0) YP=-YP
C
C SAVE CURRENT DERIVATIVE (= TANGENT VECTOR) IN YPOLD .
YPOLD=YP
RETURN
END SUBROUTINE FODE
SUBROUTINE FODEDS(S,Y,YP,N,IFLAG,YPOLD,A,NDIMA,LENQR,MODE,NFE)
C
C SUBROUTINE FODEDS IS USED BY SUBROUTINE STEPDS TO SPECIFY THE
C ORDINARY DIFFERENTIAL EQUATION DY/DS = G(S,Y) , WHOSE SOLUTION
C IS THE ZERO CURVE OF THE HOMOTOPY MAP. S = ARC LENGTH,
C YP = DY/DS, AND Y(S) = (X(S), LAMBDA(S)) .
C
USE HOMOTOPY, QR => QRSPARSE
USE REAL_PRECISION
REAL (KIND=R8):: LAMBDA,S,YPNORM
INTEGER:: IFLAG,J,JPOS,LENQR,MODE,N,NDIMA,NFE,NP1
REAL (KIND=R8):: A(NDIMA),Y(N+1),YP(N+1),YPOLD(N+1)
INTERFACE
SUBROUTINE PCGDS(N,LENQR,IFLAG,YP,RHO)
USE REAL_PRECISION
INTEGER, INTENT(IN):: LENQR,N
INTEGER, INTENT(IN OUT):: IFLAG
REAL (KIND=R8), INTENT(IN OUT):: YP(N+1)
REAL (KIND=R8), OPTIONAL, INTENT(IN):: RHO(N)
END SUBROUTINE PCGDS
SUBROUTINE GMRILUDS(N,LENQR,IFLAG,YP,RHO)
USE REAL_PRECISION
INTEGER, INTENT(IN):: LENQR,N
INTEGER, INTENT(IN OUT):: IFLAG
REAL (KIND=R8), INTENT(IN OUT):: YP(N+1)
REAL (KIND=R8), OPTIONAL, INTENT(IN):: RHO(N)
END SUBROUTINE GMRILUDS
FUNCTION DNRM2(N,X,STRIDE)
USE REAL_PRECISION
INTEGER:: N,STRIDE
REAL (KIND=R8):: DNRM2,X(N)
END FUNCTION DNRM2
END INTERFACE
C
C ***** END OF SPECIFICATION INFORMATION. *****
C
NP1=N+1
NFE=NFE+1
C NFE CONTAINS THE NUMBER OF JACOBIAN EVALUATIONS.
LAMBDA=Y(NP1)
ROWPOS(NP1)=LENQR+1
C * * * * * * * * * * * * * * * * *
C MODE = 1 STORAGE FORMAT.
C * * * * * * * * * * * * * * * * *
C
IF (MODE .EQ. 1) THEN
C COMPUTE THE JACOBIAN MATRIX, STORE IT IN [QR | -PP] .
C
IF (IFLAG .EQ. -2) THEN
C
C [QR | -PP] = [ D RHO(A,X,LAMBDA)/DX | D RHO(A,X,LAMBDA)/D LAMBDA ] .
C
C PP = - (D RHO(A,X,LAMBDA)/D LAMBDA) .
CALL RHOJS(A,LAMBDA,Y(1:N))
C
ELSE
CALL F(Y(1:N),PP)
IF (IFLAG .EQ. 0) THEN
C
C [QR | -PP] = [ I - LAMBDA*DF(X) | A - F(X) ] .
C
PP = PP - A(1:N)
CALL FJACS(Y(1:N))
QR = (-LAMBDA)*QR
QR(ROWPOS(1:N)) = QR(ROWPOS(1:N)) + 1.0
ELSE
C
C [QR | -PP] = [ LAMBDA*DF(X) + (1 - LAMBDA)*I | F(X) - X + A ] .
C
PP = Y(1:N) - A(1:N) - PP
CALL FJACS(Y(1:N))
QR = LAMBDA*QR
QR(ROWPOS(1:N)) = QR(ROWPOS(1:N)) + 1.0 - LAMBDA
ENDIF
ENDIF
ELSE
C * * * * * * * * * * * * * * * * *
C MODE = 2 STORAGE FORMAT.
C * * * * * * * * * * * * * * * * *
C
IF (IFLAG .EQ. -2) THEN
C
C [QR] = [ D RHO(A,X,LAMBDA)/DX , D RHO(A,X,LAMBDA)/D LAMBDA ] .
C
CALL RHOJS(A,LAMBDA,Y(1:N))
C
ELSE
CALL F(Y(1:N),PP)
IF (IFLAG .EQ. 0) THEN
C
C [QR | -PP] = [ I - LAMBDA*DF(X) | A - F(X) ] .
C
PP = PP - A(1:N)
CALL FJACS(Y(1:N))
QR = (-LAMBDA)*QR
C FIND INDEX JPOS OF DIAGONAL ELEMENT IN JTH ROW OF QR.
DO J=1,N
JPOS=ROWPOS(J)
DO
IF (COLPOS(JPOS) .EQ. J) EXIT
JPOS=JPOS+1
IF (JPOS < ROWPOS(J+1)) CYCLE
IFLAG=4
RETURN
END DO
QR(JPOS) = QR(JPOS) + 1.0
END DO
ELSE
C
C [QR | -PP] = [ LAMBDA*DF(X) + (1 - LAMBDA)*I | F(X) - X + A ] .
C
PP = Y(1:N) - A(1:N) - PP
CALL FJACS(Y(1:N))
QR = LAMBDA*QR
C FIND INDEX JPOS OF DIAGONAL ELEMENT IN JTH ROW OF QR.
DO J=1,N
JPOS=ROWPOS(J)
DO
IF (COLPOS(JPOS) .EQ. J) EXIT
JPOS=JPOS+1
IF (JPOS < ROWPOS(J+1)) CYCLE
IFLAG=4
RETURN
END DO
QR(JPOS) = QR(JPOS) + 1.0 - LAMBDA
END DO
ENDIF
ENDIF
ENDIF
C
C * * * * * * * * * * * * * * * * *
YP=YPOLD
C COMPUTE KERNEL OF JACOBIAN, WHICH SPECIFIES YP=DY/DS, USING A
C PRECONDITIONED CONJUGATE GRADIENT ALGORITHM.
SELECT CASE (MODE)
CASE (1)
CALL PCGDS(N,LENQR,IFLAG,YP)
CASE (2)
CALL GMRILUDS(N,LENQR,IFLAG,YP)
END SELECT
IF (IFLAG .GT. 0) RETURN
C
C NORMALIZE TANGENT VECTOR YP.
YPNORM=DNRM2(NP1,YP,1)
YP = (1.0/YPNORM)*YP
C
C CHOOSE UNIT TANGENT VECTOR DIRECTION TO MAINTAIN CONTINUITY.
IF (DOT_PRODUCT(YP,YPOLD) .LT. 0.0) YP = -YP
C
C SAVE CURRENT DERIVATIVE (= TANGENT VECTOR) IN YPOLD .
YPOLD = YP
C
RETURN
END SUBROUTINE FODEDS
SUBROUTINE GFUNP(N,IDEG,PDG,QDG,X,XDGM1,XDG,PXDGM1,PXDG,G,DG)
C
C GFUNP EVALUATES THE START EQUATION "G".
C
C ON INPUT:
C
C N IS THE NUMBER OF VARIABLES.
C
C IDEG(J) IS THE DEGREE OF THE J-TH EQUATION.
C
C PDG(1,J), PDG(2,J) ARE THE REAL AND IMAGINARY PARTS
C OF THE POWERS OF P USED TO DEFINE G.
C
C QDG(1,J), QDG(2,J) ARE THE REAL AND IMAGINARY PARTS
C OF THE POWERS OF Q USED TO DEFINE G.
C
C X(1,J), X(2,J) ARE THE REAL AND IMAGINARY PARTS OF THE
C J-TH INDEPENDENT VARIABLE.
C
C XDGM1,XDG,PXDGM1,PXDG ARE WORKSPACE ARRAYS.
C
C ON OUTPUT:
C
C N,IDEG,PDG,QDG, AND X ARE UNCHANGED.
C
C G(1,J),G(2,J) ARE THE REAL AND IMAGINARY PARTS OF THE
C J-TH START EQUATION.
C
C DG(1,J),DG(2,J) ARE THE REAL AND IMAGINARY PARTS OF THE
C PARTIAL DERIVATIVES OF THE J-TH START EQUATION WITH RESPECT TO THE
C J-TH INDEPENDENT VARIABLE.
C
C SUBROUTINES: MULP, POWP.
C
USE REAL_PRECISION
C
INTERFACE
SUBROUTINE DIVP(XXXX,YYYY,ZZZZ,IERR)
USE REAL_PRECISION
REAL (KIND=R8), DIMENSION(2), INTENT(IN):: XXXX,YYYY
INTEGER, INTENT(OUT):: IERR
REAL (KIND=R8), DIMENSION(2), INTENT(OUT):: ZZZZ
END SUBROUTINE DIVP
SUBROUTINE MULP(XXXX,YYYY,ZZZZ)
USE REAL_PRECISION
REAL (KIND=R8), DIMENSION(2), INTENT(IN):: XXXX,YYYY
REAL (KIND=R8), DIMENSION(2), INTENT(OUT):: ZZZZ
END SUBROUTINE MULP
SUBROUTINE POWP(NNNN,XXXX,YYYY)
USE REAL_PRECISION
INTEGER, INTENT(IN):: NNNN
REAL (KIND=R8), DIMENSION(2), INTENT(IN):: XXXX
REAL (KIND=R8), DIMENSION(2), INTENT(IN OUT):: YYYY
END SUBROUTINE POWP
END INTERFACE
C
C DECLARATION OF INPUT AND OUTPUT:
INTEGER, INTENT(IN):: N,IDEG(N)
REAL (KIND=R8), INTENT(IN):: PDG(2,N),QDG(2,N),X(2,N)
REAL (KIND=R8), INTENT(IN OUT):: XDGM1(2,N),XDG(2,N),PXDGM1(2,N),
& PXDG(2,N)
REAL (KIND=R8), INTENT(OUT):: G(2,N),DG(2,N)
C
C DECLARATION LOCAL OF VARIABLES
INTEGER:: I,J
C
C COMPUTE THE (IDEG-1)-TH AND IDEG-TH POWER OF X
DO J=1,N
CALL POWP(IDEG(J)-1,X(1,J), XDGM1(1,J))
CALL MULP(X(1,J),XDGM1(1,J), XDG(1,J))
END DO
C
C COMPUTE THE PRODUCT OF PDG AND XDGM1
DO J=1,N
CALL MULP( PDG(1,J), XDGM1(1,J), PXDGM1(1,J) )
END DO
C
C COMPUTE THE PRODUCT OF PDG AND XDG
DO J=1,N
CALL MULP( PDG(1,J), XDG(1,J), PXDG(1,J) )
END DO
G = PXDG - QDG
DO J=1,N
DG(1:2,J)= IDEG(J)*PXDGM1(1:2,J)
END DO
RETURN
END SUBROUTINE GFUNP
SUBROUTINE GMFADS(NN,A,NWK,MAXA)
C
C This subroutine computes the LDU decomposition of a symmetric positive
C definite matrix B where only the upper triangular skyline structure
C is stored. The decomposition is done by the Gill-Murray
C strategy from P.E. Gill and W. Murray, Newton type Methods
C for Unconstrained and Linearly Constrained Optimization,
C Mathematical Programming, 7, 311-350 (1974) and gives an
C approximate decomposition in the case of a nonpositive
C definite or ill-conditioned matrix.
C
C Input variables:
C
C NN -- dimension of B.
C
C A -- one dimensional real array containing the upper
C triangular skyline portion of a symmetric matrix B in
C packed skyline storage format.
C
C NWK -- number of elements in A.
C
C MAXA -- an integer array of dimension NN+1 containing the
C locations of the diagonal elements of B in A.
C By convention, MAXA(NN+1)=NWK+1.
C
C Output variables:
C
C A -- the upper triangular skyline portion of the LDU
C decomposition of the symmetric matrix B (or B + E if B
C was not sufficiently positive definite).
C
C
C No working storage is required by this routine.
C
C
USE REAL_PRECISION
INTEGER, INTENT(IN):: NN,NWK,MAXA(NN+1)
REAL (KIND=R8), INTENT(IN OUT):: A(NWK)
INTEGER:: I,I0,I1,I2,I3,I4,J,J1,K,K1,K2,KH,KL,KN,KU,KZ,L,L1,
& L2,L3,M,M1,N1,NNN
REAL (KIND=R8):: BET,DEL,DJ,G,GAM,GAM1,PHI,
& THE,THE1,XT1,XT2,ZET,ZET1
G=0.0
GAM=0.0
DO I=1,NN
K=MAXA(I)
G=G+A(K)*A(K)
GAM1=ABS(A(K))
IF(GAM1.GT.GAM)GAM=GAM1
END DO
ZET=0.0
DO I=1,NN
K=MAXA(I)
K1=MAXA(I+1)-1
K2=K1-K
IF (K2.EQ.0) CYCLE
L=K+1
DO J=L,K1
G=G+2.0*A(J)*A(J)
ZET1=ABS(A(J))
IF(ZET1.GT.ZET)ZET=ZET1
END DO
END DO
ZET=ZET/NN
DEL=EPSILON(1.0_R8)
BET=DEL
IF (ZET .GT. BET) BET=ZET
IF (GAM .GT. BET) BET=GAM
G=SQRT(G)
IF (G .GT. 1.0) DEL=DEL*G
DO I=1,NN
N1=I-1
KN=MAXA(I)
KL=KN+1
KU=MAXA(I+1)-1
KH=KU-KL
PHI=A(KN)
IF (KH .GE. 0) THEN
K1=KN+1
K2=I
DO J=K1,KU
K2=K2-1
KZ=MAXA(K2)
PHI=PHI-A(J)*A(J)*A(KZ)
END DO
END IF
PHI=ABS(PHI)
L=I+1
THE=0.0
NNN=NN+1
IF (L .NE. NNN) THEN
DO J=L,NN
L1=MAXA(J)
L2=MAXA(J+1)
L3=L2-L1-1
M=J-I
IF (L3 .LT. M) CYCLE
M1=L1+M
IF (N1 .NE. 0) THEN
DO J1=1,N1
I0=MAXA(J1)
I1=MAXA(L)
I2=I-J1
I3=I1-KN-1
I4=J-J1
IF (I3 .LT. I2) CYCLE
IF (L3 .LT. I4) CYCLE
XT1=A(KN+I2)
XT2=A(L1+I4)
A(M1)=A(M1)-XT1*XT2*A(I0)
END DO
END IF
THE1=ABS(A(M1))
IF (THE .LT. THE1) THE=THE1
END DO
END IF
THE=THE*THE/BET
DJ=DEL
IF (PHI .GT. DJ) DJ=PHI
IF (THE .GT. DJ) DJ=THE
A(KN)=DJ
IF (L .EQ. NNN) CYCLE
DO J=L,NN
L1=MAXA(J)
L2=MAXA(J+1)
L3=L2-L1-1
M=J-I
IF (L3 .LT. M) CYCLE
M1=L1+M
A(M1)=A(M1)/A(KN)
END DO
END DO
RETURN
END SUBROUTINE GMFADS
SUBROUTINE GMRES(N, KDMAX, ITMAX, RHS, X, KVAL,
& PRECON, IFLAG, ROWPOSP, COLPOSP)
C
C THIS ROUTINE IS AN EXTENSION OF THE STANDARD RESTARTED GMRES METHOD OF
C Y. SAAD AND M. SCHULTZ, "GMRES: A GENERALIZED MINIMAL RESIDUAL METHOD
C FOR SOLVING NONSYMMETRIC LINEAR SYSTEMS", SIAM J. SCI. STAT. COMPUT., 7
C (1986), PP. 856-869, THAT ALLOWS THE MAXIMUM KRYLOV SUBSPACE DIMENSION
C TO INCREASE (UP TO A MAXIMUM SPECIFIED PARAMETER VALUE) IF STAGNATION
C OCCURS. THE ARNOLDI BASIS VECTORS ARE GENERATED USING ORTHOGONALIZATION
C WITH HOUSEHOLDER TRANSFORMATIONS.
C ON RESTARTS, RESIDUAL VECTORS ARE COMPUTED BY DIRECT EVALUATION.
C CONDITIONING OF THE GMRES LEAST-SQUARES PROBLEM IS MONITORED USING
C LAPACK ROUTINE DLAIC1 (SEE P. N. BROWN AND H. F. WALKER, "GMRES ON
C (NEARLY) SINGULAR SYSTEMS", UTAH STATE UNIV. C MATH. STAT. DEPT. RES.
C REPORT 2/94/73, FEBRUARY, 1994).
C
C
C INPUT VARIABLES:
C
C N IS THE DIMENSION OF THE LINEAR SYSTEM AA*X = RHS. THE SPARSE
C MATRIX DATA STRUCTURE FOR AA IS STORED IN THE MODULE HOMOTOPY.
C
C KDMAX IS THE MAXIMUM KRYLOV SUBSPACE DIMENSION ALLOWED BEFORE
C STAGNATION DETECTION BEGINS.
C
C ITMAX IS THE MAXIMUM OVERALL NUMBER OF GMRES ITERATIONS ALLOWED.
C
C RHS IS THE RIGHT HAND SIDE OF THE LINEAR SYSTEM.
C
C X IS THE INITIAL APPROXIMATE SOLUTION OF THE LINEAR SYSTEM.
C
C KVAL IS AN INDEX INTO X REQUIRED FOR MATRIX-VECTOR MULTIPLICATION.
C
C PRECON IS A ONE-DIMENSIONAL REAL ARRAY CONTAINING A PRECONDITIONING
C MATRIX Q FOR AA.
C
C IFLAG INDICATES WHAT DATA IS USED IN MATRIX-VECTOR MULTIPLICATION.
C
C ROWPOSP IS AN OPTIONAL INTEGER ARRAY USED IN THE DATA STRUCTURE
C DESCRIBING THE SPARSE PRECONDITIONING MATRIX Q.
C
C COLPOSP IS AN OPTIONAL INTEGER ARRAY USED IN THE DATA STRUCTURE
C DESCRIBING THE SPARSE PRECONDITIONING MATRIX Q.
C
C OUTPUT VARIABLES:
C
C KDMAX IS THE INCREASED MAXIMUM KRYLOV SUBSPACE DIMENSION ON
C TERMINATION.
C
C IFLAG IS UNCHANGED ON NORMAL TERMINATION (DESIRED TOLERANCE MET);
C = 4, IF GMRES ITERATION FAILS TO CONVERGE BECAUSE
C - DESIRED TOLERANCE NOT MET AFTER ITMAX ITERATIONS;
C - DANGEROUS NEAR-SINGULARITY DETECTED;
C - LIMITS OF NUMERICAL ACCURACY REACHED;
C - STAGNATION OCCURS.
C
C ITMAX IS THE TOTAL NUMBER OF GMRES ITERATIONS PERFORMED.
C
C X IS THE FINAL APPROXIMATE SOLUTION.
C
C
C CALLS DNRM2, ILUSOLVDS, MULTDS, MULT2DS, SOLVDS, AND INTERNAL
C SUBROUTINES MULM1 AND MULM2 AND INTERNAL FUNCTIONS APPL_HOUSE, HOUSE.
C
USE HOMOTOPY, AA=>QRSPARSE
USE REAL_PRECISION
INTEGER, INTENT(IN):: KVAL, N
INTEGER, INTENT(IN OUT):: IFLAG, ITMAX, KDMAX
REAL (KIND=R8), INTENT(IN):: PRECON(:), RHS(N)
REAL (KIND=R8), INTENT(IN OUT):: X(N)
INTEGER, INTENT(IN), OPTIONAL:: COLPOSP(:), ROWPOSP(:)
C
C IPRINT IS A PARAMETER FOR THE FORTRAN UNIT NUMBER, WHICH IF POSITIVE
C ACTIVATES TRACING OF THE GMRES ITERATIONS.
C M IS A PARAMETER USED IN THE ADAPTIVE STRATEGY FOR INCREMENTING THE
C KRYLOV SUBSPACE DIMENSION (VALUE KDMAX).
C STRONG_VERSION IS A LOGICAL PARAMETER, WHICH IF TRUE CAUSES GMRES TO
C STOP IN TWO CASES: WHEN THE ESTIMATED CONDITION NUMBER IS GREATER
C THAN 1.0/(50*EPSILON);
C WHEN THE TRUE RESIDUAL NORM IS LARGE AND IT INCREASED DURING THE LAST
C RESTART CYCLE.
C
INTEGER, PARAMETER:: IPRINT=0, M=4
LOGICAL, PARAMETER:: STRONG_VERSION=.TRUE.
C
C LOCAL VARIABLES.
C
INTEGER:: I, IFLAGC, IFLAGI, IQUIT, ITNO, KD,
& KDP1, KDLIMIT, LENAA
REAL (KIND=R8), ALLOCATABLE:: C(:), R(:,:), S(:), SVBIG(:),
& SVSML(:), V(:,:), W(:)
REAL (KIND=R8):: VTEMP(N), VTEMP2(N)
REAL (KIND=R8):: BIG, BIGCND, CC, CNDMAX,
& RSN, RSNOLD, SESTPR, SMALL, SS, TEMP, TOL
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 DEFINE MAXIMUM ALLOWED KRYLOV SUBSPACE DIMENSION.
C
KDLIMIT = MIN(MAX(200, INT(SQRT(REAL(N)))), (N+1)/2)
KDMAX = MIN(KDMAX,KDLIMIT)
C
C ALLOCATE LOCAL ARRAYS.
C
IF(.NOT. ALLOCATED(C)) ALLOCATE(C(KDLIMIT))
IF(.NOT. ALLOCATED(R)) ALLOCATE(R(KDLIMIT,KDLIMIT))
IF(.NOT. ALLOCATED(S)) ALLOCATE(S(KDLIMIT))
IF(.NOT. ALLOCATED(SVBIG)) ALLOCATE(SVBIG(KDLIMIT))
IF(.NOT. ALLOCATED(SVSML)) ALLOCATE(SVSML(KDLIMIT))
IF(.NOT. ALLOCATED(V)) ALLOCATE(V(N,KDLIMIT+1))
IF(.NOT. ALLOCATED(W)) ALLOCATE(W(KDLIMIT+1))
C
C PERFORM INITIALIZATIONS.
C
ITNO = 0
IFLAGC = 0
IFLAGI = 1
LENAA = ROWPOS(N)-1
TOL = MAX(100.0, 1.01*REAL(LENAA)/REAL(N))*EPSILON(1.0_R8)
VTEMP = X
IF (PRESENT(ROWPOSP) .AND. PRESENT( COLPOSP)) THEN
CALL MULM2(V(:,1),VTEMP) ! MODE=2
ELSE
CALL MULM1(V(:,1),VTEMP) ! MODE=1
END IF
V(:,1) = RHS - V(:,1)
RSN = DNRM2(N, V(:,1), 1)
TEMP = DNRM2(N, RHS, 1)
TOL = MAX(RSN,TEMP)*TOL
IF (RSN .LE. TOL) THEN
ITMAX = ITNO
CALL CLEANUP
RETURN
ELSE
IF (ITMAX .EQ. 0) THEN
ITMAX = ITNO
IFLAG = 4
CALL CLEANUP
RETURN
ENDIF
ENDIF
CNDMAX = 1.0/(50.0*EPSILON(1.0_R8))
IQUIT = 0
BIGCND = 0.0
C
C FOR PRINTING:
C
IF (IPRINT .GT. 0) THEN
WRITE(IPRINT, 19) ITNO, RSN
19 FORMAT(//,9X,'GMRES ITERATIONS:',//,T4,'ITERATION',TR4,
& 'RESIDUAL NORM',TR5,'CONDITION NUMBER',/,T4,I6,TR7,ES16.8)
ENDIF
C
C BEGIN THE OUTER LOOP.
C
OUTER: DO
C
C FOR PRINTING:
C
IF (IPRINT .GT. 0 ) THEN
IF (IFLAGI .NE. 0) THEN
WRITE (IPRINT, 21) KDMAX
21 FORMAT(/' RESTART WITH KRYLOV SUBSPACE DIMENSION =',I2)
ELSE
WRITE(IPRINT, 22)
22 FORMAT(/' RESTART')
ENDIF
ENDIF
C
KD = 0
RSNOLD = RSN
IQUIT = 0
IFLAGI = 0
C
C FIND HOUSEHOLDER VECTOR V(:,1).
C
W(1)=V(1,1)
CALL HREFG(N,V(2:N,1),V(1,1),W(1))
C
C BEGIN THE INNER LOOP.
C
200 INNER: DO
KD = KD + 1
KDP1 = KD + 1
ITNO = ITNO + 1
IQUIT = 0
C
C FIND THE KD-TH ORTHOGONAL VECTOR.
C
VTEMP = 0.0
VTEMP(KD) = 1.0
DO I=KD,1,-1
TEMP=V(I,I)
V(I,I)=1.0
CALL HREFX(N-I+1,V(I:N,I),TEMP,VTEMP(I:N))
V(I,I)=TEMP
END DO
C
C EVALUATE AA*(KD-TH KRYLOV SUBSPACE BASIS VECTOR).
C
IF (PRESENT(ROWPOSP) .AND. PRESENT(COLPOSP)) THEN ! MODE = 2
CALL ILUSOLVDS(N, PRECON(1:ROWPOSP(N+1)-1),
& ROWPOSP(N+1)-1, ROWPOSP(1:N+1),
& COLPOSP(1:ROWPOSP(N+1)-1), VTEMP(1:N))
CALL MULM2(VTEMP2, VTEMP(1:N))
ELSE ! MODE = 1
CALL SOLVDS(N, PRECON, ROWPOS(N+1)-1, ROWPOS(1:N+1),
& VTEMP(1:N))
CALL MULM1(VTEMP2, VTEMP(1:N))
ENDIF
C
C FIND HOUSEHOLDER VECTOR V(:,KDP1).
C
VTEMP(1:N)=VTEMP2(1:N)
DO I=1,KD
TEMP=V(I,I)
V(I,I)=1.0
CALL HREFX(N-I+1,V(I:N,I),TEMP,VTEMP(I:N))
V(I,I)=TEMP
END DO
IF ( MAXVAL(ABS(VTEMP(KDP1:N))) .NE. 0.0) THEN
V(KDP1+1:N, KDP1) = VTEMP(KDP1+1:N)
CALL HREFG(N-KDP1+1,V(KDP1+1:N, KDP1),V(KDP1,KDP1),
& VTEMP(KDP1))
ENDIF
C
C IF KD .GT. 1, APPLY THE PREVIOUS GIVENS ROTATIONS.
C
DO I = 1, KD-1
TEMP = VTEMP(I)
VTEMP(I) = C(I)*TEMP + S(I)*VTEMP(I+1)
VTEMP(I+1) = -S(I)*TEMP + C(I)*VTEMP(I+1)
END DO
C
C DETERMINE AND APPLY THE NEXT ROTATION.
C
IF (VTEMP(KDP1) .NE. 0.0) THEN
TEMP = VTEMP(KD)
R(KD,KD) = SQRT(VTEMP(KD)**2 + VTEMP(KDP1)**2)
C(KD) = TEMP/R(KD,KD)
S(KD) = VTEMP(KDP1)/R(KD,KD)
VTEMP(KD) = C(KD)*TEMP + S(KD)*VTEMP(KDP1)
END IF
R(1:KD, KD) = VTEMP(1:KD)
C
C COMPUTE AND TEST INCREMENTAL CONDITION NUMBER.
C
IF (KD .EQ. 1) THEN
BIG = R(KD,KD)
SMALL = BIG
SVBIG(1) = 1.0
SVSML(1) = 1.0
ELSE
I = 1
CALL DLAIC1(I, KD-1, SVBIG, BIG, R(1,KD), R(KD,KD),
& SESTPR, SS, CC)
BIG = SESTPR
SVBIG(1:KD-1) = SS*SVBIG(1:KD-1)
SVBIG(KD) = CC
I = 2
CALL DLAIC1(I, KD-1, SVSML, SMALL, R(1,KD), R(KD,KD),
& SESTPR, SS, CC)
SMALL = SESTPR
SVSML(1:KD-1) = SS*SVSML(1:KD-1)
SVSML(KD) = CC
ENDIF
IF (STRONG_VERSION) THEN
IF (BIG .GE. SMALL*CNDMAX) THEN
IF (IPRINT .GT. 0) THEN
WRITE (IPRINT, 230) CNDMAX
230 FORMAT(/,4X, 'IN GMRES CONDITION NUMBER .GE.',
& ES16.8)
ENDIF
ITNO = ITNO - 1
IFLAGC = 4
EXIT OUTER
ENDIF
ENDIF
TEMP = BIG/SMALL
IF (TEMP .GT. BIGCND) BIGCND = TEMP
C
C UPDATE W AND THE RESIDUAL NORM.
C
IF (VTEMP(KDP1) .NE. 0.0) THEN
W(KDP1)= -S(KD)*W(KD)
W(KD) = C(KD)*W(KD)
ELSE
W(KDP1) = 0.0
ENDIF
RSN = ABS(W(KDP1))
C
C FOR PRINTING:
C
IF (IPRINT .GT. 0 ) THEN
WRITE (IPRINT, 240) ITNO, RSN, BIG/SMALL
240 FORMAT(T4,I6, TR7,2ES16.8)
ENDIF
C
C TEST FOR TERMINATION OF THE INNER LOOP.
C
IF (RSN .LE. TOL .OR. KD .EQ. KDMAX .OR. ITNO .EQ. ITMAX)
& EXIT INNER
END DO INNER
C
C IF TERMINATING THE INNER LOOP, TEST FOR TERMINATION OF THE OUTER LOOP,
C COMPUTE THE CORRECTION, AND UPDATE THE APPROXIMATE SOLUTION.
C
C TEST FOR TERMINATION OF THE OUTER LOOP.
C
IF (RSN .LE. TOL .OR. ITNO .EQ. ITMAX ) THEN
IQUIT = 1
C
C TEST FOR STAGNATION.
C
ELSE
TEMP = KD*LOG(TOL/RSN)/
& LOG(RSN/((1.0 + 10.0*EPSILON(1.0_R8))*RSNOLD))
IF (TEMP .GE. 40.0*(ITMAX - ITNO)) THEN
IF (KDMAX .LE. KDLIMIT-M) THEN
IFLAGI = 3
KDMAX = KDMAX + M ! INCREASE KDMAX BY M
GO TO 200
ELSE IF (KDMAX .NE. KDLIMIT) THEN
IFLAGI = 3
KDMAX = KDLIMIT
GO TO 200
END IF
ENDIF
ENDIF
C
C COMPUTE THE CORRECTION IN V(:,1).
C
DO I = KD, 1, -1
W(I) = W(I)/R(I,I)
IF (I .GT. 1) W(1:I-1) = W(1:I-1)-W(I)*R(1:I-1,I)
END DO
C
C COMPUTE KD ORTHOGONAL VECTORS FROM HOUSEHOLDER VECTORS.
C
VTEMP = 0.0
VTEMP(1:KD)=W(1:KD)
DO I=KD,1,-1
TEMP=V(I,I)
V(I,I)=1.0
CALL HREFX(N-I+1,V(I:N,I), TEMP, VTEMP(I:N))
V(I,I)=TEMP
END DO
C
C ADD THE CORRECTION TO THE APPROXIMATE SOLUTION.
C
IF (PRESENT(ROWPOSP) .AND. PRESENT(COLPOSP)) THEN ! MODE = 2
CALL ILUSOLVDS(N, PRECON(1:ROWPOSP(N+1)-1),
& ROWPOSP(N+1)-1, ROWPOSP(1:N+1),
& COLPOSP(1:ROWPOSP(N+1)-1), VTEMP(1:N))
ELSE ! MODE = 1
CALL SOLVDS(N, PRECON, ROWPOS(N+1)-1, ROWPOS(1:N+1),
& VTEMP(1:N))
ENDIF
X = X+VTEMP
C
C UPDATE THE RESIDUAL VECTOR BY DIRECT EVALUATION IN (V:,1) AND RECOMPUTE
C THE RESIDUAL NORM.
C
IF (IQUIT .EQ. 0 ) THEN
VTEMP=X
IF (PRESENT(ROWPOSP) .AND. PRESENT( COLPOSP)) THEN
CALL MULM2(V(:,1),VTEMP) ! MODE=2
ELSE
CALL MULM1(V(:,1),VTEMP) ! MODE=1
END IF
V(:,1) = RHS - V(:,1)
RSN = DNRM2(N, V(:,1), 1)
ENDIF
C
C TERMINATE, OR GO TO THE TOP OF THE OUTER LOOP.
C
IF (RSN .LE. TOL) THEN
IFLAGC = 0
EXIT OUTER
ENDIF
IF (ITNO .EQ. ITMAX) THEN
IF (IPRINT .GT. 0) THEN
WRITE (IPRINT, 250) ITMAX
250 FORMAT(/,4X, 'MAXIMUM NUMBER OF ITERATIONS REACHED:',I7)
ENDIF
IFLAGC = 4
EXIT OUTER
END IF
IF (IQUIT .EQ. 0 ) THEN
IF (RSN .GT. RSNOLD ) THEN
IF (IPRINT .GT. 0) THEN
WRITE (IPRINT, 260) RSN
260 FORMAT(/,4X, '(TRUE) RESIDUAL NORM INCREASED TO',ES16.8)
ENDIF
IF (RSN .LE. TOL**(2.0/3.0)) THEN
IFLAGC = 0
ELSE IF (STRONG_VERSION) THEN
IF (IPRINT .GT. 0) THEN
WRITE (IPRINT, 270) TOL**(2.0/3.0)
270 FORMAT(/,4X,'(TRUE) RESIDUAL NORM IS LARGER THAN',ES16.8)
ENDIF
IFLAGC = 4
END IF
EXIT OUTER
END IF
C
C TEST FOR STAGNATION USING TRUE RESIDUAL NORM.
C
TEMP = KD*LOG(TOL/RSN)/
& LOG(RSN/((1.0 + 10.0*EPSILON(1.0_R8))*RSNOLD))
IF (TEMP .GE. 70.0*(ITMAX - ITNO)) THEN
IFLAGI = 3
IF ( KDMAX .LE. KDLIMIT-M ) THEN
KDMAX = KDMAX + M ! INCREASE KDMAX
ELSE IF (KDMAX .NE. KDLIMIT) THEN
KDMAX = KDLIMIT
END IF
ELSE IF (TEMP .GE. 1000.0*(ITMAX - ITNO)) THEN
IFLAGC = 4
EXIT OUTER
END IF
END IF
END DO OUTER
C
C RETURN.
C
ITMAX = ITNO
IF (IFLAGC .NE. 0) IFLAG = IFLAGC
CALL CLEANUP
RETURN
C
C INTERNAL SUBROUTINES.
C
CONTAINS
SUBROUTINE MULM1(Y,X)
C
C MATRIX-VECTOR MULTIPLY FOR MODE = 1.
C
REAL (KIND=R8), INTENT(IN):: X(N)
REAL (KIND=R8), INTENT(OUT):: Y(N)
Y = 0.0
CALL MULTDS(Y(1:N-1), AA, X(1:N-1), ROWPOS(1:N),
& N-1, ROWPOS(N)-1)
C
C RESULT MODIFIED ACCORDING TO KVAL.
C
Y(KVAL) = Y(KVAL)+X(N)
IF (KVAL .LT. N)
& Y(N) = X(KVAL) + X(N)*(1.0/ABS(AA(ROWPOS(KVAL)))+1.0)
RETURN
END SUBROUTINE MULM1
C
SUBROUTINE MULM2(Y,X)
C
C MATRIX-VECTOR MULTIPLY FOR MODE = 2.
C
REAL (KIND=R8), INTENT(IN):: X(N)
REAL (KIND=R8), INTENT(OUT):: Y(N)
INTERFACE
SUBROUTINE MULT2DS(Y, B, X, ROWPOS, COLPOS, N, LENB)
USE REAL_PRECISION
INTEGER, INTENT (IN):: LENB, N, ROWPOS(N+1),
& COLPOS(LENB)
REAL (KIND=R8), INTENT(IN):: X(:), B(LENB)
REAL (KIND=R8), INTENT (OUT):: Y(N)
END SUBROUTINE MULT2DS
END INTERFACE
Y = 0.0
IF (IFLAG .NE. -2) THEN
CALL MULT2DS(Y(1:N-1), AA, X(1:N-1), ROWPOS(1:N),
& COLPOS(1:LENAA), N-1, LENAA)
Y(1:N-1) = Y(1:N-1)-X(N)*PP(1:N-1)
ELSE
CALL MULT2DS(Y(1:N-1), AA, X(1:N), ROWPOS(1:N),
& COLPOS(1:LENAA), N-1, LENAA)
END IF
C
C RESULT MODIFIED ACCORDING TO KVAL.
C
Y(N) = X(KVAL)
RETURN
END SUBROUTINE MULM2
SUBROUTINE CLEANUP ! DEALLOCATE TEMPORARY LOCAL STORAGE.
IF(ALLOCATED(C)) DEALLOCATE(C)
IF(ALLOCATED(R)) DEALLOCATE(R)
IF(ALLOCATED(S)) DEALLOCATE(S)
IF(ALLOCATED(SVBIG)) DEALLOCATE(SVBIG)
IF(ALLOCATED(SVSML)) DEALLOCATE(SVSML)
IF(ALLOCATED(V)) DEALLOCATE(V)
IF(ALLOCATED(W)) DEALLOCATE(W)
END SUBROUTINE CLEANUP
C
SUBROUTINE HREFG(N, X, TAU, ALPHA)
C
C HOUSEHOLDER VECTOR CONSTRUCTION. HREFG IS THE LAPACK AUXILIARY
C ROUTINE DLARFG WITH TRIVIAL MODIFICATIONS FOR COMPATIBILITY WITH
C HOMPACK90.
C
INTEGER, INTENT(IN):: N
REAL (KIND=R8), INTENT(OUT):: ALPHA, TAU
REAL (KIND=R8), INTENT(IN OUT):: X(N-1)
!
! Purpose
! =======
!
! HREFG generates a real elementary reflector H of order N, such
! that
!
! H * ( ALPHA ) = ( BETA ), H' * H = I.
! ( X ) ( 0 )
!
! where ALPHA and BETA are scalars, and X is an (N-1)-element real
! vector. H is represented in the form
!
! H = I - TAU * ( 1 ) * ( 1 V' ) ,
! ( V )
!
! where TAU is a real scalar and V is a real (N-1)-element
! vector.
!
! If the elements of X are all zero, then TAU = 0 and H is taken to be
! the unit matrix.
!
! Otherwise 1 <= TAU <= 2.
!
! Arguments
! =========
!
! N (input)
! The order of the elementary reflector.
!
!
! X (input/output)
! On entry, the vector X.
! On exit, it is overwritten with the vector V.
!
! TAU (output)
! The value TAU.
!
! ALPHA (INPUT/OUTPUT)
! On entry, the value ALPHA.
! On exit, it is overwritten with the value BETA.
!
! =====================================================================
!
! .. Parameters ..
REAL (KIND=R8), PARAMETER:: ONE = 1.0_R8, ZERO = 0.0_R8
! .. Local Scalars ..
INTEGER:: J, KNT
REAL (KIND=R8):: BETA, RSAFMN, SAFMIN, XNORM
!
! .. External Functions ..
INTERFACE
FUNCTION DLAMCH(CMACH)
USE REAL_PRECISION
CHARACTER (LEN=1):: CMACH
REAL (KIND=R8):: DLAMCH
END FUNCTION DLAMCH
FUNCTION DLAPY2(X,Y)
USE REAL_PRECISION
REAL (KIND=R8):: DLAPY2,X,Y
END FUNCTION DLAPY2
FUNCTION DNRM2(N,X,STRIDE)
USE REAL_PRECISION
INTEGER:: N,STRIDE
REAL (KIND=R8):: DNRM2,X(N)
END FUNCTION DNRM2
END INTERFACE
!
! .. Executable Statements ..
!
IF( N.LE.1 ) THEN
TAU = ZERO
RETURN
END IF
XNORM = DNRM2( N-1, X, 1)
IF( XNORM.EQ.ZERO ) THEN
!
! H = I
!
TAU = ZERO
ELSE
!
! General case.
!
BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
SAFMIN = DLAMCH( 'S' )
IF( ABS( BETA ).LT.SAFMIN ) THEN
!
! XNORM, BETA may be inaccurate; scale X and recompute them.
!
RSAFMN = ONE / SAFMIN
KNT = 0
DO WHILE (ABS(BETA) < SAFMIN)
KNT = KNT + 1
X=RSAFMN*X
BETA = BETA*RSAFMN
ALPHA = ALPHA*RSAFMN
END DO
!
! New BETA is at most 1, at least SAFMIN.
!
XNORM = DNRM2( N-1, X, 1 )
BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
TAU = ( BETA-ALPHA ) / BETA
X = X/( ALPHA-BETA )
!
! If ALPHA is subnormal, it may lose relative accuracy.
!
ALPHA = BETA
DO J = 1, KNT
ALPHA = ALPHA*SAFMIN
END DO
ELSE
TAU = ( BETA-ALPHA) / BETA
X = X/( ALPHA-BETA )
ALPHA = BETA
END IF
END IF
RETURN
END SUBROUTINE HREFG
C
SUBROUTINE HREFX(M, V, TAU, C)
C
C HOUSEHOLDER VECTOR APPLICATION. HREFX IS A FORTRAN 90 VERSION OF
C THE LAPACK ROUTINE DLARFX WITH MODIFICATIONS FOR
C COMPATIBILITY WITH HOMPACK 90.
C
INTEGER, INTENT(IN):: M
REAL (KIND=R8), INTENT(IN):: TAU, V(M)
REAL (KIND=R8), INTENT(IN OUT):: C(M)
!
! Purpose
! =======
!
! HREFX applies a real elementary reflector H to a real M-element
! vector C. H is represented in the form
!
! H = I - TAU * V * V'
!
! where TAU is a real scalar and V is a real vector.
!
! If TAU = 0, then H is taken to be the unit matrix.
!
! This version uses inline code if H has order < 11.
!
! Arguments
! =========
!
! M (input)
! The order of the vector C.
!
! V (input)
! The vector V in the representation of H.
!
! TAU (input)
! The value TAU in the representation of H.
!
! C (input/output)
! On entry, the M-element vector C.
! On exit, C is overwritten by H * C.
!
! =====================================================================
!
! .. Parameters ..
REAL (KIND=R8), PARAMETER:: ONE=1.0_R8, ZERO=0.0_R8
!
! .. Local Scalars ..
REAL (KIND=R8):: SUM, T1, T10, T2, T3, T4, T5, T6, T7, T8, T9,
& V1, V10, V2, V3, V4, V5, V6, V7, V8, V9
!
! .. Executable Statements ..
!
IF( TAU.EQ.ZERO ) RETURN
!
! Form H * C, where H has order M.
!
SELECT CASE(M)
!
! Code for general M.
!
CASE (11:)
!
! C := C - TAU * V * V' * C
!
C = C - (TAU * DOT_PRODUCT(V,C)) * V
!
! Special code for 1 x 1 Householder.
!
CASE(1)
T1 = ONE - TAU*V( 1 )*V( 1 )
C(1) = T1*C(1)
!
! Special code for 2 x 2 Householder.
!
CASE(2)
V1 = V( 1 )
T1 = TAU*V1
V2 = V( 2 )
T2 = TAU*V2
SUM = V1*C(1) + V2*C(2)
C(1) = C(1) - SUM*T1
C(2) = C(2) - SUM*T2
!
! Special code for 3 x 3 Householder.
!
CASE(3)
V1 = V( 1 )
T1 = TAU*V1
V2 = V( 2 )
T2 = TAU*V2
V3 = V( 3 )
T3 = TAU*V3
SUM = V1*C(1) + V2*C(2) + V3*C(3)
C(1) = C(1) - SUM*T1
C(2) = C(2) - SUM*T2
C(3) = C(3) - SUM*T3
!
! Special code for 4 x 4 Householder.
!
CASE(4)
V1 = V( 1 )
T1 = TAU*V1
V2 = V( 2 )
T2 = TAU*V2
V3 = V( 3 )
T3 = TAU*V3
V4 = V( 4 )
T4 = TAU*V4
SUM = V1*C(1) + V2*C(2) + V3*C(3) + V4*C(4)
C( 1 ) = C( 1 ) - SUM*T1
C( 2 ) = C( 2 ) - SUM*T2
C( 3 ) = C( 3 ) - SUM*T3
C( 4 ) = C( 4 ) - SUM*T4
!
! Special code for 5 x 5 Householder.
!
CASE(5)
V1 = V( 1 )
T1 = TAU*V1
V2 = V( 2 )
T2 = TAU*V2
V3 = V( 3 )
T3 = TAU*V3
V4 = V( 4 )
T4 = TAU*V4
V5 = V( 5 )
T5 = TAU*V5
SUM = V1*C( 1 ) + V2*C( 2 ) + V3*C( 3 ) +
& V4*C( 4 ) + V5*C( 5 )
C( 1 ) = C( 1 ) - SUM*T1
C( 2 ) = C( 2 ) - SUM*T2
C( 3 ) = C( 3 ) - SUM*T3
C( 4 ) = C( 4 ) - SUM*T4
C( 5 ) = C( 5 ) - SUM*T5
!
! Special code for 6 x 6 Householder.
!
CASE(6)
V1 = V( 1 )
T1 = TAU*V1
V2 = V( 2 )
T2 = TAU*V2
V3 = V( 3 )
T3 = TAU*V3
V4 = V( 4 )
T4 = TAU*V4
V5 = V( 5 )
T5 = TAU*V5
V6 = V( 6 )
T6 = TAU*V6
SUM = V1*C( 1 ) + V2*C( 2 ) + V3*C( 3 ) +
& V4*C( 4 ) + V5*C( 5 ) + V6*C( 6 )
C( 1 ) = C( 1 ) - SUM*T1
C( 2 ) = C( 2 ) - SUM*T2
C( 3 ) = C( 3 ) - SUM*T3
C( 4 ) = C( 4 ) - SUM*T4
C( 5 ) = C( 5 ) - SUM*T5
C( 6 ) = C( 6 ) - SUM*T6
!
! Special code for 7 x 7 Householder.
!
CASE(7)
V1 = V( 1 )
T1 = TAU*V1
V2 = V( 2 )
T2 = TAU*V2
V3 = V( 3 )
T3 = TAU*V3
V4 = V( 4 )
T4 = TAU*V4
V5 = V( 5 )
T5 = TAU*V5
V6 = V( 6 )
T6 = TAU*V6
V7 = V( 7 )
T7 = TAU*V7
SUM = V1*C( 1 ) + V2*C( 2 ) + V3*C( 3 ) +
& V4*C( 4 ) + V5*C( 5 ) + V6*C( 6 ) +
& V7*C( 7 )
C( 1 ) = C( 1 ) - SUM*T1
C( 2 ) = C( 2 ) - SUM*T2
C( 3 ) = C( 3 ) - SUM*T3
C( 4 ) = C( 4 ) - SUM*T4
C( 5 ) = C( 5 ) - SUM*T5
C( 6 ) = C( 6 ) - SUM*T6
C( 7 ) = C( 7 ) - SUM*T7
!
! Special code for 8 x 8 Householder.
!
CASE(8)
V1 = V( 1 )
T1 = TAU*V1
V2 = V( 2 )
T2 = TAU*V2
V3 = V( 3 )
T3 = TAU*V3
V4 = V( 4 )
T4 = TAU*V4
V5 = V( 5 )
T5 = TAU*V5
V6 = V( 6 )
T6 = TAU*V6
V7 = V( 7 )
T7 = TAU*V7
V8 = V( 8 )
T8 = TAU*V8
SUM = V1*C( 1 ) + V2*C( 2 ) + V3*C( 3 ) +
& V4*C( 4 ) + V5*C( 5 ) + V6*C( 6 ) +
& V7*C( 7 ) + V8*C( 8 )
C( 1 ) = C( 1 ) - SUM*T1
C( 2 ) = C( 2 ) - SUM*T2
C( 3 ) = C( 3 ) - SUM*T3
C( 4 ) = C( 4 ) - SUM*T4
C( 5 ) = C( 5 ) - SUM*T5
C( 6 ) = C( 6 ) - SUM*T6
C( 7 ) = C( 7 ) - SUM*T7
C( 8 ) = C( 8 ) - SUM*T8
!
! Special code for 9 x 9 Householder.
!
CASE(9)
V1 = V( 1 )
T1 = TAU*V1
V2 = V( 2 )
T2 = TAU*V2
V3 = V( 3 )
T3 = TAU*V3
V4 = V( 4 )
T4 = TAU*V4
V5 = V( 5 )
T5 = TAU*V5
V6 = V( 6 )
T6 = TAU*V6
V7 = V( 7 )
T7 = TAU*V7
V8 = V( 8 )
T8 = TAU*V8
V9 = V( 9 )
T9 = TAU*V9
SUM = V1*C( 1 ) + V2*C( 2 ) + V3*C( 3 ) +
& V4*C( 4 ) + V5*C( 5 ) + V6*C( 6 ) +
& V7*C( 7 ) + V8*C( 8 ) + V9*C( 9 )
C( 1 ) = C( 1 ) - SUM*T1
C( 2 ) = C( 2 ) - SUM*T2
C( 3 ) = C( 3 ) - SUM*T3
C( 4 ) = C( 4 ) - SUM*T4
C( 5 ) = C( 5 ) - SUM*T5
C( 6 ) = C( 6 ) - SUM*T6
C( 7 ) = C( 7 ) - SUM*T7
C( 8 ) = C( 8 ) - SUM*T8
C( 9 ) = C( 9 ) - SUM*T9
!
! Special code for 10 x 10 Householder.
!
CASE (10)
V1 = V( 1 )
T1 = TAU*V1
V2 = V( 2 )
T2 = TAU*V2
V3 = V( 3 )
T3 = TAU*V3
V4 = V( 4 )
T4 = TAU*V4
V5 = V( 5 )
T5 = TAU*V5
V6 = V( 6 )
T6 = TAU*V6
V7 = V( 7 )
T7 = TAU*V7
V8 = V( 8 )
T8 = TAU*V8
V9 = V( 9 )
T9 = TAU*V9
V10 = V( 10 )
T10 = TAU*V10
SUM = V1*C( 1 ) + V2*C( 2 ) + V3*C( 3 ) +
& V4*C( 4 ) + V5*C( 5 ) + V6*C( 6 ) +
& V7*C( 7 ) + V8*C( 8 ) + V9*C( 9 ) +
& V10*C( 10 )
C( 1 ) = C( 1 ) - SUM*T1
C( 2 ) = C( 2 ) - SUM*T2
C( 3 ) = C( 3 ) - SUM*T3
C( 4 ) = C( 4 ) - SUM*T4
C( 5 ) = C( 5 ) - SUM*T5
C( 6 ) = C( 6 ) - SUM*T6
C( 7 ) = C( 7 ) - SUM*T7
C( 8 ) = C( 8 ) - SUM*T8
C( 9 ) = C( 9 ) - SUM*T9
C( 10 ) = C( 10 ) - SUM*T10
END SELECT
RETURN
END SUBROUTINE HREFX
END SUBROUTINE GMRES
SUBROUTINE GMRILUDS(NN,LENAA,IFLAG,START,RHS)
C
C This subroutine solves a system of equations using a
C preconditioned adaptive GMRES(k).
C GMRILUDS is the MODE=2 equivalent of subroutine PCGDS, which
C is only for MODE=1 storage.
C
C The linear system to be solved is Mx=b.
C If IFLAG = -1 or 0,
C
C +-- --+
C | | |
C | AA | -PP |
C M = | | | ,
C +--------+-----+
C | E(k)**t |
C +-- --+
C
C where AA is an (NN x NN) matrix stored in compressed-row format,
C PP is an (NN x 1) vector.
C It is assumed that rank [AA,-PP]=NN and M is invertible.
C If IFLAG = -2,
C
C +-- --+
C | |
C | AA |
C M = | | ,
C +--------+-----+
C | E(k)**t |
C +-- --+
C
C where AA is an (NN x (NN+1)) matrix stored in compressed-row
C format. It is assumed that rank [AA]=NN and M is invertible.
C
C +- -+ +- -+
C | 0 | | |
C | ... | | -RHS |
C b = | 0 |, or b = | |,
C +-----+ +------+
C | T | | T |
C +- -+ +- -+
C
C T = START(k), where |START(k)|= max|START(i)|
C 1<=i<=NN+1
C
C b is of length NN+1 and E(k)**t is the ( 1 x (NN+1) ) vector
C consisting of all zeros, except for a '1' in its k-th position.
C
C Input variables:
C
C NN -- row dimension of the matrix packed in AA.
C
C LENAA -- number of elements in the packed array AA.
C
C IFLAG -- indicator of how M is assembled.
C
C START -- vector of length NN+1, normally the solution to the
C previous linear system; used to determine the index k .
C
C RHS -- optional vector of length NN, used to define right hand
C side for normal flow correction calculation. It is
C assumed that GMRILUDS is called without RHS present before
C it is called with RHS present. An ILU factorization based
C preconditioner is computed only when RHS is not present.
C
C Input variables defined in module HOMOTOPY:
C
C AA -- one dimensional real array containing a submatrix of M
C in compressed-row storage form. The logical dimensions
C of AA depend on IFLAG.
C
C ROWPOS -- integer array used for specifying information about AA.
C Using compressed-row storage, it has length NN+2, and
C stores the indices of where each row begins within AA.
C ROWPOS(NN+1) = LENAA + 1 and ROWPOS(NN+2) = LENAA + 2.
C (NOTE: The value of ROWPOS(NN+2) is set by this
C subroutine when the preconditioning matrix Q is
C initialized.)
C
C COLPOS -- integer array used for specifying information about AA.
C Using compressed-row storage, it has length LENAA,
C and contains the column indices of the corresponding
C elements in AA.
C
C For example, using the compressed-row storage scheme
C with IFLAG = -2, a 5 x 6 matrix of the form
C
C +-- --+
C | 1 3 0 0 0 10 |
C | 3 2 0 7 0 11 |
C | 0 0 4 6 0 12 |
C | 0 7 6 5 9 13 |
C | 0 0 0 9 8 14 |
C +-- --+
C
C would result in NN=5, LENAA=18, ROWPOS=(1,4,8,11,16,19,*),
C AA=(1,3,10,3,2,7,11,4,6,12,7,6,5,9,13,9,8,14),
C COLPOS=(1,2,6,1,2,6,3,4,6,2,3,4,5,6,4,5,6).
C
C PP -- vector of length NN, used for (NN+1)st column of
C augmented matrix M when IFLAG = -1 or 0 .
C
C Output variables:
C
C START -- solution vector x of M x = b (defined above).
C
C IFLAG -- normally unchanged on output. If the GMRES
C iteration fails to converge in 10*(NN+1) iterations (most
C likely due to a singular Jacobian matrix), GMRILUDS returns
C with IFLAG = 4 , and does not compute x.
C
C Calls subroutines ILUFDS and GMRES.
C
USE HOMOTOPY, AA => QRSPARSE, WORK => PAR, IWORK => IPAR
USE REAL_PRECISION
C
INTEGER, INTENT(IN):: NN,LENAA
INTEGER, INTENT(IN OUT):: IFLAG
REAL (KIND=R8), INTENT(IN OUT):: START(NN+1)
REAL (KIND=R8), INTENT(IN), OPTIONAL :: RHS(NN)
C
C LOCAL VARIABLES.
C
INTEGER:: I, ITMAX, K, KD, NP1, QIND, ZBIND
INTEGER:: CIND, RIND, ROWL, STRT
REAL (KIND=R8):: STARTK
REAL (KIND=R8):: RHSC(NN+1)
C
C GMRES PARAMETER.
C
INTEGER, PARAMETER:: SUBSPACE=8 ! KRYLOV SUBSPACE VALUE.
C
INTERFACE
SUBROUTINE GMRES(N, KDMAX, ITMAX, RHSC, X, KVAL,
& Q, IFLAG, ROWPOSP, COLPOSP)
USE HOMOTOPY
USE REAL_PRECISION
INTEGER, INTENT(IN):: KVAL, N
INTEGER, INTENT(IN OUT):: IFLAG, ITMAX, KDMAX
REAL (KIND=R8), INTENT(IN):: Q(:), RHSC(N)
REAL (KIND=R8), INTENT(IN OUT):: X(N)
INTEGER, INTENT(IN), OPTIONAL:: COLPOSP(:), ROWPOSP(:)
END SUBROUTINE GMRES
SUBROUTINE ILUFDS(NN, Q, LENQ, ROWPOS, COLPOS)
USE REAL_PRECISION
INTEGER, INTENT(IN):: LENQ, NN, COLPOS(LENQ), ROWPOS(NN+1)
REAL (KIND=R8), INTENT(IN OUT):: Q(LENQ)
END SUBROUTINE ILUFDS
END INTERFACE
C
NP1=NN+1
C
C INITIALIZE START POSITIONS WITHIN WORK AND IWORK.
C
ZBIND = 1
QIND = NP1+1
RIND = 1
CIND = NP1+2
C
IF (.NOT. ALLOCATED(WORK)) THEN
IF (IFLAG .EQ. -2) THEN
ALLOCATE(WORK(NP1+LENAA+2))
ELSE
ALLOCATE(WORK(NP1+LENAA+NN+2))
END IF
WORK(1:NP1) = 0.0
END IF
IF (.NOT. ALLOCATED(IWORK)) THEN
IF (IFLAG .EQ. -2) THEN
ALLOCATE(IWORK(NP1+1+LENAA+2))
ELSE
ALLOCATE(IWORK(NP1+1+LENAA+NN+2))
END IF
END IF
C
C FIND THE ELEMENT OF LARGEST MAGNITUDE IN THE INITIAL VECTOR, AND
C RECORD ITS POSITION IN K.
C
K = MAXVAL(MAXLOC(ABS(START)))
STARTK = START(K)
C
C SET VALUES OF ROWPOS(NN+1) AND ROWPOS(NN+2), AND
C COMPUTE THE PRECONDITIONER Q FOR M.
C
IF (.NOT. PRESENT(RHS)) THEN
ROWPOS(NP1) = LENAA+1
ROWPOS(NN+2) = LENAA+2
IF (IFLAG .EQ. -2) THEN
WORK(QIND:QIND+LENAA-1) = AA(1:LENAA)
IWORK(RIND:RIND+NP1) = ROWPOS
IWORK(CIND:CIND+LENAA-1) = COLPOS
ELSE
C MERGE AA AND -PP INTO Q ONLY FOR IFLAG >= -1.
IWORK(RIND) = 1
DO I=1,NN
STRT = IWORK(RIND+I-1)
ROWL = ROWPOS(I+1)-ROWPOS(I)
WORK(QIND+STRT-1:QIND+STRT+ROWL-2) =
& AA(ROWPOS(I):ROWPOS(I+1)-1)
WORK(QIND+STRT+ROWL-1) = -PP(I)
IWORK(CIND+STRT-1:CIND+STRT+ROWL-2) =
& COLPOS(ROWPOS(I):ROWPOS(I+1)-1)
IWORK(CIND+STRT+ROWL-1) = NN+1
IWORK(RIND+I) = ROWPOS(I+1)+I
END DO
IWORK(RIND+NP1) = ROWPOS(NN+2)+NN
END IF
WORK(QIND+IWORK(RIND+NN)-1) = 1.0
IWORK(CIND+IWORK(RIND+NN)-1) = K
IF (K. LT. NP1) THEN
WORK(QIND+IWORK(RIND+NN)) = 0.0
IWORK(CIND+IWORK(RIND+NN)) = NP1
IWORK(RIND+NP1) = IWORK(RIND+NP1)+1
END IF
CALL ILUFDS(NP1, WORK(QIND:QIND+IWORK(RIND+NP1)-2),
& IWORK(RIND+NP1)-1, IWORK(RIND:RIND+NP1),
& IWORK(CIND:CIND+IWORK(RIND+NP1)-2))
END IF
C
C COMPUTE RIGHT HAND SIDE FOR Mx=b.
C
RHSC(NP1) = STARTK
IF (PRESENT(RHS)) THEN
RHSC(1:NN) = -RHS
ELSE
RHSC(1:NN) = 0.0
END IF
C
C CALL GMRES FOR THE Mx=b SYSTEM.
C
ITMAX=15*NP1
KD=SUBSPACE
CALL GMRES(NP1, KD, ITMAX, RHSC, WORK(ZBIND:ZBIND+NN),
& K, WORK(QIND:QIND+IWORK(RIND+NP1)-2), IFLAG,
& ROWPOSP=IWORK(RIND:RIND+NP1),
& COLPOSP=IWORK(CIND:CIND+IWORK(RIND+NP1)-2))
IF ( IFLAG .GT. 0) RETURN
START(1:NP1) = WORK(ZBIND:ZBIND+NN)
RETURN
END SUBROUTINE GMRILUDS
SUBROUTINE HFUN1P(QDG,LAMBDA,X,
& PDG,CL,COEF,RHO,
& DRHOX,DRHOL,XDGM1,XDG,
& G,DG,PXDGM1,PXDG,
& F,DF,XX,TRM,
& DTRM,CLX,DXNP1,
& N,MAXT,IDEG,
& NUMT,KDEG)
C
C HFUN1P EVALUATES THE CONTINUATION EQUATION "RHO".
C
C NOTE THAT:
C DRHOX IS THE "REALIFICATION" OF DCRHOX, WHERE
C DCRHOX DENOTES THE (COMPLEX) PARTIAL
C DERIVATIVE MATRIX OF THE CONTINUATION SYSTEM
C WITH RESPECT TO X, AND
C DRHOL IS THE "REALIFICATION" OF DCRHOL, WHERE
C DCRHOL DENOTES THE (COMPLEX) PARTIAL
C DERIVATIVE MATRIX OF THE CONTINUATION SYSTEM
C WITH RESPECT TO LAMBDA. THUS
C DRHOX(2J-1,2K-1) = DCRHOX(1,J,K)
C DRHOX(2J ,2K ) = DCRHOX(1,J,K)
C DRHOX(2J-1,2K ) =-DCRHOX(2,J,K)
C DRHOX(2J ,2K-1) = DCRHOX(2,J,K)
C DRHOL(2J-1,N2P1) = DCRHOL(1,J)
C DRHOL(2J ,N2P1) = DCRHOL(2,J)
C RHO(2J-1) = CRHO(1,J)
C RHO(2J ) = CRHO(2,J)
C WHERE CRHO DENOTES THE (COMPLEX) CONTINUATION SYSTEM,
C THE INITIAL "1" OR "2" DENOTES REAL OR IMAGINARY PARTS,
C RESPECTIVELY, "J" INDEXES THE EQUATION, "K" INDEXES THE PARTIAL
C DERIVATIVE, AND NEITHER DCRHOX NOR DCRHOL ARE PROGRAM VARIABLES.
C
C ON INPUT:
C
C QDG IS THE "RANDOM" PARAMETER "A".
C
C LAMBDA IS THE CONTINUATION PARAMETER.
C
C X IS THE INDEPENDENT VARIABLE.
C
C PDG IS ONE OF THE PARAMETERS THAT DEFINES G (SEE SUBROUTINE
C GFUNP).
C
C CL IS ONE OF THE PARAMETERS THAT DEFINES F (SEE SUBROUTINE
C FFUNP).
C
C COEF IS ONE OF THE PARAMETERS THAT DEFINES F (SEE SUBROUTINE
C FFUNP).
C
C ON OUTPUT:
C
C RHO IS THE HOMOTOPY.
C
C DRHOX CONTAINS THE PARTIAL DERIVATIVES OF RHO WITH RESPECT
C TO X.
C
C DRHOL CONTAINS THE PARTIAL DERIVATIVES OF RHO WITH RESPECT
C TO LAMBDA.
C
C THE FOLLOWING ARE VARIABLES WHOSE WORKSPACE IS PASSED FROM HFUNP:
C XDGM1
C XDG
C G
C DG
C PXDGM1
C PXDG
C F
C DF
C XX
C TRM
C DTRM
C CLX
C DXNP1
C N
C MAXT
C IDEG
C NUMT
C KDEG
C
C OTHER VARIABLES:
C ONEML
C
C SUBROUTINES: GFUNP, FFUNP.
C
USE REAL_PRECISION
C DECLARATION OF INPUT, WORKSPACE, AND OUTPUT:
INTEGER, INTENT(IN):: N,MAXT,IDEG(N),NUMT(N),KDEG(N,N+1,MAXT)
REAL (KIND=R8), INTENT(IN):: QDG(2,N),LAMBDA,X(2,N),
& PDG(2,N),CL(2,N+1),COEF(N,MAXT)
REAL (KIND=R8), INTENT(OUT):: RHO(2*N),DRHOX(2*N,2*N),DRHOL(2*N)
REAL (KIND=R8), INTENT(IN OUT):: XDGM1(2,N),XDG(2,N),
& G(2,N),DG(2,N),PXDGM1(2,N),PXDG(2,N),
& F(2,N), DF(2,N,N+1),XX(2,N,N+1,MAXT),TRM(2,N,MAXT),
& DTRM(2,N,N+1,MAXT),CLX(2,N),DXNP1(2,N)
C
C DECLARATION OF LOCAL VARIABLES:
INTEGER:: J,J2,J2M1,K,K2,K2M1
REAL (KIND=R8):: ONEML
C
CALL GFUNP(N,IDEG,PDG,QDG,X,XDGM1,XDG,PXDGM1,PXDG,G,DG)
CALL FFUNP(N,NUMT,MAXT,KDEG,COEF,CL,X,XX,TRM,DTRM,CLX,DXNP1,F,DF)
ONEML=1.0 - LAMBDA
DO J=1,N
J2=2*J
J2M1=J2-1
DO K=1,N
K2=2*K
K2M1=K2-1
DRHOX(J2M1,K2M1)= LAMBDA*DF(1,J,K)
DRHOX(J2 ,K2 )= DRHOX(J2M1,K2M1)
DRHOX(J2 ,K2M1)= LAMBDA*DF(2,J,K)
DRHOX(J2M1,K2 )=-DRHOX(J2 ,K2M1)
END DO
DRHOX(J2M1,J2M1)= DRHOX(J2M1,J2M1) + ONEML*DG(1,J)
DRHOX(J2 ,J2 )= DRHOX(J2M1,J2M1)
DRHOX(J2 ,J2M1)= DRHOX(J2 ,J2M1) + ONEML*DG(2,J)
DRHOX(J2M1,J2 )=-DRHOX(J2 ,J2M1)
DRHOL(J2M1) = F(1,J) - G(1,J)
DRHOL(J2) = F(2,J) - G(2,J)
RHO(J2M1) = LAMBDA*F(1,J) + ONEML* G(1,J)
RHO(J2 ) = LAMBDA*F(2,J) + ONEML* G(2,J)
END DO
RETURN
END SUBROUTINE HFUN1P
SUBROUTINE HFUNP(N,QDG,LAMBDA,X)
C
C HFUNP ALLOCATES STORAGE FOR SUBROUTINE HFUN1P FROM THE WORK ARRAYS
C PAR AND IPAR, AS FOLLOWS:
C
C DOUBLE PRECISION VARIABLES AND ARRAYS PASSED IN PAR
C
C PAR INDEX VARIABLE NAME LENGTH
C ---------- ------------- -----------------
C 1 PDG 2*N
C 2 CL 2*(N+1)
C 3 COEF N*MAXT
C 4 RHO N2
C 5 DRHOX N2*N2
C 6 DRHOL N2
C 7 XDGM1 2*N
C 8 XDG 2*N
C 9 G 2*N
C 10 DG 2*N
C 11 PXDGM1 2*N
C 12 PXDG 2*N
C 13 F 2*N
C 14 DF 2*N*(N+1)
C 15 XX 2*N*(N+1)*MAXT
C 16 TRM 2*N*MAXT
C 17 DTRM 2*N*(N+1)*MAXT
C 18 CLX 2*N
C 19 DXNP1 2*N
C
C INTEGER VARIABLES AND ARRAYS PASSED IN IPAR
C
C IPAR INDEX VARIABLE NAME LENGTH OFFSET
C ---------- ------------- -----------------
C 1 N 1 1
C 2 MAXT 1 2
C 3 PROFF 25 3
C 4 IPROFF 15 28
C 5 IDEG N 43
C 6 NUMT N 43+N
C 7 KDEG N*(N+1)*MAXT 43+N2
C
C ON INPUT:
C
C N IS THE NUMBER OF EQUATIONS AND VARIABLES.
C
C QDG IS THE "RANDOM" VECTOR DENOTED "A" IN HOMPACK DOCUMENTATION.
C
C LAMBDA IS THE CONTINUATION PARAMETER.
C
C X IS THE INDEPENDENT VARIABLE.
C
C ON OUTPUT:
C
C THE GLOBAL WORK ARRAYS PAR AND IPAR HAVE BEEN UPDATED.
C
C SUBROUTINES: HFUN1P.
C
USE HOMOTOPY
USE REAL_PRECISION
INTEGER, INTENT(IN):: N
REAL (KIND=R8), INTENT(IN):: QDG(2,N),LAMBDA,X(2,N)
C
CALL HFUN1P(QDG,LAMBDA,X,
& PAR( IPAR(3 + ( 1-1))), PAR( IPAR(3 + ( 2-1))),
& PAR( IPAR(3 + ( 3-1))), PAR( IPAR(3 + ( 4-1))),
& PAR( IPAR(3 + ( 5-1))), PAR( IPAR(3 + ( 6-1))),
& PAR( IPAR(3 + ( 7-1))), PAR( IPAR(3 + ( 8-1))),
& PAR( IPAR(3 + ( 9-1))), PAR( IPAR(3 + (10-1))),
& PAR( IPAR(3 + (11-1))), PAR( IPAR(3 + (12-1))),
& PAR( IPAR(3 + (13-1))), PAR( IPAR(3 + (14-1))),
& PAR( IPAR(3 + (15-1))), PAR( IPAR(3 + (16-1))),
& PAR( IPAR(3 + (17-1))), PAR( IPAR(3 + (18-1))),
& PAR( IPAR(3 + (19-1))),
&IPAR( IPAR(28+ ( 1-1))),IPAR( IPAR(28+ ( 2-1))),
&IPAR( IPAR(28+ ( 5-1))),IPAR( IPAR(28+ ( 6-1))),
&IPAR( IPAR(28+ ( 7-1))) )
C
RETURN
END SUBROUTINE HFUNP
SUBROUTINE ILUFDS(NP1, B, LENB, ROWPOSP, COLPOSP)
c
C Computes the incomplete LU factorization of the matrix B,
C where B is NP1 x NP1. B is assumed to be stored in the general
C sparse row scheme.
C
C The method used is that found in TR 89-41, Department of Computer
C Science, VPI&SU, Blacksburg, VA, 1989: 'Preconditioned
C conjugate gradient algorithms for homotopy curve tracking',
C page 10.
C---------------------------------------------------------------------
C
C Input variables:
C B matrix to be factored.
C ROWPOSP indices of row-starts within B.
C COLPOSP column indices for matrix B stored in the general
C sparse row scheme.
C LENB number of entries in B.
C NP1 the dimension of B.
C
C Output variables:
C B the ILU factors of input matrix B.
C-------------------------------------------------------------
C
USE REAL_PRECISION
INTEGER, INTENT(IN):: LENB, NP1, ROWPOSP(NP1+1), COLPOSP(LENB)
REAL (KIND=R8), INTENT(IN OUT):: B(LENB)
C
C LOCAL VARIABLES
REAL (KIND=R8):: SIJ, LIT, LII
INTEGER:: I, J, COUNT, ISTRT, IFIN, TMAX, K, T, M
C
C
DO I = 1, NP1
ISTRT = ROWPOSP(I)
IFIN = ROWPOSP(I+1) - 1
C--------------------------------------- FOR EACH ELEMENT IN ROW I,
C COMPUTE THE COLUMN NUMBER
C AND T.
DO COUNT = ISTRT, IFIN
J = COLPOSP(COUNT)
TMAX = MIN0(I,J) - 1
SIJ = B(COUNT)
C--------------------------------------- COMPUTE THE CORRESPONDING
C SUM OF PRODUCTS OF ELEMENTS
C OF L AND U.
K = ISTRT
42 T = COLPOSP(K)
IF (T .LE. TMAX) THEN
LIT = B(K)
M = ROWPOSP(T)
C--------------------------------------- FIND VALUE OF U_{TJ}.
20 IF (COLPOSP(M) .LT. J) THEN
M = M + 1
GOTO 20
ENDIF
IF (COLPOSP(M) .EQ. J) SIJ = SIJ - LIT*B(M)
K = K + 1
GOTO 42
ENDIF
C--------------------------------------- END OF 'T' LOOP.
K = ISTRT
C--------------------------------------- FIND VALUE OF L_{II}.
30 IF (COLPOSP(K) .LT. I) THEN
K = K+1
GOTO 30
ENDIF
IF (COLPOSP(K) .EQ. I) THEN
LII = B(K)
IF (DABS(LII) .EQ. 0.0) THEN
LII = 0.00001
B(K) = 0.00001
ENDIF
ELSE
LII = 0.00001
ENDIF
C--------------------------------------- UPDATE L OR U, AS NEEDED.
IF (I .GE. J) THEN
B(COUNT) = SIJ
ELSE
B(COUNT) = SIJ/LII
ENDIF
END DO
C--------------------------------------- END OF 'COUNT' LOOP.
END DO
C
RETURN
END SUBROUTINE ILUFDS
SUBROUTINE ILUSOLVDS(NN, Q, LENQ, ROWPOSP, COLPOSP, B)
C
C Computes Q^{-1}*B -- returns result as B.
C---------------------------------------------------------------------
C
C Input variables:
C
C Q triangular factors of preconditioning matrix, stored
C in the general sparse row scheme.
C ROWPOSP indices of row-starts within B.
C COLPOSP column indices for matrix B stored in the general
C sparse row scheme.
C NN logical row dimension of Q.
C LENQ number of data entries in Q.
C B right hand side -- should have dimension NN.
C
C Output variables:
C
C B solution of Qx = B.
C
C---------------------------------------------------------------------
C
USE REAL_PRECISION
INTEGER, INTENT(IN) :: NN, LENQ, ROWPOSP(NN+1), COLPOSP(LENQ)
REAL (KIND=R8), INTENT(IN) :: Q(LENQ)
REAL (KIND=R8), INTENT(IN OUT) :: B(NN)
C
C LOCAL VARIABLES
INTEGER:: DIAG(NN), I, K, J
C------------------------------------------------ COMPUTE B = INV(L)*B
B(1) = B(1)/Q(1)
DIAG(1) = 1
DO I = 2, NN
K = ROWPOSP(I)
42 J = COLPOSP(K)
IF (J .LT. I) THEN
B(I) = B(I) - Q(K)*B(J)
K = K + 1
GOTO 42
ELSE
DIAG(I) = K
B(I) = B(I)/Q(K)
ENDIF
END DO
C------------------------------------------------ COMPUTE B = INV(U)*B
DO I = NN-1, 1, -1
DO K = DIAG(I)+1, ROWPOSP(I+1)-1
B(I) = B(I) - Q(K)*B(COLPOSP(K))
END DO
END DO
RETURN
END SUBROUTINE ILUSOLVDS
SUBROUTINE INITP(IFLG1,N,NUMT,KDEG,COEF,
& IDEG,FACV,CL,PDG,QDG,R)
C
C INITP INITIALIZES THE CONSTANTS THAT DEFINE THE POLSYS HOMOTOPY,
C INITIALIZES THE CONSTANTS THAT DEFINE THE PROJECTIVE TRANSFORMATION,
C AND SCALES THE COEFFICIENTS (IF SCALING IS SPECIFIED).
C
C ON INPUT:
C
C IFLG1 IS A FLAG THAT SPECIFIES WHETHER THE COEFFICIENTS ARE TO
C BE SCALED OR NOT AND WHETHER THE PROJECTIVE TRANSFORMATION IS TO
C BE USED OR NOT. IFLG1=A*10+B. SCALING IS SPECIFIED WHEN B=1. THE
C PROJECTIVE TRANSFORMATION IS SPECIFIED WHEN A=1. OTHERWISE, A AND/OR
C B =0. SCALING IS EVOKED BY A CALL TO THE SUBROUTINE SCLGNP. THE
C PROJECTIVE TRANSFORMATION IS EVOKED BY SETTING THE CL ARRAY EQUAL
C TO RANDOM COMPLEX NUMBERS. OTHERWISE, CL IS SET TO NOMINAL VALUES.
C
C N IS THE NUMBER OF EQUATIONS AND VARIABLES.
C
C NUMT(J) IS THE NUMBER OF TERMS IN EQUATION J, FOR J=1 TO N.
C
C KDEG(J,L,K) IS THE DEGREE OF THE L-TH VARIABLE, X(L), IN THE K-TH
C TERM OF THE J-TH EQUATION, WHERE J=1 TO N, L=1 TO N+1, AND K=1 TO
C NUMT(J). THE CASE "L=N+1" IS SPECIAL, AND KDEG IS NOT AN INPUT
C VALUE TO POLSYS , BUT RATHER IS COMPUTED IN THIS SUBROUTINE.
C
C COEF(J,K) IS THE COEFFICIENT OF THE K-TH TERM FOR THE J-TH
C EQUATION, WHERE J=1 TO N AND K=1 TO NUMT(J).
C
C
C ON OUTPUT:
C
C IDEG(J) IS THE DEGREE OF THE J-TH EQUATION FOR J=1 TO N.
C
C FACV(J) IS THE SCALE FACTOR FOR THE J-TH VARIABLE.
C
C CL(2,1:N+1) IS AN ARRAY USED TO DEFINE THE PROJECTIVE
C TRANSFORMATION. IT IS USED IN SUBROUTINES FFUNP AND OTPUTP
C TO DEFINE THE PROJECTIVE COORDINATE, XNP1.
C
C PDG IS USED IN SUBROUTINE GFUNP TO DEFINE THE INITIAL SYSTEM,
C G(X)=0.
C
C QDG IS USED IN SUBROUTINE GFUNP TO DEFINE THE INITIAL SYSTEM,
C G(X)=0.
C
C R IS USED IN SUBROUTINE STRPTP TO GENERATE SOLUTIONS TO G(X)=0.
C
C
USE HOMOTOPY
USE REAL_PRECISION
C DECLARATIONS OF INPUT AND OUTPUT:
INTEGER, INTENT(IN):: IFLG1,N,NUMT(:)
INTEGER, INTENT(IN OUT):: KDEG(:,:,:)
REAL (KIND=R8), INTENT(IN OUT):: COEF(:,:)
INTEGER, INTENT(OUT):: IDEG(N)
REAL (KIND=R8), INTENT(OUT)::
& FACV(N),CL(2,N+1),PDG(2,N),QDG(2,N),R(2,N)
C
C DECLARATIONS OF LOCAL VARIABLES:
INTEGER:: IERR,J,JJ,MAXT,N2,NP1
REAL (KIND=R8):: CCL(2,11),P(2,10),Q(2,10),ZERO
C
INTERFACE
SUBROUTINE DIVP(XXXX,YYYY,ZZZZ,IERR)
USE REAL_PRECISION
REAL (KIND=R8), DIMENSION(2), INTENT(IN):: XXXX,YYYY
INTEGER, INTENT(OUT):: IERR
REAL (KIND=R8), DIMENSION(2), INTENT(OUT):: ZZZZ
END SUBROUTINE DIVP
SUBROUTINE MULP(XXXX,YYYY,ZZZZ)
USE REAL_PRECISION
REAL (KIND=R8), DIMENSION(2), INTENT(IN):: XXXX,YYYY
REAL (KIND=R8), DIMENSION(2), INTENT(OUT):: ZZZZ
END SUBROUTINE MULP
SUBROUTINE POWP(NNNN,XXXX,YYYY)
USE REAL_PRECISION
INTEGER, INTENT(IN):: NNNN
REAL (KIND=R8), DIMENSION(2), INTENT(IN):: XXXX
REAL (KIND=R8), DIMENSION(2), INTENT(IN OUT):: YYYY
END SUBROUTINE POWP
SUBROUTINE SCLGNP(N,MAXT,NUMT,DEG,MODE,EPS0,COEF,
& NNUMT,DDEG,CCOEF,ALPHA,BETA,RWORK,XWORK,
& FACV,FACE,COESCL,IERR)
USE REAL_PRECISION
INTEGER, INTENT(IN):: N,MAXT,NUMT(:),DEG(:,:,:),MODE
REAL (KIND=R8), INTENT(IN):: EPS0,COEF(:,:)
INTEGER, INTENT(IN OUT):: NNUMT(N),DDEG(N,N+1,MAXT)
REAL (KIND=R8), INTENT(IN OUT):: CCOEF(N,MAXT),ALPHA(2*N,2*N),
& BETA(2*N),RWORK(N*(2*N+1)),XWORK(2*N)
REAL (KIND=R8), INTENT(OUT):: FACV(N),FACE(N),COESCL(N,MAXT)
INTEGER, INTENT(OUT):: IERR
END SUBROUTINE SCLGNP