SUBROUTINE PITCON(NVAR, LIM, IT, XIT, KSTEP, IPC, IPCFIX, DIRIPC, & HTANCF, IRET, MODCON, IPIVOT, HMAX, HMIN, HFACT, ABSERR, RELERR, & RWORK, ISIZE, NROW, NCOL, FXNAME, FPNAME, SLNAME, LUNIT) c*********************************************************************72 c cc pitcon() is the driver routine for the continuation package. 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 Reference: c c Werner Rheinboldt, c Solution Field of Nonlinear Equations and Continuation Methods, c SIAM Journal on Numerical Analysis, c Volume 17, Number 2, April 1980, pages 221-237. c c Cor den Heijer, Werner Rheinboldt, c On Steplength Algorithms for a Class of Continuation Methods, c SIAM Journal on Numerical Analysis c Volume 18, Number 5, October 1981, pages 925-947. c c Werner Rheinboldt, c Numerical Analysis of Continuation Methods for Nonlinear c Structural Problems, c Computers and Structures, c Volume 13, 1981, pages 103-114. c c Werner Rheinboldt, John Burkardt, c A Locally Parameterized Continuation Process, c ACM Transactions on Mathematical Software, c Volume 9, Number 2, June 1983, pages 215-235. c c Werner Rheinboldt, John Burkardt, c Algorithm 596: c A Program for a Locally Parameterized c Continuation Process, c ACM Transactions on Mathematical Software, c Volume 9, Number 2, June 1983, pages 236-241. 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 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 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 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 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 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 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 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 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 real acof(12) real alfmin double precision cosalf double precision dadjus double precision dtcipc double precision dtlipc real epmach real epsate real epsqrt external fpname external fxname integer idone integer ierr integer ifp integer iret integer ipivot(nvar) integer itc integer itl integer ixc integer ixf integer ixr integer jfp integer jtc integer jtl integer jxc integer jxf integer jxr integer kstep integer kstepl real rmach3 real rwork(isize) external slname real tcos real tenm1 real tenm2 real tenm3 real tsin real wrge(8) 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 save save acof save alfmin save idone save ifp save itc save itl save ixc save ixf save ixr save jfp save jtc save jtl save jxc save jxf save jxr save kstepl save rmach3 save tcos save tenm1 save tenm2 save tenm3 save tsin save wrge data idone /0/ data rmach3 /0.59604644775E-07/ data tenm1 /0.1/ data tenm2 /0.01/ data tenm3 /0.001/ 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 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 continue 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 continue IF (KSTEP.EQ.0) then CALL SCOPY(NVAR, RWORK(IXR), 1, RWORK(IXC), 1) CALL SCOPY(NVAR, RWORK(IXR), 1, RWORK(IXF), 1) end if 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 30 continue 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 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 40 continue 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 continue 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 continue 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 continue 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 continue IF (KSTEP.LE.0) GO TO 200 COSALF = 0.0D0 DO 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) end do 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 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 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 continue CALL SCOPY(NVAR, RWORK(IXC), 1, RWORK(IXR), 1) GO TO 350 110 continue ATCLIM = ABS(TCLIM) IF (ATCLIM.GT.0.5*ABSERR) GO TO 130 120 continue 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 continue 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 continue 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 continue 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 continue 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 continue 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 continue MODLIM = 1 190 continue 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 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 200 continue c c CHECK IF DEFAULT STEP MUST BE USED c IF PREVIOUS STEP WAS OF SIZE ZERO, USE STEPSIZE HDEF=(HMIN+HMAX)/2 c IF (KSTEP.GT.0 .AND. HSECLC.GT.0.0) GO TO 210 IF (KSTEP.GT.0) HTANCF = HDEF GO TO 230 210 continue 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 continue 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 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 230 continue KSTEPL = KSTEP KSTEP = KSTEP + 1 NREDXF = 0 240 continue CALL SCOPY(NVAR, RWORK(IXF), 1, RWORK(IXR), 1) CALL SAXPY(NVAR, HTANCF, RWORK(ITC), 1, RWORK(IXR), 1) 250 continue 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 continue 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 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 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 RETURNS. SET VALUE OF IRET. IF AN ERROR OCCURRED, PRINT c A MESSAGE AS WELL. 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 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, & NROW, NCOL, IPIVOT, ABSERR, RELERR, XSTEP, NEQN, FNRM, FXNAME, & FPNAME, SLNAME, LUNIT, DET, IEX) c*********************************************************************72 c cc CORECT corrects an approximate solution. 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 EXTERNAL FXNAME, FPNAME, SLNAME REAL X(NVAR), WORK(NVAR), FPRYM(NROW,NCOL) INTEGER IPIVOT(NVAR) save 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 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 end do 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 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, & NEQN, DET, IEX, IPCFIX, ABSERR, RELERR, NROW, NCOL, FPNAME, & SLNAME) c*********************************************************************72 c cc TANGNT computes a unit tangent vector at a point on the solution curve. 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 EXTERNAL FPNAME, SLNAME REAL X(NVAR), TAN(NVAR), FPRYM(NROW,NCOL) INTEGER IPIVOT(NVAR) save 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 I=1,NEQN TAN(I) = 0.0 end do 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 END SUBROUTINE ROOT(A, FA, B, FB, U, FU, KOUNT, IFLAG) c*********************************************************************72 c cc ROOT seeks a root of a scalar nonlinear equation in a change-of-sign interval. c c ROOT is 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 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 END SUBROUTINE FSOLVE(NVAR, X, Y, IP, DET, IEX, IERR, ICALL, MODNEW, & FPRYM, NROW, NCOL, IPIVOT, FPNAME) c*********************************************************************72 c cc FSOLVE solves a linear system related to the continuation process. 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 that 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 Parameters: c c NVAR THE NUMBER OF VARIABLES IN THE NONLINEAR SYSTEM. c c X THE POINT AT WHICH TO EVALUATE FPRYM. c c Y THE RIGHT HAND SIDE ON INPUT, THE SOLUTION c c Input, integer IP. THE VARIABLE USED IN THE AUGMENTING EQUATION THAT c IS OF THE 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 c DET BINARY MANTISSA OF THE DETERMINANT OF JACOBIAN DFA(X,IP). c c IEX BINARY EXPONENT OF THE DETERMINANT OF JACOBIAN DFA(X,IP). c c IERR RETURN FLAG, 0 MEANS SUCCESSFUL SOLUTION, 1 MEANS FAILURE c c ICALL SET UP FLAG. c IF (ICALL.EQ.0.AND.MODNEW.NE.0) DO NOT RE-EVALUATE JACOBIAN. c c MODNEW NEWTON METHOD FLAG. c 0, JACOBIAN IS TO BE EVALUATED EVERY CORRECTOR STEP c AND EVERY TANGENT CALCULATION c 1, JACOBIAN IS TO BE EVALUATED FOR FIRST CORRECTOR c STEP ONLY, AND EVERY TANGENT CALCULATION. c c FPRYM ARRAY WHERE DFA(X,IP) IS TO BE STORED. c c NROW NUMBER OF ROWS USED TO STORE FPRYM. c c NCOL NUMBER OF COLUMNS USED TO STORE FPRYM. c c IPIVOT INTEGER WORK SPACE FOR PIVOT ROW SWITCHES DEMANDED BY SGEFA c c FPNAME EXTERNAL NAME OF JACOBIAN EVALUATOR. c implicit none integer ncol integer nrow integer nvar real det external fpname real fprym(nrow,ncol) integer i integer icall integer ierr integer iex integer ifeval integer info integer ip integer ipeval integer ipivot(nvar) integer isolve integer job integer modnew integer ncorxr integer ncrsum integer nrdsum integer nredxf real two real x(nvar) real y(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 I=1,NVAR FPRYM(NVAR,I) = 0.0 end do 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 END INTEGER FUNCTION ISAMAX(N, SX, INCX) c*********************************************************************72 c cc ISAMAX returns the index of the vector element of maximum absolute value. 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 I=2,N IF (ABS(SX(IX)).LE.SMAX) GO TO 10 ISAMAX = I SMAX = ABS(SX(IX)) 10 IX = IX + INCX end do 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 END REAL FUNCTION SASUM(N, SX, INCX) c*********************************************************************72 c cc SASUM returns the sum of the absolute values of the entries of a vector. 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 I=1,NINCX,INCX STEMP = STEMP + ABS(SX(I)) end do 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 I=1,M STEMP = STEMP + ABS(SX(I)) end do IF (N.LT.6) GO TO 60 40 MP1 = M + 1 DO 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)) end do 60 SASUM = STEMP RETURN END SUBROUTINE SAXPY(N, SA, SX, INCX, SY, INCY) c*********************************************************************72 c cc SAXPY adds a multiple of one vector to another. 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 I=1,N SY(IY) = SY(IY) + SA*SX(IX) IX = IX + INCX IY = IY + INCY end do 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 I=1,M SY(I) = SY(I) + SA*SX(I) end do IF (N.LT.4) RETURN 40 MP1 = M + 1 DO 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) end do RETURN END SUBROUTINE SCOPY(N, SX, INCX, SY, INCY) c*********************************************************************72 c cc SCOPY copies a vector. 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 I=1,N SY(IY) = SX(IX) IX = IX + INCX IY = IY + INCY end do 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 I=1,M SY(I) = SX(I) end do IF (N.LT.7) RETURN 40 MP1 = M + 1 DO 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) end do RETURN END REAL FUNCTION SDOT(N, SX, INCX, SY, INCY) c*********************************************************************72 c cc SDOT computes the dot product of two vectors. 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 I=1,N STEMP = STEMP + SX(IX)*SY(IY) IX = IX + INCX IY = IY + INCY end do 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 I=1,M STEMP = STEMP + SX(I)*SY(I) end do IF (N.LT.5) GO TO 60 40 MP1 = M + 1 DO 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) end do 60 SDOT = STEMP RETURN END REAL FUNCTION SNRM2(N, SX, INCX) c*********************************************************************72 c cc SNRM2 computes the Euclidean norm of a vector. 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 J=I,NN,INCX IF (ABS(SX(J)).GE.HITEST) GO TO 50 SUM = SUM + SX(J)**2 end do SNRM2 = SQRT(SUM) GO TO 140 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 END SUBROUTINE SSCAL(N, SA, SX, INCX) c*********************************************************************72 c cc SSCAL multiplies a vector by a scale factor. 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 I=1,NINCX,INCX SX(I) = SA*SX(I) end do 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 I=1,M SX(I) = SA*SX(I) end do IF (N.LT.5) RETURN 40 MP1 = M + 1 DO 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) end do RETURN END SUBROUTINE SGEFA(A, LDA, N, IPVT, INFO) c*********************************************************************72 c cc SGEFA factors a real matrix by Gaussian elimination. 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 END SUBROUTINE SGESL(A, LDA, N, IPVT, B, JOB) c*********************************************************************72 c cc SGESL solves a linear system that was factored by SGEFA. 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 K=1,N T = SDOT(K-1,A(1,K),1,B(1),1) B(K) = (B(K)-T)/A(K,K) end do 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 END