C ALGORITHM 652 (NEW VERSION), COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 13, NO. 3, PP. 281-310. C C HOMPACK is a suite of FORTRAN 77 subroutines for solving nonlinear C systems of equations by homotopy methods. There are subroutines for C fixed point, zero finding, and general homotopy curve tracking problems, C utilizing both dense and sparse Jacobian matrices, and implementing C three different algorithms: ODE-based, normal flow, and augmented C Jacobian. The (driver) subroutines called by the user are given in the C table below, and are well documented internally. The user need not C be concerned with any other subroutines in HOMPACK. C C C Problem type C --------|--------|--------|--------|--------|--------| C x = f(x) | F(x) = 0 |rho(a,lambda,x)=0| C --------|--------|--------|--------|--------|--------| C dense | sparse | dense | sparse | dense | sparse | Algorithm C --------|--------|--------|--------|--------|--------|--------------------- C FIXPDF | FIXPDS | FIXPDF | FIXPDS | FIXPDF | FIXPDS | ODE based C --------|--------|--------|--------|--------|--------|--------------------- C FIXPNF | FIXPNS | FIXPNF | FIXPNS | FIXPNF | FIXPNS | normal flow C --------|--------|--------|--------|--------|--------|--------------------- C FIXPQF | FIXPQS | FIXPQF | FIXPQS | FIXPQF | FIXPQS | augmented Jacobian C --------|--------|--------|--------|--------|--------|--------------------- C C C The sparse subroutines use the packed skyline storage scheme standard in C structural mechanics, but any sparse storage scheme can be used by C replacing some of the low-level HOMPACK routines with user-written C routines. The stepping subroutines STEP?? may be of interest to some C users with special curve tracking needs. C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ORGANIZATIONAL DETAILS. HOMPACK is organized in two different ways: by algorithm/problem type and by subroutine level. There are three levels of subroutines. The top level consists of drivers, one for each problem type and algorithm type. Normally these drivers are called by the user, and the user need know nothing beyond them. They allocate storage for the lower level routines, and all the arrays are variable dimension, so there is no limit on problem size. The second subroutine level implements the major components of the algorithms such as stepping along the homotopy zero curve, computing tangents, and the end game for the solution at lambda = 1 . A sophisticated user might call these routines directly to have complete control of the algorithm, or for some other task such as tracking an arbitrary parametrized curve over an arbitrary parameter range. The lowest subroutine level handles the numerical linear algebra, and includes some BLAS routines. All the linear algebra and associated data structure handling are concentrated in these routines, so a user could incorporate his own data structures by writing his own versions of these low level routines. The organization of HOMPACK by algorithm/problem type is shown in the above table, which lists the driver name for each algorithm and problem type. Using brackets to indicate the three subroutine levels described above, the natural grouping of the HOMPACK routines is: [FIXPDF] [FODE, ROOT, SINTRP, STEPS] [DCPOSE] [FIXPDS] [FODEDS, ROOT, SINTRP, STEPDS] [GMFADS, MFACDS, MULTDS, PCGDS, QIMUDS, SOLVDS] [FIXPNF] [ROOTNF, STEPNF, [TANGNF]] [ROOT] [FIXPNS] [ROOTNS, STEPNS, TANGNS] [GMFADS, MFACDS, MULTDS, PCGDS, PCGNS, QIMUDS, ROOT, SOLVDS] [FIXPQF] [ROOTQF, STEPQF, TANGQF] [QRFAQF, QRSLQF, R1UPQF, UPQRQF] [FIXPQS] [ROOTQS, STEPQS, TANGQS] [GMFADS, MULTDS, PCGQS, SOLVDS] [POLSYS] [POLYNF, POLYP, ROOTNF, STEPNF, TANGNF] [DIVP, FFUNP, GFUNP, HFUNP, HFUN1P, INITP, MULP, OTPUTP, POWP, RHO, RHOJAC, ROOT, SCLGNP, STRPTP] The BLAS subroutines used by HOMPACK are DAXPY, DCOPY, DDOT, DNRM2, DSCAL, D1MACH, IDAMAX. The user written subroutines, of which exactly two must be supplied depending on the driver chosen, are F, FJAC, FJACS, RHO, RHOA, RHOJAC, RHOJS. For testing, there are three main test programs MAINF, MAINP, and MAINS, and one data file INNHP.DAT (read by MAINP). Inquiries should be directed to Layne T. Watson, Department of Computer Science, VPI & SU, Blacksburg, VA 24061; (703) 961-7540; watson@cs.vt.edu ltw@vtopus.cs.vt.edu C MAIN PROGRAM TO TEST FIXPNF, FIXPQF, AND FIXPDF C BROWN'S FUNCTION, ZERO FINDING. C C THIS PROGRAM TESTS THE HOMPACK ROUTINES FIXPNF, FIXPQF, AND C FIXPDF. THE USER MAY INSERT CALLS TO A SYSTEM TIMER AT THE C DESIGNATED LOCATIONS IN ORDER TO GET EXECUTION TIME FOR THESE C ROUTINES. C C THE MODIFICATIONS TO BE MADE FOR THE SYSTEM TIMER ARE INDICATED C BY A LINE OF M'S, E.G. CMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM C C C THE OUTPUT FROM THIS ROUTINE SHOULD BE AS FOLLOWS, WITH THE C EXECUTION TIMES CORRESPONDING TO A VAX 11/785. C C TESTING FIXPQF C C LAMBDA = 1.00000000 FLAG = 1 6 JACOBIAN EVALUATIONS C EXECUTION TIME(SECS) = 0.44 ARCLEN = 2.693 C 1.00000000E+00 1.00000000E+00 1.00000000E+00 C 1.00000000E+00 1.00000000E+00 C C TESTING FIXPNF C C LAMBDA = 1.00000000 FLAG = 1 22 JACOBIAN EVALUATIONS C EXECUTION TIME(SECS) = 0.19 ARCLEN = 2.682 C 1.00000000E+00 1.00000000E+00 1.00000000E+00 C 1.00000000E+00 1.00000000E+00 C C TESTING FIXPDF C C LAMBDA = 1.00000000 FLAG = 1 71 JACOBIAN EVALUATIONS C EXECUTION TIME(SECS) = 0.57 ARCLEN = 2.712 C 1.00000000E+00 1.00000000E+00 1.00000000E+00 C 1.00000000E+00 1.00000000E+00 C C C PROGRAM TEST1 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION WT(101),PHI(101,16),P(101) DOUBLE PRECISION ARCLEN,QT(101,101),R(101*52),F0(101) DOUBLE PRECISION F1(101),DZ(101),T(101) DOUBLE PRECISION Y(101),W(101),WP(101),Z0(101),Z1(101), + YP(101),YOLD(101),YPOLD(101),A(100),QR(101,102), + ALPHA(100),TZ(101),SSPAR(8),YSAV(101),PAR(1) INTEGER PIVOT(101),CODE,TIME,IPAR(1),N,NDIMA,NFE,TRACE, + IFLAG,II,J,NP1 CHARACTER*6 NAME REAL DTIME COMMON /SIZE/ N C C TEST EACH OF THE THREE ALGORITHMS. C DO 60 II=1,3 C C INITIALIZE TIMER VARIABLES. C CODE=2 TIME=0 DTIME=0.0 C C DEFINE ARGUMENTS FOR CALL TO HOMPACK PROCEDURE. C N=5 NP1=N+1 ARCRE=0.5D-4 ARCAE=0.5D-4 ANSRE=1.0D-10 ANSAE=1.0D-10 TRACE=0 DO 30 J=1,8 30 SSPAR(J)=0.0 IFLAG=-1 DO 40 J=2,NP1 40 Y(J)=0.0D0 C CMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM C C INSERT CALL TO INITIALIZE SYSTEM TIMER HERE. FOR EXAMPLE, FOR C THE VAX, THE FOLLOWING STATEMENT IS USED. C C CALL LIB$INIT_TIMER C CMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM C C CALL TO HOMPACK ROUTINE. C IF (II .EQ. 1) THEN NAME='FIXPQF' CALL FIXPQF(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,A,NFE, + ARCLEN,YP,YOLD,YPOLD,QT,R,F0,F1,Z0,DZ,W,T,YSAV, + SSPAR,PAR,IPAR) ELSE IF (II .EQ. 2) THEN NAME='FIXPNF' CALL FIXPNF(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,A,NFE, + ARCLEN,YP,YOLD,YPOLD,QR,ALPHA,TZ,PIVOT,W,WP,Z0,Z1, + SSPAR,PAR,IPAR) ELSE NAME='FIXPDF' CALL FIXPDF(N,Y,IFLAG,ARCRE,ANSRE,TRACE,A,NDIMA,NFE, + ARCLEN,YP,YPOLD,QR,ALPHA,TZ,PIVOT,WT,PHI,P,PAR,IPAR) END IF C CMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM C C INSERT CALL TO RETURN EXECUTION TIME IN SECONDS IN DTIME. C FOR EXAMPLE, THE VAX STATEMENTS ARE AS FOLLOWS. C CALL LIB$STAT_TIMER(CODE,TIME) C DTIME=TIME/100.0 C CMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM C WRITE (6,45) NAME 45 FORMAT (//,8X,'TESTING',1X,6A) WRITE (6,50) Y(1),IFLAG,NFE,DTIME,ARCLEN,(Y(J),J=2,NP1) 50 FORMAT(//' LAMBDA =',F11.8,' FLAG =',I2,I8,' JACOBIAN ', + 'EVALUATIONS',/,1X,'EXECUTION TIME(SECS) =',F10.2,4X, + 'ARCLEN =',F10.3/(1X,1P,3E16.8)) 60 CONTINUE 400 STOP END SUBROUTINE F(X,V) C******************************************************************** C C SUBROUTINE F(X,V) -- EVALUATES BROWN'S FUNCTION AT THE POINT C X, AND RETURNS THE VALUE IN V. C C******************************************************************** C DOUBLE PRECISION X(1),V(1),PROD,SUM INTEGER J,N COMMON /SIZE/ N PROD=1.0D0 DO 10 J=1,N 10 PROD=PROD*X(J) V(1)=PROD-1.0D0 SUM=0.0D0 DO 20 J=1,N 20 SUM=SUM+X(J) SUM=SUM-DBLE(N+1) DO 30 J=2,N 30 V(J)=SUM+X(J) RETURN END SUBROUTINE FJAC(X,V,K) C******************************************************************** C C SUBROUTINE FJAC(X,V,K) -- EVALUATES THE K-TH COLUMN OF C THE JACOBIAN MATRIX FOR BROWN'S FUNCTION EVALUATED AT C THE POINT X, RETURNING THE VALUE IN V. C C******************************************************************** C DOUBLE PRECISION X(1),V(1),PROD INTEGER J,K,N COMMON /SIZE/ N PROD=1.0D0 DO 10 J=1,K-1 10 PROD=PROD*X(J) DO 15 J=K+1,N 15 PROD=PROD*X(J) V(1)=PROD DO 20 J=2,N 20 V(J)=1.0D0 IF (K .GT. 1) V(K)=V(K)+1.0D0 RETURN END C C MAIN ROUTINE TO TEST POLSYS C C THIS ROUTINE REQUIRES ONE INPUT FILE, READ AS UNIT FOR007. C C A SAMPLE INPUT FILE AND ASSOCIATED OUTPUT ARE GIVEN C IN THE COMMENTS THAT FOLLOW. THIS SAMPLE PROBLEM IS C CITED IN THE HOMPACK REPORT. C C***** SAMPLE INPUT DATA: C TWO QUADRICS, NO SOLUTIONS AT INFINITY, TWO REAL SOLUTIONS. C 00001 IFLGHM C 00001 IFLGSC C 4 ITOTDG C 1.D-04 EPSBIG C 1.D-14 EPSSML C 1.D-00 SSPAR(5) C 10 NUMRR C 2 N C 00006 NUMTRM(1) C 00002 DEG(1,1,1) C 00000 DEG(1,2,1) C -.00098D 00 C 00000 DEG(1,1,2) C 00002 DEG(1,2,2) C 978000.D 00 C 00001 DEG(1,1,3) C 00001 DEG(1,2,3) C -9.8D 00 C 00001 DEG(1,1,4) C 00000 DEG(1,2,4) C -235.0D 00 C 00000 DEG(1,1,5) C 00001 DEG(1,2,5) C 88900.0D 00 C 00000 DEG(1,1,6) C 00000 DEG(1,2,6) C -1.000D 00 C 00006 NUMTRM(2) C 00002 DEG(2,1,1) C 00000 DEG(2,2,1) C -.0100D 00 C 00000 DEG(2,1,2) C 00002 DEG(2,2,2) C -.9840D 00 C 00001 DEG(2,1,3) C 00001 DEG(2,2,3) C -29.70D 00 C 00001 DEG(2,1,4) C 00000 DEG(2,2,4) C .00987D 00 C 00000 DEG(2,1,5) C 00001 DEG(2,2,5) C -.1240D 00 C 00000 DEG(2,1,6) C 00000 DEG(2,2,6) C -.2500D 00 C***** END OF SAMPLE INPUT DATA. C C***** ASSOCIATED SAMPLE OUTPUT: C C C POLYS TEST ROUTINE 5/20/85 C C C TWO QUADRICS PBHP0403, NO SOLUTIONS AT INFINITY ......... C C IF IFLGHM=1,HOMOGENEOUS;IF IFLGHM=2,INHOMOGENEOUS;IFLGHM= 1 C C IF IFLGSC=1,SCLGEN USED; IF IFLGSC=2, NO SCALING; IFLGSC= 1 C C ITOTDG= 4 C C EPSBIG,EPSSML = 0.100000000000000D-03 0.100000000000000D-13 C NUMBER OF EQUATIONS = 2 C C C NUMBER OF RECALLS WHEN IFLAG=3: 40 C C C C ****** COEFFICIENT TABLEAU ****** C C C NUMT( 1)= 6 C KDEG( 1, 1, 1)= 2 C KDEG( 1, 2, 1)= 0 C COEF( 1, 1)=-0.980000000000000D-03 C KDEG( 1, 1, 2)= 0 C KDEG( 1, 2, 2)= 2 C COEF( 1, 2)= 0.978000000000000D+06 C KDEG( 1, 1, 3)= 1 C KDEG( 1, 2, 3)= 1 C COEF( 1, 3)=-0.980000000000000D+01 C KDEG( 1, 1, 4)= 1 C KDEG( 1, 2, 4)= 0 C COEF( 1, 4)=-0.235000000000000D+03 C KDEG( 1, 1, 5)= 0 C KDEG( 1, 2, 5)= 1 C COEF( 1, 5)= 0.889000000000000D+05 C KDEG( 1, 1, 6)= 0 C KDEG( 1, 2, 6)= 0 C COEF( 1, 6)=-0.100000000000000D+01 C C C NUMT( 2)= 6 C KDEG( 2, 1, 1)= 2 C KDEG( 2, 2, 1)= 0 C COEF( 2, 1)=-0.100000000000000D-01 C KDEG( 2, 1, 2)= 0 C KDEG( 2, 2, 2)= 2 C COEF( 2, 2)=-0.984000000000000D+00 C KDEG( 2, 1, 3)= 1 C KDEG( 2, 2, 3)= 1 C COEF( 2, 3)=-0.297000000000000D+02 C KDEG( 2, 1, 4)= 1 C KDEG( 2, 2, 4)= 0 C COEF( 2, 4)= 0.987000000000000D-02 C KDEG( 2, 1, 5)= 0 C KDEG( 2, 2, 5)= 1 C COEF( 2, 5)=-0.124000000000000D+00 C KDEG( 2, 1, 6)= 0 C KDEG( 2, 2, 6)= 0 C COEF( 2, 6)=-0.250000000000000D+00 C C C C C PATH NUMBER = 1 C C FINAL VALUES FOR PATH C C ARCLEN = 0.100553311312353D+02 C NFE = 53 C IFLG2 = 1 C T = 0.100000000000000D+01 C X = 0.234233851959126D+04 0.791152831437911D-11 C X =-0.788344824094138D+00-0.268347762088076D-14 C X =-0.949359459408658D-02-0.106447550900261D-02 C X = C C C PATH NUMBER = 2 C C FINAL VALUES FOR PATH C C ARCLEN = 0.172112868960496D+01 C NFE = 37 C IFLG2 = 1 C T = 0.100000000000000D+01 C X = 0.161478579234367D-01 0.168496955498881D+01 C X = 0.267994739614462D-03 0.442802993973661D-02 C X =-0.381948972942403D+00 0.372068943457283D+00 C X = C C C PATH NUMBER = 3 C C FINAL VALUES FOR PATH C C ARCLEN = 0.202329539135269D+01 C NFE = 35 C IFLG2 = 1 C T = 0.100000000000000D+01 C X = 0.161478579234362D-01-0.168496955498881D+01 C X = 0.267994739614461D-03-0.442802993973661D-02 C X =-0.329370493847660D+00 0.556619775523013D+00 C X = C C C PATH NUMBER = 4 C C FINAL VALUES FOR PATH C C ARCLEN = 0.416327291917901D+01 C NFE = 46 C IFLG2 = 1 C T = 0.100000000000000D+01 C X = 0.908921229615394D-01-0.111985846294633D-14 C X =-0.911497098197500D-01 0.117962440099502D-17 C X =-0.573673395727962D-01 0.136243663709219D+00 C X = C C C TOTAL NFE OVER ALL PATHS = 171 C C***** END OF ASSOCIATED SAMPLE OUTPUT. C C ************************************************************* C C PROGRAM DESCRIPTION: 1. READS IN DATA. C 2. GENERATES POLSYS INPUT. C 3. CALLS POLSYS. C 4. WRITES POLSYS OUTPUT. C C DIMENSIONS SHOULD BE SET AS FOLLOWS: C C DIMENSION NUMT(NN),COEF(NN,MMAXT),KDEG(NN,NN+1,MMAXT) C DIMENSION IFLG2(TTOTDG) C DIMENSION LAMBDA(TTOTDG),ROOTS(2,NN+1,TTOTDG),ARCLEN(TTOTDG), C + NFE(TTOTDG) C DIMENSION WK(LENWK),IWK(LENIWK) C WHERE: C N IS THE NUMBER OF EQUATIONS. NN .GE. N. C MAXT IS THE MAXIMUM NUMBER OF TERMS IN ANY ONE EQUATION. C MMAXT .GE. MAXT. C TOTDG IS THE TOTAL DEGREE OF THE SYSTEM. TTOTDG .GE. TOTDG. 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 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 THIS TEST CODE HAS DIMENSIONS SET AS FOLLOWS: C C NN=10, MMAXT=30, TTOTDG=999 C LENWK = 21 + 610 + 1000 + 2100 + 12000 = 15731 C LENIWK = 43 + 70 + 3300 = 3413 C PROGRAM TESTP INTEGER IFLG1,IFLG2,IFLGHM,IFLGSC,ITEST,ITOTIT,IWK,J,K,KDEG, + L,LENIWK,LENWK,M,MMAXT,N,NFE,NN,NP1,NT,NUMRR,NUMT,TTOTDG DOUBLE PRECISION ARCLEN,COEF,EPSBIG,EPSSML,LAMBDA,ROOTS, + SSPAR,WK CHARACTER*72 TITLE C DIMENSION ARCLEN(999),COEF(10,30),IFLG2(999),IWK(3413), + KDEG(10,11,30),LAMBDA(999),NFE(999),NUMT(10),ROOTS(2,11,999), + SSPAR(8),WK(15731) C NN=10 MMAXT=30 TTOTDG=999 LENWK=15731 LENIWK=3413 C C OPEN (UNIT=7,FILE='INNHP.DAT',STATUS='UNKNOWN') OPEN (UNIT=6,FILE='OUTHP.DAT',STATUS='UNKNOWN') C SSPAR(1)=.0 SSPAR(2)=.0 SSPAR(3)=.0 SSPAR(4)=.0 SSPAR(6)=.0 SSPAR(7)=.0 SSPAR(8)=.0 C 1000 FORMAT(I5) 2000 FORMAT(D22.15) C WRITE(6,10) 10 FORMAT( ' POLSYS TEST ROUTINE 8/12/85',//) C READ(7,*) TITLE WRITE(6,21) TITLE 21 FORMAT(' ',A72) C READ(7,1000) IFLGHM READ(7,1000) IFLGSC READ(7,1000) TTOTDG C READ(7,2000) EPSBIG READ(7,2000) EPSSML READ(7,2000) SSPAR(5) READ(7,1000) NUMRR READ(7,1000) N C WRITE(6,100) IFLGHM 100 FORMAT(/ +' IF IFLGHM=1,HOMOGENEOUS;IF IFLGHM=0,INHOMOGENEOUS;IFLGHM=',I2) WRITE(6,102) IFLGSC 102 FORMAT(/ +' IF IFLGSC=1,SCLGNP USED; IF IFLGSC=0, NO SCALING; IFLGSC=',I2) WRITE(6,104) TTOTDG 104 FORMAT(/,' TTOTDG=',I5) C C WRITE(6,106) EPSBIG,EPSSML,SSPAR(5),N 106 FORMAT(/,' EPSBIG,EPSSML =',2D22.15, + //,' SSPAR(5) =',D22.15, + //,' NUMBER OF EQUATIONS =',I5) WRITE(6,108) NUMRR 108 FORMAT(/,' NUMBER OF RECALLS WHEN IFLAG=3: ',I5) C NP1=N+1 C C NOTE THAT THE DEGREES OF VARIABLES IN EACH TERM OF EACH EQUATION C ARE DEFINED BY THE FOLLOWING INDEXING SCHEME: C C KDEG(J, L, K) C C ^ ^ ^ C C E V T C Q A E C U R R C A I M C T A C I B C O L C N E C C WRITE(6,200) 200 FORMAT(//,' ****** COEFFICIENT TABLEAU ******') C DO 202 J=1,N C WRITE(6,205) 205 FORMAT(/) C READ(7,1000) NUMT(J) WRITE(6,210) J,NUMT(J) 210 FORMAT(' NUMT(',I2,')=',I5) C NT=NUMT(J) C DO 215 K=1,NT C DO 218 L=1,N READ(7,1000) KDEG(J,L,K) WRITE(6,220) J,L,K,KDEG(J,L,K) 220 FORMAT(' KDEG(',I2,',',I2,',',I2,')=',I5) C 218 CONTINUE C READ(7,2000) COEF(J,K) WRITE(6,230) J,K,COEF(J,K) 230 FORMAT(' COEF(',I2,',',I2,')=',D22.15) C 215 CONTINUE C 202 CONTINUE C WRITE(6,205) WRITE(6,205) C IFLG1=10*IFLGHM+IFLGSC C DO 235 M=1,TTOTDG IFLG2(M)=-2 235 CONTINUE C CALL POLSYS(N,NUMT,COEF,KDEG,IFLG1,IFLG2, + EPSBIG,EPSSML,SSPAR, + NUMRR,NN,MMAXT,TTOTDG,LENWK,LENIWK, + LAMBDA,ROOTS,ARCLEN,NFE,WK,IWK) WRITE(6,240) IFLG1 240 FORMAT(/,' IFLG1=',I5,/) C ITOTIT=0 DO 250 M=1,TTOTDG C ITOTIT=ITOTIT+NFE(M) C WRITE(6,260) M 260 FORMAT(' PATH NUMBER =',I5) WRITE(6,270) 270 FORMAT(/' FINAL VALUES FOR PATH'/) C WRITE(6,280) ARCLEN(M) 280 FORMAT(' ARCLEN =',D22.15) WRITE(6,290) NFE(M) 290 FORMAT(' NFE =',I5) WRITE(6,300) IFLG2(M) 300 FORMAT(' IFLG2 =',I5) C C******************************* C C DESIGNATE SOLUTIONS "REAL" OR "COMPLEX" C ITEST=0 DO 310 J=1,N IF(ABS(ROOTS(2,J,M)).GE.1.E-4) ITEST=1 310 CONTINUE IF( ITEST.EQ.1) THEN WRITE(6,779) 779 FORMAT(' COMPLEX SOLUTION ') ELSE WRITE(6,780) 780 FORMAT(' REAL SOLUTION ') END IF C C******************************* C C C DESIGNATE SOLUTION "FINITE" OR "INFINITE" C IF( ABS(ROOTS(1,NP1,M))+ABS(ROOTS(2,NP1,M)) .LT. 1.E-6) THEN WRITE(6,781) 781 FORMAT(' INFINITE SOLUTION ') ELSE WRITE(6,782) 782 FORMAT(' FINITE SOLUTION ') END IF C C******************************* C WRITE(6,320) LAMBDA(M),(ROOTS(1,J,M),ROOTS(2,J,M),J=1,N) 320 FORMAT(' LAMBDA =',D22.15,/,10(' X =',2D22.15,/)) WRITE(6,330) ROOTS(1,NP1,M),ROOTS(2,NP1,M) 330 FORMAT(/,' XNP1 =',2D22.15,/) C WRITE(6,205) C 250 CONTINUE C WRITE(6,400) ITOTIT 400 FORMAT(' TOTAL NFE OVER ALL PATHS = ',I10) C C STOP END C MAIN PROGRAM TO TEST FIXPQS, FIXPNS, AND FIXPDS C C THIS PROGRAM TESTS THE HOMPACK ROUTINES FIXPNS, FIXPQS, AND C FIXPDS. THE USER MAY INSERT CALLS TO A SYSTEM TIMER AT THE C DESIGNATED LOCATIONS IN ORDER TO GET EXECUTION TIME FOR THESE C ROUTINES. C C THE MODIFICATIONS TO BE MADE FOR THE SYSTEM TIMER ARE INDICATED C BY A LINE OF M'S, E.G. CMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM C C C THE OUTPUT FROM THIS ROUTINE SHOULD BE AS FOLLOWS, WITH THE C EXECUTION TIMES CORRESPONDING TO A VAX 11/785. C C TESTING FIXPQS C C LAMBDA = 1.00000000 FLAG = 1 33 JACOBIAN EVALUATIONS C ARC LENGTH = 1.274 EXECUTION TIME(SECS) = 2.31 C 4.00864019E-01 2.65454893E-01 8.40421103E-02 4.83042527E-01 C 3.01797132E-01 2.32508994E-01 4.96639853E-01 3.00908894E-01 C C TESTING FIXPNS C C LAMBDA = 1.00000000 FLAG = 1 20 JACOBIAN EVALUATIONS C ARC LENGTH = 1.275 EXECUTION TIME(SECS) = 1.04 C 4.00864019E-01 2.65454893E-01 8.40421103E-02 4.83042527E-01 C 3.01797132E-01 2.32508994E-01 4.96639853E-01 3.00908894E-01 C C TESTING FIXPDS C C LAMBDA = 1.00000000 FLAG = 1 70 JACOBIAN EVALUATIONS C ARC LENGTH = 1.281 EXECUTION TIME(SECS) = 1.78 C 4.00864019E-01 2.65454893E-01 8.40421103E-02 4.83042527E-01 C 3.01797132E-01 2.32508994E-01 4.96639853E-01 3.00908894E-01 C C PROGRAM TEST1 IMPLICIT DOUBLE PRECISION(A-H,O-Z) DOUBLE PRECISION Y(9), + YP(9),YOLD(9),YPOLD(9),A(8),QR(18),WORK(200), + SSPAR(8),PAR(1),PP(8),RHOVEC(9),Z0(9),DZ(9),T(9), + WT(9),PHI(9,16),P(9) INTEGER PIVOT(10),IPAR(1) INTEGER IFLAG,II,J,LENQR,N,NFE,NP1,NDIMA,TRACE DOUBLE PRECISION ARCRE,ARCAE,ANSRE,ANSAE,ARCLEN CHARACTER*6 NAME INTEGER TIME,CODE REAL DTIME C C TEST EACH OF THE THREE ALGORITHMS. C DO 60 II=1,3 C C INITIALIZE TIMER VARIABLES. C CODE=2 TIME=0 DTIME=0.0 C C DFEFINE ARGUMENTS FOR CALL TO HOMPACK PROCEDURE. C N=8 DO 7 J=1,8 7 SSPAR(J)=0.0D0 ARCRE=.5D-4 ARCAE=.5D-4 ANSRE=1.0D-12 ANSAE=1.0D-12 TRACE=0 IFLAG=-1 LENQR=18 NP1=N+1 DO 40 J=1,N 40 Y(J)=0.5D0 C CMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM C C INSERT CALL TO INITIALIZE SYSTEM TIMER HERE. FOR EXAMPLE, FOR C THE VAX, THE FOLLOWING STATEMENT IS USED. C C CALL LIB$INIT_TIMER C CMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM C C CALL TO HOMPACK ROUTINE. C C IF (II .EQ. 1) THEN NAME='FIXPQS' CALL FIXPQS(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE, + A,NFE,ARCLEN,YP,YOLD,YPOLD,QR,LENQR,PIVOT,PP,RHOVEC, + Z0,DZ,T,WORK,SSPAR,PAR,IPAR) ELSE IF (II .EQ. 2) THEN NAME='FIXPNS' CALL FIXPNS(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,A, + NFE,ARCLEN,YP,YOLD,YPOLD,QR,LENQR,PIVOT,WORK, + SSPAR,PAR,IPAR) ELSE NAME='FIXPDS' CALL FIXPDS(N,Y,IFLAG,ARCRE,ANSRE,TRACE,A,NDIMA,NFE, + ARCLEN,YP,YPOLD,QR,LENQR,PIVOT,PP,WORK,WT,PHI,P, + PAR,IPAR) END IF C CMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM C C INSERT CALL TO RETURN EXECUTION TIME IN SECONDS IN DTIME. C FOR EXAMPLE, THE VAX STATEMENTS ARE AS FOLLOWS. C CALL LIB$STAT_TIMER(CODE,TIME) C DTIME=TIME/100.0 C CMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM C WRITE (6,45) NAME 45 FORMAT(//,8X,'TESTING',1X,6A) WRITE (6,50) Y(NP1),IFLAG,NFE,ARCLEN,DTIME,(Y(J),J=1,N) 50 FORMAT(/' LAMBDA =',F11.8,' FLAG =',I2,I8,' JACOBIAN ', + 'EVALUATIONS',/,1X,' ARC LENGTH =',F8.3, + ' EXECUTION TIME(SECS) =',F10.2/(1X,1P,4E16.8)) 60 CONTINUE STOP END SUBROUTINE F(X,V) C C**************************************************************** C C SUBROUTINE F(X,V) -- COMPUTES F AT THE POINT X, C RETURNING THE VALUE IN V. C C**************************************************************** DOUBLE PRECISION X(8),V(8) V(1)=X(1)**3+6.0*X(2)*X(3)-1+2.0*X(1) V(2)=6.0*X(1)*X(3)+X(2)**4*X(5)-1+3.0*X(2) V(3)=6.0*X(1)*X(2)+X(3)*X(5)-1+4.0*X(3) V(4)=X(4)**3*X(8)-1+2.0*X(4) V(5)=X(2)**5/5.0 + X(3)**2/2.0 + X(8)*X(5)-1+3.0*X(5) V(6)=X(6)*X(8)-1+4.0*X(6) V(7)=X(7)**2*X(8)**3-1+2.0*X(7) V(8)=X(4)**4/4.0 + X(5)**2/2.0 + X(6)**2/2.0 + X(7)**3* + X(8)**2-1+3.0*X(8) RETURN END SUBROUTINE FJACS(X,QR,LENQR,PIVOT) C****************************************************************** C C SUBROUTINE FJACS(X,QR,LENQR,PIVOT) C C -- COMPUTES THE JACOBIAN OF F AT THE POINT X, RETURNING C THE JACOBIAN MATRIX IN PACKED SKYLINE FORM IN THE C ARRAYS QR, AND PIVOT. C C***************************************************************** DOUBLE PRECISION X(8),QR(LENQR) INTEGER LENQR,PIVOT(9) PIVOT(1)=1 PIVOT(2)=2 PIVOT(3)=4 PIVOT(4)=7 PIVOT(5)=8 PIVOT(6)=12 PIVOT(7)=13 PIVOT(8)=14 PIVOT(9)=19 QR(1)=3.0*X(1)**2+2.0 QR(2)=4.0*X(2)**3*X(5)+3.0 QR(3)=6.0*X(3) QR(4)=X(5)+4.0 QR(5)=6.0*X(1) QR(6)=6.0*X(2) QR(7)=3.0*X(4)**2*X(8)+2.0 QR(8)=X(8)+3.0 QR(9)=.0 QR(10)=X(3) QR(11)=X(2)**4 QR(12)=X(8)+4.0 QR(13)=2.0*X(7)*X(8)**3+2.0 QR(14)=2.0*X(7)**3*X(8)+3.0 QR(15)=3.0*X(7)**2*X(8)**2 QR(16)=X(6) QR(17)=X(5) QR(18)=X(4)**3 RETURN END C C HOMPACK is a suite of FORTRAN 77 subroutines for solving nonlinea C systems of equations by homotopy methods. There are subroutines for C fixed point, zero finding, and general homotopy curve tracking problem C utilizing both dense and sparse Jacobian matrices, and implementing C three different algorithms: ODE-based, normal flow, and augmented C Jacobian. The (driver) subroutines called by the user are given in th C table below, and are well documented internally. The user need not C be concerned with any other subroutines in HOMPACK. C C C Problem type C --------|--------|--------|--------|--------|--------| C x = f(x) | F(x) = 0 |rho(a,lambda,x)=0| C --------|--------|--------|--------|--------|--------| C dense | sparse | dense | sparse | dense | sparse | Algorithm C --------|--------|--------|--------|--------|--------|---------------- C FIXPDF | FIXPDS | FIXPDF | FIXPDS | FIXPDF | FIXPDS | ODE based C --------|--------|--------|--------|--------|--------|---------------- C FIXPNF | FIXPNS | FIXPNF | FIXPNS | FIXPNF | FIXPNS | normal flow C --------|--------|--------|--------|--------|--------|---------------- C FIXPQF | FIXPQS | FIXPQF | FIXPQS | FIXPQF | FIXPQS | augmented Jacob C --------|--------|--------|--------|--------|--------|---------------- C C C The sparse subroutines use the packed skyline storage scheme standard C structural mechanics, but any sparse storage scheme can be used by C replacing some of the low-level HOMPACK routines with user-written C routines. The stepping subroutines STEP?? may be of interest to some C users with special curve tracking needs. C C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C DOUBLE PRECISION FUNCTION D1MACH(I) C***BEGIN PROLOGUE D1MACH C***DATE WRITTEN 750101 (YYMMDD) C***REVISION DATE 870717 (YYMMDD) C***CATEGORY NO. Q3 C***KEYWORDS MACHINE CONSTANTS C***AUTHOR FOX, P. A., (BELL LABS) C HALL, A. D., (BELL LABS) C SCHRYER, N. L., (BELL LABS) C WATSON, L. T., (VPI & SU) C***PURPOSE Returns double precision machine dependent constants C***DESCRIPTION C D1MACH can be used to obtain machine-dependent parameters C for the local machine environment. It is a function C subprogram with one (input) argument, and can be called C as follows, for example C C D = D1MACH(I) C C where I=1,...,5. The (output) value of D above is C determined by the (input) value of I. The results for C various values of I are discussed below. C C Double-precision machine constants C D1MACH( 1) = B**(EMIN-1), the smallest positive magnitude. C D1MACH( 2) = B**EMAX*(1 - B**(-T)), the largest magnitude. C D1MACH( 3) = B**(-T), the smallest relative spacing. C D1MACH( 4) = B**(1-T), the largest relative spacing. C D1MACH( 5) = LOG10(B) C***REFERENCES FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A C PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188. C***END PROLOGUE D1MACH C INTEGER I INTEGER SMALL(4) INTEGER LARGE(4) INTEGER RIGHT(4) INTEGER DIVER(4) INTEGER LOG10(4) C DOUBLE PRECISION DMACH(5) C EQUIVALENCE (DMACH(1),SMALL(1)) EQUIVALENCE (DMACH(2),LARGE(1)) EQUIVALENCE (DMACH(3),RIGHT(1)) EQUIVALENCE (DMACH(4),DIVER(1)) EQUIVALENCE (DMACH(5),LOG10(1)) C SAVE DMACH C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C DATA SMALL(1) / ZC00800000 / C DATA SMALL(2) / Z000000000 / C C DATA LARGE(1) / ZDFFFFFFFF / C DATA LARGE(2) / ZFFFFFFFFF / C C DATA RIGHT(1) / ZCC5800000 / C DATA RIGHT(2) / Z000000000 / C C DATA DIVER(1) / ZCC6800000 / C DATA DIVER(2) / Z000000000 / C C DATA LOG10(1) / ZD00E730E7 / C DATA LOG10(2) / ZC77800DC0 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. C C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O0000000000000000 / C C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O0007777777777777 / C C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 / C C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. C C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O7770000000000000 / C C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O7777777777777777 / C C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES. C C DATA SMALL(1) / 00604000000000000000B / C DATA SMALL(2) / 00000000000000000000B / C DATA LARGE(1) / 37767777777777777777B / C DATA LARGE(2) / 37167777777777777777B / C DATA RIGHT(1) / 15604000000000000000B / C DATA RIGHT(2) / 15000000000000000000B / C DATA DIVER(1) / 15614000000000000000B / C DATA DIVER(2) / 15010000000000000000B / C DATA LOG10(1) / 17164642023241175717B / C DATA LOG10(2) / 16367571421742254654B / C C MACHINE CONSTANTS FOR THE CDC CYBER SERIES. C C DATA SMALL(1) / O"00604000000000000000" / C DATA SMALL(2) / O"00000000000000000000" / C DATA LARGE(1) / O"37767777777777777777" / C DATA LARGE(2) / O"37167777777777777777" / C DATA RIGHT(1) / O"15604000000000000000" / C DATA RIGHT(2) / O"15000000000000000000" / C DATA DIVER(1) / O"15614000000000000000" / C DATA DIVER(2) / O"15010000000000000000" / C DATA LOG10(1) / O"17164642023241175717" / C DATA LOG10(2) / O"16367571421742254654" / C C MACHINE CONSTANTS FOR CONVEX C-1 C C DATA SMALL(1),SMALL(2) / '00100000'X, '00000000'X / C DATA LARGE(1),LARGE(2) / '7FFFFFFF'X, 'FFFFFFFF'X / C DATA RIGHT(1),RIGHT(2) / '3CC00000'X, '00000000'X / C DATA DIVER(1),DIVER(2) / '3CD00000'X, '00000000'X / C DATA LOG10(1),LOG10(2) / '3FF34413'X, '509F79FF'X / C C MACHINE CONSTANTS FOR THE CRAY 1 (ERIC GROSSE) C C DATA SMALL(1) / 201354000000000000000B / C DATA SMALL(2) / 000000000000000000000B / C DATA LARGE(1) / 577767777777777777777B / C DATA LARGE(2) / 000007777777777777776B / C DATA RIGHT(1) / 376434000000000000000B / C DATA RIGHT(2) / 000000000000000000000B / C DATA DIVER(1) / 376444000000000000000B / C DATA DIVER(2) / 000000000000000000000B / C DATA LOG10(1) / 377774642023241175717B / C DATA LOG10(2) / 000007571421742254654B / C C MACHINE CONSTANTS FOR THE CRAY 1 (SLATEC LIBRARY) C C DATA SMALL(1) / 200004000000000000000B / C DATA SMALL(2) / 000000000000000000000B / C DATA LARGE(1) / 577777777777777777777B / C DATA LARGE(2) / 000007777777777777777B / C DATA RIGHT(1) / 377214000000000000000B / C DATA RIGHT(2) / 000000000000000000000B / C DATA DIVER(1) / 377224000000000000000B / C DATA DIVER(2) / 000000000000000000000B / C DATA LOG10(1) / 377774642023241175717B / C DATA LOG10(2) / 000007571421742254654B / C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD - C STATIC DMACH(5) C C DATA SMALL/20K,3*0/,LARGE/77777K,3*177777K/ C DATA RIGHT/31420K,3*0/,DIVER/32020K,3*0/ C DATA LOG10/40423K,42023K,50237K,74776K/ C C MACHINE CONSTANTS FOR THE DATA GENERAL MV/8000, 10000 C C DATA SMALL(1),SMALL(2) / 1048576, 0 / C DATA LARGE(1),LARGE(2) / 2147483647, -1 / C DATA RIGHT(1),RIGHT(2) / 856686592, 0 / C DATA DIVER(1),DIVER(2) / 873463808, 0 / C DATA LOG10(1),LOG10(2) / 1091781651, 1352628734 / C C MACHINE CONSTANTS FOR THE HARRIS 220 C C DATA SMALL(1),SMALL(2) / '20000000, '00000201 / C DATA LARGE(1),LARGE(2) / '37777777, '37777577 / C DATA RIGHT(1),RIGHT(2) / '20000000, '00000333 / C DATA DIVER(1),DIVER(2) / '20000000, '00000334 / C DATA LOG10(1),LOG10(2) / '23210115, '10237777 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES. C C DATA SMALL(1),SMALL(2) / O402400000000, O000000000000 / C DATA LARGE(1),LARGE(2) / O376777777777, O777777777777 / C DATA RIGHT(1),RIGHT(2) / O604400000000, O000000000000 / C DATA DIVER(1),DIVER(2) / O606400000000, O000000000000 / C DATA LOG10(1),LOG10(2) / O776464202324, O117571775714 / C C MACHINE CONSTANTS FOR THE HP 2100 C THREE WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA SMALL(1), SMALL(2), SMALL(3) / 40000B, 0, 1 / C DATA LARGE(1), LARGE(2), LARGE(3) / 77777B, 177777B, 177776B / C DATA RIGHT(1), RIGHT(2), RIGHT(3) / 40000B, 0, 265B / C DATA DIVER(1), DIVER(2), DIVER(3) / 40000B, 0, 276B / C DATA LOG10(1), LOG10(2), LOG10(3) / 46420B, 46502B, 77777B / C C C MACHINE CONSTANTS FOR THE HP 2100 C FOUR WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA SMALL(1), SMALL(2) / 40000B, 0 / C DATA SMALL(3), SMALL(4) / 0, 1 / C DATA LARGE(1), LARGE(2) / 77777B, 177777B / C DATA LARGE(3), LARGE(4) / 177777B, 177776B / C DATA RIGHT(1), RIGHT(2) / 40000B, 0 / C DATA RIGHT(3), RIGHT(4) / 0, 225B / C DATA DIVER(1), DIVER(2) / 40000B, 0 / C DATA DIVER(3), DIVER(4) / 0, 227B / C DATA LOG10(1), LOG10(2) / 46420B, 46502B / C DATA LOG10(3), LOG10(4) / 76747B, 176377B / C C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND C THE PERKIN ELMER (INTERDATA) 7/32. C C DATA SMALL(1),SMALL(2) / Z00100000, Z00000000 / C DATA LARGE(1),LARGE(2) / Z7FFFFFFF, ZFFFFFFFF / C DATA RIGHT(1),RIGHT(2) / Z33100000, Z00000000 / C DATA DIVER(1),DIVER(2) / Z34100000, Z00000000 / C DATA LOG10(1),LOG10(2) / Z41134413, Z509F79FF / C C MACHINE CONSTANTS FOR THE INTEL 8087, 80287, SEQUENT BALANCE. C ASSUMES INTEGER*2 AS THE DEFAULT FOR TYPE INTEGER. C C DATA SMALL(1), SMALL(2) / 0, 0/ C DATA SMALL(3), SMALL(4) / 0, 16/ C DATA LARGE(1), LARGE(2) / -1, -1/ C DATA LARGE(3), LARGE(4) / -1, 32751/ C DATA RIGHT(1), RIGHT(2) / 0, 0/ C DATA RIGHT(3), RIGHT(4) / 0, 15520/ C DATA DIVER(1), DIVER(2) / 0, 0/ C DATA DIVER(3), DIVER(4) / 0, 15536/ C DATA LOG10(1), LOG10(2) / 31231, 20639/ C DATA LOG10(3), LOG10(4) / 17427, 16339/ C C MACHINE CONSTANTS FOR THE INTEL 8087, 80287, SEQUENT BALANCE. C ASSUMES INTEGER*4 AS THE DEFAULT FOR TYPE INTEGER. C C DATA SMALL(1),SMALL(2) / 0, 1048576 / C DATA LARGE(1),LARGE(2) / -1, 2146435071 / C DATA RIGHT(1),RIGHT(2) / 0, 1017118720 / C DATA DIVER(1),DIVER(2) / 0, 1018167296 / C DATA LOG10(1),LOG10(2) / 1352628735, 1070810131 / C C MACHINE CONSTANTS FOR THE MOTOROLA 68000 SERIES, AT&T 3B SERIES. C ASSUMES INTEGER*2 AS THE DEFAULT FOR TYPE INTEGER. C C DATA SMALL(1), SMALL(2) / 16, 0/ C DATA SMALL(3), SMALL(4) / 0, 0/ C DATA LARGE(1), LARGE(2) / 32751, -1/ C DATA LARGE(3), LARGE(4) / -1, -1/ C DATA RIGHT(1), RIGHT(2) / 15520, 0/ C DATA RIGHT(3), RIGHT(4) / 0, 0/ C DATA DIVER(1), DIVER(2) / 15536, 0/ C DATA DIVER(3), DIVER(4) / 0, 0/ C DATA LOG10(1), LOG10(2) / 16339, 17427/ C DATA LOG10(3), LOG10(4) / 20639, 31231/ C C MACHINE CONSTANTS FOR THE MOTOROLA 68000 SERIES, AT&T 3B SERIES. C ASSUMES INTEGER*4 AS THE DEFAULT FOR TYPE INTEGER. C C DATA SMALL(1),SMALL(2) / 1048576, 0 / C DATA LARGE(1),LARGE(2) / 2146435071, -1 / C DATA RIGHT(1),RIGHT(2) / 1017118720, 0 / C DATA DIVER(1),DIVER(2) / 1018167296, 0 / C DATA LOG10(1),LOG10(2) / 1070810131, 1352628735 / C C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR). C C DATA SMALL(1),SMALL(2) / "033400000000, "000000000000 / C DATA LARGE(1),LARGE(2) / "377777777777, "344777777777 / C DATA RIGHT(1),RIGHT(2) / "113400000000, "000000000000 / C DATA DIVER(1),DIVER(2) / "114400000000, "000000000000 / C DATA LOG10(1),LOG10(2) / "177464202324, "144117571776 / C C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR). C C DATA SMALL(1),SMALL(2) / "000400000000, "000000000000 / C DATA LARGE(1),LARGE(2) / "377777777777, "377777777777 / C DATA RIGHT(1),RIGHT(2) / "103400000000, "000000000000 / C DATA DIVER(1),DIVER(2) / "104400000000, "000000000000 / C DATA LOG10(1),LOG10(2) / "177464202324, "476747767461 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN'S SUPPORTING C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1),SMALL(2) / 8388608, 0 / C DATA LARGE(1),LARGE(2) / 2147483647, -1 / C DATA RIGHT(1),RIGHT(2) / 612368384, 0 / C DATA DIVER(1),DIVER(2) / 620756992, 0 / C DATA LOG10(1),LOG10(2) / 1067065498, -2063872008 / C C DATA SMALL(1),SMALL(2) / O00040000000, O00000000000 / C DATA LARGE(1),LARGE(2) / O17777777777, O37777777777 / C DATA RIGHT(1),RIGHT(2) / O04440000000, O00000000000 / C DATA DIVER(1),DIVER(2) / O04500000000, O00000000000 / C DATA LOG10(1),LOG10(2) / O07746420232, O20476747770 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN'S SUPPORTING C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1),SMALL(2) / 128, 0 / C DATA SMALL(3),SMALL(4) / 0, 0 / C C DATA LARGE(1),LARGE(2) / 32767, -1 / C DATA LARGE(3),LARGE(4) / -1, -1 / C C DATA RIGHT(1),RIGHT(2) / 9344, 0 / C DATA RIGHT(3),RIGHT(4) / 0, 0 / C C DATA DIVER(1),DIVER(2) / 9472, 0 / C DATA DIVER(3),DIVER(4) / 0, 0 / C C DATA LOG10(1),LOG10(2) / 16282, 8346 / C DATA LOG10(3),LOG10(4) / -31493, -12296 / C C DATA SMALL(1),SMALL(2) / O000200, O000000 / C DATA SMALL(3),SMALL(4) / O000000, O000000 / C C DATA LARGE(1),LARGE(2) / O077777, O177777 / C DATA LARGE(3),LARGE(4) / O177777, O177777 / C C DATA RIGHT(1),RIGHT(2) / O022200, O000000 / C DATA RIGHT(3),RIGHT(4) / O000000, O000000 / C C DATA DIVER(1),DIVER(2) / O022400, O000000 / C DATA DIVER(3),DIVER(4) / O000000, O000000 / C C DATA LOG10(1),LOG10(2) / O037632, O020232 / C DATA LOG10(3),LOG10(4) / O102373, O147770 / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. C C DATA SMALL(1),SMALL(2) / O000040000000, O000000000000 / C DATA LARGE(1),LARGE(2) / O377777777777, O777777777777 / C DATA RIGHT(1),RIGHT(2) / O170540000000, O000000000000 / C DATA DIVER(1),DIVER(2) / O170640000000, O000000000000 / C DATA LOG10(1),LOG10(2) / O177746420232, O411757177572 / C C C MACHINE CONSTANTS FOR VAX 11/780 C (EXPRESSED IN INTEGER AND HEXADECIMAL) C ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS*** C *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS*** C DATA SMALL(1), SMALL(2) / 128, 0 / DATA LARGE(1), LARGE(2) / -32769, -1 / DATA RIGHT(1), RIGHT(2) / 9344, 0 / DATA DIVER(1), DIVER(2) / 9472, 0 / DATA LOG10(1), LOG10(2) / 546979738, -805665541 / C C DATA SMALL(1), SMALL(2) / Z00000080, Z00000000 / C DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / C DATA RIGHT(1), RIGHT(2) / Z00002480, Z00000000 / C DATA DIVER(1), DIVER(2) / Z00002500, Z00000000 / C DATA LOG10(1), LOG10(2) / Z209A3F9A, ZCFFA84FB / C C MACHINE CONSTANTS FOR VAX 11/780 (G-FLOATING) C (EXPRESSED IN INTEGER AND HEXADECIMAL) C ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS*** C *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS*** C C DATA SMALL(1), SMALL(2) / 16, 0 / C DATA LARGE(1), LARGE(2) / -32769, -1 / C DATA RIGHT(1), RIGHT(2) / 15552, 0 / C DATA DIVER(1), DIVER(2) / 15568, 0 / C DATA LOG10(1), LOG10(2) / 1142112243, 2046775455 / C C DATA SMALL(1), SMALL(2) / Z00000010, Z00000000 / C DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / C DATA RIGHT(1), RIGHT(2) / Z00003CC0, Z00000000 / C DATA DIVER(1), DIVER(2) / Z00003CD0, Z00000000 / C DATA LOG10(1), LOG10(2) / Z44133FF3, Z79FF509F / C C***FIRST EXECUTABLE STATEMENT D1MACH C D1MACH = DMACH(I) RETURN C END SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY) C C CONSTANT TIMES A VECTOR PLUS A VECTOR. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1),DA INTEGER I,INCX,INCY,IX,IY,M,MP1,N C IF(N.LE.0)RETURN IF (DA .EQ. 0.0D0) RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DY(IY) = DY(IY) + DA*DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,4) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DY(I) = DY(I) + DA*DX(I) 30 CONTINUE IF( N .LT. 4 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,4 DY(I) = DY(I) + DA*DX(I) DY(I + 1) = DY(I + 1) + DA*DX(I + 1) DY(I + 2) = DY(I + 2) + DA*DX(I + 2) DY(I + 3) = DY(I + 3) + DA*DX(I + 3) 50 CONTINUE RETURN END SUBROUTINE DCOPY(N,DX,INCX,DY,INCY) C C COPIES A VECTOR, X, TO A VECTOR, Y. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1) INTEGER I,INCX,INCY,IX,IY,M,MP1,N C IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DY(IY) = DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,7) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DY(I) = DX(I) 30 CONTINUE IF( N .LT. 7 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,7 DY(I) = DX(I) DY(I + 1) = DX(I + 1) DY(I + 2) = DX(I + 2) DY(I + 3) = DX(I + 3) DY(I + 4) = DX(I + 4) DY(I + 5) = DX(I + 5) DY(I + 6) = DX(I + 6) 50 CONTINUE RETURN END SUBROUTINE DCPOSE(NDIM,N,QR,ALPHA,PIVOT,IERR,Y,SUM) C C SUBROUTINE DCPOSE IS A MODIFICATION OF THE ALGOL PROCEDURE C DECOMPOSE IN P. BUSINGER AND G. H. GOLUB, LINEAR LEAST C SQUARES SOLUTIONS BY HOUSEHOLDER TRANSFORMATIONS, C NUMER. MATH. 7 (1965) 269-276. C INTEGER NDIM,N,PIVOT(1) DOUBLE PRECISION QR(NDIM,1),ALPHA(N) INTEGER IERR,I,J,JBAR,K,KP1,NP1 DOUBLE PRECISION BETA,SIGMA,ALPHAK,QRKK,Y(1),SUM(1) DOUBLE PRECISION DDOT IERR=0 NP1=N+1 DO 20 J=1,NP1 SUM(J)=DDOT(N,QR(1,J),1,QR(1,J),1) 20 PIVOT(J)=J DO 500 K=1,N SIGMA=SUM(K) JBAR=K KP1=K+1 DO 40 J=KP1,NP1 IF (SIGMA .GE. SUM(J)) GO TO 40 SIGMA=SUM(J) JBAR=J 40 CONTINUE IF (JBAR .EQ. K) GO TO 70 I=PIVOT(K) PIVOT(K)=PIVOT(JBAR) PIVOT(JBAR)=I SUM(JBAR)=SUM(K) SUM(K)=SIGMA DO 50 I=1,N SIGMA=QR(I,K) QR(I,K)=QR(I,JBAR) QR(I,JBAR)=SIGMA 50 CONTINUE C END OF COLUMN INTERCHANGE. 70 SIGMA=DDOT(N-K+1,QR(K,K),1,QR(K,K),1) IF (SIGMA .NE. 0.0) GO TO 60 IERR=1 RETURN 60 IF (K .EQ. N) GO TO 500 QRKK=QR(K,K) ALPHAK=-SQRT(SIGMA) IF (QRKK .LT. 0.0) ALPHAK=-ALPHAK ALPHA(K)=ALPHAK BETA=1.0/(SIGMA-QRKK*ALPHAK) QR(K,K)=QRKK-ALPHAK DO 80 J=KP1,NP1 80 Y(J)=BETA*DDOT(N-K+1,QR(K,K),1,QR(K,J),1) DO 100 J=KP1,NP1 DO 90 I=K,N QR(I,J)=QR(I,J)-QR(I,K)*Y(J) 90 CONTINUE SUM(J)=SUM(J)-QR(K,J)**2 100 CONTINUE 500 CONTINUE ALPHA(N)=QR(N,N) RETURN END DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY) C C FORMS THE DOT PRODUCT OF TWO VECTORS. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1),DTEMP INTEGER I,INCX,INCY,IX,IY,M,MP1,N C DDOT = 0.0D0 DTEMP = 0.0D0 IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DTEMP = DTEMP + DX(IX)*DY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE DDOT = DTEMP RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DTEMP = DTEMP + DX(I)*DY(I) 30 CONTINUE IF( N .LT. 5 ) GO TO 60 40 MP1 = M + 1 DO 50 I = MP1,N,5 DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) + * DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4) 50 CONTINUE 60 DDOT = DTEMP RETURN END 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 D1MACH . C C 0 IF DIVISION DOES NOT CAUSE OVERFLOW. C C DECLARATION OF INPUT DOUBLE PRECISION XXXX,YYYY DIMENSION XXXX(2),YYYY(2) C C DECLARATION OF OUTPUT INTEGER IERR DOUBLE PRECISION ZZZZ DIMENSION ZZZZ(2) C C DECLARATION OF VARIABLES DOUBLE PRECISION DENOM,XNUM,D1MACH 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)/D1MACH(2) .LT. ABS(DENOM) ) ) THEN ZZZZ(1) = XNUM/DENOM ELSE ZZZZ(1) = D1MACH(2) 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)/D1MACH(2) .LT. ABS(DENOM) ) ) THEN ZZZZ(2) = XNUM/DENOM ELSE ZZZZ(2) = D1MACH(2) IERR =1 END IF RETURN END * DOUBLE PRECISION FUNCTION DNRM2 ( N, DX, INCX) INTEGER I,INCX,J,N,NEXT,NN DOUBLE PRECISION DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX DOUBLE PRECISION ONE,ZERO PARAMETER (ZERO=0.0D0, ONE=1.0D0) C C EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE C INCREMENT INCX . C IF N .LE. 0 RETURN WITH RESULT = 0. C IF N .GE. 1 THEN INCX MUST BE .GE. 1 C C C.L.LAWSON, 1978 JAN 08 C C FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE C HOPEFULLY APPLICABLE TO ALL MACHINES. C CUTLO = MAXIMUM OF DSQRT(U/EPS) OVER ALL KNOWN MACHINES. C CUTHI = MINIMUM OF DSQRT(V) OVER ALL KNOWN MACHINES. C WHERE C EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1. C U = SMALLEST POSITIVE NO. (UNDERFLOW LIMIT) C V = LARGEST NO. (OVERFLOW LIMIT) C C BRIEF OUTLINE OF ALGORITHM.. C C PHASE 1 SCANS ZERO COMPONENTS. C MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO C MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO C MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M C WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX. C C VALUES FOR CUTLO AND CUTHI.. C FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER C DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS.. C CUTLO, S.P. U/EPS = 2**(-102) FOR HONEYWELL. CLOSE SECONDS ARE C UNIVAC AND DEC AT 2**(-103) C THUS CUTLO = 2**(-51) = 4.44089E-16 C CUTHI, S.P. V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC. C THUS CUTHI = 2**(63.5) = 1.30438E19 C CUTLO, D.P. U/EPS = 2**(-67) FOR HONEYWELL AND DEC. C THUS CUTLO = 2**(-33.5) = 8.23181D-11 C CUTHI, D.P. SAME AS S.P. CUTHI = 1.30438D19 C DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / C DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 / DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / C IF(N .GT. 0) GO TO 10 DNRM2 = ZERO GO TO 300 C 10 ASSIGN 30 TO NEXT SUM = ZERO NN = N * INCX C BEGIN MAIN LOOP I = 1 20 GO TO NEXT,(30, 50, 70, 110) 30 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85 ASSIGN 50 TO NEXT XMAX = ZERO C C PHASE 1. SUM IS ZERO C 50 IF( DX(I) .EQ. ZERO) GO TO 200 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85 C C PREPARE FOR PHASE 2. ASSIGN 70 TO NEXT GO TO 105 C C PREPARE FOR PHASE 4. C 100 I = J ASSIGN 110 TO NEXT SUM = (SUM / DX(I)) / DX(I) 105 XMAX = DABS(DX(I)) GO TO 115 C C PHASE 2. SUM IS SMALL. C SCALE TO AVOID DESTRUCTIVE UNDERFLOW. C 70 IF( DABS(DX(I)) .GT. CUTLO ) GO TO 75 C C COMMON CODE FOR PHASES 2 AND 4. C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW. C 110 IF( DABS(DX(I)) .LE. XMAX ) GO TO 115 SUM = ONE + SUM * (XMAX / DX(I))**2 XMAX = DABS(DX(I)) GO TO 200 C 115 SUM = SUM + (DX(I)/XMAX)**2 GO TO 200 C C C PREPARE FOR PHASE 3. C 75 SUM = (SUM * XMAX) * XMAX C C C FOR REAL OR D.P. SET HITEST = CUTHI/N C FOR COMPLEX SET HITEST = CUTHI/(2*N) C 85 HITEST = CUTHI/FLOAT( N ) C C PHASE 3. SUM IS MID-RANGE. NO SCALING. C DO 95 J =I,NN,INCX IF(DABS(DX(J)) .GE. HITEST) GO TO 100 95 SUM = SUM + DX(J)**2 DNRM2 = DSQRT( SUM ) GO TO 300 C 200 CONTINUE I = I + INCX IF ( I .LE. NN ) GO TO 20 C C END OF MAIN LOOP. C C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING. C DNRM2 = XMAX * DSQRT(SUM) 300 CONTINUE RETURN END SUBROUTINE DSCAL(N,DA,DX,INCX) C C SCALES A VECTOR BY A CONSTANT. C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DA,DX(1) INTEGER I,INCX,M,MP1,N,NINCX C IF(N.LE.0)RETURN IF(INCX.EQ.1)GO TO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C NINCX = N*INCX DO 10 I = 1,NINCX,INCX DX(I) = DA*DX(I) 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DX(I) = DA*DX(I) 30 CONTINUE IF( N .LT. 5 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 DX(I) = DA*DX(I) DX(I + 1) = DA*DX(I + 1) DX(I + 2) = DA*DX(I + 2) DX(I + 3) = DA*DX(I + 3) DX(I + 4) = DA*DX(I + 4) 50 CONTINUE RETURN END SUBROUTINE F(X,V) DOUBLE PRECISION X(1),V(1) C C EVALUATE F(X) AND RETURN IN THE VECTOR V . C RETURN END SUBROUTINE FFUNP(N,NUMT,MMAXT,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 MMAXT 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 C C DECLARATION OF INPUT AND OUTPUT: INTEGER N,NUMT,MMAXT,KDEG DOUBLE PRECISION COEF,CL,X,XX,TRM,DTRM,CLX,DXNP1,F,DF DIMENSION NUMT(N),KDEG(N,N+1,MMAXT), $ COEF(N,MMAXT),CL(2,N+1),X(2,N), $ XX(2,N,N+1,MMAXT),TRM(2,N,MMAXT),DTRM(2,N,N+1,MMAXT), $ CLX(2,N),DXNP1(2,N),F(2,N),DF(2,N,N+1) C C DECLARATION OF VARIABLES: INTEGER I,IERR,J,K,L,M,NNNN,NP1 DOUBLE PRECISION TEMP1,TEMP2,XNP1 DIMENSION TEMP1(2),TEMP2(2),XNP1(2) C NP1=N+1 C C GENERATE XNP1, THE PROJECTIVE COORDINATE, AND ITS DERIVATIVES. DO 40 J=1,N CALL MULP(CL(1,J),X(1,J),CLX(1,J)) 40 CONTINUE C DO 60 I=1,2 XNP1(I)=CL(I,NP1) DO 50 J=1,N XNP1(I) = XNP1(I) + CLX(I,J) DXNP1(I,J)=CL(I,J) 50 CONTINUE 60 CONTINUE 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 100 J=1,N DO 100 K=1,NUMT(J) CALL POWP(KDEG(J,NP1,K),XNP1, XX(1,J,NP1,K)) DO 100 L=1,N CALL POWP(KDEG(J, L,K),X(1,L),XX(1,J, L,K)) 100 CONTINUE DO 200 J=1,N DO 200 K=1,NUMT(J) TRM(1,J,K)=COEF(J,K) TRM(2,J,K)=0.0 DO 120 L=1,NP1 CALL MULP(XX(1,J,L,K), TRM(1,J,K),TEMP1) TRM(1,J,K )=TEMP1(1) TRM(2,J,K )=TEMP1(2) 120 CONTINUE 200 CONTINUE DO 300 J=1,N F(1,J)=0.0 F(2,J)=0.0 DO 220 I=1,2 DO 220 K=1,NUMT(J) F(I,J)= F(I,J) + TRM(I,J,K) 220 CONTINUE 300 CONTINUE C C END OF "COMPUTE F" C C BEGIN "COMPUTE DF" C DO 400 J=1,N DO 400 K=1,NUMT(J) DO 400 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,J,M,K)=0.0 DTRM(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,J,M,K)=KDEG(J,M,K)*DTRM(1,J,M,K) DTRM(2,J,M,K)=KDEG(J,M,K)*DTRM(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 320 L=1,NP1 IF (L .EQ. M) GOTO 320 CALL MULP(XX(1,J,L,K),DTRM(1,J,M,K),TEMP1) DTRM(1,J,M,K)=TEMP1(1) DTRM(2,J,M,K)=TEMP1(2) 320 CONTINUE 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,J,M,K)=KDEG(J,M,K)*DTRM(1,J,M,K) DTRM(2,J,M,K)=KDEG(J,M,K)*DTRM(2,J,M,K) END IF END IF 400 CONTINUE DO 600 J=1,N DO 600 M=1,NP1 DF(1,J,M)=0.0 DF(2,J,M)=0.0 DO 420 I=1,2 DO 420 K=1,NUMT(J) DF(I,J,M)= DF(I,J,M) + DTRM(I,J,M,K) 420 CONTINUE 600 CONTINUE 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 700 J=1,N DO 700 K=1,N CALL MULP(DF(1,J,NP1),DXNP1(1,K),TEMP1) DO 700 I=1,2 DF(I,J,K)=DF(I,J,K)+TEMP1(I) 700 CONTINUE RETURN END * SUBROUTINE FIXPDF(N,Y,IFLAG,ARCTOL,EPS,TRACE,A,NDIMA,NFE, $ ARCLEN,YP,YPOLD,QR,ALPHA,TZ,PIVOT,WT,PHI,P,PAR,IPAR) C C SUBROUTINE FIXPDF FINDS A FIXED POINT OR ZERO OF THE C N-DIMENSIONAL VECTOR FUNCTION F(X), OR TRACKS A ZERO CURVE C OF A GENERAL HOMOTOPY MAP RHO(A,LAMBDA,X). FOR THE FIXED C POINT PROBLEM F(X) IS ASSUMED TO BE A C2 MAP OF SOME BALL C INTO ITSELF. THE EQUATION X = F(X) IS SOLVED BY C FOLLOWING THE ZERO CURVE OF THE HOMOTOPY MAP C C LAMBDA*(X - F(X)) + (1 - LAMBDA)*(X - A) , C C STARTING FROM LAMBDA = 0, X = A. THE CURVE IS PARAMETERIZED C BY ARC LENGTH S, AND IS FOLLOWED BY SOLVING THE ORDINARY C DIFFERENTIAL EQUATION D(HOMOTOPY MAP)/DS = 0 FOR C Y(S) = (LAMBDA(S), X(S)). C C FOR THE ZERO FINDING PROBLEM F(X) IS ASSUMED TO BE A C2 MAP C SUCH THAT FOR SOME R > 0, X*F(X) >= 0 WHENEVER NORM(X) = R. C THE EQUATION F(X) = 0 IS SOLVED BY FOLLOWING THE ZERO CURVE C OF THE HOMOTOPY MAP C C LAMBDA*F(X) + (1 - LAMBDA)*(X - A) C C EMANATING FROM LAMBDA = 0, X = A. C C A MUST BE AN INTERIOR POINT OF THE ABOVE MENTIONED BALLS. C C FOR THE CURVE TRACKING PROBLEM RHO(A,LAMBDA,X) IS ASSUMED TO C BE A C2 MAP FROM E**M X [0,1) X E**N INTO E**N, WHICH FOR C ALMOST ALL PARAMETER VECTORS A IN SOME NONEMPTY OPEN SUBSET C OF E**M SATISFIES C C RANK [D RHO(A,LAMBDA,X)/D LAMBDA , D RHO(A,LAMBDA,X)/DX] = N C C FOR ALL POINTS (LAMBDA,X) SUCH THAT RHO(A,LAMBDA,X)=0. IT IS C FURTHER ASSUMED THAT C C RANK [ D RHO(A,0,X0)/DX ] = N . C C WITH A FIXED, THE ZERO CURVE OF RHO(A,LAMBDA,X) EMANATING C FROM LAMBDA = 0, X = X0 IS TRACKED UNTIL LAMBDA = 1 BY C SOLVING THE ORDINARY DIFFERENTIAL EQUATION C D RHO(A,LAMBDA(S),X(S))/DS = 0 FOR Y(S) = (LAMBDA(S), X(S)), C WHERE S IS ARC LENGTH ALONG THE ZERO CURVE. ALSO THE HOMOTOPY C MAP RHO(A,LAMBDA,X) IS ASSUMED TO BE CONSTRUCTED SUCH THAT C C D LAMBDA(0)/DS > 0 . C C THIS CODE IS BASED ON THE ALGORITHM IN L. T. WATSON, A C GLOBALLY CONVERGENT ALGORITHM FOR COMPUTING FIXED POINTS OF C C2 MAPS, APPL. MATH. COMPUT., 5 (1979) 297-311. C C C FOR THE FIXED POINT AND ZERO FINDING PROBLEMS, THE USER C MUST SUPPLY A SUBROUTINE F(X,V) WHICH EVALUATES F(X) AT X C AND RETURNS THE VECTOR F(X) IN V, AND A SUBROUTINE FJAC(X,V,K) C WHICH RETURNS IN V THE KTH COLUMN OF THE JACOBIAN MATRIX OF C F(X) EVALUATED AT X. FOR THE CURVE TRACKING PROBLEM, THE USER MUST C SUPPLY A SUBROUTINE RHOA(V,LAMBDA,X,PAR,IPAR) WHICH GIVEN C (LAMBDA,X) RETURNS A PARAMETER VECTOR A IN V SUCH THAT C RHO(A,LAMBDA,X)=0, AND A SUBROUTINE RHOJAC(A,LAMBDA,X,V,K,PAR,IPAR) C WHICH RETURNS IN V THE KTH COLUMN OF THE N X (N+1) JACOBIAN C MATRIX [D RHO/D LAMBDA, D RHO/DX] EVALUATED AT (A,LAMBDA,X). C FIXPDF DIRECTLY OR INDIRECTLY USES THE SUBROUTINES C STEPS , SINTRP , ROOT , FODE , F (OR RHOA ), C FJAC (OR RHOJAC ), DCPOSE , D1MACH , AND THE BLAS FUNCTIONS C DDOT AND DNRM2 . ONLY D1MACH CONTAINS MACHINE C DEPENDENT CONSTANTS. NO OTHER MODIFICATIONS BY THE USER ARE C REQUIRED. C C ***WARNING: THIS SUBROUTINE IS GENERALLY MORE ROBUST THAN FIXPNF C AND FIXPQF , BUT MAY BE SLOWER THAN THOSE SUBROUTINES BY A C FACTOR OF TWO. C C C ON INPUT: C C N IS THE DIMENSION OF X, F(X), AND RHO(A,LAMBDA,X). C C Y IS AN ARRRAY OF LENGTH N + 1. (Y(2),...,Y(N+1)) = A IS THE C STARTING POINT FOR THE ZERO CURVE FOR THE FIXED POINT AND C ZERO FINDING PROBLEMS. (Y(2),...,Y(N+1)) = X0 FOR THE CURVE C TRACKING PROBLEM. C C IFLAG CAN BE -2, -1, 0, 2, OR 3. IFLAG SHOULD BE 0 ON THE C FIRST CALL TO FIXPDF FOR THE PROBLEM X=F(X), -1 FOR THE C PROBLEM F(X)=0, AND -2 FOR THE PROBLEM RHO(A,LAMBDA,X)=0. C IN CERTAIN SITUATIONS IFLAG IS SET TO 2 OR 3 BY FIXPDF, C AND FIXPDF CAN BE CALLED AGAIN WITHOUT CHANGING IFLAG. C C ARCTOL IS THE LOCAL ERROR ALLOWED THE ODE SOLVER WHEN C FOLLOWING THE ZERO CURVE. IF ARCTOL .LE. 0.0 ON INPUT C IT IS RESET TO .5*DSQRT(EPS). NORMALLY ARCTOL SHOULD C BE CONSIDERABLY LARGER THAN EPS. C C EPS IS THE LOCAL ERROR ALLOWED THE ODE SOLVER WHEN VERY C NEAR THE FIXED POINT(ZERO). EPS IS APPROXIMATELY THE C MIXED ABSOLUTE AND RELATIVE ERROR IN THE COMPUTED FIXED C POINT(ZERO). C C TRACE IS AN INTEGER SPECIFYING THE LOGICAL I/O UNIT FOR C INTERMEDIATE OUTPUT. IF TRACE .GT. 0 THE POINTS COMPUTED ON C THE ZERO CURVE ARE WRITTEN TO I/O UNIT TRACE . C C A(1:NDIMA) CONTAINS THE PARAMETER VECTOR A . FOR THE FIXED POINT C AND ZERO FINDING PROBLEMS, A NEED NOT BE INITIALIZED BY THE C USER, AND IS ASSUMED TO HAVE LENGTH N. FOR THE CURVE C TRACKING PROBLEM, A HAS LENGTH NDIMA AND MUST BE INITIALIZED C BY THE USER. C C NDIMA IS THE DIMENSION OF A, AND IS USED ONLY FOR THE CURVE C TRACKING PROBLEM. C C YP(1:N+1) IS A WORK ARRAY CONTAINING THE CURRENT TANGENT C VECTOR TO THE ZERO CURVE. C C YPOLD(1:N+1) IS A WORK ARRAY CONTAINING THE PREVIOUS TANGENT C VECTOR TO THE ZERO CURVE. C C QR(1:N,1:N+1), ALPHA(1:N), TZ(1:N+1), AND PIVOT(1:N+1) ARE C ALL WORK ARRAYS USED BY FODE TO CALCULATE THE TANGENT C VECTOR YP. C C WT(1:N+1), PHI(1:N+1,1:16), AND P(1:N+1) ARE ALL WORK ARRAYS C USED BY THE ODE SUBROUTINE STEPS . C C PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS, C WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES C RHOA, RHOJAC. C C Y, ARCTOL, EPS, ARCLEN, NFE, AND IFLAG SHOULD ALL BE C VARIABLES IN THE CALLING PROGRAM. C C C ON OUTPUT: C C N AND TRACE ARE UNCHANGED. C C Y(1) = LAMBDA, (Y(2),...,Y(N+1)) = X, AND Y IS AN APPROXIMATE C ZERO OF THE HOMOTOPY MAP. NORMALLY LAMBDA = 1 AND X IS A C FIXED POINT(ZERO) OF F(X). IN ABNORMAL SITUATIONS LAMBDA C MAY ONLY BE NEAR 1 AND X IS NEAR A FIXED POINT(ZERO). C C IFLAG = C -2 CAUSES FIXPDF TO INITIALIZE EVERYTHING FOR THE PROBLEM C RHO(A,LAMBDA,X) = 0 (USE ON FIRST CALL). C C -1 CAUSES FIXPDF TO INITIALIZE EVERYTHING FOR THE PROBLEM C F(X) = 0 (USE ON FIRST CALL). C C 0 CAUSES FIXPDF TO INITIALIZE EVERYTHING FOR THE PROBLEM C X = F(X) (USE ON FIRST CALL). C C 1 NORMAL RETURN. C C 2 SPECIFIED ERROR TOLERANCE CANNOT BE MET. EPS HAS BEEN C INCREASED TO A SUITABLE VALUE. TO CONTINUE, JUST CALL C FIXPDF AGAIN WITHOUT CHANGING ANY PARAMETERS. C C 3 STEPS HAS BEEN CALLED 1000 TIMES. TO CONTINUE, CALL C FIXPDF AGAIN WITHOUT CHANGING ANY PARAMETERS. C C 4 JACOBIAN MATRIX DOES NOT HAVE FULL RANK. THE ALGORITHM C HAS FAILED (THE ZERO CURVE OF THE HOMOTOPY MAP CANNOT BE C FOLLOWED ANY FURTHER). C C 5 EPS (OR ARCTOL ) IS TOO LARGE. THE PROBLEM SHOULD BE C RESTARTED BY CALLING FIXPDF WITH A SMALLER EPS (OR C ARCTOL ) AND IFLAG = 0 (-1, -2). C C 6 I - DF(X) IS NEARLY SINGULAR AT THE FIXED POINT (DF(X) IS C NEARLY SINGULAR AT THE ZERO, OR D RHO(A,LAMBDA,X)/DX IS C NEARLY SINGULAR AT LAMBDA = 1 ). ANSWER MAY NOT BE C ACCURATE. C C 7 ILLEGAL INPUT PARAMETERS, A FATAL ERROR. C C ARCTOL = EPS AFTER A NORMAL RETURN (IFLAG = 1). C C EPS IS UNCHANGED AFTER A NORMAL RETURN (IFLAG = 1). IT IS C INCREASED TO AN APPROPRIATE VALUE ON THE RETURN IFLAG = 2. C C A WILL (NORMALLY) HAVE BEEN MODIFIED. C C NFE IS THE NUMBER OF FUNCTION EVALUATIONS (= NUMBER OF C JACOBIAN EVALUATIONS). C C ARCLEN IS THE LENGTH OF THE PATH FOLLOWED. C C DOUBLE PRECISION AOLD,ARCLEN,ARCTOL,CURSW,CURTOL,EPS, 1 EPSSTP,EPST,H,HOLD,S,S99,SA,SB,SOUT,SQNP1,XOLD,Y1SOUT INTEGER IFLAG,IFLAGC,ITER,IVC,J,JUDY,JW,K,KGI,KOLD, 1 KSTEPS,LCODE,LIMIT,LIMITD,N,NDIMA,NFE,NFEC,NP1,TRACE LOGICAL START,CRASH,ST99 C C ***** ARRAY DECLARATIONS. ***** C C ARRAYS NEEDED BY THE ODE SUBROUTINE STEPS . DOUBLE PRECISION Y(N+1),WT(N+1),PHI(N+1,16),P(N+1),YP(N+1), 1 ALPHAS(12),W(12),G(13),GI(11) INTEGER IV(10) C C ARRAYS NEEDED BY FIXPDF , FODE , AND DCPOSE . DOUBLE PRECISION YPOLD(N+1),A(N),QR(N,N+1),ALPHA(N),TZ(N+1), $ PAR(1) INTEGER PIVOT(N+1),IPAR(1) C C ***** END OF DIMENSIONAL INFORMATION. ***** C SAVE EXTERNAL FODE C C LIMITD IS AN UPPER BOUND ON THE NUMBER OF STEPS. IT MAY BE C CHANGED BY CHANGING THE FOLLOWING PARAMETER STATEMENT: PARAMETER (LIMITD=1000) C C C : : : : : : : : : : : : : : : : : : : : : IF (N .LE. 0 .OR. EPS .LE. 0.0 ) IFLAG=7 IF (IFLAG .GE. -2 .AND. IFLAG .LE. 0) GO TO 10 IF (IFLAG .EQ. 2) GO TO 35 IF (IFLAG .EQ. 3) GO TO 30 C ONLY VALID INPUT FOR IFLAG IS -2, -1, 0, 2, 3. IFLAG=7 RETURN C C ***** INITIALIZATION BLOCK. ***** C 10 ARCLEN=0.0 S=0.0 IF (ARCTOL .LE. 0.0) ARCTOL=.5*SQRT(EPS) NFEC=0 IFLAGC=IFLAG NP1=N+1 SQNP1=SQRT(DBLE(NP1)) C C SWITCH FROM THE TOLERANCE ARCTOL TO THE (FINER) TOLERANCE EPS IF C THE CURVATURE OF ANY COMPONENT OF Y EXCEEDS CURSW. C CURSW=10.0 C ST99=.FALSE. START=.TRUE. CRASH=.FALSE. HOLD=1.0 H=.1 EPSSTP=ARCTOL KSTEPS=0 C SET INITIAL CONDITIONS FOR ORDINARY DIFFERENTIAL EQUATION. YPOLD(1)=1.0 YP(1)=1.0 Y(1)=0.0 DO 20 J=2,NP1 YPOLD(J)=0.0 YP(J)=0.0 20 CONTINUE C LOAD A FOR THE FIXED POINT AND ZERO FINDING PROBLEMS. IF (IFLAGC .GE. -1) THEN DO 23 J=2,NP1 A(J-1)=Y(J) 23 CONTINUE ENDIF 30 LIMIT=LIMITD C C ***** END OF INITIALIZATION BLOCK. ***** C C C ***** MAIN LOOP. ***** C 35 DO 150 ITER=1,LIMIT IF (Y(1) .LT. 0.0) THEN 40 ARCLEN=ARCLEN+S IFLAG=5 RETURN ENDIF 50 IF (S .LE. 7.0*SQNP1) GO TO 80 C ARC LENGTH IS GETTING TOO LONG, THE PROBLEM WILL BE C RESTARTED WITH A DIFFERENT A VECTOR. ARCLEN=ARCLEN+S S=0.0 60 START=.TRUE. CRASH=.FALSE. C COMPUTE A NEW A VECTOR. IF (IFLAGC .EQ. -2) THEN DO 63 JW=1,NDIMA QR(JW,1)=A(JW) 63 CONTINUE CALL RHOA(A,Y(1),Y(2),PAR,IPAR) DO 65 JW=1,NDIMA AOLD=QR(JW,1) C IF NEW AND OLD A DIFFER BY TOO MUCH, TRACKING SHOULD NOT CONTINUE. IF (ABS(A(JW)-AOLD) .GT. 1.0+ABS(AOLD)) THEN ARCLEN=ARCLEN+S IFLAG=5 RETURN ENDIF 65 CONTINUE ELSE CALL F(Y(2),YP) DO 70 JW=1,N AOLD=A(JW) IF (IFLAGC .EQ. -1) THEN A(JW)=Y(1)*YP(JW)/(1.0 - Y(1)) + Y(JW+1) ELSE A(JW)=(Y(JW+1) - Y(1)*YP(JW))/(1.0 - Y(1)) ENDIF C IF NEW AND OLD A DIFFER BY TOO MUCH, TRACKING SHOULD NOT CONTINUE. IF (ABS(A(JW)-AOLD) .GT. 1.0+ABS(AOLD)) THEN ARCLEN=ARCLEN+S IFLAG=5 RETURN ENDIF 70 CONTINUE ENDIF GO TO 100 80 IF (Y(1) .LE. .99 .OR. ST99) GO TO 100 C WHEN LAMBDA REACHES .99, THE PROBLEM WILL BE RESTARTED WITH C A NEW A VECTOR. 90 ST99=.TRUE. EPSSTP=EPS ARCTOL=EPS GO TO 60 C C SET DIFFERENT ERROR TOLERANCE FOR HIGH CURVATURE COMPONENTS OF THE C TRAJECTORY Y(S). 100 CURTOL=CURSW*HOLD EPST=EPS/EPSSTP DO 110 JW=1,NP1 IF (ABS(YP(JW)-YPOLD(JW)) .LE. CURTOL) THEN WT(JW)=(ABS(Y(JW))+1.0) ELSE WT(JW)=(ABS(Y(JW))+1.0)*EPST ENDIF 110 CONTINUE C C TAKE A STEP ALONG THE CURVE. CALL STEPS(FODE,NP1,Y,S,H,EPSSTP,WT,START,HOLD,K,KOLD,CRASH, + PHI,P,YP,ALPHAS,W,G,KSTEPS,XOLD,IVC,IV,KGI,GI, + YPOLD,A,QR,ALPHA,TZ,PIVOT,NFEC,IFLAGC,PAR,IPAR) C PRINT LATEST POINT ON CURVE IF REQUESTED. IF (TRACE .GT. 0) THEN WRITE (TRACE,117) ITER,NFEC,S,Y(1),(Y(JW),JW=2,NP1) 117 FORMAT(/' STEP',I5,3X,'NFE =',I5,3X,'ARC LENGTH =',F9.4,3X, $ 'LAMBDA =',F7.4,5X,'X vector:'/1P,(1X,6E12.4)) ENDIF NFE=NFEC C CHECK IF THE STEP WAS SUCCESSFUL. IF (IFLAGC .EQ. 4) THEN ARCLEN=ARCLEN+S IFLAG=4 RETURN ENDIF 120 IF (CRASH) THEN C RETURN CODE FOR ERROR TOLERANCE TOO SMALL. IFLAG=2 C CHANGE ERROR TOLERANCES. EPS=EPSSTP IF (ARCTOL .LT. EPS) ARCTOL=EPS C CHANGE LIMIT ON NUMBER OF ITERATIONS. LIMIT=LIMIT-ITER RETURN ENDIF C 130 IF (Y(1) .GE. 1.0) THEN IF (ST99) GO TO 160 C C IF LAMBDA .GE. 1.0 BUT THE PROBLEM HAS NOT BEEN RESTARTED C WITH A NEW A VECTOR, BACK UP AND RESTART. C S99=S-.5*HOLD C GET AN APPROXIMATE ZERO Y(S) WITH Y(1)=LAMBDA .LT. 1.0 . 135 CALL SINTRP(S,Y,S99,WT,YP,NP1,KOLD,PHI,IVC,IV,KGI,GI, $ ALPHAS,G,W,XOLD,P) IF (WT(1) .LT. 1.0) GO TO 140 S99=.5*(S-HOLD+S99) GO TO 135 C 140 DO 144 JUDY=1,NP1 Y(JUDY)=WT(JUDY) YPOLD(JUDY)=YP(JUDY) 144 CONTINUE S=S99 GO TO 90 ENDIF C 150 CONTINUE C C ***** END OF MAIN LOOP. ***** C C LAMBDA HAS NOT REACHED 1 IN 1000 STEPS. IFLAG=3 RETURN C C C USE INVERSE INTERPOLATION TO GET THE ANSWER AT LAMBDA = 1.0 . C 160 SA=S-HOLD SB=S LCODE=1 170 CALL ROOT(SOUT,Y1SOUT,SA,SB,EPS,EPS,LCODE) C ROOT FINDS S SUCH THAT Y(1)(S) = LAMBDA = 1 . IF (LCODE .GT. 0) GO TO 190 CALL SINTRP(S,Y,SOUT,WT,YP,NP1,KOLD,PHI,IVC,IV,KGI,GI, $ ALPHAS,G,W,XOLD,P) Y1SOUT=WT(1)-1.0 GO TO 170 190 IFLAG=1 C SET IFLAG = 6 IF ROOT COULD NOT GET LAMBDA = 1.0 . IF (LCODE .GT. 2) IFLAG=6 ARCLEN=ARCLEN+SA C LAMBDA(SA) = 1.0 . CALL SINTRP(S,Y,SA,WT,YP,NP1,KOLD,PHI,IVC,IV,KGI,GI, $ ALPHAS,G,W,XOLD,P) C DO 210 J=1,NP1 210 Y(J)=WT(J) RETURN END SUBROUTINE FIXPDS(N,Y,IFLAG,ARCTOL,EPS,TRACE,A,NDIMA,NFE, $ ARCLEN,YP,YPOLD,QR,LENQR,PIVOT,PP,WORK,WT,PHI,P, $ PAR,IPAR) C C SUBROUTINE FIXPDS FINDS A FIXED POINT OR ZERO OF THE C N-DIMENSIONAL VECTOR FUNCTION F(X), OR TRACKS A ZERO CURVE C OF A GENERAL HOMOTOPY MAP RHO(A,X,LAMBDA). FOR THE FIXED C POINT PROBLEM F(X) IS ASSUMED TO BE A C2 MAP OF SOME BALL C INTO ITSELF. THE EQUATION X = F(X) IS SOLVED BY C FOLLOWING THE ZERO CURVE OF THE HOMOTOPY MAP C C LAMBDA*(X - F(X)) + (1 - LAMBDA)*(X - A) , C C STARTING FROM LAMBDA = 0, X = A. THE CURVE IS PARAMETERIZED C BY ARC LENGTH S, AND IS FOLLOWED BY SOLVING THE ORDINARY C DIFFERENTIAL EQUATION D(HOMOTOPY MAP)/DS = 0 FOR C Y(S) = (X(S), LAMBDA(S)). C C FOR THE ZERO FINDING PROBLEM F(X) IS ASSUMED TO BE A C2 MAP C SUCH THAT FOR SOME R > 0, X*F(X) >= 0 WHENEVER NORM(X) = R. C THE EQUATION F(X) = 0 IS SOLVED BY FOLLOWING THE ZERO CURVE C OF THE HOMOTOPY MAP C C LAMBDA*F(X) + (1 - LAMBDA)*(X - A) C C EMANATING FROM LAMBDA = 0, X = A. C C A MUST BE AN INTERIOR POINT OF THE ABOVE MENTIONED BALLS. C C FOR THE CURVE TRACKING PROBLEM RHO(A,X,LAMBDA) IS ASSUMED TO C BE A C2 MAP FROM E**M X E**N X [0,1) INTO E**N, WHICH FOR C ALMOST ALL PARAMETER VECTORS A IN SOME NONEMPTY OPEN SUBSET C OF E**M SATISFIES C C RANK [D RHO(A,X,LAMBDA)/D LAMBDA , D RHO(A,X,LAMBDA)/DX] = N C C FOR ALL POINTS (X,LAMBDA) SUCH THAT RHO(A,X,LAMBDA)=0. IT IS C FURTHER ASSUMED THAT C C RANK [ D RHO(A,X0,0)/DX ] = N . C C WITH A FIXED, THE ZERO CURVE OF RHO(A,X,LAMBDA) EMANATING C FROM LAMBDA = 0, X = X0 IS TRACKED UNTIL LAMBDA = 1 BY C SOLVING THE ORDINARY DIFFERENTIAL EQUATION C D RHO(A,X(S),LAMBDA(S))/DS = 0 FOR Y(S) = (X(S), LAMBDA(S)), C WHERE S IS ARC LENGTH ALONG THE ZERO CURVE. ALSO THE HOMOTOPY C MAP RHO(A,X,LAMBDA) IS ASSUMED TO BE CONSTRUCTED SUCH THAT C C D LAMBDA(0)/DS > 0 . C C THIS CODE IS BASED ON THE ALGORITHM IN L. T. WATSON, A C GLOBALLY CONVERGENT ALGORITHM FOR COMPUTING FIXED POINTS OF C C2 MAPS, APPL. MATH. COMPUT., 5 (1979) 297-311. C C C FOR THE FIXED POINT AND ZERO FINDING PROBLEMS, THE USER C MUST SUPPLY A SUBROUTINE F(X,V) WHICH EVALUATES F(X) AT X C AND RETURNS THE VECTOR F(X) IN V, AND A SUBROUTINE C FJACS(X,QR,LENQR,PIVOT) WHICH EVALUATES THE (SYMMETRIC) C JACOBIAN MATRIX OF F(X) AT X, AND RETURNS THE SYMMETRIC C JACOBIAN MATRIX IN PACKED SKYLINE STORAGE FORMAT IN QR. LENQR C AND PIVOT DESCRIBE THE DATA STRUCTURE IN QR. FOR THE CURVE C TRACKING PROBLEM, THE USER MUST SUPPLY A SUBROUTINE C RHOA(V,LAMBDA,X,PAR,IPAR) WHICH GIVEN (X,LAMBDA) RETURNS A C PARAMETER VECTOR A IN V SUCH THAT RHO(A,X,LAMBDA)=0, AND A C SUBROUTINE RHOJS(A,LAMBDA,X,QR,LENQR,PIVOT,PP,PAR,IPAR) WHICH C RETURNS IN QR THE SYMMETRIC N X N JACOBIAN MATRIX [D RHO/DX] C EVALUATED AT (A,X,LAMBDA) AND STORED IN PACKED SKYLINE FORMAT, C AND RETURNS IN PP THE VECTOR -(D RHO/D LAMBDA) EVALUATED AT C (A,X,LAMBDA). LENQR AND PIVOT DESCRIBE THE DATA STRUCTURE IN QR. C *** NOTE THE MINUS SIGN IN THE DEFINITION OF PP. *** C C C FUNCTIONS AND SUBROUTINES DIRECTLY OR INDIRECTLY CALLED BY FIXPDS: C D1MACH , F (OR RHOA ), FJACS (OR RHOJS ), FODEDS , GMFADS , C MFACDS , MULTDS , PCGDS , QIMUDS , ROOT , SINTRP , SOLVDS , C STEPDS , AND THE BLAS FUNCTIONS DAXPY , DCOPY , DDOT , DNRM2 , C DSCAL , IDAMAX . ONLY D1MACH CONTAINS MACHINE DEPENDENT C CONSTANTS. NO OTHER MODIFICATIONS BY THE USER ARE REQUIRED. C C ***WARNING: THIS SUBROUTINE IS GENERALLY MORE ROBUST THAN FIXPNS C AND FIXPQS , BUT MAY BE SLOWER THAN THOSE SUBROUTINES BY A C FACTOR OF TWO. C C C ON INPUT: C C N IS THE DIMENSION OF X, F(X), AND RHO(A,X,LAMBDA). C C Y IS AN ARRRAY OF LENGTH N + 1. (Y(1),...,Y(N)) = A IS THE C STARTING POINT FOR THE ZERO CURVE FOR THE FIXED POINT AND C ZERO FINDING PROBLEMS. (Y(1),...,Y(N)) = X0 FOR THE CURVE C TRACKING PROBLEM. C C IFLAG CAN BE -2, -1, 0, 2, OR 3. IFLAG SHOULD BE 0 ON THE C FIRST CALL TO FIXPDS FOR THE PROBLEM X=F(X), -1 FOR THE C PROBLEM F(X)=0, AND -2 FOR THE PROBLEM RHO(A,X,LAMBDA)=0. C IN CERTAIN SITUATIONS IFLAG IS SET TO 2 OR 3 BY FIXPDS, C AND FIXPDS CAN BE CALLED AGAIN WITHOUT CHANGING IFLAG. C C ARCTOL IS THE LOCAL ERROR ALLOWED THE ODE SOLVER WHEN C FOLLOWING THE ZERO CURVE. IF ARCTOL .LE. 0.0 ON INPUT C IT IS RESET TO .5*DSQRT(EPS). NORMALLY ARCTOL SHOULD C BE CONSIDERABLY LARGER THAN EPS. C C EPS IS THE LOCAL ERROR ALLOWED THE ODE SOLVER WHEN VERY C NEAR THE FIXED POINT(ZERO). EPS IS APPROXIMATELY THE C MIXED ABSOLUTE AND RELATIVE ERROR IN THE COMPUTED FIXED C POINT(ZERO). C C TRACE IS AN INTEGER SPECIFYING THE LOGICAL I/O UNIT FOR C INTERMEDIATE OUTPUT. IF TRACE .GT. 0 THE POINTS COMPUTED ON C THE ZERO CURVE ARE WRITTEN TO I/O UNIT TRACE . C C A(1:NDIMA) CONTAINS THE PARAMETER VECTOR A . FOR THE FIXED POINT C AND ZERO FINDING PROBLEMS, A NEED NOT BE INITIALIZED BY THE C USER, AND IS ASSUMED TO HAVE LENGTH N. FOR THE CURVE C TRACKING PROBLEM, A HAS LENGTH NDIMA AND MUST BE INITIALIZED C BY THE USER. C C NDIMA IS THE DIMENSION OF A , AND IS USED ONLY FOR THE CURVE C TRACKING PROBLEM. C C YP(1:N+1) IS A WORK ARRAY CONTAINING THE CURRENT TANGENT C VECTOR TO THE ZERO CURVE. C C YPOLD(1:N+1) IS A WORK ARRAY CONTAINING THE PREVIOUS TANGENT C VECTOR TO THE ZERO CURVE. C C QR(1:LENQR) IS A WORK ARRAY CONTAINING THE (SYMMETRIC) JACOBIAN C MATRIX WITH RESPECT TO X, IN THE PACKED SKYLINE STORAGE FORMAT. C C LENQR IS THE DIMENSION OF QR . C C PIVOT(1:N+2), PP(1:N), AND WORK(1:6*(N+1)+LENQR) ARE ALL WORK C ARRAYS USED BY FODEDS TO CALCULATE THE TANGENT VECTOR YP. C C WT(1:N+1), PHI(1:N+1,1:16), AND P(1:N+1) ARE ALL WORK ARRAYS C USED BY THE ODE SUBROUTINE STEPDS . C C PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS, C WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES C RHOA, RHOJS. C C Y, ARCTOL, EPS, ARCLEN, NFE, AND IFLAG SHOULD ALL BE C VARIABLES IN THE CALLING PROGRAM. C C C ON OUTPUT: C C N AND TRACE ARE UNCHANGED. C C (Y(1),...,Y(N)) = X, Y(N+1) = LAMBDA, AND Y IS AN APPROXIMATE C ZERO OF THE HOMOTOPY MAP. NORMALLY LAMBDA = 1 AND X IS A C FIXED POINT(ZERO) OF F(X). IN ABNORMAL SITUATIONS LAMBDA C MAY ONLY BE NEAR 1 AND X IS NEAR A FIXED POINT(ZERO). C C IFLAG = C -2 CAUSES FIXPDS TO INITIALIZE EVERYTHING FOR THE PROBLEM C RHO(A,X,LAMBDA) = 0 (USE ON FIRST CALL). C C -1 CAUSES FIXPDS TO INITIALIZE EVERYTHING FOR THE PROBLEM C F(X) = 0 (USE ON FIRST CALL). C C 0 CAUSES FIXPDS TO INITIALIZE EVERYTHING FOR THE PROBLEM C X = F(X) (USE ON FIRST CALL). C C 1 NORMAL RETURN. C C 2 SPECIFIED ERROR TOLERANCE CANNOT BE MET. EPS HAS BEEN C INCREASED TO A SUITABLE VALUE. TO CONTINUE, JUST CALL C FIXPDS AGAIN WITHOUT CHANGING ANY PARAMETERS. C C 3 STEPDS HAS BEEN CALLED 1000 TIMES. TO CONTINUE, CALL C FIXPDS AGAIN WITHOUT CHANGING ANY PARAMETERS. C C 4 JACOBIAN MATRIX DOES NOT HAVE FULL RANK AND/OR THE CONJUGATE C GRADIENT ITERATION FOR THE KERNEL OF THE JACOBIAN MATRIX C FAILED TO CONVERGE. THE ALGORITHM HAS FAILED (THE ZERO C CURVE OF THE HOMOTOPY MAP CANNOT BE FOLLOWED ANY FURTHER). C C 5 EPS (OR ARCTOL ) IS TOO LARGE. THE PROBLEM SHOULD BE C RESTARTED BY CALLING FIXPDS WITH A SMALLER EPS (OR C ARCTOL ) AND IFLAG = 0 (-1, -2). C C 6 I - DF(X) IS NEARLY SINGULAR AT THE FIXED POINT (DF(X) IS C NEARLY SINGULAR AT THE ZERO, OR D RHO(A,X,LAMBDA)/DX IS C NEARLY SINGULAR AT LAMBDA = 1 ). ANSWER MAY NOT BE C ACCURATE. C C 7 ILLEGAL INPUT PARAMETERS, A FATAL ERROR. C C ARCTOL = EPS AFTER A NORMAL RETURN (IFLAG = 1). C C EPS IS UNCHANGED AFTER A NORMAL RETURN (IFLAG = 1). IT IS C INCREASED TO AN APPROPRIATE VALUE ON THE RETURN IFLAG = 2. C C A WILL (NORMALLY) HAVE BEEN MODIFIED. C C NFE IS THE NUMBER OF FUNCTION EVALUATIONS (= NUMBER OF C JACOBIAN EVALUATIONS). C C ARCLEN IS THE LENGTH OF THE PATH FOLLOWED. C C DOUBLE PRECISION AOLD,ARCLEN,ARCTOL,CURSW,CURTOL,EPS, 1 EPSSTP,EPST,H,HOLD,S,S99,SA,SB,SOUT,SQNP1,XOLD,Y1SOUT INTEGER IFLAG,IFLAGC,ITER,IVC,J,JW,K,KGI,KOLD, 1 KSTEPS,LCODE,LENQR,LIMIT,LIMITD,N,NDIMA,NFE,NFEC,NP1,TRACE LOGICAL START,CRASH,ST99 C C ***** ARRAY DECLARATIONS. ***** C C ARRAYS NEEDED BY THE ODE SUBROUTINE STEPDS . DOUBLE PRECISION Y(N+1),WT(N+1),PHI(N+1,16),P(N+1),YP(N+1), 1 ALPHAS(12),W(12),G(13),GI(11) INTEGER IV(10) C C ARRAYS NEEDED BY FIXPDS , FODEDS , AND PCGDS . DOUBLE PRECISION YPOLD(N+1),A(N),QR(LENQR),PP(N), 1 WORK(6*(N+1)+LENQR),PAR(1) INTEGER PIVOT(N+2),IPAR(1) C C ***** END OF DIMENSIONAL INFORMATION. ***** C SAVE EXTERNAL FODEDS C C LIMITD IS AN UPPER BOUND ON THE NUMBER OF STEPS. IT MAY BE C CHANGED BY CHANGING THE FOLLOWING PARAMETER STATEMENT: PARAMETER (LIMITD=1000) C C C : : : : : : : : : : : : : : : : : : : : : IF (N .LE. 0 .OR. EPS .LE. 0.0 ) IFLAG=7 IF (IFLAG .GE. -2 .AND. IFLAG .LE. 0) GO TO 10 IF (IFLAG .EQ. 2) GO TO 35 IF (IFLAG .EQ. 3) GO TO 30 C ONLY VALID INPUT FOR IFLAG IS -2, -1, 0, 2, 3. IFLAG=7 RETURN C C ***** INITIALIZATION BLOCK. ***** C 10 ARCLEN=0.0 S=0.0 IF (ARCTOL .LE. 0.0) ARCTOL=.5*SQRT(EPS) NFEC=0 IFLAGC=IFLAG NP1=N+1 SQNP1=SQRT(DBLE(NP1)) C C SWITCH FROM THE TOLERANCE ARCTOL TO THE (FINER) TOLERANCE EPS IF C THE CURVATURE OF ANY COMPONENT OF Y EXCEEDS CURSW. C CURSW=10.0 C ST99=.FALSE. START=.TRUE. CRASH=.FALSE. HOLD=1.0 H=.1 EPSSTP=ARCTOL KSTEPS=0 C SET INITIAL CONDITIONS FOR ORDINARY DIFFERENTIAL EQUATION. YPOLD(NP1)=1.0 YP(NP1)=1.0 Y(NP1)=0.0 WORK(2*NP1)=0.0 WORK(3*NP1)=0.0 DO 20 J=1,N YPOLD(J)=0.0 YP(J)=0.0 WORK(NP1+J)=0.0 WORK(2*NP1+J)=0.0 20 CONTINUE C LOAD A FOR THE FIXED POINT AND ZERO FINDING PROBLEMS. IF (IFLAGC .GE. -1) THEN CALL DCOPY(N,Y,1,A,1) ENDIF 30 LIMIT=LIMITD C C ***** END OF INITIALIZATION BLOCK. ***** C C C ***** MAIN LOOP. ***** C 35 DO 150 ITER=1,LIMIT IF (Y(NP1) .LT. 0.0) THEN 40 ARCLEN=ARCLEN+S IFLAG=5 RETURN ENDIF 50 IF (S .LE. 7.0*SQNP1) GO TO 80 C ARC LENGTH IS GETTING TOO LONG, THE PROBLEM WILL BE C RESTARTED WITH A DIFFERENT A VECTOR. ARCLEN=ARCLEN+S S=0.0 60 START=.TRUE. CRASH=.FALSE. C COMPUTE A NEW A VECTOR. IF (IFLAGC .EQ. -2) THEN DO 63 JW=1,NDIMA QR(JW)=A(JW) 63 CONTINUE CALL RHOA(A,Y(NP1),Y,PAR,IPAR) DO 65 JW=1,NDIMA AOLD=QR(JW) C IF NEW AND OLD A DIFFER BY TOO MUCH, TRACKING SHOULD NOT CONTINUE. IF (ABS(A(JW)-AOLD) .GT. 1.0+ABS(AOLD)) THEN ARCLEN=ARCLEN+S IFLAG=5 RETURN ENDIF 65 CONTINUE ELSE CALL F(Y,YP) DO 70 JW=1,N AOLD=A(JW) IF (IFLAGC .EQ. -1) THEN A(JW)=Y(NP1)*YP(JW)/(1.0 - Y(NP1)) + Y(JW) ELSE A(JW)=(Y(JW) - Y(NP1)*YP(JW))/(1.0 - Y(NP1)) ENDIF C IF NEW AND OLD A DIFFER BY TOO MUCH, TRACKING SHOULD NOT CONTINUE. IF (ABS(A(JW)-AOLD) .GT. 1.0+ABS(AOLD)) THEN ARCLEN=ARCLEN+S IFLAG=5 RETURN ENDIF 70 CONTINUE ENDIF GO TO 100 80 IF (Y(NP1) .LE. .99 .OR. ST99) GO TO 100 C WHEN LAMBDA REACHES .99, THE PROBLEM WILL BE RESTARTED WITH C A NEW A VECTOR. 90 ST99=.TRUE. EPSSTP=EPS ARCTOL=EPS GO TO 60 C C SET DIFFERENT ERROR TOLERANCE FOR HIGH CURVATURE COMPONENTS OF THE C TRAJECTORY Y(S). 100 CURTOL=CURSW*HOLD EPST=EPS/EPSSTP DO 110 JW=1,NP1 IF (ABS(YP(JW)-YPOLD(JW)) .LE. CURTOL) THEN WT(JW)=(ABS(Y(JW))+1.0) ELSE WT(JW)=(ABS(Y(JW))+1.0)*EPST ENDIF 110 CONTINUE C C TAKE A STEP ALONG THE CURVE. CALL STEPDS(FODEDS,NP1,Y,S,H,EPSSTP,WT,START,HOLD,K,KOLD,CRASH, + PHI,P,YP,ALPHAS,W,G,KSTEPS,XOLD,IVC,IV,KGI,GI, + YPOLD,A,QR,LENQR,PIVOT,PP,WORK,NFEC,IFLAGC,PAR,IPAR) C PRINT LATEST POINT ON CURVE IF REQUESTED. IF (TRACE .GT. 0) THEN WRITE (TRACE,117) ITER,NFEC,S,Y(NP1),(Y(JW),JW=1,N) 117 FORMAT(/' STEP',I5,3X,'NFE =',I5,3X,'ARC LENGTH =',F9.4,3X, $ 'LAMBDA =',F7.4,5X,'X vector:'/1P,(1X,6E12.4)) ENDIF NFE=NFEC C CHECK IF THE STEP WAS SUCCESSFUL. IF (IFLAGC .EQ. 4) THEN ARCLEN=ARCLEN+S IFLAG=4 RETURN ENDIF 120 IF (CRASH) THEN C RETURN CODE FOR ERROR TOLERANCE TOO SMALL. IFLAG=2 C CHANGE ERROR TOLERANCES. EPS=EPSSTP IF (ARCTOL .LT. EPS) ARCTOL=EPS C CHANGE LIMIT ON NUMBER OF ITERATIONS. LIMIT=LIMIT-ITER RETURN ENDIF C 130 IF (Y(NP1) .GE. 1.0) THEN IF (ST99) GO TO 160 C C IF LAMBDA .GE. 1.0 BUT THE PROBLEM HAS NOT BEEN RESTARTED C WITH A NEW A VECTOR, BACK UP AND RESTART. C S99=S-.5*HOLD C GET AN APPROXIMATE ZERO Y(S) WITH Y(NP1)=LAMBDA .LT. 1.0 . 135 CALL SINTRP(S,Y,S99,WT,YP,NP1,KOLD,PHI,IVC,IV,KGI,GI, $ ALPHAS,G,W,XOLD,P) IF (WT(NP1) .LT. 1.0) GO TO 140 S99=.5*(S-HOLD+S99) GO TO 135 C 140 CALL DCOPY(NP1,WT,1,Y,1) CALL DCOPY(NP1,YP,1,YPOLD,1) S=S99 GO TO 90 ENDIF C 150 CONTINUE C C ***** END OF MAIN LOOP. ***** C C LAMBDA HAS NOT REACHED 1 IN 1000 STEPS. IFLAG=3 RETURN C C C USE INVERSE INTERPOLATION TO GET THE ANSWER AT LAMBDA = 1.0 . C 160 SA=S-HOLD SB=S LCODE=1 170 CALL ROOT(SOUT,Y1SOUT,SA,SB,EPS,EPS,LCODE) C ROOT FINDS S SUCH THAT Y(NP1)(S) = LAMBDA = 1 . IF (LCODE .GT. 0) GO TO 190 CALL SINTRP(S,Y,SOUT,WT,YP,NP1,KOLD,PHI,IVC,IV,KGI,GI, $ ALPHAS,G,W,XOLD,P) Y1SOUT=WT(NP1)-1.0 GO TO 170 190 IFLAG=1 C SET IFLAG = 6 IF ROOT COULD NOT GET LAMBDA = 1.0 . IF (LCODE .GT. 2) IFLAG=6 ARCLEN=ARCLEN+SA C LAMBDA(SA) = 1.0 . CALL SINTRP(S,Y,SA,WT,YP,NP1,KOLD,PHI,IVC,IV,KGI,GI, $ ALPHAS,G,W,XOLD,P) C CALL DCOPY(NP1,WT,1,Y,1) C RETURN END SUBROUTINE FIXPNF(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,A,NFE, $ ARCLEN,YP,YOLD,YPOLD,QR,ALPHA,TZ,PIVOT,W,WP,Z0,Z1,SSPAR, $ PAR,IPAR) C C SUBROUTINE FIXPNF FINDS A FIXED POINT OR ZERO OF THE C N-DIMENSIONAL VECTOR FUNCTION F(X), OR TRACKS A ZERO CURVE C OF A GENERAL HOMOTOPY MAP RHO(A,LAMBDA,X). FOR THE FIXED C POINT PROBLEM F(X) IS ASSUMED TO BE A C2 MAP OF SOME BALL C INTO ITSELF. THE EQUATION X = F(X) IS SOLVED BY C FOLLOWING THE ZERO CURVE OF THE HOMOTOPY MAP C C LAMBDA*(X - F(X)) + (1 - LAMBDA)*(X - A) , C C STARTING FROM LAMBDA = 0, X = A. THE CURVE IS PARAMETERIZED C BY ARC LENGTH S, AND IS FOLLOWED BY SOLVING THE ORDINARY C DIFFERENTIAL EQUATION D(HOMOTOPY MAP)/DS = 0 FOR C Y(S) = (LAMBDA(S), X(S)) USING A HERMITE CUBIC PREDICTOR AND A C CORRECTOR WHICH RETURNS TO THE ZERO CURVE ALONG THE FLOW NORMAL C TO THE DAVIDENKO FLOW (WHICH CONSISTS OF THE INTEGRAL CURVES OF C D(HOMOTOPY MAP)/DS ). C C FOR THE ZERO FINDING PROBLEM F(X) IS ASSUMED TO BE A C2 MAP C SUCH THAT FOR SOME R > 0, X*F(X) >= 0 WHENEVER NORM(X) = R. C THE EQUATION F(X) = 0 IS SOLVED BY FOLLOWING THE ZERO CURVE C OF THE HOMOTOPY MAP C C LAMBDA*F(X) + (1 - LAMBDA)*(X - A) C C EMANATING FROM LAMBDA = 0, X = A. C C A MUST BE AN INTERIOR POINT OF THE ABOVE MENTIONED BALLS. C C FOR THE CURVE TRACKING PROBLEM RHO(A,LAMBDA,X) IS ASSUMED TO C BE A C2 MAP FROM E**M X [0,1) X E**N INTO E**N, WHICH FOR C ALMOST ALL PARAMETER VECTORS A IN SOME NONEMPTY OPEN SUBSET C OF E**M SATISFIES C C RANK [D RHO(A,LAMBDA,X)/D LAMBDA , D RHO(A,LAMBDA,X)/DX] = N C C FOR ALL POINTS (LAMBDA,X) SUCH THAT RHO(A,LAMBDA,X)=0. IT IS C FURTHER ASSUMED THAT C C RANK [ D RHO(A,0,X0)/DX ] = N . C C WITH A FIXED, THE ZERO CURVE OF RHO(A,LAMBDA,X) EMANATING C FROM LAMBDA = 0, X = X0 IS TRACKED UNTIL LAMBDA = 1 BY C SOLVING THE ORDINARY DIFFERENTIAL EQUATION C D RHO(A,LAMBDA(S),X(S))/DS = 0 FOR Y(S) = (LAMBDA(S), X(S)), C WHERE S IS ARC LENGTH ALONG THE ZERO CURVE. ALSO THE HOMOTOPY C MAP RHO(A,LAMBDA,X) IS ASSUMED TO BE CONSTRUCTED SUCH THAT C C D LAMBDA(0)/DS > 0 . C C C FOR THE FIXED POINT AND ZERO FINDING PROBLEMS, THE USER MUST SUPPLY C A SUBROUTINE F(X,V) WHICH EVALUATES F(X) AT X AND RETURNS THE C VECTOR F(X) IN V, AND A SUBROUTINE FJAC(X,V,K) WHICH RETURNS IN V C THE KTH COLUMN OF THE JACOBIAN MATRIX OF F(X) EVALUATED AT X. FOR C THE CURVE TRACKING PROBLEM, THE USER MUST SUPPLY A SUBROUTINE C RHO(A,LAMBDA,X,V,PAR,IPAR) WHICH EVALUATES THE HOMOTOPY MAP RHO AT C (A,LAMBDA,X) AND RETURNS THE VECTOR RHO(A,LAMBDA,X) IN V, AND A C SUBROUTINE RHOJAC(A,LAMBDA,X,V,K,PAR,IPAR) WHICH RETURNS IN V THE KT C COLUMN OF THE N X (N+1) JACOBIAN MATRIX [D RHO/D LAMBDA, D RHO/DX] C EVALUATED AT (A,LAMBDA,X). FIXPNF DIRECTLY OR INDIRECTLY USES C THE SUBROUTINES STEPNF , TANGNF , ROOTNF , ROOT , F (OR RHO ), C FJAC (OR RHOJAC ), D1MACH , AND THE BLAS FUNCTIONS DDOT AND C DNRM2 . ONLY D1MACH CONTAINS MACHINE DEPENDENT CONSTANTS. C NO OTHER MODIFICATIONS BY THE USER ARE REQUIRED. C C C ON INPUT: C C N IS THE DIMENSION OF X, F(X), AND RHO(A,LAMBDA,X). C C Y IS AN ARRRAY OF LENGTH N + 1. (Y(2),...,Y(N+1)) = A IS THE C STARTING POINT FOR THE ZERO CURVE FOR THE FIXED POINT AND C ZERO FINDING PROBLEMS. (Y(2),...,Y(N+1)) = X0 FOR THE CURVE C TRACKING PROBLEM. C C IFLAG CAN BE -2, -1, 0, 2, OR 3. IFLAG SHOULD BE 0 ON THE C FIRST CALL TO FIXPNF FOR THE PROBLEM X=F(X), -1 FOR THE C PROBLEM F(X)=0, AND -2 FOR THE PROBLEM RHO(A,LAMBDA,X)=0. C IN CERTAIN SITUATIONS IFLAG IS SET TO 2 OR 3 BY FIXPNF, C AND FIXPNF CAN BE CALLED AGAIN WITHOUT CHANGING IFLAG. C C ARCRE , ARCAE ARE THE RELATIVE AND ABSOLUTE ERRORS, RESPECTIVELY, C ALLOWED THE NORMAL FLOW ITERATION ALONG THE ZERO CURVE. IF C ARC?E .LE. 0.0 ON INPUT IT IS RESET TO .5*SQRT(ANS?E) . C NORMALLY ARC?E SHOULD BE CONSIDERABLY LARGER THAN ANS?E . C C ANSRE , ANSAE ARE THE RELATIVE AND ABSOLUTE ERROR VALUES USED FOR C THE ANSWER AT LAMBDA = 1. THE ACCEPTED ANSWER Y = (LAMBDA, X) C SATISFIES C C |Y(1) - 1| .LE. ANSRE + ANSAE .AND. C C ||Z|| .LE. ANSRE*||X|| + ANSAE WHERE C C (.,Z) IS THE NEWTON STEP TO Y. C C TRACE IS AN INTEGER SPECIFYING THE LOGICAL I/O UNIT FOR C INTERMEDIATE OUTPUT. IF TRACE .GT. 0 THE POINTS COMPUTED ON C THE ZERO CURVE ARE WRITTEN TO I/O UNIT TRACE . C C A(1:*) CONTAINS THE PARAMETER VECTOR A . FOR THE FIXED POINT C AND ZERO FINDING PROBLEMS, A NEED NOT BE INITIALIZED BY THE C USER, AND IS ASSUMED TO HAVE LENGTH N. FOR THE CURVE C TRACKING PROBLEM, A MUST BE INITIALIZED BY THE USER. C C YP(1:N+1) IS A WORK ARRAY CONTAINING THE TANGENT VECTOR TO C THE ZERO CURVE AT THE CURRENT POINT Y . C C YOLD(1:N+1) IS A WORK ARRAY CONTAINING THE PREVIOUS POINT FOUND C ON THE ZERO CURVE. C C YPOLD(1:N+1) IS A WORK ARRAY CONTAINING THE TANGENT VECTOR TO C THE ZERO CURVE AT YOLD . C C QR(1:N,1:N+2), ALPHA(1:N), TZ(1:N+1), PIVOT(1:N+1) , W(1:N+1) , C WP(1:N+1) , Z0(1:N+1) , Z1(1:N+1) ARE ALL WORK ARRAYS USED BY C STEPNF TO CALCULATE THE TANGENT VECTORS AND NEWTON STEPS. C C SSPAR(1:8) = (LIDEAL, RIDEAL, DIDEAL, HMIN, HMAX, BMIN, BMAX, P) IS C A VECTOR OF PARAMETERS USED FOR THE OPTIMAL STEP SIZE ESTIMATION. C IF SSPAR(J) .LE. 0.0 ON INPUT, IT IS RESET TO A DEFAULT VALUE C BY FIXPNF . OTHERWISE THE INPUT VALUE OF SSPAR(J) IS USED. C SEE THE COMMENTS BELOW AND IN STEPNF FOR MORE INFORMATION ABOUT C THESE CONSTANTS. C C PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS, C WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES C RHO, RHOJAC. C C C ON OUTPUT: C C N , TRACE , A ARE UNCHANGED. C C Y(1) = LAMBDA, (Y(2),...,Y(N+1)) = X, AND Y IS AN APPROXIMATE C ZERO OF THE HOMOTOPY MAP. NORMALLY LAMBDA = 1 AND X IS A C FIXED POINT(ZERO) OF F(X). IN ABNORMAL SITUATIONS LAMBDA C MAY ONLY BE NEAR 1 AND X IS NEAR A FIXED POINT(ZERO). C C IFLAG = C -2 CAUSES FIXPNF TO INITIALIZE EVERYTHING FOR THE PROBLEM C RHO(A,LAMBDA,X) = 0 (USE ON FIRST CALL). C C -1 CAUSES FIXPNF TO INITIALIZE EVERYTHING FOR THE PROBLEM C F(X) = 0 (USE ON FIRST CALL). C C 0 CAUSES FIXPNF TO INITIALIZE EVERYTHING FOR THE PROBLEM C X = F(X) (USE ON FIRST CALL). C C 1 NORMAL RETURN. C C 2 SPECIFIED ERROR TOLERANCE CANNOT BE MET. SOME OR ALL OF C ARCRE , ARCAE , ANSRE , ANSAE HAVE BEEN INCREASED TO C SUITABLE VALUES. TO CONTINUE, JUST CALL FIXPNF AGAIN C WITHOUT CHANGING ANY PARAMETERS. C C 3 STEPNF HAS BEEN CALLED 1000 TIMES. TO CONTINUE, CALL C FIXPNF AGAIN WITHOUT CHANGING ANY PARAMETERS. C C 4 JACOBIAN MATRIX DOES NOT HAVE FULL RANK. THE ALGORITHM C HAS FAILED (THE ZERO CURVE OF THE HOMOTOPY MAP CANNOT BE C FOLLOWED ANY FURTHER). C C 5 THE TRACKING ALGORITHM HAS LOST THE ZERO CURVE OF THE C HOMOTOPY MAP AND IS NOT MAKING PROGRESS. THE ERROR TOLERANCES C ARC?E AND ANS?E WERE TOO LENIENT. THE PROBLEM SHOULD BE C RESTARTED BY CALLING FIXPNF WITH SMALLER ERROR TOLERANCES C AND IFLAG = 0 (-1, -2). C C 6 THE NORMAL FLOW NEWTON ITERATION IN STEPNF OR ROOTNF C FAILED TO CONVERGE. THE ERROR TOLERANCES ANS?E MAY BE TOO C STRINGENT. C C 7 ILLEGAL INPUT PARAMETERS, A FATAL ERROR. C C ARCRE , ARCAE , ANSRE , ANSAE ARE UNCHANGED AFTER A NORMAL RETURN C (IFLAG = 1). THEY ARE INCREASED TO APPROPRIATE VALUES ON THE C RETURN IFLAG = 2 . C C NFE IS THE NUMBER OF FUNCTION EVALUATIONS (= NUMBER OF C JACOBIAN EVALUATIONS). C C ARCLEN IS THE LENGTH OF THE PATH FOLLOWED. C C C DOUBLE PRECISION ABSERR,ANSAE,ANSRE,ARCAE,ARCLEN,ARCRE, 1 CURSW,CURTOL,D1MACH,DNRM2,H,HOLD,RELERR,S INTEGER IFLAG,IFLAGC,ITER,JW,LIMIT,LIMITD,N,NC,NFE,NFEC,NP1, 1 TRACE LOGICAL CRASH,POLSYS,START C C ***** ARRAY DECLARATIONS. ***** C DOUBLE PRECISION Y(N+1),YP(N+1),YOLD(N+1),YPOLD(N+1),A(N), $ QR(N,N+2),ALPHA(N),TZ(N+1),W(N+1),WP(N+1),Z0(N+1), $ Z1(N+1),SSPAR(8),PAR(1) INTEGER PIVOT(N+1),IPAR(1) C C ***** END OF DIMENSIONAL INFORMATION. ***** C SAVE C C LIMITD IS AN UPPER BOUND ON THE NUMBER OF STEPS. IT MAY BE C CHANGED BY CHANGING THE FOLLOWING PARAMETER STATEMENT: PARAMETER (LIMITD=1000) C C SWITCH FROM THE TOLERANCE ARC?E TO THE (FINER) TOLERANCE ANS?E IF C THE CURVATURE OF ANY COMPONENT OF Y EXCEEDS CURSW. PARAMETER (CURSW=10.0) C C C C : : : : : : : : : : : : : : : : : : : : : : : : C SET LOGICAL SWITCH TO REFLECT ENTRY POINT. POLSYS=.FALSE. GO TO 11 ENTRY POLYNF(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,A,NFE, $ ARCLEN,YP,YOLD,YPOLD,QR,ALPHA,TZ,PIVOT,W,WP,Z0,Z1,SSPAR, $ PAR,IPAR) POLSYS=.TRUE. 11 CONTINUE C IF (N .LE. 0 .OR. ANSRE .LE. 0.0 .OR. ANSAE .LT. 0.0) $ IFLAG=7 IF (IFLAG .GE. -2 .AND. IFLAG .LE. 0) GO TO 20 IF (IFLAG .EQ. 2) GO TO 120 IF (IFLAG .EQ. 3) GO TO 90 C ONLY VALID INPUT FOR IFLAG IS -2, -1, 0, 2, 3. IFLAG=7 RETURN C C ***** INITIALIZATION BLOCK. ***** C 20 ARCLEN=0.0 IF (ARCRE .LE. 0.0) ARCRE=.5*SQRT(ANSRE) IF (ARCAE .LE. 0.0) ARCAE=.5*SQRT(ANSAE) NC=N NFEC=0 IFLAGC=IFLAG NP1=N+1 C SET INITIAL CONDITIONS FOR FIRST CALL TO STEPNF . START=.TRUE. CRASH=.FALSE. HOLD=1.0 H=.1 S=0.0 YPOLD(1)=1.0 YP(1)=1.0 Y(1)=0.0 DO 40 JW=2,NP1 YPOLD(JW)=0.0 YP(JW)=0.0 40 CONTINUE C SET OPTIMAL STEP SIZE ESTIMATION PARAMETERS. C LET Z[K] DENOTE THE NEWTON ITERATES ALONG THE FLOW NORMAL TO THE C DAVIDENKO FLOW AND Y THEIR LIMIT. C IDEAL CONTRACTION FACTOR: ||Z[2] - Z[1]|| / ||Z[1] - Z[0]|| IF (SSPAR(1) .LE. 0.0) SSPAR(1)= .5 C IDEAL RESIDUAL FACTOR: ||RHO(A, Z[1])|| / ||RHO(A, Z[0])|| IF (SSPAR(2) .LE. 0.0) SSPAR(2)= .01 C IDEAL DISTANCE FACTOR: ||Z[1] - Y|| / ||Z[0] - Y|| IF (SSPAR(3) .LE. 0.0) SSPAR(3)= .5 C MINIMUM STEP SIZE HMIN . IF (SSPAR(4) .LE. 0.0) SSPAR(4)= (SQRT(N+1.0)+4.0)*D1MACH(4) C MAXIMUM STEP SIZE HMAX . IF (SSPAR(5) .LE. 0.0) SSPAR(5)= 1.0 C MINIMUM STEP SIZE REDUCTION FACTOR BMIN . IF (SSPAR(6) .LE. 0.0) SSPAR(6)= .1 C MAXIMUM STEP SIZE EXPANSION FACTOR BMAX . IF (SSPAR(7) .LE. 0.0) SSPAR(7)= 3.0 C ASSUMED OPERATING ORDER P . IF (SSPAR(8) .LE. 0.0) SSPAR(8)= 2.0 C C LOAD A FOR THE FIXED POINT AND ZERO FINDING PROBLEMS. IF (IFLAGC .GE. -1) THEN DO 60 JW=2,NP1 A(JW-1)=Y(JW) 60 CONTINUE ENDIF 90 LIMIT=LIMITD C C ***** END OF INITIALIZATION BLOCK. ***** C C C ***** MAIN LOOP. ***** C 120 DO 400 ITER=1,LIMIT IF (Y(1) .LT. 0.0) THEN ARCLEN=S IFLAG=5 RETURN ENDIF C C SET DIFFERENT ERROR TOLERANCE IF THE TRAJECTORY Y(S) HAS ANY HIGH C CURVATURE COMPONENTS. 140 CURTOL=CURSW*HOLD RELERR=ARCRE ABSERR=ARCAE DO 160 JW=1,NP1 IF (ABS(YP(JW)-YPOLD(JW)) .GT. CURTOL) THEN RELERR=ANSRE ABSERR=ANSAE GO TO 200 ENDIF 160 CONTINUE C C TAKE A STEP ALONG THE CURVE. 200 CALL STEPNF(NC,NFEC,IFLAGC,START,CRASH,HOLD,H,RELERR,ABSERR, + S,Y,YP,YOLD,YPOLD,A,QR,ALPHA,TZ,PIVOT,W,WP,Z0,Z1,SSPAR, + PAR,IPAR) C PRINT LATEST POINT ON CURVE IF REQUESTED. IF (TRACE .GT. 0) THEN WRITE (TRACE,217) ITER,NFEC,S,Y(1),(Y(JW),JW=2,NP1) 217 FORMAT(/' STEP',I5,3X,'NFE =',I5,3X,'ARC LENGTH =',F9.4,3X, $ 'LAMBDA =',F7.4,5X,'X vector:'/1P,(1X,6E12.4)) ENDIF NFE=NFEC C CHECK IF THE STEP WAS SUCCESSFUL. IF (IFLAGC .GT. 0) THEN ARCLEN=S IFLAG=IFLAGC RETURN ENDIF IF (CRASH) THEN C RETURN CODE FOR ERROR TOLERANCE TOO SMALL. IFLAG=2 C CHANGE ERROR TOLERANCES. IF (ARCRE .LT. RELERR) ARCRE=RELERR IF (ANSRE .LT. RELERR) ANSRE=RELERR IF (ARCAE .LT. ABSERR) ARCAE=ABSERR IF (ANSAE .LT. ABSERR) ANSAE=ABSERR C CHANGE LIMIT ON NUMBER OF ITERATIONS. LIMIT=LIMIT-ITER RETURN ENDIF C IF (Y(1) .GE. 1.0) THEN C C USE HERMITE CUBIC INTERPOLATION AND NEWTON ITERATION TO GET THE C ANSWER AT LAMBDA = 1.0 . C C SAVE YOLD FOR ARC LENGTH CALCULATION LATER. DO 260 JW=1,NP1 Z0(JW)=YOLD(JW) 260 CONTINUE CALL ROOTNF(NC,NFEC,IFLAGC,ANSRE,ANSAE,Y,YP,YOLD,YPOLD, $ A,QR,ALPHA,TZ,PIVOT,W,WP,PAR,IPAR) C NFE=NFEC IFLAG=1 C SET ERROR FLAG IF ROOTNF COULD NOT GET THE POINT ON THE ZERO C CURVE AT LAMBDA = 1.0 . IF (IFLAGC .GT. 0) IFLAG=IFLAGC C CALCULATE FINAL ARC LENGTH. DO 290 JW=1,NP1 W(JW)=Y(JW) - Z0(JW) 290 CONTINUE ARCLEN=S - HOLD + DNRM2(NP1,W,1) RETURN ENDIF C C FOR POLYNOMIAL SYSTEMS AND THE POLSYS HOMOTOPY MAP, C D LAMBDA/DS .GE. 0 NECESSARILY. THIS CONDITION IS FORCED HERE IF C THE ENTRY POINT WAS POLYNF . C IF (POLSYS) THEN IF (YP(1) .LT. 0.0) THEN C REVERSE TANGENT DIRECTION SO D LAMBDA/DS = YP(1) > 0 . DO 310 JW=1,NP1 YP(JW)=-YP(JW) YPOLD(JW)=YP(JW) 310 CONTINUE C FORCE STEPNF TO USE THE LINEAR PREDICTOR FOR THE NEXT STEP ONLY. START=.TRUE. ENDIF ENDIF C 400 CONTINUE C C ***** END OF MAIN LOOP. ***** C C LAMBDA HAS NOT REACHED 1 IN 1000 STEPS. IFLAG=3 ARCLEN=S RETURN C END SUBROUTINE FIXPNS(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,A, $ NFE,ARCLEN,YP,YOLD,YPOLD,QR,LENQR,PIVOT,WORK,SSPAR, $ PAR,IPAR) C C SUBROUTINE FIXPNS FINDS A FIXED POINT OR ZERO OF THE C N-DIMENSIONAL VECTOR FUNCTION F(X), OR TRACKS A ZERO CURVE C OF A GENERAL HOMOTOPY MAP RHO(A,X,LAMBDA). FOR THE FIXED C POINT PROBLEM F(X) IS ASSUMED TO BE A C2 MAP OF SOME BALL C INTO ITSELF. THE EQUATION X = F(X) IS SOLVED BY C FOLLOWING THE ZERO CURVE OF THE HOMOTOPY MAP C C LAMBDA*(X - F(X)) + (1 - LAMBDA)*(X - A) , C C STARTING FROM LAMBDA = 0, X = A. THE CURVE IS PARAMETERIZED C BY ARC LENGTH S, AND IS FOLLOWED BY SOLVING THE ORDINARY C DIFFERENTIAL EQUATION D(HOMOTOPY MAP)/DS = 0 FOR C Y(S) = (X(S), LAMBDA(S)) USING A HERMITE CUBIC PREDICTOR AND A C CORRECTOR WHICH RETURNS TO THE ZERO CURVE ALONG THE FLOW NORMAL C TO THE DAVIDENKO FLOW (WHICH CONSISTS OF THE INTEGRAL CURVES OF C D(HOMOTOPY MAP)/DS ). C C FOR THE ZERO FINDING PROBLEM F(X) IS ASSUMED TO BE A C2 MAP C SUCH THAT FOR SOME R > 0, X*F(X) >= 0 WHENEVER NORM(X) = R. C THE EQUATION F(X) = 0 IS SOLVED BY FOLLOWING THE ZERO CURVE C OF THE HOMOTOPY MAP C C LAMBDA*F(X) + (1 - LAMBDA)*(X - A) C C EMANATING FROM LAMBDA = 0, X = A. C C A MUST BE AN INTERIOR POINT OF THE ABOVE MENTIONED BALLS. C C FOR THE CURVE TRACKING PROBLEM RHO(A,X,LAMBDA) IS ASSUMED TO C BE A C2 MAP FROM E**M X E**N X [0,1) INTO E**N, WHICH FOR C ALMOST ALL PARAMETER VECTORS A IN SOME NONEMPTY OPEN SUBSET C OF E**M SATISFIES C C RANK [D RHO(A,X,LAMBDA)/DX , D RHO(A,X,LAMBDA)/D LAMBDA] = N C C FOR ALL POINTS (X,LAMBDA) SUCH THAT RHO(A,X,LAMBDA)=0. IT IS C FURTHER ASSUMED THAT C C RANK [ D RHO(A,X0,0)/DX ] = N . C C WITH A FIXED, THE ZERO CURVE OF RHO(A,X,LAMBDA) EMANATING C FROM LAMBDA = 0, X = X0 IS TRACKED UNTIL LAMBDA = 1 BY C SOLVING THE ORDINARY DIFFERENTIAL EQUATION C D RHO(A,X(S),LAMBDA(S))/DS = 0 FOR Y(S) = (X(S), LAMBDA(S)), C WHERE S IS ARC LENGTH ALONG THE ZERO CURVE. ALSO THE HOMOTOPY C MAP RHO(A,X,LAMBDA) IS ASSUMED TO BE CONSTRUCTED SUCH THAT C C D LAMBDA(0)/DS > 0 . C C C FOR THE FIXED POINT AND ZERO FINDING PROBLEMS, THE USER C MUST SUPPLY A SUBROUTINE F(X,V) WHICH EVALUATES F(X) AT X C AND RETURNS THE VECTOR F(X) IN V, AND A SUBROUTINE C FJACS(X,QR,LENQR,PIVOT) WHICH EVALUATES THE (SYMMETRIC) C JACOBIAN MATRIX OF F(X) AT X, AND RETURNS THE SYMMETRIC C JACOBIAN MATRIX IN PACKED SKYLINE STORAGE FORMAT IN QR. LENQR C AND PIVOT DESCRIBE THE DATA STRUCTURE IN QR. FOR THE CURVE C TRACKING PROBLEM, THE USER MUST SUPPLY A SUBROUTINE C RHO(A,LAMBDA,X,V,PAR,IPAR) WHICH EVALUATES THE HOMOTOPY MAP RHO C AT (A,X,LAMBDA) AND RETURNS THE VECTOR RHO(A,X,LAMBDA) IN V, AND C A SUBROUTINE RHOJS(A,LAMBDA,X,QR,LENQR,PIVOT,PP,PAR,IPAR) WHICH C RETURNS IN QR THE SYMMETRIC N X N JACOBIAN MATRIX [D RHO/DX] C EVALUATED AT (A,X,LAMBDA) AND STORED IN PACKED SKYLINE FORMAT, AND C RETURNS IN PP THE VECTOR -(D RHO/D LAMBDA) EVALUATED AT C (A,X,LAMBDA). LENQR AND PIVOT DESCRIBE THE DATA STRUCTURE C IN QR. C *** NOTE THE MINUS SIGN IN THE DEFINITION OF PP. *** C C C FUNCTIONS AND SUBROUTINES DIRECTLY OR INDIRECTLY CALLED BY FIXPDS: C D1MACH , F (OR RHO ), FJACS (OR RHOJS ), GMFADS , MFACDS , C MULTDS , PCGDS , PCGNS , QIMUDS , ROOT , ROOTNS , SOLVDS , C STEPNS , TANGNS , AND THE BLAS FUNCTIONS DAXPY , DCOPY , DDOT , C DNRM2 , DSCAL , IDAMAX . ONLY D1MACH CONTAINS MACHINE DEPENDENT C CONSTANTS. NO OTHER MODIFICATIONS BY THE USER ARE REQUIRED. C C C ON INPUT: C C N IS THE DIMENSION OF X, F(X), AND RHO(A,X,LAMBDA). C C Y IS AN ARRRAY OF LENGTH N + 1. (Y(1),...,Y(N)) = A IS THE C STARTING POINT FOR THE ZERO CURVE FOR THE FIXED POINT AND C ZERO FINDING PROBLEMS. (Y(1),...,Y(N)) = X0 FOR THE CURVE C TRACKING PROBLEM. C C IFLAG CAN BE -2, -1, 0, 2, OR 3. IFLAG SHOULD BE 0 ON THE C FIRST CALL TO FIXPNS FOR THE PROBLEM X=F(X), -1 FOR THE C PROBLEM F(X)=0, AND -2 FOR THE PROBLEM RHO(A,X,LAMBDA)=0. C IN CERTAIN SITUATIONS IFLAG IS SET TO 2 OR 3 BY FIXPNS, C AND FIXPNS CAN BE CALLED AGAIN WITHOUT CHANGING IFLAG. C C ARCRE , ARCAE ARE THE RELATIVE AND ABSOLUTE ERRORS, RESPECTIVELY, C ALLOWED THE NORMAL FLOW ITERATION ALONG THE ZERO CURVE. IF C ARC?E .LE. 0.0 ON INPUT IT IS RESET TO .5*SQRT(ANS?E) . C NORMALLY ARC?E SHOULD BE CONSIDERABLY LARGER THAN ANS?E . C C ANSRE , ANSAE ARE THE RELATIVE AND ABSOLUTE ERROR VALUES USED FOR C THE ANSWER AT LAMBDA = 1. THE ACCEPTED ANSWER Y = (X, LAMBDA) C SATISFIES C C |Y(NP1) - 1| .LE. ANSRE + ANSAE .AND. C C ||Z|| .LE. ANSRE*||X|| + ANSAE WHERE C C (Z,.) IS THE NEWTON STEP TO Y. C C TRACE IS AN INTEGER SPECIFYING THE LOGICAL I/O UNIT FOR C INTERMEDIATE OUTPUT. IF TRACE .GT. 0 THE POINTS COMPUTED ON C THE ZERO CURVE ARE WRITTEN TO I/O UNIT TRACE . C C A(1:*) CONTAINS THE PARAMETER VECTOR A . FOR THE FIXED POINT C AND ZERO FINDING PROBLEMS, A NEED NOT BE INITIALIZED BY THE C USER, AND IS ASSUMED TO HAVE LENGTH N. FOR THE CURVE C TRACKING PROBLEM, A MUST BE INITIALIZED BY THE USER. C C YP(1:N+1) IS A WORK ARRAY CONTAINING THE TANGENT VECTOR TO C THE ZERO CURVE AT THE CURRENT POINT Y . C C YOLD(1:N+1) IS A WORK ARRAY CONTAINING THE PREVIOUS POINT FOUND C ON THE ZERO CURVE. C C YPOLD(1:N+1) IS A WORK ARRAY CONTAINING THE TANGENT VECTOR TO C THE ZERO CURVE AT YOLD . C C QR(1:LENQR) IS A WORK ARRAY CONTAINING THE N X N SYMMETRIC C JACOBIAN MATRIX WITH RESPECT TO X STORED IN PACKED SKYLINE C STORAGE FORMAT. LENQR AND PIVOT DESCRIBE THE DATA C STRUCTURE IN QR . C C LENQR IS THE LENGTH OF THE ONE-DIMENSIONAL ARRAY QR . C C PIVOT(1:N+2) IS A WORK ARRAY CONTAINING THE INDICES OF THE C DIAGONAL ELEMENTS OF THE N X N SYMMETRIC JACOBIAN MATRIX C (WITH RESPECT TO X) WITHIN QR . C C WORK(1:13*(N+1)+2*N+LENQR) IS A WORK ARRAY SPLIT UP AND USED C FOR THE CALCULATION OF THE JACOBIAN MATRIX KERNEL, THE C NEWTON STEP, INTERPOLATION, AND THE ESTIMATION OF THE OPTIMAL C STEP SIZE H . 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 FIXPNS . OTHERWISE THE INPUT VALUE OF SSPAR(J) IS USED. C SEE THE COMMENTS BELOW AND IN STEPNS FOR MORE INFORMATION ABOUT C THESE CONSTANTS. C C PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS, C WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES C RHO, RHOJS. C C C ON OUTPUT: C C N , TRACE , A ARE UNCHANGED. C C (Y(1),...,Y(N)) = X, Y(NP1) = LAMBDA, AND Y IS AN APPROXIMATE C ZERO OF THE HOMOTOPY MAP. NORMALLY LAMBDA = 1 AND X IS A C FIXED POINT(ZERO) OF F(X). IN ABNORMAL SITUATIONS LAMBDA C MAY ONLY BE NEAR 1 AND X IS NEAR A FIXED POINT(ZERO). C C IFLAG = C -2 CAUSES FIXPNS TO INITIALIZE EVERYTHING FOR THE PROBLEM C RHO(A,X,LAMBDA) = 0 (USE ON FIRST CALL). C C -1 CAUSES FIXPNS TO INITIALIZE EVERYTHING FOR THE PROBLEM C F(X) = 0 (USE ON FIRST CALL). C C 0 CAUSES FIXPNS TO INITIALIZE EVERYTHING FOR THE PROBLEM C X = F(X) (USE ON FIRST CALL). C C 1 NORMAL RETURN. C C 2 SPECIFIED ERROR TOLERANCE CANNOT BE MET. SOME OR ALL OF C ARCRE , ARCAE , ANSRE , ANSAE HAVE BEEN INCREASED TO C SUITABLE VALUES. TO CONTINUE, JUST CALL FIXPNS AGAIN C WITHOUT CHANGING ANY PARAMETERS. C C 3 STEPNS HAS BEEN CALLED 1000 TIMES. TO CONTINUE, CALL C FIXPNS AGAIN WITHOUT CHANGING ANY PARAMETERS. C C 4 THE PRECONDITIONED CONJUGATE GRADIENT ITERATION FAILED TO C CONVERGE (PROBABLY BECAUSE THE JACOBIAN MATRIX DID NOT HAVE C FULL RANK). THE ALGORITHM HAS FAILED (THE ZERO CURVE OF C THE HOMOTOPY MAP CANNOT BE FOLLOWED ANY FURTHER). C C 5 THE TRACKING ALGORITHM HAS LOST THE ZERO CURVE OF THE C HOMOTOPY MAP AND IS NOT MAKING PROGRESS. THE ERROR TOLERANCES C ARC?E AND ANS?E WERE TOO LENIENT. THE PROBLEM SHOULD BE C RESTARTED BY CALLING FIXPNS WITH SMALLER ERROR TOLERANCES C AND IFLAG = 0 (-1, -2). C C 6 THE NORMAL FLOW NEWTON ITERATION IN STEPNS OR ROOTNS C FAILED TO CONVERGE. THE ERROR TOLERANCES ANS?E MAY BE TOO C STRINGENT. C C 7 ILLEGAL INPUT PARAMETERS, A FATAL ERROR. C C ARCRE , ARCAE , ANSRE , ANSAE ARE UNCHANGED AFTER A NORMAL RETURN C (IFLAG = 1). THEY ARE INCREASED TO APPROPRIATE VALUES ON THE C RETURN IFLAG = 2 . C C NFE IS THE NUMBER OF FUNCTION EVALUATIONS (= NUMBER OF C JACOBIAN EVALUATIONS). C C ARCLEN IS THE LENGTH OF THE PATH FOLLOWED. C C C DOUBLE PRECISION ABSERR,ANSAE,ANSRE,ARCAE,ARCLEN,ARCRE, 1 CURSW,CURTOL,D1MACH,DNRM2,H,HOLD,RELERR,S INTEGER IFLAG,IFLAGC,IPP,IRHO,ITANGW,ITER,ITZ,IW,IWP, 1 IZ0,IZ1,JW,LENQR,LIMIT,LIMITD,N,NC,NFE,NFEC,NP1,TRACE LOGICAL START,CRASH C C ***** ARRAY DECLARATIONS. ***** C DOUBLE PRECISION Y(N+1),YP(N+1),YOLD(N+1),YPOLD(N+1),A(N), $ QR(LENQR),WORK(13*(N+1)+2*N+LENQR),SSPAR(8),PAR(1) INTEGER PIVOT(N+2),IPAR(1) C C ***** END OF DIMENSIONAL INFORMATION. ***** C SAVE C C LIMITD IS AN UPPER BOUND ON THE NUMBER OF STEPS. IT MAY BE C CHANGED BY CHANGING THE FOLLOWING PARAMETER STATEMENT: PARAMETER (LIMITD=1000) C C SWITCH FROM THE TOLERANCE ARC?E TO THE (FINER) TOLERANCE ANS?E IF C THE CURVATURE OF ANY COMPONENT OF Y EXCEEDS CURSW. PARAMETER (CURSW=10.0) C C C C : : : : : : : : : : : : : : : : : : : : : : : : IF (N .LE. 0 .OR. ANSRE .LE. 0.0 .OR. ANSAE .LT. 0.0) $ IFLAG=7 IF (IFLAG .GE. -2 .AND. IFLAG .LE. 0) GO TO 20 IF (IFLAG .EQ. 2) GO TO 120 IF (IFLAG .EQ. 3) GO TO 90 C ONLY VALID INPUT FOR IFLAG IS -2, -1, 0, 2, 3. IFLAG=7 RETURN C C ***** INITIALIZATION BLOCK. ***** C 20 ARCLEN=0.0 IF (ARCRE .LE. 0.0) ARCRE=.5*SQRT(ANSRE) IF (ARCAE .LE. 0.0) ARCAE=.5*SQRT(ANSAE) NC=N NFEC=0 IFLAGC=IFLAG NP1=N+1 C SET INDICES FOR SPLITTING UP WORK ARRAY. IPP=1 IRHO=N+1 IW=IRHO+N IWP=IW+NP1 ITZ=IWP+NP1 IZ0=ITZ+NP1 IZ1=IZ0+NP1 ITANGW=IZ1+NP1 C SET INITIAL CONDITIONS FOR FIRST CALL TO STEPNS . START=.TRUE. CRASH=.FALSE. HOLD=1.0 H=.1 S=0.0 YPOLD(NP1)=1.0 YP(NP1)=1.0 Y(NP1)=0.0 DO 40 JW=1,N YPOLD(JW)=0.0 YP(JW)=0.0 40 CONTINUE DO 50 JW=ITANGW,ITANGW+NP1+N WORK(JW)=0.0 50 CONTINUE C SET OPTIMAL STEP SIZE ESTIMATION PARAMETERS. C LET Z[K] DENOTE THE NEWTON ITERATES ALONG THE FLOW NORMAL TO THE C DAVIDENKO FLOW AND Y THEIR LIMIT. C IDEAL CONTRACTION FACTOR: ||Z[2] - Z[1]|| / ||Z[1] - Z[0]|| IF (SSPAR(1) .LE. 0.0) SSPAR(1)= .5 C IDEAL RESIDUAL FACTOR: ||RHO(A, Z[1])|| / ||RHO(A, Z[0])|| IF (SSPAR(2) .LE. 0.0) SSPAR(2)= .01 C IDEAL DISTANCE FACTOR: ||Z[1] - Y|| / ||Z[0] - Y|| IF (SSPAR(3) .LE. 0.0) SSPAR(3)= .5 C MINIMUM STEP SIZE HMIN . IF (SSPAR(4) .LE. 0.0) SSPAR(4)= (SQRT(N+1.0)+4.0)*D1MACH(4) C MAXIMUM STEP SIZE HMAX . IF (SSPAR(5) .LE. 0.0) SSPAR(5)= 1.0 C MINIMUM STEP SIZE REDUCTION FACTOR BMIN . IF (SSPAR(6) .LE. 0.0) SSPAR(6)= .1 C MAXIMUM STEP SIZE EXPANSION FACTOR BMAX . IF (SSPAR(7) .LE. 0.0) SSPAR(7)= 3.0 C ASSUMED OPERATING ORDER P . IF (SSPAR(8) .LE. 0.0) SSPAR(8)= 2.0 C C LOAD A FOR THE FIXED POINT AND ZERO FINDING PROBLEMS. IF (IFLAGC .GE. -1) THEN CALL DCOPY(N,Y,1,A,1) ENDIF 90 LIMIT=LIMITD C C ***** END OF INITIALIZATION BLOCK. ***** C C C ***** MAIN LOOP. ***** C 120 DO 400 ITER=1,LIMIT IF (Y(NP1) .LT. 0.0) THEN ARCLEN=S IFLAG=5 RETURN ENDIF C C SET DIFFERENT ERROR TOLERANCE IF THE TRAJECTORY Y(S) HAS ANY HIGH C CURVATURE COMPONENTS. 140 CURTOL=CURSW*HOLD RELERR=ARCRE ABSERR=ARCAE DO 160 JW=1,NP1 IF (ABS(YP(JW)-YPOLD(JW)) .GT. CURTOL) THEN RELERR=ANSRE ABSERR=ANSAE GO TO 200 ENDIF 160 CONTINUE C C TAKE A STEP ALONG THE CURVE. 200 CALL STEPNS(NC,NFEC,IFLAGC,START,CRASH,HOLD,H,RELERR,ABSERR, + S,Y,YP,YOLD,YPOLD,A,QR,LENQR,PIVOT,WORK,SSPAR,PAR,IPAR) C PRINT LATEST POINT ON CURVE IF REQUESTED. IF (TRACE .GT. 0) THEN WRITE (TRACE,217) ITER,NFEC,S,Y(NP1),(Y(JW),JW=1,NC) 217 FORMAT(/' STEP',I5,3X,'NFE =',I5,3X,'ARC LENGTH =',F9.4,3X, $ 'LAMBDA =',F7.4,5X,'X vector:'/1P,(1X,6E12.4)) ENDIF NFE=NFEC C CHECK IF THE STEP WAS SUCCESSFUL. IF (IFLAGC .GT. 0) THEN ARCLEN=S IFLAG=IFLAGC RETURN ENDIF IF (CRASH) THEN C RETURN CODE FOR ERROR TOLERANCE TOO SMALL. IFLAG=2 C CHANGE ERROR TOLERANCES. IF (ARCRE .LT. RELERR) ARCRE=RELERR IF (ANSRE .LT. RELERR) ANSRE=RELERR IF (ARCAE .LT. ABSERR) ARCAE=ABSERR IF (ANSAE .LT. ABSERR) ANSAE=ABSERR C CHANGE LIMIT ON NUMBER OF ITERATIONS. LIMIT=LIMIT-ITER RETURN ENDIF C IF (Y(NP1) .GE. 1.0) THEN C C USE HERMITE CUBIC INTERPOLATION AND NEWTON ITERATION TO GET THE C ANSWER AT LAMBDA = 1.0 . C C SAVE YOLD FOR ARC LENGTH CALCULATION LATER. CALL DCOPY(NP1,YOLD,1,WORK(IZ0),1) C CALL ROOTNS(NC,NFEC,IFLAGC,ANSRE,ANSAE,Y,YP,YOLD,YPOLD, $ A,QR,LENQR,PIVOT,WORK,PAR,IPAR) C NFE=NFEC IFLAG=1 C SET ERROR FLAG IF ROOTNS COULD NOT GET THE POINT ON THE ZERO C CURVE AT LAMBDA = 1.0 . IF (IFLAGC .GT. 0) IFLAG=IFLAGC C CALCULATE FINAL ARC LENGTH. DO 290 JW=1,NP1 WORK(JW)=Y(JW) - WORK(IZ0+JW-1) 290 CONTINUE ARCLEN=S - HOLD + DNRM2(NP1,WORK,1) RETURN ENDIF C 400 CONTINUE C C ***** END OF MAIN LOOP. ***** C C LAMBDA HAS NOT REACHED 1 IN 1000 STEPS. IFLAG=3 ARCLEN=S RETURN C C END SUBROUTINE FIXPQF(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,A, $ NFE,ARCLEN,YP,YOLD,YPOLD,QT,R,F0,F1,Z0,DZ,W,T,YSAV, $ SSPAR,PAR,IPAR) C C SUBROUTINE FIXPQF FINDS A FIXED POINT OR ZERO OF THE C N-DIMENSIONAL VECTOR FUNCTION F(X), OR TRACKS A ZERO CURVE OF A C GENERAL HOMOTOPY MAP RHO(A,LAMBDA,X). FOR THE FIXED POINT PROBLEM C F(X) IS ASSUMED TO BE A C2 MAP OF SOME BALL INTO ITSELF. THE C EQUATION X=F(X) IS SOLVED BY FOLLOWING THE ZERO CURVE OF THE C HOMOTOPY MAP C C LAMBDA*(X - F(X)) + (1 - LAMBDA)*(X - A), C C STARTING FROM LAMBDA = 0, X = A. THE CURVE IS PARAMETERIZED C BY ARC LENGTH S, AND IS FOLLOWED BY SOLVING THE ORDINARY C DIFFERENTIAL EQUATION D(HOMOTOPY MAP)/DS = 0 FOR C Y(S) = (LAMBDA(S), X(S)). THIS IS DONE BY USING A HERMITE CUBIC C PREDICTOR AND A CORRECTOR WHICH RETURNS TO THE ZERO CURVE IN A C HYPERPLANE PERPENDICULAR TO THE TANGENT TO THE ZERO CURVE AT THE C MOST RECENT POINT. C C FOR THE ZERO FINDING PROBLEM F(X) IS ASSUMED TO BE A C2 MAP C SUCH THAT FOR SOME R > 0, X*F(X) >= 0 WHENEVER NORM(X) = R. C THE EQUATION F(X) = 0 IS SOLVED BY FOLLOWING THE ZERO CURVE OF C THE HOMOTOPY MAP C C LAMBDA*F(X) + (1 - LAMBDA)*(X - A) C C EMANATING FROM LAMBDA = 0, X = A. C C A MUST BE AN INTERIOR POINT OF THE ABOVE MENTIONED BALLS. C C FOR THE CURVE TRACKING PROBLEM RHO(A,LAMBDA,X) IS ASSUMED TO C BE A C2 MAP FROM E**M X [0,1) X E**N INTO E**N, WHICH FOR C ALMOST ALL PARAMETER VECTORS A IN SOME NONEMPTY OPEN SUBSET C OF E**M SATISFIES C C RANK [D RHO(A,LAMBDA,X)/D LAMBDA, D RHO(A,LAMBDA,X)/DX] = N C C FOR ALL POINTS (LAMBDA,X) SUCH THAT RHO(A,LAMBDA,X) = 0. IT IS C FURTHER ASSUMED THAT C C RANK [ D RHO(A,0,X0)/DX ] = N. C C WITH A FIXED, THE ZERO CURVE OF RHO(A,LAMBDA,X) EMANATING FROM C LAMBDA = 0, X = X0 IS TRACKED UNTIL LAMBDA = 1 BY SOLVING THE C ORDINARY DIFFERENTIAL EQUATION D RHO(A,LAMBDA(S),X(S))/DS = 0 C FOR Y(S) = (LAMBDA(S), X(S)), WHERE S IS ARC LENGTH ALONG THE C ZERO CURVE. ALSO THE HOMOTOPY MAP RHO(A,LAMBDA,X) IS ASSUMED TO C BE CONSTRUCTED SUCH THAT C C D LAMBDA(0)/DS > 0. C C FOR THE FIXED POINT AND ZERO FINDING PROBLEMS, THE USER MUST SUPPLY C A SUBROUTINE F(X,V) WHICH EVALUATES F(X) AT X AND RETURNS THE C VECTOR F(X) IN V, AND A SUBROUTINE FJAC(X,V,K) WHICH RETURNS IN V C THE KTH COLUMN OF THE JACOBIAN MATRIX OF F(X) EVALUATED AT X. FOR C THE CURVE TRACKING PROBLEM, THE USER MUST SUPPLY A SUBROUTINE C RHO(A,LAMBDA,X,V,PAR,IPAR) WHICH EVALUATES THE HOMOTOPY MAP RHO AT C (A,LAMBDA,X) AND RETURNS THE VECTOR RHO(A,LAMBDA,X) IN V, AND C A SUBROUTINE RHOJAC(A,LAMBDA,X,V,K,PAR,IPAR) WHICH RETURNS IN V C THE KTH COLUMN OF THE N X (N+1) JACOBIAN MATRIX C [D RHO/D LAMBDA, D RHO/DX] EVALUATED AT (A,LAMBDA,X). FIXPQF C DIRECTLY OR INDIRECTLY USES THE SUBROUTINES D1MACH, F (OR RHO), C FJAC (OR RHOJAC), QRFAQF, QRSLQF, ROOT, ROOTQF, STEPQF, TANGQF, C UPQRQF AND THE BLAS ROUTINES DAXPY, DCOPY, DDOT, DNRM2, AND DSCAL. C ONLY D1MACH CONTAINS MACHINE DEPENDENT CONSTANTS. NO OTHER C MODIFICATIONS BY THE USER ARE REQUIRED. C C C ON INPUT: C C N IS THE DIMENSION OF X, F(X), AND RHO(A,LAMBDA,X). C C Y(1:N+1) CONTAINS THE STARTING POINT FOR TRACKING THE HOMOTOPY MAP. C (Y(2),...,Y(N+1)) = A FOR THE FIXED POINT AND ZERO FINDING C PROBLEMS. (Y(2),...,Y(N+1)) = X0 FOR THE CURVE TRACKING PROBLEM. C Y(1) NEED NOT BE DEFINED BY THE USER. C C IFLAG CAN BE -2, -1, 0, 2, OR 3. IFLAG SHOULD BE 0 ON THE FIRST C CALL TO FIXPQF FOR THE PROBLEM X=F(X), -1 FOR THE PROBLEM C F(X)=0, AND -2 FOR THE PROBLEM RHO(A,LAMBDA,X)=0. IN CERTAIN C SITUATIONS IFLAG IS SET TO 2 OR 3 BY FIXPQF, AND FIXPQF CAN C BE CALLED AGAIN WITHOUT CHANGING IFLAG. C C ARCRE, ARCAE ARE THE RELATIVE AND ABSOLUTE ERRORS, RESPECTIVELY, C ALLOWED THE QUASI-NEWTON ITERATION ALONG THE ZERO CURVE. IF C ARC?E .LE. 0.0 ON INPUT, IT IS RESET TO .5*SQRT(ANS?E). C NORMALLY ARC?E SHOULD BE CONSIDERABLY LARGER THAN ANS?E. C C ANSRE, ANSAE ARE THE RELATIVE AND ABSOLUTE ERROR VALUES USED FOR C THE ANSWER AT LAMBDA = 1. THE ACCEPTED ANSWER Y = (LAMBDA, X) C SATISFIES C C |Y(1) - 1| .LE. ANSRE + ANSAE .AND. C C ||DZ|| .LE. ANSRE*||Y|| + ANSAE WHERE C C DZ IS THE QUASI-NEWTON STEP TO Y. C C TRACE IS AN INTEGER SPECIFYING THE LOGICAL I/O UNIT FOR C INTERMEDIATE OUTPUT. IF TRACE .GT. 0 THE POINTS COMPUTED ON C THE ZERO CURVE ARE WRITTEN TO I/O UNIT TRACE . C C A(1:*) CONTAINS THE PARAMETER VECTOR A. FOR THE FIXED POINT C AND ZERO FINDING PROBLEMS, A NEED NOT BE INITIALIZED BY THE C USER, AND IS ASSUMED TO HAVE LENGTH N. FOR THE CURVE C TRACKING PROBLEM, A MUST BE INITIALIZED BY THE USER. C C YP(1:N+1) IS A WORK ARRAY CONTAINING THE TANGENT VECTOR TO THE C ZERO CURVE AT THE CURRENT POINT Y. C C YOLD(1:N+1) IS A WORK ARRAY CONTAINING THE PREVIOUS POINT FOUND C ON THE ZERO CURVE. C C YPOLD(1:N+1) IS A WORK ARRAY CONTAINING THE TANGENT VECTOR TO C THE ZERO CURVE AT YOLD. C C QT(1:N+1,1:N+1), R((N+1)*(N+2)/2), F0(1:N+1), F1(1:N+1), Z0(1:N+1), C DZ(1:N+1), W(1:N+1), T(1:N+1), YSAV(1:N+1) ARE ALL WORK ARRAYS C USED BY STEPQF, TANGQF AND ROOTQF TO CALCULATE THE TANGENT C VECTORS AND QUASI-NEWTON STEPS. C C SSPAR(1:4) = (HMIN, HMAX, BMIN, BMAX) IS A VECTOR OF PARAMETERS C USED FOR THE OPTIMAL STEP SIZE ESTIMATION. A DEFAULT VALUE C CAN BE SPECIFIED FOR ANY OF THESE FOUR PARAMETERS BY SETTING IT C .LE. 0.0 ON INPUT. SEE THE COMMENTS IN STEPQF FOR MORE C INFORMATION ABOUT THESE PARAMETERS. C C PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS, C WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES C RHO, RHOJAC. C C C ON OUTPUT: C C N , TRACE , A ARE UNCHANGED. C C Y(1) = LAMBDA, (Y(2),...,Y(N+1)) = X, AND Y IS AN APPROXIMATE C ZERO OF THE HOMOTOPY MAP. NORMALLY LAMBDA = 1 AND X IS A C FIXED POINT OR ZERO OF F(X). IN ABNORMAL SITUATIONS, LAMBDA C MAY ONLY BE NEAR 1 AND X NEAR A FIXED POINT OR ZERO. C C IFLAG = C C 1 NORMAL RETURN C C 2 SPECIFIED ERROR TOLERANCE CANNOT BE MET. SOME OR ALL OF C ARCRE, ARCAE, ANSRE, ANSAE HAVE BEEN INCREASED TO C SUITABLE VALUES. TO CONTINUE, JUST CALL FIXPQF AGAIN C WITHOUT CHANGING ANY PARAMETERS. C C 3 STEPQF HAS BEEN CALLED 1000 TIMES. TO CONTINUE, CALL C FIXPQF AGAIN WITHOUT CHANGING ANY PARAMETERS. C C 4 JACOBIAN MATRIX DOES NOT HAVE FULL RANK. THE ALGORITHM C HAS FAILED (THE ZERO CURVE OF THE HOMOTOPY MAP CANNOT BE C FOLLOWED ANY FURTHER). C C 5 THE TRACKING ALGORITHM HAS LOST THE ZERO CURVE OF THE C HOMOTOPY MAP AND IS NOT MAKING PROGRESS. THE ERROR C TOLERANCES ARC?E AND ANS?E WERE TOO LENIENT. THE PROBLEM C SHOULD BE RESTRARTED BY CALLING FIXPQF WITH SMALLER ERROR C TOLERANCES AND IFLAG = 0 (-1, -2). C C 6 THE QUASI-NEWTON ITERATION IN STEPQF OR ROOTQF FAILED TO C CONVERGE. THE ERROR TOLERANCES ANS?E MAY BE TOO STRINGENT. C C 7 ILLEGAL INPUT PARAMETERS, A FATAL ERROR. C C ARCRE, ARCAE, ANSRE, ANSAE ARE UNCHANGED AFTER A NORMAL RETURN C (IFLAG = 1). THEY ARE INCREASED TO APPROPRIATE VALUES ON THE C RETURN IFLAG = 2. C C NFE IS THE NUMBER OF JACOBIAN EVALUATIONS. C C ARCLEN IS THE APPROXIMATE LENGTH OF THE ZERO CURVE. C C ***** DECLARATIONS ***** C C FUNCTION DECLARATIONS C DOUBLE PRECISION D1MACH, DNRM2 C C LOCAL VARIABLES C DOUBLE PRECISION ABSERR, H, HOLD, RELERR, S, WK INTEGER IFLAGC, ITER, JW, LIMITD, LIMIT, NP1 LOGICAL CRASH, START C C SCALAR ARGUMENTS C DOUBLE PRECISION ARCRE, ARCAE, ANSRE, ANSAE, ARCLEN INTEGER N,IFLAG,TRACE,NFE C C ARRAY DECLARATIONS C DOUBLE PRECISION A(N), Y(N+1), YP(N+1), YOLD(N+1), YPOLD(N+1), $ QT(N+1,N+1), R((N+1)*(N+2)/2), F0(N+1), F1(N+1), Z0(N+1), $ DZ(N+1), W(N+1), T(N+1), YSAV(N+1), SSPAR(4), PAR(1) INTEGER IPAR(1) C SAVE C C ***** END OF DECLARATIONS ***** C C LIMITD IS AN UPPER BOUND ON THE NUMBER OF STEPS. IT MAY BE C CHANGED BY CHANGING THE FOLLOWING PARAMETER STATEMENT: PARAMETER (LIMITD =1000) C C ***** FIRST EXECUTABLE STATEMENT ***** C C CHECK IFLAG C IF (N .LE. 0 .OR. ANSRE .LE. 0.0 .OR. ANSAE .LT. 0.0) $ IFLAG = 7 IF (IFLAG .GE. -2 .AND. IFLAG .LE. 0) GO TO 10 IF (IFLAG .EQ. 2) GO TO 50 IF (IFLAG .EQ. 3) GO TO 40 C C ONLY VALID INPUT FOR IFLAG IS -2, -1, 0, 2, 3. C IFLAG = 7 RETURN C C ***** INITIALIZATION BLOCK ***** C 10 ARCLEN = 0.0 IF (ARCRE .LE. 0.0) ARCRE = .5*SQRT(ANSRE) IF (ARCAE .LE. 0.0) ARCAE = .5*SQRT(ANSAE) NFE=0 IFLAGC = IFLAG NP1=N+1 C C SET INITIAL CONDITIONS FOR FIRST CALL TO STEPQF. C START=.TRUE. CRASH=.FALSE. RELERR = ARCRE ABSERR = ARCAE HOLD=1.0 H=0.1 S=0.0 YPOLD(1) = 1.0 Y(1) = 0.0 DO 20 JW=2,NP1 YPOLD(JW)=0.0 20 CONTINUE C C SET OPTIMAL STEP SIZE ESTIMATION PARAMETERS. C C MINIMUM STEP SIZE HMIN IF (SSPAR(1) .LE. 0.0) SSPAR(1)= (SQRT(N+1.0)+4.0)*D1MACH(4) C MAXIMUM STEP SIZE HMAX IF (SSPAR(2) .LE. 0.0) SSPAR(2)= 1.0 C MINIMUM STEP REDUCTION FACTOR BMIN IF (SSPAR(3) .LE. 0.0) SSPAR(3)= 0.1 C MAXIMUM STEP EXPANSION FACTOR BMAX IF (SSPAR(4) .LE. 0.0) SSPAR(4)= 7.0 C C LOAD A FOR THE FIXED POINT AND ZERO FINDING PROBLEMS. C IF (IFLAGC .GE. -1) THEN CALL DCOPY(N,Y(2),1,A,1) ENDIF C 40 LIMIT=LIMITD C C ***** END OF INITIALIZATION BLOCK. ***** C C ***** MAIN LOOP. ***** C 50 DO 400 ITER=1,LIMIT IF (Y(1) .LT. 0.0) THEN ARCLEN = S IFLAG = 5 RETURN END IF C C TAKE A STEP ALONG THE CURVE. C CALL STEPQF(N,NFE,IFLAGC,START,CRASH,HOLD,H,WK, $ RELERR,ABSERR,S,Y,YP,YOLD,YPOLD,A,QT,R,F0,F1,Z0,DZ, $ W,T,SSPAR,PAR,IPAR) C C PRINT LATEST POINT ON CURVE IF REQUESTED. C IF (TRACE .GT. 0) THEN WRITE (TRACE,217) ITER,NFE,S,Y(1),(Y(JW),JW=2,NP1) 217 FORMAT(/' STEP',I5,3X,'NFE =',I5,3X,'ARC LENGTH =',F9.4,3X, $ 'LAMBDA =',F7.4,5X,'X vector:'/1P,(1X,6E12.4)) ENDIF C C CHECK IF THE STEP WAS SUCCESSFUL. C IF (IFLAGC .GT. 0) THEN ARCLEN=S IFLAG=IFLAGC RETURN END IF C IF (CRASH) THEN C C RETURN CODE FOR ERROR TOLERANCE TOO SMALL. C IFLAG=2 C C CHANGE ERROR TOLERANCES. C IF (ARCRE .LT. RELERR) THEN ARCRE=RELERR ANSRE=RELERR END IF IF (ARCAE .LT. ABSERR) ARCAE=ABSERR C C CHANGE LIMIT ON NUMBER OF ITERATIONS. C LIMIT = LIMIT - ITER RETURN END IF C C IF LAMBDA >= 1.0, USE ROOTQF TO FIND SOLUTION. C IF (Y(1) .GE. 1.0) GOTO 500 C 400 CONTINUE C C ***** END OF MAIN LOOP ***** C C DID NOT CONVERGE IN LIMIT ITERATIONS, SET IFLAG AND RETURN. C ARCLEN = S IFLAG = 3 RETURN C C ***** FINAL STEP -- FIND SOLUTION AT LAMBDA=1 ***** C C SAVE YOLD FOR ARC LENGTH CALCULATION LATER. C 500 CALL DCOPY(NP1,YOLD,1,YSAV,1) C C FIND SOLUTION. C CALL ROOTQF(N,NFE,IFLAGC,ANSRE,ANSAE,Y,YP,YOLD, $ YPOLD,A,QT,R,DZ,Z0,W,T,F0,F1,PAR,IPAR) C C CHECK IF SOLUTION WAS FOUND AND SET IFLAG ACCORDINGLY. C IFLAG=1 C C SET ERROR FLAG IF ROOTQF COULD NOT GET THE POINT ON THE ZERO C CURVE AT LAMBDA = 1.0. C