SUBROUTINE POLSYS(N,NUMT,COEF,KDEG,IFLG1,IFLG2,EPSBIG,EPSSML,
     $     SSPAR,NUMRR,NN,MMAXT,TTOTDG,LENWK,LENIWK,
     $     LAMBDA,ROOTS,ARCLEN,NFE,WK,IWK)
C
C POLSYS FINDS ALL (COMPLEX) SOLUTIONS TO A SYSTEM
C F(X)=0 OF N POLYNOMIAL EQUATIONS IN N UNKNOWNS
C WITH REAL COEFFICIENTS. IF IFLG=10 OR IFLG=11, POLSYS
C RETURNS THE SOLUTIONS AT INFINITY ALSO.
C
C THE SYSTEM F(X)=0 IS DESCRIBED VIA THE COEFFICENTS,
C "COEF", AND THE PARAMETERS "N, NUMT, KDEG", AS FOLLOWS.
C
C
C       NUMT(J)
C
C F(J) = SUM  COEF(J,K) * X(1)**KDEG(J,1,K)...X(N)**KDEG(J,N,K)
C
C        K=1
C
C FOR J=1, ..., N.
C
C
C POLSYS HAS TWO MAIN RUN OPTIONS:  AUTOMATIC SCALING AND
C THE PROJECTIVE TRANSFORMATION.  THESE ARE EVOKED VIA THE
C FLAG "IFLG1", AS DESCRIBED BELOW.  THE OTHER INPUT
C PARAMETERS ARE THE SAME WHETHER ONE OR BOTH OF THESE OPTIONS
C ARE SPECIFIED, AND THE OUTPUT IS ALWAYS RETURNED UNSCALED
C AND UNTRANSFORMED.
C
C IF AUTOMATIC SCALING IS SPECIFIED, THEN THE INPUT
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 M-TH PATH,
C
C ROOTS(2,J,M) = IMAGINARY PART OF X(J) FOR THE M-TH PATH,
C
C FOR J=1, ..., N, AND
C
C ROOTS(1,N+1,M) = REAL PART OF Z(N+1) FOR THE M-TH PATH,
C
C ROOTS(2,N+1,M) = IMAGINARY PART OF Z(N+1) FOR THE M-TH 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 POLSYS.
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 THE 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  POLSYS  AND WHETHER THE PROJECTIVE TRANSFORMATION
C OF THE SYSTEM IS TO BE AUTOMATICALLY EVOKED BY POLSYS.  SEE
C "ON INPUT" BELOW.
C
C IFLG2, EPSBIG, EPSSML, AND  SSPAR  TELL THE PATH TRACKER
C POLYNF  WHICH PATHS TO TRACK AND SET PARAMETERS FOR THE PATH
C TRACKER.
C
C NUMRR  TELLS  POLSYS  HOW MANY MULTIPLES OF 1000 STEPS TO TRY
C BEFORE ABANDONING A PATH.
C
C NN, MMAXT, TTOTDG, LENWK, LENIWK  GIVE THE DIMENSIONS OF ARRAYS.
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 THE FOLLOWING SUBROUTINES ARE USED DIRECTLY OR INDIRECTLY BY
C POLSYS:
C         SPECIAL FOR POLSYS:
C           POLYP , INITP , STRPTP ,
C           OTPUTP , RHO , RHOJAC ,
C           HFUNP , HFUN1P , GFUNP , FFUNP ,
C           MULP , POWP , DIVP,
C           SCLGNP.
C         FROM THE GENERAL HOMPACK ROUTINES:
C           POLYNF , STEPNF , TANGNF , ROOTNF , ROOT ,
C           QRFAQF, QRSLQF , D1MACH , DDOT , DNRM2.
C
C ON INPUT:
C
C N  IS THE NUMBER OF EQUATIONS AND VARIABLES.
C
C NUMT(1:NN)  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:NN,1:MMAXT)  IS A REAL ARRAY.  COEF(J,K)  IS
C   THE K-TH COEFFICIENT OF THE J-TH EQUATION FOR J=1 TO N,
C   K=1 TO NUMT(J).
C
C KDEG(1:NN,1:NN+1,1:MMAXT)  IS AN INTEGER ARRAY.
C   KDEG(J,L,K)  IS THE DEGREE OF THE L-TH VARIABLE IN THE K-TH
C   TERM OF THE J-TH EQUATION FOR  J=1 TO N, L=1 TO N, K=1 TO NUMT(J).
C
C IFLG1 =
C   00  IF THE PROBLEM IS TO BE SOLVED WITHOUT
C       CALLING POLSYS' 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:TTOTDG)  IS AN INTEGER ARRAY.  IF IFLG2(M) = -2, THEN THE
C   M-TH PATH IS TRACKED.  OTHERWISE THE M-TH 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 IFLG(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  POLYNF ) ARE SET TO  EPSBIG.
C
C EPSSML  IS THE ACCURACY DESIRED FOR THE FINAL SOLUTION.  ANSRE AND
C   ANSAE (IN  POLYNF ) 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  POLYNF .  OTHERWISE THE INPUT VALUE OF  SSPAR(J)  IS USED.
C    SEE THE COMMENTS IN  POLYNF  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 NN  IS THE DECLARED DIMENSION OF  NUMT  AND OF THE
C   FIRST INDEX OF  COEF  AND  KDEG .  THE SECOND INDEX OF
C   KDEG  AND  ROOTS  IS DIMENSIONED NN+1.  NN  MUST
C   BE GREATER THAN OR EQUAL TO N.
C
C MMAXT  IS THE DECLARED DIMENSION OF THE SECOND INDEX OF
C   COEF  AND THE THIRD INDEX OF  KDEG.  MMAXT  MUST BE
C   GREATER THAN OR EQUAL TO THE MAXIMUM NUMBER OF
C   TERMS IN EACH EQUATION.  IN OTHER WORDS,
C   MMAXT .GE. MAX {NUMT(J) | J=1, ..., N} .
C
C TTOTDG  IS THE DECLARED DIMENSION OF  IFLG2 , LAMBDA , ARCLEN ,
C   NFE , AND OF THE THIRD INDEX OF  ROOTS.  TTOTDG
C   MUST BE GREATER THAN OR EQUAL TO TOTDG, THE TOTAL
C   DEGREE OF THE SYSTEM.
C
C LENWK  IS THE DIMENSION OF THE WORKSPACE  WK .  LENWK  MUST
C   BE GREATER THAN OR EQUAL TO
C     21 + 61*N + 10*N**2 + 7*N*MMAXT + 4*N**2*MMAXT.
C
C LENIWK  IS THE DIMENSION OF THE WORKSPACE  IWK .  LENIWK  MUST BE
C   GREATER THAN OR EQUAL TO  43 + 7*N + N*(N+1)*MMAXT.
C
C
C ON OUTPUT:
C
C N, NUMT, COEF, KDEG, NN, MMAXT, TTOTDG, LENWK, AND  LENIWK
C   ARE UNCHANGED.
C
C IFLG1=
C   -1  IF  NN  IS TOO SMALL.
C   -2  IF  MMAXT  IS TOO SMALL.
C   -3  IF  TTOTDG  IS TOO SMALL.
C   -4  IF  LENWK  IS TOO SMALL.
C   -5  IF  LENIWK  IS TOO SMALL.
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 M-TH 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 TOLERANCE  EPSBIG  MAY BE TOO
C       STRINGENT.
C
C   7   ILLEGAL INPUT PARAMETERS, A FATAL ERROR.
C
C LAMBDA(M)  IS THE FINAL LAMBDA VALUE FOR THE M-TH 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 M-TH 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 M-TH PATH FOR M = 1, ..., TOTDG.
C
C NFE(M)  IS THE NUMBER OF JACOBIAN MATRIX EVALUATIONS REQUIRED TO
C   TRACK THE M-TH PATH FOR M =1, ..., TOTDG.
C
C ----------------------------------------------------------------------
C TYPE DECLARATIONS FOR INPUT AND OUTPUT
      INTEGER N,NUMT,KDEG,IFLG1,IFLG2,NUMRR,NN,MMAXT,
     $  TTOTDG,LENWK,LENIWK,NFE,IWK
      DOUBLE PRECISION COEF,EPSBIG,EPSSML,SSPAR,LAMBDA,ROOTS,
     $  ARCLEN,WK
C
C ARRAY DECLARATIONS FOR INPUT AND OUTPUT
      DIMENSION NUMT(NN),KDEG(NN,NN+1,MMAXT),IFLG2(TTOTDG),
     $  NFE(TTOTDG),IWK(LENIWK)
      DIMENSION COEF(NN,MMAXT),SSPAR(8),LAMBDA(TTOTDG),
     $  ROOTS(2,NN+1,TTOTDG), ARCLEN(TTOTDG),WK(LENWK)
C
C TYPE DECLARATIONS FOR VARIABLES
      INTEGER I,IDEG,IIDEG,IWKOFF,J,K,L,LENIWW,LENWKK,LIWK,LWK,
     $  MAXT,N2,TOTDG,WKOFF
C
C ARRAY DECLARATIONS FOR VARIABLES
      DIMENSION LWK(19),LIWK(4),WKOFF(19),IWKOFF(4)
C
C CHECK THAT BASIC DIMENSIONAL PARAMETERS ARE BIG ENOUGH
C
      IF(NN.LT.N) THEN
        IFLG1=-1
        RETURN
      END IF
      MAXT=0
      DO 50 J=1,N
        IF(MAXT.LT.NUMT(J))MAXT=NUMT(J)
 50   CONTINUE
      IF(MMAXT.LT.MAXT) THEN
        IFLG1=-2
        RETURN
      END IF
      TOTDG=1
      DO 80 J=1,N
         IDEG=0
         DO 70 K=1,NUMT(J)
             IIDEG=0
             DO 60 L=1,N
                IIDEG=IIDEG+KDEG(J,L,K)
 60          CONTINUE
             IF(IIDEG.GT.IDEG)IDEG=IIDEG
 70      CONTINUE
         TOTDG=TOTDG*IDEG
 80   CONTINUE
      IF(TTOTDG.LT.TOTDG) THEN
        IFLG1=-3
        RETURN
      END IF
      LENWKK = 21 + 61*N + 10*N**2 + 7*N*MMAXT + 4*N**2*MMAXT
      IF(LENWK.LT.LENWKK) THEN
        IFLG1=-4
        RETURN
      END IF
      LENIWW = 43 +  7*N + N*(N+1)*MMAXT
      IF(LENIWK.LT.LENIWW) THEN
        IFLG1=-5
        RETURN
      END IF
      IF(IFLG1.NE.0.AND.IFLG1.NE.1.AND.IFLG1.NE.10.AND.IFLG1.NE.11) THEN
        IFLG1=-6
        RETURN
      END IF
C
C VARIABLES THAT ARE PASSED IN ARRAY WK: (LENGTHS ARE IN THE
C INTEGER ARRAY LWK.)
C
C    VARIABLE NAME       LENGTH
C
C    1   PDG               N2
C    2   QDG               N2
C    3   R                 N2
C    4   FACV              N
C    5   CL                2*(N+1)
C    6   Y                 N2+1
C    7   YP                N2+1
C    8   YOLD              N2+1
C    9   YPOLD             N2+1
C   10   QR                N2*(N2+2)
C   11   ALPHA             N2
C   12   TZ                N2+1
C   13   W                 N2+1
C   14   WP                N2+1
C   15   Z0                N2+1
C   16   Z1                N2+1
C   17   SSPAR             8
C   18   PAR               2 + 28*N + 6*N**2 + 7*N*MMAXT + 4*N**2*MMAXT
C
C VARIABLES THAT ARE PASSED IN ARRAY IWK: (LENGTHS ARE IN THE
C INTEGER ARRAY LIWK.)
C
C    VARIABLE NAME       LENGTH
C
C    1   IDEG              N
C    2   ICOUNT            N
C    3   PIVOT             N2+1
C    4   IPAR              42 + 2*N + N*(N+1)*MMAXT
C
      N2=2*N
      LWK(1)= N2
      LWK(2)= N2
      LWK(3)= N2
      LWK(4)= N
      LWK(5)= 2*(N+1)
      LWK(6)= N2+1
      LWK(7)= N2+1
      LWK(8)= N2+1
      LWK(9)= N2+1
      LWK(10)=N2*(N2+2)
      LWK(11)=N2
      LWK(12)=N2+1
      LWK(13)=N2+1
      LWK(14)=N2+1
      LWK(15)=N2+1
      LWK(16)=N2+1
      LWK(17)=8
      LWK(18)= 2 + 28*N + 6* N**2 + 7*N*MMAXT + 4* N**2 *MMAXT
C
      LIWK(1)=N
      LIWK(2)=N
      LIWK(3)=2*N+1
      LIWK(4)= 42 + 2*N + N*(N+1)*MMAXT
C
C WKOFF AND IWKOFF ARE OFFSETS THAT DEFINE THE VARIABLES LISTED ABOVE
C
      WKOFF(1)=1
      DO 100 I=2,18
          WKOFF(I)=WKOFF(I-1)+LWK(I-1)
 100  CONTINUE
      IWKOFF(1)=1
      DO 200 I=2,4
          IWKOFF(I)=IWKOFF(I-1)+LIWK(I-1)
 200  CONTINUE
      DO 300 J=1,8
        WK(WKOFF(17) + (J-1))=SSPAR(J)
 300  CONTINUE
C
      CALL POLYP(N,NUMT,COEF,KDEG,IFLG1,IFLG2,EPSBIG,EPSSML,
     $ NUMRR,NN,MMAXT,TTOTDG,LAMBDA,ROOTS,ARCLEN,NFE,TOTDG,
     $ WK( WKOFF( 1)),WK( WKOFF( 2)),WK( WKOFF( 3)),WK( WKOFF( 4)),
     $ WK( WKOFF( 5)),WK( WKOFF( 6)),WK( WKOFF( 7)),WK( WKOFF( 8)),
     $ WK( WKOFF( 9)),WK( WKOFF(10)),WK( WKOFF(11)),WK( WKOFF(12)),
     $ WK( WKOFF(13)),WK( WKOFF(14)),WK( WKOFF(15)),WK( WKOFF(16)),
     $ WK( WKOFF(17)),WK( WKOFF(18)),
     $IWK(IWKOFF( 1)),IWK(IWKOFF( 2)),IWK(IWKOFF( 3)),IWK(IWKOFF( 4)))
C
      RETURN
      END

.