SUBROUTINE DG7QTS(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 DOUBLE PRECISION D(P), DIG(P), DIHDI(1), L(1), V(21), STEP(P), 1 W(1) 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 (DG7QTS 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 DD7TPR - RETURNS INNER PRODUCT OF TWO VECTORS. C DL7ITV - APPLIES INVERSE-TRANSPOSE OF COMPACT LOWER TRIANG. MATRIX. C DL7IVM - APPLIES INVERSE OF COMPACT LOWER TRIANG. MATRIX. C DL7SRT - FINDS CHOLESKY FACTOR (OF COMPACTLY STORED LOWER TRIANG.). C DL7SVN - RETURNS APPROX. TO MIN. SING. VALUE OF LOWER TRIANG. MATRIX. C DR7MDC - RETURNS MACHINE-DEPENDENT CONSTANTS. C DV2NRM - 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 DOUBLE PRECISION 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 *** DOUBLE PRECISION BIG, DGXFAC, EPSFAC, FOUR, HALF, KAPPA, NEGONE, 1 ONE, P001, SIX, THREE, TWO, ZERO C C *** INTRINSIC FUNCTIONS *** C/+ DOUBLE PRECISION DSQRT C/ C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C DOUBLE PRECISION DD7TPR, DL7SVN, DR7MDC, DV2NRM EXTERNAL DD7TPR, DL7ITV, DL7IVM,DL7SRT, DL7SVN, DR7MDC, DV2NRM C C *** SUBSCRIPTS FOR V *** C INTEGER DGNORM, DSTNRM, DST0, EPSLON, GTSTEP, STPPAR, NREDUC, 1 PHMNFC, PHMXFC, PREDUC, RADIUS, RAD0 C/6 C DATA DGNORM/1/, DSTNRM/2/, DST0/3/, EPSLON/19/, GTSTEP/4/, C 1 NREDUC/6/, PHMNFC/20/, PHMXFC/21/, PREDUC/7/, RADIUS/8/, C 2 RAD0/9/, STPPAR/5/ C/7 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) C/ C C/6 C DATA EPSFAC/50.0D+0/, FOUR/4.0D+0/, HALF/0.5D+0/, C 1 KAPPA/2.0D+0/, NEGONE/-1.0D+0/, ONE/1.0D+0/, P001/1.0D-3/, C 2 SIX/6.0D+0/, THREE/3.0D+0/, TWO/2.0D+0/, ZERO/0.0D+0/ C/7 PARAMETER (EPSFAC=50.0D+0, FOUR=4.0D+0, HALF=0.5D+0, 1 KAPPA=2.0D+0, NEGONE=-1.0D+0, ONE=1.0D+0, P001=1.0D-3, 2 SIX=6.0D+0, THREE=3.0D+0, TWO=2.0D+0, ZERO=0.0D+0) SAVE DGXFAC C/ DATA BIG/0.D+0/, DGXFAC/0.D+0/ C C *** BODY *** C IF (BIG .LE. ZERO) BIG = DR7MDC(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*DMIN1(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) = DV2NRM(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 = DABS(DIHDI(I)) IF (T1 .LT. T) T1 = T 20 CONTINUE W(DGGDMX) = T1 C C *** TRY ALPHA = 0 *** C 30 CALL DL7SRT(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 DL7ITV(IRC, W, L, W) T1 = DV2NRM(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 = DL7SVN(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 DL7IVM(P, W(Q), L, DIG) GTSTA = DD7TPR(P, W(Q), W(Q)) V(NREDUC) = HALF * GTSTA CALL DL7ITV(P, W(Q), L, W(Q)) DST = DV2NRM(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 = DABS(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 = DABS(DIHDI(K1)) SI = W(I) J = DIAG0 + I T1 = HALF * (AKK - W(J) + SI - AKI) T1 = T1 + DSQRT(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 = DABS(DIHDI(K1)) SI = W(I) J = DIAG0 + I T1 = HALF * (W(J) + SI - AKI - AKK) T1 = T1 + DSQRT(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 = DMAX1(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 = DABS(V(STPPAR)) * V(RAD0)/RAD ALPHAK = DMIN1(UK, DMAX1(ALPHAK, LK)) C IF (IRC .NE. 0) GO TO 210 C C *** COMPUTE L0 FOR POSITIVE DEFINITE H *** C CALL DL7IVM(P, W, L, W(Q)) T = DV2NRM(P, W) W(PHIPIN) = RAD / T / T LK = DMAX1(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 * DMAX1(P001, DSQRT(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 DL7SRT(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 DL7ITV(IRC, W, L, W) T1 = DV2NRM(IRC, W) LK = ALPHAK - T/T1/T1 V(DST0) = -LK IF (UK .LT. LK) UK = 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 DL7IVM(P, W(Q), L, DIG) GTSTA = DD7TPR(P, W(Q), W(Q)) CALL DL7ITV(P, W(Q), L, W(Q)) DST = DV2NRM(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 DMIN1 IS NECESSARY BECAUSE OF RESTARTS *** IF (PHI .LT. ZERO) UK = DMIN1(UK, ALPHAK) C *** KAMIN = 0 ONLY IFF THE GRADIENT VANISHES *** IF (KAMIN .EQ. 0) GO TO 210 CALL DL7IVM(P, W, L, W(Q)) C *** THE FOLLOWING, COMMENTED CALCULATION OF ALPHAK IS SOMETIMES C *** SAFER BUT WORSE IN PERFORMANCE... C T1 = DST / DV2NRM(P, W) C ALPHAK = ALPHAK + T1 * (PHI/RAD) * T1 T1 = DV2NRM(P, W) ALPHAK = ALPHAK + (PHI/T1) * (DST/T1) * (DST/RAD) LK = DMAX1(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 * (DABS(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 = DABS(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 = DMAX1(LK, T - W(EMAX)) IF (V(DST0) .GT. ZERO) LK = DMAX1(LK, (V(DST0)-RAD)*W(PHIPIN)) GO TO 250 C C *** BIGGER RADIUS *** 320 IF (ALPHAK .GT. ZERO) UK = DMIN1(UK, W(UK0)) LK = DMAX1(ZERO, -V(DST0), T - W(EMAX)) IF (V(DST0) .GT. ZERO) LK = DMAX1(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 + DMIN1(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 DL7SVN TO OBTAIN APPROXIMATE EIGENVECTOR CORRESPONDING C *** TO SMALLEST EIGENVALUE OF (D**-1)*H*(D**-1). DL7SVN RETURNS C *** X AND W WITH L*W = X. C 340 T = DL7SVN(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 DL7ITV(P, W, L, W) T2 = ONE/DV2NRM(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 = DD7TPR(P, W(Q), W) T1 = (RAD + DST) * (RAD - DST) ROOT = DSQRT(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) = DMIN1(V(DST0), T2**2 - ALPHAK) LK = DMAX1(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 * DR7MDC(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*DD7TPR(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) = DD7TPR(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 999 RETURN C C *** LAST CARD OF DG7QTS FOLLOWS *** END