share/hompack/fortran/polsys.f

   1       SUBROUTINE POLSYS(N,NUMT,COEF,KDEG,IFLG1,IFLG2,EPSBIG,EPSSML,
   2      $     SSPAR,NUMRR,NN,MMAXT,TTOTDG,LENWK,LENIWK,
   3      $     LAMBDA,ROOTS,ARCLEN,NFE,WK,IWK)
   4 C
   5 C POLSYS FINDS ALL (COMPLEX) SOLUTIONS TO A SYSTEM
   6 C F(X)=0 OF N POLYNOMIAL EQUATIONS IN N UNKNOWNS
   7 C WITH REAL COEFFICIENTS. IF IFLG=10 OR IFLG=11, POLSYS
   8 C RETURNS THE SOLUTIONS AT INFINITY ALSO.
   9 C
  10 C THE SYSTEM F(X)=0 IS DESCRIBED VIA THE COEFFICENTS,
  11 C "COEF", AND THE PARAMETERS "N, NUMT, KDEG", AS FOLLOWS.
  12 C
  13 C
  14 C       NUMT(J)
  15 C
  16 C F(J) = SUM  COEF(J,K) * X(1)**KDEG(J,1,K)...X(N)**KDEG(J,N,K)
  17 C
  18 C        K=1
  19 C
  20 C FOR J=1, ..., N.
  21 C
  22 C
  23 C POLSYS HAS TWO MAIN RUN OPTIONS:  AUTOMATIC SCALING AND
  24 C THE PROJECTIVE TRANSFORMATION.  THESE ARE EVOKED VIA THE
  25 C FLAG "IFLG1", AS DESCRIBED BELOW.  THE OTHER INPUT
  26 C PARAMETERS ARE THE SAME WHETHER ONE OR BOTH OF THESE OPTIONS
  27 C ARE SPECIFIED, AND THE OUTPUT IS ALWAYS RETURNED UNSCALED
  28 C AND UNTRANSFORMED.
  29 C
  30 C IF AUTOMATIC SCALING IS SPECIFIED, THEN THE INPUT
  31 C COEFFICIENTS ARE MODIFIED BY SUBROUTINE  SCLGNP . THE PROBLEM
  32 C IS SOLVED WITH THE SCALED COEFFICIENTS AND SCALED VARIABLES.
  33 C THE COEFFICIENTS ARE RETURNED SCALED.
  34 C
  35 C IF THE PROJECTIVE TRANSFORMATION IS SPECIFIED, THEN
  36 C ESSENTIALLY THE SYSTEM IS REFORMULATED IN HOMOGENEOUS
  37 C COORDINATES, Z(1), ..., Z(N+1), AND SOLVED IN COMPLEX
  38 C PROJECTIVE SPACE.  THE RESULTING SOLUTIONS ARE
  39 C UNTRANSFORMED VIA
  40 C
  41 C X(J) = Z(J)/Z(N+1)   J=1, ..., N.
  42 C
  43 C ON RETURN,
  44 C
  45 C ROOTS(1,J,M) = REAL PART OF X(J) FOR THE M-TH PATH,
  46 C
  47 C ROOTS(2,J,M) = IMAGINARY PART OF X(J) FOR THE M-TH PATH,
  48 C
  49 C FOR J=1, ..., N, AND
  50 C
  51 C ROOTS(1,N+1,M) = REAL PART OF Z(N+1) FOR THE M-TH PATH,
  52 C
  53 C ROOTS(2,N+1,M) = IMAGINARY PART OF Z(N+1) FOR THE M-TH PATH.
  54 C
  55 C IF ROOTS(*,N+1,M) IS SMALL, THEN THE ASSOCIATED SOLUTION
  56 C SHOULD BE REGARDED AS BEING "NEAR INFINITY".  NOTE THAT,
  57 C WHEN THE PROJECTIVE TRANSFORMATION HAS BEEN SPECIFIED, THE
  58 C ROOTS VALUES HAVE BEEN UNTRANSFORMED -- THAT IS, DIVIDED
  59 C THROUGH BY Z(N+1) -- UNLESS SUCH DIVISION WOULD HAVE CAUSED
  60 C OVERFLOW.  IN THIS LATTER CASE, THE AFFECTED COMPONENTS OF
  61 C ROOTS ARE SET TO THE LARGEST FLOATING POINT NUMBER (MACHINE
  62 C INFINITY).
  63 C
  64 C THE CODE CAN BE MODIFIED EASILY TO SOLVE SYSTEMS WITH COMPLEX
  65 C COEFFICIENTS,  COEF .  ONLY THE SUBROUTINES  INITP  AND  FFUNP
  66 C NEED BE CHANGED.
  67 C
  68 C THE FORTRAN COMPLEX DECLARATION IS NOT USED IN POLSYS.
  69 C COMPLEX VARIABLES ARE REPRESENTED BY REAL ARRAYS WITH FIRST
  70 C INDEX DIMENSIONED 2 AND COMPLEX OPERATIONS ARE EVOKED BY
  71 C SUBROUTINE CALLS.
  72 C
  73 C THE TOTAL NUMBER OF PATHS THAT WILL THE TRACKED (IF
  74 C IFLG2(M)=-2 FOR ALL M) IS EQUAL TO THE "TOTAL DEGREE" OF THE
  75 C SYSTEM, TOTDG.   TOTDG IS EQUAL TO THE PRODUCTS OF THE
  76 C DEGREES OF ALL THE EQUATIONS IN THE SYSTEM.  THE DEGREE OF
  77 C AN EQUATION IS THE MAXIMUM OF THE DEGREES OF ITS TERMS.  THE
  78 C DEGREE OF A TERM IS THE SUM OF THE DEGREES OF THE VARIABLES.
  79 C THUS, TOTDG = IDEG(1) * ... * IDEG(N) WHERE IDEG(J) =
  80 C MAX {JDEG(J,K) | K=1,...,NUMT(J)} WHERE JDEG(J,K) = KDEG(J,1,K) +
  81 C ... + KDEG(J,N,K).
  82 C
  83 C IFLG1  DETERMINES WHETHER THE SYSTEM IS TO BE AUTOMATICALLY
  84 C SCALED BY  POLSYS  AND WHETHER THE PROJECTIVE TRANSFORMATION
  85 C OF THE SYSTEM IS TO BE AUTOMATICALLY EVOKED BY POLSYS.  SEE
  86 C "ON INPUT" BELOW.
  87 C
  88 C IFLG2, EPSBIG, EPSSML, AND  SSPAR  TELL THE PATH TRACKER
  89 C POLYNF  WHICH PATHS TO TRACK AND SET PARAMETERS FOR THE PATH
  90 C TRACKER.
  91 C
  92 C NUMRR  TELLS  POLSYS  HOW MANY MULTIPLES OF 1000 STEPS TO TRY
  93 C BEFORE ABANDONING A PATH.
  94 C
  95 C NN, MMAXT, TTOTDG, LENWK, LENIWK  GIVE THE DIMENSIONS OF ARRAYS.
  96 C
  97 C THE OUTPUT CONSISTS OF  IFLG1, AND OF  LAMBDA, ROOTS, ARCLEN, AND
  98 C NFE  FOR EACH PATH.  IFLG1  RETURNS INPUT DATA ERROR INFORMATION.
  99 C ROOTS  GIVES THE SOLUTIONS THEMSELVES, WHILE  LAMBDA, ARCLEN,
 100 C AND  NFE  GIVE INFORMATION ABOUT THE ASSOCIATED PATHS.
 101 C
 102 C THE FOLLOWING SUBROUTINES ARE USED DIRECTLY OR INDIRECTLY BY
 103 C POLSYS:
 104 C         SPECIAL FOR POLSYS:
 105 C           POLYP , INITP , STRPTP ,
 106 C           OTPUTP , RHO , RHOJAC ,
 107 C           HFUNP , HFUN1P , GFUNP , FFUNP ,
 108 C           MULP , POWP , DIVP,
 109 C           SCLGNP.
 110 C         FROM THE GENERAL HOMPACK ROUTINES:
 111 C           POLYNF , STEPNF , TANGNF , ROOTNF , ROOT ,
 112 C           QRFAQF, QRSLQF , D1MACH , DDOT , DNRM2.
 113 C
 114 C ON INPUT:
 115 C
 116 C N  IS THE NUMBER OF EQUATIONS AND VARIABLES.
 117 C
 118 C NUMT(1:NN)  IS AN INTEGER ARRAY.  NUMT(J)  IS THE NUMBER OF TERMS
 119 C   IN THE JTH EQUATION FOR J=1 TO N.
 120 C
 121 C COEF(1:NN,1:MMAXT)  IS A REAL ARRAY.  COEF(J,K)  IS
 122 C   THE K-TH COEFFICIENT OF THE J-TH EQUATION FOR J=1 TO N,
 123 C   K=1 TO NUMT(J).
 124 C
 125 C KDEG(1:NN,1:NN+1,1:MMAXT)  IS AN INTEGER ARRAY.
 126 C   KDEG(J,L,K)  IS THE DEGREE OF THE L-TH VARIABLE IN THE K-TH
 127 C   TERM OF THE J-TH EQUATION FOR  J=1 TO N, L=1 TO N, K=1 TO NUMT(J).
 128 C
 129 C IFLG1 =
 130 C   00  IF THE PROBLEM IS TO BE SOLVED WITHOUT
 131 C       CALLING POLSYS' SCALING ROUTINE, SCLGNP, AND
 132 C       WITHOUT USING THE PROJECTIVE TRANSFORMTION.
 133 C
 134 C   01  IF SCALING BUT NO PROJECTIVE TRANSFORMATION IS TO BE USED.
 135 C
 136 C   10  IF NO SCALING BUT PROJECTIVE TRANSFORMATION IS TO BE USED.
 137 C
 138 C   11  IF BOTH SCALING AND PROJECTIVE TRANSFORMATION ARE TO BE USED.
 139 C
 140 C IFLG2(1:TTOTDG)  IS AN INTEGER ARRAY.  IF IFLG2(M) = -2, THEN THE
 141 C   M-TH PATH IS TRACKED.  OTHERWISE THE M-TH PATH IS SKIPPED.
 142 C   THUS, TO FIND ALL SOLUTIONS SET IFLG2(M) = -2 FOR M=1,...,TOTDG.
 143 C   SELECTED PATHS CAN BE RERUN BY SETTING IFLG2(M)=-2 FOR
 144 C   THE PATHS TO BE RERUN AND IFLG(M).NE.-2 FOR THE OTHERS.
 145 C
 146 C EPSBIG  IS THE LOCAL ERROR TOLERANCE ALLOWED THE PATH TRACKER ALONG
 147 C   THE PATH.  ARCRE AND ARCAE (IN  POLYNF ) ARE SET TO  EPSBIG.
 148 C
 149 C EPSSML  IS THE ACCURACY DESIRED FOR THE FINAL SOLUTION.  ANSRE AND
 150 C   ANSAE (IN  POLYNF ) ARE SET TO  EPSSML.
 151 C
 152 C SSPAR(1:8) = (LIDEAL, RIDEAL, DIDEAL, HMIN, HMAX, BMIN, BMAX, P)  IS
 153 C    A VECTOR OF PARAMETERS USED FOR THE OPTIMAL STEP SIZE ESTIMATION.
 154 C    IF  SSPAR(J) .LE. 0.0  ON INPUT, IT IS RESET TO A DEFAULT VALUE
 155 C    BY  POLYNF .  OTHERWISE THE INPUT VALUE OF  SSPAR(J)  IS USED.
 156 C    SEE THE COMMENTS IN  POLYNF  AND IN  STEPNF  FOR MORE INFORMATION
 157 C    ABOUT THESE CONSTANTS.
 158 C
 159 C NUMRR  IS THE NUMBER OF MULTIPLES OF 1000 STEPS THAT WILL BE TRIED
 160 C   BEFORE ABANDONING A PATH.
 161 C
 162 C NN  IS THE DECLARED DIMENSION OF  NUMT  AND OF THE
 163 C   FIRST INDEX OF  COEF  AND  KDEG .  THE SECOND INDEX OF
 164 C   KDEG  AND  ROOTS  IS DIMENSIONED NN+1.  NN  MUST
 165 C   BE GREATER THAN OR EQUAL TO N.
 166 C
 167 C MMAXT  IS THE DECLARED DIMENSION OF THE SECOND INDEX OF
 168 C   COEF  AND THE THIRD INDEX OF  KDEG.  MMAXT  MUST BE
 169 C   GREATER THAN OR EQUAL TO THE MAXIMUM NUMBER OF
 170 C   TERMS IN EACH EQUATION.  IN OTHER WORDS,
 171 C   MMAXT .GE. MAX {NUMT(J) | J=1, ..., N} .
 172 C
 173 C TTOTDG  IS THE DECLARED DIMENSION OF  IFLG2 , LAMBDA , ARCLEN ,
 174 C   NFE , AND OF THE THIRD INDEX OF  ROOTS.  TTOTDG
 175 C   MUST BE GREATER THAN OR EQUAL TO TOTDG, THE TOTAL
 176 C   DEGREE OF THE SYSTEM.
 177 C
 178 C LENWK  IS THE DIMENSION OF THE WORKSPACE  WK .  LENWK  MUST
 179 C   BE GREATER THAN OR EQUAL TO
 180 C     21 + 61*N + 10*N**2 + 7*N*MMAXT + 4*N**2*MMAXT.
 181 C
 182 C LENIWK  IS THE DIMENSION OF THE WORKSPACE  IWK .  LENIWK  MUST BE
 183 C   GREATER THAN OR EQUAL TO  43 + 7*N + N*(N+1)*MMAXT.
 184 C
 185 C
 186 C ON OUTPUT:
 187 C
 188 C N, NUMT, COEF, KDEG, NN, MMAXT, TTOTDG, LENWK, AND  LENIWK
 189 C   ARE UNCHANGED.
 190 C
 191 C IFLG1=
 192 C   -1  IF  NN  IS TOO SMALL.
 193 C   -2  IF  MMAXT  IS TOO SMALL.
 194 C   -3  IF  TTOTDG  IS TOO SMALL.
 195 C   -4  IF  LENWK  IS TOO SMALL.
 196 C   -5  IF  LENIWK  IS TOO SMALL.
 197 C   -6  IF  IFLG1  ON INPUT IS NOT 00 OR 01 OR 10 OR 11.
 198 C   UNCHANGED OTHERWISE.
 199 C
 200 C IFLG2(1:TOTDG)  GIVES INFORMATION ABOUT HOW THE M-TH PATH TERMINATED:
 201 C IFLG2(M) =
 202 C   1   NORMAL RETURN.
 203 C
 204 C   2   SPECIFIED ERROR TOLERANCE CANNOT BE MET.  INCREASE  EPSBIG
 205 C       AND  EPSSML  AND RERUN.
 206 C
 207 C   3   MAXIMUM NUMBER OF STEPS EXCEEDED.  TO TRACK THE PATH FURTHER,
 208 C       INCREASE  NUMRR  AND RERUN THE PATH.  HOWEVER, THE PATH MAY
 209 C       BE DIVERGING, IF THE  LAMBDA  VALUE IS NEAR 1 AND THE  ROOTS
 210 C       VALUES ARE LARGE.
 211 C
 212 C   4   JACOBIAN MATRIX DOES NOT HAVE FULL RANK.  THE ALGORITHM
 213 C       HAS FAILED (THE ZERO CURVE OF THE HOMOTOPY MAP CANNOT BE
 214 C       FOLLOWED ANY FURTHER).
 215 C
 216 C   5   THE TRACKING ALGORITHM HAS LOST THE ZERO CURVE OF THE
 217 C       HOMOTOPY MAP AND IS NOT MAKING PROGRESS.  THE ERROR TOLERANCES
 218 C       EPSBIG  AND  EPSSML  WERE TOO LENIENT.  THE PROBLEM SHOULD BE
 219 C       RESTARTED WITH SMALLER ERROR TOLERANCES.
 220 C
 221 C   6   THE NORMAL FLOW NEWTON ITERATION IN  STEPNF  OR  ROOTNF
 222 C       FAILED TO CONVERGE.  THE ERROR TOLERANCE  EPSBIG  MAY BE TOO
 223 C       STRINGENT.
 224 C
 225 C   7   ILLEGAL INPUT PARAMETERS, A FATAL ERROR.
 226 C
 227 C LAMBDA(M)  IS THE FINAL LAMBDA VALUE FOR THE M-TH PATH, M = 1, ...,
 228 C   TOTDG, WHERE LAMBDA IS THE CONTINUATION PARAMETER.
 229 C
 230 C ROOTS(1,J,M), ROOTS(2,J,M)  ARE THE REAL AND IMAGINARY PARTS
 231 C   OF THE JTH VARIABLE RESPECTIVELY, FOR J = 1,...,N, FOR
 232 C   THE M-TH PATH, FOR M = 1,...,TOTDG.  IF  IFLG1 = 10 OR 11, THEN
 233 C   ROOTS(1,N+1,M)  AND  ROOTS(2,N+1,M)  ARE THE REAL AND
 234 C   IMAGINARY PARTS RESPECTIVELY OF THE PROJECTIVE
 235 C   COORDINATE OF THE SOLUTION.
 236 C
 237 C ARCLEN(M)  IS THE ARC LENGTH OF THE M-TH PATH FOR M = 1, ..., TOTDG.
 238 C
 239 C NFE(M)  IS THE NUMBER OF JACOBIAN MATRIX EVALUATIONS REQUIRED TO
 240 C   TRACK THE M-TH PATH FOR M =1, ..., TOTDG.
 241 C
 242 C ----------------------------------------------------------------------
 243 C TYPE DECLARATIONS FOR INPUT AND OUTPUT
 244       INTEGER N,NUMT,KDEG,IFLG1,IFLG2,NUMRR,NN,MMAXT,
 245      $  TTOTDG,LENWK,LENIWK,NFE,IWK
 246       DOUBLE PRECISION COEF,EPSBIG,EPSSML,SSPAR,LAMBDA,ROOTS,
 247      $  ARCLEN,WK
 248 C
 249 C ARRAY DECLARATIONS FOR INPUT AND OUTPUT
 250       DIMENSION NUMT(NN),KDEG(NN,NN+1,MMAXT),IFLG2(TTOTDG),
 251      $  NFE(TTOTDG),IWK(LENIWK)
 252       DIMENSION COEF(NN,MMAXT),SSPAR(8),LAMBDA(TTOTDG),
 253      $  ROOTS(2,NN+1,TTOTDG), ARCLEN(TTOTDG),WK(LENWK)
 254 C
 255 C TYPE DECLARATIONS FOR VARIABLES
 256       INTEGER I,IDEG,IIDEG,IWKOFF,J,K,L,LENIWW,LENWKK,LIWK,LWK,
 257      $  MAXT,N2,TOTDG,WKOFF
 258 C
 259 C ARRAY DECLARATIONS FOR VARIABLES
 260       DIMENSION LWK(19),LIWK(4),WKOFF(19),IWKOFF(4)
 261 C
 262 C CHECK THAT BASIC DIMENSIONAL PARAMETERS ARE BIG ENOUGH
 263 C
 264       IF(NN.LT.N) THEN
 265         IFLG1=-1
 266         RETURN
 267       END IF
 268       MAXT=0
 269       DO 50 J=1,N
 270         IF(MAXT.LT.NUMT(J))MAXT=NUMT(J)
 271  50   CONTINUE
 272       IF(MMAXT.LT.MAXT) THEN
 273         IFLG1=-2
 274         RETURN
 275       END IF
 276       TOTDG=1
 277       DO 80 J=1,N
 278          IDEG=0
 279          DO 70 K=1,NUMT(J)
 280              IIDEG=0
 281              DO 60 L=1,N
 282                 IIDEG=IIDEG+KDEG(J,L,K)
 283  60          CONTINUE
 284              IF(IIDEG.GT.IDEG)IDEG=IIDEG
 285  70      CONTINUE
 286          TOTDG=TOTDG*IDEG
 287  80   CONTINUE
 288       IF(TTOTDG.LT.TOTDG) THEN
 289         IFLG1=-3
 290         RETURN
 291       END IF
 292       LENWKK = 21 + 61*N + 10*N**2 + 7*N*MMAXT + 4*N**2*MMAXT
 293       IF(LENWK.LT.LENWKK) THEN
 294         IFLG1=-4
 295         RETURN
 296       END IF
 297       LENIWW = 43 +  7*N + N*(N+1)*MMAXT
 298       IF(LENIWK.LT.LENIWW) THEN
 299         IFLG1=-5
 300         RETURN
 301       END IF
 302       IF(IFLG1.NE.0.AND.IFLG1.NE.1.AND.IFLG1.NE.10.AND.IFLG1.NE.11) THEN
 303         IFLG1=-6
 304         RETURN
 305       END IF
 306 C
 307 C VARIABLES THAT ARE PASSED IN ARRAY WK: (LENGTHS ARE IN THE
 308 C INTEGER ARRAY LWK.)
 309 C
 310 C    VARIABLE NAME       LENGTH
 311 C
 312 C    1   PDG               N2
 313 C    2   QDG               N2
 314 C    3   R                 N2
 315 C    4   FACV              N
 316 C    5   CL                2*(N+1)
 317 C    6   Y                 N2+1
 318 C    7   YP                N2+1
 319 C    8   YOLD              N2+1
 320 C    9   YPOLD             N2+1
 321 C   10   QR                N2*(N2+2)
 322 C   11   ALPHA             N2
 323 C   12   TZ                N2+1
 324 C   13   W                 N2+1
 325 C   14   WP                N2+1
 326 C   15   Z0                N2+1
 327 C   16   Z1                N2+1
 328 C   17   SSPAR             8
 329 C   18   PAR               2 + 28*N + 6*N**2 + 7*N*MMAXT + 4*N**2*MMAXT
 330 C
 331 C VARIABLES THAT ARE PASSED IN ARRAY IWK: (LENGTHS ARE IN THE
 332 C INTEGER ARRAY LIWK.)
 333 C
 334 C    VARIABLE NAME       LENGTH
 335 C
 336 C    1   IDEG              N
 337 C    2   ICOUNT            N
 338 C    3   PIVOT             N2+1
 339 C    4   IPAR              42 + 2*N + N*(N+1)*MMAXT
 340 C
 341       N2=2*N
 342       LWK(1)= N2
 343       LWK(2)= N2
 344       LWK(3)= N2
 345       LWK(4)= N
 346       LWK(5)= 2*(N+1)
 347       LWK(6)= N2+1
 348       LWK(7)= N2+1
 349       LWK(8)= N2+1
 350       LWK(9)= N2+1
 351       LWK(10)=N2*(N2+2)
 352       LWK(11)=N2
 353       LWK(12)=N2+1
 354       LWK(13)=N2+1
 355       LWK(14)=N2+1
 356       LWK(15)=N2+1
 357       LWK(16)=N2+1
 358       LWK(17)=8
 359       LWK(18)= 2 + 28*N + 6* N**2 + 7*N*MMAXT + 4* N**2 *MMAXT
 360 C
 361       LIWK(1)=N
 362       LIWK(2)=N
 363       LIWK(3)=2*N+1
 364       LIWK(4)= 42 + 2*N + N*(N+1)*MMAXT
 365 C
 366 C WKOFF AND IWKOFF ARE OFFSETS THAT DEFINE THE VARIABLES LISTED ABOVE
 367 C
 368       WKOFF(1)=1
 369       DO 100 I=2,18
 370           WKOFF(I)=WKOFF(I-1)+LWK(I-1)
 371  100  CONTINUE
 372       IWKOFF(1)=1
 373       DO 200 I=2,4
 374           IWKOFF(I)=IWKOFF(I-1)+LIWK(I-1)
 375  200  CONTINUE
 376       DO 300 J=1,8
 377         WK(WKOFF(17) + (J-1))=SSPAR(J)
 378  300  CONTINUE
 379 C
 380       CALL POLYP(N,NUMT,COEF,KDEG,IFLG1,IFLG2,EPSBIG,EPSSML,
 381      $ NUMRR,NN,MMAXT,TTOTDG,LAMBDA,ROOTS,ARCLEN,NFE,TOTDG,
 382      $ WK( WKOFF( 1)),WK( WKOFF( 2)),WK( WKOFF( 3)),WK( WKOFF( 4)),
 383      $ WK( WKOFF( 5)),WK( WKOFF( 6)),WK( WKOFF( 7)),WK( WKOFF( 8)),
 384      $ WK( WKOFF( 9)),WK( WKOFF(10)),WK( WKOFF(11)),WK( WKOFF(12)),
 385      $ WK( WKOFF(13)),WK( WKOFF(14)),WK( WKOFF(15)),WK( WKOFF(16)),
 386      $ WK( WKOFF(17)),WK( WKOFF(18)),
 387      $IWK(IWKOFF( 1)),IWK(IWKOFF( 2)),IWK(IWKOFF( 3)),IWK(IWKOFF( 4)))
 388 C
 389       RETURN
 390       END