C ALGORITHM 596, COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.9, NO. 2, C JUN., 1983, P. 236-241; also P. 215-235. C***********************************************************************MAN 10 C MAN 20 C TEST MAIN PROGRAM FOR CONTINUATION CODE MAN 30 C AIRCRAFT STABILITY PROBLEM MAN 40 C MAN 50 C RUN WITH X(6) = VAL1 = -.05, -.008, 0.0, 0.1 MAN 60 C MAN 70 C SEEKING LIMIT POINTS IN X(7). MAN 80 C MAN 90 C X(1) = ROLL MAN 100 C X(2) = PITCH MAN 110 C X(3) = YAW MAN 120 C X(4) = ANGLE OF ATTACK MAN 130 C X(5) = SIDESLIP MAN 140 C X(6) = ELEVATOR MAN 150 C X(7) = AILERON MAN 160 C X(8) = RUDDER MAN 170 C MAN 180 C***********************************************************************MAN 190 C MAN 200 EXTERNAL FXAIR, FPAIR, FSOLVE MAN 210 REAL RWORK(104) MAN 220 INTEGER IPIVOT(8) MAN 230 COMMON /COEF/ BARRAY(5,8) MAN 240 COMMON /CONTRL/ IFIX1, IFIX2, VAL1, VAL2 MAN 250 COMMON /COUNT1/ ICRSL, ITNSL, NSTCR, NCNCR, NTRCR, NLMCR, NLMRT MAN 260 COMMON /COUNT2/ IFEVAL, IPEVAL, ISOLVE, NREDXF, NRDSUM, NCORXR, MAN 270 * NCRSUM MAN 280 COMMON /DETEXP/ IEXTXL, IEXTXC, IEXTXR, IEXNXR MAN 290 COMMON /DETMAN/ DETTXL, DETTXC, DETTXR, DETNXR MAN 300 COMMON /OUTPUT/ IWRITE MAN 310 COMMON /POINT/ CURVCF, CORDXF, ALPHLC, HSECLC, FNRM MAN 320 IWRITE = 1 MAN 330 LUNIT = 6 MAN 340 MAXSTP = 25 MAN 350 C MAN 360 C IFIX1 AND IFIX2 RECORD THE INDICES OF X WHICH ARE TO BE HELD FIXED. MAN 370 C VAL1 AND VAL2 ARE THE VALUES AT WHICH THESE QUANTITIES ARE HELD. MAN 380 C MAN 390 IFIX1 = 6 MAN 400 IFIX2 = 8 MAN 410 VAL1 = -.05 MAN 420 VAL2 = 0.0 MAN 430 C MAN 440 C SET PITCON ARGUMENTS MAN 450 C MAN 460 NVAR = 8 MAN 470 LIM = 7 MAN 480 IT = 0 MAN 490 XIT = 0.0 MAN 500 KSTEP2 = -3 MAN 510 KSTEP1 = -2 MAN 520 KSTEP = -1 MAN 530 IPC = 7 MAN 540 IPCFIX = 0 MAN 550 DIRIPC = -1.0 MAN 560 HTANCF = .3 MAN 570 IRET = 0 MAN 580 MODCON = 0 MAN 590 HMAX = 100.0 MAN 600 HMIN = .001 MAN 610 HFACT = 3.0 MAN 620 ABSERR = .00001 MAN 630 RELERR = .00001 MAN 640 RWORK(1) = 0.0 MAN 650 RWORK(2) = 0.0 MAN 660 RWORK(3) = 0.0 MAN 670 RWORK(4) = 0.0 MAN 680 RWORK(5) = 0.0 MAN 690 RWORK(6) = 0.0 MAN 700 RWORK(7) = 0.0 MAN 710 RWORK(8) = 0.0 MAN 720 RWORK(IFIX1) = VAL1 MAN 730 RWORK(IFIX2) = VAL2 MAN 740 ISIZE = 104 MAN 750 NROW = NVAR MAN 760 NCOL = NVAR MAN 770 C MAN 780 C SET LOCATION POINTERS MAN 790 C MAN 800 JXR = 0 MAN 810 IXR = 1 MAN 820 JXC = JXR + NVAR MAN 830 IXC = IXR + NVAR MAN 840 JXF = JXC + NVAR MAN 850 IXF = IXC + NVAR MAN 860 JTL = JXF + NVAR MAN 870 ITL = IXF + NVAR MAN 880 JTC = JTL + NVAR MAN 890 ITC = ITL + NVAR MAN 900 JFP = JTC + NVAR MAN 910 IFP = ITC + NVAR MAN 920 C MAN 930 C SET B MATRIX MAN 940 C MAN 950 CALL STORE(BARRAY) MAN 960 WRITE (LUNIT,99987) RWORK(6), RWORK(7), RWORK(8) MAN 970 WRITE (LUNIT,99988) ABSERR, RELERR, HTANCF, HFACT, HMIN, HMAX, MAN 980 * MODCON MAN 990 IF (LIM.NE.0) WRITE (LUNIT,99984) LIM MAN 1000 IF (IPCFIX.NE.0) WRITE (LUNIT,99983) IPC MAN 1010 C MAN 1020 C BEGIN TRACING OUT THE CURVE MAN 1030 C MAN 1040 DO 20 I=1,MAXSTP MAN 1050 CALL PITCON(NVAR, LIM, IT, XIT, KSTEP, IPC, IPCFIX, DIRIPC, MAN 1060 * HTANCF, IRET, MODCON, IPIVOT, HMAX, HMIN, HFACT, ABSERR, MAN 1070 * RELERR, RWORK, ISIZE, NROW, NCOL, FXAIR, FPAIR, FSOLVE, LUNIT) MAN 1080 IF (IRET.LT.(-1)) GO TO 30 MAN 1090 IF (IRET.LT.0) GO TO 20 MAN 1100 IF (IRET.EQ.1) GO TO 40 MAN 1110 IF (IRET.EQ.2) GO TO 10 MAN 1120 C MAN 1130 C NORMAL CONTINUATION POINT RETURN MAN 1140 C MAN 1150 KSTEP1 = KSTEP - 1 MAN 1160 KSTEP2 = KSTEP1 - 1 MAN 1170 IF (KSTEP.GT.1) WRITE (LUNIT,99996) KSTEP2, KSTEP1, HSECLC MAN 1180 IHI = JTC + NVAR MAN 1190 WRITE (LUNIT,99995) KSTEP1 MAN 1200 WRITE (LUNIT,99982) (RWORK(J),J=ITC,JFP) MAN 1210 WRITE (LUNIT,99994) KSTEP, HTANCF MAN 1220 WRITE (LUNIT,99993) KSTEP MAN 1230 WRITE (LUNIT,99982) (RWORK(J),J=IXR,JXC) MAN 1240 WRITE (LUNIT,99992) FNRM MAN 1250 GO TO 20 MAN 1260 C MAN 1270 C LIMIT POINT FOUND, PRINT OUT AND CONTINUE MAN 1280 C MAN 1290 10 WRITE (LUNIT,99991) MAN 1300 WRITE (LUNIT,99986) MAN 1310 WRITE (LUNIT,99982) (RWORK(J),J=IXR,JXC) MAN 1320 WRITE (LUNIT,99989) MAN 1330 WRITE (LUNIT,99982) (RWORK(J),J=ITL,JTC) MAN 1340 WRITE (LUNIT,99992) FNRM MAN 1350 20 CONTINUE MAN 1360 C MAN 1370 C LOOP COMPLETED WITHOUT REACHING TARGET POINT MAN 1380 C MAN 1390 WRITE (LUNIT,99999) MAN 1400 GO TO 50 MAN 1410 C MAN 1420 C ERROR RETURN FROM CONTINUATION CODE MAN 1430 C MAN 1440 30 WRITE (LUNIT,99998) IRET MAN 1450 GO TO 50 MAN 1460 C MAN 1470 C TARGET POINT REACHED MAN 1480 C MAN 1490 40 WRITE (LUNIT,99997) MAN 1500 WRITE (LUNIT,99986) MAN 1510 WRITE (LUNIT,99982) (RWORK(J),J=IXR,JXC) MAN 1520 WRITE (LUNIT,99992) FNRM MAN 1530 50 WRITE (LUNIT,99990) NROW, NCOL, KSTEP, NCRSUM, NRDSUM, IFEVAL, MAN 1540 * IPEVAL, ISOLVE MAN 1550 WRITE (LUNIT,99985) ICRSL, ITNSL, NSTCR, NCNCR, NTRCR, NLMCR, MAN 1560 * NLMRT MAN 1570 STOP MAN 1580 C MAN 1590 C***********************************************************************MAN 1600 C MAN 1610 99999 FORMAT (47H0PROGRAM TOOK FULL NUMBER OF CONTINUATION STEPS) MAN 1620 99998 FORMAT (24H0ERROR RETURN WITH IRET=, I4) MAN 1630 99997 FORMAT (19H0TARGET FOUND, STOP) MAN 1640 99996 FORMAT (17H0SECANT STEP FROM, I3, 4H TO , I3, 3H = , G14.7) MAN 1650 99995 FORMAT (18H TANGENT AT POINT , I3, 5H WAS ) MAN 1660 99994 FORMAT (16H TO REACH POINT , I3, 22H TANGENT STEPSIZE WAS , G14.7)MAN 1670 99993 FORMAT (7H POINT , I3, 4H IS ) MAN 1680 99992 FORMAT (15H FUNCTION NORM=, G14.7) MAN 1690 99991 FORMAT (18H0LIMIT POINT FOUND) MAN 1700 99990 FORMAT (17H0MATRIX SIZE , I4, 1HX, I4/17H STEPS TAKEN , MAN 1710 * I4/17H NEWTON STEPS , I4/17H STEP REDUCTIONS , I4/9H FUNCTION,MAN 1720 * 8H CALLS , I4/17H PRIME CALLS , I4/17H SOLVE CALLS , I4)MAN 1730 99989 FORMAT (9H TANGENT ) MAN 1740 99988 FORMAT (43H0UNIVERSITY OF PITTSBURGH CONTINUATION CODE/8H ABSERR=,MAN 1750 * G14.7, 8H RELERR=, G14.7/8H H= , G14.7, 8H HFACT= , MAN 1760 * G14.7/8H HMIN= , G14.7, 8H HMAX= , G14.7/18H NEWTON METHOD OPT,MAN 1770 * 4HION=, I2) MAN 1780 99987 FORMAT (27H1AIRCRAFT STABILITY PROBLEM/10H ELEVATOR=, MAN 1790 * G14.7/10H AILERON =, G14.7/10H RUDDER =, G14.7) MAN 1800 99986 FORMAT (9H POINT ) MAN 1810 99985 FORMAT (8H0CALLS -/40H SOLVE CALLED FOR CORRECTION , MAN 1820 * I3/40H SOLVE CALLED FOR TANGENT , I3/10H0CORRECTOR,MAN 1830 * 30H CALLED FOR STARTING POINT , I3/23H CORRECTOR CALLED FOR C,MAN 1840 * 17HONTINUATION POINT, I3/36H CORRECTOR CALLED FOR TARGET POINT ,MAN 1850 * 4H , I3/40H CORRECTOR CALLED FOR LIMIT POINT , MAN 1860 * I3/40H0ROOTFINDER CALLED FOR LIMIT POINT , I3) MAN 1870 99984 FORMAT (33H0LIMIT POINTS SOUGHT FOR VARIABLE, I3) MAN 1880 99983 FORMAT (36H CONTINUATION VARIABLE FORCED TO BE , I3) MAN 1890 99982 FORMAT (1X, 5G14.7) MAN 1900 END MAN 1910 SUBROUTINE STORE(BARRAY) STO 10 C C*********************************************************************** C REAL BARRAY(5,8) BARRAY(1,1) = -3.933 BARRAY(2,1) = 0.0 BARRAY(3,1) = 0.002 BARRAY(4,1) = 0.0 BARRAY(5,1) = 0.0 BARRAY(1,2) = 0.107 BARRAY(2,2) = -0.987 BARRAY(3,2) = 0.0 BARRAY(4,2) = 1.0 BARRAY(5,2) = 0.0 BARRAY(1,3) = 0.126 BARRAY(2,3) = 0.0 BARRAY(3,3) = -0.235 BARRAY(4,3) = 0.0 BARRAY(5,3) = -1.0 BARRAY(1,4) = 0.0 BARRAY(2,4) = -22.95 BARRAY(3,4) = 0.0 BARRAY(4,4) = -1.0 BARRAY(5,4) = 0.0 BARRAY(1,5) = -9.99 BARRAY(2,5) = 0.0 BARRAY(3,5) = 5.67 BARRAY(4,5) = 0.0 BARRAY(5,5) = -0.196 BARRAY(1,6) = 0.0 BARRAY(2,6) = -28.37 BARRAY(3,6) = 0.0 BARRAY(4,6) = -0.168 BARRAY(5,6) = 0.0 BARRAY(1,7) = -45.83 BARRAY(2,7) = 0.0 BARRAY(3,7) = -0.921 BARRAY(4,7) = 0.0 BARRAY(5,7) = -0.0071 BARRAY(1,8) = -7.64 BARRAY(2,8) = 0.0 BARRAY(3,8) = -6.51 BARRAY(4,8) = 0.0 BARRAY(5,8) = 0.0 RETURN C C*********************************************************************** C END SUBROUTINE FXAIR(NVAR, X, FX) FXA 10 C C*********************************************************************** C C THE FUNCTION IS OF THE FOLLOWING FORM C C FOR INDICES I=1 THROUGH 5, C C F(I)=SUM ON J AND K OF B(I,J)*X(J)+PHI(I,J,K)*X(J)*X(K) C C F(6)=X(IFIX1)-VAL1 C F(7)=X(IFIX2)-VAL2 C C*********************************************************************** C REAL X(NVAR), FX(NVAR) COMMON /COEF/ BARRAY(5,8) COMMON /CONTRL/ IFIX1, IFIX2, VAL1, VAL2 C C COMPUTE LINEAR TERMS C DO 20 I=1,5 FX(I) = 0.0 DO 10 J=1,8 FX(I) = FX(I) + BARRAY(I,J)*X(J) 10 CONTINUE 20 CONTINUE C C COMPUTE NONLINEAR TERMS C PHI = -.727*X(2)*X(3) + 8.39*X(3)*X(4) - 684.4*X(4)*X(5) + * 63.5*X(4)*X(7) FX(1) = FX(1) + PHI PHI = .949*X(1)*X(3) + .173*X(1)*X(5) FX(2) = FX(2) + PHI PHI = -.716*X(1)*X(2) - 1.578*X(1)*X(4) + 1.132*X(4)*X(7) FX(3) = FX(3) + PHI PHI = -X(1)*X(5) FX(4) = FX(4) + PHI PHI = X(1)*X(4) FX(5) = FX(5) + PHI C C SET FUNCTION VALUES FOR TWO FIXED VARIABLES C FX(6) = X(IFIX1) - VAL1 FX(7) = X(IFIX2) - VAL2 RETURN C C*********************************************************************** C END SUBROUTINE FPAIR(NVAR, X, FPRYM, NROW, NCOL) FPA 10 C C*********************************************************************** C REAL X(NVAR), FPRYM(NROW,NCOL) COMMON /COEF/ BARRAY(5,8) COMMON /CONTRL/ IFIX1, IFIX2, VAL1, VAL2 C C COMPUTE TERMS FROM LINEAR PART OF FUNCTION C DO 20 I=1,5 DO 10 J=1,8 FPRYM(I,J) = BARRAY(I,J) 10 CONTINUE 20 CONTINUE DO 40 I=6,7 DO 30 J=1,8 FPRYM(I,J) = 0.0 30 CONTINUE 40 CONTINUE FPRYM(6,IFIX1) = 1.0 FPRYM(7,IFIX2) = 1.0 C C COMPUTE TERMS FROM NONLINEAR PART OF FUNCTION C FPRYM(1,2) = FPRYM(1,2) - .727*X(3) FPRYM(1,3) = FPRYM(1,3) - .727*X(2) + 8.39*X(4) FPRYM(1,4) = FPRYM(1,4) + 8.39*X(3) - 684.4*X(5) + 63.5*X(7) FPRYM(1,5) = FPRYM(1,5) - 684.4*X(4) FPRYM(1,7) = FPRYM(1,7) + 63.5*X(4) FPRYM(2,1) = FPRYM(2,1) + .949*X(3) + .173*X(5) FPRYM(2,3) = FPRYM(2,3) + .949*X(1) FPRYM(2,5) = FPRYM(2,5) + .173*X(1) FPRYM(3,1) = FPRYM(3,1) - .716*X(2) - 1.578*X(4) FPRYM(3,2) = FPRYM(3,2) - .716*X(1) FPRYM(3,4) = FPRYM(3,4) - 1.578*X(1) + 1.132*X(7) FPRYM(3,7) = FPRYM(3,7) + 1.132*X(4) FPRYM(4,1) = FPRYM(4,1) - X(5) FPRYM(4,5) = FPRYM(4,5) - X(1) FPRYM(5,1) = FPRYM(5,1) + X(4) FPRYM(5,4) = FPRYM(5,4) + X(1) RETURN C C*********************************************************************** C END SUBROUTINE PITCON(NVAR, LIM, IT, XIT, KSTEP, IPC, IPCFIX, DIRIPC, PIT 10 * HTANCF, IRET, MODCON, IPIVOT, HMAX, HMIN, HFACT, ABSERR, RELERR, * RWORK, ISIZE, NROW, NCOL, FXNAME, FPNAME, SLNAME, LUNIT) C C*********************************************************************** C C 1. INTRODUCTION C C THIS IS THE 30 JUNE 1982 VERSION OF PITCON, C THE UNIVERSITY OF PITTSBURGH CONTINUATION PACKAGE. C THIS VERSION USES SINGLE PRECISION AND FULL MATRIX STORAGE. C C THIS PACKAGE WAS PREPARED WITH THE PARTIAL SUPPORT OF C THE NATIONAL SCIENCE FOUNDATION, UNDER GRANT MCS-78-05299, C BY WERNER C. RHEINBOLDT AND JOHN V. BURKARDT, C UNIVERSITY OF PITTSBURGH, PITTSBURGH, PA 15261. C C SUBROUTINE PITCON COMPUTES POINTS ALONG A SOLUTION CURVE OF AN C UNDERDETERMINED SYSTEM OF NONLINEAR EQUATIONS OF THE FORM FX=0. C THE CURVE IS SPECIFIED TO BEGIN AT A GIVEN STARTING SOLUTION C X OF THE SYSTEM. HERE X DENOTES A REAL VECTOR OF NVAR C COMPONENTS AND FX A REAL VECTOR OF NVAR-1 COMPONENTS. C NORMALLY EACH CALL TO PITCON PRODUCES A NEW POINT FURTHER ALONG C THE SOLUTION CURVE IN A USER-SPECIFIED DIRECTION. C C AN OPTION ALLOWS THE SEARCH FOR AND COMPUTATION OF TARGET POINTS, C THAT IS, SOLUTION POINTS X FOR WHICH X(IT) = XIT FOR SOME USER C SPECIFIED VALUES OF IT AND XIT. C C A FURTHER OPTION ALLOWS THE SEARCH FOR AND COMPUTATION OF LIMIT C POINTS FOR SPECIFIED COORDINATE LIM, THAT IS, SOLUTION POINTS FOR C WHICH THE LIM-TH COMPONENT OF THE TANGENT VECTOR IS ZERO. C C EXPLANATIONS OF THE ALGORITHMS USED IN THIS PACKAGE MAY C BE FOUND IN C C WERNER RHEINBOLDT, C SOLUTION FIELD OF NONLINEAR EQUATIONS AND CONTINUATION METHODS C SIAM JOURNAL OF NUMERICAL ANALYSIS, 17, 1980, PP 221-237 C C COR DEN HEIJER AND WERNER RHEINBOLDT, C ON STEPLENGTH ALGORITHMS FOR A CLASS OF CONTINUATION METHODS, C SIAM JOURNAL OF NUMERICAL ANALYSIS 18, 1981, PP 925-947 C C WERNER RHEINBOLDT, C NUMERICAL ANALYSIS OF CONTINUATION METHODS FOR NONLINEAR C STRUCTURAL PROBLEMS, C COMPUTERS AND STRUCTURES, 13, 1981, PP 103-114 C C*********************************************************************** C C 2. SUBROUTINES USED C C C 2A. USER SPECIFIED SUBROUTINES C C THE FOLLOWING THREE SUBROUTINES MUST BE SPECIFIED BY THE C USER. THE ACTUAL NAMES OF THESE THREE ROUTINES MUST C BE DECLARED EXTERNAL IN THE CALLING PROGRAM, AND PASSED C TO THE CONTINUATION PROGRAM. THE SUBROUTINES ARE LISTED C HERE WITH THE DUMMY EXTERNAL NAMES USED BY THE C CONTINUATION PACKAGE. NOTE THAT A FULL STORAGE SOLVER C IS AVAILABLE IN PITCON, AND WILL BE USED IF THE USER C PASSES THE NAME FSOLVE FOR SLNAME. C C FXNAME EVALUATES THE NVAR-1 COMPONENT FUNCTION FX GIVEN X, AN NVAR C COMPONENT VECTOR. THIS FUNCTION DESCRIBES THE NONLINEAR C SYSTEM. THE AUGMENTING EQUATION IS HANDLED BY THE C CONTINUATION PACKAGE. THE FUNCTION MUST BE DECLARED C EXTERNAL IN THE CALLING PROGRAM, FOR EXAMPLE, C EXTERNAL FCTN, AND MUST HAVE THE FORM - C C SUBROUTINE FCTN(NVAR,X,FX) C REAL X(NVAR),FX(NVAR) C USING INPUT X, EVALUATE NVAR-1 COMPONENTS OF F(X). C C FPNAME EVALUATES THE NVAR-1 BY NVAR JACOBIAN FPRYM (X) AT X C AND STORES IT IN THE NROW BY NCOL ARRAY FPRYM, WITH THE C ACTUAL STORAGE OF THE I,J ENTRY BEING DETERMINED BY THE C SOLVER USED. NOTE THAT THE LAST ROW OF C FPRYM (FOR THE AUGMENTING EQUATION) IS INSERTED BY THE C SOLVING SUBROUTINE. THE JACOBIAN ROUTINE MUST BE DECLARED C EXTERNAL IN THE CALLING PROGRAM. FOR EXAMPLE, C EXTERNAL FPRIME, AND MUST HAVE THE FORM- C C SUBROUTINE FPRIME(NVAR,X,FPRYM,NROW,NCOL) C REAL X(NVAR),FPRYM(NROW,NCOL) C AND MUST STORE THE JACOBIAN ELEMENTS D FI/D XJ C IN THE MATRIX FPRYM USING A FORMAT COMPATIBLE WITH THE C SOLVER CHOSEN. C C SLNAME THE SOLVING ROUTINE FOR THE CONTINUATION CODE. C SLNAME MUST BE DECLARED EXTERNAL IN THE CALLING PROGRAM, C FOR EXAMPLE, EXTERNAL SOLVE. SLNAME MAY BE FSOLVE C IF THE FULL SOLVER SUPPLIED WITH THIS CODE IS USED. C THE SOLVER MUST HAVE THE FORM- C C SUBROUTINE SOLVE(NVAR,X,Y,IHOLD,DET,IEX,IERR,ICALL, C MODNEW,FPRYM,NROW,NCOL,IPIVOT,FPNAME) C EXTERNAL FPNAME C REAL X(NVAR),Y(NVAR),FPRYM(NROW,NCOL) C INTEGER IPIVOT(NVAR) C C THE SOLVER MUST PERFORM THE FOLLOWING OPERATIONS- C C IF (ICALL.EQ.1.OR.MODNEW.EQ.0) REEVALUATE THE JACOBIAN BY C CALLING THE JACOBIAN SUBROUTINE NAMED BY FPNAME. C OTHERWISE, USE THE JACOBIAN AS STORED IN FPRYM. C C MODIFY FPRYM SO THAT THE NVAR-TH ROW IS ALL 0 EXCEPT C THAT THE (NVAR,IHOLD) ENTRY IS 1.0. C C STORE THE DETERMINANT OF THE RESULTING MATRIX IN TWO C PARTS, DET AND IEX, SO THAT DETERMINANT=DET*2**IEX. C DETA SHOULD SATISFY (1.0.LE.ABS(DETA).AND.ABS(DETA).LT.2.0) C C SOLVE FPRYM*Y=Y, OVERWRITING RIGHT HAND SIDE WITH SOLUTION. C SET IERR=0 FOR SATISFACTORY SOLUTION. C SET IERR=1 IF ANY FATAL ERROR WAS DETECTED (IN PARTICULAR, C SINGULARITY). C C IPIVOT IS ONLY USED BY THE SOLVING ROUTINE, AND IS SUPPLIED C FOR PIVOT STORAGE. A USER DESIGNED SOLVER MAY NOT C NEED IPIVOT, IN WHICH CASE THE CODE MAY BE CHANGED TO C DISPENSE WITH IT ENTIRELY. C C C 2B. SUBROUTINES IN THE PITCON PACKAGE ARE C C C PITCON DRIVING ROUTINE OF CONTINUATION CODE. INITIALIZES C INFORMATION, DETERMINES WHETHER LIMIT, TARGET OR C CONTINUATION POINT WILL BE SOUGHT THIS STEP, COMPUTES C STEPLENGTHS, CONTROLS CORRECTOR PROCESS, AND HANDLES ERROR C RETURNS. C C *NOTE* SUBROUTINE PITCON CONTAINS MACHINE DEPENDENT CONSTANTS. C C CORECT USES A FORM OF NEWTON'S METHOD TO SOLVE THE AUGMENTED C NONLINEAR SYSTEM FA(X)=0 WITH AUGMENTING EQUATION C X(IHOLD)=B, THAT IS, X(IHOLD) IS HELD FIXED DURING THE C CORRECTION PROCESS. C C TANGNT SETS UP THE TANGENT SYSTEM DFA(XC,IPL)*TC=E(NVAR), C SOLVING FOR A VECTOR IN THE TANGENT PLANE. THE VECTOR C IS NORMALIZED, AND ITS MAXIMUM ENTRY DETERMINES THE C NEW CONTINUATION PARAMETER UNLESS THE CONTINUATION C PARAMETER HAS BEEN FIXED BY THE USER OR THE PROGRAM. C C ROOT ROOT FINDER USED TO LOCATE LIMIT POINTS. THIS ROUTINE IS A C MODIFIED VERSION OF THE FORTRAN FUNCTION ZERO GIVEN IN C ALGORITHMS FOR MINIMIZATION WITHOUT DERIVATIVES C BY RICHARD P BRENT, PRENTICE HALL, 1973. C C FSOLVE A FULL STORAGE LINEAR EQUATION SOLVER, WHICH MAY BE USED C AS THE EXTERNAL PARAMETER SLNAME WHEN CALLING PITCON. C IT SETS UP AND SOLVES DFA(X,IP)*Y(OUTPUT)=Y(INPUT) C WHERE DFA(X,IP) IS THE JACOBIAN OF FA AT X, C AND Y IS A SUPPLIED RIGHT HAND SIDE. C C *NOTE* SUBROUTINE FSOLVE USES FULL MATRIX STORAGE TO SOLVE THE C SYSTEM. THE USER MAY WISH TO REPLACE THIS ROUTINE WITH ONE C MORE SUITABLE. JUST PASS THE NAME OF THE SOLVER IN SLNAME. C C C 2C. SUBROUTINES FROM THE LINPACK LIBRARY ARE C C C ISAMAX INTEGER FUNCTION RETURNS THE POSITION OF LARGEST ELEMENT OF C THE VECTOR SX. C C SAXPY SETS VECTOR SY(I) = SA*SX(I)+SY(I). C C SCOPY COPIES SY(I)=SX(I). C C SDOT SDOT = SUM(I=1 TO N) SX(I)*SY(I). C C SNRM2 SNRM2 = EUCLIDEAN NORM OF VECTOR SX. C C *NOTE* SNRM2 HAS MACHINE DEPENDENT CUTOFF CONSTANTS. C C SSCAL SETS SX(I)=SA*SX(I). C C SGEFA FACTORS MATRIX A WHOSE LEADING DIMENSION IS LDA C AND WHOSE ACTUAL USED DIMENSION IS N, SETS UP PIVOT C INFORMATION IN VECTOR IPIVOT AND WARNS OF ZERO PIVOTS. C C SGESL ACCEPTS OUTPUT OF SGEFA, AND A RIGHT HAND SIDE B, AND C SOLVES SYSTEM A*X=B, RETURNING X IN B. FOR MODIFIED C NEWTON'S METHOD, ONCE MATRIX IS FACTORED BY SGEFA, ONLY C SGESL IS CALLED FOR SUCCESSIVE RIGHT HAND SIDES. C C *NOTE* IF SUBROUTINE FSOLVE IS REPLACED TO TAKE ADVANTAGE C OF THE SPARSITY OF THE JACOBIAN, THEN THE SUBROUTINES C SGEFA AND SGESL, WHICH ARE ONLY CALLED BY FSOLVE, C ARE NO LONGER NEEDED. C C C LINPACK REFERENCE C C LINPACK USER'S GUIDE, C J J DONGARRA, J R BUNCH, C B MOLER AND G W STEWART, C SOCIETY FOR INDUSTRIAL AND APPLIED MATHEMATICS, C PHILADELPHIA, 1979. C C C 2D. FORTRAN LIBRARY SUBROUTINES USED ARE C C C ABS, ALOG, AMAX1, AMIN1, ATAN2, DBLE, FLOAT, SIGN, SIN, SNGL, SQRT C C*********************************************************************** C C 3. PROGRAMMING NOTES C C THE USER MUST - C C 3A. WRITE AND PASS SUBROUTINES C C THE CALLING PROGRAM MUST DECLARE EXTERNAL THE NAMES TO BE PASSED C AS THE ARGUMENTS FXNAME, FPNAME AND SLNAME. THE FUNCTION AND C JACOBIAN ROUTINES MUST BE WRITTEN BY THE USER. THE SOLVER MUST C BE WRITTEN BY THE USER UNLESS FSOLVE IS C TO BE USED, IN WHICH CASE FSOLVE SHOULD BE PASSED AS C THE SLNAME ARGUMENT. C C 3B. DIMENSION AND PASS STORAGE AREAS C C DECLARE A REAL VECTOR RWORK OF SIZE ISIZE, C WHERE ISIZE.GE.(5*NVAR+NROW*NCOL) (DEPENDS ON SOLVER USED) C AND AN INTEGER VECTOR IPIVOT OF SIZE NVAR. C C 3C. PASS CERTAIN NON-DEFAULTABLE VALUES C C PASS NVAR GREATER THAN ONE, ISIZE.GE.(5*NVAR+NROW*NCOL) AND SET C IRET=0, KSTEP=-1 OR KSTEP=0 ON FIRST CALL FOR A NEW PROBLEM. C SET VALUE OF LUNIT TO FORTRAN OUTPUT UNIT NUMBER. C C 3D. STORE A STARTING POINT C C A POINT XR SHOULD BE STORED IN THE FIRST NVAR ENTRIES OF RWORK C WHICH IS AN APPROXIMATE SOLUTION TO THE EQUATION. THE USER MAY C REQUEST THAT THIS POINT BE IMPROVED TO THE TOLERANCE ABSERR C BY ALSO SETTING KSTEP=-1. C C THE USER SHOULD - C C 3E. RESPOND TO ERROR RETURNS FROM PITCON C C CERTAIN ERROR RETURNS FROM PITCON INDICATE THAT THE RUN SHOULD C BE ABORTED. OTHERS IMPLY THAT SOME INFORMATION RETURNING MAY C BE FAULTY OR UNRELIABLE. SOME ERROR RETURNS MAY POINT AS WELL C TO AN INCORRECT JACOBIAN, POOR STEPSIZE LIMITS OR INAPPROPRIATE C ERROR TOLERANCES. C C THE USER MAY - C C 3F. OBSERVE THE PASSING OF BIFURCATION POINTS C C DETTXL AND DETTXC ARE THE WEIGHTED MANTISSAS OF THE C TANGENT SYSTEMS AT XL AND XC RESPECTIVELY. ON RETURN WITH XF, C IF DETTXL AND DETTXC DIFFER IN SIGN, THEN A SIMPLE BIFURCATION C POINT LIES BETWEEN XL AND XC. DETTXL AND DETTXC ARE AVAILABLE IN C COMMON /DETMAN/. C C 3G. NOTE WORK C C BY EXAMINING THE CONTENTS OF COMMON BLOCKS /COUNT1/ AND /COUNT2/ C THE USER MAY SEE WHERE THE GREATEST WORK IS BEING EXPENDED. C IN PARTICULAR, THIS INFORMATION IS USEFUL WHEN DECIDING WHETHER C TO EMPLOY FULL OR MODIFIED NEWTON'S METHOD ON A SECOND RUN OF C THE PROBLEM. C C 3H. RESET PARAMETERS C C IT IS POSSIBLE TO CHANGE THE PITCON PARAMETERS (ERROR TOLERANCES, C TARGET OR LIMIT POINT INDICES, STEPSIZES AND SO FORTH) IN THE MIDST C OF TRACING A CURVE. IT IS SUGGESTED, HOWEVER, THAT THE CODE BE C WARNED OF THIS CHANGE BY PASSING IN THE VALUE -1 FOR KSTEP C SO THAT CERTAIN AUXILLIARY DATA MAY BE RESET. C C 3I. RUN SEVERAL PROBLEMS C C BECAUSE THE NAMES OF THE FUNCTION, JACOBIAN, AND SOLVER ARE C EXTERNAL PARAMETERS, IT IS FAIRLY EASY TO RUN SEVERAL PROBLEMS C INVOLVING SYSTEMS OF DIFFERENT SIZES OR FORMS, TO SOLVE SOME C SYSTEMS WITH A FULL SOLVER AND OTHERS WITH A BANDED SOLVER, C AND SO FORTH. AGAIN, THE CODE SHOULD BE WARNED OF A SWITCH C TO A NEW PROBLEM BY SETTING KSTEP=-1 BEFORE PROCEEDING. C C C*********************************************************************** C C 4. PITCON PARAMETERS, THEIR DEFINITIONS, TYPES AND DEFAULTS. C C C NVAR = THE NUMBER OF VARIABLES IN THE NONLINEAR SYSTEM. NVAR IS C THE DIMENSION OF THE PIVOT VECTOR IPIVOT, AND OF THE C VECTORS XR, XC, XF, TL AND TC WHICH ARE CONTAINED IN RWORK. C NVAR MUST BE GREATER THAN 1, AND MUST NOT BE CHANGED DURING C THE COURSE OF A PROBLEM RUN. NVAR HAS NO DEFAULT VALUE. C C LIM = LIMIT POINT FLAG AND INDEX. IF (LIM.EQ.0), LIMIT POINTS C ARE NOT TO BE LOOKED FOR. OTHERWISE, THE USER SHOULD SET C LIM TO THE INDEX OF THE VARIABLE FOR WHICH LIMIT C POINTS ARE TO BE SOUGHT. LIM DEFAULTS TO ZERO. C LIM MUST SATISFY 0.LE.LIM.LE.NVAR. C C IT = TARGET POINT FLAG AND INDEX. IF (IT.EQ.0), TARGET POINTS C ARE NOT TO BE LOOKED FOR. OTHERWISE, THE USER SHOULD SET C IT TO THE INDEX OF THE VARIABLE FOR WHICH TARGET C VALUES XIT ARE DESIRED. IT HAS THE DEFAULT VALUE ZERO. C IT MAY BE RESET BY THE USER AT ANY TIME DURING A RUN. C IT MUST SATISFY 0.LE.IT.LE.NVAR. C C XIT = THE VALUE OF THE TARGET VECTOR COMPONENT SOUGHT, IF IT.NE.0. C TARGET POINTS XR SATISFY XR(IT)=XIT. XIT HAS NO DEFAULT. C C KSTEP = THE NUMBER OF CONTINUATION STEPS TAKEN. THIS DOES NOT C INCLUDE FAILURES, OR TARGET AND LIMIT POINTS. THE PROGRAM C INCREMENTS KSTEP EACH TIME A NEW POINT XF IS COMPUTED. C ON THE FIRST CALL TO PITCON FOR A PARTICULAR PROBLEM, THE C USER SHOULD SET KSTEP TO 0 OR -1. IF KSTEP=-1, THE PROGRAM C WILL CHECK THE ACCURACY OF THE STARTING POINT XR, AND IF C NECESSARY, ATTEMPT TO CORRECT IT USING NEWTON'S METHOD. C IF KSTEP=0, THE PROGRAM PERFORMS NO CHECK ON THE STARTING C POINT, AND PROCEEDS TO THE CONTINUATION LOOP. IF THE USER C WISHES TO RUN A DIFFERENT PROBLEM THEN A CALL TO PITCON C WITH KSTEP=-1 OR 0 WILL RESET THE CODE, DESTROYING THE C INFORMATION FROM THE PREVIOUS RUN. KSTEP DEFAULTS TO -1. C C IPC = THE COMPONENT OF XC TO BE USED AS CONTINUATION PARAMETER. C ON INITIAL CALLS WITH KSTEP.EQ.-1 OR KSTEP.EQ.0, THE USER C SUPPLIES A VALUE FOR IPC OR THE DEFAULT IS USED. NOTE C THAT POOR CHOICES OF IPC CAN CAUSE PROGRAM FAILURE. IN C PARTICULAR, ON INITIAL CALL, IF KSTEP.EQ.-1, THE PROGRAM C TRIES TO IMPROVE THE APPROXIMATE INPUT SOLUTION XR TO AN C ACCEPTABLE SOLUTION, WITH THE IPC-TH COMPONENT OF XR C HELD FIXED DURING THE NEWTON ITERATION BUT FOR SOME C VALUES OF IPC, NO SUCH SOLUTION MAY EXIST. FURTHERMORE C IF KSTEP.EQ.0, THE PROGRAM SOLVES A TANGENT SYSTEM ASSUMING C THAT THE TANGENT AT THE POINT HAS A NONZERO IPC-TH C COMPONENT. A POOR CHOICE OF IPC MAY THEN CAUSE THE TANGENT C CALCULATION TO FAIL. IPC DEFAULTS TO NVAR. C C IPCFIX= A FLAG WHICH CAN FORCE THE CONTINUATION PARAMETER IPC TO C REMAIN UNCHANGED. IF (IPCFIX.NE.1) THE PROGRAM IS FREE C TO CHOOSE THE CONTINUATION PARAMETER, GENERALLY C AS THE INDEX OF THE MAXIMUM ENTRY OF THE TANGENT. C IF (IPCFIX.EQ.1) THE INPUT VALUE OF IPC IS USED FOR THE C CONTINUATION PARAMETER ON THE CURRENT STEP. NOTE THAT THE C CHOICE IPCFIX=1 ASSERTS THAT THE CURVE MAY BE GLOBALLY C PARAMETERIZED BY A SINGLE FIXED VARIABLE. THE PROGRAM C WILL FAIL IF THIS PARAMETERIZATION DETERIORATES OR C FAILS, PARTICULARLY AT A LIMIT POINT IN THE PARAMETER. C IPCFIX DEFAULTS TO 0. C C DIRIPC= LOCAL CONTINUATION DIRECTION WITH RESPECT TO CURRENT C CONTINUATION PARAMETER IPC. DIRIPC IS PLUS OR MINUS C ONE, DENOTING CURRENT CONTINUATION IN THE DIRECTION C OF INCREASING OR DECREASING X(IPC) RESPECTIVELY. C ON INITIAL CALL, DIRIPC SHOULD BE PASSED BY THE USER. C THEREAFTER, DIRIPC IS COMPUTED BY THE CODE TO PRESERVE C THE SENSE IN WHICH THE CURVE IS BEING FOLLOWED. C DIRIPC DEFAULTS TO +1.0. C C HTANCF= STEPSIZE TO BE USED ON NEXT PREDICTOR STEP ALONG THE C UNIT TANGENT. ON FIRST CALL, THE USER SHOULD SUPPLY A C VALUE OF HTANCF. THEREAFTER, THE VALUE OF HTANCF IS C DETERMINED BY THE PROGRAM, SUBJECT TO THE CONSTRAINTS C IMPOSED BY HMAX, HMIN, AND HFACT, AS WELL AS INFORMATION C FROM RECENT STEPS TAKEN BY THE PROGRAM. HTANCF SHOULD C BE POSITIVE. HTANCF DEFAULTS TO .5*(HMAX+HMIN). C C IRET = A RETURN FLAG TO INDICATE ERRORS OR THE TYPE OF POINT C RETURNED IN XR. NONNEGATIVE VALUES OF IRET REPRESENT C NORMAL RETURNS. NEGATIVE VALUES OF IRET INDICATE THAT C SOME ERROR OR DIFFICULTY HAS BEEN ENCOUNTERED. VALUES C OF IRET BETWEEN -1 AND -4 ARE SIMPLY REPORTS THAT AN C ATTEMPT TO COMPUTE A LIMIT OR TARGET POINT FAILED. THESE C DO NOT AFFECT FURTHER CONTINUATION STEPS, AND THE USER NEED C NOT MODIFY ANY VARIABLES BEFORE PROCEEDING. C VALUES OF IRET OF -5 AND -6 REFER TO DANGEROUS SITUATIONS C THAT MAY BE CORRECTABLE. C VALUES OF IRET FROM -7 TO -10 ARE SERIOUS, FATAL ERRORS. C THE USER SHOULD HALT THE PROGRAM AND EXAMINE THE INPUT C AND THE INTERIM RESULTS. C IRET SHOULD BE ZERO ON THE FIRST CALL FOR A PROBLEM. C THE SPECIFIC VALUES OF IRET AND THEIR MEANINGS ARE C C IRET=2 NORMAL RETURN WITH LIMIT POINT IN XR AND TANGENT C AT XR CONTAINED IN TL. C C IRET=1 NORMAL RETURN WITH TARGET POINT IN XR. C C IRET=0 NORMAL RETURN WITH NEW CONTINUATION POINT IN XR. C C IRET=-1 AN ERROR OCCURRED DURING COMPUTATION OF C A LIMIT POINT. C C IRET=-2 CORECT CALLED FOR TARGET POINT CALCULATION FAILED C TO CONVERGE IN MAXCOR ITERATIONS. C C IRET=-3 THE SOLVER WAS CALLED BY CORECT FOR TARGET POINT C CALCULATION, AND FAILED. C C IRET=-4 UNACCEPTABLE CORRECTOR STEP IN TARGET POINT C CALCULATION. C C IRET=-5 PREDICTION STEP HTANCF IS LESS THAN HMIN, PERHAPS C BECAUSE OF REPEATED FAILURE OF CORRECTOR, AND C CONSEQUENT STEPSIZE REDUCTION. USER MIGHT REDUCE C HMIN, OR SWITCH FROM MODCON=1 TO MODCON=0, OR C INCREASE ABSERR AND RELERR. BUT BE AWARE THAT C REPEATED STEPSIZE REDUCTIONS MAY INDICATE AN C INTRACTABLE FUNCTION OR AN INCORRECT JACOBIAN. C C IRET=-6 FUNCTION VALUE FNRM OF INPUT XR IS TOO LARGE C AND COULD NOT BE IMPROVED BY CORRECTOR. USER C SHOULD CHECK IF THE VALUE OF IPC IS APPROPRIATE, C OR WHETHER ABSERR AND RELERR SHOULD BE INCREASED, C OR WHETHER THE INPUT POINT XR MIGHT BE IMPROVED C BEFORE CALLING PROGRAM, AS WELL AS WHETHER THE C FUNCTION AND JACOBIAN EVALUATORS ARE CORRECT. C C IRET=-7 THE SOLVER FAILED IN A CALL FROM TANGNT. C THIS RETURN MAY MEAN THAT THE PROGRAM'S CHOICE C OF IPC WAS FAULTY, IN WHICH CASE A RESTART C WITH A NEW VALUE OF IPC AND KSTEP.EQ.0 MAY C ALLOW THE USER TO RECOVER. C C IRET=-8 THE SOLVER FAILED IN A CALL FROM CORECT. C THIS RETURN MAY MEAN THAT THE JACOBIAN IS C TOO ILL CONDITIONED FOR THE SOLVER, OR THAT C A SINGULARITY OF THE JACOBIAN IS NEARBY. C C IRET=-9 THE TANGENT VECTOR TC AT XC IS ZERO. C THIS RETURN INDICATES AN ERROR IN THE ROUTINE C WHICH SOLVES THE LINEAR SYSTEMS. C C IRET=-10 IMPROPER INPUT, NVAR.LE.1, OR C ISIZE.LT.5*NVAR+NROW*NCOL, OR C PROGRAM HAS BEEN CALLED AGAIN AFTER FATAL ERROR. C C MODCON= FLAG FOR TYPE OF NEWTON METHOD TO BE USED IN CORRECTING C CONTINUATION AND TARGET POINTS. (NOTE THAT FOR LIMIT C POINTS, MODIFIED NEWTON'S METHOD IS TRIED FIRST, AND C FULL NEWTON'S METHOD IS APPLIED ONLY IF THE MODIFIED C METHOD DOES NOT CONVERGE.) C C MODCON=0 UPDATE JACOBIAN FOR TANGENT CALCULATION, C UPDATE JACOBIAN FOR EACH CORRECTION STEP C OF CONTINUATION POINTS. C C MODCON=1 UPDATE JACOBIAN FOR TANGENT CALCULATION, C EVALUATE JACOBIAN ONLY AT FIRST CORRECTION STEP C OF CONTINUATION POINTS. C C IPIVOT= INTEGER VECTOR USER DECLARED TO BE OF SIZE NVAR, C USED DURING THE MATRIX FACTORIZATION TO STORE PIVOT C INFORMATION. C C HMAX = THE MAXIMUM PREDICTOR STEP SIZE ALLOWED. HTANCF WILL C ALWAYS BE FORCED TO BE LESS THAN OR EQUAL TO HMAX. C IF (HMAX.LE.0.0) ON INITIAL CALL, THE DEFAULT VALUE IS USED. C HMAX DEFAULTS TO SQRT(NVAR). C C HMIN = THE MINIMUM PREDICTOR STEP SIZE ALLOWED. HTANCF WILL C ALWAYS BE FORCED TO BE GREATER THAN OR EQUAL TO HMIN. C IF HTANCF FALLS BELOW HMIN BECAUSE OF STEP SIZE REDUCTIONS C CAUSED BY FAILURE OF NEWTON'S METHOD, AN ERROR RETURN IS C FORCED. IF (HMIN.LE.EPSQRT) ON INITIAL CALL, THE DEFAULT C VALUE IS USED. THE DEFAULT VALUE OF HMIN IS C EPSQRT=SQRT(EPMACH) WHERE EPMACH IS THE MACHINE PRECISION C CONSTANT, THE SMALLEST NUMBER X SO THAT COMPUTATIONALLY, C (1.0.LT.(1.0+X)). C C HFACT = MAXIMUM GROWTH FACTOR FOR PREDICTOR STEP SIZE HTANCF C BASED ON PREVIOUS SECANT STEP SIZE HSECLC. THE BOUNDS C SATISFIED ARE (HSECLC/HFACT).LE.HTANCF.LE.(HSECLC*HFACT). C IF THE CORRECTOR STEP FAILS, THEN HFACT IS ALSO USED TO C REDUCE THE PREDICTOR STEP HTANCF TO (HTANCF/HFACT). C IF (HFACT.LE.1.0) ON INITIAL CALL, HFACT DEFAULTS TO 3.0. C C ABSERR= ABSOLUTE ERROR CONTROL. A POINT X IS ACCEPTED AS A C SOLUTION OF FX=0 IF THE LARGEST COMPONENT C OF F(X) IS STILL LESS THAN ABSERR IN ABSOLUTE VALUE. C NOTE THAT POINTS WHICH DO NOT SATISFY THIS CRITERION C MAY NONETHELESS BE ACCEPTED BY THE CORRECTOR. C IF (ABSERR.LE.0.0) ON FIRST CALL, ABSERR IS SET TO C THE DEFAULT VALUE OF SQRT(EPMACH). C C RELERR= RELATIVE ERROR CONTROL. DURING THE CORRECTOR ITERATION, C A CONVERGENCE TEST IS APPLIED TO THE SIZE OF THE CORRECTOR C STEPS. RELERR IS USED IN COMPUTING THE STEPSIZE TOLERANCE C TEST=ABSERR+RELERR*XNRM, WHERE XNRM IS THE MAXIMUM NORM C OF THE CURRENT ITERATE. NOTE THAT STEPS MAY BE ACCEPTED C WHICH DO NOT SATISFY THIS CRITERION. IF (RELERR.LE.0.0) C ON FIRST CALL, RELERR IS SET TO THE DEFAULT VALUE. C RELERR DEFAULTS TO SQRT(EPMACH) WHERE EPMACH IS THE C MACHINE PRECISION CONSTANT. C C RWORK = USER DECLARED VECTOR OF SIZE ISIZE=(5*NVAR+NROW*NCOL). C RWORK STORES FIVE VECTORS AND AN ARRAY IN THE ORDER C (XR,XC,XF,TL,TC,FPRYM). THEIR BEGINNING LOCATIONS C ARE IXR=1, IXC=NVAR+1, IXF=2*NVAR+1, ITL=3*NVAR+1, C ITC=4*NVAR+1, IFP=5*NVAR+1. FPRYM IS AN ARRAY OF C SIZE NROW X NCOL. THE MEANINGS OF THESE COMPONENTS OF C RWORK ARE DESCRIBED BELOW. THE USER SHOULD SET XR C ON FIRST CALL, BUT NO OTHER PORTIONS OF RWORK SHOULD BE C SET. AFTER THE FIRST CALL, NO ENTRIES OF RWORK C SHOULD BE ALTERED. C C (XR) = REAL VECTOR OF LENGTH NVAR, FIRST NVAR ENTRIES OF RWORK. C ON FIRST CALL, A USER SUPPLIED STARTING POINT, WHICH MAY BE C IMPROVED BY THE PROGRAM IF KSTEP=-1. ON NORMAL RETURN FROM C PITCON, XR WILL HOLD THE MOST RECENTLY FOUND POINT, WHETHER C A CONTINUATION POINT, TARGET POINT, OR LIMIT POINT. C C (XC) = REAL VECTOR OF LENGTH NVAR, ENTRIES NVAR+1 TO 2*NVAR OF C RWORK. THE PREVIOUS CONTINUATION POINT. C C (XF) = REAL VECTOR OF LENGTH NVAR, ENTRIES 2*NVAR+1 TO 3*NVAR OF C RWORK. THE CURRENT CONTINUATION POINT. C C (TL) = REAL VECTOR OF LENGTH NVAR, ENTRIES 3*NVAR+1 TO 4*NVAR OF C RWORK. THE CONTENTS OF TL ON RETURN VARY. C ON RETURN WITH CONTINUATION POINT, TL CONTAINS THE PREVIOUS C VALUE OF TC, A TANGENT VECTOR CORRESPONDING TO THE C CONTINUATION POINT XL WHICH IS NOT STORED. C ON RETURN WITH A TARGET POINT, TL CONTAINS THE FUNCTION C VALUES AT THE TARGET POINT. C ON RETURN WITH A LIMIT POINT, TL CONTAINS THE VALUE C OF THE TANGENT VECTOR AT THE LIMIT POINT. C C (TC) = REAL VECTOR OF LENGTH NVAR, ENTRIES 4*NVAR+1 TO 5*NVAR OF C RWORK. THE TANGENT VECTOR AT THE PREVIOUS CONTINUATION C POINT XC. C C (FPRYM)= REAL ARRAY OF DIMENSIONS (NROW,NCOL). ENTRIES 5*NVAR+1 TO C 5*NVAR+NROW*NCOL OF RWORK. FPRYM IS USED TO STORE THE C AUGMENTED JACOBIAN OF THE SYSTEM. THE METHOD OF STORAGE C USED BY THE JACOBIAN ROUTINE MUST BE COMPATIBLE WITH THE C METHOD OF RETRIEVAL AND SOLUTION USED BY THE SOLVER. C IF USING THE PREPROGRAMMED ROUTINE FSOLVE, THE C MATRIX FPRYM IS SQUARE OF DIMENSIONS (NVAR,NVAR) AND C FPRYM(I,J)=D FI(X)/D X(J). C C ISIZE = USER SET DIMENSION FOR VECTOR RWORK, WHICH MUST BE AT C LEAST OF SIZE (5*NVAR+NROW*NCOL). C FULL MATRIX STORAGE USING SLNAME=FSOLVE REQUIRES NROW=NVAR, C NCOL=NVAR. C C NROW = NUMBER OF ROWS OF MATRIX STORAGE TO BE USED IN RWORK. C NROW=NVAR FOR FULL MATRIX STORAGE AND SLNAME=FSOLVE. C C NCOL = NUMBER OF COLUMNS OF MATRIX STORAGE TO BE USED IN RWORK. C NCOL=NVAR FOR FULL MATRIX STORAGE AND SLNAME=FSOLVE. C C FXNAME= EXTERNAL QUANTITY, THE NAME OF THE FUNCTION ROUTINE C WHICH WAS DECLARED EXTERNAL IN THE CALLING PROGRAM. C FOR DETAILS OF THE FORMAT OF THE ROUTINE, SEE PART 2A. C C FPNAME= EXTERNAL QUANTITY, THE NAME OF THE DERIVATIVE ROUTINE C WHICH WAS DECLARED EXTERNAL IN THE CALLING PROGRAM. C FOR DETAILS ON THE FORMAT OF THIS ROUTINE, SEE PART 2A. C C SLNAME= EXTERNAL QUANTITY, THE NAME OF THE SOLVING ROUTINE C WHICH WAS DECLARED EXTERNAL IN THE CALLING PROGRAM C FOR DETAILS ON THE FORMAT OF THIS ROUTINE, SEE PART 2A. C C LUNIT = INTEGER VARIABLE WHICH CONTAINS THE LOGICAL OUTPUT UNIT C FOR ERROR MESSAGES AND DIAGNOSTICS FROM THE PROGRAM. C ALL WRITE STATEMENTS IN THE PROGRAM HAVE THE FORM C WRITE(LUNIT,FORMAT NUMBER)QUANTITIES C C*********************************************************************** C C 5. LABELED COMMON BLOCKS C C /COUNT1/ COUNTS THE NUMBER OF CALLS BY ONE ROUTINE OF ANOTHER. C C ICRSL = NUMBER OF CALLS BY CORRECTOR TO SOLVER. C ITNSL = NUMBER OF CALLS BY TANGNT TO SOLVER. C NSTCR = NUMBER OF CORRECTOR STEPS TAKEN FOR IMPROVING C STARTING POINTS. C NCNCR = NUMBER OF CORRECTOR STEPS TAKEN FOR COMPUTING C CONTINUATION POINTS. C NTRCR = NUMBER OF CORRECTOR STEPS TAKEN FOR COMPUTING C TARGET POINTS. C NLMCR = NUMBER OF CORRECTOR STEPS TAKEN FOR COMPUTING C LIMIT POINTS. C NLMRT = NUMBER OF ROOT FINDING STEPS TAKEN FOR COMPUTING C LIMIT POINTS. C C /COUNT2/ KEEPS PERFORMANCE AND WORK STATISTICS. C C IFEVAL = NUMBER OF CALLS TO FUNCTION EVALUATOR. C IPEVAL = NUMBER OF CALLS TO JACOBIAN EVALUATOR. C ISOLVE = NUMBER OF CALLS TO SOLVER. C NREDXF = NUMBER OF STEPSIZE REDUCTIONS MADE BEFORE PREDICTOR C POINT CONVERGED TO THE NEW CONTINUATION POINT. C NRDSUM = TOTAL NUMBER OF STEPSIZE REDUCTIONS. C NCORXR = NUMBER OF CORRECTOR ITERATION STEPS TAKEN IN MOST RECENT C CALL TO CORECT. C NCRSUM = TOTAL NUMBER OF CORRECTOR ITERATION STEPS. C C /DETEXP/ CONTAINS THE BINARY EXPONENTS OF RECENT DETERMINANTS. C LET DET1 DENOTE THE WEIGHTED MANTISSA OF THE DETERMINANT C AT X1 AND IEX1 DENOTE THE BINARY EXPONENT OF THE C DETERMINANT AT X1. THEN THE DETERMINANT MAY BE C RECOVERED AS DETERMINANT=DET1*2.0**IEX1. C C IEXTXL = EXPONENT FOR DETERMINANT OF TANGENT SYSTEM AT XL. C IEXTXC = EXPONENT FOR DETERMINANT OF TANGENT SYSTEM AT XC. C IEXTXR = EXPONENT FOR DETERMINANT OF TANGENT SYSTEM AT XR C WHERE XR IS A LIMIT POINT. C IEXNXR = EXPONENT FOR DETERMINANT OF NEWTON SYSTEM WHICH C PRODUCED XR, WHERE XR MAY BE XF OR A TARGET OR C LIMIT POINT. C C /DETMAN/ CONTAINS THE MANTISSAS OF RECENT DETERMINANTS. C C DETTXL = MANTISSA FOR DETERMINANT OF TANGENT SYSTEM AT XL. C DETTXC = MANTISSA FOR DETERMINANT OF TANGENT SYSTEM AT XC. C DETTXR = MANTISSA FOR DETERMINANT OF TANGENT SYSTEM AT XR C WHERE XR IS A LIMIT POINT. C DETNXR = MANTISSA FOR DETERMINANT OF NEWTON SYSTEM WHICH C PRODUCED XR, WHERE XR MAY BE XF OR A TARGET OR C LIMIT POINT. C C /OUTPUT/ C C IWRITE = USER ACCESSIBLE OUTPUT INDICATOR. C IWRITE=0, NO OUTPUT PRINTED BY PITCON. C IWRITE=1, ERROR MESSAGES PRINTED BY PITCON. C IWRITE=2, CERTAIN OUTPUT WILL BE PRINTED BY PITCON. C C /POINT/ CONTAINS DATA ABOUT THE SOLUTION CURVE. C C CURVCF = ESTIMATED CURVATURE BETWEEN XC AND XF. C CORDXF = NORM OF THE CORRECTOR STEP FROM PREDICTED POINT TO C CORRECTED POINT, WITH MAXIMUM ABSOLUTE VALUE AS THE NORM. C THIS QUANTITY IS MODIFIED TO AN OPTIMAL VALUE. C ALPHLC = ANGLE BETWEEN OLD AND NEW TANGENTS TL AND TC C HSECLC = EUCLIDEAN NORM OF SECANT BETWEEN XL AND XC. C FNRM = MAXIMUM NORM OF FUNCTION VALUE AT CURRENT POINT XR. C C /TOL/ CONTAINS MACHINE DEPENDENT TOLERANCES. C THE QUANTITY EPMACH IS THE MACHINE PRECISION CONSTANT OR C SMALLEST RELATIVE SPACING. THE VALUES FOR EPMACH IN THE C CODE FOR DIFFERENT MACHINES ARE TAKEN FROM THE BELL LABS C PORT FUNCTION R1MACH AND CORRESPOND TO R1MACH(3). C C EPMACH= SMALLEST NUMBER SUCH THAT 1.0+EPMACH.GT.EPMACH C .5*BETA**(1-TAU) FOR ROUNDED, TAU-DIGIT ARITHMETIC C BASE BETA. TWICE THIS VALUE FOR TRUNCATED ARITHMETIC. C THIS IS THE RELATIVE MACHINE PRECISION. C EPMACH=2**(-27) FOR DEC-10. C EPSATE= 8*EPMACH. C EPSQRT= SQUARE ROOT OF EPMACH. C C REFERENCE FOR MACHINE CONSTANTS- C C MACHINE CONSTANTS FOR PORTABLE FORTRAN LIBRARIES, C PHYLLIS FOX, A D HALL, N L SCHRYER, C COMPUTING SCIENCE TECHNICAL REPORT 37, C BELL LABORATORIES, C MURRAY HILL, NEW JERSEY, AUGUST 1975. C C*********************************************************************** C C 6. NOMENCLATURE FOR STEP DEPENDENT VARIABLES C C C THE PROGRAM ACCUMULATES INFORMATION THAT IS ASSOCIATED WITH C SEVERAL PREVIOUS CONTINUATION POINTS OR THE STEPS MADE BETWEEN C THEM. IN INTERPRETING THE CODE OR ITS OUTPUT, IT IS IMPORTANT C TO KNOW WHERE SUCH QUANTITIES APPLY. C C THE POINTS XLL AND XL WILL HAVE BEEN DISCARDED BY THE PROGRAM, C BUT SOME QUANTITIES ASSOCIATED WITH THEM STILL SURVIVE. C C QUANTITIES ASSOCIATED WITH STEP FROM XLL TO XL C C HSECLL = SIZE OF SECANT FROM XLL TO XL, EUCLIDEAN NORM(XLL-XL) C C QUANTITIES ASSOCIATED WITH THE POINT XL C C DETTXL = BINARY MANTISSA OF DETERMINANT OF DFA(XL,IPLL), DIVIDED C BY IPLL-TH COMPONENT OF TANGENT AT XL. C IEXTXL = BINARY EXPONENT OF DETERMINANT OF TANGENT SYSTEM AT XL. C IPL = THE LOCATION OF THE FIRST OR SECOND LARGEST COMPONENT C OF THE TANGENT VECTOR AT XL. C TLLIM = VALUE OF LIM-TH COMPONENT OF TANGENT VECTOR AT XL. C TL = TANGENT VECTOR AT XL, ALTHOUGH LIMIT OR TARGET POINT C CALCULATIONS COULD HAVE OVERWRITTEN THIS VECTOR. C C QUANTITIES ASSOCIATED WITH INTERVAL FROM XL TO XC C C ALPHLC = THE ANGLE BETWEEN THE TANGENTS AT XL AND XC. C CURVLC = ESTIMATED CURVATURE BETWEEN XL AND XC. C HSECLC = SIZE OF SECANT BETWEEN XL AND XC, EUCLIDEAN NORM(XL-XC). C C QUANTITIES ASSOCIATED WITH THE POINT XC C C DETTXC = BINARY MANTISSA OF DETERMINANT OF TANGENT SYSTEM AT XC. C DIRIPC = SIGN OF DETTXC, DETERMINES SENSE OF CONTINUATION. C IEXTXC = BINARY EXPONENT OF DETERMINANT OF TANGENT SYSTEM AT XC. C IPC = LOCATION OF FIRST OR SECOND LARGEST COMPONENT OF TANGENT C VECTOR AT XC. C TC = TANGENT VECTOR AT XC. C TCLIM = VALUE OF LIM-TH COMPONENT OF TANGENT AT XC. C TCIPC = VALUE OF TC(IPC). C C QUANTITIES ASSOCIATED WITH THE INTERVAL FROM XC TO XF C C CURVCF = ESTIMATED CURVATURE BETWEEN XC AND XF. C HTANCF = STEPSIZE USED ALONG TANGENT TO GET PREDICTED POINT C WHICH WAS CORRECTED TO SOLUTION POINT XF. C C QUANTITIES ASSOCIATED WITH THE POINT XF C C CORDXF = SIZE OF THE TOTAL CORRECTION FROM PREDICTED POINT C X=XC+HTANCF*TC TO CORRECTED POINT XF. C NOTE THAT THIS HAS BEEN MODIFIED TO AN OPTIMAL VALUE. C CURVXF = A PREDICTED VALUE OF THE CURVATURE AT XF. C NREDXF = NUMBER OF STEPSIZE REDUCTIONS MADE BEFORE THE CORRECTOR C ITERATION WOULD CONVERGE TO XF. C C QUANTITIES ASSOCIATED WITH THE POINT XR, WHICH MAY BE THE CURRENT C CONTINUATION POINT XF, OR A LIMIT OR TARGET POINT C C DETNXR = BINARY MANTISSA OF DETERMINANT OF NEWTON SYSTEM WHICH C CONVERGED TO XR. C DETTXR = BINARY MANTISSA OF DETERMINANT OF TANGENT SYSTEM AT XR, C IF XR IS A LIMIT POINT. C FNRM = MAXIMUM NORM OF FUNCTION VALUE AT XR. C FPRYM = CONTAINS THE FACTORED NEWTON SYSTEM AT THE PENULTIMATE C NEWTON ITERATE FOR CONTINUATION OR TARGET POINT RETURNS, C OR THE FACTORED TANGENT SYSTEM FOR LIMIT POINT RETURNS. C IEXNXR = BINARY EXPONENT OF DETERMINANT OF NEWTON SYSTEM WHICH C CONVERGED TO XR. C IEXTXR = BINARY EXPONENT OF DETERMINANT OF TANGENT SYSTEM AT XR, C IF XR IS A LIMIT POINT. C NCORXR = NUMBER OF NEWTON STEPS TAKEN IN LAST CORRECTOR ITERATION. C XSTEP = SIZE OF THE LAST CORRECTOR STEP TAKEN IN CONVERGING TO XR. C C*********************************************************************** C EXTERNAL FXNAME, FPNAME, SLNAME INTEGER IPIVOT(NVAR) REAL RWORK(ISIZE) REAL WRGE(8), ACOF(12) DOUBLE PRECISION DTLIPC, DTCIPC, DADJUS, COSALF COMMON /COUNT1/ ICRSL, ITNSL, NSTCR, NCNCR, NTRCR, NLMCR, NLMRT COMMON /COUNT2/ IFEVAL, IPEVAL, ISOLVE, NREDXF, NRDSUM, NCORXR, * NCRSUM COMMON /DETEXP/ IEXTXL, IEXTXC, IEXTXR, IEXNXR COMMON /DETMAN/ DETTXL, DETTXC, DETTXR, DETNXR COMMON /OUTPUT/ IWRITE COMMON /POINT/ CURVCF, CORDXF, ALPHLC, HSECLC, FNRM COMMON /TOL/ EPMACH, EPSATE, EPSQRT DATA IDONE /0/ DATA TENM1 /0.1/ DATA TENM2 /0.01/ DATA TENM3 /0.001/ C C*********************************************************************** C C RMACH3 IS A MACHINE DEPENDENT QUANTITY. C THE APPROPRIATE DATA STATEMENT SHOULD BE ACTIVATED. C C*********************************************************************** C C RMACH3 FOR HONEYWELL 6000 SERIES FOLLOWS C C DATA RMACH3 /0714400000000/ C C RMACH3 FOR IBM 360 AND 370 AND SEL SYSTEMS 85/86 FOLLOWS C C DATA RMACH3 /Z3B100000/ C C RMACH3 FOR CDC 6000 AND 7000 SERIES FOLLOWS C C DATA RMACH3 /16404000000000000000B/ C C RMACH3 FOR PDP-10 KA OR KI PROCESSOR FOLLOWS C DATA RMACH3 /"146400000000/ C C RMACH3 FOR THE PDP-11 FOLLOWS C C DATA RMACH3 /0.59604644775E-07/ C C RMACH3 FOR THE UNIVAC 1100 SERIES FOLLOWS C C DATA RMACH3 /0146400000000/ C C*********************************************************************** C C 1. PREPARATIONS. C ON FIRST CALL FOR THIS PROBLEM, INITIALIZE COUNTERS AND VARIABLES, C CHECK USER INFORMATION AND SET DEFAULTS, AND IF (KSTEP.EQ.-1), C CHECK NORM OF F(XR) AND CORRECT XR IF NECESSARY. C ON EACH CALL, IF INPUT IRET HAS NONFATAL VALUE, RESET IRET C SO THAT CONTINUATION LOOP PICKS UP WHERE IT WAS HALTED. C C*********************************************************************** C IERR = 0 IF (IRET.EQ.(-1)) IRET = 2 IF (IRET.EQ.(-2) .OR. IRET.EQ.(-3) .OR. IRET.EQ.(-4)) IRET = 1 IF (IRET.EQ.(-5) .OR. IRET.EQ.(-6)) IRET = 0 C C IF CODE WAS CALLED AGAIN AFTER FATAL VALUE OF IRET, C THEN RETURN WITH ERROR VALUE IRET=-10. C IF (IRET.LT.0) GO TO 470 C C PERFORM ONE-TIME ONLY INITIALIZATIONS C IF (IDONE.NE.0) GO TO 10 WRGE(1) = 0.8735115 WRGE(2) = 0.1531947 WRGE(3) = 0.03191815 WRGE(4) = 0.3339946E-10 WRGE(5) = 0.4677788 WRGE(6) = 0.6970123E-03 WRGE(7) = 0.1980863E-05 WRGE(8) = 0.1122789E-08 ACOF(1) = 0.9043128 ACOF(2) = -0.7075675 ACOF(3) = -4.667383 ACOF(4) = -3.677482 ACOF(5) = 0.8516099 ACOF(6) = -0.1953119 ACOF(7) = -4.830636 ACOF(8) = -0.9770528 ACOF(9) = 1.040061 ACOF(10) = 0.03793395 ACOF(11) = 1.042177 ACOF(12) = 0.04450706 C C SET THE MACHINE DEPENDENT VARIABLE EPMACH, THE SMALLEST NUMBER C SO THAT (1.0+EPMACH.GT.1.0) C EPMACH = RMACH3 C C SET EPSATE=8*EPMACH, EPSQRT=SQRT(EPMACH) C EPSATE = 8.0*EPMACH EPSQRT = SQRT(EPMACH) TCOS = 1.0 - EPMACH TSIN = SQRT(1.0-TCOS*TCOS) ALFMIN = 2.0*ATAN2(TSIN,TCOS) IF (KSTEP.LT.(-1) .OR. KSTEP.GT.0) KSTEP = -1 KSTEPL = -2 IDONE = 1 C C PERFORM INITIALIZATIONS AND CHECKS FOR NEW PROBLEM ONLY C 10 IF (KSTEP.GT.0) GO TO 30 IF (KSTEPL.EQ.(-1) .AND. KSTEP.EQ.0) GO TO 20 IF (NVAR.LE.1) GO TO 470 IF (ISIZE.LT.(5*NVAR+NROW*NCOL)) GO TO 470 IXR = 1 JXR = 0 IXC = IXR + NVAR JXC = JXR + NVAR IXF = IXC + NVAR JXF = JXC + NVAR ITL = IXF + NVAR JTL = JXF + NVAR ITC = ITL + NVAR JTC = JTL + NVAR IFP = ITC + NVAR JFP = JTC + NVAR DETTXL = 0.0 DETTXC = 0.0 DETTXR = 0.0 DETNXR = 0.0 IEXTXL = 0 IEXTXC = 0 IEXTXR = 0 IEXNXR = 0 TCIPC = 0.0 CORDXF = 0.0 CURVCF = 0.0 HSECLL = 0.0 HSECLC = 0.0 XITO = 0.0 ITO = 0 NEQN = NVAR - 1 NREDXF = 0 NCRSUM = 0 NRDSUM = 0 ICRSL = 0 ITNSL = 0 NSTCR = 0 NCNCR = 0 NTRCR = 0 NLMCR = 0 NLMRT = 0 IFEVAL = 0 ISOLVE = 0 IPEVAL = 0 IF (HMIN.LE.EPSQRT) HMIN = EPSQRT IF (HMAX.LE.0.0) HMAX = SQRT(FLOAT(NVAR)) HDEF = .5*(HMAX+HMIN) IF (HFACT.LE.1.0) HFACT = 3.0 HRED = 1.0/HFACT IF (ABSERR.LE.0.0) ABSERR = EPSQRT IF (RELERR.LE.0.0) RELERR = EPSQRT IF (IPCFIX.NE.1) IPCFIX = 0 IF (IPC.LE.0 .OR. IPC.GT.NVAR) IPC = NVAR IF (LIM.LT.0 .OR. LIM.GT.NVAR) LIM = 0 IF (HTANCF.LE.0.0) HTANCF = HDEF IF (DIRIPC.NE.(-1.0)) DIRIPC = 1.0 C C IF (KSTEP.LT.0) CHECK NORM OF F(XR) AT STARTING POINT, C IF ACCEPTABLE, RETURN IMMEDIATELY WITH KSTEP=0, C OTHERWISE APPLY NEWTON'S METHOD, HOLDING VALUE OF C IPC-TH COMPONENT FIXED. C IF (KSTEP.GE.0) GO TO 20 CALL CORECT(NVAR, RWORK(IXR), IPC, RWORK(ITL), IERR, MODCON, * RWORK(IFP), NROW, NCOL, IPIVOT, ABSERR, RELERR, XSTEP, NEQN, * FNRM, FXNAME, FPNAME, SLNAME, LUNIT, DETNXR, IEXNXR) NSTCR = NSTCR + NCORXR C C IF NO ACCEPTABLE POINT FOUND, ERROR RETURN C IF (IERR.NE.0) GO TO 430 KSTEPL = -1 KSTEP = 0 GO TO 370 20 IF (KSTEP.EQ.0) CALL SCOPY(NVAR, RWORK(IXR), 1, RWORK(IXC), 1) IF (KSTEP.EQ.0) CALL SCOPY(NVAR, RWORK(IXR), 1, RWORK(IXF), 1) C C*********************************************************************** C C 2. TARGET POINT CHECK. IF (IT.NE.0) TARGET POINTS ARE SOUGHT. C CHECK TO SEE IF TARGET COMPONENT IT HAS VALUE XIT LYING C BETWEEN XC(IT) AND XF(IT). IF SO, GET LINEARLY INTERPOLATED C STARTING POINT, AND USE NEWTON'S METHOD TO GET TARGET POINT C C*********************************************************************** C 30 IF (IT.LT.0 .OR. IT.GT.NVAR) IT = 0 IF (IT.EQ.0) GO TO 40 IF (IRET.EQ.1 .AND. XIT.EQ.XITO .AND. IT.EQ.ITO) GO TO 40 INDEX = JXC + IT XCIT = RWORK(INDEX) INDEX = JXF + IT XFIT = RWORK(INDEX) IF ((XIT.LE.XCIT) .AND. (XIT.LT.XFIT)) GO TO 40 IF ((XIT.GE.XCIT) .AND. (XIT.GT.XFIT)) GO TO 40 DEL = XFIT - XCIT RAT = 0.0 IF (ABS(DEL).GT.EPSQRT) RAT = (XIT-XCIT)/DEL CALL SCOPY(NVAR, RWORK(IXF), 1, RWORK(IXR), 1) CALL SAXPY(NVAR, -1.0, RWORK(IXC), 1, RWORK(IXR), 1) CALL SSCAL(NVAR, RAT, RWORK(IXR), 1) CALL SAXPY(NVAR, 1.0, RWORK(IXC), 1, RWORK(IXR), 1) INDEX = JXR + IT RWORK(INDEX) = XIT CALL CORECT(NVAR, RWORK(IXR), IT, RWORK(ITL), IERR, MODCON, * RWORK(IFP), NROW, NCOL, IPIVOT, ABSERR, RELERR, XSTEP, NEQN, * FNRM, FXNAME, FPNAME, SLNAME, LUNIT, DETNXR, IEXNXR) NTRCR = NTRCR + NCORXR ITO = IT XITO = XIT IF (IERR.EQ.0) GO TO 360 IF (IERR.EQ.(-1)) GO TO 400 IF (IERR.EQ.(-2)) GO TO 390 IF (IERR.EQ.(-3)) GO TO 410 C C*********************************************************************** C C 3. TANGENT AND LOCAL CONTINUATION PARAMETER CALCULATION. UNLESS C THE TANGENT AND LIMIT POINT CALCULATIONS WERE ALREADY PERFORMED C (BECAUSE THE LOOP WAS INTERRUPTED FOR A LIMIT POINT), SET UP AND C SOLVE THE EQUATION FOR THE TANGENT VECTOR. FORCE THE TANGENT C TO BE OF UNIT LENGTH, AND FORCE THE IPL-COMPONENT TO HAVE THE SAME C SIGN AS THE IPL-TH COMPONENT OF THE PREVIOUS TANGENT VECTOR, OR (ON C FIRST STEP) THE SAME SIGN AS THE INPUT DIRECTION DIRIPC. SET THE C LOCAL PARAMETER IPC TO THE LOCATION OF THE LARGEST COMPONENT C OF THE TANGENT VECTOR, UNLESS A LIMIT POINT IN THAT DIRECTION C APPEARS TO BE APPROACHING AND ANOTHER CHOICE IS AVAILABLE. C C*********************************************************************** C 40 IF (IRET.NE.2) GO TO 50 IRET = 0 GO TO 200 C C STORE OLD TANGENT IN TL, COMPUTE NEW TANGENT FOR XC C 50 IPL = IPC IF (KSTEP.GT.0) CALL SCOPY(NVAR, RWORK(ITC), 1, RWORK(ITL), 1) ICALL = 1 DETTXL = DETTXC IEXTXL = IEXTXC CALL TANGNT(NVAR, RWORK(IXF), IPC, RWORK(ITC), IRET, ICALL, * RWORK(IFP), IPIVOT, NEQN, DETTXC, IEXTXC, IPCFIX, ABSERR, * RELERR, NROW, NCOL, FPNAME, SLNAME) IF (IRET.EQ.(-2)) GO TO 460 IF (IRET.EQ.(-1)) GO TO 440 C C TANGNT RETURNED IPC, THE LOCATION OF THE LARGEST COMPONENT C OF THE TANGENT VECTOR. THIS WILL BE USED FOR THE LOCAL C PARAMETERIZATION OF THE CURVE UNLESS A LIMIT POINT IN IPC SEEMS C TO BE COMING. TO CHECK THIS, WE COMPARE TCIPC =TC(IPC) AND THE C SECOND LARGEST COMPONENT TCJPC =TC(JPC). IF TCJPC IS NO LESS C THAN .1 OF TCIPC, AND TC(JPC) IS LARGER THAN TL(JPC), C WHEREAS TC(IPC) IS LESS THAN TL(IPC), WE WILL RESET THE C LOCAL PARAMETER IPC =JPC. C C BUT IF IPCFIX.EQ.1, IGNORE THE CHECK C TLIPL = TCIPC INDEX = JTC + IPC TCIPC = RWORK(INDEX) JPC = IPC INDEX = JTL + IPC IF (IPCFIX.EQ.1) GO TO 60 IF (ABS(TCIPC).GE.ABS(RWORK(INDEX))) GO TO 60 IF (TLIPL.EQ.0.0) GO TO 60 INDEX = JTC + IPC RWORK(INDEX) = 0.0 JPC = ISAMAX(NVAR,RWORK(ITC),1) INDEX = JTC + JPC TCJPC = RWORK(INDEX) INDEX = JTC + IPC RWORK(INDEX) = TCIPC IF (ABS(TCJPC).LT.TENM1*ABS(TCIPC)) GO TO 60 INDEX = JTL + JPC IF (ABS(TCJPC).LT.ABS(RWORK(INDEX))) GO TO 60 IF (IWRITE.GE.2) WRITE (LUNIT,99981) IPC IPC = JPC 60 INDEX = JTC + IPL TCIPL = RWORK(INDEX) INDEX = JTL + IPC DTLIPC = DBLE(RWORK(INDEX)) DETTXC = DETTXC/TCIPL C C ADJUST SIGN OF TANGENT. C COMPARE THE SIGN OF TC(IPL) WITH SIGN OF TL(IPL), EXCEPT ON C FIRST STEP, COMPARE SIGN OF TC(IPL) WITH USER INPUT DIRIPC. C IF THESE SIGNS DIFFER, CHANGE THE SIGN OF TC TO FORCE AGREEMENT, C AND THE SIGN OF DETTXC. C THEN RECORD DIRIPC = SIGN OF DETERMINANT = SIGN(DETTXC). C STLIPL = DIRIPC IF (TLIPL.NE.0.0) STLIPL = SIGN(1.0,TLIPL) IF (SIGN(1.0,TCIPL).EQ.STLIPL) GO TO 70 CALL SSCAL(NVAR, -1.0, RWORK(ITC), 1) DETTXC = -DETTXC TCIPL = -TCIPL 70 INDEX = JTC + IPC TCIPC = RWORK(INDEX) INDEX = JTC + JPC TCJPC = RWORK(INDEX) DTCIPC = DBLE(TCIPC) DIRIPC = SIGN(1.0,DETTXC) IF (LIM.EQ.0) GO TO 80 TLLIM = TCLIM INDEX = JTC + LIM TCLIM = RWORK(INDEX) C C COMPUTE ALPHLC, THE ANGLE BETWEEN TANGENT AT XL AND TANGENT AT XC C AND HSECLC, THE EUCLIDEAN NORM OF SECANT FROM XL TO XC. C 80 IF (KSTEP.LE.0) GO TO 200 COSALF = 0.0D0 DO 90 I=1,NVAR INDEX = JTC + I JNDEX = JTL + I COSALF = COSALF + DBLE(RWORK(INDEX))*DBLE(RWORK(JNDEX)) INDEX = JXR + I JNDEX = JXF + I KNDEX = JXC + I RWORK(INDEX) = RWORK(JNDEX) - RWORK(KNDEX) 90 CONTINUE HSECLL = HSECLC HSECLC = SNRM2(NVAR,RWORK(IXR),1) TCOS = SNGL(COSALF) IF (TCOS.GT.1.0) TCOS = 1.0 IF (TCOS.LT.(-1.0)) TCOS = -1.0 TSIN = SQRT(1.0-TCOS*TCOS) ALPHLC = ATAN2(TSIN,TCOS) IF (IWRITE.GE.2) WRITE (LUNIT,99989) ALPHLC C C*********************************************************************** C C 4. LIMIT POINT CHECK. IF (LIM.NE.0) CHECK TO SEE IF C OLD AND NEW TANGENTS DIFFER IN SIGN OF LIM-TH COMPONENT. C IF SO, ATTEMPT TO COMPUTE A POINT XR BETWEEN XC AND XF C FOR WHICH TANGENT COMPONENT VANISHES C C*********************************************************************** C IF (LIM.LE.0 .OR. KSTEP.LE.0) GO TO 200 C C CHECK FOR LIMIT INTERVAL C IF (SIGN(1.0,TCLIM).EQ.SIGN(1.0,TLLIM)) GO TO 200 C C TEST FOR ACCEPTABLE ENDPOINTS C ATLLIM = ABS(TLLIM) IF (ATLLIM.GT.0.5*ABSERR) GO TO 110 C C IF XC IS LIMIT POINT, TL ALREADY CONTAINS TANGENT AT XC C 100 CALL SCOPY(NVAR, RWORK(IXC), 1, RWORK(IXR), 1) GO TO 350 110 ATCLIM = ABS(TCLIM) IF (ATCLIM.GT.0.5*ABSERR) GO TO 130 120 CALL SCOPY(NVAR, RWORK(IXF), 1, RWORK(IXR), 1) CALL SCOPY(NVAR, RWORK(ITC), 1, RWORK(ITL), 1) GO TO 350 C C TEST FOR SMALL INTERVAL C 130 INDEX = JXC + LIM XCLIM = RWORK(INDEX) INDEX = JXF + LIM XFLIM = RWORK(INDEX) DEL = ABS(XFLIM-XCLIM) XABS = AMAX1(ABS(XCLIM),ABS(XFLIM)) IF (DEL.GT.(ABSERR+RELERR*XABS)) GO TO 140 IF (ATLLIM.GT.ATCLIM) GO TO 120 GO TO 100 C C BEGIN ROOT-FINDING ITERATION WITH INTERVAL (0,1) AND C FUNCTION VALUES TL(LIM), TC(LIM). C 140 KOUNT = 0 A = 0.0 FA = TLLIM B = 1.0 FB = TCLIM C C SET LPC TO THE INDEX OF MAXIMUM ENTRY OF SECANT C (EXCEPT THAT LPC MUST NOT EQUAL LIM) C AND SAVE THE SIGN OF THE MAXIMUM COMPONENT IN DIRLPC C SO THAT NEW TANGENTS MAY BE PROPERLY SIGNED. C CALL SCOPY(NVAR, RWORK(IXF), 1, RWORK(IXR), 1) CALL SAXPY(NVAR, -1.0, RWORK(IXC), 1, RWORK(IXR), 1) INDEX = JXR + LIM RWORK(INDEX) = 0.0 LPC = ISAMAX(NVAR,RWORK(IXR),1) INDEX = JXR + LPC DIRLPC = SIGN(1.0,RWORK(INDEX)) C C SET FIRST APPROXIMATION TO ROOT TO WHICHEVER ENDPOINT C HAS SMALLEST LIM-TH COMPONENT OF TANGENT C IF (ATCLIM.LT.ATLLIM) GO TO 150 SN = 0.0 TSN = TLLIM CALL SCOPY(NVAR, RWORK(IXC), 1, RWORK(IXR), 1) GO TO 160 150 SN = 1.0 TSN = TCLIM CALL SCOPY(NVAR, RWORK(IXF), 1, RWORK(IXR), 1) C C CALL ROOTFINDER FOR APPROXIMATE ROOT SN. USE LINEAR COMBINATION C OF PREVIOUS ROOT SNL, AND ONE OF 0.0 AND 1.0 TO GET A STARTING C POINT FOR CORRECTION. RETURN TO CURVE ON LINE X(LPC)=CONSTANT, C COMPUTE TANGENT THERE, AND SET FUNCTION VALUE TO TANGENT(LIM) C 160 SNL = SN CALL ROOT(A, FA, B, FB, SN, TSN, KOUNT, IFLAG) NLMRT = NLMRT + 1 IF (IFLAG.LT.(-1)) GO TO 380 IF (IFLAG.EQ.(-1) .OR. IFLAG.EQ.0) GO TO 350 C C FIND WHETHER SN LIES IN (0.0,SNL) OR (SNL,1.0). C USE APPROPRIATE LINEAR COMBINATION TO GET STARTING POINT. C IF (SN.GT.SNL) GO TO 170 C C SET X(SN)=(SNL-SN)/(SNL-0.0)*X(0.0)+(SN-0.0)/(SNL-0.0)*X(SNL) C IF (SNL.LE.0.0) SCALER = 0.0 IF (SNL.GT.0.0) SCALER = SN/SNL CALL SSCAL(NVAR, SCALER, RWORK(IXR), 1) SCALER = 1.0 - SCALER CALL SAXPY(NVAR, SCALER, RWORK(IXC), 1, RWORK(IXR), 1) GO TO 180 C C SET X(SN)=(SN-SNL)/(1.0-SNL)*X(1.0)+(1.0-SN)/(1.0-SNL)*X(SNL) C 170 IF (SNL.GE.1.0) SCALER = 0.0 IF (SNL.LT.1.0) SCALER = (1.0-SN)/(1.0-SNL) CALL SSCAL(NVAR, SCALER, RWORK(IXR), 1) SCALER = 1.0 - SCALER CALL SAXPY(NVAR, SCALER, RWORK(IXF), 1, RWORK(IXR), 1) C C NOTE THAT CORRECTION IS ALWAYS DONE WITH MODIFIED NEWTON PROCESS. C MODLIM=1 UNLESS CORRECTOR RETURNS A NON FATAL ERROR, WHEN WE SET C MODLIM=0 AND TRY AGAIN. C 180 MODLIM = 1 190 CALL CORECT(NVAR, RWORK(IXR), LPC, RWORK(ITL), IERR, MODLIM, * RWORK(IFP), NROW, NCOL, IPIVOT, ABSERR, RELERR, XSTPLM, NEQN, * FNRM, FXNAME, FPNAME, SLNAME, LUNIT, DETNXR, IEXNXR) NLMCR = NLMCR + NCORXR IF (IERR.NE.0 .AND. MODLIM.EQ.0) GO TO 380 IF (IERR.NE.0) MODLIM = 0 IF (IERR.NE.0) GO TO 190 ICALL = 1 LPCFIX = 1 CALL TANGNT(NVAR, RWORK(IXR), LPC, RWORK(ITL), IRET, ICALL, * RWORK(IFP), IPIVOT, NEQN, DETTXR, IEXTXR, LPCFIX, ABSERR, * RELERR, NROW, NCOL, FPNAME, SLNAME) IF (IRET.LT.0) GO TO 380 C C ADJUST THE SIGN OF THE TANGENT SO THAT THE LPC-TH COMPONENT C HAS THE SAME SIGN AS THE LPC-TH COMPONENT OF THE SECANT C INDEX = JTL + LPC IF (DIRLPC.NE.SIGN(1.0,RWORK(INDEX))) CALL SSCAL(NVAR, -1.0, * RWORK(ITL), 1) C C SEE IF WE CAN ACCEPT THE NEW POINT BECAUSE TANGENT(LIM) IS SMALL C OR IF WE MUST GO ON. C INDEX = JTL + LIM TSN = RWORK(INDEX) IF (ABS(TSN).LE.0.5*ABSERR) GO TO 350 GO TO 160 C C*********************************************************************** C C 5. STEP LENGTH COMPUTATION. COMPUTE STEPLENGTH HTANCF C C THE FORMULAS UNDERLYING THE ALGORITHM ARE C C LET C C ALPHLC = MAXIMUM OF ARCCOS(TL,TC) AND ALFMIN = 2*ARCCOS(1-EPMACH) C HSECLC = NORM(XL-XC) C HSECLL = NORM(XL-XLL) C ABSIN = ABS(SIN(.5*ALPHLC)) C CURVLC = LAST VALUE OF CURVCF C CURVCF = 2*ABSIN/HSECLC C CORDXF = OPTIMIZED CORRECTOR DISTANCE TO CURRENT CONTINUATION C POINT. C BUT CORDXF FORCED TO LIE BETWEEN .01*HSECLC AND HSECLC. C UNLESS (CORDXF.EQ.0.0), BECAUSE THE PREDICTED POINT WAS C ACCEPTED. IN SUCH A CASE, SET HTANCF=HFACT*HSECLC C INSTEAD OF USING FIRST ESTIMATE FOR HTANCF. C C THEN C C CURVXF = CURVCF + HSECLC*(CURVCF-CURVLC)/(HSECLC+HSECLL) C A SIMPLER FORMULA IS USED IF WE DO NOT HAVE DATA AT TWO PREVIOUS C POINTS. C C IF (HSECLC.GT.0.0) CURVXF MUST BE AT LEAST AS LARGE AS .01/HSECLC C IF (HSECLC.EQ.0.0) CURVXF MUST BE AT LEAST .001 C C FIRST ESTIMATE (UNLESS (CORDXF.EQ.0.0) ) C C HTANCF = SQRT(2*CORDXF/CURVXF) C C ADJUSTED VALUE C C HTANCF = HTANCF*(1.0 + HTANCF*(TC(IPC)-TL(IPC))/(2*HSECLC*TC(IPC))) C C READJUSTMENT AND TRUNCATIONS C C IF STEPSIZE REDUCTION OCCURRED DURING LAST CORRECTOR PROCESS, C HTANCF IS FORCED TO BE LESS THAN (HFACT-1)*HSECLC/2. C C HTANCF MUST LIE BETWEEN (HSECLC/HFACT) AND (HSECLC*HFACT). C C HTANCF IS ALWAYS FORCED TO LIE BETWEEN HMIN AND HMAX. C C*********************************************************************** C C CHECK IF DEFAULT STEP MUST BE USED C IF PREVIOUS STEP WAS OF SIZE ZERO, USE STEPSIZE HDEF=(HMIN+HMAX)/2 C 200 IF (KSTEP.GT.0 .AND. HSECLC.GT.0.0) GO TO 210 IF (KSTEP.GT.0) HTANCF = HDEF GO TO 230 210 IF (ALPHLC.LT.ALFMIN) ALPHLC = ALFMIN ABSIN = ABS(SIN(.5*ALPHLC)) C C COMPUTE NEW CURVATURE DATA C CURVLC = CURVCF CURVCF = 2.0*ABSIN/HSECLC CURVXF = CURVCF IF (HSECLL.NE.0.0) CURVXF = CURVCF + HSECLC*(CURVCF-CURVLC)/ * (HSECLC+HSECLL) TEMP = TENM3 IF (HSECLC.GT.0.0) TEMP = TENM2/HSECLC CURVXF = AMAX1(CURVXF,TEMP) IF (IWRITE.GE.2) WRITE (LUNIT,99988) CURVCF IF (IWRITE.GE.2) WRITE (LUNIT,99987) CURVXF C C IF CONVERGENCE DISTANCE IS ZERO, SET ESTIMATE TO HFACT*HSECLC. C OTHERWISE, TRUNCATE CORDXF TO LIE BETWEEN .01*HSECLC AND HSECLC. C HTANCF = HFACT*HSECLC IF (CORDXF.EQ.0.0) GO TO 220 TEMP = TENM2*HSECLC CORDXF = AMAX1(CORDXF,TEMP) CORDXF = AMIN1(CORDXF,HSECLC) C C SET HTANCF, THEN ADJUST FOR CURVATURE IN CONTINUATION C PARAMETER DIRECTION, THEN TRUNCATE C HTANCF = SQRT(2.0*CORDXF/CURVXF) 220 IF (NREDXF.GT.0) HTANCF = AMIN1(HTANCF,(HFACT-1.0)*HSECLC*.5) DADJUS = 1.0D0 + (1.0D0-DTLIPC/DTCIPC)*DBLE(.5*HTANCF)/ * DBLE(HSECLC) HTANCF = HTANCF*SNGL(DADJUS) HTANCF = AMAX1(HTANCF,HSECLC*HRED) HTANCF = AMIN1(HTANCF,HSECLC*HFACT) HTANCF = AMAX1(HTANCF,HMIN) HTANCF = AMIN1(HTANCF,HMAX) C C*********************************************************************** C C 6. PREDICTION AND CORRECTION STEPS. USING XR=XC+HTANCF*TC C AS STARTING POINT, CORRECT XR WITH A FULL OR MODIFIED NEWTON C ITERATION. IF CORECT FAILS, REDUCE STEPSIZE USED FOR PREDICTOR C POINT, AND TRY AGAIN. CORRECTION WILL BE ABANDONED IF STEPSIZE C FALLS BELOW HMIN. C C*********************************************************************** C 230 KSTEPL = KSTEP KSTEP = KSTEP + 1 NREDXF = 0 240 CALL SCOPY(NVAR, RWORK(IXF), 1, RWORK(IXR), 1) CALL SAXPY(NVAR, HTANCF, RWORK(ITC), 1, RWORK(IXR), 1) 250 IF (IWRITE.GE.2) WRITE (LUNIT,99986) HTANCF INDEX = JXR + 1 JNDEX = JXR + NVAR IF (IWRITE.GE.3) WRITE (LUNIT,99985) (RWORK(I),I=INDEX,JNDEX) CALL CORECT(NVAR, RWORK(IXR), IPC, RWORK(ITL), IERR, MODCON, * RWORK(IFP), NROW, NCOL, IPIVOT, ABSERR, RELERR, XSTEP, NEQN, * FNRM, FXNAME, FPNAME, SLNAME, LUNIT, DETNXR, IEXNXR) NCNCR = NCNCR + NCORXR IF (IERR.EQ.0) GO TO 270 IF (IERR.EQ.(-1)) GO TO 450 C C NO CONVERGENCE, TRY A SMALLER STEPSIZE C HTANCF = HRED*HTANCF IF (HTANCF.LT.HMIN) GO TO 420 NREDXF = NREDXF + 1 IF (IERR.EQ.(-2)) GO TO 260 GO TO 240 260 CALL SAXPY(NVAR, -1.0, RWORK(IXF), 1, RWORK(IXR), 1) CALL SSCAL(NVAR, HRED, RWORK(IXR), 1) CALL SAXPY(NVAR, 1.0, RWORK(IXF), 1, RWORK(IXR), 1) GO TO 250 C C*********************************************************************** C C 7. SUCCESSFUL STEP. STORE INFORMATION BEFORE RETURN. C UPDATE OLD AND CURRENT CONTINUATION POINTS. C COMPUTE CORDXF, THE SIZE OF THE CORRECTOR STEP. COMPUTE C A FACTOR THETA WHICH RESCALES CORDXF TO A VALUE WHICH WOULD C CORRESPOND TO A DESIRABLE NUMBER OF CORRECTOR STEPS C (4 FOR FULL NEWTON, 10 FOR MODIFIED NEWTON). C SEE REFERENCE DEN HEIJER AND RHEINBOLDT, LOC. CIT. C C*********************************************************************** C 270 NRDSUM = NRDSUM + NREDXF IF (NREDXF.NE.0 .AND. IWRITE.GE.2) WRITE (LUNIT,99984) NREDXF C C COMPUTE CORRECTOR STEP = XC+HTANCF*TC-XF C SET CORDXF = MAX NORM OF CORRECTOR STEP C CALL SCOPY(NVAR, RWORK(IXF), 1, RWORK(IXC), 1) CALL SAXPY(NVAR, -1.0, RWORK(IXR), 1, RWORK(IXC), 1) CALL SAXPY(NVAR, HTANCF, RWORK(ITC), 1, RWORK(IXC), 1) IMAX = ISAMAX(NVAR,RWORK(IXC),1) INDEX = JXC + IMAX CORDXF = ABS(RWORK(INDEX)) IF (NCORXR.EQ.0) CORDXF = 0.0 C C MODIFY CORDXF TO A VALUE THAT WOULD CORRESPOND TO THE C DESIRED NUMBER OF CORRECTOR STEPS C IF (CORDXF.EQ.0.0) GO TO 340 OMEGA = XSTEP/CORDXF THETA = 0.0 IF (MODCON.EQ.1) GO TO 300 C C FULL NEWTON METHOD FOR CORRECTOR STEPS C IF (NCORXR.LE.1) THETA = 8.0 IF (NCORXR.EQ.4) THETA = 1.0 IF (THETA.NE.0.0) GO TO 330 IF (NCORXR.GT.4) GO TO 280 LK = 4*NCORXR - 7 THETA = 1.0 IF (OMEGA.GE.WRGE(LK)) GO TO 330 LST = LK INDEX = LK + 1 IF (OMEGA.GE.WRGE(INDEX)) GO TO 290 LST = LK + 2 IF (OMEGA.GE.WRGE(LST)) GO TO 290 THETA = 8.0 GO TO 330 280 THETA = 0.125 IF (NCORXR.GE.7) GO TO 330 LK = 4*NCORXR - 16 IF (OMEGA.LE.WRGE(LK)) GO TO 330 LST = 2*NCORXR - 1 290 INDEX = LST + 1 THETA = ACOF(LST) + ACOF(INDEX)*ALOG(OMEGA) GO TO 330 C C MODIFIED NEWTON METHOD FOR CORRECTOR STEPS C 300 IF (NCORXR.LE.1) THETA = 8.0 IF (NCORXR.EQ.10) THETA = 1.0 IF (THETA.NE.0.0) GO TO 330 EXPON = FLOAT(NCORXR-1)/FLOAT(NCORXR-10) C C AVOID OVERFLOW OR UNDERFLOW BY ANTICIPATING C CUTOFF VALUES OF THETA C TEST = 8.0**EXPON IF (NCORXR.GT.10) GO TO 310 IF (TEST.GT.OMEGA) THETA = 8.0 IF (1.0.LT.(TEST*OMEGA)) THETA = .125 IF (THETA.NE.0.0) GO TO 330 GO TO 320 310 IF (TEST.LT.OMEGA) THETA = 8.0 IF (1.0.GT.(TEST*OMEGA)) THETA = .125 IF (THETA.NE.0.0) GO TO 330 320 EXPON = 1.0/EXPON THETA = OMEGA**EXPON THETA = AMAX1(THETA,0.125) THETA = AMIN1(THETA,8.0) C C SET THE MODIFIED VALUE OF CORDXF C 330 CORDXF = THETA*CORDXF IF (IWRITE.GE.3) WRITE (LUNIT,99983) OMEGA, THETA IF (IWRITE.GE.2) WRITE (LUNIT,99982) CORDXF 340 CALL SCOPY(NVAR, RWORK(IXF), 1, RWORK(IXC), 1) CALL SCOPY(NVAR, RWORK(IXR), 1, RWORK(IXF), 1) GO TO 370 C C*********************************************************************** C C RETURNS. SET VALUE OF IRET. IF AN ERROR OCCURRED, PRINT C A MESSAGE AS WELL. C C*********************************************************************** C C C RETURN LIMIT POINT C 350 IRET = 2 RETURN C C RETURN WITH TARGET POINT C 360 IRET = 1 RETURN C C RETURN WITH CONTINUATION POINT C CALL SCOPY(NVAR, RWORK(IXF), 1, RWORK(IXR), 1) 370 IRET = 0 RETURN C C ERROR RETURNS C 380 IRET = -1 IF (IWRITE.GE.1) WRITE (LUNIT,99999) RETURN 390 IRET = -2 IF (IWRITE.GE.1) WRITE (LUNIT,99998) RETURN 400 IRET = -3 IF (IWRITE.GE.1) WRITE (LUNIT,99997) RETURN 410 IRET = -4 IF (IWRITE.GE.1) WRITE (LUNIT,99996) RETURN 420 IRET = -5 IF (IWRITE.GE.1) WRITE (LUNIT,99995) HTANCF, HMIN RETURN 430 IRET = -6 IF (IWRITE.GE.1) WRITE (LUNIT,99994) RETURN 440 IRET = -7 IF (IWRITE.GE.1) WRITE (LUNIT,99993) RETURN 450 IRET = -8 IF (IWRITE.GE.1) WRITE (LUNIT,99992) RETURN 460 IRET = -9 IF (IWRITE.GE.1) WRITE (LUNIT,99991) RETURN 470 IRET = -10 IF (IWRITE.GE.1) WRITE (LUNIT,99990) NVAR, ISIZE RETURN C C*********************************************************************** C 99999 FORMAT (26H0LIMIT POINT FINDER FAILED) 99998 FORMAT (53H0CORRECTOR, SEEKING TARGET POINT, TOOK TOO MANY STEPS) 99997 FORMAT (54H0CORRECTOR, SEEKING TARGET, CALLED SOLVER WHICH FAILED) 99996 FORMAT (48H0CORRECTOR, SEEKING TARGET, FAILED WITH BAD STEP) 99995 FORMAT (10H0STEPSIZE , G14.7, 15H LESS THAN HMIN, G14.7) 99994 FORMAT (42H0NORM OF F(X) IS TOO LARGE ON INITIAL CALL) 99993 FORMAT (34H0SOLVER FAILED IN CALL FROM TANGNT) 99992 FORMAT (34H0SOLVER FAILED IN CALL FROM CORECT) 99991 FORMAT (23H0TANGENT VECTOR IS ZERO) 99990 FORMAT (26H0UNACCEPTABLE INPUT NVAR=, I10, 7H ISIZE=, I10) 99989 FORMAT (36H ANGLE BETWEEN OLD AND NEW TANGENTS=, G14.7, 7H RADIAN, * 1HS) 99988 FORMAT (21H OLD CURVATURE , G14.7) 99987 FORMAT (21H PREDICTED CURVATURE , G14.7) 99986 FORMAT (29H PREDICTOR IS USING STEPSIZE=, G14.7) 99985 FORMAT (12H PREDICTED X/1X, 5G14.7) 99984 FORMAT (1H , I2, 20H STEPSIZE REDUCTIONS) 99983 FORMAT (7H OMEGA=, G14.7, 7H THETA=, G14.7) 99982 FORMAT (43H ROUGH RADIUS OF CONVERGENCE FOR CORRECTOR , G14.7) 99981 FORMAT (43H TANGNT ANTICIPATES LIMIT POINT IN VARIABLE, I3) END SUBROUTINE CORECT(NVAR, X, IHOLD, WORK, IERR, MODNEW, FPRYM, COR 10 * NROW, NCOL, IPIVOT, ABSERR, RELERR, XSTEP, NEQN, FNRM, FXNAME, * FPNAME, SLNAME, LUNIT, DET, IEX) C C*********************************************************************** C C SUBROUTINE CORECT PERFORMS THE CORRECTION OF AN APPROXIMATE C SOLUTION OF THE EQUATION F(X)=0.0. THE CORRECTION METHOD IS EITHER C FULL OR MODIFIED NEWTON'S METHOD. FOR MODIFIED NEWTON'S METHOD, C THE JACOBIAN IS TO BE EVALUATED ONLY AT THE STARTING POINT. C IF B IS THE VALUE OF X(IHOLD) FOR THE INPUT STARTING POINT, C THEN THE AUGMENTING EQUATION IS X(IHOLD)=B, THAT IS, THE IHOLD-TH C COMPONENT OF X IS TO BE HELD FIXED. THE AUGMENTED SYSTEM TO BE C SOLVED IS THEN DFA(X,IHOLD)*(-DELX)=FA(X) C C ACCEPTANCE AND REJECTION CRITERIA- C C LET NCORXR BE THE CURRENT ITERATION STEP. NCORXR=0 FOR INPUT C PREDICTED POINT. C AFTER EACH ITERATION, THE CURRENT POINT IS ACCEPTED C IF ANY OF THE FOLLOWING HOLD- C C STRONG ACCEPTANCE CRITERION C C 1. ABS(F(X)).LE.ABSERR AND XSTEP.LE.(ABSERR+RELERR*XNRM)) C C WEAK ACCEPTANCE CRITERIA C C 2. (NCORXR.EQ.0) AND C ABS(F(X)).LE.(.5*ABSERR) C 3. ABS(F(X)).LE.EPSATE C 4. (NCORXR.GE.2) AND C (ABS(F(X))+ABS(FOLD(X)).LE.ABSERR AND C XSTEP.LE.8*(ABSERR+RELERR*XNRM) C 5. (NCORXR.GE.2) AND C ABS(F(X).LE.8.0*ABSERR AND C (XSTEP+XSTEPOLD).LE.(ABSERR+RELERR*XNRM) C C AFTER EACH ITERATION, THE ITERATION IS TO BE ABORTED IF C ANY OF THE FOLLOWING HOLD C C 1. FNRM.GT.(FMP*FNRML+ABSERR) C 2. (NCORXR.GE.2) AND C (XSTEP.GT.(FMP*XSTEPL+ABSERR)) C 3. NCORXR.GE.MAXCOR (NUMBER OF ITERATIONS EXCEEDED). C C (FMP=2.0 FOR NCORXR.EQ.1, FMP=1.05 FOR NCORXR.GT.1) C C C INPUT C NVAR = NUMBER OF VARIABLES. C X = THE STARTING POINT FOR THE CORRECTOR ITERATION. C IHOLD = COMPONENT OF X THAT WILL NOT BE CHANGED DURING ITERATION. C WORK = WORK SPACE OF DIMENSION NVAR. C MODNEW= FLAG FOR TYPE OF NEWTON'S METHOD TO BE USED. C IF (MODNEW.EQ.0), JACOBIAN IS TO BE EVALUATED AT EVERY C CORRECTOR ITERATE. MAXCOR IS SET TO 10. C IF (MODNEW.EQ.1), THE JACOBIAN IS ONLY EVALUATED AT THE C STARTING POINT, AND MAXCOR IS SET TO 20. C FPRYM = STORAGE SPACE FOR AUGMENTED JACOBIAN OF SIZE (NROW,NCOL). C NROW = NUMBER OF ROWS USED FOR STORING FPRYM. C NCOL = NUMBER OF COLUMNS USED FOR STORING FPRYM. C IPIVOT= STORAGE SPACE FOR PIVOT INDICES OF SIZE NVAR. C ABSERR= ABSOLUTE ERROR TOLERANCE. C RELERR= RELATIVE ERROR TOLERANCE. C NEQN = NUMBER OF EQUATIONS = NVAR-1. C FXNAME= EXTERNAL NAME OF FUNCTION EVALUATOR. C FPNAME= EXTERNAL NAME OF JACOBIAN EVALUATOR. C SLNAME= EXTERNAL NAME OF SOLVER. C LUNIT = LOGICAL FORTRAN OUTPUT UNIT. C C OUTPUT C X = SOLUTION VECTOR ON A SUCCESSFUL CALL TO CORECT. C WORK = THE RESIDUAL F(X), AFTER A SUCCESSFUL CALL TO CORECT. C IERR = THE RETURN FLAG WITH THE FOLLOWING VALUES C 0 SUCCESSFUL CORRECTION. VECTOR X RETURNED USUALLY C SATISIFIES ABS(F(X)).LE.ABSERR. C -1 ERROR RETURN FROM SOLVER CALLED BY CORECT. C -2 MAXIMUM NUMBER OF CORRECTOR ITERATIONS WERE TAKEN. C -3 CORRECTOR STEP WAS UNACCEPTABLE, CORRECTION FAILED. C XSTEP = ABSOLUTE VALUE OF MAXIMUM COMPONENT OF LAST CORRECTION C TO X. C FNRM = ABSOLUTE VALUE OF MAXIMUM COMPONENT OF FUNCTION VALUE C OF CORRECTED X. C DET = MANTISSA OF DETERMINANT OF NEWTON SYSTEM ON LAST C EVALUATION. C IEX = BINARY EXPONENT OF DETERMINANT OF NEWTON SYSTEM ON LAST C EVALUATION. C C OUTPUT THROUGH COMMON /COUNT2/ C NCORXR= THE NUMBER OF CORRECTOR ITERATIONS TAKEN ON THIS CALL. C C THIS SUBROUTINE IS CALLED BY C PITCON C AND CALLS C SOLVER (SLNAME) C FUNCTION (FXNAME) C FORTRAN ABS C LINPAK ISAMAX C LINPAK SAXPY C C*********************************************************************** C EXTERNAL FXNAME, FPNAME, SLNAME REAL X(NVAR), WORK(NVAR), FPRYM(NROW,NCOL) INTEGER IPIVOT(NVAR) COMMON /COUNT1/ ICRSL, ITNSL, NSTCR, NCNCR, NTRCR, NLMCR, NLMRT COMMON /COUNT2/ IFEVAL, IPEVAL, ISOLVE, NREDXF, NRDSUM, NCORXR, * NCRSUM COMMON /OUTPUT/ IWRITE COMMON /TOL/ EPMACH, EPSATE, EPSQRT C C INITIALIZE C NCORXR = 0 MAXCOR = 10 IF (MODNEW.EQ.1) MAXCOR = 20 IERR = 0 FMP = 2.0 ICALL = 1 XSTEP = 0.0 C C GET INITIAL FUNCTION VALUE C CALL FXNAME(NVAR, X, WORK) IFEVAL = IFEVAL + 1 IFMAX = ISAMAX(NEQN,WORK,1) FNRM = ABS(WORK(IFMAX)) IF (IWRITE.GE.2) WRITE (LUNIT,99992) FNRM WORK(NVAR) = 0.0 C C CHECK WEAK ACCEPTANCE OF PREDICTED POINT C IF (FNRM.LE.0.5*ABSERR) GO TO 70 C C NEWTON ITERATION LOOP BEGINS C DO 30 I=1,MAXCOR NCORXR = I C C SOLVE SYSTEM FPRYM*(-DELX)=FX C CALL SLNAME(NVAR, X, WORK, IHOLD, DET, IEX, IERR, ICALL, * MODNEW, FPRYM, NROW, NCOL, IPIVOT, FPNAME) ICRSL = ICRSL + 1 IF (MODNEW.EQ.1) ICALL = 0 ISOLVE = ISOLVE + 1 IF (IERR.NE.0) GO TO 60 C C STORE OLD INFORMATION, THEN GET NORMS OF C FUNCTION VALUE, SOLUTION POINT X, AND STEP FROM C PREVIOUS X TO CURRENT ONE. C FNRML = FNRM XSTEPL = XSTEP CALL SAXPY(NVAR, -1.0, WORK, 1, X, 1) ISMAX = ISAMAX(NVAR,WORK,1) XSTEP = ABS(WORK(ISMAX)) IF (IWRITE.GE.2) WRITE (LUNIT,99993) XSTEP IXMAX = ISAMAX(NVAR,X,1) XNRM = ABS(X(IXMAX)) CALL FXNAME(NVAR, X, WORK) IFEVAL = IFEVAL + 1 IFMAX = ISAMAX(NEQN,WORK,1) FNRM = ABS(WORK(IFMAX)) IF (IWRITE.GE.2) WRITE (LUNIT,99992) FNRM WORK(NVAR) = 0.0 C C CHECK FOR STRONG ACCEPTANCE CRITERION C IF (FNRM.LE.ABSERR .AND. XSTEP.LE.(ABSERR+RELERR*XNRM)) GO TO 80 C C CHECK FOR WEAK ACCEPTANCE CRITERIA C IF (FNRM.LE.EPSATE) GO TO 70 IF (NCORXR.LE.1) GO TO 20 IF ((FNRM+FNRML).LE.ABSERR .AND. XSTEP.LE.8.0*(ABSERR+RELERR* * XNRM)) GO TO 70 IF (FNRM.LE.8.0*ABSERR .AND. (XSTEP+XSTEPL).LE.(ABSERR+RELERR* * XNRM)) GO TO 70 GO TO 10 C C SEE IF NEWTON ITERATION SHOULD BE ABORTED C 10 IF (XSTEP.GT.(FMP*XSTEPL+ABSERR)) GO TO 40 20 IF (FNRM.GT.(FMP*FNRML+ABSERR)) GO TO 40 FMP = 1.05 30 CONTINUE C C FALL THROUGH LOOP IF MAXIMUM NUMBER OF STEPS TAKEN C WITHOUT ACCEPTANCE OR REJECTION C GO TO 50 C C UNSUCCESSFUL STEP C 40 IERR = -3 IF (IWRITE.GE.2) WRITE (LUNIT,99995) GO TO 90 C C MAXIMUM NUMBER OF CORRECTOR STEPS REACHED C 50 IERR = -2 IF (IWRITE.GE.2) WRITE (LUNIT,99996) GO TO 90 C C ERROR RETURN FROM SOLVING ROUTINE C 60 IERR = -1 IF (IWRITE.GE.2) WRITE (LUNIT,99997) GO TO 90 C C ACCEPTED POINT SATISFIED A WEAK CRITERION C 70 IF (IWRITE.GE.1) WRITE (LUNIT,99994) GO TO 80 C C ACCEPTED POINT SATISFIED THE STRONG CRITERION C 80 IERR = 0 90 NCRSUM = NCRSUM + NCORXR IF (IWRITE.GE.2) WRITE (LUNIT,99999) NCORXR IF (IWRITE.GE.2) WRITE (LUNIT,99998) IHOLD RETURN C C*********************************************************************** C 99999 FORMAT (16H CORRECTOR TOOK , I2, 6H STEPS) 99998 FORMAT (24H CORRECTOR HELD VARIABLE, I3, 6H FIXED) 99997 FORMAT (35H0SOLVER FAILED, CALLED BY CORRECTOR) 99996 FORMAT (25H0TOO MANY CORRECTOR STEPS) 99995 FORMAT (24H0CORRECTOR STEP REJECTED) 99994 FORMAT (40H ACCEPTED POINT SATISFIED WEAK CRITERION) 99993 FORMAT (23H CORRECTOR STEP LENGTH=, G14.7) 99992 FORMAT (15H FUNCTION NORM=, G14.7) END SUBROUTINE TANGNT(NVAR, X, IP, TAN, IRET, ICALL, FPRYM, IPIVOT, TAN 10 * NEQN, DET, IEX, IPCFIX, ABSERR, RELERR, NROW, NCOL, FPNAME, * SLNAME) C C*********************************************************************** C C SUBROUTINE TANGNT COMPUTES A UNIT TANGENT VECTOR TO THE SOLUTION C CURVE OF THE UNDERDETERMINED NONLINEAR SYSTEM FX = 0. THE C TANGENT VECTOR TAN IS THE SOLUTION OF THE LINEAR SYSTEM C C DFA(X,IPL)*TAN = E(NVAR) C C WHERE DFA(X,IPL) IS THE AUGMENTED JACOBIAN WHOSE FIRST NEQN ROWS C ARE DFX/DX (X), THE DERIVATIVE OF FX EVALUATED AT X, AND WHOSE LAST C ROW IS (E(IPL)) TRANSPOSE, THE NVAR COMPONENT EUCLIDEAN COORDINATE C VECTOR WITH 1 IN THE IPL-TH ENTRY AND ZEROS ELSEWHERE. E(NVAR) IS C THE NVAR COMPONENT EUCLIDEAN COORDINATE VECTOR WITH ONE IN THE C LAST COMPONENT. C C THE TANGENT VECTOR IS THEN NORMALIZED. ADJUSTMENT OF SIGN C IS DONE OUTSIDE THIS ROUTINE. C C INPUT C C NVAR = THE NUMBER OF VARIABLES. C X = POINT AT WHICH THE TANGENT SYSTEM IS TO BE SOLVED. C IP = INDEX USED FOR LAST EQUATION IN TANGENT SYSTEM. C TAN = WORK SPACE OF SIZE NVAR. C ICALL = JACOBIAN EVALUATION SWITCH, SET TO 1 TO FORCE SOLVER C TO EVALUATE AUGMENTED JACOBIAN AT X. C FPRYM = WORK SPACE TO HOLD AUGMENTED JACOBIAN, OF SIZE (NROW,NCOL). C IPIVOT= WORK SPACE TO HOLD PIVOT INDICES, OF SIZE NVAR. C NEQN = NUMBER OF EQUATIONS, NVAR-1. C IPCFIX= SWITCH WHICH CAN OVERRIDE CHOICE OF NEW CONTINUATION C VARIABLE, FORCING OLD ONE TO BE USED INSTEAD. C IPCFIX.NE.1 MEANS CHOOSE NEW CONTINUATION VARIABLE. C IPCFIX.EQ.1 MEANS MUST USE OLD CONTINUATION VARIABLE. C ABSERR= ABSOLUTE ERROR TOLERANCE. C RELERR= RELATIVE ERROR TOLERANCE. C NROW = NUMBER OF ROWS USED IN STORING FPRYM. C NCOL = NUMBER OF COLUMNS USED IN STORING FPRYM. C FPNAME= EXTERNAL NAME OF JACOBIAN EVALUATION ROUTINE. C SLNAME= EXTERNAL NAME OF SOLVING ROUTINE. C C OUTPUT C C IP = UNCHANGED IF IPCFIX.EQ.1. OTHERWISE, IP IS THE C LOCATION OF LARGEST COMPONENT OF TANGENT VECTOR TAN, C THE CANDIDATE FOR NEW CONTINUATION COMPONENT. C TAN = A UNIT TANGENT VECTOR AT X WHOSE SIGN MUST STILL BE SET. C IRET = ERROR FLAG. C 0 NO ERROR. C -1 ERROR RETURN FROM SOLVER. C -2 NORM OF TANGENT VECTOR IS ZERO. C DET = BINARY MANTISSA OF DETERMINANT OF JACOBIAN DFA(X,IPL) C IEX = BINARY EXPONENT OF THE DETERMINANT OF JACOBIAN DFA(X,IPL) C C THIS SUBROUTINE IS CALLED BY C PITCON C AND CALLS C SOLVER (SLNAME) C LINPAK ISAMAX C LINPAK SNRM2 C LINPAK SSCAL C C*********************************************************************** C EXTERNAL FPNAME, SLNAME REAL X(NVAR), TAN(NVAR), FPRYM(NROW,NCOL) INTEGER IPIVOT(NVAR) COMMON /COUNT1/ ICRSL, ITNSL, NSTCR, NCNCR, NTRCR, NLMCR, NLMRT COMMON /COUNT2/ IFEVAL, IPEVAL, ISOLVE, NREDXF, NRDSUM, NCORXR, * NCRSUM COMMON /OUTPUT/ IWRITE COMMON /TOL/ EPMACH, EPSATE, EPSQRT C C COMPUTE TANGENT VECTOR C DO 10 I=1,NEQN TAN(I) = 0.0 10 CONTINUE TAN(NVAR) = 1.0 IERR = 0 MODTAN = 0 CALL SLNAME(NVAR, X, TAN, IP, DET, IEX, IERR, ICALL, MODTAN, * FPRYM, NROW, NCOL, IPIVOT, FPNAME) ITNSL = ITNSL + 1 ISOLVE = ISOLVE + 1 IF (IERR.NE.0) IRET = -1 IF (IRET.LT.0) RETURN C C UNLESS IPCFIX.EQ.1, SET IP TO MAXIMUM ENTRY OF TAN C IF (IPCFIX.NE.1) IP = ISAMAX(NVAR,TAN,1) C C OBTAIN EUCLIDEAN NORM OF TANGENT VECTOR C TNORM = SNRM2(NVAR,TAN,1) IF (TNORM.EQ.0.0) IRET = -2 IF (IRET.LT.0) RETURN C C NORMALIZE THE VECTOR C SCALER = 1.0/TNORM CALL SSCAL(NVAR, SCALER, TAN, 1) RETURN C C*********************************************************************** C END SUBROUTINE ROOT(A, FA, B, FB, U, FU, KOUNT, IFLAG) ROO 10 C C*********************************************************************** C C SUBROUTINE ROOT SEEKS A ROOT OF THE EQUATION F(X)=0.0, C GIVEN A STARTING INTERVAL (A,B) ON WHICH F CHANGES SIGN. C ON FIRST CALL TO ROOT, THE INTERVAL AND FUNCTION VALUES FA AND FB C ARE INPUT AND AN APPROXIMATION U FOR THE ROOT IS RETURNED. C BEFORE EACH SUBSEQUENT CALL, THE USER EVALUATES FU=F(U), AND THE C PROGRAM TRIES TO RETURN A BETTER APPROXIMATION U. C C THIS PROGRAM IS BASED ON THE FORTRAN FUNCTION ZERO C GIVEN IN THE BOOK C ALGORITHMS FOR MINIMIZATION WITHOUT DERIVATIVES C BY RICHARD P. BRENT, PRENTICE HALL,INC, 1973 C C THE MODIFICATIONS WERE DONE BY JOHN BURKARDT. C C ON INPUT C C A - IS ONE ENDPOINT OF AN INTERVAL IN WHICH F CHANGES SIGN. C FA - THE VALUE OF F(A). THE USER MUST EVALUATE F(A) BEFORE C FIRST CALL ONLY. THEREAFTER THE PROGRAM SETS FA. C B - IS THE OTHER ENDPOINT OF THE INTERVAL IN WHICH C F CHANGES SIGN. NOTE THAT THE PROGRAM WILL RETURN C IMMEDIATELY WITH AN ERROR FLAG IF FB*FA.GT.0.0. C FB - THE VALUE OF F(B). THE USER MUST EVALUATE F(B) BEFORE C FIRST CALL ONLY. THERAFTER THE PROGRAM SETS FB. C U - ON FIRST CALL, U SHOULD NOT BE SET BY THE USER. C ON SUBSEQUENT CALLS, U SHOULD NOT BE CHANGED C FROM ITS OUTPUT VALUE, THE CURRENT APPROXIMANT C TO THE ROOT. C FU - ON FIRST CALL, FU SHOULD NOT BE SET BY THE USER. C THEREAFTER, THE USER SHOULD EVALUATE THE FUNCTION C AT THE OUTPUT VALUE U, AND RETURN FU=F(U). C KOUNT - A COUNTER FOR THE NUMBER OF CALLS TO ROOT. KOUNT C SHOULD BE SET TO ZERO ON THE FIRST CALL FOR A GIVEN C ROOT PROBLEM. C IFLAG - AN ERROR RETURN FLAG WHOSE INPUT VALUE IS IMMATERIAL. C C ON RETURN FROM A CALL TO ROOT C C A - ONE ENDPOINT OF CHANGE OF SIGN INTERVAL. C FA - THE VALUE OF F(A). C B - OTHER ENDPOINT OF CHANGE IN SIGN INTERVAL. C FB - THE VALUE OF F(B). C U - CURRENT APPROXIMATION TO THE ROOT. BEFORE ANOTHER C CALL TO ROOT, EVALUATE F(U). C FU - FU WILL BE OVERWRITTEN BY THE USER BEFORE ANOTHER C CALL. ITS VALUE ON RETURN IS ONE OF FA, FU OR FB. C KOUNT - CURRENT NUMBER OF CALLS TO ROOT. C IFLAG - PROGRAM RETURN FLAG C IFLAG=-2 MEANS THAT ON FIRST CALL, FA*FB.GT.0.0. C THIS IS AN ERROR RETURN, SINCE A BRACKETING C INTERVAL SHOULD BE SUPPLIED ON FIRST CALL. C IFLAG=-1 MEANS THAT THE CURRENT BRACKETING INTERVAL C WHOSE ENDPOINTS ARE STORED IN A AND B C IS SO SMALL (LESS THAN 4*EPMACH*ABS(U)+EPMACH) C THAT U SHOULD BE ACCEPTED AS THE ROOT. C THE FUNCTION VALUE F(U) IS STORED IN FU. C IFLAG= 0 MEANS THAT THE INPUT VALUE FU IS EXACTLY C ZERO, AND U SHOULD BE ACCEPTED AS THE ROOT. C C IFLAG.GT.0 MEANS THAT THE CURRENT APPROXIMATION TO C THE ROOT IS CONTAINED IN U. IF A BETTER C APPROXIMATION IS DESIRED, SET FU=F(U) C AND CALL ROOT AGAIN. IFLAG INDICATES C THE METHOD THAT WAS USED TO PRODUCE U. C C IFLAG= 1 BISECTION WAS USED. C IFLAG= 2 LINEAR INTERPOLATION (SECANT METHOD). C IFLAG= 3 INVERSE QUADRATIC INTERPOLATION. C C LOCAL VARIABLES INCLUDE C C EPMACH- SMALLEST POSITIVE NUMBER SUCH THAT 1.0+EPMACH.GT.1.0 C .5*BETA**(1-TAU) FOR ROUNDED, TAU-DIGIT ARITHMETIC C BASE BETA. TWICE THAT VALUE FOR TRUNCATED ARITHMETIC. C THIS IS THE RELATIVE MACHINE PRECISION. C HALFUB- SIGNED HALFWIDTH OF INTERVAL. DURING SEGMENT 3, THE C CHANGE OF SIGN INTERVAL IS (U,B) OR (B,U). THE MIDPOINT C IS XMID=U+HALFUB, REGARDLESS OF ORIENTATION. C SDEL1 - SIZE OF CHANGE IN SIGN INTERVAL. C SDEL2 - PREVIOUS VALUE OF SDEL1. C SDEL3 - PREVIOUS VALUE OF SDEL2. C SDEL4 - PREVIOUS VALUE OF SDEL3. C STEP - THE NEW ROOT IS COMPUTED AS A CORRECTION TO U OF THE C FORM U(NEW)=U(OLD)+STEP. C TOLER - A NUMBER WE ACCEPT AS SMALL WHEN EXAMINING INTERVAL C SIZE OR STEP SIZE. TOLER=2.0*EPMACH*ABS(U) + EPMACH IS C A MINIMUM BELOW WHICH WE DO NOT ALLOW SUCH VALUES TO FALL. C C THIS SUBROUTINE IS CALLED BY C PITCON C AND CALLS C FORTRAN ABS C FORTRAN SIGN C C*********************************************************************** C REAL A, B, U, FA, FB, FU, STEP, TOLER, P, Q, R, S COMMON /TOL/ EPMACH, EPSATE, EPSQRT C C SEGMENT 1 FIRST CALL HANDLED SPECIALLY. DO BOOKKEEPING. C C SET CERTAIN VALUES ONLY FOR INITIAL CALL WITH KOUNT=0 C IF (KOUNT.GT.0) GO TO 10 IF (FA.GT.0.0 .AND. FB.GT.0.0) GO TO 110 IF (FA.LT.0.0 .AND. FB.LT.0.0) GO TO 110 KOUNT = 1 SDEL1 = 2.0*ABS(B-A) SDEL2 = 2.0*SDEL1 SDEL3 = 2.0*SDEL2 U = B FU = FB GO TO 20 C C ON EVERY CALL, INCREMENT COUNTER C 10 KOUNT = KOUNT + 1 C C RETURN IF HIT MACHINE ZERO FOR F(U) C IF (FU.EQ.0.0) GO TO 90 C C SEGMENT 2 REARRANGE POINTS AND FUNCTION VALUES IF C NECESSARY SO THAT FU*FB.LT.0.0, AND SO THAT C ABS(FU).LT.ABS(FB) C IF ((FU.LE.0.0) .AND. (FB.GT.0.0)) GO TO 30 IF ((FU.GT.0.0) .AND. (FB.LE.0.0)) GO TO 30 C C FU AND FB ARE OF SAME SIGN. C (ROOT CHANGED SIGN) C OVERWRITE B WITH VALUE OF A C 20 B = A FB = FA C C IF NECESSARY, SET A =U, U =B, B =U C TO ENSURE THAT ABS(FU).LE.ABS(FB) C 30 IF (ABS(FB).GE.ABS(FU)) GO TO 40 A = U U = B B = A FA = FU FU = FB FB = FA C C SEGMENT 3 CHECK FOR ACCEPTANCE BECAUSE OF SMALL INTERVAL C CURRENT CHANGE IN SIGN INTERVAL IS (B,U) OR (U,B). C 40 TOLER = 2.0*EPMACH*ABS(U) + EPMACH HALFUB = 0.5*(B-U) SDEL4 = SDEL3 SDEL3 = SDEL2 SDEL2 = SDEL1 SDEL1 = ABS(B-U) IF (ABS(HALFUB).LE.TOLER) GO TO 100 C C SEGMENT 4 COMPUTE NEW APPROXIMANT TO ROOT OF THE FORM C U(NEW)=U(OLD)+STEP. C METHODS AVAILABLE ARE LINEAR INTERPOLATION C INVERSE QUADRATIC INTERPOLATION C AND BISECTION. C IF (ABS(FU).GE.ABS(FA)) GO TO 70 IF (A.NE.B) GO TO 50 C C ATTEMPT LINEAR INTERPOLATION IF ONLY TWO POINTS AVAILABLE C COMPUTE P AND Q FOR APPROXIMATION U(NEW)=U(OLD)+P/Q C IFLAG = 2 S = FU/FA P = 2.0*HALFUB*S Q = 1.0 - S GO TO 60 C C ATTEMPT INVERSE QUADRATIC INTERPOLATION IF THREE POINTS AVAILABLE C COMPUTE P AND Q FOR APPROXIMATION U(NEW)=U(OLD)+P/Q C 50 IFLAG = 3 S = FU/FA Q = FA/FB R = FU/FB P = S*(2.0*HALFUB*Q*(Q-R)-(U-A)*(R-1.0)) Q = (Q-1.0)*(R-1.0)*(S-1.0) C C CORRECT THE SIGNS OF P AND Q C 60 IF (P.GT.0.0) Q = -Q P = ABS(P) C C IF P/Q IS TOO LARGE, GO BACK TO BISECTION C IF (8.0*SDEL1.GT.SDEL4) GO TO 70 IF (P.GE.1.5*ABS(HALFUB*Q)-ABS(TOLER*Q)) GO TO 70 STEP = P/Q GO TO 80 C C PERFORM BISECTION C IF ABS(FU).GE.ABS(FA) C OR INTERPOLATION IS UNSAFE (P/Q IS LARGE) C OR IF THREE CONSECUTIVE STEPS HAVE NOT DECREASED C THE SIZE OF THE INTERVAL BY A FACTOR OF 8.0 C 70 IFLAG = 1 STEP = HALFUB GO TO 80 C C SEGMENT 5 VALUE OF STEP HAS BEEN COMPUTED. C UPDATE INFORMATION A =U, FA =FU, U =U+STEP. C CHANGE IN SIGN INTERVAL IS NOW (A,B) OR (B,A). C 80 A = U FA = FU IF (ABS(STEP).LE.TOLER) STEP = SIGN(TOLER,HALFUB) U = U + STEP RETURN C C SPECIAL RETURNS C C INPUT POINT U IS EXACT ROOT C 90 IFLAG = 0 RETURN C C CHANGE IN SIGN INTERVAL IS OF SIZE LESS THAN 4*EPMACH*ABS(U)+EPMACH C INTERVAL RETURNED AS (U,B) OR (B,U). C ACCEPT U AS ROOT WITH RESIDUAL F(U) STORED IN FU. C 100 IFLAG = -1 A = U FA = FU RETURN C C CHANGE OF SIGN CONDITION VIOLATED C 110 IFLAG = -2 KOUNT = 0 RETURN C C*********************************************************************** C END SUBROUTINE FSOLVE(NVAR, X, Y, IP, DET, IEX, IERR, ICALL, MODNEW, FSO 10 * FPRYM, NROW, NCOL, IPIVOT, FPNAME) C C*********************************************************************** C C THIS SUBROUTINE IS CALLED BY C CORECT C TANGNT C AND CALLS C JACOBIAN (FPNAME) C FORTRAN ABS C LINPAK SGEFA C LINPAK SGESL C C THIS SUBROUTINE SOLVES THE LINEAR SYSTEM DFA(X,IP)*Y = Y C WHERE DFA(X,IP) IS THE (NVAR)X(NVAR) MATRIX WHOSE FIRST NVAR - 1 C ROWS ARE THE JACOBIAN OF THE FUNCTION, AND WHOSE LAST C ROW IS ALL 0 EXCEPT FOR A 1 IN THE IP-TH COMPONENT. C C ON INPUT, Y CONTAINS THE RIGHT HAND SIDE. ON OUTPUT, C Y SHOULD CONTAIN THE SOLUTION TO THE SYSTEM. C C **NOTE** SUBROUTINE FSOLVE USES FULL MATRIX STORAGE TO SOLVE THE C LINEAR SYSTEM. THE USER MAY WISH TO REPLACE THIS ROUTINE WITH C ONE MORE SUITABLE. C C DET BINARY MANTISSA OF THE DETERMINANT OF JACOBIAN DFA(X,IP) C IEX BINARY EXPONENT OF THE DETERMINANT OF JACOBIAN DFA(X,IP) C MODNEW NEWTON METHOD FLAG. C MODNEW=0, JACOBIAN IS TO BE EVALUATED EVERY CORRECTOR STEP C AND EVERY TANGENT CALCULATION C MODNEW=1, JACOBIAN IS TO BE EVALUATED FOR FIRST CORRECTOR C STEP ONLY, AND EVERY TANGENT CALCULATION. C ICALL SET UP FLAG. C IF (ICALL.EQ.0.AND.MODNEW.NE.0) DO NOT RE-EVALUATE JACOBIAN C INFO OUTPUT FROM SGEFA. IF INFO.NE.0, SGEFA FOUND A ZERO C PIVOT WHEN ELIMINATING INFO-TH VARIABLE. C IERR RETURN FLAG, 0 MEANS SUCCESSFUL SOLUTION, 1 MEANS FAILURE C NVAR THE NUMBER OF VARIABLES IN THE NONLINEAR SYSTEM C X THE POINT AT WHICH TO EVALUATE FPRYM C Y THE RIGHT HAND SIDE ON INPUT, THE SOLUTION C ON OUTPUT C FPRYM ARRAY WHERE DFA(X,IP) IS TO BE STORED. C IPIVOT INTEGER WORK SPACE FOR PIVOT ROW SWITCHES DEMANDED BY SGEFA C IP THE VARIABLE USED IN THE AUGMENTING EQUATION THAT IS OF THE C FORM X(IP)=B. HENCE THE LAST ROW OF DFA(X,IP) IS ALL C ZERO EXCEPT FOR A 1.0 IN THE IP-TH COLUMN. C NROW NUMBER OF ROWS USED TO STORE FPRYM. C NCOL NUMBER OF COLUMNS USED TO STORE FPRYM. C FPNAME EXTERNAL NAME OF JACOBIAN EVALUATOR. C C*********************************************************************** C EXTERNAL FPNAME REAL X(NVAR), Y(NVAR), FPRYM(NROW,NCOL) INTEGER IPIVOT(NVAR) COMMON /COUNT2/ IFEVAL, IPEVAL, ISOLVE, NREDXF, NRDSUM, NCORXR, * NCRSUM C C DEPENDING ON VALUES OF ICALL AND MODNEW, EITHER SET UP C AUGMENTED JACOBIAN, DECOMPOSE INTO L-U FACTORS, AND GET DETERMINANT, C OR USE CURRENT FACTORED JACOBIAN WITH NEW RIGHT HAND SIDE. C IF (ICALL.EQ.0 .AND. MODNEW.NE.0) GO TO 50 CALL FPNAME(NVAR, X, FPRYM, NROW, NCOL) IPEVAL = IPEVAL + 1 DO 10 I=1,NVAR FPRYM(NVAR,I) = 0.0 10 CONTINUE FPRYM(NVAR,IP) = 1.0 C C CARRY OUT IN CORE LU DECOMPOSITION OF NVAR BY NVAR MATRIX C AND USE PIVOT INFORMATION TO COMPUTE DETERMINANT C CALL SGEFA(FPRYM, NROW, NROW, IPIVOT, INFO) DET = 1.0 IEX = 0 TWO = 2.0 DO 40 I=1,NVAR IF (IPIVOT(I).NE.I) DET = -DET DET = FPRYM(I,I)*DET IF (DET.EQ.0.0) GO TO 60 20 IF (ABS(DET).GE.1.0) GO TO 30 DET = DET*TWO IEX = IEX - 1 GO TO 20 30 IF (ABS(DET).LT.TWO) GO TO 40 DET = DET/TWO IEX = IEX + 1 GO TO 30 40 CONTINUE IF (INFO.NE.0) GO TO 60 C C USING L-U FACTORED MATRIX, SOLVE SYSTEM USING FORWARD-BACKWARD C ELIMINATION, AND OVERWRITE RIGHT HAND SIDE WITH SOLUTION C 50 JOB = 0 CALL SGESL(FPRYM, NROW, NROW, IPIVOT, Y, JOB) IERR = 0 RETURN 60 IERR = 1 INFO = 0 RETURN C C*********************************************************************** C END INTEGER FUNCTION ISAMAX(N, SX, INCX) ISA 10 C C*********************************************************************** C C C FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE. C JACK DONGARRA, LINPACK, 3/11/78. C REAL SX(1), SMAX INTEGER I, INCX, IX, N C ISAMAX = 0 IF (N.LT.1) RETURN ISAMAX = 1 IF (N.EQ.1) RETURN IF (INCX.EQ.1) GO TO 30 C C CODE FOR INCREMENT NOT EQUAL TO 1 C IX = 1 SMAX = ABS(SX(1)) IX = IX + INCX DO 20 I=2,N IF (ABS(SX(IX)).LE.SMAX) GO TO 10 ISAMAX = I SMAX = ABS(SX(IX)) 10 IX = IX + INCX 20 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C 30 SMAX = ABS(SX(1)) DO 40 I=2,N IF (ABS(SX(I)).LE.SMAX) GO TO 40 ISAMAX = I SMAX = ABS(SX(I)) 40 CONTINUE RETURN C C*********************************************************************** C END REAL FUNCTION SASUM(N, SX, INCX) SAS 10 C C*********************************************************************** C C C TAKES THE SUM OF THE ABSOLUTE VALUES. C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C REAL SX(1), STEMP INTEGER I, INCX, M, MP1, N, NINCX C SASUM = 0.0E0 STEMP = 0.0E0 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 STEMP = STEMP + ABS(SX(I)) 10 CONTINUE SASUM = STEMP RETURN C C CODE FOR INCREMENT EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,6) IF (M.EQ.0) GO TO 40 DO 30 I=1,M STEMP = STEMP + ABS(SX(I)) 30 CONTINUE IF (N.LT.6) GO TO 60 40 MP1 = M + 1 DO 50 I=MP1,N,6 STEMP = STEMP + ABS(SX(I)) + ABS(SX(I+1)) + ABS(SX(I+2)) + * ABS(SX(I+3)) + ABS(SX(I+4)) + ABS(SX(I+5)) 50 CONTINUE 60 SASUM = STEMP RETURN C C*********************************************************************** C END SUBROUTINE SAXPY(N, SA, SX, INCX, SY, INCY) SAX 10 C C*********************************************************************** C C C CONSTANT TIMES A VECTOR PLUS A VECTOR. C USES UNROLLED LOOP FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C REAL SX(1), SY(1), SA INTEGER I, INCX, INCY, IX, IY, M, MP1, N C IF (N.LE.0) RETURN IF (SA.EQ.0.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 SY(IY) = SY(IY) + SA*SX(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 SY(I) = SY(I) + SA*SX(I) 30 CONTINUE IF (N.LT.4) RETURN 40 MP1 = M + 1 DO 50 I=MP1,N,4 SY(I) = SY(I) + SA*SX(I) SY(I+1) = SY(I+1) + SA*SX(I+1) SY(I+2) = SY(I+2) + SA*SX(I+2) SY(I+3) = SY(I+3) + SA*SX(I+3) 50 CONTINUE RETURN C C*********************************************************************** C END SUBROUTINE SCOPY(N, SX, INCX, SY, INCY) SCO 10 C C*********************************************************************** C C C COPIES A VECTOR, X, TO A VECTOR, Y. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO 1. C JACK DONGARRA, LINPACK, 3/11/78. C REAL SX(1), SY(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 SY(IY) = SX(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 SY(I) = SX(I) 30 CONTINUE IF (N.LT.7) RETURN 40 MP1 = M + 1 DO 50 I=MP1,N,7 SY(I) = SX(I) SY(I+1) = SX(I+1) SY(I+2) = SX(I+2) SY(I+3) = SX(I+3) SY(I+4) = SX(I+4) SY(I+5) = SX(I+5) SY(I+6) = SX(I+6) 50 CONTINUE RETURN C C*********************************************************************** C END REAL FUNCTION SDOT(N, SX, INCX, SY, INCY) SDO 10 C C*********************************************************************** C 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 REAL SX(1), SY(1), STEMP INTEGER I, INCX, INCY, IX, IY, M, MP1, N C STEMP = 0.0E0 SDOT = 0.0E0 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 STEMP = STEMP + SX(IX)*SY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE SDOT = STEMP 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 STEMP = STEMP + SX(I)*SY(I) 30 CONTINUE IF (N.LT.5) GO TO 60 40 MP1 = M + 1 DO 50 I=MP1,N,5 STEMP = STEMP + SX(I)*SY(I) + SX(I+1)*SY(I+1) + SX(I+2)*SY(I+2) * + SX(I+3)*SY(I+3) + SX(I+4)*SY(I+4) 50 CONTINUE 60 SDOT = STEMP RETURN C C*********************************************************************** C END REAL FUNCTION SNRM2(N, SX, INCX) SNR 10 C C*********************************************************************** C INTEGER NEXT REAL SX(1), CUTLO, CUTHI, HITEST, SUM, XMAX, ZERO, ONE DATA ZERO, ONE /0.0E0,1.0E0/ C C EUCLIDEAN NORM OF THE N-VECTOR STORED IN SX() 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 SQRT(U/EPS) OVER ALL KNOWN MACHINES. C CUTHI = MINIMUM OF SQRT(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 /4.441E-16,1.304E19/ C IF (N.GT.0) GO TO 10 SNRM2 = ZERO GO TO 140 C 10 ASSIGN 30 TO NEXT SUM = ZERO NN = N*INCX C BEGIN MAIN LOOP I = 1 20 GO TO NEXT, (30, 40, 70, 80) 30 IF (ABS(SX(I)).GT.CUTLO) GO TO 110 ASSIGN 40 TO NEXT XMAX = ZERO C C PHASE 1. SUM IS ZERO C 40 IF (SX(I).EQ.ZERO) GO TO 130 IF (ABS(SX(I)).GT.CUTLO) GO TO 110 C C PREPARE FOR PHASE 2. ASSIGN 70 TO NEXT GO TO 60 C C PREPARE FOR PHASE 4. C 50 I = J ASSIGN 80 TO NEXT SUM = (SUM/SX(I))/SX(I) 60 XMAX = ABS(SX(I)) GO TO 90 C C PHASE 2. SUM IS SMALL. C SCALE TO AVOID DESTRUCTIVE UNDERFLOW. C 70 IF (ABS(SX(I)).GT.CUTLO) GO TO 100 C C COMMON CODE FOR PHASES 2 AND 4. C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW. C 80 IF (ABS(SX(I)).LE.XMAX) GO TO 90 SUM = ONE + SUM*(XMAX/SX(I))**2 XMAX = ABS(SX(I)) GO TO 130 C 90 SUM = SUM + (SX(I)/XMAX)**2 GO TO 130 C C C PREPARE FOR PHASE 3. C 100 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 110 HITEST = CUTHI/FLOAT(N) C C PHASE 3. SUM IS MID-RANGE. NO SCALING. C DO 120 J=I,NN,INCX IF (ABS(SX(J)).GE.HITEST) GO TO 50 SUM = SUM + SX(J)**2 120 CONTINUE SNRM2 = SQRT(SUM) GO TO 140 C 130 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 SNRM2 = XMAX*SQRT(SUM) 140 CONTINUE RETURN C C*********************************************************************** C END SUBROUTINE SSCAL(N, SA, SX, INCX) SSC 10 C C*********************************************************************** C C C SCALES A VECTOR BY A CONSTANT. C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO 1. C JACK DONGARRA, LINPACK, 3/11/78. C REAL SA, SX(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 SX(I) = SA*SX(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 SX(I) = SA*SX(I) 30 CONTINUE IF (N.LT.5) RETURN 40 MP1 = M + 1 DO 50 I=MP1,N,5 SX(I) = SA*SX(I) SX(I+1) = SA*SX(I+1) SX(I+2) = SA*SX(I+2) SX(I+3) = SA*SX(I+3) SX(I+4) = SA*SX(I+4) 50 CONTINUE RETURN C C*********************************************************************** C END SUBROUTINE SGEFA(A, LDA, N, IPVT, INFO) SGE 10 C C*********************************************************************** C INTEGER LDA, N, IPVT(1), INFO REAL A(LDA,1) C C SGEFA FACTORS A REAL MATRIX BY GAUSSIAN ELIMINATION. C C SGEFA IS USUALLY CALLED BY SGECO, BUT IT CAN BE CALLED C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. C (TIME FOR SGECO) = (1 + 9/N)*(TIME FOR SGEFA) . C C ON ENTRY C C A REAL(LDA, N) C THE MATRIX TO BE FACTORED. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C ON RETURN C C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS C WHICH WERE USED TO OBTAIN IT. C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. C C IPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C INFO INTEGER C = 0 NORMAL VALUE. C = K IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR C CONDITION FOR THIS SUBROUTINE, BUT IT DOES C INDICATE THAT SGESL OR SGEDI WILL DIVIDE BY ZERO C IF CALLED. USE RCOND IN SGECO FOR A RELIABLE C INDICATION OF SINGULARITY. C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS SAXPY,SSCAL,ISAMAX C C INTERNAL VARIABLES C REAL T INTEGER ISAMAX, J, K, KP1, L, NM1 C C C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING C INFO = 0 NM1 = N - 1 IF (NM1.LT.1) GO TO 70 DO 60 K=1,NM1 KP1 = K + 1 C C FIND L = PIVOT INDEX C L = ISAMAX(N-K+1,A(K,K),1) + K - 1 IPVT(K) = L C C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED C IF (A(L,K).EQ.0.0E0) GO TO 40 C C INTERCHANGE IF NECESSARY C IF (L.EQ.K) GO TO 10 T = A(L,K) A(L,K) = A(K,K) A(K,K) = T 10 CONTINUE C C COMPUTE MULTIPLIERS C T = -1.0E0/A(K,K) CALL SSCAL(N-K, T, A(K+1,K), 1) C C ROW ELIMINATION WITH COLUMN INDEXING C DO 30 J=KP1,N T = A(L,J) IF (L.EQ.K) GO TO 20 A(L,J) = A(K,J) A(K,J) = T 20 CONTINUE CALL SAXPY(N-K, T, A(K+1,K), 1, A(K+1,J), 1) 30 CONTINUE GO TO 50 40 CONTINUE INFO = K 50 CONTINUE 60 CONTINUE 70 CONTINUE IPVT(N) = N IF (A(N,N).EQ.0.0E0) INFO = N RETURN C C*********************************************************************** C END SUBROUTINE SGESL(A, LDA, N, IPVT, B, JOB) SGE 10 C C*********************************************************************** C INTEGER LDA, N, IPVT(1), JOB REAL A(LDA,1), B(1) C C SGESL SOLVES THE REAL SYSTEM C A * X = B OR TRANS(A) * X = B C USING THE FACTORS COMPUTED BY SGECO OR SGEFA. C C ON ENTRY C C A REAL(LDA, N) C THE OUTPUT FROM SGECO OR SGEFA. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C IPVT INTEGER(N) C THE PIVOT VECTOR FROM SGECO OR SGEFA. C C B REAL(N) C THE RIGHT HAND SIDE VECTOR. C C JOB INTEGER C = 0 TO SOLVE A*X = B , C = NONZERO TO SOLVE TRANS(A)*X = B WHERE C TRANS(A) IS THE TRANSPOSE. C C ON RETURN C C B THE SOLUTION VECTOR X . C C ERROR CONDITION C C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER C SETTING OF LDA . IT WILL NOT OCCUR IF THE SUBROUTINES ARE C CALLED CORRECTLY AND IF SGECO HAS SET RCOND .GT. 0.0 C OR SGEFA HAS SET INFO .EQ. 0 . C C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX C WITH P COLUMNS C CALL SGECO(A,LDA,N,IPVT,RCOND,Z) C IF (RCOND IS TOO SMALL) GO TO ... C DO 10 J = 1, P C CALL SGESL(A,LDA,N,IPVT,C(1,J),0) C 10 CONTINUE C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS SAXPY,SDOT C C INTERNAL VARIABLES C REAL SDOT, T INTEGER K, KB, L, NM1 C NM1 = N - 1 IF (JOB.NE.0) GO TO 50 C C JOB = 0 , SOLVE A * X = B C FIRST SOLVE L*Y = B C IF (NM1.LT.1) GO TO 30 DO 20 K=1,NM1 L = IPVT(K) T = B(L) IF (L.EQ.K) GO TO 10 B(L) = B(K) B(K) = T 10 CONTINUE CALL SAXPY(N-K, T, A(K+1,K), 1, B(K+1), 1) 20 CONTINUE 30 CONTINUE C C NOW SOLVE U*X = Y C DO 40 KB=1,N K = N + 1 - KB B(K) = B(K)/A(K,K) T = -B(K) CALL SAXPY(K-1, T, A(1,K), 1, B(1), 1) 40 CONTINUE GO TO 100 50 CONTINUE C C JOB = NONZERO, SOLVE TRANS(A) * X = B C FIRST SOLVE TRANS(U)*Y = B C DO 60 K=1,N T = SDOT(K-1,A(1,K),1,B(1),1) B(K) = (B(K)-T)/A(K,K) 60 CONTINUE C C NOW SOLVE TRANS(L)*X = Y C IF (NM1.LT.1) GO TO 90 DO 80 KB=1,NM1 K = N - KB B(K) = B(K) + SDOT(N-K,A(K+1,K),1,B(K+1),1) L = IPVT(K) IF (L.EQ.K) GO TO 70 T = B(L) B(L) = B(K) B(K) = T 70 CONTINUE 80 CONTINUE 90 CONTINUE 100 CONTINUE RETURN C C*********************************************************************** C END 1AIRCRAFT STABILITY PROBLEM ELEVATOR=-0.5000000E-01 AILERON = 0.0000000E+00 RUDDER = 0.0000000E+00 0UNIVERSITY OF PITTSBURGH CONTINUATION CODE ABSERR= 0.1000000E-04 RELERR= 0.1000000E-04 H= 0.3000000 HFACT= 3.000000 HMIN= 0.1000000E-02 HMAX= 100.0000 NEWTON METHOD OPTION= 0 0LIMIT POINTS SOUGHT FOR VARIABLE 7 ACCEPTED POINT SATISFIED WEAK CRITERION TANGENT AT POINT -1 WAS 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 TO REACH POINT 0 TANGENT STEPSIZE WAS 0.3000000 POINT 0 IS 0.1060012E-02 0.5120612E-01 0.5799534E-04 0.5960608E-01 0.2646838E-04 -0.5000000E-01 0.0000000E+00 0.0000000E+00 FUNCTION NORM= 0.2621267E-07 TANGENT AT POINT 0 WAS 0.9928986 0.3971849E-04 0.5807782E-01 0.3491732E-05 0.9383239E-02 0.0000000E+00-0.1033984 0.0000000E+00 TO REACH POINT 1 TANGENT STEPSIZE WAS 0.3000000 POINT 1 IS 0.2989296 0.5224776E-01 0.1752600E-01 0.5978443E-01 0.2888085E-02 -0.5000000E-01-0.3108907E-01 0.0000000E+00 FUNCTION NORM= 0.1303852E-07 ACCEPTED POINT SATISFIED WEAK CRITERION 0SECANT STEP FROM 0 TO 1 = 0.3000121 TANGENT AT POINT 1 WAS 0.9927709 0.6996797E-02 0.5851547E-01 0.1185679E-02 0.9848173E-02 0.0000000E+00-0.1040892 0.0000000E+00 TO REACH POINT 2 TANGENT STEPSIZE WAS 0.7085235 POINT 2 IS 0.9812587 0.6401174E-01 0.5894670E-01 0.6152920E-01 0.1109039E-01 -0.5000000E-01-0.1048387 0.0000000E+00 FUNCTION NORM= 0.2980232E-07 0SECANT STEP FROM 1 TO 2 = 0.6877037 TANGENT AT POINT 2 WAS 0.9911383 0.2978522E-01 0.6259765E-01 0.3870231E-02 0.1520790E-01 0.0000000E+00-0.1122183 0.0000000E+00 TO REACH POINT 3 TANGENT STEPSIZE WAS 2.063111 POINT 3 IS 2.440764 0.2531042 0.1521706 0.6774194E-01 0.7938589E-01 -0.5000000E-01-0.3363575 0.0000000E+00 FUNCTION NORM= 0.5960464E-07 0SECANT STEP FROM 2 TO 3 = 1.494291 TANGENT AT POINT 3 WAS 0.8966243 0.3538527 0.3001525E-01-0.3850312E-02 0.1173910 0.0000000E+00-0.2369742 0.0000000E+00 TO REACH POINT 4 TANGENT STEPSIZE WAS 1.496469 POINT 4 IS 2.952565 0.7826338 0.8401720E-01 0.4403858E-01 0.2529988 -0.5000000E-01-0.5039722 0.0000000E+00 FUNCTION NORM= 0.1192093E-06 0SECANT STEP FROM 3 TO 4 = 0.7783215 TANGENT AT POINT 4 WAS 0.2754016 0.8876920 -0.2103540 -0.5578997E-01 0.2959480 -0.3963379E-10-0.3478882E-01 0.0000000E+00 TO REACH POINT 5 TANGENT STEPSIZE WAS 1.444062 POINT 5 IS 2.916958 2.064516 -0.2698410 -0.4357284E-01 0.7255810 -0.5000000E-01 0.7422988E-01 0.0000000E+00 FUNCTION NORM= 0.5960464E-07 0LIMIT POINT FOUND POINT 2.964867 0.8255649 0.7366088E-01 0.4130952E-01 0.2673494 -0.5000000E-01-0.5048105 0.0000000E+00 TANGENT 0.2362643 0.8957867 -0.2197420 -0.5755886E-01 0.3002430 0.0000000E+00-0.4514060E-07 0.0000000E+00 FUNCTION NORM= 0.3576279E-06 0SECANT STEP FROM 4 TO 5 = 1.528081 TANGENT AT POINT 5 WAS -0.1237267 0.7033410 -0.1898283 -0.4611499E-01 0.2877072 0.0000000E+00 0.6075083 0.0000000E+00 TO REACH POINT 6 TANGENT STEPSIZE WAS 4.584244 POINT 6 IS 2.496304 4.058408 -0.7461970 -0.1579089 1.692389 -0.5000000E-01 2.859196 0.0000000E+00 FUNCTION NORM= 0.5960464E-07 0SECANT STEP FROM 5 TO 6 = 3.617084 TANGENT AT POINT 6 WAS -0.9773980E-01 0.4392467 -0.9534489E-01-0.2223779E-01 0.2511305 0.0000000E+00 0.8513857 0.0000000E+00 TO REACH POINT 7 TANGENT STEPSIZE WAS 10.85125 POINT 7 IS 1.807083 7.067870 -1.413138 -0.2922075 4.077554 -0.5000000E-01 12.09780 0.0000000E+00 FUNCTION NORM= 0.7629395E-05 0SECANT STEP FROM 6 TO 7 = 10.05164 TANGENT AT POINT 7 WAS -0.4952314E-01 0.2046322 -0.5205488E-01-0.8192124E-02 0.2295179 0.0000000E+00 0.9487981 0.0000000E+00 TO REACH POINT 8 TANGENT STEPSIZE WAS 30.15492 POINT 8 IS 0.9593399 10.01302 -2.803647 -0.4014672 10.86465 -0.5000000E-01 40.70873 0.0000000E+00 FUNCTION NORM= 0.1192093E-06 0SECANT STEP FROM 7 TO 8 = 29.59710 TANGENT AT POINT 8 WAS -0.1635089E-01 0.4301921E-01-0.4714328E-01-0.1489147E-02 0.2315707 0.0000000E+00 0.9705835 0.0000000E+00 TO REACH POINT 9 TANGENT STEPSIZE WAS 88.79130 POINT 9 IS 0.3628370 11.12638 -7.312350 -0.4391644 31.89847 -0.5000000E-01 126.8881 0.0000000E+00 FUNCTION NORM= 0.9536743E-06 0SECANT STEP FROM 8 TO 9 = 88.83260 TANGENT AT POINT 9 WAS -0.2632842E-02 0.2777145E-02-0.5267731E-01-0.9206447E-04 0.2393715 0.0000000E+00 0.9694905 0.0000000E+00 TO REACH POINT 10 TANGENT STEPSIZE WAS 100.0000 POINT 10 IS 0.2096539 11.25436 -12.62353 -0.4433980 55.83562 -0.5000000E-01 223.4916 0.0000000E+00 FUNCTION NORM= 0.7450581E-07 0SECANT STEP FROM 9 TO 10 = 99.66683 TANGENT AT POINT 10 WAS -0.8934377E-03 0.5492248E-03-0.5365569E-01-0.1813405E-04 0.2406490 0.0000000E+00 0.9691275 0.0000000E+00 TO REACH POINT 11 TANGENT STEPSIZE WAS 100.0000 POINT 11 IS 0.1469338 11.28720 -18.00007 -0.4444822 79.90052 -0.5000000E-01 320.3178 0.0000000E+00 FUNCTION NORM= 0.3814697E-05 0SECANT STEP FROM 10 TO 11 = 99.91665 TANGENT AT POINT 11 WAS -0.4407084E-03 0.1906852E-03-0.5391949E-01-0.6291352E-05 0.2409916 0.0000000E+00 0.9690281 0.0000000E+00 TO REACH POINT 12 TANGENT STEPSIZE WAS 100.0000 POINT 12 IS 0.1130149 11.30020 -23.39625 -0.4449112 103.9997 -0.5000000E-01 417.1873 0.0000000E+00 FUNCTION NORM= 0.7629395E-05 0SECANT STEP FROM 11 TO 12 = 99.96792 TANGENT AT POINT 12 WAS -0.2611652E-03 0.8714148E-04-0.5402515E-01-0.2877278E-05 0.2411286 0.0000000E+00 0.9689882 0.0000000E+00 TO REACH POINT 13 TANGENT STEPSIZE WAS 33.33333 POINT 13 IS 0.1049311 11.30280 -25.19735 -0.4449971 112.0373 -0.5000000E-01 449.4848 0.0000000E+00 FUNCTION NORM= 0.3814697E-05 0SECANT STEP FROM 12 TO 13 = 33.33136 TANGENT AT POINT 13 WAS -0.2252165E-03 0.6982322E-04-0.5404638E-01-0.2306299E-05 0.2411561 0.0000000E+00 0.9689802 0.0000000E+00 TO REACH POINT 14 TANGENT STEPSIZE WAS 33.33327 POINT 14 IS 0.9792452E-01 11.30490 -26.99910 -0.4450665 120.0758 -0.5000000E-01 481.7825 0.0000000E+00 FUNCTION NORM= 0.5722046E-05 0SECANT STEP FROM 13 TO 14 = 33.33168 TANGENT AT POINT 14 WAS -0.1961974E-03 0.5683755E-04-0.5406353E-01-0.1876989E-05 0.2411783 0.0000000E+00 0.9689737 0.0000000E+00 TO REACH POINT 15 TANGENT STEPSIZE WAS 47.14118 POINT 15 IS 0.8947291E-01 11.30724 -29.54805 -0.4451439 131.4452 -0.5000000E-01 527.4584 0.0000000E+00 FUNCTION NORM= 0.3129244E-06 0SECANT STEP FROM 14 TO 15 = 47.13864 TANGENT AT POINT 15 WAS -0.1638424E-03 0.4340454E-04-0.5408267E-01-0.1434397E-05 0.2412031 0.0000000E+00 0.9689665 0.0000000E+00 TO REACH POINT 16 TANGENT STEPSIZE WAS 66.66896 POINT 16 IS 0.7973719E-01 11.30969 -33.15418 -0.4452246 147.5260 -0.5000000E-01 592.0545 0.0000000E+00 FUNCTION NORM= 0.7629395E-05 0SECANT STEP FROM 15 TO 16 = 66.66518 TANGENT AT POINT 16 WAS -0.1301677E-03 0.3078870E-04-0.5410261E-01-0.1017784E-05 0.2412290 0.0000000E+00 0.9689590 0.0000000E+00 TO REACH POINT 17 TANGENT STEPSIZE WAS 94.28594 POINT 17 IS 0.6910028E-01 11.31204 -38.25599 -0.4453023 170.2705 -0.5000000E-01 683.4083 0.0000000E+00 FUNCTION NORM= 0.1028180E-05 0SECANT STEP FROM 16 TO 17 = 94.28075 TANGENT AT POINT 17 WAS -0.9778527E-04 0.2009472E-04-0.5412181E-01-0.6646108E-06 0.2412538 0.0000000E+00 0.9689517 0.0000000E+00 TO REACH POINT 18 TANGENT STEPSIZE WAS 33.33333 POINT 18 IS 0.6598751E-01 11.31266 -40.06011 -0.4453230 178.3123 -0.5000000E-01 715.7062 0.0000000E+00 FUNCTION NORM= 0.3814697E-05 0SECANT STEP FROM 17 TO 18 = 33.33286 TANGENT AT POINT 18 WAS -0.8918105E-04 0.1751778E-04-0.5412691E-01-0.5794362E-06 0.2412604 0.0000000E+00 0.9689498 0.0000000E+00 TO REACH POINT 19 TANGENT STEPSIZE WAS 33.33332 POINT 19 IS 0.6314289E-01 11.31321 -41.86439 -0.4453411 186.3543 -0.5000000E-01 748.0041 0.0000000E+00 FUNCTION NORM= 0.3814697E-05 0SECANT STEP FROM 18 TO 19 = 33.33291 TANGENT AT POINT 19 WAS -0.8166372E-04 0.1535104E-04-0.5413137E-01-0.5082582E-06 0.2412662 0.0000000E+00 0.9689481 0.0000000E+00 TO REACH POINT 20 TANGENT STEPSIZE WAS 47.14066 POINT 20 IS 0.5951433E-01 11.31387 -44.41628 -0.4453630 197.7278 -0.5000000E-01 793.6802 0.0000000E+00 FUNCTION NORM= 0.1902226E-06 0SECANT STEP FROM 19 TO 20 = 47.13995 TANGENT AT POINT 20 WAS -0.7255367E-04 0.1291316E-04-0.5413678E-01-0.4262529E-06 0.2412732 0.0000000E+00 0.9689460 0.0000000E+00 TO REACH POINT 21 TANGENT STEPSIZE WAS 66.66732 POINT 21 IS 0.5504079E-01 11.31464 -48.02558 -0.4453884 213.8128 -0.5000000E-01 858.2760 0.0000000E+00 FUNCTION NORM= 0.3814697E-05 0SECANT STEP FROM 20 TO 21 = 66.66616 TANGENT AT POINT 21 WAS -0.6206258E-04 0.1019965E-04-0.5414301E-01-0.3379401E-06 0.2412813 0.0000000E+00 0.9689437 0.0000000E+00 TO REACH POINT 22 TANGENT STEPSIZE WAS 31.42747 POINT 22 IS 0.5315705E-01 11.31494 -49.72719 -0.4453984 221.3957 -0.5000000E-01 888.7272 0.0000000E+00 FUNCTION NORM= 0.3814697E-05 0SECANT STEP FROM 21 TO 22 = 31.42726 TANGENT AT POINT 22 WAS -0.5788949E-04 0.9184560E-05-0.5414549E-01-0.3047140E-06 0.2412845 0.0000000E+00 0.9689428 0.0000000E+00 TO REACH POINT 23 TANGENT STEPSIZE WAS 31.42746 POINT 23 IS 0.5139791E-01 11.31522 -51.42887 -0.4454076 228.9786 -0.5000000E-01 919.1784 0.0000000E+00 FUNCTION NORM= 0.3814697E-05 0PROGRAM TOOK FULL NUMBER OF CONTINUATION STEPS 0MATRIX SIZE 8X 8 STEPS TAKEN 23 NEWTON STEPS 77 STEP REDUCTIONS 3 FUNCTION CALLS 107 PRIME CALLS 100 SOLVE CALLS 103 0CALLS - SOLVE CALLED FOR CORRECTION 77 SOLVE CALLED FOR TANGENT 26 0CORRECTOR CALLED FOR STARTING POINT 2 CORRECTOR CALLED FOR CONTINUATION POINT 69 CORRECTOR CALLED FOR TARGET POINT 0 CORRECTOR CALLED FOR LIMIT POINT 6 0ROOTFINDER CALLED FOR LIMIT POINT 3