* File: SIVSET.for 28 subrs used by all main subrs of the * David Gay & Linda Kaufman nonlinear LS package. c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. c ALL RIGHTS RESERVED. c Based on Government Sponsored Research NAS7-03001. C>> 2001-06-18 SIVSET Krogh Replaced ". AND." with " .AND." C>> 1998-10-29 SIVSET Krogh Moved external statements up for mangle. C>> 1996-06-18 SIVSET Krogh More work for for C conversion. c>> 1995-11-15 SIVSET Krogh Moved formats up for C conversion. c>> 1994-11-02 SIVSET Krogh Changes to use M77CON c>> 1992-04-27 SIVSET CLL Removed unreferenced stmt labels. c>> 1992-04-13 CLL Change from Hollerith to '...' syntax in formats. c>> 1991-05-21 CLL JPL Changing (1) to (*) in declarations. c>> 1990-06-15 CLL JPL c>> 1990-03-21 CLL JPL c>> 1990-02-16 CLL JPL SUBROUTINE SIVSET(ALG, IV, LIV, LV, V) C C *** SUPPLY ***SOL (VERSION 2.3) DEFAULT VALUES TO IV AND V *** C C *** ALG = 1 MEANS REGRESSION CONSTANTS. C *** ALG = 2 MEANS GENERAL UNCONSTRAINED OPTIMIZATION CONSTANTS. c ------------------------------------------------------------------ c 1990-06-25 CLL changed default settings for printing so the c default is no printing. This involved setting COVPRT, OUTLEV, C PARPRT, SOLPRT, STATPR, and X0PRT to zero. Previously these were c all set to 1 except that COVPRT = 3. c ------------------------------------------------------------------ C INTEGER LIV, LV INTEGER ALG, IV(LIV) REAL V(LV) C EXTERNAL SV7DFL, I1MACH INTEGER I1MACH c--S replaces "?": ?IVSET,?V7DFL,?R7MDC,?L7SVN,?D7TPR,?V2AXY,?V2NRM c--& ?Q7APL,?D7UPD,?V7SCP,?Q7RAD,?V7SCL,?PARCK,?V7CPY,?S7LVM,?S7LUP c--& ?L7MST,?L7ITV,?L7IVM,?G7QTS,?L7SRT,?L7SVX,?A7SST,?ITSUM,?L7SQR c--& ?L7TVM,?L7VML,?RLDST C I1MACH... RETURNS MACHINE-DEPENDENT INTEGER CONSTANTS. C SV7DFL.... PROVIDES DEFAULT VALUES TO V. C INTEGER ALG1, MIV, MV INTEGER MINIV(4), MINV(4) C C *** SUBSCRIPTS FOR IV *** C INTEGER ALGSAV, COVPRT, COVREQ, DRADPR, DTYPE, HC, IERR, INITH, 1 INITS, IPIVOT, IVNEED, LASTIV, LASTV, LMAT, MXFCAL, 2 MXITER, NFCOV, NGCOV, NVDFLT, NVSAVE, OUTLEV, PARPRT, 3 PARSAV, PERM, PRUNIT, QRTYP, RDREQ, RMAT, SOLPRT, STATPR, 4 VNEED, VSAVE, X0PRT C C *** IV SUBSCRIPT VALUES *** C PARAMETER (ALGSAV=51, COVPRT=14, COVREQ=15, DRADPR=101, DTYPE=16, 1 HC=71, IERR=75, INITH=25, INITS=25, IPIVOT=76, 2 IVNEED=3, LASTIV=44, LASTV=45, LMAT=42, MXFCAL=17, 3 MXITER=18, NFCOV=52, NGCOV=53, NVDFLT=50, NVSAVE=9, 4 OUTLEV=19, PARPRT=20, PARSAV=49, PERM=58, PRUNIT=21, 5 QRTYP=80, RDREQ=57, RMAT=78, SOLPRT=22, STATPR=23, 6 VNEED=4, VSAVE=60, X0PRT=24) DATA MINIV(1)/82/, MINIV(2)/59/, MINIV(3)/103/, MINIV(4)/103/, 1 MINV(1)/98/, MINV(2)/71/, MINV(3)/101/, MINV(4)/85/ C C ------------------------------ BODY -------------------------------- C C I1MACH(2) returns the unit number for standard output. c IF (PRUNIT .LE. LIV) IV(PRUNIT) = I1MACH(2) IF (ALGSAV .LE. LIV) IV(ALGSAV) = ALG IF (ALG .LT. 1 .OR. ALG .GT. 4) GO TO 40 MIV = MINIV(ALG) IF (LIV .LT. MIV) GO TO 20 MV = MINV(ALG) IF (LV .LT. MV) GO TO 30 ALG1 = MOD(ALG-1,2) + 1 CALL SV7DFL(ALG1, LV, V) IV(1) = 12 IF (ALG .GT. 2) IV(DRADPR) = 1 IV(IVNEED) = 0 IV(LASTIV) = MIV IV(LASTV) = MV IV(LMAT) = MV + 1 IV(MXFCAL) = 200 IV(MXITER) = 150 IV(OUTLEV) = 0 IV(PARPRT) = 0 IV(PERM) = MIV + 1 IV(SOLPRT) = 0 IV(STATPR) = 0 IV(VNEED) = 0 IV(X0PRT) = 0 C IF (ALG1 .GE. 2) GO TO 10 C C *** REGRESSION VALUES C IV(COVPRT) = 0 IV(COVREQ) = 1 IV(DTYPE) = 1 IV(HC) = 0 IV(IERR) = 0 IV(INITS) = 0 IV(IPIVOT) = 0 IV(NVDFLT) = 32 IV(VSAVE) = 58 IF (ALG .GT. 2) IV(VSAVE) = IV(VSAVE) + 3 IV(PARSAV) = IV(VSAVE) + NVSAVE IV(QRTYP) = 1 IV(RDREQ) = 3 IV(RMAT) = 0 GO TO 999 C C *** GENERAL OPTIMIZATION VALUES C 10 IV(DTYPE) = 0 IV(INITH) = 1 IV(NFCOV) = 0 IV(NGCOV) = 0 IV(NVDFLT) = 25 IV(PARSAV) = 47 IF (ALG .GT. 2) IV(PARSAV) = 61 GO TO 999 C 20 IV(1) = 15 GO TO 999 C 30 IV(1) = 16 GO TO 999 C 40 IV(1) = 67 C 999 RETURN C *** LAST CARD OF SIVSET FOLLOWS *** END SUBROUTINE SV7DFL(ALG, LV, V) C C *** SUPPLY ***SOL (VERSION 2.3) DEFAULT VALUES TO V *** C C *** ALG = 1 MEANS REGRESSION CONSTANTS. C *** ALG = 2 MEANS GENERAL UNCONSTRAINED OPTIMIZATION CONSTANTS. C INTEGER ALG, LV REAL V(LV) C EXTERNAL SR7MDC REAL SR7MDC C SR7MDC... RETURNS MACHINE-DEPENDENT CONSTANTS C REAL MACHEP, MEPCRT, ONE, SQTEPS, THREE C C *** SUBSCRIPTS FOR V *** C INTEGER AFCTOL, BIAS, COSMIN, DECFAC, DELTA0, DFAC, DINIT, DLTFDC, 1 DLTFDJ, DTINIT, D0INIT, EPSLON, ETA0, FUZZ, HUBERC, 2 INCFAC, LMAX0, LMAXS, PHMNFC, PHMXFC, RDFCMN, RDFCMX, 3 RFCTOL, RLIMIT, RSPTOL, SCTOL, SIGMIN, TUNER1, TUNER2, 4 TUNER3, TUNER4, TUNER5, XCTOL, XFTOL C PARAMETER (ONE=1.E+0, THREE=3.E+0) C C *** V SUBSCRIPT VALUES *** C PARAMETER (AFCTOL=31, BIAS=43, COSMIN=47, DECFAC=22, DELTA0=44, 1 DFAC=41, DINIT=38, DLTFDC=42, DLTFDJ=43, DTINIT=39, 2 D0INIT=40, EPSLON=19, ETA0=42, FUZZ=45, HUBERC=48, 3 INCFAC=23, LMAX0=35, LMAXS=36, PHMNFC=20, PHMXFC=21, 4 RDFCMN=24, RDFCMX=25, RFCTOL=32, RLIMIT=46, RSPTOL=49, 5 SCTOL=37, SIGMIN=50, TUNER1=26, TUNER2=27, TUNER3=28, 6 TUNER4=29, TUNER5=30, XCTOL=33, XFTOL=34) C ------------------------------ BODY -------------------------------- C MACHEP = SR7MDC(3) V(AFCTOL) = 1.E-20 IF (MACHEP .GT. 1.E-10) V(AFCTOL) = MACHEP**2 V(DECFAC) = 0.5E+0 SQTEPS = SR7MDC(4) V(DFAC) = 0.6E+0 V(DTINIT) = 1.E-6 MEPCRT = MACHEP ** (ONE/THREE) V(D0INIT) = 1.E+0 V(EPSLON) = 0.1E+0 V(INCFAC) = 2.E+0 V(LMAX0) = 1.E+0 V(LMAXS) = 1.E+0 V(PHMNFC) = -0.1E+0 V(PHMXFC) = 0.1E+0 V(RDFCMN) = 0.1E+0 V(RDFCMX) = 4.E+0 V(RFCTOL) = max(1.E-10, MEPCRT**2) V(SCTOL) = V(RFCTOL) V(TUNER1) = 0.1E+0 V(TUNER2) = 1.E-4 V(TUNER3) = 0.75E+0 V(TUNER4) = 0.5E+0 V(TUNER5) = 0.75E+0 V(XCTOL) = SQTEPS V(XFTOL) = 1.E+2 * MACHEP C IF (ALG .GE. 2) GO TO 10 C C *** REGRESSION VALUES C V(COSMIN) = max(1.E-6, 1.E+2 * MACHEP) V(DINIT) = 0.E+0 V(DELTA0) = SQTEPS V(DLTFDC) = MEPCRT V(DLTFDJ) = SQTEPS V(FUZZ) = 1.5E+0 V(HUBERC) = 0.7E+0 V(RLIMIT) = SR7MDC(5) V(RSPTOL) = 1.E-3 V(SIGMIN) = 1.E-4 GO TO 999 C C *** GENERAL OPTIMIZATION VALUES C 10 V(BIAS) = 0.8E+0 V(DINIT) = -1.0E+0 V(ETA0) = 1.0E+3 * MACHEP C 999 RETURN C *** LAST CARD OF SV7DFL FOLLOWS *** END REAL FUNCTION SR7MDC(K) C C *** RETURN MACHINE DEPENDENT CONSTANTS USED BY NL2SOL *** C INTEGER K C C *** THE CONSTANT RETURNED DEPENDS ON K... C C *** K = 1... SMALLEST POS. ETA SUCH THAT -ETA EXISTS. C *** K = 2... SQUARE ROOT OF ETA. C *** K = 3... UNIT ROUNDOFF = SMALLEST POS. NO. MACHEP SUCH C *** THAT 1 + MACHEP .GT. 1 .AND. 1 - MACHEP .LT. 1. C *** K = 4... SQUARE ROOT OF MACHEP. C *** K = 5... SQUARE ROOT OF BIG (SEE K = 6). C *** K = 6... LARGEST MACHINE NO. BIG SUCH THAT -BIG EXISTS. C REAL BIG, ETA, MACHEP C EXTERNAL R1MACH REAL R1MACH, ZERO DATA BIG/0.E+0/, ETA/0.E+0/, MACHEP/0.E+0/, ZERO/0.E+0/ IF (BIG .GT. ZERO) GO TO 1 BIG = R1MACH(2) ETA = R1MACH(1) MACHEP = R1MACH(4) 1 CONTINUE C C ------------------------------ BODY -------------------------------- C GO TO (10, 20, 30, 40, 50, 60), K C 10 SR7MDC = ETA GO TO 999 C 20 SR7MDC = sqrt(256.E+0*ETA)/16.E+0 GO TO 999 C 30 SR7MDC = MACHEP GO TO 999 C 40 SR7MDC = sqrt(MACHEP) GO TO 999 C 50 SR7MDC = sqrt(BIG/256.E+0)*16.E+0 GO TO 999 C 60 SR7MDC = BIG C 999 RETURN C *** LAST CARD OF SR7MDC FOLLOWS *** END REAL FUNCTION SL7SVN(P, L, X, Y) C C *** ESTIMATE SMALLEST SING. VALUE OF PACKED LOWER TRIANG. MATRIX L C C *** PARAMETER DECLARATIONS *** C INTEGER P REAL L(*), X(P), Y(P) C DIMENSION L(P*(P+1)/2) C C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C *** PURPOSE *** C C THIS FUNCTION RETURNS A GOOD OVER-ESTIMATE OF THE SMALLEST C SINGULAR VALUE OF THE PACKED LOWER TRIANGULAR MATRIX L. C C *** PARAMETER DESCRIPTION *** C C P (IN) = THE ORDER OF L. L IS A P X P LOWER TRIANGULAR MATRIX. C L (IN) = ARRAY HOLDING THE ELEMENTS OF L IN ROW ORDER, I.E. C L(1,1), L(2,1), L(2,2), L(3,1), L(3,2), L(3,3), ETC. C X (OUT) IF SL7SVN RETURNS A POSITIVE VALUE, THEN X IS A NORMALIZED C APPROXIMATE LEFT SINGULAR VECTOR CORRESPONDING TO THE C SMALLEST SINGULAR VALUE. THIS APPROXIMATION MAY BE VERY C CRUDE. IF SL7SVN RETURNS ZERO, THEN SOME COMPONENTS OF X C ARE ZERO AND THE REST RETAIN THEIR INPUT VALUES. C Y (OUT) IF SL7SVN RETURNS A POSITIVE VALUE, THEN Y = (L**-1)*X IS AN C UNNORMALIZED APPROXIMATE RIGHT SINGULAR VECTOR CORRESPOND- C ING TO THE SMALLEST SINGULAR VALUE. THIS APPROXIMATION C MAY BE CRUDE. IF SL7SVN RETURNS ZERO, THEN Y RETAINS ITS C INPUT VALUE. THE CALLER MAY PASS THE SAME VECTOR FOR X C AND Y (NONSTANDARD FORTRAN USAGE), IN WHICH CASE Y OVER- C WRITES X (FOR NONZERO SL7SVN RETURNS). C C *** ALGORITHM NOTES *** C C THE ALGORITHM IS BASED ON (1), WITH THE ADDITIONAL PROVISION THAT C SL7SVN = 0 IS RETURNED IF THE SMALLEST DIAGONAL ELEMENT OF L C (IN MAGNITUDE) IS NOT MORE THAN THE UNIT ROUNDOFF TIMES THE C LARGEST. THE ALGORITHM USES A RANDOM NUMBER GENERATOR PROPOSED C IN (4), WHICH PASSES THE SPECTRAL TEST WITH FLYING COLORS -- SEE C (2) AND (3). C C *** SUBROUTINES AND FUNCTIONS CALLED *** C C SV2NRM - FUNCTION, RETURNS THE 2-NORM OF A VECTOR. C C *** REFERENCES *** C C (1) CLINE, A., MOLER, C., STEWART, G., AND WILKINSON, J.H.(1977), C AN ESTIMATE FOR THE CONDITION NUMBER OF A MATRIX, REPORT C TM-310, APPLIED MATH. DIV., ARGONNE NATIONAL LABORATORY. C C (2) HOAGLIN, D.C. (1976), THEORETICAL PROPERTIES OF CONGRUENTIAL C RANDOM-NUMBER GENERATORS -- AN EMPIRICAL VIEW, C MEMORANDUM NS-340, DEPT. OF STATISTICS, HARVARD UNIV. C C (3) KNUTH, D.E. (1969), THE ART OF COMPUTER PROGRAMMING, VOL. 2 C (SEMINUMERICAL ALGORITHMS), ADDISON-WESLEY, READING, MASS. C C (4) SMITH, C.S. (1971), MULTIPLICATIVE PSEUDO-RANDOM NUMBER C GENERATORS WITH PRIME MODULUS, J. ASSOC. COMPUT. MACH. 18, C PP. 586-593. C C *** HISTORY *** C C DESIGNED AND CODED BY DAVID M. GAY (WINTER 1977/SUMMER 1978). C C *** GENERAL *** C C THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH C SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS C MCS-7600324, DCR75-10143, 76-14311DSS, AND MCS76-11989. C C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C *** LOCAL VARIABLES *** C INTEGER I, II, IX, J, JI, JJ, JJJ, JM1, J0, PM1 REAL B, SMINUS, SPLUS, T, XMINUS, XPLUS C C *** CONSTANTS *** C REAL HALF, ONE, R9973, ZERO C C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C EXTERNAL SD7TPR, SV2NRM,SV2AXY REAL SD7TPR, SV2NRM C PARAMETER (HALF=0.5E+0, ONE=1.E+0, R9973=9973.E+0, ZERO=0.E+0) C C *** BODY *** C IX = 2 PM1 = P - 1 C C *** FIRST CHECK WHETHER TO RETURN SL7SVN = 0 AND INITIALIZE X *** C II = 0 J0 = P*PM1/2 JJ = J0 + P IF (L(JJ) .EQ. ZERO) GO TO 110 IX = MOD(3432*IX, 9973) B = HALF*(ONE + real(IX)/R9973) XPLUS = B / L(JJ) X(P) = XPLUS IF (P .LE. 1) GO TO 60 DO 10 I = 1, PM1 II = II + I IF (L(II) .EQ. ZERO) GO TO 110 JI = J0 + I X(I) = XPLUS * L(JI) 10 CONTINUE C C *** SOLVE (L**T)*X = B, WHERE THE COMPONENTS OF B HAVE RANDOMLY C *** CHOSEN MAGNITUDES IN (.5,1) WITH SIGNS CHOSEN TO MAKE X LARGE. C C DO J = P-1 TO 1 BY -1... DO 50 JJJ = 1, PM1 J = P - JJJ C *** DETERMINE X(J) IN THIS ITERATION. NOTE FOR I = 1,2,...,J C *** THAT X(I) HOLDS THE CURRENT PARTIAL SUM FOR ROW I. IX = MOD(3432*IX, 9973) B = HALF*(ONE + real(IX)/R9973) XPLUS = (B - X(J)) XMINUS = (-B - X(J)) SPLUS = abs(XPLUS) SMINUS = abs(XMINUS) JM1 = J - 1 J0 = J*JM1/2 JJ = J0 + J XPLUS = XPLUS/L(JJ) XMINUS = XMINUS/L(JJ) IF (JM1 .EQ. 0) GO TO 30 DO 20 I = 1, JM1 JI = J0 + I SPLUS = SPLUS + abs(X(I) + L(JI)*XPLUS) SMINUS = SMINUS + abs(X(I) + L(JI)*XMINUS) 20 CONTINUE 30 IF (SMINUS .GT. SPLUS) XPLUS = XMINUS X(J) = XPLUS C *** UPDATE PARTIAL SUMS *** IF (JM1 .GT. 0) CALL SV2AXY(JM1, X, XPLUS, L(J0+1), X) 50 CONTINUE C C *** NORMALIZE X *** C 60 T = ONE/SV2NRM(P, X) DO 70 I = 1, P 70 X(I) = T*X(I) C C *** SOLVE L*Y = X AND RETURN SL7SVN = 1/TWONORM(Y) *** C DO 100 J = 1, P JM1 = J - 1 J0 = J*JM1/2 JJ = J0 + J T = ZERO IF (JM1 .GT. 0) T = SD7TPR(JM1, L(J0+1), Y) Y(J) = (X(J) - T) / L(JJ) 100 CONTINUE C SL7SVN = ONE/SV2NRM(P, Y) GO TO 999 C 110 SL7SVN = ZERO 999 RETURN C *** LAST CARD OF SL7SVN FOLLOWS *** END SUBROUTINE SQ7APL(NN, N, P, J, R, IERR) C *****PARAMETERS. INTEGER NN, N, P, IERR REAL J(NN,P), R(N) C C .................................................................. C .................................................................. C C *****PURPOSE. C THIS SUBROUTINE APPLIES TO R THE ORTHOGONAL TRANSFORMATIONS C STORED IN J BY QRFACT C C *****PARAMETER DESCRIPTION. C ON INPUT. C C NN IS THE ROW DIMENSION OF THE MATRIX J AS DECLARED IN C THE CALLING PROGRAM DIMENSION STATEMENT C C N IS THE NUMBER OF ROWS OF J AND THE SIZE OF THE VECTOR R C C P IS THE NUMBER OF COLUMNS OF J AND THE SIZE OF SIGMA C C J CONTAINS ON AND BELOW ITS DIAGONAL THE COLUMN VECTORS C U WHICH DETERMINE THE HOUSEHOLDER TRANSFORMATIONS C IDENT - U*U.TRANSPOSE C C R IS THE RIGHT HAND SIDE VECTOR TO WHICH THE ORTHOGONAL C TRANSFORMATIONS WILL BE APPLIED C C IERR IF NON-ZERO INDICATES THAT NOT ALL THE TRANSFORMATIONS C WERE SUCCESSFULLY DETERMINED AND ONLY THE FIRST C ABS(IERR) - 1 TRANSFORMATIONS WILL BE USED C C ON OUTPUT. C C R HAS BEEN OVERWRITTEN BY ITS TRANSFORMED IMAGE C C *****APPLICATION AND USAGE RESTRICTIONS. C NONE C C *****ALGORITHM NOTES. C THE VECTORS U WHICH DETERMINE THE HOUSEHOLDER TRANSFORMATIONS C ARE NORMALIZED SO THAT THEIR 2-NORM SQUARED IS 2. THE USE OF C THESE TRANSFORMATIONS HERE IS IN THE SPIRIT OF (1). C C *****SUBROUTINES AND FUNCTIONS CALLED. C C SD7TPR - FUNCTION, RETURNS THE INNER PRODUCT OF VECTORS C C *****REFERENCES. C (1) BUSINGER, P. A., AND GOLUB, G. H. (1965), LINEAR LEAST SQUARES C SOLUTIONS BY HOUSEHOLDER TRANSFORMATIONS, NUMER. MATH. 7, C PP. 269-276. C C *****HISTORY. C DESIGNED BY DAVID M. GAY, CODED BY STEPHEN C. PETERS (WINTER 1977) C CALL ON SV2AXY SUBSTITUTED FOR DO LOOP, FALL 1983. C C *****GENERAL. C C THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH C SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS C MCS-7600324, DCR75-10143, 76-14311DSS, AND MCS76-11989. C C .................................................................. C .................................................................. C C *****LOCAL VARIABLES. INTEGER K, L, NL1 C *****FUNCTIONS. EXTERNAL SD7TPR,SV2AXY REAL SD7TPR c ------------------------------------------------------------------ C C *** BODY *** C K = P IF (IERR .NE. 0) K = abs(IERR) - 1 IF ( K .EQ. 0) GO TO 999 C DO 20 L = 1, K NL1 = N - L + 1 CALL SV2AXY(NL1, R(L), -SD7TPR(NL1,J(L,L),R(L)), J(L,L), R(L)) 20 CONTINUE C 999 RETURN C *** LAST LINE OF SQ7APL FOLLOWS *** END REAL FUNCTION SV2NRM(P, X) C C *** RETURN THE 2-NORM OF THE P-VECTOR X, TAKING *** C *** CARE TO AVOID THE MOST LIKELY UNDERFLOWS. *** C INTEGER P REAL X(P) C INTEGER I, J REAL ONE, R, SCALE, SQTETA, T, XI, ZERO EXTERNAL SR7MDC REAL SR7MDC c ------------------------------------------------------------------ PARAMETER (ONE=1.E+0, ZERO=0.E+0) SAVE SQTETA DATA SQTETA/0.E+0/ c ------------------------------------------------------------------ IF (P .GT. 0) GO TO 10 SV2NRM = ZERO GO TO 999 10 DO 20 I = 1, P IF (X(I) .NE. ZERO) GO TO 30 20 CONTINUE SV2NRM = ZERO GO TO 999 C 30 SCALE = abs(X(I)) IF (I .LT. P) GO TO 40 SV2NRM = SCALE GO TO 999 40 T = ONE IF (SQTETA .EQ. ZERO) SQTETA = SR7MDC(2) C C *** SQTETA IS (SLIGHTLY LARGER THAN) THE SQUARE ROOT OF THE C *** SMALLEST POSITIVE FLOATING POINT NUMBER ON THE MACHINE. C *** THE TESTS INVOLVING SQTETA ARE DONE TO PREVENT UNDERFLOWS. C J = I + 1 DO 60 I = J, P XI = abs(X(I)) IF (XI .GT. SCALE) GO TO 50 R = XI / SCALE IF (R .GT. SQTETA) T = T + R*R GO TO 60 50 R = SCALE / XI IF (R .LE. SQTETA) R = ZERO T = ONE + T * R*R SCALE = XI 60 CONTINUE C SV2NRM = SCALE * sqrt(T) 999 RETURN C *** LAST LINE OF SV2NRM FOLLOWS *** END REAL FUNCTION SD7TPR(P, X, Y) C C *** RETURN THE INNER PRODUCT OF THE P-VECTORS X AND Y. *** C INTEGER P REAL X(P), Y(P) C INTEGER I REAL ONE, SQTETA, T, ZERO EXTERNAL SR7MDC REAL SR7MDC C C *** SR7MDC(2) RETURNS A MACHINE-DEPENDENT CONSTANT, SQTETA, WHICH C *** IS SLIGHTLY LARGER THAN THE SMALLEST POSITIVE NUMBER THAT C *** CAN BE SQUARED WITHOUT UNDERFLOWING. C PARAMETER (ONE=1.E+0, ZERO=0.E+0) DATA SQTETA/0.E+0/ C SD7TPR = ZERO IF (P .LE. 0) GO TO 999 IF (SQTETA .EQ. ZERO) SQTETA = SR7MDC(2) DO 20 I = 1, P T = max(abs(X(I)), abs(Y(I))) IF (T .GT. ONE) GO TO 10 IF (T .LT. SQTETA) GO TO 20 T = (X(I)/SQTETA)*Y(I) IF (abs(T) .LT. SQTETA) GO TO 20 10 SD7TPR = SD7TPR + X(I)*Y(I) 20 CONTINUE C 999 RETURN C *** LAST LINE OF SD7TPR FOLLOWS *** END SUBROUTINE SD7UPD(D, DR, IV, LIV, LV, N, ND, NN, N2, P, V) C C *** UPDATE SCALE VECTOR D FOR NL2IT *** C C *** PARAMETER DECLARATIONS *** C INTEGER LIV, LV, N, ND, NN, N2, P INTEGER IV(LIV) REAL D(P), DR(ND,P), V(LV) C C *** LOCAL VARIABLES *** C INTEGER D0, I, JCN0, JCN1, JCNI, JTOL0, JTOLI, K, SII REAL T, VDFAC C C *** CONSTANTS *** C REAL ZERO C C *** EXTERNAL SUBROUTINE *** C EXTERNAL SV7SCP C C SV7SCP... SETS ALL COMPONENTS OF A VECTOR TO A SCALAR. C C *** SUBSCRIPTS FOR IV AND V *** C INTEGER DFAC, DTYPE, JCN, JTOL, NITER, S PARAMETER (DFAC=41, DTYPE=16, JCN=66, JTOL=59, NITER=31, S=62) PARAMETER (ZERO=0.E+0) C ------------------------------ BODY -------------------------------- C IF (IV(DTYPE) .NE. 1 .AND. IV(NITER) .GT. 0) GO TO 999 JCN1 = IV(JCN) JCN0 = abs(JCN1) - 1 IF (JCN1 .LT. 0) GO TO 10 IV(JCN) = -JCN1 CALL SV7SCP(P, V(JCN1), ZERO) 10 DO 30 I = 1, P JCNI = JCN0 + I T = V(JCNI) DO 20 K = 1, NN 20 T = max(T, abs(DR(K,I))) V(JCNI) = T 30 CONTINUE IF (N2 .LT. N) GO TO 999 VDFAC = V(DFAC) JTOL0 = IV(JTOL) - 1 D0 = JTOL0 + P SII = IV(S) - 1 DO 50 I = 1, P SII = SII + I JCNI = JCN0 + I T = V(JCNI) IF (V(SII) .GT. ZERO) T = max(sqrt(V(SII)), T) JTOLI = JTOL0 + I D0 = D0 + 1 IF (T .LT. V(JTOLI)) T = max(V(D0), V(JTOLI)) D(I) = max(VDFAC*D(I), T) 50 CONTINUE C 999 RETURN C *** LAST CARD OF SD7UPD FOLLOWS *** END SUBROUTINE SQ7RAD(N, NN, P, QTR, QTRSET, RMAT, W, Y) C C *** ADD ROWS W TO QR FACTORIZATION WITH R MATRIX RMAT AND C *** Q**T * RESIDUAL = QTR. Y = NEW COMPONENTS OF RESIDUAL C *** CORRESPONDING TO W. QTR, Y REFERENCED ONLY IF QTRSET = .TRUE. C LOGICAL QTRSET INTEGER N, NN, P REAL QTR(P), RMAT(*), W(NN,P), Y(N) C DIMENSION RMAT(P*(P+1)/2) EXTERNAL SD7TPR, SR7MDC,SV2AXY, SV7SCL, SV2NRM REAL SD7TPR, SR7MDC, SV2NRM C C *** LOCAL VARIABLES *** C INTEGER I, II, IJ, IP1, J, K, NK REAL ARI, QRI, RI, SS, T, WI REAL BIG, BIGRT, ONE, TINY, TINYRT, ZERO SAVE BIGRT, TINY, TINYRT DATA BIG/-1.E+0/, BIGRT/-1.E+0/, ONE/1.E+0/, TINY/0.E+0/, 1 TINYRT/0.E+0/, ZERO/0.E+0/ C C ----------------------------- BODY ----------------------------------- C IF (TINY .GT. ZERO) GO TO 10 TINY = SR7MDC(1) BIG = SR7MDC(6) IF (TINY*BIG .LT. ONE) TINY = ONE / BIG 10 K = 1 NK = N II = 0 DO 180 I = 1, P II = II + I IP1 = I + 1 IJ = II + I IF (NK .LE. 1) T = abs(W(K,I)) IF (NK .GT. 1) T = SV2NRM(NK, W(K,I)) IF (T .LT. TINY) GOTO 180 RI = RMAT(II) IF (RI .NE. ZERO) GO TO 100 IF (NK .GT. 1) GO TO 30 IJ = II DO 20 J = I, P RMAT(IJ) = W(K,J) IJ = IJ + J 20 CONTINUE IF (QTRSET) QTR(I) = Y(K) W(K,I) = ZERO GO TO 999 30 WI = W(K,I) IF (BIGRT .GT. ZERO) GO TO 40 BIGRT = SR7MDC(5) TINYRT = SR7MDC(2) 40 IF (T .LE. TINYRT) GO TO 50 IF (T .GE. BIGRT) GO TO 50 IF (WI .LT. ZERO) T = -T WI = WI + T SS = sqrt(T * WI) GO TO 70 50 SS = sqrt(T) IF (WI .LT. ZERO) GO TO 60 WI = WI + T SS = SS * sqrt(WI) GO TO 70 60 T = -T WI = WI + T SS = SS * sqrt(-WI) 70 W(K,I) = WI CALL SV7SCL(NK, W(K,I), ONE/SS, W(K,I)) RMAT(II) = -T IF (.NOT. QTRSET) GO TO 80 CALL SV2AXY(NK, Y(K), -SD7TPR(NK,Y(K),W(K,I)), W(K,I), Y(K)) QTR(I) = Y(K) 80 IF (IP1 .GT. P) GO TO 999 DO 90 J = IP1, P CALL SV2AXY(NK, W(K,J), -SD7TPR(NK,W(K,J),W(K,I)), 1 W(K,I), W(K,J)) RMAT(IJ) = W(K,J) IJ = IJ + J 90 CONTINUE IF (NK .LE. 1) GO TO 999 K = K + 1 NK = NK - 1 GO TO 180 C 100 ARI = abs(RI) IF (ARI .GT. T) GO TO 110 T = T * sqrt(ONE + (ARI/T)**2) GO TO 120 110 T = ARI * sqrt(ONE + (T/ARI)**2) 120 IF (RI .LT. ZERO) T = -T RI = RI + T RMAT(II) = -T SS = -RI / T IF (NK .LE. 1) GO TO 150 CALL SV7SCL(NK, W(K,I), ONE/RI, W(K,I)) IF (.NOT. QTRSET) GO TO 130 QRI = QTR(I) T = SS * ( QRI + SD7TPR(NK, Y(K), W(K,I)) ) QTR(I) = QRI + T 130 IF (IP1 .GT. P) GO TO 999 IF (QTRSET) CALL SV2AXY(NK, Y(K), T, W(K,I), Y(K)) DO 140 J = IP1, P RI = RMAT(IJ) T = SS * ( RI + SD7TPR(NK, W(K,J), W(K,I)) ) CALL SV2AXY(NK, W(K,J), T, W(K,I), W(K,J)) RMAT(IJ) = RI + T IJ = IJ + J 140 CONTINUE GO TO 180 C 150 WI = W(K,I) / RI W(K,I) = WI IF (.NOT. QTRSET) GO TO 160 QRI = QTR(I) T = SS * ( QRI + Y(K)*WI ) QTR(I) = QRI + T 160 IF (IP1 .GT. P) GO TO 999 IF (QTRSET) Y(K) = T*WI + Y(K) DO 170 J = IP1, P RI = RMAT(IJ) T = SS * (RI + W(K,J)*WI) W(K,J) = W(K,J) + T*WI RMAT(IJ) = RI + T IJ = IJ + J 170 CONTINUE 180 CONTINUE C 999 RETURN C *** LAST LINE OF SQ7RAD FOLLOWS *** END SUBROUTINE SPARCK(ALG, D, IV, LIV, LV, N, V) C C *** CHECK ***SOL (VERSION 2.3) PARAMETERS, PRINT CHANGED VALUES *** C C *** ALG = 1 FOR REGRESSION, ALG = 2 FOR GENERAL UNCONSTRAINED OPT. C INTEGER ALG, LIV, LV, N INTEGER IV(LIV) REAL D(N), V(LV) C EXTERNAL SIVSET, SR7MDC,SV7CPY,SV7DFL REAL SR7MDC C SIVSET -- SUPPLIES DEFAULT VALUES TO BOTH IV AND V. C SR7MDC -- RETURNS MACHINE-DEPENDENT CONSTANTS. C SV7CPY -- COPIES ONE VECTOR TO ANOTHER. C SV7DFL -- SUPPLIES DEFAULT PARAMETER VALUES TO V ALONE. C C *** LOCAL VARIABLES *** C INTEGER ALG1, I, II, IV1, J, K, L, M, MIV1, MIV2, NDFALT, PARSV1, 1 PU INTEGER IJMP, JLIM(4), MINIV(4), NDFLT(4) CHARACTER*1 VARNM(2), SH(2) CHARACTER*4 CNGD(3), DFLT(3), VN(2,34), WHICH(3) REAL BIG, MACHEP, TINY, VK, VM(34), VX(34), ZERO C C *** IV AND V SUBSCRIPTS *** C INTEGER ALGSAV, DINIT, DTYPE, DTYPE0, EPSLON, INITS, IVNEED, 1 LASTIV, LASTV, LMAT, NEXTIV, NEXTV, NVDFLT, OLDN, 2 PARPRT, PARSAV, PERM, PRUNIT, VNEED C C PARAMETER (ALGSAV=51, DINIT=38, DTYPE=16, DTYPE0=54, EPSLON=19, 1 INITS=25, IVNEED=3, LASTIV=44, LASTV=45, LMAT=42, 2 NEXTIV=46, NEXTV=47, NVDFLT=50, OLDN=38, PARPRT=20, 3 PARSAV=49, PERM=58, PRUNIT=21, VNEED=4) SAVE BIG, MACHEP, TINY C DATA BIG/0.E+0/, MACHEP/-1.E+0/, TINY/1.E+0/, ZERO/0.E+0/ DATA VN(1,1),VN(2,1)/'EPSL','ON..'/ DATA VN(1,2),VN(2,2)/'PHMN','FC..'/ DATA VN(1,3),VN(2,3)/'PHMX','FC..'/ DATA VN(1,4),VN(2,4)/'DECF','AC..'/ DATA VN(1,5),VN(2,5)/'INCF','AC..'/ DATA VN(1,6),VN(2,6)/'RDFC','MN..'/ DATA VN(1,7),VN(2,7)/'RDFC','MX..'/ DATA VN(1,8),VN(2,8)/'TUNE','R1..'/ DATA VN(1,9),VN(2,9)/'TUNE','R2..'/ DATA VN(1,10),VN(2,10)/'TUNE','R3..'/ DATA VN(1,11),VN(2,11)/'TUNE','R4..'/ DATA VN(1,12),VN(2,12)/'TUNE','R5..'/ DATA VN(1,13),VN(2,13)/'AFCT','OL..'/ DATA VN(1,14),VN(2,14)/'RFCT','OL..'/ DATA VN(1,15),VN(2,15)/'XCTO','L...'/ DATA VN(1,16),VN(2,16)/'XFTO','L...'/ DATA VN(1,17),VN(2,17)/'LMAX','0...'/ DATA VN(1,18),VN(2,18)/'LMAX','S...'/ DATA VN(1,19),VN(2,19)/'SCTO','L...'/ DATA VN(1,20),VN(2,20)/'DINI','T...'/ DATA VN(1,21),VN(2,21)/'DTIN','IT..'/ DATA VN(1,22),VN(2,22)/'D0IN','IT..'/ DATA VN(1,23),VN(2,23)/'DFAC','....'/ DATA VN(1,24),VN(2,24)/'DLTF','DC..'/ DATA VN(1,25),VN(2,25)/'DLTF','DJ..'/ DATA VN(1,26),VN(2,26)/'DELT','A0..'/ DATA VN(1,27),VN(2,27)/'FUZZ','....'/ DATA VN(1,28),VN(2,28)/'RLIM','IT..'/ DATA VN(1,29),VN(2,29)/'COSM','IN..'/ DATA VN(1,30),VN(2,30)/'HUBE','RC..'/ DATA VN(1,31),VN(2,31)/'RSPT','OL..'/ DATA VN(1,32),VN(2,32)/'SIGM','IN..'/ DATA VN(1,33),VN(2,33)/'ETA0','....'/ DATA VN(1,34),VN(2,34)/'BIAS','....'/ C DATA VM(1)/1.0E-3/, VM(2)/-0.99E+0/, VM(3)/1.0E-3/, VM(4)/1.0E-2/, 1 VM(5)/1.2E+0/, VM(6)/1.E-2/, VM(7)/1.2E+0/, VM(8)/0.E+0/, 2 VM(9)/0.E+0/, VM(10)/1.E-3/, VM(11)/-1.E+0/, VM(13)/0.E+0/, 3 VM(15)/0.E+0/, VM(16)/0.E+0/, VM(19)/0.E+0/, VM(20)/-10.E+0/, 4 VM(21)/0.E+0/, VM(22)/0.E+0/, VM(23)/0.E+0/, VM(27)/1.01E+0/, 5 VM(28)/1.E+10/, VM(30)/0.E+0/, VM(31)/0.E+0/, VM(32)/0.E+0/, 6 VM(34)/0.E+0/ DATA VX(1)/0.9E+0/, VX(2)/-1.E-3/, VX(3)/1.E+1/, VX(4)/0.8E+0/, 1 VX(5)/1.E+2/, VX(6)/0.8E+0/, VX(7)/1.E+2/, VX(8)/0.5E+0/, 2 VX(9)/0.5E+0/, VX(10)/1.E+0/, VX(11)/1.E+0/, VX(14)/0.1E+0/, 3 VX(15)/1.E+0/, VX(16)/1.E+0/, VX(19)/1.E+0/, VX(23)/1.E+0/, 4 VX(24)/1.E+0/, VX(25)/1.E+0/, VX(26)/1.E+0/, VX(27)/1.E+10/, 5 VX(29)/1.E+0/, VX(31)/1.E+0/, VX(32)/1.E+0/, VX(33)/1.E+0/, 6 VX(34)/1.E+0/ C DATA VARNM(1)/'P'/, VARNM(2)/'P'/, SH(1)/'S'/, SH(2)/'H'/ DATA CNGD(1),CNGD(2),CNGD(3)/'---C','HANG','ED V'/, 1 DFLT(1),DFLT(2),DFLT(3)/'NOND','EFAU','LT V'/ DATA IJMP/33/, JLIM(1)/0/, JLIM(2)/24/, JLIM(3)/0/, JLIM(4)/24/, 1 NDFLT(1)/32/, NDFLT(2)/25/, NDFLT(3)/32/, NDFLT(4)/25/ DATA MINIV(1)/82/, MINIV(2)/59/, MINIV(3)/103/, MINIV(4)/103/ C C............................... BODY ................................ C 10 FORMAT(/' THE FIRST PARAMETER TO SIVSET SHOULD BE',I3, 1 ' RATHER THAN',I3) 40 FORMAT(/' /// BAD',A1,' =',I5) 70 FORMAT(/' /// ',1A1,' CHANGED FROM ',I5,' TO ',I5) 90 FORMAT(/' /// IV(1) =',I5,' SHOULD BE BETWEEN 0 AND 14.') 130 FORMAT(/' /// ',2A4,'.. V(',I2,') =',g11.3,' should', 1 ' be between',g11.3,' and',g11.3) 160 FORMAT(/' IV(NVDFLT) =',I5,' RATHER THAN ',I5) 180 FORMAT(/' /// D(',I3,') =',g11.3,' SHOULD BE POSITIVE') 220 FORMAT(/' NONDEFAULT VALUES....'/' INIT',A1,'..... IV(25) =',I3) 240 FORMAT(' DTYPE..... IV(16) =',I3) 260 FORMAT(/' ',3A4,'ALUES....'/) 270 FORMAT(1X,2A4,'.. V(',I2,') =',g15.7) 310 FORMAT(/' /// LIV =',I5,' MUST BE AT LEAST',I5) 330 FORMAT(/' /// LV =',I5,' MUST BE AT LEAST',I5) 350 FORMAT(/' /// ALG =',I5,' MUST BE 1 2, 3, OR 4') 370 FORMAT(/' /// LIV =',I5,' MUST BE AT LEAST',I5, 1 ' TO COMPUTE TRUE MIN. LIV AND MIN. LV') c c++(~.C.) Default UNITNO='(PU,' c++(.C.) Default UNITNO='(*,' c++ Replace "(PU," = UNITNO PU = 0 IF (PRUNIT .LE. LIV) PU = IV(PRUNIT) IF (ALGSAV .GT. LIV) GO TO 20 IF (ALG .EQ. IV(ALGSAV)) GO TO 20 IF (PU .NE. 0) WRITE(PU,10) ALG, IV(ALGSAV) IV(1) = 67 GO TO 999 20 IF (ALG .LT. 1 .OR. ALG .GT. 4) GO TO 340 MIV1 = MINIV(ALG) IF (IV(1) .EQ. 15) GO TO 360 ALG1 = MOD(ALG-1,2) + 1 IF (IV(1) .EQ. 0) CALL SIVSET(ALG, IV, LIV, LV, V) IV1 = IV(1) IF (IV1 .NE. 13 .AND. IV1 .NE. 12) GO TO 30 IF (PERM .LE. LIV) MIV1 = max(MIV1, IV(PERM) - 1) IF (IVNEED .LE. LIV) MIV2 = MIV1 + max(IV(IVNEED), 0) IF (LASTIV .LE. LIV) IV(LASTIV) = MIV2 IF (LIV .LT. MIV1) GO TO 300 IV(IVNEED) = 0 IV(LASTV) = max(IV(VNEED), 0) + IV(LMAT) - 1 IV(VNEED) = 0 IF (LIV .LT. MIV2) GO TO 300 IF (LV .LT. IV(LASTV)) GO TO 320 30 IF (IV1 .LT. 12 .OR. IV1 .GT. 14) GO TO 60 IF (N .GE. 1) GO TO 50 IV(1) = 81 IF (PU .EQ. 0) GO TO 999 WRITE(PU,40) VARNM(ALG1), N GO TO 999 50 IF (IV1 .NE. 14) IV(NEXTIV) = IV(PERM) IF (IV1 .NE. 14) IV(NEXTV) = IV(LMAT) IF (IV1 .EQ. 13) GO TO 999 K = IV(PARSAV) - EPSLON CALL SV7DFL(ALG1, LV-K, V(K+1)) IV(DTYPE0) = 2 - ALG1 IV(OLDN) = N WHICH(1) = DFLT(1) WHICH(2) = DFLT(2) WHICH(3) = DFLT(3) GO TO 110 60 IF (N .EQ. IV(OLDN)) GO TO 80 IV(1) = 17 IF (PU .EQ. 0) GO TO 999 WRITE(PU,70) VARNM(ALG1), IV(OLDN), N GO TO 999 C 80 IF (IV1 .LE. 11 .AND. IV1 .GE. 1) GO TO 100 IV(1) = 80 IF (PU .NE. 0) WRITE(PU,90) IV1 GO TO 999 C 100 WHICH(1) = CNGD(1) WHICH(2) = CNGD(2) WHICH(3) = CNGD(3) C 110 IF (IV1 .EQ. 14) IV1 = 12 IF (BIG .GT. TINY) GO TO 120 TINY = SR7MDC(1) MACHEP = SR7MDC(3) BIG = SR7MDC(6) VM(12) = MACHEP VX(12) = BIG VX(13) = BIG VM(14) = MACHEP VM(17) = TINY VX(17) = BIG VM(18) = TINY VX(18) = BIG VX(20) = BIG VX(21) = BIG VX(22) = BIG VM(24) = MACHEP VM(25) = MACHEP VM(26) = MACHEP VX(28) = SR7MDC(5) VM(29) = MACHEP VX(30) = BIG VM(33) = MACHEP 120 M = 0 I = 1 J = JLIM(ALG1) K = EPSLON NDFALT = NDFLT(ALG1) DO 150 L = 1, NDFALT VK = V(K) IF (VK .GE. VM(I) .AND. VK .LE. VX(I)) GO TO 140 M = K IF (PU .NE. 0) WRITE(PU,130) VN(1,I), VN(2,I), K, VK, 1 VM(I), VX(I) 140 K = K + 1 I = I + 1 IF (I .EQ. J) I = IJMP 150 CONTINUE C IF (IV(NVDFLT) .EQ. NDFALT) GO TO 170 IV(1) = 51 IF (PU .EQ. 0) GO TO 999 WRITE(PU,160) IV(NVDFLT), NDFALT GO TO 999 170 IF ((IV(DTYPE) .GT. 0 .OR. V(DINIT) .GT. ZERO) .AND. IV1 .EQ. 12) 1 GO TO 200 DO 190 I = 1, N IF (D(I) .GT. ZERO) GO TO 190 M = 18 IF (PU .NE. 0) WRITE(PU,180) I, D(I) 190 CONTINUE 200 IF (M .EQ. 0) GO TO 210 IV(1) = M GO TO 999 C 210 IF (PU .EQ. 0 .OR. IV(PARPRT) .EQ. 0) GO TO 999 IF (IV1 .NE. 12 .OR. IV(INITS) .EQ. ALG1-1) GO TO 230 M = 1 WRITE(PU,220) SH(ALG1), IV(INITS) 230 IF (IV(DTYPE) .EQ. IV(DTYPE0)) GO TO 250 IF (M .EQ. 0) WRITE(PU,260) WHICH M = 1 WRITE(PU,240) IV(DTYPE) 250 I = 1 J = JLIM(ALG1) K = EPSLON L = IV(PARSAV) NDFALT = NDFLT(ALG1) DO 290 II = 1, NDFALT IF (V(K) .EQ. V(L)) GO TO 280 IF (M .EQ. 0) WRITE(PU,260) WHICH M = 1 WRITE(PU,270) VN(1,I), VN(2,I), K, V(K) 280 K = K + 1 L = L + 1 I = I + 1 IF (I .EQ. J) I = IJMP 290 CONTINUE C IV(DTYPE0) = IV(DTYPE) PARSV1 = IV(PARSAV) CALL SV7CPY(IV(NVDFLT), V(PARSV1), V(EPSLON)) GO TO 999 C 300 IV(1) = 15 IF (PU .EQ. 0) GO TO 999 WRITE(PU,310) LIV, MIV2 IF (LIV .LT. MIV1) GO TO 999 IF (LV .LT. IV(LASTV)) GO TO 320 GO TO 999 C 320 IV(1) = 16 IF (PU .NE. 0) WRITE(PU,330) LV, IV(LASTV) GO TO 999 C 340 IV(1) = 67 IF (PU .NE. 0) WRITE(PU,350) ALG GO TO 999 360 IF (PU .NE. 0) WRITE(PU,370) LIV, MIV1 IF (LASTIV .LE. LIV) IV(LASTIV) = MIV1 IF (LASTV .LE. LIV) IV(LASTV) = 0 C 999 RETURN C *** LAST LINE OF SPARCK FOLLOWS *** END SUBROUTINE SS7LVM(P, Y, SS, X) C C *** SET Y = SS * X, SS = P X P SYMMETRIC MATRIX. *** C *** LOWER TRIANGLE OF SS STORED ROWWISE. *** C C *** PARAMETER DECLARATIONS *** C INTEGER P REAL SS(*), X(P), Y(P) C DIMENSION SS(P*(P+1)/2) C C *** LOCAL VARIABLES *** C INTEGER I, IM1, J, K REAL XI C C *** EXTERNAL FUNCTION *** C EXTERNAL SD7TPR REAL SD7TPR C C ---------------------------------------------------------------------- C J = 1 DO 10 I = 1, P Y(I) = SD7TPR(I, SS(J), X) J = J + I 10 CONTINUE C IF (P .LE. 1) GO TO 999 J = 1 DO 40 I = 2, P XI = X(I) IM1 = I - 1 J = J + 1 DO 30 K = 1, IM1 Y(K) = Y(K) + SS(J)*XI J = J + 1 30 CONTINUE 40 CONTINUE C 999 RETURN C *** LAST CARD OF SS7LVM FOLLOWS *** END SUBROUTINE SS7LUP(A, COSMIN, P, SIZE, STEP, U, W, WCHMTD, WSCALE, 1 Y) C C *** UPDATE SYMMETRIC A SO THAT A * STEP = Y *** C *** (LOWER TRIANGLE OF A STORED ROWWISE *** C C *** PARAMETER DECLARATIONS *** C INTEGER P REAL A(*), COSMIN, SIZE, STEP(P), U(P), W(P), 1 WCHMTD(P), WSCALE, Y(P) C DIMENSION A(P*(P+1)/2) C C *** LOCAL VARIABLES *** C INTEGER I, J, K REAL DENMIN, SDOTWM, T, UI, WI C C *** CONSTANTS *** REAL HALF, ONE, ZERO C C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C EXTERNAL SD7TPR, SS7LVM, SV2NRM REAL SD7TPR, SV2NRM PARAMETER (HALF=0.5E+0, ONE=1.E+0, ZERO=0.E+0) C ------------------------------------------------------------------ C SDOTWM = SD7TPR(P, STEP, WCHMTD) DENMIN = COSMIN * SV2NRM(P,STEP) * SV2NRM(P,WCHMTD) WSCALE = ONE IF (DENMIN .NE. ZERO) WSCALE = min(ONE, abs(SDOTWM/DENMIN)) T = ZERO IF (SDOTWM .NE. ZERO) T = WSCALE / SDOTWM DO 10 I = 1, P 10 W(I) = T * WCHMTD(I) CALL SS7LVM(P, U, A, STEP) T = HALF * (SIZE * SD7TPR(P, STEP, U) - SD7TPR(P, STEP, Y)) DO 20 I = 1, P 20 U(I) = T*W(I) + Y(I) - SIZE*U(I) C C *** SET A = A + U*(W**T) + W*(U**T) *** C K = 1 DO 40 I = 1, P UI = U(I) WI = W(I) DO 30 J = 1, I A(K) = SIZE*A(K) + UI*W(J) + WI*U(J) K = K + 1 30 CONTINUE 40 CONTINUE C C *** LAST CARD OF SS7LUP FOLLOWS *** END SUBROUTINE SL7MST(D, G, IERR, IPIVOT, KA, P, QTR, R, STEP, V, W) C C *** COMPUTE LEVENBERG-MARQUARDT STEP USING MORE-HEBDEN TECHNIQUE ** C *** NL2SOL VERSION 2.2. *** C C *** PARAMETER DECLARATIONS *** C INTEGER IERR, KA, P INTEGER IPIVOT(P) REAL D(P), G(P), QTR(P), R(*), STEP(P), V(21), W(*) C DIMENSION R(*), W(P*(P+5)/2 + 4) C C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C *** PURPOSE *** C C GIVEN THE R MATRIX FROM THE QR DECOMPOSITION OF A JACOBIAN C MATRIX, J, AS WELL AS Q-TRANSPOSE TIMES THE CORRESPONDING C RESIDUAL VECTOR, RESID, THIS SUBROUTINE COMPUTES A LEVENBERG- C MARQUARDT STEP OF APPROXIMATE LENGTH V(RADIUS) BY THE MORE- C TECHNIQUE. C C *** PARAMETER DESCRIPTION *** C C D (IN) = THE SCALE VECTOR. C G (IN) = THE GRADIENT VECTOR (J**T)*R. C IERR (I/O) = RETURN CODE FROM QRFACT OR DQ7RGS -- 0 MEANS R HAS C FULL RANK. C IPIVOT (I/O) = PERMUTATION ARRAY FROM QRFACT OR DQ7RGS, WHICH COMPUTE C QR DECOMPOSITIONS WITH COLUMN PIVOTING. C KA (I/O). KA .LT. 0 ON INPUT MEANS THIS IS THE FIRST CALL ON C SL7MST FOR THE CURRENT R AND QTR. ON OUTPUT KA CON- C TAINS THE NUMBER OF HEBDEN ITERATIONS NEEDED TO DETERMINE C STEP. KA = 0 MEANS A GAUSS-NEWTON STEP. C P (IN) = NUMBER OF PARAMETERS. C QTR (IN) = (Q**T)*RESID = Q-TRANSPOSE TIMES THE RESIDUAL VECTOR. C R (IN) = THE R MATRIX, STORED COMPACTLY BY COLUMNS. C STEP (OUT) = THE LEVENBERG-MARQUARDT STEP COMPUTED. C V (I/O) CONTAINS VARIOUS CONSTANTS AND VARIABLES DESCRIBED BELOW. C W (I/O) = WORKSPACE OF LENGTH P*(P+5)/2 + 4. C C *** ENTRIES IN V *** C C V(DGNORM) (I/O) = 2-NORM OF (D**-1)*G. C V(DSTNRM) (I/O) = 2-NORM OF D*STEP. C V(DST0) (I/O) = 2-NORM OF GAUSS-NEWTON STEP (FOR NONSING. J). C V(EPSLON) (IN) = MAX. REL. ERROR ALLOWED IN TWONORM(R)**2 MINUS C TWONORM(R - J*STEP)**2. (SEE ALGORITHM NOTES BELOW.) C V(GTSTEP) (OUT) = INNER PRODUCT BETWEEN G AND STEP. C V(NREDUC) (OUT) = HALF THE REDUCTION IN THE SUM OF SQUARES PREDICTED C FOR A GAUSS-NEWTON STEP. C V(PHMNFC) (IN) = TOL. (TOGETHER WITH V(PHMXFC)) FOR ACCEPTING STEP C (MORE*S SIGMA). THE ERROR V(DSTNRM) - V(RADIUS) MUST LIE C BETWEEN V(PHMNFC)*V(RADIUS) AND V(PHMXFC)*V(RADIUS). C V(PHMXFC) (IN) (SEE V(PHMNFC).) C V(PREDUC) (OUT) = HALF THE REDUCTION IN THE SUM OF SQUARES PREDICTED C BY THE STEP RETURNED. C V(RADIUS) (IN) = RADIUS OF CURRENT (SCALED) TRUST REGION. C V(RAD0) (I/O) = VALUE OF V(RADIUS) FROM PREVIOUS CALL. C V(STPPAR) (I/O) = MARQUARDT PARAMETER (OR ITS NEGATIVE IF THE SPECIAL C CASE MENTIONED BELOW IN THE ALGORITHM NOTES OCCURS). C C NOTE -- SEE DATA STATEMENT BELOW FOR VALUES OF ABOVE SUBSCRIPTS. C C *** USAGE NOTES *** C C IF IT IS DESIRED TO RECOMPUTE STEP USING A DIFFERENT VALUE OF C V(RADIUS), THEN THIS ROUTINE MAY BE RESTARTED BY CALLING IT C WITH ALL PARAMETERS UNCHANGED EXCEPT V(RADIUS). (THIS EXPLAINS C WHY MANY PARAMETERS ARE LISTED AS I/O). ON AN INTIIAL CALL (ONE C WITH KA = -1), THE CALLER NEED ONLY HAVE INITIALIZED D, G, KA, P, C QTR, R, V(EPSLON), V(PHMNFC), V(PHMXFC), V(RADIUS), AND V(RAD0). C C *** APPLICATION AND USAGE RESTRICTIONS *** C C THIS ROUTINE IS CALLED AS PART OF THE NL2SOL (NONLINEAR LEAST- C SQUARES) PACKAGE (REF. 1). C C *** ALGORITHM NOTES *** C C THIS CODE IMPLEMENTS THE STEP COMPUTATION SCHEME DESCRIBED IN C REFS. 2 AND 4. FAST GIVENS TRANSFORMATIONS (SEE REF. 3, PP. 60- C 62) ARE USED TO COMPUTE STEP WITH A NONZERO MARQUARDT PARAMETER. C A SPECIAL CASE OCCURS IF J IS (NEARLY) SINGULAR AND V(RADIUS) C IS SUFFICIENTLY LARGE. IN THIS CASE THE STEP RETURNED IS SUCH C THAT TWONORM(R)**2 - TWONORM(R - J*STEP)**2 DIFFERS FROM ITS C OPTIMAL VALUE BY LESS THAN V(EPSLON) TIMES THIS OPTIMAL VALUE, C WHERE J AND R DENOTE THE ORIGINAL JACOBIAN AND RESIDUAL. (SEE C REF. 2 FOR MORE DETAILS.) C C *** FUNCTIONS AND SUBROUTINES CALLED *** C C SD7TPR - RETURNS INNER PRODUCT OF TWO VECTORS. C SL7ITV - APPLY INVERSE-TRANSPOSE OF COMPACT LOWER TRIANG. MATRIX. C SL7IVM - APPLY INVERSE OF COMPACT LOWER TRIANG. MATRIX. C SV7CPY - COPIES ONE VECTOR TO ANOTHER. C SV2NRM - RETURNS 2-NORM OF A VECTOR. C C *** REFERENCES *** C C 1. DENNIS, J.E., GAY, D.M., AND WELSCH, R.E. (1981), AN ADAPTIVE C NONLINEAR LEAST-SQUARES ALGORITHM, ACM TRANS. MATH. C SOFTWARE, VOL. 7, NO. 3. C 2. GAY, D.M. (1981), COMPUTING OPTIMAL LOCALLY CONSTRAINED STEPS, C SIAM J. SCI. STATIST. COMPUTING, VOL. 2, NO. 2, PP. C 186-197. C 3. LAWSON, C.L., AND HANSON, R.J. (1974), SOLVING LEAST SQUARES C PROBLEMS, PRENTICE-HALL, ENGLEWOOD CLIFFS, N.J. C 4. MORE, J.J. (1978), THE LEVENBERG-MARQUARDT ALGORITHM, IMPLEMEN- C TATION AND THEORY, PP.105-116 OF SPRINGER LECTURE NOTES C IN MATHEMATICS NO. 630, EDITED BY G.A. WATSON, SPRINGER- C VERLAG, BERLIN AND NEW YORK. C C *** GENERAL *** C C CODED BY DAVID M. GAY. C THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH C SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS C MCS-7600324, DCR75-10143, 76-14311DSS, MCS76-11989, AND C MCS-7906671. C C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C *** LOCAL VARIABLES *** C INTEGER DSTSAV, I, IP1, I1, J1, K, KALIM, L, LK0, PHIPIN, 1 PP1O2, RES, RES0, RMAT, RMAT0, UK0 REAL A, ADI, ALPHAK, B, DFACSQ, DST, DTOL, D1, D2, 1 LK, OLDPHI, PHI, PHIMAX, PHIMIN, PSIFAC, RAD, 2 SI, SJ, SQRTAK, T, TWOPSI, UK, WL C C *** CONSTANTS *** REAL DFAC, EIGHT, HALF, NEGONE, ONE, P001, THREE, 1 TTOL, ZERO REAL BIG C C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C EXTERNAL SD7TPR, SL7ITV, SL7IVM, SL7SVN, SR7MDC,SV7CPY, SV2NRM REAL SD7TPR, SL7SVN, SR7MDC, SV2NRM C C *** SUBSCRIPTS FOR V *** C INTEGER DGNORM, DSTNRM, DST0, EPSLON, GTSTEP, NREDUC, PHMNFC, 1 PHMXFC, PREDUC, RADIUS, RAD0, STPPAR PARAMETER (DGNORM=1, DSTNRM=2, DST0=3, EPSLON=19, GTSTEP=4, 1 NREDUC=6, PHMNFC=20, PHMXFC=21, PREDUC=7, RADIUS=8, 2 RAD0=9, STPPAR=5) PARAMETER (DFAC=256.E+0, EIGHT=8.E+0, HALF=0.5E+0, NEGONE=-1.E+0, 1 ONE=1.E+0, P001=1.E-3, THREE=3.E+0, TTOL=2.5E+0, 2 ZERO=0.E+0) SAVE BIG DATA BIG/0.E+0/ C C *** BODY *** C C *** FOR USE IN RECOMPUTING STEP, THE FINAL VALUES OF LK AND UK, C *** THE INVERSE DERIVATIVE OF MORE*S PHI AT 0 (FOR NONSING. J) C *** AND THE VALUE RETURNED AS V(DSTNRM) ARE STORED AT W(LK0), C *** W(UK0), W(PHIPIN), AND W(DSTSAV) RESPECTIVELY. LK0 = P + 1 PHIPIN = LK0 + 1 UK0 = PHIPIN + 1 DSTSAV = UK0 + 1 RMAT0 = DSTSAV C *** A COPY OF THE R-MATRIX FROM THE QR DECOMPOSITION OF J IS C *** STORED IN W STARTING AT W(RMAT), AND A COPY OF THE RESIDUAL C *** VECTOR IS STORED IN W STARTING AT W(RES). THE LOOPS BELOW C *** THAT UPDATE THE QR DECOMP. FOR A NONZERO MARQUARDT PARAMETER C *** WORK ON THESE COPIES. RMAT = RMAT0 + 1 PP1O2 = P * (P + 1) / 2 RES0 = PP1O2 + RMAT0 RES = RES0 + 1 RAD = V(RADIUS) IF (RAD .GT. ZERO) 1 PSIFAC = V(EPSLON)/((EIGHT*(V(PHMNFC) + ONE) + THREE) * RAD**2) IF (BIG .LE. ZERO) BIG = SR7MDC(6) PHIMAX = V(PHMXFC) * RAD PHIMIN = V(PHMNFC) * RAD C *** DTOL, DFAC, AND DFACSQ ARE USED IN RESCALING THE FAST GIVENS C *** REPRESENTATION OF THE UPDATED QR DECOMPOSITION. DTOL = ONE/DFAC DFACSQ = DFAC*DFAC C *** OLDPHI IS USED TO DETECT LIMITS OF NUMERICAL ACCURACY. IF C *** WE RECOMPUTE STEP AND IT DOES NOT CHANGE, THEN WE ACCEPT IT. OLDPHI = ZERO LK = ZERO UK = ZERO KALIM = KA + 12 C C *** START OR RESTART, DEPENDING ON KA *** C IF (KA) 10, 20, 370 C C *** FRESH START -- COMPUTE V(NREDUC) *** C 10 KA = 0 KALIM = 12 K = P IF (IERR .NE. 0) K = abs(IERR) - 1 V(NREDUC) = HALF*SD7TPR(K, QTR, QTR) C C *** SET UP TO TRY INITIAL GAUSS-NEWTON STEP *** C 20 V(DST0) = NEGONE IF (IERR .NE. 0) GO TO 90 T = SL7SVN(P, R, STEP, W(RES)) IF (T .GE. ONE) GO TO 30 IF (SV2NRM(P, QTR) .GE. BIG*T) GO TO 90 C C *** COMPUTE GAUSS-NEWTON STEP *** C C *** NOTE -- THE R-MATRIX IS STORED COMPACTLY BY COLUMNS IN C *** R(1), R(2), R(3), ... IT IS THE TRANSPOSE OF A C *** LOWER TRIANGULAR MATRIX STORED COMPACTLY BY ROWS, AND WE C *** TREAT IT AS SUCH WHEN USING SL7ITV AND SL7IVM. 30 CALL SL7ITV(P, W, R, QTR) C *** TEMPORARILY STORE PERMUTED -D*STEP IN STEP. DO 60 I = 1, P J1 = IPIVOT(I) STEP(I) = D(J1)*W(I) 60 CONTINUE DST = SV2NRM(P, STEP) V(DST0) = DST PHI = DST - RAD IF (PHI .LE. PHIMAX) GO TO 410 C *** IF THIS IS A RESTART, GO TO 110 *** IF (KA .GT. 0) GO TO 110 C C *** GAUSS-NEWTON STEP WAS UNACCEPTABLE. COMPUTE L0 *** C DO 70 I = 1, P J1 = IPIVOT(I) STEP(I) = D(J1)*(STEP(I)/DST) 70 CONTINUE CALL SL7IVM(P, STEP, R, STEP) T = ONE / SV2NRM(P, STEP) W(PHIPIN) = (T/RAD)*T LK = PHI*W(PHIPIN) C C *** COMPUTE U0 *** C 90 DO 100 I = 1, P 100 W(I) = G(I)/D(I) V(DGNORM) = SV2NRM(P, W) UK = V(DGNORM)/RAD IF (UK .LE. ZERO) GO TO 390 C C *** ALPHAK WILL BE USED AS THE CURRENT MARQUARDT PARAMETER. WE C *** USE MORE*S SCHEME FOR INITIALIZING IT. C ALPHAK = abs(V(STPPAR)) * V(RAD0)/RAD ALPHAK = min(UK, max(ALPHAK, LK)) C C C *** TOP OF LOOP -- INCREMENT KA, COPY R TO RMAT, QTR TO RES *** C 110 KA = KA + 1 CALL SV7CPY(PP1O2, W(RMAT), R) CALL SV7CPY(P, W(RES), QTR) C C *** SAFEGUARD ALPHAK AND INITIALIZE FAST GIVENS SCALE VECTOR. *** C IF (ALPHAK .LE. ZERO .OR. ALPHAK .LT. LK .OR. ALPHAK .GE. UK) 1 ALPHAK = UK * max(P001, sqrt(LK/UK)) IF (ALPHAK .LE. ZERO) ALPHAK = HALF * UK SQRTAK = sqrt(ALPHAK) DO 120 I = 1, P 120 W(I) = ONE C C *** ADD ALPHAK*D AND UPDATE QR DECOMP. USING FAST GIVENS TRANS. *** C DO 270 I = 1, P C *** GENERATE, APPLY 1ST GIVENS TRANS. FOR ROW I OF ALPHAK*D. C *** (USE STEP TO STORE TEMPORARY ROW) *** L = I*(I+1)/2 + RMAT0 WL = W(L) D2 = ONE D1 = W(I) J1 = IPIVOT(I) ADI = SQRTAK*D(J1) IF (ADI .GE. abs(WL)) GO TO 150 130 A = ADI/WL B = D2*A/D1 T = A*B + ONE IF (T .GT. TTOL) GO TO 150 W(I) = D1/T D2 = D2/T W(L) = T*WL A = -A DO 140 J1 = I, P L = L + J1 STEP(J1) = A*W(L) 140 CONTINUE GO TO 170 C 150 B = WL/ADI A = D1*B/D2 T = A*B + ONE IF (T .GT. TTOL) GO TO 130 W(I) = D2/T D2 = D1/T W(L) = T*ADI DO 160 J1 = I, P L = L + J1 WL = W(L) STEP(J1) = -WL W(L) = A*WL 160 CONTINUE C 170 IF (I .EQ. P) GO TO 280 C C *** NOW USE GIVENS TRANS. TO ZERO ELEMENTS OF TEMP. ROW *** C IP1 = I + 1 DO 260 I1 = IP1, P L = I1*(I1+1)/2 + RMAT0 WL = W(L) SI = STEP(I1-1) D1 = W(I1) C C *** RESCALE ROW I1 IF NECESSARY *** C IF (D1 .GE. DTOL) GO TO 190 D1 = D1*DFACSQ WL = WL/DFAC K = L DO 180 J1 = I1, P K = K + J1 W(K) = W(K)/DFAC 180 CONTINUE C C *** USE GIVENS TRANS. TO ZERO NEXT ELEMENT OF TEMP. ROW C 190 IF (abs(SI) .GT. abs(WL)) GO TO 220 IF (SI .EQ. ZERO) GO TO 260 200 A = SI/WL B = D2*A/D1 T = A*B + ONE IF (T .GT. TTOL) GO TO 220 W(L) = T*WL W(I1) = D1/T D2 = D2/T DO 210 J1 = I1, P L = L + J1 WL = W(L) SJ = STEP(J1) W(L) = WL + B*SJ STEP(J1) = SJ - A*WL 210 CONTINUE GO TO 240 C 220 B = WL/SI A = D1*B/D2 T = A*B + ONE IF (T .GT. TTOL) GO TO 200 W(I1) = D2/T D2 = D1/T W(L) = T*SI DO 230 J1 = I1, P L = L + J1 WL = W(L) SJ = STEP(J1) W(L) = A*WL + SJ STEP(J1) = B*SJ - WL 230 CONTINUE C C *** RESCALE TEMP. ROW IF NECESSARY *** C 240 IF (D2 .GE. DTOL) GO TO 260 D2 = D2*DFACSQ DO 250 K = I1, P 250 STEP(K) = STEP(K)/DFAC 260 CONTINUE 270 CONTINUE C C *** COMPUTE STEP *** C 280 CALL SL7ITV(P, W(RES), W(RMAT), W(RES)) C *** RECOVER STEP AND STORE PERMUTED -D*STEP AT W(RES) *** DO 290 I = 1, P J1 = IPIVOT(I) K = RES0 + I T = W(K) STEP(J1) = -T W(K) = T*D(J1) 290 CONTINUE DST = SV2NRM(P, W(RES)) PHI = DST - RAD IF (PHI .LE. PHIMAX .AND. PHI .GE. PHIMIN) GO TO 430 IF (OLDPHI .EQ. PHI) GO TO 430 OLDPHI = PHI C C *** CHECK FOR (AND HANDLE) SPECIAL CASE *** C IF (PHI .GT. ZERO) GO TO 310 IF (KA .GE. KALIM) GO TO 430 TWOPSI = ALPHAK*DST*DST - SD7TPR(P, STEP, G) IF (ALPHAK .GE. TWOPSI*PSIFAC) GO TO 310 V(STPPAR) = -ALPHAK GO TO 440 C C *** UNACCEPTABLE STEP -- UPDATE LK, UK, ALPHAK, AND TRY AGAIN *** C 300 IF (PHI .LT. ZERO) UK = min(UK, ALPHAK) GO TO 320 310 IF (PHI .LT. ZERO) UK = ALPHAK 320 DO 330 I = 1, P J1 = IPIVOT(I) K = RES0 + I STEP(I) = D(J1) * (W(K)/DST) 330 CONTINUE CALL SL7IVM(P, STEP, W(RMAT), STEP) DO 340 I = 1, P 340 STEP(I) = STEP(I) / sqrt(W(I)) T = ONE / SV2NRM(P, STEP) ALPHAK = ALPHAK + T*PHI*T/RAD LK = max(LK, ALPHAK) ALPHAK = LK GO TO 110 C C *** RESTART *** C 370 LK = W(LK0) UK = W(UK0) IF (V(DST0) .GT. ZERO .AND. V(DST0) - RAD .LE. PHIMAX) GO TO 20 ALPHAK = abs(V(STPPAR)) DST = W(DSTSAV) PHI = DST - RAD T = V(DGNORM)/RAD IF (RAD .GT. V(RAD0)) GO TO 380 C C *** SMALLER RADIUS *** UK = T IF (ALPHAK .LE. ZERO) LK = ZERO IF (V(DST0) .GT. ZERO) LK = max(LK, (V(DST0)-RAD)*W(PHIPIN)) GO TO 300 C C *** BIGGER RADIUS *** 380 IF (ALPHAK .LE. ZERO .OR. UK .GT. T) UK = T LK = ZERO IF (V(DST0) .GT. ZERO) LK = max(LK, (V(DST0)-RAD)*W(PHIPIN)) GO TO 300 C C *** SPECIAL CASE -- RAD .LE. 0 OR (G = 0 AND J IS SINGULAR) *** C 390 V(STPPAR) = ZERO DST = ZERO LK = ZERO UK = ZERO V(GTSTEP) = ZERO V(PREDUC) = ZERO DO 400 I = 1, P 400 STEP(I) = ZERO GO TO 450 C C *** ACCEPTABLE GAUSS-NEWTON STEP -- RECOVER STEP FROM W *** C 410 ALPHAK = ZERO DO 420 I = 1, P J1 = IPIVOT(I) STEP(J1) = -W(I) 420 CONTINUE C C *** SAVE VALUES FOR USE IN A POSSIBLE RESTART *** C 430 V(STPPAR) = ALPHAK 440 V(GTSTEP) = min(SD7TPR(P,STEP,G), ZERO) V(PREDUC) = HALF * (ALPHAK*DST*DST - V(GTSTEP)) 450 V(DSTNRM) = DST W(DSTSAV) = DST W(LK0) = LK W(UK0) = UK V(RAD0) = RAD C C C *** LAST CARD OF SL7MST FOLLOWS *** END SUBROUTINE SG7QTS(D, DIG, DIHDI, KA, L, P, STEP, V, W) C C *** COMPUTE GOLDFELD-QUANDT-TROTTER STEP BY MORE-HEBDEN TECHNIQUE *** C *** (NL2SOL VERSION 2.2), MODIFIED A LA MORE AND SORENSEN *** C C *** PARAMETER DECLARATIONS *** C INTEGER KA, P REAL D(P), DIG(P), DIHDI(*), L(*), V(21), STEP(P), 1 W(*) C DIMENSION DIHDI(P*(P+1)/2), L(P*(P+1)/2), W(4*P+7) C C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C *** PURPOSE *** C C GIVEN THE (COMPACTLY STORED) LOWER TRIANGLE OF A SCALED C HESSIAN (APPROXIMATION) AND A NONZERO SCALED GRADIENT VECTOR, C THIS SUBROUTINE COMPUTES A GOLDFELD-QUANDT-TROTTER STEP OF C APPROXIMATE LENGTH V(RADIUS) BY THE MORE-HEBDEN TECHNIQUE. IN C OTHER WORDS, STEP IS COMPUTED TO (APPROXIMATELY) MINIMIZE C PSI(STEP) = (G**T)*STEP + 0.5*(STEP**T)*H*STEP SUCH THAT THE C 2-NORM OF D*STEP IS AT MOST (APPROXIMATELY) V(RADIUS), WHERE C G IS THE GRADIENT, H IS THE HESSIAN, AND D IS A DIAGONAL C SCALE MATRIX WHOSE DIAGONAL IS STORED IN THE PARAMETER D. C (SG7QTS ASSUMES DIG = D**-1 * G AND DIHDI = D**-1 * H * D**-1.) C C *** PARAMETER DESCRIPTION *** C C D (IN) = THE SCALE VECTOR, I.E. THE DIAGONAL OF THE SCALE C MATRIX D MENTIONED ABOVE UNDER PURPOSE. C DIG (IN) = THE SCALED GRADIENT VECTOR, D**-1 * G. IF G = 0, THEN C STEP = 0 AND V(STPPAR) = 0 ARE RETURNED. C DIHDI (IN) = LOWER TRIANGLE OF THE SCALED HESSIAN (APPROXIMATION), C I.E., D**-1 * H * D**-1, STORED COMPACTLY BY ROWS., I.E., C IN THE ORDER (1,1), (2,1), (2,2), (3,1), (3,2), ETC. C KA (I/O) = THE NUMBER OF HEBDEN ITERATIONS (SO FAR) TAKEN TO DETER- C MINE STEP. KA .LT. 0 ON INPUT MEANS THIS IS THE FIRST C ATTEMPT TO DETERMINE STEP (FOR THE PRESENT DIG AND DIHDI) C -- KA IS INITIALIZED TO 0 IN THIS CASE. OUTPUT WITH C KA = 0 (OR V(STPPAR) = 0) MEANS STEP = -(H**-1)*G. C L (I/O) = WORKSPACE OF LENGTH P*(P+1)/2 FOR CHOLESKY FACTORS. C P (IN) = NUMBER OF PARAMETERS -- THE HESSIAN IS A P X P MATRIX. C STEP (I/O) = THE STEP COMPUTED. C V (I/O) CONTAINS VARIOUS CONSTANTS AND VARIABLES DESCRIBED BELOW. C W (I/O) = WORKSPACE OF LENGTH 4*P + 6. C C *** ENTRIES IN V *** C C V(DGNORM) (I/O) = 2-NORM OF (D**-1)*G. C V(DSTNRM) (OUTPUT) = 2-NORM OF D*STEP. C V(DST0) (I/O) = 2-NORM OF D*(H**-1)*G (FOR POS. DEF. H ONLY), OR C OVERESTIMATE OF SMALLEST EIGENVALUE OF (D**-1)*H*(D**-1). C V(EPSLON) (IN) = MAX. REL. ERROR ALLOWED FOR PSI(STEP). FOR THE C STEP RETURNED, PSI(STEP) WILL EXCEED ITS OPTIMAL VALUE C BY LESS THAN -V(EPSLON)*PSI(STEP). SUGGESTED VALUE = 0.1. C V(GTSTEP) (OUT) = INNER PRODUCT BETWEEN G AND STEP. C V(NREDUC) (OUT) = PSI(-(H**-1)*G) = PSI(NEWTON STEP) (FOR POS. DEF. C H ONLY -- V(NREDUC) IS SET TO ZERO OTHERWISE). C V(PHMNFC) (IN) = TOL. (TOGETHER WITH V(PHMXFC)) FOR ACCEPTING STEP C (MORE*S SIGMA). THE ERROR V(DSTNRM) - V(RADIUS) MUST LIE C BETWEEN V(PHMNFC)*V(RADIUS) AND V(PHMXFC)*V(RADIUS). C V(PHMXFC) (IN) (SEE V(PHMNFC).) C SUGGESTED VALUES -- V(PHMNFC) = -0.25, V(PHMXFC) = 0.5. C V(PREDUC) (OUT) = PSI(STEP) = PREDICTED OBJ. FUNC. REDUCTION FOR STEP. C V(RADIUS) (IN) = RADIUS OF CURRENT (SCALED) TRUST REGION. C V(RAD0) (I/O) = VALUE OF V(RADIUS) FROM PREVIOUS CALL. C V(STPPAR) (I/O) IS NORMALLY THE MARQUARDT PARAMETER, I.E. THE ALPHA C DESCRIBED BELOW UNDER ALGORITHM NOTES. IF H + ALPHA*D**2 C (SEE ALGORITHM NOTES) IS (NEARLY) SINGULAR, HOWEVER, C THEN V(STPPAR) = -ALPHA. C C *** USAGE NOTES *** C C IF IT IS DESIRED TO RECOMPUTE STEP USING A DIFFERENT VALUE OF C V(RADIUS), THEN THIS ROUTINE MAY BE RESTARTED BY CALLING IT C WITH ALL PARAMETERS UNCHANGED EXCEPT V(RADIUS). (THIS EXPLAINS C WHY STEP AND W ARE LISTED AS I/O). ON AN INITIAL CALL (ONE WITH C KA .LT. 0), STEP AND W NEED NOT BE INITIALIZED AND ONLY COMPO- C NENTS V(EPSLON), V(STPPAR), V(PHMNFC), V(PHMXFC), V(RADIUS), AND C V(RAD0) OF V MUST BE INITIALIZED. C C *** ALGORITHM NOTES *** C C THE DESIRED G-Q-T STEP (REF. 2, 3, 4, 6) SATISFIES C (H + ALPHA*D**2)*STEP = -G FOR SOME NONNEGATIVE ALPHA SUCH THAT C H + ALPHA*D**2 IS POSITIVE SEMIDEFINITE. ALPHA AND STEP ARE C COMPUTED BY A SCHEME ANALOGOUS TO THE ONE DESCRIBED IN REF. 5. C ESTIMATES OF THE SMALLEST AND LARGEST EIGENVALUES OF THE HESSIAN C ARE OBTAINED FROM THE GERSCHGORIN CIRCLE THEOREM ENHANCED BY A C SIMPLE FORM OF THE SCALING DESCRIBED IN REF. 7. CASES IN WHICH C H + ALPHA*D**2 IS NEARLY (OR EXACTLY) SINGULAR ARE HANDLED BY C THE TECHNIQUE DISCUSSED IN REF. 2. IN THESE CASES, A STEP OF C (EXACT) LENGTH V(RADIUS) IS RETURNED FOR WHICH PSI(STEP) EXCEEDS C ITS OPTIMAL VALUE BY LESS THAN -V(EPSLON)*PSI(STEP). THE TEST C SUGGESTED IN REF. 6 FOR DETECTING THE SPECIAL CASE IS PERFORMED C ONCE TWO MATRIX FACTORIZATIONS HAVE BEEN DONE -- DOING SO SOONER C SEEMS TO DEGRADE THE PERFORMANCE OF OPTIMIZATION ROUTINES THAT C CALL THIS ROUTINE. C C *** FUNCTIONS AND SUBROUTINES CALLED *** C C SD7TPR - RETURNS INNER PRODUCT OF TWO VECTORS. C SL7ITV - APPLIES INVERSE-TRANSPOSE OF COMPACT LOWER TRIANG. MATRIX. C SL7IVM - APPLIES INVERSE OF COMPACT LOWER TRIANG. MATRIX. C SL7SRT - FINDS CHOLESKY FACTOR (OF COMPACTLY STORED LOWER TRIANG.). C SL7SVN - RETURNS APPROX. TO MIN. SING. VALUE OF LOWER TRIANG. MATRIX. C SR7MDC - RETURNS MACHINE-DEPENDENT CONSTANTS. C SV2NRM - RETURNS 2-NORM OF A VECTOR. C C *** REFERENCES *** C C 1. DENNIS, J.E., GAY, D.M., AND WELSCH, R.E. (1981), AN ADAPTIVE C NONLINEAR LEAST-SQUARES ALGORITHM, ACM TRANS. MATH. C SOFTWARE, VOL. 7, NO. 3. C 2. GAY, D.M. (1981), COMPUTING OPTIMAL LOCALLY CONSTRAINED STEPS, C SIAM J. SCI. STATIST. COMPUTING, VOL. 2, NO. 2, PP. C 186-197. C 3. GOLDFELD, S.M., QUANDT, R.E., AND TROTTER, H.F. (1966), C MAXIMIZATION BY QUADRATIC HILL-CLIMBING, ECONOMETRICA 34, C PP. 541-551. C 4. HEBDEN, M.D. (1973), AN ALGORITHM FOR MINIMIZATION USING EXACT C SECOND DERIVATIVES, REPORT T.P. 515, THEORETICAL PHYSICS C DIV., A.E.R.E. HARWELL, OXON., ENGLAND. C 5. MORE, J.J. (1978), THE LEVENBERG-MARQUARDT ALGORITHM, IMPLEMEN- C TATION AND THEORY, PP.105-116 OF SPRINGER LECTURE NOTES C IN MATHEMATICS NO. 630, EDITED BY G.A. WATSON, SPRINGER- C VERLAG, BERLIN AND NEW YORK. C 6. MORE, J.J., AND SORENSEN, D.C. (1981), COMPUTING A TRUST REGION C STEP, TECHNICAL REPORT ANL-81-83, ARGONNE NATIONAL LAB. C 7. VARGA, R.S. (1965), MINIMAL GERSCHGORIN SETS, PACIFIC J. MATH. 15, C PP. 719-729. C C *** GENERAL *** C C CODED BY DAVID M. GAY. C THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH C SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS C MCS-7600324, DCR75-10143, 76-14311DSS, MCS76-11989, AND C MCS-7906671. C C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C *** LOCAL VARIABLES *** C LOGICAL RESTRT INTEGER DGGDMX, DIAG, DIAG0, DSTSAV, EMAX, EMIN, I, IM1, INC, IRC, 1 J, K, KALIM, KAMIN, K1, LK0, PHIPIN, Q, Q0, UK0, X REAL ALPHAK, AKI, AKK, DELTA, DST, EPS, GTSTA, LK, 1 OLDPHI, PHI, PHIMAX, PHIMIN, PSIFAC, RAD, RADSQ, 2 ROOT, SI, SK, SW, T, TWOPSI, T1, T2, UK, WI C C *** CONSTANTS *** REAL BIG, DGXFAC, EPSFAC, FOUR, HALF, KAPPA, NEGONE, 1 ONE, P001, SIX, THREE, TWO, ZERO C C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C EXTERNAL SD7TPR, SL7ITV, SL7IVM,SL7SRT, SL7SVN, SR7MDC, SV2NRM REAL SD7TPR, SL7SVN, SR7MDC, SV2NRM C C *** SUBSCRIPTS FOR V *** C INTEGER DGNORM, DSTNRM, DST0, EPSLON, GTSTEP, STPPAR, NREDUC, 1 PHMNFC, PHMXFC, PREDUC, RADIUS, RAD0 PARAMETER (DGNORM=1, DSTNRM=2, DST0=3, EPSLON=19, GTSTEP=4, 1 NREDUC=6, PHMNFC=20, PHMXFC=21, PREDUC=7, RADIUS=8, 2 RAD0=9, STPPAR=5) PARAMETER (EPSFAC=50.0E+0, FOUR=4.0E+0, HALF=0.5E+0, 1 KAPPA=2.0E+0, NEGONE=-1.0E+0, ONE=1.0E+0, P001=1.0E-3, 2 SIX=6.0E+0, THREE=3.0E+0, TWO=2.0E+0, ZERO=0.0E+0) SAVE DGXFAC DATA BIG/0.E+0/, DGXFAC/0.E+0/ C C *** BODY *** C IF (BIG .LE. ZERO) BIG = SR7MDC(6) C C *** STORE LARGEST ABS. ENTRY IN (D**-1)*H*(D**-1) AT W(DGGDMX). DGGDMX = P + 1 C *** STORE GERSCHGORIN OVER- AND UNDERESTIMATES OF THE LARGEST C *** AND SMALLEST EIGENVALUES OF (D**-1)*H*(D**-1) AT W(EMAX) C *** AND W(EMIN) RESPECTIVELY. EMAX = DGGDMX + 1 EMIN = EMAX + 1 C *** FOR USE IN RECOMPUTING STEP, THE FINAL VALUES OF LK, UK, DST, C *** AND THE INVERSE DERIVATIVE OF MORE*S PHI AT 0 (FOR POS. DEF. C *** H) ARE STORED IN W(LK0), W(UK0), W(DSTSAV), AND W(PHIPIN) C *** RESPECTIVELY. LK0 = EMIN + 1 PHIPIN = LK0 + 1 UK0 = PHIPIN + 1 DSTSAV = UK0 + 1 C *** STORE DIAG OF (D**-1)*H*(D**-1) IN W(DIAG),...,W(DIAG0+P). DIAG0 = DSTSAV DIAG = DIAG0 + 1 C *** STORE -D*STEP IN W(Q),...,W(Q0+P). Q0 = DIAG0 + P Q = Q0 + 1 C *** ALLOCATE STORAGE FOR SCRATCH VECTOR X *** X = Q + P RAD = V(RADIUS) RADSQ = RAD**2 C *** PHITOL = MAX. ERROR ALLOWED IN DST = V(DSTNRM) = 2-NORM OF C *** D*STEP. PHIMAX = V(PHMXFC) * RAD PHIMIN = V(PHMNFC) * RAD PSIFAC = BIG T1 = TWO * V(EPSLON) / (THREE * (FOUR * (V(PHMNFC) + ONE) * 1 (KAPPA + ONE) + KAPPA + TWO) * RAD) IF (T1 .LT. BIG*min(RAD,ONE)) PSIFAC = T1 / RAD C *** OLDPHI IS USED TO DETECT LIMITS OF NUMERICAL ACCURACY. IF C *** WE RECOMPUTE STEP AND IT DOES NOT CHANGE, THEN WE ACCEPT IT. OLDPHI = ZERO EPS = V(EPSLON) IRC = 0 RESTRT = .FALSE. KALIM = KA + 50 C C *** START OR RESTART, DEPENDING ON KA *** C IF (KA .GE. 0) GO TO 290 C C *** FRESH START *** C K = 0 UK = NEGONE KA = 0 KALIM = 50 V(DGNORM) = SV2NRM(P, DIG) V(NREDUC) = ZERO V(DST0) = ZERO KAMIN = 3 IF (V(DGNORM) .EQ. ZERO) KAMIN = 0 C C *** STORE DIAG(DIHDI) IN W(DIAG0+1),...,W(DIAG0+P) *** C J = 0 DO 10 I = 1, P J = J + I K1 = DIAG0 + I W(K1) = DIHDI(J) 10 CONTINUE C C *** DETERMINE W(DGGDMX), THE LARGEST ELEMENT OF DIHDI *** C T1 = ZERO J = P * (P + 1) / 2 DO 20 I = 1, J T = abs(DIHDI(I)) IF (T1 .LT. T) T1 = T 20 CONTINUE W(DGGDMX) = T1 C C *** TRY ALPHA = 0 *** C 30 CALL SL7SRT(1, P, L, DIHDI, IRC) IF (IRC .EQ. 0) GO TO 50 C *** INDEF. H -- UNDERESTIMATE SMALLEST EIGENVALUE, USE THIS C *** ESTIMATE TO INITIALIZE LOWER BOUND LK ON ALPHA. J = IRC*(IRC+1)/2 T = L(J) L(J) = ONE DO 40 I = 1, IRC 40 W(I) = ZERO W(IRC) = ONE CALL SL7ITV(IRC, W, L, W) T1 = SV2NRM(IRC, W) LK = -T / T1 / T1 V(DST0) = -LK IF (RESTRT) GO TO 210 GO TO 70 C C *** POSITIVE DEFINITE H -- COMPUTE UNMODIFIED NEWTON STEP. *** 50 LK = ZERO T = SL7SVN(P, L, W(Q), W(Q)) IF (T .GE. ONE) GO TO 60 IF (V(DGNORM) .GE. T*T*BIG) GO TO 70 60 CALL SL7IVM(P, W(Q), L, DIG) GTSTA = SD7TPR(P, W(Q), W(Q)) V(NREDUC) = HALF * GTSTA CALL SL7ITV(P, W(Q), L, W(Q)) DST = SV2NRM(P, W(Q)) V(DST0) = DST PHI = DST - RAD IF (PHI .LE. PHIMAX) GO TO 260 IF (RESTRT) GO TO 210 C C *** PREPARE TO COMPUTE GERSCHGORIN ESTIMATES OF LARGEST (AND C *** SMALLEST) EIGENVALUES. *** C 70 K = 0 DO 100 I = 1, P WI = ZERO IF (I .EQ. 1) GO TO 90 IM1 = I - 1 DO 80 J = 1, IM1 K = K + 1 T = abs(DIHDI(K)) WI = WI + T W(J) = W(J) + T 80 CONTINUE 90 W(I) = WI K = K + 1 100 CONTINUE C C *** (UNDER-)ESTIMATE SMALLEST EIGENVALUE OF (D**-1)*H*(D**-1) *** C K = 1 T1 = W(DIAG) - W(1) IF (P .LE. 1) GO TO 120 DO 110 I = 2, P J = DIAG0 + I T = W(J) - W(I) IF (T .GE. T1) GO TO 110 T1 = T K = I 110 CONTINUE C 120 SK = W(K) J = DIAG0 + K AKK = W(J) K1 = K*(K-1)/2 + 1 INC = 1 T = ZERO DO 150 I = 1, P IF (I .EQ. K) GO TO 130 AKI = abs(DIHDI(K1)) SI = W(I) J = DIAG0 + I T1 = HALF * (AKK - W(J) + SI - AKI) T1 = T1 + sqrt(T1*T1 + SK*AKI) IF (T .LT. T1) T = T1 IF (I .LT. K) GO TO 140 130 INC = I 140 K1 = K1 + INC 150 CONTINUE C W(EMIN) = AKK - T UK = V(DGNORM)/RAD - W(EMIN) IF (V(DGNORM) .EQ. ZERO) UK = UK + P001 + P001*UK IF (UK .LE. ZERO) UK = P001 C C *** COMPUTE GERSCHGORIN (OVER-)ESTIMATE OF LARGEST EIGENVALUE *** C K = 1 T1 = W(DIAG) + W(1) IF (P .LE. 1) GO TO 170 DO 160 I = 2, P J = DIAG0 + I T = W(J) + W(I) IF (T .LE. T1) GO TO 160 T1 = T K = I 160 CONTINUE C 170 SK = W(K) J = DIAG0 + K AKK = W(J) K1 = K*(K-1)/2 + 1 INC = 1 T = ZERO DO 200 I = 1, P IF (I .EQ. K) GO TO 180 AKI = abs(DIHDI(K1)) SI = W(I) J = DIAG0 + I T1 = HALF * (W(J) + SI - AKI - AKK) T1 = T1 + sqrt(T1*T1 + SK*AKI) IF (T .LT. T1) T = T1 IF (I .LT. K) GO TO 190 180 INC = I 190 K1 = K1 + INC 200 CONTINUE C W(EMAX) = AKK + T LK = max(LK, V(DGNORM)/RAD - W(EMAX)) C C *** ALPHAK = CURRENT VALUE OF ALPHA (SEE ALG. NOTES ABOVE). WE C *** USE MORE*S SCHEME FOR INITIALIZING IT. ALPHAK = abs(V(STPPAR)) * V(RAD0)/RAD ALPHAK = min(UK, max(ALPHAK, LK)) C IF (IRC .NE. 0) GO TO 210 C C *** COMPUTE L0 FOR POSITIVE DEFINITE H *** C CALL SL7IVM(P, W, L, W(Q)) T = SV2NRM(P, W) W(PHIPIN) = RAD / T / T LK = max(LK, PHI*W(PHIPIN)) C C *** SAFEGUARD ALPHAK AND ADD ALPHAK*I TO (D**-1)*H*(D**-1) *** C 210 KA = KA + 1 IF (-V(DST0) .GE. ALPHAK .OR. ALPHAK .LT. LK .OR. ALPHAK .GE. UK) 1 ALPHAK = UK * max(P001, sqrt(LK/UK)) IF (ALPHAK .LE. ZERO) ALPHAK = HALF * UK IF (ALPHAK .LE. ZERO) ALPHAK = UK K = 0 DO 220 I = 1, P K = K + I J = DIAG0 + I DIHDI(K) = W(J) + ALPHAK 220 CONTINUE C C *** TRY COMPUTING CHOLESKY DECOMPOSITION *** C CALL SL7SRT(1, P, L, DIHDI, IRC) IF (IRC .EQ. 0) GO TO 240 C C *** (D**-1)*H*(D**-1) + ALPHAK*I IS INDEFINITE -- OVERESTIMATE C *** SMALLEST EIGENVALUE FOR USE IN UPDATING LK *** C J = (IRC*(IRC+1))/2 T = L(J) L(J) = ONE DO 230 I = 1, IRC 230 W(I) = ZERO W(IRC) = ONE CALL SL7ITV(IRC, W, L, W) T1 = SV2NRM(IRC, W) LK = ALPHAK - T/T1/T1 V(DST0) = -LK IF (ALPHAK .LT. LK) GO TO 210 C C *** NASTY CASE -- EXACT GERSCHGORIN BOUNDS. FUDGE LK, UK... C T = P001 * ALPHAK IF (T .LE. ZERO) T = P001 LK = ALPHAK + T IF (UK .LE. LK) UK = LK + T GO TO 210 C C *** ALPHAK MAKES (D**-1)*H*(D**-1) POSITIVE DEFINITE. C *** COMPUTE Q = -D*STEP, CHECK FOR CONVERGENCE. *** C 240 CALL SL7IVM(P, W(Q), L, DIG) GTSTA = SD7TPR(P, W(Q), W(Q)) CALL SL7ITV(P, W(Q), L, W(Q)) DST = SV2NRM(P, W(Q)) PHI = DST - RAD IF (PHI .LE. PHIMAX .AND. PHI .GE. PHIMIN) GO TO 270 IF (PHI .EQ. OLDPHI) GO TO 270 OLDPHI = PHI IF (PHI .LT. ZERO) GO TO 330 C C *** UNACCEPTABLE ALPHAK -- UPDATE LK, UK, ALPHAK *** C 250 IF (KA .GE. KALIM) GO TO 270 C *** THE FOLLOWING min IS NECESSARY BECAUSE OF RESTARTS *** IF (PHI .LT. ZERO) UK = min(UK, ALPHAK) C *** KAMIN = 0 ONLY IFF THE GRADIENT VANISHES *** IF (KAMIN .EQ. 0) GO TO 210 CALL SL7IVM(P, W, L, W(Q)) C *** THE FOLLOWING, COMMENTED CALCULATION OF ALPHAK IS SOMETIMES C *** SAFER BUT WORSE IN PERFORMANCE... C T1 = DST / SV2NRM(P, W) C ALPHAK = ALPHAK + T1 * (PHI/RAD) * T1 T1 = SV2NRM(P, W) ALPHAK = ALPHAK + (PHI/T1) * (DST/T1) * (DST/RAD) LK = max(LK, ALPHAK) ALPHAK = LK GO TO 210 C C *** ACCEPTABLE STEP ON FIRST TRY *** C 260 ALPHAK = ZERO C C *** SUCCESSFUL STEP IN GENERAL. COMPUTE STEP = -(D**-1)*Q *** C 270 DO 280 I = 1, P J = Q0 + I STEP(I) = -W(J)/D(I) 280 CONTINUE V(GTSTEP) = -GTSTA V(PREDUC) = HALF * (abs(ALPHAK)*DST*DST + GTSTA) GO TO 410 C C C *** RESTART WITH NEW RADIUS *** C 290 IF (V(DST0) .LE. ZERO .OR. V(DST0) - RAD .GT. PHIMAX) GO TO 310 C C *** PREPARE TO RETURN NEWTON STEP *** C RESTRT = .TRUE. KA = KA + 1 K = 0 DO 300 I = 1, P K = K + I J = DIAG0 + I DIHDI(K) = W(J) 300 CONTINUE UK = NEGONE GO TO 30 C 310 KAMIN = KA + 3 IF (V(DGNORM) .EQ. ZERO) KAMIN = 0 IF (KA .EQ. 0) GO TO 50 C DST = W(DSTSAV) ALPHAK = abs(V(STPPAR)) PHI = DST - RAD T = V(DGNORM)/RAD UK = T - W(EMIN) IF (V(DGNORM) .EQ. ZERO) UK = UK + P001 + P001*UK IF (UK .LE. ZERO) UK = P001 IF (RAD .GT. V(RAD0)) GO TO 320 C C *** SMALLER RADIUS *** LK = ZERO IF (ALPHAK .GT. ZERO) LK = W(LK0) LK = max(LK, T - W(EMAX)) IF (V(DST0) .GT. ZERO) LK = max(LK, (V(DST0)-RAD)*W(PHIPIN)) GO TO 250 C C *** BIGGER RADIUS *** 320 IF (ALPHAK .GT. ZERO) UK = min(UK, W(UK0)) LK = max(max(ZERO, -V(DST0)), T - W(EMAX)) IF (V(DST0) .GT. ZERO) LK = max(LK, (V(DST0)-RAD)*W(PHIPIN)) GO TO 250 C C *** DECIDE WHETHER TO CHECK FOR SPECIAL CASE... IN PRACTICE (FROM C *** THE STANDPOINT OF THE CALLING OPTIMIZATION CODE) IT SEEMS BEST C *** NOT TO CHECK UNTIL A FEW ITERATIONS HAVE FAILED -- HENCE THE C *** TEST ON KAMIN BELOW. C 330 DELTA = ALPHAK + min(ZERO, V(DST0)) TWOPSI = ALPHAK*DST*DST + GTSTA IF (KA .GE. KAMIN) GO TO 340 C *** IF THE TEST IN REF. 2 IS SATISFIED, FALL THROUGH TO HANDLE C *** THE SPECIAL CASE (AS SOON AS THE MORE-SORENSEN TEST DETECTS C *** IT). IF (PSIFAC .GE. BIG) GO TO 340 IF (DELTA .GE. PSIFAC*TWOPSI) GO TO 370 C C *** CHECK FOR THE SPECIAL CASE OF H + ALPHA*D**2 (NEARLY) C *** SINGULAR. USE ONE STEP OF INVERSE POWER METHOD WITH START C *** FROM SL7SVN TO OBTAIN APPROXIMATE EIGENVECTOR CORRESPONDING C *** TO SMALLEST EIGENVALUE OF (D**-1)*H*(D**-1). SL7SVN RETURNS C *** X AND W WITH L*W = X. C 340 T = SL7SVN(P, L, W(X), W) C C *** NORMALIZE W *** DO 350 I = 1, P 350 W(I) = T*W(I) C *** COMPLETE CURRENT INV. POWER ITER. -- REPLACE W BY (L**-T)*W. CALL SL7ITV(P, W, L, W) T2 = ONE/SV2NRM(P, W) DO 360 I = 1, P 360 W(I) = T2*W(I) T = T2 * T C C *** NOW W IS THE DESIRED APPROXIMATE (UNIT) EIGENVECTOR AND C *** T*X = ((D**-1)*H*(D**-1) + ALPHAK*I)*W. C SW = SD7TPR(P, W(Q), W) T1 = (RAD + DST) * (RAD - DST) ROOT = sqrt(SW*SW + T1) IF (SW .LT. ZERO) ROOT = -ROOT SI = T1 / (SW + ROOT) C C *** THE ACTUAL TEST FOR THE SPECIAL CASE... C IF ((T2*SI)**2 .LE. EPS*(DST**2 + ALPHAK*RADSQ)) GO TO 380 C C *** UPDATE UPPER BOUND ON SMALLEST EIGENVALUE (WHEN NOT POSITIVE) C *** (AS RECOMMENDED BY MORE AND SORENSEN) AND CONTINUE... C IF (V(DST0) .LE. ZERO) V(DST0) = min(V(DST0), T2**2 - ALPHAK) LK = max(LK, -V(DST0)) C C *** CHECK WHETHER WE CAN HOPE TO DETECT THE SPECIAL CASE IN C *** THE AVAILABLE ARITHMETIC. ACCEPT STEP AS IT IS IF NOT. C C *** IF NOT YET AVAILABLE, OBTAIN MACHINE DEPENDENT VALUE DGXFAC. 370 IF (DGXFAC .EQ. ZERO) DGXFAC = EPSFAC * SR7MDC(3) C IF (DELTA .GT. DGXFAC*W(DGGDMX)) GO TO 250 GO TO 270 C C *** SPECIAL CASE DETECTED... NEGATE ALPHAK TO INDICATE SPECIAL CASE C 380 ALPHAK = -ALPHAK V(PREDUC) = HALF * TWOPSI C C *** ACCEPT CURRENT STEP IF ADDING SI*W WOULD LEAD TO A C *** FURTHER RELATIVE REDUCTION IN PSI OF LESS THAN V(EPSLON)/3. C T1 = ZERO T = SI*(ALPHAK*SW - HALF*SI*(ALPHAK + T*SD7TPR(P,W(X),W))) IF (T .LT. EPS*TWOPSI/SIX) GO TO 390 V(PREDUC) = V(PREDUC) + T DST = RAD T1 = -SI 390 DO 400 I = 1, P J = Q0 + I W(J) = T1*W(I) - W(J) STEP(I) = W(J) / D(I) 400 CONTINUE V(GTSTEP) = SD7TPR(P, DIG, W(Q)) C C *** SAVE VALUES FOR USE IN A POSSIBLE RESTART *** C 410 V(DSTNRM) = DST V(STPPAR) = ALPHAK W(LK0) = LK W(UK0) = UK V(RAD0) = RAD W(DSTSAV) = DST C C *** RESTORE DIAGONAL OF DIHDI *** C J = 0 DO 420 I = 1, P J = J + I K = DIAG0 + I DIHDI(J) = W(K) 420 CONTINUE C C C *** LAST CARD OF SG7QTS FOLLOWS *** END SUBROUTINE SL7IVM(N, X, L, Y) C C *** SOLVE L*X = Y, WHERE L IS AN N X N LOWER TRIANGULAR C *** MATRIX STORED COMPACTLY BY ROWS. X AND Y MAY OCCUPY THE SAME C *** STORAGE. *** C INTEGER N REAL X(N), L(*), Y(N) EXTERNAL SD7TPR REAL SD7TPR INTEGER I, J, K REAL T, ZERO PARAMETER (ZERO=0.E+0) C DO 10 K = 1, N IF (Y(K) .NE. ZERO) GO TO 20 X(K) = ZERO 10 CONTINUE GO TO 999 20 J = K*(K+1)/2 X(K) = Y(K) / L(J) IF (K .GE. N) GO TO 999 K = K + 1 DO 30 I = K, N T = SD7TPR(I-1, L(J+1), X) J = J + I X(I) = (Y(I) - T)/L(J) 30 CONTINUE 999 RETURN C *** LAST CARD OF SL7IVM FOLLOWS *** END REAL FUNCTION SL7SVX(P, L, X, Y) C C *** ESTIMATE LARGEST SING. VALUE OF PACKED LOWER TRIANG. MATRIX L C C *** PARAMETER DECLARATIONS *** C INTEGER P REAL L(*), X(P), Y(P) C DIMENSION L(P*(P+1)/2) C C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C *** PURPOSE *** C C THIS FUNCTION RETURNS A GOOD UNDER-ESTIMATE OF THE LARGEST C SINGULAR VALUE OF THE PACKED LOWER TRIANGULAR MATRIX L. C C *** PARAMETER DESCRIPTION *** C C P (IN) = THE ORDER OF L. L IS A P X P LOWER TRIANGULAR MATRIX. C L (IN) = ARRAY HOLDING THE ELEMENTS OF L IN ROW ORDER, I.E. C L(1,1), L(2,1), L(2,2), L(3,1), L(3,2), L(3,3), ETC. C X (OUT) IF SL7SVX RETURNS A POSITIVE VALUE, THEN X = (L**T)*Y IS AN C (UNNORMALIZED) APPROXIMATE RIGHT SINGULAR VECTOR C CORRESPONDING TO THE LARGEST SINGULAR VALUE. THIS C APPROXIMATION MAY BE CRUDE. C Y (OUT) IF SL7SVX RETURNS A POSITIVE VALUE, THEN Y = L*X IS A C NORMALIZED APPROXIMATE LEFT SINGULAR VECTOR CORRESPOND- C ING TO THE LARGEST SINGULAR VALUE. THIS APPROXIMATION C MAY BE VERY CRUDE. THE CALLER MAY PASS THE SAME VECTOR C FOR X AND Y (NONSTANDARD FORTRAN USAGE), IN WHICH CASE X C OVER-WRITES Y. C C *** ALGORITHM NOTES *** C C THE ALGORITHM IS BASED ON ANALOGY WITH (1). IT USES A C RANDOM NUMBER GENERATOR PROPOSED IN (4), WHICH PASSES THE C SPECTRAL TEST WITH FLYING COLORS -- SEE (2) AND (3). C C *** SUBROUTINES AND FUNCTIONS CALLED *** C C SV2NRM - FUNCTION, RETURNS THE 2-NORM OF A VECTOR. C C *** REFERENCES *** C C (1) CLINE, A., MOLER, C., STEWART, G., AND WILKINSON, J.H.(1977), C AN ESTIMATE FOR THE CONDITION NUMBER OF A MATRIX, REPORT C TM-310, APPLIED MATH. DIV., ARGONNE NATIONAL LABORATORY. C C (2) HOAGLIN, D.C. (1976), THEORETICAL PROPERTIES OF CONGRUENTIAL C RANDOM-NUMBER GENERATORS -- AN EMPIRICAL VIEW, C MEMORANDUM NS-340, DEPT. OF STATISTICS, HARVARD UNIV. C C (3) KNUTH, D.E. (1969), THE ART OF COMPUTER PROGRAMMING, VOL. 2 C (SEMINUMERICAL ALGORITHMS), ADDISON-WESLEY, READING, MASS. C C (4) SMITH, C.S. (1971), MULTIPLICATIVE PSEUDO-RANDOM NUMBER C GENERATORS WITH PRIME MODULUS, J. ASSOC. COMPUT. MACH. 18, C PP. 586-593. C C *** HISTORY *** C C DESIGNED AND CODED BY DAVID M. GAY (WINTER 1977/SUMMER 1978). C C *** GENERAL *** C C THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH C SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS C MCS-7600324, DCR75-10143, 76-14311DSS, AND MCS76-11989. C C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C *** LOCAL VARIABLES *** C INTEGER I, IX, J, JI, JJ, JJJ, JM1, J0, PM1, PPLUS1 REAL B, BLJI, SMINUS, SPLUS, T, YI C C *** CONSTANTS *** C REAL HALF, ONE, R9973, ZERO C C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C EXTERNAL SD7TPR, SV2NRM,SV2AXY REAL SD7TPR, SV2NRM PARAMETER (HALF=0.5E+0, ONE=1.E+0, R9973=9973.E+0, ZERO=0.E+0) C C *** BODY *** C IX = 2 PPLUS1 = P + 1 PM1 = P - 1 C C *** FIRST INITIALIZE X TO PARTIAL SUMS *** C J0 = P*PM1/2 JJ = J0 + P IX = MOD(3432*IX, 9973) B = HALF*(ONE + real(IX)/R9973) X(P) = B * L(JJ) IF (P .LE. 1) GO TO 40 DO 10 I = 1, PM1 JI = J0 + I X(I) = B * L(JI) 10 CONTINUE C C *** COMPUTE X = (L**T)*B, WHERE THE COMPONENTS OF B HAVE RANDOMLY C *** CHOSEN MAGNITUDES IN (.5,1) WITH SIGNS CHOSEN TO MAKE X LARGE. C C DO J = P-1 TO 1 BY -1... DO 30 JJJ = 1, PM1 J = P - JJJ C *** DETERMINE X(J) IN THIS ITERATION. NOTE FOR I = 1,2,...,J C *** THAT X(I) HOLDS THE CURRENT PARTIAL SUM FOR ROW I. IX = MOD(3432*IX, 9973) B = HALF*(ONE + real(IX)/R9973) JM1 = J - 1 J0 = J*JM1/2 SPLUS = ZERO SMINUS = ZERO DO 20 I = 1, J JI = J0 + I BLJI = B * L(JI) SPLUS = SPLUS + abs(BLJI + X(I)) SMINUS = SMINUS + abs(BLJI - X(I)) 20 CONTINUE IF (SMINUS .GT. SPLUS) B = -B X(J) = ZERO C *** UPDATE PARTIAL SUMS *** CALL SV2AXY(J, X, B, L(J0+1), X) 30 CONTINUE C C *** NORMALIZE X *** C 40 T = SV2NRM(P, X) IF (T .LE. ZERO) GO TO 80 T = ONE / T DO 50 I = 1, P 50 X(I) = T*X(I) C C *** COMPUTE L*X = Y AND RETURN SVMAX = TWONORM(Y) *** C DO 60 JJJ = 1, P J = PPLUS1 - JJJ JI = J*(J-1)/2 + 1 Y(J) = SD7TPR(J, L(JI), X) 60 CONTINUE C C *** NORMALIZE Y AND SET X = (L**T)*Y *** C T = ONE / SV2NRM(P, Y) JI = 1 DO 70 I = 1, P YI = T * Y(I) X(I) = ZERO CALL SV2AXY(I, X, YI, L(JI), X) JI = JI + I 70 CONTINUE SL7SVX = SV2NRM(P, X) GO TO 999 C 80 SL7SVX = ZERO C 999 RETURN C *** LAST CARD OF SL7SVX FOLLOWS *** END c ================================================================== SUBROUTINE SA7SST(IV, LIV, LV, V) C C *** ASSESS CANDIDATE STEP (***SOL VERSION 2.3) *** C c ------------------------------------------------------------------ INTEGER LIV, LV INTEGER IV(LIV) REAL V(LV) C C *** PURPOSE *** C C THIS SUBROUTINE IS CALLED BY AN UNCONSTRAINED MINIMIZATION C ROUTINE TO ASSESS THE NEXT CANDIDATE STEP. IT MAY RECOMMEND ONE C OF SEVERAL COURSES OF ACTION, SUCH AS ACCEPTING THE STEP, RECOM- C PUTING IT USING THE SAME OR A NEW QUADRATIC MODEL, OR HALTING DUE C TO CONVERGENCE OR FALSE CONVERGENCE. SEE THE RETURN CODE LISTING C BELOW. C C ------------------------- PARAMETER USAGE -------------------------- C C IV (I/O) INTEGER PARAMETER AND SCRATCH VECTOR -- SEE DESCRIPTION C BELOW OF IV VALUES REFERENCED. C LIV (IN) LENGTH OF IV ARRAY. C LV (IN) LENGTH OF V ARRAY. C V (I/O) REAL PARAMETER AND SCRATCH VECTOR -- SEE DESCRIPTION C BELOW OF V VALUES REFERENCED. C C *** IV VALUES REFERENCED *** C C IV(IRC) (I/O) ON INPUT FOR THE FIRST STEP TRIED IN A NEW ITERATION, C IV(IRC) SHOULD BE SET TO 3 OR 4 (THE VALUE TO WHICH IT IS C SET WHEN STEP IS DEFINITELY TO BE ACCEPTED). ON INPUT C AFTER STEP HAS BEEN RECOMPUTED, IV(IRC) SHOULD BE C UNCHANGED SINCE THE PREVIOUS RETURN OF SA7SST. C ON OUTPUT, IV(IRC) IS A RETURN CODE HAVING ONE OF THE C FOLLOWING VALUES... C 1 = SWITCH MODELS OR TRY SMALLER STEP. C 2 = SWITCH MODELS OR ACCEPT STEP. C 3 = ACCEPT STEP AND DETERMINE V(RADFAC) BY GRADIENT C TESTS. C 4 = ACCEPT STEP, V(RADFAC) HAS BEEN DETERMINED. C 5 = RECOMPUTE STEP (USING THE SAME MODEL). C 6 = RECOMPUTE STEP WITH RADIUS = V(LMAXS) BUT DO NOT C EVAULATE THE OBJECTIVE FUNCTION. C 7 = X-CONVERGENCE (SEE V(XCTOL)). C 8 = RELATIVE FUNCTION CONVERGENCE (SEE V(RFCTOL)). C 9 = BOTH X- AND RELATIVE FUNCTION CONVERGENCE. C 10 = ABSOLUTE FUNCTION CONVERGENCE (SEE V(AFCTOL)). C 11 = SINGULAR CONVERGENCE (SEE V(LMAXS)). C 12 = FALSE CONVERGENCE (SEE V(XFTOL)). C 13 = IV(IRC) WAS OUT OF RANGE ON INPUT. C RETURN CODE I HAS PRECDENCE OVER I+1 FOR I = 9, 10, 11. C IV(MLSTGD) (I/O) SAVED VALUE OF IV(MODEL). C IV(MODEL) (I/O) ON INPUT, IV(MODEL) SHOULD BE AN INTEGER IDENTIFYING C THE CURRENT QUADRATIC MODEL OF THE OBJECTIVE FUNCTION. C IF A PREVIOUS STEP YIELDED A BETTER FUNCTION REDUCTION, C THEN IV(MODEL) WILL BE SET TO IV(MLSTGD) ON OUTPUT. C IV(NFCALL) (IN) INVOCATION COUNT FOR THE OBJECTIVE FUNCTION. C IV(NFGCAL) (I/O) VALUE OF IV(NFCALL) AT STEP THAT GAVE THE BIGGEST C FUNCTION REDUCTION THIS ITERATION. IV(NFGCAL) REMAINS C UNCHANGED UNTIL A FUNCTION REDUCTION IS OBTAINED. C IV(RADINC) (I/O) THE NUMBER OF RADIUS INCREASES (OR MINUS THE NUMBER C OF DECREASES) SO FAR THIS ITERATION. C IV(RESTOR) (OUT) SET TO 1 IF V(F) HAS BEEN RESTORED AND X SHOULD BE C RESTORED TO ITS INITIAL VALUE, TO 2 IF X SHOULD BE SAVED, C TO 3 IF X SHOULD BE RESTORED FROM THE SAVED VALUE, AND TO C 0 OTHERWISE. C IV(STAGE) (I/O) COUNT OF THE NUMBER OF MODELS TRIED SO FAR IN THE C CURRENT ITERATION. C IV(STGLIM) (IN) MAXIMUM NUMBER OF MODELS TO CONSIDER. C IV(SWITCH) (OUT) SET TO 0 UNLESS A NEW MODEL IS BEING TRIED AND IT C GIVES A SMALLER FUNCTION VALUE THAN THE PREVIOUS MODEL, C IN WHICH CASE SA7SST SETS IV(SWITCH) = 1. C IV(TOOBIG) (IN) IS NONZERO IF STEP WAS TOO BIG (E.G. IF IT CAUSED C OVERFLOW). C IV(XIRC) (I/O) VALUE THAT IV(IRC) WOULD HAVE IN THE ABSENCE OF C CONVERGENCE, FALSE CONVERGENCE, AND OVERSIZED STEPS. C C *** V VALUES REFERENCED *** C C V(AFCTOL) (IN) ABSOLUTE FUNCTION CONVERGENCE TOLERANCE. IF THE C ABSOLUTE VALUE OF THE CURRENT FUNCTION VALUE V(F) IS LESS C THAN V(AFCTOL), THEN SA7SST RETURNS WITH IV(IRC) = 10. C V(DECFAC) (IN) FACTOR BY WHICH TO DECREASE RADIUS WHEN IV(TOOBIG) IS C NONZERO. C V(DSTNRM) (IN) THE 2-NORM OF D*STEP. C V(DSTSAV) (I/O) VALUE OF V(DSTNRM) ON SAVED STEP. C V(DST0) (IN) THE 2-NORM OF D TIMES THE NEWTON STEP (WHEN DEFINED, C I.E., FOR V(NREDUC) .GE. 0). C V(F) (I/O) ON BOTH INPUT AND OUTPUT, V(F) IS THE OBJECTIVE FUNC- C TION VALUE AT X. IF X IS RESTORED TO A PREVIOUS VALUE, C THEN V(F) IS RESTORED TO THE CORRESPONDING VALUE. C V(FDIF) (OUT) THE FUNCTION REDUCTION V(F0) - V(F) (FOR THE OUTPUT C VALUE OF V(F) IF AN EARLIER STEP GAVE A BIGGER FUNCTION C DECREASE, AND FOR THE INPUT VALUE OF V(F) OTHERWISE). C V(FLSTGD) (I/O) SAVED VALUE OF V(F). C V(F0) (IN) OBJECTIVE FUNCTION VALUE AT START OF ITERATION. C V(GTSLST) (I/O) VALUE OF V(GTSTEP) ON SAVED STEP. C V(GTSTEP) (IN) INNER PRODUCT BETWEEN STEP AND GRADIENT. C V(INCFAC) (IN) MINIMUM FACTOR BY WHICH TO INCREASE RADIUS. C V(LMAXS) (IN) MAXIMUM REASONABLE STEP SIZE (AND INITIAL STEP BOUND). C IF THE ACTUAL FUNCTION DECREASE IS NO MORE THAN TWICE C WHAT WAS PREDICTED, IF A RETURN WITH IV(IRC) = 7, 8, 9, C OR 10 DOES NOT OCCUR, IF V(DSTNRM) .GT. V(LMAXS), AND IF C V(PREDUC) .LE. V(SCTOL) * ABS(V(F0)), THEN SA7SST RE- C TURNS WITH IV(IRC) = 11. IF SO DOING APPEARS WORTHWHILE, C THEN SA7SST REPEATS THIS TEST WITH V(PREDUC) COMPUTED FOR C A STEP OF LENGTH V(LMAXS) (BY A RETURN WITH IV(IRC) = 6). C V(NREDUC) (I/O) FUNCTION REDUCTION PREDICTED BY QUADRATIC MODEL FOR C NEWTON STEP. IF SA7SST IS CALLED WITH IV(IRC) = 6, I.E., C IF V(PREDUC) HAS BEEN COMPUTED WITH RADIUS = V(LMAXS) FOR C USE IN THE SINGULAR CONVERVENCE TEST, THEN V(NREDUC) IS C SET TO -V(PREDUC) BEFORE THE LATTER IS RESTORED. C V(PLSTGD) (I/O) VALUE OF V(PREDUC) ON SAVED STEP. C V(PREDUC) (I/O) FUNCTION REDUCTION PREDICTED BY QUADRATIC MODEL FOR C CURRENT STEP. C V(RADFAC) (OUT) FACTOR TO BE USED IN DETERMINING THE NEW RADIUS, C WHICH SHOULD BE V(RADFAC)*DST, WHERE DST IS EITHER THE C OUTPUT VALUE OF V(DSTNRM) OR THE 2-NORM OF C DIAG(NEWD)*STEP FOR THE OUTPUT VALUE OF STEP AND THE C UPDATED VERSION, NEWD, OF THE SCALE VECTOR D. FOR C IV(IRC) = 3, V(RADFAC) = 1.0 IS RETURNED. C V(RDFCMN) (IN) MINIMUM VALUE FOR V(RADFAC) IN TERMS OF THE INPUT C VALUE OF V(DSTNRM) -- SUGGESTED VALUE = 0.1. C V(RDFCMX) (IN) MAXIMUM VALUE FOR V(RADFAC) -- SUGGESTED VALUE = 4.0. C V(RELDX) (IN) SCALED RELATIVE CHANGE IN X CAUSED BY STEP, COMPUTED C (E.G.) BY FUNCTION SRLDST AS C MAX (D(I)*ABS(X(I)-X0(I)), 1 .LE. I .LE. P) / C MAX (D(I)*(ABS(X(I))+ABS(X0(I))), 1 .LE. I .LE. P). C V(RFCTOL) (IN) RELATIVE FUNCTION CONVERGENCE TOLERANCE. IF THE C ACTUAL FUNCTION REDUCTION IS AT MOST TWICE WHAT WAS PRE- C DICTED AND V(NREDUC) .LE. V(RFCTOL)*ABS(V(F0)), THEN C SA7SST RETURNS WITH IV(IRC) = 8 OR 9. C V(STPPAR) (IN) MARQUARDT PARAMETER -- 0 MEANS FULL NEWTON STEP. C V(TUNER1) (IN) TUNING CONSTANT USED TO DECIDE IF THE FUNCTION C REDUCTION WAS MUCH LESS THAN EXPECTED. SUGGESTED C VALUE = 0.1. C V(TUNER2) (IN) TUNING CONSTANT USED TO DECIDE IF THE FUNCTION C REDUCTION WAS LARGE ENOUGH TO ACCEPT STEP. SUGGESTED C VALUE = 10**-4. C V(TUNER3) (IN) TUNING CONSTANT USED TO DECIDE IF THE RADIUS C SHOULD BE INCREASED. SUGGESTED VALUE = 0.75. C V(XCTOL) (IN) X-CONVERGENCE CRITERION. IF STEP IS A NEWTON STEP C (V(STPPAR) = 0) HAVING V(RELDX) .LE. V(XCTOL) AND GIVING C AT MOST TWICE THE PREDICTED FUNCTION DECREASE, THEN C SA7SST RETURNS IV(IRC) = 7 OR 9. C V(XFTOL) (IN) FALSE CONVERGENCE TOLERANCE. IF STEP GAVE NO OR ONLY C A SMALL FUNCTION DECREASE AND V(RELDX) .LE. V(XFTOL), C THEN SA7SST RETURNS WITH IV(IRC) = 12. C C ------------------------------ NOTES ------------------------------- C C *** APPLICATION AND USAGE RESTRICTIONS *** C C THIS ROUTINE IS CALLED AS PART OF THE NL2SOL (NONLINEAR C LEAST-SQUARES) PACKAGE. IT MAY BE USED IN ANY UNCONSTRAINED C MINIMIZATION SOLVER THAT USES DOGLEG, GOLDFELD-QUANDT-TROTTER, C OR LEVENBERG-MARQUARDT STEPS. C C *** ALGORITHM NOTES *** C C SEE (1) FOR FURTHER DISCUSSION OF THE ASSESSING AND MODEL C SWITCHING STRATEGIES. WHILE NL2SOL CONSIDERS ONLY TWO MODELS, C SA7SST IS DESIGNED TO HANDLE ANY NUMBER OF MODELS. C C *** USAGE NOTES *** C C ON THE FIRST CALL OF AN ITERATION, ONLY THE I/O VARIABLES C STEP, X, IV(IRC), IV(MODEL), V(F), V(DSTNRM), V(GTSTEP), AND C V(PREDUC) NEED HAVE BEEN INITIALIZED. BETWEEN CALLS, NO I/O C VALUES EXECPT STEP, X, IV(MODEL), V(F) AND THE STOPPING TOLER- C ANCES SHOULD BE CHANGED. C AFTER A RETURN FOR CONVERGENCE OR FALSE CONVERGENCE, ONE CAN C CHANGE THE STOPPING TOLERANCES AND CALL SA7SST AGAIN, IN WHICH C CASE THE STOPPING TESTS WILL BE REPEATED. C C *** REFERENCES *** C C (1) DENNIS, J.E., JR., GAY, D.M., AND WELSCH, R.E. (1981), C AN ADAPTIVE NONLINEAR LEAST-SQUARES ALGORITHM, C ACM TRANS. MATH. SOFTWARE, VOL. 7, NO. 3. C C (2) POWELL, M.J.D. (1970) A FORTRAN SUBROUTINE FOR SOLVING C SYSTEMS OF NONLINEAR ALGEBRAIC EQUATIONS, IN NUMERICAL C METHODS FOR NONLINEAR ALGEBRAIC EQUATIONS, EDITED BY C P. RABINOWITZ, GORDON AND BREACH, LONDON. C C *** HISTORY *** C C JOHN DENNIS DESIGNED MUCH OF THIS ROUTINE, STARTING WITH C IDEAS IN (2). ROY WELSCH SUGGESTED THE MODEL SWITCHING STRATEGY. C DAVID GAY AND STEPHEN PETERS CAST THIS SUBROUTINE INTO A MORE C PORTABLE FORM (WINTER 1977), AND DAVID GAY CAST IT INTO ITS C PRESENT FORM (FALL 1978). C C *** GENERAL *** C C THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH C SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS C MCS-7600324, DCR75-10143, 76-14311DSS, MCS76-11989, AND C MCS-7906671. C C ----------------------- EXTERNAL QUANTITIES ------------------------ C C *** NO EXTERNAL FUNCTIONS AND SUBROUTINES *** C C ------------------------- LOCAL VARIABLES -------------------------- C LOGICAL GOODX INTEGER I, NFC REAL EMAX, EMAXS, GTS, RFAC1, XMAX REAL HALF, ONE, ONEP2, TWO, ZERO C C *** SUBSCRIPTS FOR IV AND V *** C INTEGER AFCTOL, DECFAC, DSTNRM, DSTSAV, DST0, F, FDIF, FLSTGD, F0, 1 GTSLST, GTSTEP, INCFAC, IRC, LMAXS, MLSTGD, MODEL, NFCALL, 2 NFGCAL, NREDUC, PLSTGD, PREDUC, RADFAC, RADINC, RDFCMN, 3 RDFCMX, RELDX, RESTOR, RFCTOL, SCTOL, STAGE, STGLIM, 4 STPPAR, SWITCH, TOOBIG, TUNER1, TUNER2, TUNER3, XCTOL, 5 XFTOL, XIRC C C *** DATA INITIALIZATIONS *** C PARAMETER (HALF=0.5E+0, ONE=1.E+0, ONEP2=1.2E+0, TWO=2.E+0, 1 ZERO=0.E+0) PARAMETER (IRC=29, MLSTGD=32, MODEL=5, NFCALL=6, NFGCAL=7, 1 RADINC=8, RESTOR=9, STAGE=10, STGLIM=11, SWITCH=12, 2 TOOBIG=2, XIRC=13) PARAMETER (AFCTOL=31, DECFAC=22, DSTNRM=2, DST0=3, DSTSAV=18, 1 F=10, FDIF=11, FLSTGD=12, F0=13, GTSLST=14, GTSTEP=4, 2 INCFAC=23, LMAXS=36, NREDUC=6, PLSTGD=15, PREDUC=7, 3 RADFAC=16, RDFCMN=24, RDFCMX=25, RELDX=17, RFCTOL=32, 4 SCTOL=37, STPPAR=5, TUNER1=26, TUNER2=27, TUNER3=28, 5 XCTOL=33, XFTOL=34) C C ++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ C NFC = IV(NFCALL) IV(SWITCH) = 0 IV(RESTOR) = 0 RFAC1 = ONE GOODX = .TRUE. I = IV(IRC) IF (I .GE. 1 .AND. I .LE. 12) 1 GO TO (20,30,10,10,40,280,220,220,220,220,220,170), I IV(IRC) = 13 GO TO 999 C C *** INITIALIZE FOR NEW ITERATION *** C 10 IV(STAGE) = 1 IV(RADINC) = 0 V(FLSTGD) = V(F0) IF (IV(TOOBIG) .EQ. 0) GO TO 110 IV(STAGE) = -1 IV(XIRC) = I GO TO 60 C C *** STEP WAS RECOMPUTED WITH NEW MODEL OR SMALLER RADIUS *** C *** FIRST DECIDE WHICH *** C 20 IF (IV(MODEL) .NE. IV(MLSTGD)) GO TO 30 C *** OLD MODEL RETAINED, SMALLER RADIUS TRIED *** C *** DO NOT CONSIDER ANY MORE NEW MODELS THIS ITERATION *** IV(STAGE) = IV(STGLIM) IV(RADINC) = -1 GO TO 110 C C *** A NEW MODEL IS BEING TRIED. DECIDE WHETHER TO KEEP IT. *** C 30 IV(STAGE) = IV(STAGE) + 1 C C *** NOW WE ADD THE POSSIBILTIY THAT STEP WAS RECOMPUTED WITH *** C *** THE SAME MODEL, PERHAPS BECAUSE OF AN OVERSIZED STEP. *** C 40 IF (IV(STAGE) .GT. 0) GO TO 50 C C *** STEP WAS RECOMPUTED BECAUSE IT WAS TOO BIG. *** C IF (IV(TOOBIG) .NE. 0) GO TO 60 C C *** RESTORE IV(STAGE) AND PICK UP WHERE WE LEFT OFF. *** C IV(STAGE) = -IV(STAGE) I = IV(XIRC) GO TO (20, 30, 110, 110, 70), I C 50 IF (IV(TOOBIG) .EQ. 0) GO TO 70 C C *** HANDLE OVERSIZE STEP *** C IF (IV(RADINC) .GT. 0) GO TO 80 IV(STAGE) = -IV(STAGE) IV(XIRC) = IV(IRC) C 60 V(RADFAC) = V(DECFAC) IV(RADINC) = IV(RADINC) - 1 IV(IRC) = 5 IV(RESTOR) = 1 GO TO 999 C 70 IF (V(F) .LT. V(FLSTGD)) GO TO 110 C C *** THE NEW STEP IS A LOSER. RESTORE OLD MODEL. *** C IF (IV(MODEL) .EQ. IV(MLSTGD)) GO TO 80 IV(MODEL) = IV(MLSTGD) IV(SWITCH) = 1 C C *** RESTORE STEP, ETC. ONLY IF A PREVIOUS STEP DECREASED V(F). C 80 IF (V(FLSTGD) .GE. V(F0)) GO TO 110 IV(RESTOR) = 1 V(F) = V(FLSTGD) V(PREDUC) = V(PLSTGD) V(GTSTEP) = V(GTSLST) IF (IV(SWITCH) .EQ. 0) RFAC1 = V(DSTNRM) / V(DSTSAV) V(DSTNRM) = V(DSTSAV) NFC = IV(NFGCAL) GOODX = .FALSE. C 110 V(FDIF) = V(F0) - V(F) IF (V(FDIF) .GT. V(TUNER2) * V(PREDUC)) GO TO 140 IF (IV(RADINC) .GT. 0) GO TO 140 C C *** NO (OR ONLY A TRIVIAL) FUNCTION DECREASE C *** -- SO TRY NEW MODEL OR SMALLER RADIUS C IF (V(F) .LT. V(F0)) GO TO 120 IV(MLSTGD) = IV(MODEL) V(FLSTGD) = V(F) V(F) = V(F0) IV(RESTOR) = 1 GO TO 130 120 IV(NFGCAL) = NFC 130 IV(IRC) = 1 IF (IV(STAGE) .LT. IV(STGLIM)) GO TO 160 IV(IRC) = 5 IV(RADINC) = IV(RADINC) - 1 GO TO 160 C C *** NONTRIVIAL FUNCTION DECREASE ACHIEVED *** C 140 IV(NFGCAL) = NFC RFAC1 = ONE V(DSTSAV) = V(DSTNRM) IF (V(FDIF) .GT. V(PREDUC)*V(TUNER1)) GO TO 190 C C *** DECREASE WAS MUCH LESS THAN PREDICTED -- EITHER CHANGE MODELS C *** OR ACCEPT STEP WITH DECREASED RADIUS. C IF (IV(STAGE) .GE. IV(STGLIM)) GO TO 150 C *** CONSIDER SWITCHING MODELS *** IV(IRC) = 2 GO TO 160 C C *** ACCEPT STEP WITH DECREASED RADIUS *** C 150 IV(IRC) = 4 C C *** SET V(RADFAC) TO FLETCHER*S DECREASE FACTOR *** C 160 IV(XIRC) = IV(IRC) EMAX = V(GTSTEP) + V(FDIF) V(RADFAC) = HALF * RFAC1 IF (EMAX .LT. V(GTSTEP)) V(RADFAC) = RFAC1 * max(V(RDFCMN), 1 HALF * V(GTSTEP)/EMAX) C C *** DO FALSE CONVERGENCE TEST *** C 170 IF (V(RELDX) .LE. V(XFTOL)) GO TO 180 IV(IRC) = IV(XIRC) IF (V(F) .LT. V(F0)) GO TO 200 GO TO 230 C 180 IV(IRC) = 12 GO TO 240 C C *** HANDLE GOOD FUNCTION DECREASE *** C 190 IF (V(FDIF) .LT. (-V(TUNER3) * V(GTSTEP))) GO TO 210 C C *** INCREASING RADIUS LOOKS WORTHWHILE. SEE IF WE JUST C *** RECOMPUTED STEP WITH A DECREASED RADIUS OR RESTORED STEP C *** AFTER RECOMPUTING IT WITH A LARGER RADIUS. C IF (IV(RADINC) .LT. 0) GO TO 210 IF (IV(RESTOR) .EQ. 1) GO TO 210 C C *** WE DID NOT. TRY A LONGER STEP UNLESS THIS WAS A NEWTON C *** STEP. C V(RADFAC) = V(RDFCMX) GTS = V(GTSTEP) IF (V(FDIF) .LT. (HALF/V(RADFAC) - ONE) * GTS) 1 V(RADFAC) = max(V(INCFAC), HALF*GTS/(GTS + V(FDIF))) IV(IRC) = 4 IF (V(STPPAR) .EQ. ZERO) GO TO 230 IF (V(DST0) .GE. ZERO .AND. (V(DST0) .LT. TWO*V(DSTNRM) 1 .OR. V(NREDUC) .LT. ONEP2*V(FDIF))) GO TO 230 C *** STEP WAS NOT A NEWTON STEP. RECOMPUTE IT WITH C *** A LARGER RADIUS. IV(IRC) = 5 IV(RADINC) = IV(RADINC) + 1 C C *** SAVE VALUES CORRESPONDING TO GOOD STEP *** C 200 V(FLSTGD) = V(F) IV(MLSTGD) = IV(MODEL) IF (IV(RESTOR) .NE. 1) IV(RESTOR) = 2 V(DSTSAV) = V(DSTNRM) IV(NFGCAL) = NFC V(PLSTGD) = V(PREDUC) V(GTSLST) = V(GTSTEP) GO TO 230 C C *** ACCEPT STEP WITH RADIUS UNCHANGED *** C 210 V(RADFAC) = ONE IV(IRC) = 3 GO TO 230 C C *** COME HERE FOR A RESTART AFTER CONVERGENCE *** C 220 IV(IRC) = IV(XIRC) IF (V(DSTSAV) .GE. ZERO) GO TO 240 IV(IRC) = 12 GO TO 240 C C *** PERFORM CONVERGENCE TESTS *** C 230 IV(XIRC) = IV(IRC) 240 IF (IV(RESTOR) .EQ. 1 .AND. V(FLSTGD) .LT. V(F0)) IV(RESTOR) = 3 IF (abs(V(F)) .LT. V(AFCTOL)) IV(IRC) = 10 IF (HALF * V(FDIF) .GT. V(PREDUC)) GO TO 999 EMAX = V(RFCTOL) * abs(V(F0)) EMAXS = V(SCTOL) * abs(V(F0)) IF (V(PREDUC) .LE. EMAXS .AND. (V(DSTNRM) .GT. V(LMAXS) .OR. 1 V(STPPAR) .EQ. ZERO)) IV(IRC) = 11 IF (V(DST0) .LT. ZERO) GO TO 250 I = 0 IF ((V(NREDUC) .GT. ZERO .AND. V(NREDUC) .LE. EMAX) .OR. 1 (V(NREDUC) .EQ. ZERO .AND. V(PREDUC) .EQ. ZERO)) I = 2 IF (V(STPPAR) .EQ. ZERO .AND. V(RELDX) .LE. V(XCTOL) 1 .AND. GOODX) I = I + 1 IF (I .GT. 0) IV(IRC) = I + 6 C C *** CONSIDER RECOMPUTING STEP OF LENGTH V(LMAXS) FOR SINGULAR C *** CONVERGENCE TEST. C 250 IF (IV(IRC) .GT. 5 .AND. IV(IRC) .NE. 12) GO TO 999 IF (V(STPPAR) .EQ. ZERO) GO TO 999 IF (V(DSTNRM) .GT. V(LMAXS)) GO TO 260 IF (V(PREDUC) .GE. EMAXS) GO TO 999 IF (V(DST0) .LE. ZERO) GO TO 270 IF (HALF * V(DST0) .LE. V(LMAXS)) GO TO 999 GO TO 270 260 IF (HALF * V(DSTNRM) .LE. V(LMAXS)) GO TO 999 XMAX = V(LMAXS) / V(DSTNRM) IF (XMAX * (TWO - XMAX) * V(PREDUC) .GE. EMAXS) GO TO 999 270 IF (V(NREDUC) .LT. ZERO) GO TO 290 C C *** RECOMPUTE V(PREDUC) FOR USE IN SINGULAR CONVERGENCE TEST *** C V(GTSLST) = V(GTSTEP) V(DSTSAV) = V(DSTNRM) IF (IV(IRC) .EQ. 12) V(DSTSAV) = -V(DSTSAV) V(PLSTGD) = V(PREDUC) I = IV(RESTOR) IV(RESTOR) = 2 IF (I .EQ. 3) IV(RESTOR) = 0 IV(IRC) = 6 GO TO 999 C C *** PERFORM SINGULAR CONVERGENCE TEST WITH RECOMPUTED V(PREDUC) *** C 280 V(GTSTEP) = V(GTSLST) V(DSTNRM) = abs(V(DSTSAV)) IV(IRC) = IV(XIRC) IF (V(DSTSAV) .LE. ZERO) IV(IRC) = 12 V(NREDUC) = -V(PREDUC) V(PREDUC) = V(PLSTGD) IV(RESTOR) = 3 290 IF (-V(NREDUC) .LE. V(RFCTOL) * abs(V(F0))) IV(IRC) = 11 C 999 RETURN C C *** LAST CARD OF SA7SST FOLLOWS *** END c ================================================================== SUBROUTINE SITSUM(D, G, IV, LIV, LV, P, V, X) C C *** PRINT ITERATION SUMMARY FOR ***SOL (VERSION 2.3) *** C c 6/28/90 CLL Added test before the GO TO (...), IV1 c to avoid printing the termination message when IV1 = 1 to 4 and c no other printing has been requested. C *** PARAMETER DECLARATIONS *** C INTEGER LIV, LV, P INTEGER IV(LIV) REAL D(P), G(P), V(LV), X(P) c ------------------------------------------------------------------ C *** LOCAL VARIABLES *** C INTEGER ALG, I, IV1, M, NF, NG, OL, PU CHARACTER*4 MODEL1(6), MODEL2(6) REAL NRELDF, OLDF, PRELDF, RELDF, ZERO C C *** NO EXTERNAL FUNCTIONS OR SUBROUTINES *** C C *** SUBSCRIPTS FOR IV AND V *** C INTEGER ALGSAV, DSTNRM, F, FDIF, F0, NEEDHD, NFCALL, NFCOV, NGCOV, 1 NGCALL, NITER, NREDUC, OUTLEV, PREDUC, PRNTIT, PRUNIT, 2 RELDX, SOLPRT, STATPR, STPPAR, SUSED, X0PRT integer COVPRT, PARPRT C C *** IV SUBSCRIPT VALUES *** C PARAMETER (ALGSAV=51, NEEDHD=36, NFCALL=6, NFCOV=52, NGCALL=30, 1 NGCOV=53, NITER=31, OUTLEV=19, PRNTIT=39, PRUNIT=21, 2 SOLPRT=22, STATPR=23, SUSED=64, X0PRT=24) parameter (COVPRT=14, PARPRT=20) C C *** V SUBSCRIPT VALUES *** C PARAMETER (DSTNRM=2, F=10, F0=13, FDIF=11, NREDUC=6, PREDUC=7, 1 RELDX=17, STPPAR=5) C PARAMETER (ZERO=0.E+0) DATA MODEL1/' ',' ',' ',' ',' G ',' S '/, 1 MODEL2/' G ',' S ','G-S ','S-G ','-S-G','-G-S'/ C C ------------------------------ BODY -------------------------------- C 30 FORMAT(/' IT NF',6X,'F',7X,'RELDF',3X,'PRELDF',3X,'RELDX', 1 2X,'MODEL STPPAR') 40 FORMAT(/' IT NF',7X,'F',8X,'RELDF',4X,'PRELDF',4X,'RELDX', 1 3X,'STPPAR') 70 FORMAT(/' IT NF',6X,'F',7X,'RELDF',3X,'PRELDF',3X,'RELDX', 1 2X,'MODEL STPPAR',2X,'D*STEP',2X,'NPRELDF') 80 FORMAT(/' IT NF',7X,'F',8X,'RELDF',4X,'PRELDF',4X,'RELDX', 1 3X,'STPPAR',3X,'D*STEP',3X,'NPRELDF') 100 FORMAT(I6,I5,g10.3,2g9.2,g8.1,A3,A4,2g8.1,g9.2) 110 FORMAT(I6,I5,g11.3,2g10.2,3g9.1,g10.2) 140 FORMAT(/' ***** COEFFICIENT CONVERGENCE *****') 160 FORMAT(/' ***** RELATIVE FUNCTION CONVERGENCE *****') 200 FORMAT(/' ***** ABSOLUTE FUNCTION CONVERGENCE *****') 220 FORMAT(/' ***** SINGULAR CONVERGENCE *****') 240 FORMAT(/' ***** FALSE CONVERGENCE *****') 260 FORMAT(/' ***** FUNCTION EVALUATION LIMIT *****') 280 FORMAT(/' ***** ITERATION LIMIT *****') 300 FORMAT(/' ***** STOPX *****') 320 FORMAT(/' ***** INITIAL F(X) CANNOT BE COMPUTED *****') 340 FORMAT(/' ***** BAD PARAMETERS (Internal error) *****') C340 FORMAT(/' ***** BAD PARAMETERS TO ASSESS *****') 360 FORMAT(/' ***** GRADIENT COULD NOT BE COMPUTED *****') 380 FORMAT(/' ***** IV(1) =',I5,' *****') 410 FORMAT(/' 0',I5,g10.3) 420 FORMAT(/' 0',I5,g11.3) 450 FORMAT(/' FUNCTION',g17.6,' RELDX',g17.3/' FUNC. EVALS', 1 I8,9X,'GRAD. EVALS',I8/' PRELDF',g16.3,6X,'NPRELDF',g15.3) 470 FORMAT(/' I FINAL X(I)',8X,'D(I)',10X,'G(I)'/) 490 FORMAT(1X,I5,g16.6,2g14.3) 510 FORMAT(/' INCONSISTENT DIMENSIONS') c PU = IV(PRUNIT) IF (PU .EQ. 0) GO TO 999 IV1 = IV(1) IF (IV1 .GT. 62) IV1 = IV1 - 51 OL = IV(OUTLEV) ALG = MOD(IV(ALGSAV)-1,2) + 1 IF (IV1 .LT. 2 .OR. IV1 .GT. 15) GO TO 370 IF (IV1 .GE. 12) GO TO 120 IF (IV1 .EQ. 2 .AND. IV(NITER) .EQ. 0) GO TO 390 IF (OL .EQ. 0) GO TO 120 IF (IV1 .GE. 10 .AND. IV(PRNTIT) .EQ. 0) GO TO 120 IF (IV1 .GT. 2) GO TO 10 IV(PRNTIT) = IV(PRNTIT) + 1 IF (IV(PRNTIT) .LT. abs(OL)) GO TO 999 10 NF = IV(NFCALL) - abs(IV(NFCOV)) IV(PRNTIT) = 0 RELDF = ZERO PRELDF = ZERO OLDF = max(abs(V(F0)), abs(V(F))) IF (OLDF .LE. ZERO) GO TO 20 RELDF = V(FDIF) / OLDF PRELDF = V(PREDUC) / OLDF 20 IF (OL .GT. 0) GO TO 60 C C *** PRINT SHORT SUMMARY LINE *** C IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 1) WRITE(PU,30) IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 2) WRITE(PU,40) IV(NEEDHD) = 0 IF (ALG .EQ. 2) GO TO 50 M = IV(SUSED) WRITE(PU,100) IV(NITER), NF, V(F), RELDF, PRELDF, V(RELDX), 1 MODEL1(M), MODEL2(M), V(STPPAR) GO TO 120 C 50 WRITE(PU,110) IV(NITER), NF, V(F), RELDF, PRELDF, V(RELDX), 1 V(STPPAR) GO TO 120 C C *** PRINT LONG SUMMARY LINE *** C 60 IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 1) WRITE(PU,70) IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 2) WRITE(PU,80) IV(NEEDHD) = 0 NRELDF = ZERO IF (OLDF .GT. ZERO) NRELDF = V(NREDUC) / OLDF IF (ALG .EQ. 2) GO TO 90 M = IV(SUSED) WRITE(PU,100) IV(NITER), NF, V(F), RELDF, PRELDF, V(RELDX), 1 MODEL1(M), MODEL2(M), V(STPPAR), V(DSTNRM), NRELDF GO TO 120 C 90 WRITE(PU,110) IV(NITER), NF, V(F), RELDF, PRELDF, 1 V(RELDX), V(STPPAR), V(DSTNRM), NRELDF C 120 IF (IV1 .LE. 2) GO TO 999 I = IV(STATPR) IF (I .EQ. (-1)) GO TO 460 IF (I + IV1 .LT. 0) GO TO 460 c c The following test skips printing the termination message c when convergence is satisfactory and no other printing has been c requested. c if(IV1 .ge. 3 .and. IV1 .le. 6 .and. * IV(COVPRT) .eq. 0 .and. IV(OUTLEV) .eq. 0 .and. * IV(PARPRT) .eq. 0 .and. IV(SOLPRT) .eq. 0 .and. * IV(STATPR) .eq. 0 .and. IV(X0PRT) .eq. 0) go to 430 GO TO (999, 999, 130, 150, 170, 190, 210, 230, 250, 270, 290, 310, 1 330, 350, 500), IV1 C 130 WRITE(PU,140) GO TO 430 C 150 WRITE(PU,160) GO TO 430 C 170 WRITE(PU,'(1x/a)') * ' ***** COEFFICIENT AND RELATIVE FUNCTION CONVERGENCE *****' GO TO 430 C 190 WRITE(PU,200) GO TO 430 C 210 WRITE(PU,220) GO TO 430 C 230 WRITE(PU,240) GO TO 430 C 250 WRITE(PU,260) GO TO 430 C 270 WRITE(PU,280) GO TO 430 C 290 WRITE(PU,300) GO TO 430 C 310 WRITE(PU,320) C GO TO 390 C 330 WRITE(PU,340) GO TO 999 C 350 WRITE(PU,360) IF (IV(NITER) .GT. 0) GO TO 460 GO TO 390 C 370 WRITE(PU,380) IV(1) GO TO 999 C C *** INITIAL CALL ON SITSUM *** C 390 IF (IV(X0PRT) .NE. 0) then WRITE(PU,'(/'' I INITIAL X(I)'',8X,''D(I)''/)') DO 400 I = 1, P 400 WRITE(PU,'(1X,I5,g17.6,g14.3)') I, X(I), D(I) END IF C *** THE FOLLOWING ARE TO AVOID UNDEFINED VARIABLES WHEN THE C *** FUNCTION EVALUATION LIMIT IS 1... V(DSTNRM) = ZERO V(FDIF) = ZERO V(NREDUC) = ZERO V(PREDUC) = ZERO V(RELDX) = ZERO IF (IV1 .GE. 12) GO TO 999 IV(NEEDHD) = 0 IV(PRNTIT) = 0 IF (OL .EQ. 0) GO TO 999 IF (OL .LT. 0 .AND. ALG .EQ. 1) WRITE(PU,30) IF (OL .LT. 0 .AND. ALG .EQ. 2) WRITE(PU,40) IF (OL .GT. 0 .AND. ALG .EQ. 1) WRITE(PU,70) IF (OL .GT. 0 .AND. ALG .EQ. 2) WRITE(PU,80) IF (ALG .EQ. 1) WRITE(PU,410) IV(NFCALL), V(F) IF (ALG .EQ. 2) WRITE(PU,420) IV(NFCALL), V(F) GO TO 999 C C *** PRINT VARIOUS INFORMATION REQUESTED ON SOLUTION *** C 430 IV(NEEDHD) = 1 IF (IV(STATPR) .LE. 0) GO TO 460 OLDF = max(abs(V(F0)), abs(V(F))) PRELDF = ZERO NRELDF = ZERO IF (OLDF .LE. ZERO) GO TO 440 PRELDF = V(PREDUC) / OLDF NRELDF = V(NREDUC) / OLDF 440 NF = IV(NFCALL) - IV(NFCOV) NG = IV(NGCALL) - IV(NGCOV) WRITE(PU,450) V(F), V(RELDX), NF, NG, PRELDF, NRELDF C 460 IF (IV(SOLPRT) .EQ. 0) GO TO 999 IV(NEEDHD) = 1 IF (IV(ALGSAV) .GT. 2) GO TO 999 WRITE(PU,470) DO 480 I = 1, P 480 WRITE(PU,490) I, X(I), D(I), G(I) GO TO 999 C 500 WRITE(PU,510) 999 RETURN C *** LAST CARD OF SITSUM FOLLOWS *** END c ================================================================== SUBROUTINE SL7ITV(N, X, L, Y) C C *** SOLVE (L**T)*X = Y, WHERE L IS AN N X N LOWER TRIANGULAR C *** MATRIX STORED COMPACTLY BY ROWS. X AND Y MAY OCCUPY THE SAME C *** STORAGE. *** c ------------------------------------------------------------------ INTEGER N REAL X(N), L(*), Y(N) INTEGER I, II, IJ, IM1, I0, J, NP1 REAL XI, ZERO PARAMETER (ZERO=0.E+0) C DO 10 I = 1, N 10 X(I) = Y(I) NP1 = N + 1 I0 = N*(N+1)/2 DO 30 II = 1, N I = NP1 - II XI = X(I)/L(I0) X(I) = XI IF (I .LE. 1) GO TO 999 I0 = I0 - I IF (XI .EQ. ZERO) GO TO 30 IM1 = I - 1 DO 20 J = 1, IM1 IJ = I0 + J X(J) = X(J) - XI*L(IJ) 20 CONTINUE 30 CONTINUE 999 RETURN C *** LAST CARD OF SL7ITV FOLLOWS *** END c ================================================================== SUBROUTINE SL7SQR(N, A, L) C C *** COMPUTE A = LOWER TRIANGLE OF L*(L**T), WITH BOTH C *** L AND A STORED COMPACTLY BY ROWS. (BOTH MAY OCCUPY THE C *** SAME STORAGE. C C *** PARAMETERS *** C c ------------------------------------------------------------------ INTEGER N REAL A(*), L(*) C DIMENSION A(N*(N+1)/2), L(N*(N+1)/2) C C *** LOCAL VARIABLES *** C INTEGER I, II, IJ, IK, IP1, I0, J, JJ, JK, J0, K, NP1 REAL T C NP1 = N + 1 I0 = N*(N+1)/2 DO 30 II = 1, N I = NP1 - II IP1 = I + 1 I0 = I0 - I J0 = I*(I+1)/2 DO 20 JJ = 1, I J = IP1 - JJ J0 = J0 - J T = 0.0E0 DO 10 K = 1, J IK = I0 + K JK = J0 + K T = T + L(IK)*L(JK) 10 CONTINUE IJ = I0 + J A(IJ) = T 20 CONTINUE 30 CONTINUE END c ================================================================== SUBROUTINE SL7SRT(N1, N, L, A, IRC) C C *** COMPUTE ROWS N1 THROUGH N OF THE CHOLESKY FACTOR L OF C *** A = L*(L**T), WHERE L AND THE LOWER TRIANGLE OF A ARE BOTH C *** STORED COMPACTLY BY ROWS (AND MAY OCCUPY THE SAME STORAGE). C *** IRC = 0 MEANS ALL WENT WELL. IRC = J MEANS THE LEADING C *** PRINCIPAL J X J SUBMATRIX OF A IS NOT POSITIVE DEFINITE -- C *** AND L(J*(J+1)/2) CONTAINS THE (NONPOS.) REDUCED J-TH DIAGONAL. C C *** PARAMETERS *** C c ------------------------------------------------------------------ INTEGER N1, N, IRC REAL L(*), A(*) C DIMENSION L(N*(N+1)/2), A(N*(N+1)/2) C C *** LOCAL VARIABLES *** C INTEGER I, IJ, IK, IM1, I0, J, JK, JM1, J0, K REAL T, TD, ZERO C PARAMETER (ZERO=0.E+0) C C *** BODY *** C I0 = N1 * (N1 - 1) / 2 DO 50 I = N1, N TD = ZERO IF (I .EQ. 1) GO TO 40 J0 = 0 IM1 = I - 1 DO 30 J = 1, IM1 T = ZERO IF (J .EQ. 1) GO TO 20 JM1 = J - 1 DO 10 K = 1, JM1 IK = I0 + K JK = J0 + K T = T + L(IK)*L(JK) 10 CONTINUE 20 IJ = I0 + J J0 = J0 + J T = (A(IJ) - T) / L(J0) L(IJ) = T TD = TD + T*T 30 CONTINUE 40 I0 = I0 + I T = A(I0) - TD IF (T .LE. ZERO) GO TO 60 L(I0) = sqrt(T) * DNSGB (9 of 10) 50 CONTINUE C IRC = 0 GO TO 999 C 60 L(I0) = T IRC = I C 999 RETURN C C *** LAST CARD OF SL7SRT *** END c ================================================================== SUBROUTINE SL7TVM(N, X, L, Y) C C *** COMPUTE X = (L**T)*Y, WHERE L IS AN N X N LOWER C *** TRIANGULAR MATRIX STORED COMPACTLY BY ROWS. X AND Y MAY C *** OCCUPY THE SAME STORAGE. *** c ------------------------------------------------------------------ INTEGER N REAL X(N), L(*), Y(N) C DIMENSION L(N*(N+1)/2) INTEGER I, IJ, I0, J REAL YI, ZERO PARAMETER (ZERO=0.E+0) C I0 = 0 DO 20 I = 1, N YI = Y(I) X(I) = ZERO DO 10 J = 1, I IJ = I0 + J X(J) = X(J) + YI*L(IJ) 10 CONTINUE I0 = I0 + I 20 CONTINUE C *** LAST CARD OF SL7TVM FOLLOWS *** END c ================================================================== SUBROUTINE SL7VML(N, X, L, Y) C C *** COMPUTE X = L*Y, WHERE L IS AN N X N LOWER TRIANGULAR C *** MATRIX STORED COMPACTLY BY ROWS. X AND Y MAY OCCUPY THE SAME C *** STORAGE. *** C c ------------------------------------------------------------------ INTEGER N REAL X(N), L(*), Y(N) C DIMENSION L(N*(N+1)/2) INTEGER I, II, IJ, I0, J, NP1 REAL T, ZERO PARAMETER (ZERO=0.E+0) C NP1 = N + 1 I0 = N*(N+1)/2 DO 20 II = 1, N I = NP1 - II I0 = I0 - I T = ZERO DO 10 J = 1, I IJ = I0 + J T = T + L(IJ)*Y(J) 10 CONTINUE X(I) = T 20 CONTINUE C *** LAST CARD OF SL7VML FOLLOWS *** END c ================================================================== REAL FUNCTION SRLDST(P, D, X, X0) C C *** COMPUTE AND RETURN RELATIVE DIFFERENCE BETWEEN X AND X0 *** C *** NL2SOL VERSION 2.2 *** c ------------------------------------------------------------------ INTEGER P REAL D(P), X(P), X0(P) C INTEGER I REAL EMAX, T, XMAX, ZERO PARAMETER (ZERO=0.E+0) C C *** BODY *** C EMAX = ZERO XMAX = ZERO DO 10 I = 1, P T = abs(D(I) * (X(I) - X0(I))) IF (EMAX .LT. T) EMAX = T T = D(I) * (abs(X(I)) + abs(X0(I))) IF (XMAX .LT. T) XMAX = T 10 CONTINUE SRLDST = ZERO IF (XMAX .GT. ZERO) SRLDST = EMAX / XMAX C *** LAST CARD OF SRLDST FOLLOWS *** END c ================================================================== SUBROUTINE SV2AXY(P, W, A, X, Y) C C *** SET W = A*X + Y -- W, X, Y = P-VECTORS, A = SCALAR *** C c ------------------------------------------------------------------ INTEGER P REAL A, W(P), X(P), Y(P) C INTEGER I C DO 10 I = 1, P 10 W(I) = A*X(I) + Y(I) RETURN END c ================================================================== SUBROUTINE SV7CPY(P, Y, X) C C *** SET Y = X, WHERE X AND Y ARE P-VECTORS *** C c ------------------------------------------------------------------ INTEGER P REAL X(P), Y(P) C INTEGER I C DO 10 I = 1, P 10 Y(I) = X(I) RETURN END c ================================================================== SUBROUTINE SV7SCL(N, X, A, Y) C C *** SET X(I) = A*Y(I), I = 1(1)N *** C c ------------------------------------------------------------------ INTEGER N REAL A, X(N), Y(N) C INTEGER I C DO 10 I = 1, N 10 X(I) = A * Y(I) C *** LAST LINE OF SV7SCL FOLLOWS *** END c ================================================================== SUBROUTINE SV7SCP(P, Y, SS) C C *** SET P-VECTOR Y TO SCALAR SS *** C c ------------------------------------------------------------------ INTEGER P REAL SS, Y(P) C INTEGER I C DO 10 I = 1, P 10 Y(I) = SS RETURN END