PROGRAM ABCALL PARAMETER(MAXDEG=13,MXDM=51) C N, THE NUMBER OF EQUATIONS IN THE SYSTEM F(U,R) = 0 C MAXDEG, MAXIMUM DEGREE PREDICTOR FOR WHICH STORAGE IS ALLOCATED C MXDM, THE LEADING DIMENSION OF TWO DIMENSIONAL ARRAYS D,PHI,PHIA C STORAGE IS ALLOCATED FOR SYSTEMS OF UP TO MXDM-1 EQUATIONS C******************************************************************** C C ** TYPE AND DIMENSION** REAL WC(MXDM),WCP(MXDM),WP(MXDM),DIR(MXDM),U(MXDM),FFIN(MXDM), + V(MXDM),STO(MXDM),STO2(MXDM),D(MXDM,MXDM),TANGNT(MXDM), C ARRAYS FOR PREDICTOR AND SCALED DIVIDED DIFFERENCE STORAGE + G(0:MAXDEG,1:MAXDEG),PSYP(0:MAXDEG),PSY(0:MAXDEG), + A(0:MAXDEG),B(0:MAXDEG),PHI(MXDM,0:MAXDEG), + PSYARC(0:MAXDEG),PHIA(MXDM,0:MAXDEG), C REAL SPECIFICATION VARIABLES + EREL,EABS,ECOR,EFIN,HMIN,HMAX,RLOW,RHIGH,ORIENT,DS C INTEGER IPIVOT(MXDM),IFLAG,KNTROL,LEX,MXDG,MXSTEP, + NPRNT,M,NEWTON,NFAC,IFIN,NP1,MPREF,COUNT,ITOK C EXTERNAL F1,DZF1,DRF1 EXTERNAL F2,DZF2,DRF2 EXTERNAL F3,DZF3,DRF3 EXTERNAL F4,DZF4,DRF4 EXTERNAL F5,DZF5,DRF5 EXTERNAL F6,DZF6,DRF6 EXTERNAL F7,DZF7,DRF7 EXTERNAL F8,DZF8,DRF8 C C INITIALIZE PARAMETERS AND INPUT VARIABLES C *********** PROBLEM 1 WATSTON N=10 ***************************** PRINT*,'PROBLEM 1 WATSON N = 10' N=10 C ENTER INITIAL SOLUTION POINT TO F(W) = 0. DO 29 L=1,N+1 WC(L) = 0.E0 29 CONTINUE C** PARAMETERS ******************** C THE BOUNDS OF CONTINUATION PARAMETER: THE CONTINUATION WILL C STOP WITH A POINT WC WITH WC(N+1) = RHIGH, (OR RLOW). RHIGH = 1. RLOW = -20. C MINIMUM PREDICTOR STEPSIZE HMIN = 1.E-6 C MAXIMUM PREDICTOR STEPSIZE HMAX = 100. C ABSOLUTE AND RELATIVE PREDICTOR ERROR PARAMETERS EABS = 1.E-2 EREL = 1.E-2 C ERROR PARAMETER FOR CORRECTED POINTS ECOR = 1.E-4 C ERROR PARAMETER FOR FINAL TARGET POINT: C IF EFIN .GT. 0, USE A CHORD CORRECTOR, STOPPING ERROR EFIN C IF EFIN .LT. 0, USE A NEWTON CORRECTOR, STOPPING ERROR -EFIN C IF EFIN .EQ. 0, USE NO CORRECTOR EFIN = +1.0E-4 C CONTROL PARAMETER FOR FINAL POINT CORRECTION C IFIN .GT. 0, FINAL CORRECTION DONE WITH A NEWTON OR CHORD ITERATION C WITH TARGET PARAMETER FIXED AND NO BORDERING AND .LE IFIN C ITERATIONS. C IFIN .LT. 0, USUAL BORDERING CORRECTOR USED BUT WITH DIR RESET C TO GIVE CORRECTION WITH TARGET PARAMETER FIXED C AND .LE. -IFIN ITERATIONS C IFIN .EQ 0, USUAL CORRECTION FOR POINTS; BORDERING ALGORITHM C WITH ITERATES IN A HYP. PLANE ORTH. TO USUAL PREDICTOR. C MAXIMUM ITERATES DETERMINED AS WITH POINTS (SEE NEWTON) IFIN = -1 C SET MAXIMUM ITERATIONS IN CORRECTOR: C IF NEWTON = 1, 2, ... CORRECTOR ABORTS AFTER NEWTON CHORD ITERATES C IF NEWTON = -1, -2, ... CORRECTOR USES UP TO "NEWTON" NEWTON ITERATES NEWTON = 7 C SET MAX CORRECTOR ITERATES FOR INCREASING STEPSIZE AFTER FAILURES ITOK = 7 C SET PRINT CODE: C IF NPRNT=0, NO PRINTING OF INFORMATION C IF NPRNT=1, INITIAL HEADING WITH PARAMETER SETTINGS PRINTED C FIRST TWO STEPS PRINTED ALONG WITH SINGULAR POINT C INFO AND FINAL POINT INFO C IF NPRNT=2, ABOVE PLUS CORRECTED POINTS (FIRST TWO, EVERY KSKIP STEP) C IF NPRNT=3, ABOVE PLUS HEADING FOR STEPS AND OTHER STEP INFO. NPRNT = -3 C SET SKIPPING PARAMETER: PRINT EVERY KSKIP STEP. KSKIP = 0 C SET THE MINIMUM NUMBER OF STEPS BEFORE FINAL POINT STOPPING C WILL BE INVOKED. THIS CAN BE USED TO PICK OUT SUBSEQUENT C PASSINGS OF A PARAMETER VALUE FOR ATTENTION. C IF KTARG IS A NEGATIVE INTEGER, FINAL PREDICTED POINT RETURNED IN WC KTARG = 1 C SET ERROR TOLERANCE FOR MAX NORM OF F(W) AT INITIAL POINT C IF THIS TOLERANCE IS NOT MET, THE CODE WILL ABORT. TOLINT = EABS C SET ORIENTATION FOR INITIAL DIRECTION OF CONTINUATION C IF ORIENT = +1, THE PARAMETER WC(N+1) INITIALLY IS INCREASING C IF ORIENT = -1, THE PARAMETER INITIALLY DECREASES ORIENT = +1. C SET THE INDEX OF THE TARGET COMPONENT FOR STOPPING THE CONTINUATION: C IF NCOMP = 1, ..., N+1, THE CONTINUATION WILL STOP AND LOCATE A POINT C WITH WC(NCOMP) = VCOMP ONCE WC(NCOMP) .GT. VCOMP C IF NCOMP = -1, ...,-(N+1) THE CONTINUATION WILL STOP AND LOCATE A C POINT WITH WC(-NCOMP) = VCOMP ONCE WC(-NCOMP) .LT. VCOMP C IF NCOMP = 0, THIS FEATURE IS DISABLED (SEE RLOW, RHIGH HOWEVER) NCOMP = N+1 C SET TARGET VALUE FOR WC(ABS(MCOMP)) (SEE NCOMP) VCOMP = 1.0 C SET MAXIMUM NUMBER OF CONTINUATION STEPS (POINTS) MXSTEP = 600 C AUTOMATIC PREDICTOR ORDER INCREASE UNTIL ORDER > MPREF MPREF = 0 C MAXIMUM DEGREE OF PREDICTOR POLYNOMIAL MXDG = 3 C NFAC = 0 FOR TANGENT COMPUTED AT PREDICTED POINT C NFAC = 1 FOR TANGENT COMPUTED AT CORRECTED POINT NFAC = 0 C INITIAL STEPSIZE. IF DS .LT.0, DS IS COMPUTED BY ABCON DS = -8.E-2 C ISING = 0 TURNS OFF BIFURCATION CHECKING (FOLDS ONLY) C ISING = 1 TURNS ON BUFURCATION POINT CHECKING ISING = 0 C CALL ABCON(N,WC,FFIN,F8,DZF8,DRF8,NPRNT,KSKIP,ORIENT,NCOMP, + VCOMP,RHIGH,RLOW,KTARG,IFIN,EABS,EREL,ECOR,EFIN,DS,HMIN,HMAX, + MXSTEP,MXDG,MPREF,NFAC,NEWTON,ITOK,TOLINT,MXDM,D,IPIVOT,WCP,WP, + TANGNT,DIR,U,V,STO,STO2,PHI,PHIA,ISING,IRCODE) C CP print*,'IRCODE =',IRCODE,'wc=',(wc(i),i=1,n+1) print*,'=====================================================' C###################################################################### C GO TO 1000 C *********** PROBLEM 2 WATSON CURVE N = 12 ************** PRINT*,'*********** PROBLEM 2 WATSON CURVE N = 12 **********' N=12 C ENTER INITIAL SOLUTION POINT TO F(W) = 0. DO 229 L=1,N+1 WC(L) = 0.E0 229 CONTINUE C** PARAMETERS ******************** RHIGH = 1. RLOW = -20. HMIN = 1.E-6 HMAX = 100. EABS = 1.E-2 EREL = 1.E-2 ECOR = 1.E-4 EFIN = +1.0E-4 IFIN = -1 NEWTON = 7 ITOK = 7 NPRNT = -3 KSKIP = 0 KTARG = 1 TOLINT = EABS ORIENT = +1. NCOMP = N+1 VCOMP = 1.0 MXSTEP = 800 MPREF = 1 C MXDG = 3 NFAC = 0 DS = -8.E-2 ISING = 0 C CALL ABCON(N,WC,FFIN,F8,DZF8,DRF8,NPRNT,KSKIP,ORIENT,NCOMP, + VCOMP,RHIGH,RLOW,KTARG,IFIN,EABS,EREL,ECOR,EFIN,DS,HMIN,HMAX, + MXSTEP,MXDG,MPREF,NFAC,NEWTON,ITOK,TOLINT,MXDM,D,IPIVOT,WCP,WP, + TANGNT,DIR,U,V,STO,STO2,PHI,PHIA,ISING,IRCODE) C CP print*,'IRCODE =',IRCODE,'wc=',(wc(i),i=1,n+1) print*,'=====================================================' C###################################################################### C *********** PROBLEM 3 (SCHWETLICK 3) N = 4 ************** PRINT*,'****** PROBLEM 3 (SCHWETLICK 3) N = 4 ********' N=4 C ENTER INITIAL SOLUTION POINT TO F(W) = 0. DO 329 L=1,N+1 WC(L) = 0.E0 329 CONTINUE WC(1) = -3.E0 WC(2) = -1.E0 WC(3) = -3.E0 WC(4) = -1.E0 C** PARAMETERS ******************** RHIGH = 1. RLOW = -20. HMIN = 1.E-8 HMAX = 100. EABS = 1.E-5 EREL = 1.E-5 ECOR = 1.E-7 EFIN = +1.0E-7 IFIN = -1 NEWTON = 7 ITOK = 7 NPRNT = -3 KSKIP = 0 KTARG = 1 TOLINT = EABS ORIENT = +1. NCOMP = N+1 VCOMP = 1.0 MXSTEP = 300 MPREF = 0 C MXDG = 3 NFAC = 0 DS = -8.E-2 C CALL ABCON(N,WC,FFIN,F1,DZF1,DRF1,NPRNT,KSKIP,ORIENT,NCOMP, + VCOMP,RHIGH,RLOW,KTARG,IFIN,EABS,EREL,ECOR,EFIN,DS,HMIN,HMAX, + MXSTEP,MXDG,MPREF,NFAC,NEWTON,ITOK,TOLINT,MXDM,D,IPIVOT,WCP,WP, + TANGNT,DIR,U,V,STO,STO2,PHI,PHIA,ISING,IRCODE) C CP print*,'IRCODE =',IRCODE,'wc=',(wc(i),i=1,n+1) print*,'=====================================================' C###################################################################### C *********** PROBLEM 4 (SCHWETLICK 4) N = 6 ************** PRINT*,'********* PROBLEM 4 (SCHWETLICK 4) N = 6 ********' N=6 C ENTER INITIAL SOLUTION POINT TO F(W) = 0. DO 429 L=1,N+1 WC(L) = 0.E0 429 CONTINUE C** PARAMETERS ******************** RHIGH = 1. RLOW = -20. HMIN = 1.E-8 HMAX = 100. EABS = 1.E-3 EREL = 1.E-3 ECOR = 1.E-5 EFIN = +1.0E-5 IFIN = -1 NEWTON = 7 ITOK = 7 NPRNT = -3 KSKIP = 0 KTARG = 52 TOLINT = EABS ORIENT = +1. NCOMP = N+1 VCOMP = 1.0 MXSTEP = 100 MPREF = 0 C MXDG = 3 NFAC = 0 DS = -8.E-2 C CALL ABCON(N,WC,FFIN,F2,DZF2,DRF2,NPRNT,KSKIP,ORIENT,NCOMP, + VCOMP,RHIGH,RLOW,KTARG,IFIN,EABS,EREL,ECOR,EFIN,DS,HMIN,HMAX, + MXSTEP,MXDG,MPREF,NFAC,NEWTON,ITOK,TOLINT,MXDM,D,IPIVOT,WCP,WP, + TANGNT,DIR,U,V,STO,STO2,PHI,PHIA,ISING,IRCODE) C CP print*,'IRCODE =',IRCODE,'wc=',(wc(i),i=1,n+1) print*,'=====================================================' C###################################################################### C *********** PROBLEM 5 (WATSON 1) N = 10 ************** PRINT*,'********** PROBLEM 5 (WATSON 1) N = 10 ************' N=10 C ENTER INITIAL SOLUTION POINT TO F(W) = 0. DO 529 L=1,N+1 WC(L) = 0.E0 529 CONTINUE C** PARAMETERS ******************** RHIGH = 1. RLOW = -20. HMIN = 1.E-6 HMAX = 100. EABS = 1.E-2 EREL = 1.E-2 ECOR = 1.E-4 EFIN = +1.0E-4 IFIN = -1 NEWTON = 7 ITOK = 7 NPRNT = -3 KSKIP = 0 KTARG = 1 TOLINT = EABS ORIENT = +1. NCOMP = N+1 VCOMP = 1.0 MXSTEP = 20 MPREF = 1 C MXDG = 3 NFAC = 0 DS = -8.E-2 C CALL ABCON(N,WC,FFIN,F4,DZF4,DRF4,NPRNT,KSKIP,ORIENT,NCOMP, + VCOMP,RHIGH,RLOW,KTARG,IFIN,EABS,EREL,ECOR,EFIN,DS,HMIN,HMAX, + MXSTEP,MXDG,MPREF,NFAC,NEWTON,ITOK,TOLINT,MXDM,D,IPIVOT,WCP,WP, + TANGNT,DIR,U,V,STO,STO2,PHI,PHIA,ISING,IRCODE) C CP print*,'IRCODE =',IRCODE,'wc=',(wc(i),i=1,n+1) print*,'=====================================================' C###################################################################### C *********** PROBLEM 6 N = 10 ************** PRINT*,'*********** PROBLEM 6 N = 10 **************' N=10 C ENTER INITIAL SOLUTION POINT TO F(W) = 0. DO 629 L=1,N+1 WC(L) = 0.E0 629 CONTINUE C** PARAMETERS ******************** RHIGH = 1. RLOW = -20. HMIN = 1.E-6 HMAX = 100. EABS = 1.E-2 EREL = 1.E-2 ECOR = 1.E-4 EFIN = +1.0E-4 IFIN = -1 NEWTON = 7 ITOK = 7 NPRNT = -3 KSKIP = 0 KTARG = 1 TOLINT = EABS ORIENT = +1. NCOMP = N+1 VCOMP = 1.0 MXSTEP = 20 MPREF = 1 C MXDG = 3 NFAC = 0 DS = -8.E-2 C CALL ABCON(N,WC,FFIN,F5,DZF5,DRF5,NPRNT,KSKIP,ORIENT,NCOMP, + VCOMP,RHIGH,RLOW,KTARG,IFIN,EABS,EREL,ECOR,EFIN,DS,HMIN,HMAX, + MXSTEP,MXDG,MPREF,NFAC,NEWTON,ITOK,TOLINT,MXDM,D,IPIVOT,WCP,WP, + TANGNT,DIR,U,V,STO,STO2,PHI,PHIA,ISING,IRCODE) C CP print*,'IRCODE =',IRCODE,'wc=',(wc(i),i=1,n+1) print*,'=====================================================' C###################################################################### C *********** PROBLEM 7 (BROWN FUNCTION) N = 10 ************** PRINT*,'********* PROBLEM 7 (BROWN FUNCTION) N = 10 ********' N=10 C ENTER INITIAL SOLUTION POINT TO F(W) = 0. DO 729 L=1,N+1 WC(L) = 0.E0 729 CONTINUE C** PARAMETERS ******************** RHIGH = 1. RLOW = -20. HMIN = 1.E-6 HMAX = 100. EABS = 1.E-2 EREL = 1.E-2 ECOR = 1.E-4 EFIN = +1.0E-4 IFIN = -1 NEWTON = 7 ITOK = 7 NPRNT = -3 KSKIP = 0 KTARG = 1 TOLINT = EABS ORIENT = +1. NCOMP = N+1 VCOMP = 1.0 MXSTEP = 20 MPREF = 1 C MXDG = 3 NFAC = 0 DS = -8.E-2 C CALL ABCON(N,WC,FFIN,F6,DZF6,DRF6,NPRNT,KSKIP,ORIENT,NCOMP, + VCOMP,RHIGH,RLOW,KTARG,IFIN,EABS,EREL,ECOR,EFIN,DS,HMIN,HMAX, + MXSTEP,MXDG,MPREF,NFAC,NEWTON,ITOK,TOLINT,MXDM,D,IPIVOT,WCP,WP, + TANGNT,DIR,U,V,STO,STO2,PHI,PHIA,ISING,IRCODE) C CP print*,'IRCODE =',IRCODE,'wc=',(wc(i),i=1,n+1) print*,'=====================================================' C###################################################################### C *********** PROBLEM 8 (BROWN FUNCTION) N = 25 ************** PRINT*,'********* PROBLEM 8 (BROWN FUNCTION) N = 25 ********' N=25 C ENTER INITIAL SOLUTION POINT TO F(W) = 0. DO 829 L=1,N+1 WC(L) = 0.E0 829 CONTINUE C** PARAMETERS ******************** RHIGH = 1. RLOW = -20. HMIN = 1.E-6 HMAX = 100. EABS = 1.E-2 EREL = 1.E-2 ECOR = 1.E-4 EFIN = +1.0E-4 IFIN = -1 NEWTON = 7 ITOK = 7 NPRNT = -3 KSKIP = 0 KTARG = 1 TOLINT = EABS ORIENT = +1. NCOMP = N+1 VCOMP = 1.0 MXSTEP = 20 MPREF = 1 C MXDG = 3 NFAC = 0 DS = -8.E-2 C CALL ABCON(N,WC,FFIN,F6,DZF6,DRF6,NPRNT,KSKIP,ORIENT,NCOMP, + VCOMP,RHIGH,RLOW,KTARG,IFIN,EABS,EREL,ECOR,EFIN,DS,HMIN,HMAX, + MXSTEP,MXDG,MPREF,NFAC,NEWTON,ITOK,TOLINT,MXDM,D,IPIVOT,WCP,WP, + TANGNT,DIR,U,V,STO,STO2,PHI,PHIA,ISING,IRCODE) C CP print*,'IRCODE =',IRCODE,'wc=',(wc(i),i=1,n+1) print*,'=====================================================' C ################################################################ C *********** PROBLEM 9 (BROWN FUNCTION) N = 50 ************** PRINT*,'********* PROBLEM 9 (BROWN FUNCTION) N = 50 ********' N=50 C ENTER INITIAL SOLUTION POINT TO F(W) = 0. DO 929 L=1,N+1 WC(L) = 0.E0 929 CONTINUE C** PARAMETERS ******************** RHIGH = 1. RLOW = -20. HMIN = 1.E-6 HMAX = 100. EABS = 1.E-2 EREL = 1.E-2 ECOR = 1.E-4 EFIN = +1.0E-4 IFIN = -1 NEWTON = 7 ITOK = 7 NPRNT = -3 KSKIP = 0 KTARG = 1 TOLINT = EABS ORIENT = +1. NCOMP = N+1 VCOMP = 1.0 MXSTEP = 20 MPREF = 1 C MXDG = 3 NFAC = 0 DS = -8.E-2 C CALL ABCON(N,WC,FFIN,F6,DZF6,DRF6,NPRNT,KSKIP,ORIENT,NCOMP, + VCOMP,RHIGH,RLOW,KTARG,IFIN,EABS,EREL,ECOR,EFIN,DS,HMIN,HMAX, + MXSTEP,MXDG,MPREF,NFAC,NEWTON,ITOK,TOLINT,MXDM,D,IPIVOT,WCP,WP, + TANGNT,DIR,U,V,STO,STO2,PHI,PHIA,ISING,IRCODE) C CP print*,'IRCODE =',IRCODE,'wc=',(wc(i),i=1,n+1) print*,'=====================================================' C################################################################### C *********** PROBLEM 10 F AND R, REGULARIZING HOMOTOPY ************** PRINT*,'PROBLEM 10 F AND R, REGULARIZING HOMOTOPY' N=2 C ENTER INITIAL SOLUTION POINT TO F(W) = 0. DO 1029 L=1,N+1 WC(L) = 0.E0 1029 CONTINUE WC(1) = 15. WC(2) = -2. WC(3) = 0. C** PARAMETERS ******************** RHIGH = 1. RLOW = -20. HMIN = 1.E-6 HMAX = 100. EABS = 1.E-2 EREL = 1.E-2 ECOR = 1.E-4 EFIN = +1.0E-4 IFIN = -1 NEWTON = 7 ITOK = 7 NPRNT = -3 KSKIP = 0 KTARG = 1 TOLINT = EABS ORIENT = +1. NCOMP = N+1 VCOMP = 1.0 MXSTEP = 100 MPREF = 1 C MXDG = 3 NFAC = 0 DS = -8.E-2 C CALL ABCON(N,WC,FFIN,F7,DZF7,DRF7,NPRNT,KSKIP,ORIENT,NCOMP, + VCOMP,RHIGH,RLOW,KTARG,IFIN,EABS,EREL,ECOR,EFIN,DS,HMIN,HMAX, + MXSTEP,MXDG,MPREF,NFAC,NEWTON,ITOK,TOLINT,MXDM,D,IPIVOT,WCP,WP, + TANGNT,DIR,U,V,STO,STO2,PHI,PHIA,ISING,IRCODE) C CP print*,'IRCODE =',IRCODE,'wc=',(wc(i),i=1,n+1) print*,'=====================================================' C###################################################################### C *********** PROBLEM 11 F AND R, STANDARD HOMOTOPY ************** PRINT*,'PROBLEM 11 F AND R, STANDARD HOMOTOPY' C do 999 l=1,5 C NPRNT = l-1 C print*,'nprnt =',nprnt N=2 C ENTER INITIAL SOLUTION POINT TO F(W) = 0. WC(1) = 15. WC(2) = -2. WC(3) = 0. C PARAMETERS RHIGH = 1. RLOW = -20. HMIN = 1.E-6 HMAX = 100. EABS = 1.E-2 EREL = 1.E-2 ECOR = 1.E-4 EFIN = +1.0E-4 IFIN = -1 NEWTON = 7 ITOK = 7 NPRNT = -3 KSKIP = 0 KTARG = 1 TOLINT = EABS ORIENT = +1. NCOMP = N+1 VCOMP = 1.0 MXSTEP = 100 MPREF = 0 MXDG = 3 NFAC = 0 DS = -8.E-2 C CALL ABCON(N,WC,FFIN,F3,DZF3,DRF3,NPRNT,KSKIP,ORIENT,NCOMP, + VCOMP,RHIGH,RLOW,KTARG,IFIN,EABS,EREL,ECOR,EFIN,DS,HMIN,HMAX, + MXSTEP,MXDG,MPREF,NFAC,NEWTON,ITOK,TOLINT,MXDM,D,IPIVOT,WCP,WP, + TANGNT,DIR,U,V,STO,STO2,PHI,PHIA,ISING,IRCODE) C CP print*,'IRCODE =',IRCODE,'wc=',(wc(i),i=1,n+1) print*,'=====================================================' 999 continue C OUTPUT RESULTS 1000 CONTINUE C END C ABCON BY B. N. LUNDBERG AND A. B. POORE C SUBROUTINE ABCON(N,WC,FFIN,F,DZF,DRF,NPRNT,KSKIP,ORIENT,NCOMP, + VCOMP,RHIGH,RLOW,KTARG,IFIN,EABS,EREL,ECOR,EFIN,DS,HMIN,HMAX, + MXSTEP,MXDG,MPREF,NFAC,NEWTON,ITOK,TOLINT,MXDM,D,IPIVOT,WCP,WP, + TANGNT,DIR,U,V,STO,STO2,PHI,PHIA,ISING,IRCODE) PARAMETER(MAXDEG=13,TOLIT=1.5,TOLF=1.5,NFCOR=15) C*********************************************************************** C *** A B C O N : PREDICTOR-CORRECTOR CONTINUATION USING ADAMS-BASHFORTH C *** VARIABLE ORDER PREDICTORS WITH ERROR-STEPSIZE CONTROL C *** (SINGLE PRECISION VERSION: UPDATED 5-27-93.) C ********************************************************************** C The authors of this software are Bruce N. Lundberg and Aubrey B. Poore C Permission to use, copy, modify, and distribute this software for any C purpose without fee is hereby granted, provided that this entire notice C is included in all copies of any software which is or includes a copy C or modification of this software and in all copies of the supporting C documentation for such software. C THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED C WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY C REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY C OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. C C E-MAIL: lundberg@math.colostate.edu C C 1 ******** P U R P O S E ****************************************** C C THIS PROGRAM COMPUTES A SEQUENCE OF SOLUTIONS (UK,RK)=W(NP1) C OF F(U,R) = 0 (WHERE F IS A FUNCTION FROM NP1=N+1 SPACE C TO N SPACE) STARTING WITH A GIVEN INITIAL SOLUTION C W = (UZERO,RZERO),AND TERMINATING AT RK = RLOW OR RHIGH, C FINAL VALUES OF THE PARAMETER R. C THE DETERMINANT OF THE SYSTEM DERIVATIVE IS MONITORED FOR C CHANGES IN SIGN TO DETECT SINGULAR POINTS. C C THIS IS AN IMPLEMENTATION OF THE ARCLENGTH CONTINUATION C ALGORITHM GIVEN IN THE PAPER : C C REFERENCE: "VARIABLE ORDER ADAMS-BASHFORTH PREDICTORS WITH C ERROR-STEPSIZE CONTROL FOR CONTINUATION METHODS", C BY BRUCE N. LUNDBERG AND AUBREY B. POORE, C SIAM J. OF SCI. AND STAT. COMPUTATION (12) C P. 695-723, MAY 1991. C C 2 ******* CALL ARGUMENTS FOR ABCON ********************************** C C ** INITIAL INPUT VARIABLES, FUNCTIONS, AND CONTROL PARAMETERS ** C N -- THE SIZ OF THE PROBLEM = NUMBER OF COMPONENTS IN F, THE C CONTINUATION PROCEEDS IN THE N+1 VARIABLES W = (Z, R) C WC(N+1) -- ON INPUT, CONTAINS INITIAL SOLUTION POINT TO F(W) = 0. C ON OUTPUT, CONTAINS LAST CORRECTED SOLUTION POINT. C FFIN(N) -- ON OUTPUT, CONTAINS F AT FINAL POINT WC C F, DZF, DRF -- USER WRITTEN SUBROUTINES EVALUATION THE PROBLEM C FUNCTIONS AND DERIVATIVES. NAMES MUST BE DECLARED EXTERNAL C IN THE CALLING PROGRAM. (SEE SECTION 6 BELOW). C NPRNT -- PRINT CONTROL CODE: C NPRNT=0, NO PRINTING OF INFORMATION C NPRNT=-3, PERFORMANCE COUNTS ONLY C NPRNT=-1, PRINT CORRECTED POINTS ONLY. C NPRNT=1, INITIAL HEADING WITH PARAMETER SETTINGS PRINTED C FIRST TWO STEPS PRINTED ALONG WITH SINGULAR POINT C INFO AND FINAL POINT INFO C NPRNT=2, ABOVE PLUS CORRECTED POINTS (FIRST TWO, EVERY KSKIP STEP) C NPRNT=3, ABOVE PLUS HEADING FOR STEPS AND OTHER STEP INFO. C NPRNT=4, ABOVE PLUS STATUS MESSAGES C NPRNT>4, ABOVE PLUS FINAL JACOBIAN CALCULATED AND PRINTED C K -- THE STEP NUMBER C DS -- THE SIZE OF THE PREDICTOR STEP TAKEN IN STEP K (K=1,2,..) C WC -- THE CORRECTED K TH TERM (UK,RK) OF THE GENERATED SEQUENCE C OF SOLUTIONS. C M+1 -- PREDICTOR USED IS AN A-B (M+1) STEP METHOD. C KNTROL -- THE NUMBER OF CORRECTOR ITERATIONS USED IN CORRECTION. C C KSKIP -- PRINT SKIPPING PARAMETER: PRINT EVERY KSKIP STEP. C ORIENT -- SETS ORIENTATION FOR INITIAL DIRECTION OF CONTINUATION C ORIENT = +1, THE PARAMETER WC(N+1) INITIALLY IS INCREASING C ORIENT = -1, THE PARAMETER INITIALLY DECREASES C NCOMP -- SETS THE TARGET COMPONENT FOR STOPPING THE CONTINUATION: C NCOMP = 1, ..., N+1, THE CONTINUATION WILL STOP AND LOCATE A POINT C WITH WC(NCOMP) = VCOMP ONCE WC(NCOMP) .GT. VCOMP C NCOMP = -1, ...,-(N+1) THE CONTINUATION WILL STOP AND LOCATE A C POINT WITH WC(-NCOMP) = VCOMP ONCE WC(-NCOMP) .LT. VCOMP C NCOMP = 0, THIS FEATURE IS DISABLED (SEE RLOW, RHIGH HOWEVER) C VCOMP -- TARGET VALUE FOR WC(ABS(NCOMP)) (SEE NCOMP) C RHIGH, RLOW -- THE BOUNDS OF CONTINUATION: IF THE PREDICTED OR C CORRECTED POINT CROSSES THE HYPERPLANES R = RHIGH OR C R = RLOW, THE CONTINUATION WILL STOP WITH A POINT WC C WITH WC(N+1) = RHIGH, (OR RLOW). C KTARG -- MINIMUM NUMBER OF STEPS BEFORE FINAL POINT STOPPING WILL BE C INVOKED. THIS CAN BE USED TO PICK OUT SUBSEQUENT PASSINGS C OF A PARAMETER VALUE FOR ATTENTION. C IFIN -- CONTROL PARAMETER FOR FINAL POINT CORRECTION C IFIN .GT. 0, FINAL CORRECTION DONE WITH A NEWTON OR CHORD C ITERATION WITH TARGET PARAMETER FIXED, NO BORDERING C AND .LE IFIN ITERATIONS. C IFIN .LT. 0, USES ROUTINE COREX BUT WITH DIR RESET TO GIVE C CORRECTION WITH TARGET PARAM. FIXED AND .LE. -IFIN ITERATIONS C IFIN .EQ 0, USES ROUTINE COREX; BORDERING ALGORITHM C WITH ITERATES IN A HYP. PLANE ORTH. TO USUAL PREDICTOR. C MAXIMUM ITERATES DETERMINED AS WITH POINTS (SEE NEWTON) C EABS, EREL -- RELATIVE AND ABSOLUTE PREDICTOR ERROR PARAMETERS C (DEFAULT VALUES COMPUTED FOR EACH ONE .LE. 0 ON ENTRY) C ECOR -- ERROR PARAMETER FOR CORRECTED POINTS C (DEFAULT VALUE COMPUTED IF .LE. 0 ON ENTRY) C EFIN -- ERROR PARAMETER FOR FINAL TARGET POINT: C EFIN .GT. 0, USE A CHORD CORRECTOR, STOPPING ERROR EFIN C EFIN .LT. 0, USE A NEWTON CORRECTOR, STOPPING ERROR -EFIN C EFIN .EQ. 0, USE NO CORRECTOR C DS -- INITIAL STEPSIZE FOR FIRST PREDICTOR STEP C (DEFAULT VALUE COMPUTED IF .LE. 0 ON ENTRY) C HMIN, HMAX -- MINIMUM AND MAXIMUM PREDICTOR STEPSIZES C MXSTEP -- MAXIMUM NUMBER OF CONTINUATION STEPS (POINTS) C MXDG -- THE MAXIMUM DEGREE OF THE INTERPOLATING POLYNOMIAL ALLOWED C IN THE PREDICTION PHASE (MXDG + 1 STEP METHOD). THIS MAY C BE SPECIFIED BY THE USER FROM 0 TO 12, OR C < 0, DEFAULT MXDG COMPUTED AS FOLLOWS: C = IMIN(12,IMAX(2,NINT(-ALOG10((EREL+EABS)/2.E0))) C THAT IS THE CLOSEST INTEGER TO THE NEGATIVE OF THE C COMMON LOG OF THE AVERAGE OF THE PREDICTOR ERRORS, C BUT AT MOST 12 AND AT LEAST 2. C MPREF -- AUTOMATIC PREDICTOR ORDER INCREASE UNTIL ORDER > MPREF C NFAC -- CONTROL PARAMETER FOR TANGENT COMPUTATION: C NFAC = 0 FOR COMPUTATION OF THE TANGENT AT WP (PREDICTED POINT) C NFAC = 1 FOR COMPUTATION OF THE TANGENT AT WC (CORRECTED POINT) C NEWTON -- CONTROL PARAMETER FOR CORRECTOR: C NEWTON = 1, 2, ... CORRECTOR ABORTS AFTER NEWTON CHORD ITERATES C NEWTON = -1, -2, ... CORRECTOR USES UP TO 'NEWTON' NEWTON ITERATES C (JACOBIAN REEVALUATED FOR EACH ITERATE) C ITOK -- MAX CORRECTOR ITERATES WHEN INCREASING NEXT STEPSIZES AFTER C CORR. FAILS(ITOK .GE.ABS(NEWTON) DISABLES THIS FEATURE). C TOLINT -- ERROR TOLERANCE FOR MAX NORM OF F(W) AT INITIAL POINT C IF THIS TOLERANCE IS NOT MET, IRCODE = -6 RETURNS. C C ** STORAGE AND WORK SPACE ARGUMENTS ** C MXDM -- LEADING DIMENSION OF CERTAIN TWO DOUBLE INDEXED WORK ARRAYS C C D(MXDM,N+1) -- ARRAY FOR STORING THE NXN JACOBIAN OF F WRT Z C IPIVOT(N+1) -- STORES PIVOT INFORMATION FOR LINEAR SYSTEM SOLVER C WCP(N+1) -- STORAGE FOR PENULTIMATE CORRECTED POINT C WP(N+1) -- PREDICTED SOLUTION POINT C TANGNT(N+1) -- UNIT TANGENT TO PATH C DIR(N+1) -- PREDICTION DIRECTION USED IN COREX AND PREDIC C U(N+1), V(N+1), STO(N+1), STO2(N+1) -- STORAGE VECTORS FOR COREX C PHI(MXDM,0:MAXDEG), PHIA(MXDM,0:MAXDEG) SCALED DIVIDED DIFF. STORAGE C C ISING -- INTEGER FLAG WHICH TURNS ON BIFURCAION POINT CHECK C = 0, NO CHECKING FOR BIFURCATION, ALL SINGULARITIES ASSUMED C TO BE FOLDS C > 0, ROUTINE SINGCK USED TO ANALYSE SINGULAR POINTS C ** ERROR RETURN CODE ** C IRCODE -- RETURNS INTEGER ERROR CODE: C = 0, NO ERRORS, NORMAL RETURN C = -1, MAXIMUM NUMBER OF STEPS COMPUTED WITHOUT REACHING TARGET C = -2, FAILURE IN FINAL PREDICTION C = -3, FAILURE IN FINAL CORRECTION C = -4, SINGULAR JACOBIAN OR ZERO PIVOT ENCOUNTER C = -5, INPUT ERROR, HMAX < HMIN C = -6, INPUT ERROR, INITIAL POINT DOES NOT SATISFY F=0 TO TOLINT C = -7, ZERO PIVOT IN INITIAL PLU C = K>0, CORRECTOR FAILED NFCOR TIMES AT STEP K C C 3 ********* INTERNAL FORTRAN PARAMETERS ****************************** C (THESE ARE SET INTERNALLY BUT MAY BE CHANGED TO SAVE STORAGE, ETC.) C MAXDEG -- THE MAXIMUM DEGREE OF INTERPOLATING POLYNOMIAL IN C THE PREDICTOR. AN M TH DEGREE POLY GIVES AN ADAMS- C BASHFORTH M+1 STEP METHOD. THE DEFAULT IS 13 C TOLIT, TOLF -- MINIMUM CONTRACTION PARAMETERS FOR COREX, SEE COREX C NFCOR -- THE MAXIMUM NUMBER OF CORRECTOR FAILURES ALLOWED C PER STEP; IF A CORRECTOR FAILS MORE THAN NFCOR C TIMES IN A STEP THE RUN IS ABORTED. C C 4 ********** LOCAL VARIABLES AND CONSTANTS *************************** C NP1 -- = N + 1 AND EQUALS THE DIMENSION OF (U,R) IN F(U,R) = 0. C M -- CURRENT DEGREE OF INTERP POLY(M+1 AND M+2 STEP METHODS) C C ** LOGICAL VARIABLES ** C ORFLG -- =, .TRUE. IF DEGREE REDUCTION IS INDICATED IN SELST2. C =, .FALSE. IF DEGREE REDUCTION IS NOT INDICATED. C MXFLG -- =, .TRUE. IF CONVERGENCE IS NOT OBTAINED IN SELST2. C =, .FALSE. IF CONVERGENCE IS OBTAINED. C C 5 ********** FUNCTIONS TO BE WRITTEN BY THE USER ********************* C THESE MUST BE CODED AS SUBROUTINES AND DECLARED EXTERNAL IN THE C CALLING PROGRAM AND HAVE THERE NAMES INCLUDED IN PLACE OF THE C ARGUMENTS F, DZF, DRF C C F -- THE FUNCTION TO BE CONTINUED FOR. (F(N,VALUE,ARGUMENT)) C DZF -- THE DERIVATIVE(NXN MATRIX) OF F WITH RESPECT TO THE FIRST C N VARIABLES.(DZF(N,MXDM,RESULT,ARGUMENT)) C DRF -- THE DERIVATIVE(N - VECTOR) OF F WITH RESPECT TO THE LAST C VARIABLE, THE PARAMETER R.(DRF(N,RESULT,ARGUMENT)) C C 6 ********** SUBROUTINES USED **************************************** C ** PREDICTOR SUBROUTINES ** C PREDIC -- COMPUTES PREDICTOR , UPDATES, SELECTS ORDER AND STEPSIZE. C PREDICTOR POLYNOMIAL. C PSDIR -- PREDICTS SOLUTION AND PRED. DIRECTION. C PRDINF -- COMPUTES PREDICTOR INFORMATION. C SELST2 --- SELECTS STEPSIZE TO SPECIFIED TOLERANCE. C ARCLG2 -- COMPUTES AN ESTIMATE OF ARCLENGTH FOR UPDATING DS. C PDNRM -- USED IN ARCLG2 TO ESTIMATE TANGENT NORM. C ** CORRECTOR SUBROUTINES ** C COREX -- TO CORRECT THE RESULTS OF PREDICTION, APPROX. SOLUTION TO C F = 0. ALSO COMPUTES NEW DW/DS AND SYSTEM DETERMINANT. C SGEFA,SGESL,SGEDI, -- SQUARE LINEAR SYSTEM SOLVER FROM C LINPACK.. C DWDS -- COMPUTES THE UNIT TANGENT TO THE PATH AT REGULAR POINTS. C FINALP -- CHECKS FOR PASSAGE OF A TARGET POINT. PREDICTS AND C CORRECTS (OR SETS UP FOR A FINAL CALL OF COREX) C FINAL TARGET POINT. C ** SINGULARITY CHECKING AND ANALYSIS ** C SINGCK -- CHECKS FOR THE PASSAGE OF FOLD AND BIFURCATION POINTS. C MAKMAT -- USED IN SINGCK TO CONSTRUCT SYSTEM MATRIX. C ** LINPACK ROUTINES USED ** C SDOT -- SINGLE PRECISION DOT PRODUCT. C SNRM2 -- SINGLE PRECISION EUCLIDEAN NORM. C SGEFA -- SINGLE PRECISION PLU DECOMPOSITION C SGESL -- SINGLE PRECISION SYSTEM SOLVER (USING PLU) C SGEDI -- SINGLE PRECISION DETERMINANT COMPUTATION C SAXPY,SASUM, SSWAP, SSCAL -- ROUTINES CALLED BY ABOVE C 7 ****************** TYPE AND DIMENSION *************************** REAL WC(N+1),WCP(N+1),WP(N+1),DIR(N+1),U(N+1),FFIN(N+1), + V(N+1),STO(N+1),STO2(N+1),D(MXDM,N+1),TANGNT(N+1), C ARRAYS FOR PREDICTOR AND SCALED DIVIDED DIFFERENCE STORAGE + G(0:MAXDEG,1:MAXDEG),PSYP(0:MAXDEG),PSY(0:MAXDEG), + A(0:MAXDEG),B(0:MAXDEG),PHI(MXDM,0:MAXDEG), + PSYARC(0:MAXDEG),PHIA(MXDM,0:MAXDEG), C REAL SPECIFICATION VARIABLES + EREL,EABS,ECOR,EFIN,HMIN,HMAX,RLOW,RHIGH,ORIENT, C REAL STORAGE VARIABLES + DS,DSP,DSM,DSMM1,DSMM2,DETOLD,DETNEW,TEMP C INTEGER IPIVOT(N+1),IFLAG,KNTROL,LEX,MXDG,MXSTEP,IRCODE, + NPRNT,M,N,NEWTON,NFAC,IFIN,NP1,MPREF,COUNT,ITOK,ISING C LOGICAL ORFLG,FAIL,UPDATE,SELECT,FINAL COMMON KJAC,KFUNC,KLIN EXTERNAL F,DZF,DRF C C 8 *********** EXECUTABLE CODE BEGINS HERE ************************** C C ****************** INITIALIZATION ********************************** C C *** INITIALIZING NP1,A,B,G *** C IRCODE = 0 NP1 = N + 1 LEX = 0 K = 1 A(0) = 1.0 B(0) = 1.0 G(0,1) = 1.0 C ** TESTING ERROR PARAMETERS AND/OR COMPUTING DEFAULTS ** EMACH = R1MACH(4) C ** TESTING AND/OR COMPUTING DEFAULT STEPSIZE FOR FIRST STEP ** IF(ECOR.LE.0.D0) ECOR = SQRT(EMACH) IF(EREL.LE.0.D0) EREL = EMACH**0.25D0 IF(EABS.LE.0.D0) EABS = EMACH**0.25D0 IF(ECOR.LE.1.E1*EMACH) ECOR = 1.E1*EMACH IF((EFIN.GT.0.D0).AND.(EFIN.LE.1.E1*EMACH)) EFIN = 1.E1*EMACH IF((EFIN.LT.0.D0).AND.(EFIN.GT.-1.E1*EMACH)) EFIN = -1.E1*EMACH IF(EREL.LE.1.E1*EMACH) EREL = 1.E1*EMACH EA = EABS ER = EREL EC = ECOR C ** COMPUTING THE DEFAULT INITIAL STEPSIZE DS IF DS.LE.0 ** IF(DS.LE.0.D0) THEN TEMP = SNRM2(NP1,WC,1) DS = 0.5*(SQRT(EABS) +SQRT(EREL)*TEMP ) END IF C ** COMPUTING THE DEFAULT MAXIMUM DEGREE (FOR PREDICTOR) IF(MXDG.LT.0) THEN MXDG = 1+NINT(-LOG10(.5E0*(EREL+EABS))) MXDG = MAX(MXDG,2) END IF MXDG = MIN(MXDG,12) C COUNTERS INITIALIZED KFUNC = 0 KJAC = 0 KLIN = 0 KCUT = 0 KSING = 0 MPREF = MAX(MPREF,1) FAIL = .FALSE. FINAL = .FALSE. C C *** ECHO PRINT OF SPECIFICATIONS ************************************ C C IF NPRNT .LT.1 SKIP PRINT IF(NPRNT.GT.0) THEN PRINT*,'********************************************************' PRINT*,'****** ARCLENGTH CONTINUATION ALONG A SOLUTION *********' PRINT*,'****** CURVE OF F(U,R) = 0. *********' PRINT*,'********************************************************' PRINT*,'CONTINUATION IN THE INTERVAL (',RLOW,' , ',RHIGH,')' PRINT*,'WITH INITIAL ORIENTATION: ',ORIENT PRINT*,'AND MAXIMUM NUMBER OF STEPS: ',MXSTEP PRINT*,' PRED STEPSIZE BOUNDS: MIN = ',HMIN,' MAX = ',HMAX PRINT*,' INITIAL STEPSIZE: ',DS PRINT*,' ERROR REQUIREMENTS FOR PREDICTOR: REL,ABS',EREL,EABS PRINT*,'THE MAX STEP METH. OF PRED: ',MXDG+1,'MPREF=',MPREF PRINT*,'PREFERED MAX. # CORR ITERATIONS ITOK= ',ITOK PRINT*,'ERROR REQ. FOR CORR.: STEPS',ECOR,'FINAL POINT',ABS(EFIN) IF(NEWTON.LT.0) THEN PRINT*,'CORRECTOR USES ** NEWTONS METHOD **' ELSE PRINT*,'CORRECTOR USES ** CHORD METHOD **' END IF IF(NFAC.EQ.1) THEN PRINT*,'WITH DW/DS EVALUATED AT THE ** CORRECTED ** PT.' ELSE PRINT*,'WITH DW/DS EVALUATED AT THE ** PREDICTED ** PT.' END IF IF(NCOMP.NE.0) THEN PRINT*,'TARGET POINT: POINT WITH',ABS(NCOMP),'COMPONENT',VCOMP END IF PRINT*,'FINAL(TARGET) POINT COMPUTED WITH:' IF(NCOMP.GT.0) THEN IF(EFIN.LT.0.) PRINT*,'NEWTON ITERATION' IF(EFIN.GT.0.) PRINT*,'CHORD ITERATION' PRINT*,'WITH THE TARGET PARAMETER HELD FIXED' END IF IF(NCOMP.LT.0) THEN IF(EFIN.LT.0.) PRINT*,'NEWTON ITERATION' IF(EFIN.GT.0.) PRINT*,'CHORD ITERATION' PRINT*,'WITH BORDERING ALGORITHM AND TARGET PARAMETER FIXED' END IF IF(NCOMP.EQ.0) THEN IF(EFIN.LT.0.) PRINT*,'NEWTON ITERATION' IF(EFIN.GT.0.) PRINT*,'CHORD ITERATION' PRINT*,'WITH ITERATES IN A HYPERPLANE ORTH. PRED. DIRECTION' END IF IF(EFIN.EQ.0) THEN PRINT*,'WITH NO CORRECTION (PREDICTION RETURNED)' END IF PRINT*,'############### START OUTPUT #####################' END IF C C **************** TEST DATA ************************************ C IF(HMIN.GT.HMAX) THEN IF(NPRNT.GT.3) THEN PRINT*,' ** DATA ERROR: STEPSIZE BOUNDS **' END IF IRCODE = -5 RETURN END IF CALL F(N,STO,WC) TEMP = 0.0 DO 25 L=1,N TEMP = MAX(TEMP,ABS(STO(L))) 25 CONTINUE IF(TEMP.GT.TOLINT) THEN IF(NPRNT.GT.1) THEN PRINT*,'DATA ERROR: INITIAL POINT NOT A ZERO OF F(X,R)' PRINT*,'CHANGE ERROR TOLERANCE IN CALLING PROGRAM' PRINT*,'OR OBTAIN MORE ACCURACY IN INITIAL POINT' PRINT*,'F : ',(STO(L),L=1,N) PRINT*,'SUP-NORM(F): ',TEMP END IF IRCODE = -6 RETURN END IF C C ***************** INITIAL STEP ************************************* C M = 0 IFAILF = 0 ARCLTH = 0. ARCP = 0. ARCS = 0. CALL DZF(N,MXDM,D,WC) CALL DRF(N,STO,WC) CALL SGEFA(D,MXDM,N,IPIVOT,IFLAG) IF(IFLAG.GT.0) THEN IF(NPRNT.GE.3) THEN PRINT*,'ZERO PIVOT IN INITIAL PLU' END IF IRCODE = -7 RETURN END IF CALL SGEDI(D,MXDM,N,IPIVOT,DETOLD,V,10) CALL DWDS(D,N,MXDM,STO,TANGNT,IPIVOT,ORIENT) 28 DO 30 L=1,NP1 PHI(L,0) = TANGNT(L) WCP(L) = WC(L) 30 CONTINUE C C** FIRST PREDICTION-CORRECTION ** C COUNT = 0 CALL PSDIR(M,N,MAXDEG,MXDM,DS,WCP,WP,DIR,G,B,PHI) C C CHECK BOUNDS OF CONTINUATION AND REDUCE STEPSIZE SO FIRST STEP C IS NOT THE LAST STEP C IF(WP(NP1).LE.RLOW) DS =0.5*(RLOW-WC(NP1))/TANGNT(NP1) IF(WP(NP1).GE.RHIGH) DS =0.5*(RHIGH-WC(NP1))/TANGNT(NP1) IF(NCOMP.GT.0) THEN IF((WP(NCOMP).GE.VCOMP).AND.(WC(NCOMP).LE.VCOMP)) THEN DS=0.5*(VCOMP-WC(NCOMP))/TANGNT(NCOMP) END IF ELSE IF(NCOMP.LT.0) THEN IF((WP(-NCOMP).LE.VCOMP).AND.(WC(-NCOMP).GE.VCOMP)) THEN DS=0.5*(VCOMP-WC(-NCOMP))/TANGNT(-NCOMP) END IF END IF C 32 CALL PSDIR(M,N,MAXDEG,MXDM,DS,WCP,WP,DIR,G,B,PHI) C KNTROL = NEWTON CALL COREX(N,MXDM,WC,WP,DIR,TANGNT,D,U,V,STO,STO2,EC,ORIENT, + DETNEW,KNTROL,NFAC,IPIVOT,ENRM,TOLIT,TOLF,F,DZF,DRF,NPRNT) C C ** HALVING STEPSIZE IF CORRECTOR FAILS -- TRY NFCOR TIMES ** IF(KNTROL.EQ.0) THEN IF(NPRNT.GT.3) PRINT*,'** CORRECTOR FAILURE **' IF(COUNT.LT.NFCOR) THEN COUNT = COUNT + 1 DS = 0.5*DS GO TO 32 ELSE IRCODE = 1 RETURN END IF END IF C C ** INITIAL OUTPUT ************************************************** C IF(ABS(NPRNT).EQ.1) THEN PRINT*,(WCP(L),L=1,NP1) PRINT*,(WC(L),L=1,NP1) ELSE IF(NPRNT.GT.2) THEN PRINT*,'THE INITIAL POINT: (U,R): ',(WCP(L),L=1,NP1) PRINT*,'** STEP 1 **: (U,R): ',(WC(L),L=1,NP1) 51 PRINT*,'TANGENT: ',(TANGNT(L),L=1,NP1) PRINT*,'SUPNORM(WC-WP):',ENRM END IF C 55 CONTINUE C PSYP(0) = DS PSY(0) = DS KPRNT = KSKIP C C********************************************************************* C ** MAIN LOOP: PREDICTOR-CORRECTOR STEPS ** C********************************************************************* C DO 1000 K=2,MXSTEP C ** SHIFTING AND UPDATING ** DETOLD = DETNEW DSP = DS ARCP = ARCP + DS UPDATE = .TRUE. SELECT = .TRUE. COUNT = 0 C C ******* P R E D I C T I O N *************************************** C IF(FAIL) THEN DS = DSP SELECT = .FALSE. ELSE SELECT = .TRUE. EA = EABS ER = EREL END IF C 200 CALL PREDIC(M,N,MAXDEG,MXDM,WCP,WC,TANGNT,DSP,DS, + WP,DIR,PSY,PSYP,PSYARC,PHI,PHIA,G,A,B,STO,EA,ER, + HMIN,HMAX,SELECT,UPDATE,LEX,MPREF,MXDG,K) C IF(NPRNT.GT.1) THEN PRINT*,'PRED. STP: ',DSP,' CORR. ITER',KNTROL PRINT*,'************** ',K,' ****************************' END IF KNTROL = NEWTON C ** UPDATE ARCLENGTH ACCUMULATION ** IF(UPDATE) ARCLTH = ARCLTH + DSP C C????????? B O U N D A N D T A R G E T C H E C K I N G ?????????? C C ** CHECKING IF PREDICTOR IS OUTSIDE OF THE BOUNDS OF CONTINUATION ** C SET UP (OR PREFORM) FINAL CORRECTION AS INDICATED (FINAL = .TRUE.) IF(K.GE.KTARG) THEN CALL FINALP(N,M,MXDM,MAXDEG,D,WCP,WP,V,U,DS,DIR,EFIN, + RLOW,RHIGH,LEX,KNTROL,G,PSYP,PSY,A,B,PHI,IPIVOT,WC, + NCOMP,VCOMP,F,DZF,DRF,NPRNT,FINAL,IFIN,EC,IRCODE) IF(FINAL) THEN IF((KNTROL.EQ.0).AND.(EFIN.NE.0.)) GO TO 225 IF(EFIN.EQ.0.) THEN IF(NPRNT.GT.3) PRINT*,'NO FINAL CORRECTION' KNTROL = 0 GO TO 1890 END IF IF((IFIN.GT.0).AND.(KNTROL.GT.0)) THEN IF(NPRNT.GT.3) PRINT*,'CORRECTION IN FINALP' GO TO 889 END IF IF((IFIN.LT.0).AND.(NPRNT.GT.3)) THEN PRINT*,'VERTICAL CORRECTION IN COREX' END IF IF((IFIN.EQ.0).AND.(NPRNT.GT.3)) THEN PRINT*,'HYPERPLANE CORRECTION, FINAL' END IF END IF END IF C C ******* C O R R E C T I O N ************************************* C CALL COREX(N,MXDM,WC,WP,DIR,TANGNT,D,U,V,STO,STO2,EC,ORIENT, + DETNEW,KNTROL,NFAC,IPIVOT,ENRM,TOLIT,TOLF,F,DZF,DRF,NPRNT) C C ??????? CORRECTOR FAILURE ?????????????????????? C C ** THE FOLLOWING CODE HANDLES CORRECTOR FAILURES C REDUCING STEPSIZE IF CORRECTOR FAILS -- TRY NFCOR TIMES ** C 225 CONTINUE IF(KNTROL.EQ.0) THEN FINAL = .FALSE. EC = ECOR IF(NPRNT.GT.3) PRINT*,' ** CORRECTOR FAILURE **' KCUT = KCUT + 1 SELECT = .FALSE. UPDATE = .FALSE. IF(COUNT.LT.NFCOR) THEN COUNT = COUNT + 1 IF((COUNT.EQ.1).AND.(.NOT.FAIL)) THEN M = 0 EA = .01*EABS ER = .01*EREL CALL SELST2(M,N,MAXDEG,MXDM,DS,DSP,EA,ER,HMIN,HMAX, + WCP,WP,G,PSYP,PSY,PHI,B,A,DIR,ORFLG) DS = MIN(DS,DSP) ELSE DS = 0.5*DS END IF GO TO 200 ELSE IRCODE = K RETURN END IF END IF C ** DECIDING IF NEXT STEP SHOULD BE CAUTIOUS RESET COUNT ** IF(KNTROL.LE.ITOK) FAIL = .FALSE. IF(COUNT.GT.0) FAIL = .TRUE. C C ??????????? SINGULARITY DETECTION ???????????????????????????? C 255 IF(DETNEW*DETOLD.LE.0.0) THEN IF(ISING.GT.0) THEN CALL SINGCK(N,MXDM,WC,WCP,D,PHI(1,0),PHI(1,0),V,DZF,DRF, + INDIC,IPIVOT) ELSE INDIC = 1 END IF IF(INDIC.GT.0) THEN KSING = KSING + 1 IF(NPRNT.GT.1) THEN PRINT*,'ARCLENGTH FROM LAST SING: ',ARCLTH-ARCS END IF ARCS = ARCLTH END IF IF(INDIC.EQ.1) THEN IF(NPRNT.GT.0) THEN PRINT*,'*** A FOLD POINT HAS BEEN PASSED ***' END IF C ** CHANGE ORIENTATION AND REVERSE TANGENT ** ORIENT = -ORIENT DO 275 L=1,NP1 TANGNT(L) = -TANGNT(L) 275 CONTINUE END IF END IF IF((NPRNT.GT.0).and.(INDIC.EQ.2)) THEN PRINT*,'*** A BIFURCATION POINT HAS BEEN PASSED ***' END IF INDIC = 0 C C **** PRINTED OUTPUT *********************************************** C 889 CONTINUE IF(FINAL.AND.(NPRNT.GT.1)) THEN PRINT*,'******* FINAL CORRECTION ********' KPRNT = KSKIP END IF IF((KPRNT.LT.KSKIP).OR.(KSKIP.LT.0)) THEN KPRNT = KPRNT+1 GO TO 900 ELSE KPRNT = 0 END IF C SKIP THIS PRINTING IF KPRNT 0, NEWTON'S METHOD CORRECTOR (NO BORDERING ALGORITHM) C WITH ICOMP COMPONENT KEPT FIXED AT VTARG DO 235 L=1,N+1 WC(L) = WP(L) 235 CONTINUE WC(ICOMP) = VTARG ITMAX = MAX(ABS(KNTROL),IFIN) DO 300 I=1,ITMAX JFLAG = 0 IF((EFIN.LT.0.).OR.(I.EQ.1)) THEN C NEWTONS METHOD OPTION AND FIRST STEP CALL DZF(N,MXDM,D,WC) IF(ICOMP.LT.NP1) THEN DO 240 J=ICOMP,N-1 DO 250 L=1,N D(L,J) = D(L,J+1) 250 CONTINUE 240 CONTINUE CALL DRF(N,STO1,WC) DO 270 L=1,N D(L,N) = STO1(L) 270 CONTINUE END IF IFLAG = 1 CALL SGEFA(D,MXDM,N,IPIVOT,IFLAG) IF(IFLAG.GT.0) THEN IF(NPRNT.GT.3) THEN PRINT*,'ZERO PIVOT IN PLU IN FINAL********' END IF IRCODE = -4 KNTROL = 0 DO 275 L=1,NP1 WC(L) = WCP(L) 275 CONTINUE RETURN END IF END IF CALL F(N,STO2,WC) CALL SGESL(D,MXDM,N,IPIVOT,STO2,0) C UPDATE ALL BUT ICOMP ENTRY (WC(ICOMP) FIXED AT VTARG) DO 280 L=1,ICOMP-1 WC(L) = WC(L) - STO2(L) DIFNRM = ABS(STO2(L)) ZNRM = (ABS(WC(L)) + ABERR)*ABS(EFIN) IF(DIFNRM.GT.ZNRM) JFLAG = 1 280 CONTINUE DO 290 L=ICOMP,N WC(L+1) = WC(L+1) - STO2(L) DIFNRM = ABS(STO2(L)) ZNRM = (ABS(WC(L+1)) + ABERR)*ABS(EFIN) IF(DIFNRM.GT.ZNRM) JFLAG = 1 290 CONTINUE C TESTING FOR CONVERGENCE IF(JFLAG.EQ.0) GO TO 310 300 CONTINUE C INDICATING FAILURE TO CONVERGE KNTROL = 0 IRCODE = -3 DO 305 L=1,NP1 WC(L) = WCP(L) 305 CONTINUE RETURN 310 CONTINUE C INDICATING NUMBER OF ITERATIONS IN FINAL CORRECTION KNTROL = I RETURN END C C C C C*********** FUNCTIONS AND DERIVATIVES ******************* C SUBROUTINE F1(N,RESULT,ARG) REAL RESULT(N),ARG(N+1) INTEGER N COMMON KJAC,KFUNC C *** THIS SUBROUTINE COMPUTES F(X,R) WHERE F IS THE FUNCTION C WE WISH TO CONTINUE ALONG. ARG CONTAINS THE VALUE OF (X(N),R) C AT WHICH F IS COMPUTED. RESULT(N) CONTAINS F(X,R) ON RETURN. C ARG,N ARE UNCHANGED. C C C ** SCHWETLICK EXAMPLE #3 ************ C C KFUNC = KFUNC + 1 RESULT(1) = 200.*ARG(1)**3-200.*ARG(1)*ARG(2)+ARG(1)-1. RESULT(2) = -100.*ARG(1)**2+110.1*ARG(2)+9.9*ARG(4)-20. RESULT(3) = 180.*ARG(3)**3-180.*ARG(3)*ARG(4)+ARG(3)-1. RESULT(4) = 9.9*ARG(2)-90.*ARG(3)**2+100.1*ARG(4)-20. C CONSTRUCTING STANDARD HOMOTOPY RESULT(1) = RESULT(1) - (1.-ARG(5))*(-6004.) RESULT(2) = RESULT(2) - (1.-ARG(5))*(-1040.) RESULT(3) = RESULT(3) - (1.-ARG(5))*(-5404.) RESULT(4) = RESULT(4) - (1.-ARG(5))*(-940.) RETURN END C C C SUBROUTINE DZF1(N,MXDM,DZ,ARG) REAL DZ(MXDM,N), ARG(N+1) INTEGER N,MXDM COMMON KJAC,KFUNC C *** COMPUTES THE DERIVATIVE(NXN MATRIX) OF THE FUNCTION F C WITH RESPECT TO U AT THE VECTOR (U(N),R), WHICH IS PASSED C IN ARG(N+1). THE DERIVATIVE IS RETURNED IN DZ(N,N). C C ********* SCHWETLICK EXAMPLE # 3 ************** C KJAC = KJAC+1 C DO 20 J=1,N DO 10 I=1,N DZ(I,J) = 0. 10 CONTINUE 20 CONTINUE C DZ(1,1) =600.*ARG(1)**2-200.*ARG(2)+1. DZ(1,2) = -200.*ARG(1) C 2ND COMPONENT DZ(2,1) = -200.*ARG(1) DZ(2,2) = 110.1 DZ(2,4) = 9.9 C THIRD COMPONENT DZ(3,3) = 540.*ARG(3)**2-180.*ARG(4)+1. DZ(3,4) = -180.*ARG(3) C FOURTH COMPONENT DZ(4,2) = 9.9 DZ(4,3) = -180.*ARG(3) DZ(4,4) = 100.1 RETURN END C C C SUBROUTINE DRF1(N,DR,ARG) REAL DR(N),ARG(N+1) INTEGER N C *** COMPUTES THE DERIVATIVE(N-VECTOR) OF THE FUNCTION F C WITH RESPECT TO R AT (U,R) WHICH IS CONTAINED IN ARG(N+1) C THE DERIVATIVE IS RETURNED IN DR(N). C DR(1) = -6004. DR(2) = -1040. DR(3) = -5404. DR(4) = - 940. RETURN END C FUNCTION SCHWETLICK # 4 ******************************************* C C **** FUNCTION AND DERIVATIVES **** C C SUBROUTINE F2(N,RESULT,ARG) REAL RESULT(N),ARG(N+1) INTEGER N COMMON KJAC,KFUNC C *** THIS SUBROUTINE COMPUTES F(X,R) WHERE F IS THE FUNCTION C WE WISH TO CONTINUE ALONG. ARG CONTAINS THE VALUE OF (X(N),R) C AT WHICH F IS COMPUTED. RESULT(N) CONTAINS F(X,R) ON RETURN. C ARG,N ARE UNCHANGED. C C C ** SCHWETLICK EXAMPLE #4 ************ C FI(Y) = 5.6E-8*(EXP(25.*Y) - 1.) FU(Y) = 7.65*ATAN(1962.*Y) C KFUNC = KFUNC + 1 RESULT(1) = 1.E-4*(ARG(1) - ARG(3))+(ARG(1)-ARG(2))/39.+ 1 (ARG(1)+ARG(7))/51. RESULT(2) = .1*(ARG(2)-ARG(6))+(ARG(2)-ARG(1))/39.+FI(ARG(2)) RESULT(3) = 1.E-4*(ARG(3)-ARG(1))+(ARG(3)-ARG(4))/25.5 RESULT(4) = (ARG(4)-ARG(3))/25.5+ARG(4)/0.62+ARG(4)-ARG(5) RESULT(5) = (ARG(5)-ARG(6))/13.+ARG(5)-ARG(4)+FI(ARG(5)) RESULT(6) = (ARG(6)-ARG(2))*.1+(ARG(6)-ARG(5))/13.+ARG(6) 1 -FU(ARG(3)-ARG(1))/0.201 RETURN END C C C SUBROUTINE DZF2(N,MXDM,DZ,ARG) REAL DZ(MXDM,N), ARG(N+1),FN(1001) INTEGER N,MXDM COMMON KJAC,KFUNC C *** COMPUTES THE DERIVATIVE(NXN MATRIX) OF THE FUNCTION F C WITH RESPECT TO U AT THE VECTOR (U(N),R), WHICH IS PASSED C IN ARG(N+1). THE DERIVATIVE IS RETURNED IN DZ(N,N). C C ********* SCHWETLICK EXAMPLE # 4 ************** C DI(Y) = 25.*5.6E-8*EXP(25.*Y) DU(Y) = 1962.*7.65/(1.+(1962.*Y)**2) KJAC = KJAC+1 C DO 20 J=1,N DO 10 I=1,N DZ(I,J) = 0. 10 CONTINUE 20 CONTINUE C DZ(1,1) =1.E-4+1./39.+1./51. DZ(1,2) = -1./39. DZ(1,3) = -1.E-4 C 2ND COMPONENT DZ(2,1) = -1./39. DZ(2,2) = .1+1./39. + DI(ARG(2)) DZ(2,6) = -.1 C THIRD COMPONENT DZ(3,1) = -1.E-4 DZ(3,3) = 1.E-4+1./25.5 DZ(3,4) = -1./25.5 C FOURTH COMPONENT DZ(4,3) = -1./25.5 DZ(4,4) = 1./25.5+1./0.62+1. DZ(4,5) = -1. C FIFTH COMPONENT DZ(5,4) = -1. DZ(5,5) = 1./13.+1.+DI(ARG(5)) DZ(5,6) = -1./13. C SIXTH COMPONENT DZ(6,1) = DU(ARG(3) - ARG(1))/0.201 DZ(6,2) = -.1 DZ(6,3) = -DU(ARG(3)-ARG(1))/0.201 DZ(6,5) = -1./13. DZ(6,6) = .1+1./13.+1. RETURN END C C C SUBROUTINE DRF2(N,DR,ARG) REAL DR(N),ARG(N+1) INTEGER N C *** COMPUTES THE DERIVATIVE(N-VECTOR) OF THE FUNCTION F C WITH RESPECT TO R AT (U,R) WHICH IS CONTAINED IN ARG(N+1) C THE DERIVATIVE IS RETURNED IN DR(N). C DO 10 I=1,N DR(I) = 0.0 10 CONTINUE DR(1) = 1./51. RETURN END C C C************ F & R FUNCTION WITH STANDARD HOMOTOPY C C **** FUNCTION AND DERIVATIVES **** C SUBROUTINE F3(N,RESULT,ARG) C ** FREUDENSTIEN AND ROTH ** REAL RESULT(N),ARG(N+1) INTEGER N COMMON KJAC,KFUNC KFUNC = KFUNC+1 C *** THIS SUBROUTINE COMPUTES F(X,R) WHERE F IS THE FUNCTION C WE WISH TO CONTINUE ALONG. ARG CONTAINS THE VALUE OF (X(N),R) C AT WHICH F IS COMPUTED. RESULT(N) CONTAINS F(X,R) ON RETURN. C ARG,N ARE UNCHANGED. C C RESULT(1)=ARG(1)+5.*ARG(2)**2-ARG(2)**3-2.*ARG(2) 1 -47.+ 34.*ARG(3) RESULT(2)=ARG(1)+ARG(2)**2+ARG(2)**3-14.*ARG(2) 1 -39. +10.*ARG(3) RETURN END C C SUBROUTINE DZF3(N,MXDM,DZ,ARG) REAL DZ(MXDM,N), ARG(N+1) INTEGER N,MXDM COMMON KJAC,KFUNC KJAC = KJAC + 1 C *** COMPUTES THE DERIVATIVE(NXN MATRIX) OF THE FUNCTION F C WITH RESPECT TO U AT THE VECTOR (U(N),R), WHICH IS PASSED C IN ARG(N+1). THE DERIVATIVE IS RETURNED IN DZ(N,N). C DZ(1,1) = 1.0 DZ(2,1) = 1.0 DZ(1,2) = -3.*ARG(2)**2+10.*ARG(2)-2. DZ(2,2) = 3.*ARG(2)**2+2.*ARG(2)-14. RETURN END C C C SUBROUTINE DRF3(N,DR,ARG) REAL DR(N),ARG(N+1) INTEGER N C *** COMPUTES THE DERIVATIVE(N-VECTOR) OF THE FUNCTION F C WITH RESPECT TO R AT (U,R) WHICH IS CONTAINED IN ARG(N+1) C THE DERIVATIVE IS RETURNED IN DR(N). C DR(1) = 34. DR(2) = 10. RETURN END C ***** FUNCTION 1 FROM WATSON ***************** C C SUBROUTINE F4(NV,RESULT,ARG) REAL RESULT(NV),ARG(NV+1) INTEGER NV COMMON KJAC,KFUNC C ** WATSON 1 ** C C KFUNC = KFUNC + 1 SUM = 0. DO 10 L=1,NV SUM = SUM + ARG(L)**3 10 CONTINUE TEMP = 1./(2.*FLOAT(NV)) DO 20 L=1,NV RESULT(L) = (SUM + FLOAT(L))*TEMP 20 CONTINUE C FORMING HOMOTOPY (REG.) TEMP = ARG(NV+1) DO 30 L=1,NV RESULT(L) = TEMP*(ARG(L)-RESULT(L))+ (1.-TEMP)*ARG(L) 30 CONTINUE RETURN END C C SUBROUTINE DZF4(NV,MXDM,DZ,ARG) REAL DZ(MXDM,NV), ARG(NV+1),FN(101) INTEGER N,NV COMMON KJAC,KFUNC KJAC = KJAC+1 SUM = 0. N = NV TEMP = 1./(2.*FLOAT(N)) TEMP =-3.*TEMP*ARG(N+1) DO 40 I=1,N TEMP1 = TEMP*ARG(I)**2 DO 30 L=1,N DZ(L,I) = TEMP1 30 CONTINUE DZ(I,I) = 1. + DZ(I,I) 40 CONTINUE RETURN END SUBROUTINE DRF4(N,DR,ARG) REAL DR(N),ARG(N+1) INTEGER N C *** COMPUTES THE DERIVATIVE(N-VECTOR) OF THE FUNCTION F C WITH RESPECT TO R AT (U,R) WHICH IS CONTAINED IN ARG(N+1) C THE DERIVATIVE IS RETURNED IN DR(N). C TEMP = 1./(2.*FLOAT(N)) SUM = 0. DO 10 L=1,N SUM = SUM + ARG(L)**3 10 CONTINUE DO 20 L=1,N DR(L) =-TEMP*(SUM + FLOAT(L)) 20 CONTINUE RETURN END C******************* WAT 2 *********************************** C SUBROUTINE F5(N,V,X) REAL V(N),X(N+1) INTEGER N COMMON KJAC,KFUNC C *** THIS SUBROUTINE COMPUTES F(X,R) WHERE F IS THE FUNCTION C WE WISH TO CONTINUE ALONG. X CONTAINS THE VALUE OF (X(N),R) C AT WHICH F IS COMPUTED. V(N) CONTAINS F(X,R) ON RETURN. C X,N ARE UNCHANGED. C C******************************************************************** C C SUBROUTINE F(X,V) -- EVALUATES WAT 2 FUNCTION AT THE POINT C X, AND RETURNS THE VALUE IN V. C C******************************************************************** C KFUNC = KFUNC + 1 V(1) = 0.01*(X(1)+X(2)+1.)**3 DO 10 L=2,N-1 V(L) = 0.01*(X(L-1)+X(L)+X(L+1)+1.)**3 10 CONTINUE V(N) = 0.01*(X(N-1)+X(N)+1.)**3 XLMDA = X(N+1) DO 20 L=1,N V(L) = X(L) - XLMDA*V(L) 20 CONTINUE RETURN END C* SUBROUTINE DZF5(N,MXDM,DZ,X) REAL DZ(MXDM,N), X(N+1) INTEGER N,MXDM COMMON KJAC,KFUNC C *** COMPUTES THE DERIVATIVE(NXN MATRIX) OF THE FUNCTION F C WITH RESPECT TO U AT THE VECTOR (U(N),R), WHICH IS PASSED C IN X(N+1). THE DERIVATIVE IS RETURNED IN DZ(N,N). C C KJAC = KJAC+1 C C******************************************************************** C C THE JACOBIAN MATRIX FOR WAT 2 FUNCTION EVALUATED AT C C******************************************************************** C DO 20 I=1,N DO 10 L=1,N DZ(L,I)= 0.0 10 CONTINUE 20 CONTINUE DZ(1,1)=0.03*(X(1)+X(2)+1.)**2 DZ(1,2)=0.03*(X(1)+X(2)+1.)**2 DZ(N,N-1)=0.03*(X(N-1)+X(N)+1.)**2 DZ(N,N)=0.03*(X(N-1)+X(N)+1.)**2 DO 30 K=2,N-1 TEMP=0.03*(X(K-1)+X(K+1)+X(K)+1.)**2 DZ(K,K-1) = TEMP DZ(K,K)=TEMP DZ(K+1,K)=TEMP 30 CONTINUE XLMDA = X(N+1) DO 40 L=1,N DO 35 J=1,N DZ(L,J) = -XLMDA*DZ(L,J) 35 CONTINUE DZ(L,L) = 1.+DZ(L,L) 40 CONTINUE RETURN END C SUBROUTINE DRF5(N,DR,X) REAL DR(N),X(N+1) INTEGER N C *** COMPUTES THE DERIVATIVE(N-VECTOR) OF THE FUNCTION F C WITH RESPECT TO R AT (U,R) WHICH IS CONTAINED IN X(N+1) C THE DERIVATIVE IS RETURNED IN DR(N). DR(1) = -0.01*(X(1)+X(2)+1.)**3 DO 10 L=2,N-1 DR(L) = -0.01*(X(L-1)+X(L)+X(L+1)+1.)**3 10 CONTINUE DR(N) = -0.01*(X(N-1)+X(N)+1.)**3 RETURN END C ***************** BR ********************** SUBROUTINE F6(N,V,X) REAL V(N),X(N+1) INTEGER N COMMON KJAC,KFUNC C *** THIS SUBROUTINE COMPUTES F(X,R) WHERE F IS THE FUNCTION C WE WISH TO CONTINUE ALONG. X CONTAINS THE VALUE OF (X(N),R) C AT WHICH F IS COMPUTED. V(N) CONTAINS F(X,R) ON RETURN. C X,N ARE UNCHANGED. C C******************************************************************** C C SUBROUTINE F(X,V) -- EVALUATES BROWN'S FUNCTION AT THE POINT C X, AND RETURNS THE VALUE IN V. C C******************************************************************** C KFUNC = KFUNC + 1 PROD=1.0E0 DO 10 J=1,N 10 PROD=PROD*X(J) V(1)=PROD-1.0E0 SUM=0.0E0 DO 20 J=1,N 20 SUM=SUM+X(J) SUM=SUM-FLOAT(N+1) DO 30 J=2,N 30 V(J)=SUM+X(J) TEMP = X(N+1) DO 40 L=1,N V(L) = TEMP*V(L)+(1.-TEMP)*X(L) 40 CONTINUE RETURN END C C SUBROUTINE DZF6(N,MXDM,DZ,X) REAL DZ(MXDM,N), X(N+1) INTEGER N,MXDM COMMON KJAC,KFUNC C *** COMPUTES THE DERIVATIVE(NXN MATRIX) OF THE FUNCTION F C WITH RESPECT TO U AT THE VECTOR (U(N),R), WHICH IS PASSED C IN X(N+1). THE DERIVATIVE IS RETURNED IN DZ(N,N). C C KJAC = KJAC+1 C C******************************************************************** C C THE JACOBIAN MATRIX FOR BROWN'S FUNCTION EVALUATED AT C THE POINT X, RETURNING THE VALUE IN V. C C******************************************************************** C DO 30 K=1,N PROD=1.0E0 DO 10 J=1,K-1 10 PROD=PROD*X(J) DO 15 J=K+1,N 15 PROD=PROD*X(J) DZ(1,K)=PROD DO 20 J=2,N 20 DZ(J,K)=1.0 IF (K .GT. 1) DZ(K,K)=DZ(K,K)+1.0 30 CONTINUE TEMP = X(N+1) DO 50 I=1,N DO 40 J=1,N DZ(J,I) = TEMP*DZ(J,I) 40 CONTINUE DZ(I,I) = 1.-TEMP +DZ(I,I) 50 CONTINUE RETURN END SUBROUTINE DRF6(N,DR,X) REAL DR(N),X(N+1) INTEGER N C *** COMPUTES THE DERIVATIVE(N-VECTOR) OF THE FUNCTION F C WITH RESPECT TO R AT (U,R) WHICH IS CONTAINED IN X(N+1) C THE DERIDRATIDRE IS RETURNED IN DR(N). PROD=1.0E0 DO 10 J=1,N 10 PROD=PROD*X(J) DR(1)=PROD-1.0E0 SUM=0.0E0 DO 20 J=1,N 20 SUM=SUM+X(J) SUM=SUM-FLOAT(N+1) DO 30 J=2,N 30 DR(J)=SUM+X(J) DO 40 J=1,N DR(J) = DR(J)- X(J) 40 CONTINUE RETURN END C *********** FREUDENSTEIN AND ROTH -- REGULARIZING ******* SUBROUTINE F7(N,RESULT,ARG) C ** FREUDENSTIEN AND ROTH ** REAL RESULT(N),ARG(N+1) INTEGER N COMMON KJAC,KFUNC KFUNC = KFUNC+1 C *** THIS SUBROUTINE COMPUTES F(X,R) WHERE F IS THE FUNCTION C WE WISH TO CONTINUE ALONG. ARG CONTAINS THE VALUE OF (X(N),R) C AT WHICH F IS COMPUTED. RESULT(N) CONTAINS F(X,R) ON RETURN. C ARG,N ARE UNCHANGED. C C RESULT(1)=ARG(1)+5.*ARG(2)**2-ARG(2)**3-2.*ARG(2)-13. RESULT(2)=ARG(1)+ARG(2)**2+ARG(2)**3-14.*ARG(2)-29. XL = ARG(3) RESULT(1) = XL*RESULT(1) + (1.-XL)*(ARG(1) - 15.) RESULT(2) = XL*RESULT(2) + (1.-XL)*(ARG(2) +2.) RETURN END C C SUBROUTINE DZF7(N,MXDM,DZ,ARG) REAL DZ(MXDM,N), ARG(N+1) INTEGER N,MXDM COMMON KJAC,KFUNC KJAC = KJAC + 1 C *** COMPUTES THE DERIVATIVE(NXN MATRIX) OF THE FUNCTION F C WITH RESPECT TO U AT THE VECTOR (U(N),R), WHICH IS PASSED C IN ARG(N+1). THE DERIVATIVE IS RETURNED IN DZ(N,N). C DZ(1,1) = 1.0 DZ(2,1) = 1.0 DZ(1,2) = -3.*ARG(2)**2+10.*ARG(2)-2. DZ(2,2) = 3.*ARG(2)**2+2.*ARG(2)-14. DO 20 J=1,2 DO 10 I=1,2 DZ(I,J) = ARG(3)*DZ(I,J) 10 CONTINUE DZ(J,J) = DZ(J,J) + 1.-ARG(3) 20 CONTINUE RETURN END C C C SUBROUTINE DRF7(N,DR,ARG) REAL DR(N),ARG(N+1) INTEGER N C *** COMPUTES THE DERIVATIVE(N-VECTOR) OF THE FUNCTION F C WITH RESPECT TO R AT (U,R) WHICH IS CONTAINED IN ARG(N+1) C THE DERIVATIVE IS RETURNED IN DR(N). C DR(1)=ARG(1)+5.*ARG(2)**2-ARG(2)**3-2.*ARG(2)-13. DR(2)=ARG(1)+ARG(2)**2+ARG(2)**3-14.*ARG(2)-29. DR(1) = DR(1) - ARG(1) + 15. DR(2) = DR(2) - ARG(2) - 2. RETURN END C ************ WATSON CURVE ******************* C SUBROUTINE FF(N,RESULT,ARG) REAL RESULT(N),ARG(N+1) INTEGER N C *** THIS SUBROUTINE COMPUTES F(X,R) WHERE F IS THE FUNCTION C WE WISH TO CONTINUE ALONG. ARG CONTAINS THE VALUE OF (X(N),R) C AT WHICH F IS COMPUTED. RESULT(N) CONTAINS F(X,R) ON RETURN. C ARG,N ARE UNCHANGED. C C SUM = 0.0 DO 10 L=1,N SUM = SUM + ARG(L) 10 CONTINUE DO 20 L=1,N RESULT(L)=EXP(COS(L*SUM)) 20 CONTINUE RETURN END C SUBROUTINE F8(N,RESULT,ARG) REAL RESULT(N),ARG(N+1) INTEGER N COMMON KJAC,KFUNC KFUNC = KFUNC + 1 C *** THIS SUBROUTINE COMPUTES F(X,R) WHERE F IS THE FUNCTION C WE WISH TO CONTINUE ALONG. ARG CONTAINS THE VALUE OF (X(N),R) C AT WHICH F IS COMPUTED. RESULT(N) CONTAINS F(X,R) ON RETURN. C ARG,N ARE UNCHANGED. C C CALL FF(N,RESULT,ARG) DO 10 L=1,N RESULT(L) = ARG(L) - ARG(N+1)*RESULT(L) 10 CONTINUE RETURN END C C C SUBROUTINE DZF8(N,MXDM,DZ,ARG) REAL DZ(MXDM,N), ARG(N+1),FN(1001) INTEGER N,MXDM COMMON KJAC,KFUNC KJAC = KJAC+1 C *** COMPUTES THE DERIVATIVE(NXN MATRIX) OF THE FUNCTION F C WITH RESPECT TO U AT THE VECTOR (U(N),R), WHICH IS PASSED C IN ARG(N+1). THE DERIVATIVE IS RETURNED IN DZ(N,N). C SUM = 0. DO 10 L=1,N SUM = SUM + ARG(L) 10 CONTINUE CALL FF(N,FN,ARG) DO 20 I=1,N TEMP = I*FN(I)*SIN(I*SUM)*ARG(N+1) DO 15 J=1,N DZ(I,J) =TEMP 15 CONTINUE DZ(I,I) = DZ(I,I) + 1. 20 CONTINUE RETURN END C C C SUBROUTINE DRF8(N,DR,ARG) REAL DR(N),ARG(N+1) INTEGER N C *** COMPUTES THE DERIVATIVE(N-VECTOR) OF THE FUNCTION F C WITH RESPECT TO R AT (U,R) WHICH IS CONTAINED IN ARG(N+1) C THE DERIVATIVE IS RETURNED IN DR(N). C CALL FF(N,DR,ARG) DO 10 I=1,N DR(I) = -DR(I) 10 CONTINUE RETURN END SUBROUTINE F9(N,RESULT,ARG) PARAMETER(A=-25.,B=-25.) REAL RESULT(N),ARG(N+1) INTEGER N COMMON KJAC,KFUNC C *** THIS SUBROUTINE COMPUTES F(X,R) WHERE F IS THE FUNCTION C WE WISH TO CONTINUE ALONG. ARG CONTAINS THE VALUE OF (X(N),R) C AT WHICH F IS COMPUTED. RESULT(N) CONTAINS F(X,R) ON RETURN. C ARG,N ARE UNCHANGED. C C C KFUNC = KFUNC + 1 RESULT(1) = ARG(1)*EXP(A*ARG(3)) RESULT(2) = ARG(2)*EXP(B*ARG(3)) RETURN END C C C SUBROUTINE DZF9(N,MXDM,DZ,ARG) PARAMETER(A=-25.,B=-25.) REAL DZ(MXDM,N), ARG(N+1) INTEGER N,MXDM COMMON KJAC,KFUNC C *** COMPUTES THE DERIVATIVE(NXN MATRIX) OF THE FUNCTION F C WITH RESPECT TO U AT THE VECTOR (U(N),R), WHICH IS PASSED C IN ARG(N+1). THE DERIVATIVE IS RETURNED IN DZ(N,N). C C KJAC = KJAC+1 C DO 20 J=1,N DO 10 I=1,N DZ(I,J) = 0. 10 CONTINUE 20 CONTINUE C DZ(1,1) =EXP(A*ARG(3)) DZ(2,2) = EXP(B*ARG(3)) RETURN END C C C SUBROUTINE DRF9(N,DR,ARG) PARAMETER(A=-25.,B=-25.) REAL DR(N),ARG(N+1) INTEGER N C *** COMPUTES THE DERIVATIVE(N-VECTOR) OF THE FUNCTION F C WITH RESPECT TO R AT (U,R) WHICH IS CONTAINED IN ARG(N+1) C THE DERIVATIVE IS RETURNED IN DR(N). C DR(1) = ARG(1)*A*EXP(A*ARG(3)) DR(2) = B*ARG(2)*EXP(B*ARG(3)) RETURN END