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