*ASSESS SUBROUTINE ASSESS (D, IV, P, STEP, STLSTG, V, X, X0) C C LATEST REVISION - 03/15/90 (JRD) C C C *** ASSESS CANDIDATE STEP (NL2SOL VERSION 2.2) *** 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 C VARIABLE DECLARATIONS C C SCALAR ARGUMENTS INTEGER + P C C ARRAY ARGUMENTS REAL + D(P),STEP(P),STLSTG(P),V(*),X(P),X0(P) INTEGER + IV(*) C C LOCAL SCALARS REAL + EMAX,GTS,HALF,ONE,RELDX1,RFAC1,TEMP,TWO,XMAX,ZERO INTEGER + AFCTOL,DECFAC,DST0,DSTNRM,DSTSAV,F,F0,FDIF,FLSTGD,GTSLST, + GTSTEP,I,INCFAC,IRC,LMAX0,MLSTGD,MODEL,NFC,NFCALL,NFGCAL, + NREDUC,PLSTGD,PREDUC,RADFAC,RADINC,RDFCMN,RDFCMX,RELDX, + RESTOR,RFCTOL,STAGE,STGLIM,STPPAR,SWITCH,TOOBIG,TUNER1, + TUNER2,TUNER3,XCTOL,XFTOL,XIRC LOGICAL + GOODX C C EXTERNAL FUNCTIONS REAL + R1MACH,RELDST EXTERNAL R1MACH,RELDST C C EXTERNAL SUBROUTINES EXTERNAL VCOPY C C INTRINSIC FUNCTIONS INTRINSIC ABS,MAX C C-------------------------- PARAMETER USAGE -------------------------- C C IV (I/O) INTEGER PARAMETER AND SCRATCH VECTOR -- SEE DESCRIPTION C BELOW OF IV VALUES REFERENCED. C D (IN) SCALE VECTOR USED IN COMPUTING V(RELDX) -- SEE BELOW. C P (IN) NUMBER OF PARAMETERS BEING OPTIMIZED. C STEP (I/O) ON INPUT, STEP IS THE STEP TO BE ASSESSED. IT IS UN- C CHANGED ON OUTPUT UNLESS A PREVIOUS STEP ACHIEVED A C BETTER OBJECTIVE FUNCTION REDUCTION, IN WHICH CASE STLSTG C WILL HAVE BEEN COPIED TO STEP. C STLSTG (I/O) WHEN ASSESS RECOMMENDS RECOMPUTING STEP EVEN THOUGH THE C CURRENT (OR A PREVIOUS) STEP YIELDS AN OBJECTIVE FUNC- C TION DECREASE, IT SAVES IN STLSTG THE STEP THAT GAVE THE C BEST FUNCTION REDUCTION SEEN SO FAR (IN THE CURRENT ITERA- C TION). IF THE RECOMPUTED STEP YIELDS A LARGER FUNCTION C VALUE, THEN STEP IS RESTORED FROM STLSTG AND C X = X0 + STEP IS RECOMPUTED. C V (I/O) REAL PARAMETER AND SCRATCH VECTOR -- SEE DESCRIPTION C BELOW OF V VALUES REFERENCED. C X (I/O) ON INPUT, X = X0 + STEP IS THE POINT AT WHICH THE OBJEC- C TIVE FUNCTION HAS JUST BEEN EVALUATED. IF AN EARLIER C STEP YIELDED A BIGGER FUNCTION DECREASE, THEN X IS C RESTORED TO THE CORRESPONDING EARLIER VALUE. OTHERWISE, C IF THE CURRENT STEP DOES NOT GIVE ANY FUNCTION DECREASE, C THEN X IS RESTORED TO X0. C X0 (IN) INITIAL OBJECTIVE FUNCTION PARAMETER VECTOR (AT THE C START OF THE CURRENT ITERATION). 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 ASSESS. 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(LMAX0) 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(LMAX0)). 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 0 UNLESS X AND V(F) HAVE BEEN RESTORED, IN C WHICH CASE ASSESS SETS IV(RESTOR) = 1. 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 ASSESS 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 ASSESS 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(LMAX0) (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(LMAX0), AND IF C V(PREDUC) .LE. V(RFCTOL) * ABS(V(F0)), THEN ASSESS RE- C TURNS WITH IV(IRC) = 11. IF SO DOING APPEARS WORTHWHILE, C THEN ASSESS REPEATS THIS TEST WITH V(PREDUC) COMPUTED FOR C A STEP OF LENGTH V(LMAX0) (BY A RETURN WITH IV(IRC) = 6). C V(NREDUC) (I/O) FUNCTION REDUCTION PREDICTED BY QUADRATIC MODEL FOR C NEWTON STEP. IF ASSESS IS CALLED WITH IV(IRC) = 6, I.E., C IF V(PREDUC) HAS BEEN COMPUTED WITH RADIUS = V(LMAX0) 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) (OUT) SCALED RELATIVE CHANGE IN X CAUSED BY STEP, COMPUTED C BY FUNCTION RELDST 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 IF AN ACCEPTABLE STEP IS RETURNED, THEN V(RELDX) IS COM- C PUTED USING THE OUTPUT (POSSIBLY RESTORED) VALUES OF X C AND STEP. OTHERWISE IT IS COMPUTED USING THE INPUT C VALUES. 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 ASSESS RETURNS WITH IV(IRC) = 8 OR 9. SEE ALSO V(LMAX0). 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 ASSESS 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 ASSESS 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 ASSESS 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 ASSESS 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. (1980), C AN ADAPTIVE NONLINEAR LEAST-SQUARES ALGORITHM, C SUBMITTED TO ACM TRANS. MATH. SOFTWARE. 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 *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C C EXTERNAL RELDST, VCOPY C REAL R1MACH, RELDST C C VCOPY.... COPIES ONE VECTOR TO ANOTHER. C C/ C *** NO COMMON BLOCKS *** C C-------------------------- LOCAL VARIABLES -------------------------- C C LOGICAL GOODX C INTEGER I, NFC C REAL EMAX, GTS, HALF, ONE, RELDX1, RFAC1, C + TEMP, TWO, XMAX, ZERO C C *** SUBSCRIPTS FOR IV AND V *** C C INTEGER AFCTOL, DECFAC, DSTNRM, DSTSAV, DST0, F, FDIF, FLSTGD, F0, C 1 GTSLST, GTSTEP, INCFAC, IRC, LMAX0, MLSTGD, MODEL, NFCALL, C 2 NFGCAL, NREDUC, PLSTGD, PREDUC, RADFAC, RADINC, RDFCMN, C 3 RDFCMX, RELDX, RESTOR, RFCTOL, STAGE, STGLIM, STPPAR, C 4 SWITCH, TOOBIG, TUNER1, TUNER2, TUNER3, XCTOL, XFTOL, C 5 XIRC C C *** DATA INITIALIZATIONS *** C DATA HALF/0.5E0/, ONE/1.0E0/, TWO/2.0E0/, ZERO/0.0E0/ C DATA IRC/3/, MLSTGD/4/, MODEL/5/, NFCALL/6/, + NFGCAL/7/, RADINC/8/, RESTOR/9/, STAGE/10/, + STGLIM/11/, SWITCH/12/, TOOBIG/2/, XIRC/13/ DATA AFCTOL/31/, DECFAC/22/, DSTNRM/2/, DST0/3/, + DSTSAV/18/, F/10/, FDIF/11/, FLSTGD/12/, F0/13/, + GTSLST/14/, GTSTEP/4/, INCFAC/23/, + LMAX0/35/, NREDUC/6/, PLSTGD/15/, PREDUC/7/, + RADFAC/16/, RDFCMN/24/, RDFCMX/25/, + RELDX/17/, RFCTOL/32/, STPPAR/5/, TUNER1/26/, + TUNER2/27/, TUNER3/28/, 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) + GO TO (20,30,10,10,40,360,290,290,290,290,290,140), 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 90 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 90 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, 90, 90, 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 GO TO 999 C 70 IF (V(F) .LT. V(FLSTGD)) GO TO 90 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 90 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 C C *** COMPUTE RELATIVE CHANGE IN X BY CURRENT STEP *** C 90 RELDX1 = RELDST(P, D, X, X0) C C *** RESTORE X AND STEP IF NECESSARY *** C IF (GOODX) GO TO 105 DO 100 I = 1, P STEP(I) = STLSTG(I) X(I) = X0(I) + STLSTG(I) 100 CONTINUE C 105 V(FDIF) = V(F0) - V(F) TEMP = 0.0 IF (V(PREDUC).GT.R1MACH(1)/V(TUNER2)) TEMP = V(TUNER2) * V(PREDUC) IF (V(FDIF).GT.TEMP) GO TO 120 C C *** NO (OR ONLY A TRIVIAL) FUNCTION DECREASE C *** -- SO TRY NEW MODEL OR SMALLER RADIUS C V(RELDX) = RELDX1 IF (V(F) .LT. V(F0)) GO TO 110 IV(MLSTGD) = IV(MODEL) V(FLSTGD) = V(F) V(F) = V(F0) CALL VCOPY(P, X, X0) IV(RESTOR) = 1 GO TO 115 110 IV(NFGCAL) = NFC 115 IV(IRC) = 1 IF (IV(STAGE) .LT. IV(STGLIM)) GO TO 130 IV(IRC) = 5 IV(RADINC) = IV(RADINC) - 1 GO TO 130 C C *** NONTRIVIAL FUNCTION DECREASE ACHIEVED *** C 120 IV(NFGCAL) = NFC RFAC1 = ONE IF (GOODX) V(RELDX) = RELDX1 V(DSTSAV) = V(DSTNRM) IF (V(FDIF) .GT. V(PREDUC)*V(TUNER1)) GO TO 200 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 125 C *** CONSIDER SWITCHING MODELS *** IV(IRC) = 2 GO TO 130 C C *** ACCEPT STEP WITH DECREASED RADIUS *** C 125 IV(IRC) = 4 C C *** SET V(RADFAC) TO FLETCHER*S DECREASE FACTOR *** C 130 IV(XIRC) = IV(IRC) EMAX = V(GTSTEP) + V(FDIF) V(RADFAC) = HALF * RFAC1 IF (EMAX .LT. V(GTSTEP)) V(RADFAC) = RFAC1 * MAX(V(RDFCMN), + HALF * V(GTSTEP)/EMAX) C C *** DO FALSE CONVERGENCE TEST *** C 140 IF (V(RELDX) .LE. V(XFTOL)) GO TO 160 IV(IRC) = IV(XIRC) IF (V(F) .LT. V(F0)) GO TO 230 GO TO 300 C 160 IV(IRC) = 12 GO TO 310 C C *** HANDLE GOOD FUNCTION DECREASE *** C 200 IF (V(FDIF) .LT. (-V(TUNER3) * V(GTSTEP))) GO TO 260 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 260 IF (IV(RESTOR) .EQ. 1) GO TO 260 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) + V(RADFAC) = MAX(V(INCFAC), HALF*GTS/(GTS + V(FDIF))) IV(IRC) = 4 IF (V(STPPAR) .EQ. ZERO) GO TO 300 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 230 V(FLSTGD) = V(F) IV(MLSTGD) = IV(MODEL) CALL VCOPY(P, STLSTG, STEP) V(DSTSAV) = V(DSTNRM) IV(NFGCAL) = NFC V(PLSTGD) = V(PREDUC) V(GTSLST) = V(GTSTEP) GO TO 300 C C *** ACCEPT STEP WITH RADIUS UNCHANGED *** C 260 V(RADFAC) = ONE IV(IRC) = 3 GO TO 300 C C *** COME HERE FOR A RESTART AFTER CONVERGENCE *** C 290 IV(IRC) = IV(XIRC) IF (V(DSTSAV) .GE. ZERO) GO TO 310 IV(IRC) = 12 GO TO 310 C C *** PERFORM CONVERGENCE TESTS *** C 300 IV(XIRC) = IV(IRC) 310 IF (ABS(V(F)) .LT. V(AFCTOL)) IV(IRC) = 10 IF (HALF * V(FDIF) .GT. V(PREDUC)) GO TO 999 EMAX = 0.0 IF (ABS(V(F0)).GT.R1MACH(1)/V(RFCTOL)) + EMAX = V(RFCTOL) * ABS(V(F0)) IF (V(DSTNRM) .GT. V(LMAX0) .AND. V(PREDUC) .LE. EMAX) + IV(IRC) = 11 IF (V(DST0) .LT. ZERO) GO TO 320 I = 0 IF ((V(NREDUC) .GT. ZERO .AND. V(NREDUC) .LE. EMAX) .OR. + (V(NREDUC) .EQ. ZERO. AND. V(PREDUC) .EQ. ZERO)) I = 2 IF (V(STPPAR) .EQ. ZERO .AND. V(RELDX) .LE. V(XCTOL)) I = I + 1 IF (I .GT. 0) IV(IRC) = I + 6 C C *** CONSIDER RECOMPUTING STEP OF LENGTH V(LMAX0) FOR SINGULAR C *** CONVERGENCE TEST. C 320 IF (ABS(IV(IRC)-3) .GT. 1 .AND. IV(IRC) .NE. 12) GO TO 999 IF (V(DSTNRM) .GT. V(LMAX0)) GO TO 330 IF (V(PREDUC) .GE. EMAX) GO TO 999 IF (V(DST0) .LT. ZERO) GO TO 340 IF (HALF * V(DST0) .LE. V(LMAX0)) GO TO 999 GO TO 340 330 IF (HALF * V(DSTNRM) .LE. V(LMAX0)) GO TO 999 XMAX = V(LMAX0) / V(DSTNRM) IF (XMAX * (TWO - XMAX) * V(PREDUC) .GE. EMAX) GO TO 999 340 IF (V(NREDUC) .LT. ZERO) GO TO 370 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) IV(IRC) = 6 CALL VCOPY(P, STLSTG, STEP) GO TO 999 C C *** PERFORM SINGULAR CONVERGENCE TEST WITH RECOMPUTED V(PREDUC) *** C 360 V(GTSTEP) = V(GTSLST) V(DSTNRM) = ABS(V(DSTSAV)) CALL VCOPY(P, STEP, STLSTG) IV(IRC) = IV(XIRC) IF (V(DSTSAV) .LE. ZERO) IV(IRC) = 12 V(NREDUC) = -V(PREDUC) V(PREDUC) = V(PLSTGD) 370 IF (-V(NREDUC) .LE. V(RFCTOL) * ABS(V(F0))) IV(IRC) = 11 C 999 RETURN C C *** LAST CARD OF ASSESS FOLLOWS *** END *COVCLC SUBROUTINE COVCLC(COVIRC, D, IV, J, N, NN, P, R, V, X) C C LATEST REVISION - 03/15/90 (JRD) C C *** COMPUTE COVARIANCE MATRIX FOR NL2ITR (NL2SOL VERSION 2.2) *** C C *** LET K = ABS(IV(COVREQ). FOR K .LE. 2, A FINITE-DIFFERENCE C *** HESSIAN H IS COMPUTED (USING FUNC. AND GRAD. VALUES IF C *** IV(COVREQ) IS NONNEGATIVE, AND USING ONLY FUNC. VALUES IF C *** IV(COVREQ) IS NEGATIVE). FOR SCALE = 2*F(X) / MAX(1, N-P), C *** WHERE 2*F(X) IS THE RESIDUAL SUM OF SQUARES, COVCLC COMPUTES... C *** K = 0 OR 1... SCALE * H**-1 * (J**T * J) * H**-1. C *** K = 2... SCALE * H**-1. C *** K .GE. 3... SCALE * (J**T * J)**-1. C C C VARIABLE DECLARATIONS C C SCALAR ARGUMENTS INTEGER + COVIRC,N,NN,P C C ARRAY ARGUMENTS REAL + D(P),J(NN,P),R(N),V(1),X(P) INTEGER + IV(1) C C LOCAL SCALARS REAL + DEL,HALF,NEGPT5,ONE,T,TWO,WK,ZERO INTEGER + COV,COVMAT,COVREQ,DELTA,DELTA0,DLTFDC,F,FX,G,G1,GP,GSAVE1, + H,HC,HMI,HPI,HPM,I,IERR,IP1,IPIV0,IPIVI,IPIVK,IPIVOT,IRC, + K,KAGQT,KALM,KIND,KL,L,LMAT,M,MM1,MM1O2,MODE,NFGCAL,PP1O2, + QTR,QTR1,RD,RD1,RSAVE,SAVEI,STP0,STPI,STPM,SWITCH,TOOBIG, + W,W0,W1,WL,XMSAVE LOGICAL + HAVEJ C C EXTERNAL SUBROUTINES EXTERNAL LINVRT,LITVMU,LIVMUL,LSQRT,LTSQAR,QRFACT,VCOPY,VSCOPY C C INTRINSIC FUNCTIONS INTRINSIC ABS,MAX C C *** PARAMETER DECLARATIONS *** C C INTEGER COVIRC, IV(1), N, NN, P C REAL D(P), J(NN,P), R(N), V(1), X(P) C DIMENSION IV(*), V(*) C C *** LOCAL VARIABLES *** C C LOGICAL HAVEJ C INTEGER COV, GP, GSAVE1, G1, HC, HMI, HPI, HPM, I, IPIVI, IPIVK, C 1 IP1, IRC, K, KIND, KL, L, M, MM1, MM1O2, PP1O2, QTR1, C 2 RD1, STPI, STPM, STP0, WL, W0, W1 C REAL DEL, HALF, NEGPT5, ONE, T, TWO, WK, ZERO C C/ C *** EXTERNAL SUBROUTINES *** C C EXTERNAL LINVRT, LITVMU, LIVMUL, LSQRT, LTSQAR, QRFACT, C 1 VCOPY, VSCOPY C C LINVRT... INVERT LOWER TRIANGULAR MATRIX. C LITVMU... APPLY INVERSE-TRANSPOSE OF COMPACT LOWER TRIANG. MATRIX. C LIVMUL... APPLY INVERSE OF COMPACT LOWER TRIANG. MATRIX. C LSQRT.... COMPUTE CHOLESKY FACTOR OF (LOWER TRINAG. OF) A SYM. MATRIX. C LTSQAR... GIVEN LOWER TRIANG. MATRIX L, COMPUTE (L**T)*L. C QRFACT... COMPUTE QR DECOMPOSITION OF A MATRIX. C VCOPY.... COPY ONE VECTOR TO ANOTHER. C VSCOPY... SET ALL ELEMENTS OF A VECTOR TO A SCALAR. C C *** SUBSCRIPTS FOR IV AND V *** C C INTEGER COVMAT, COVREQ, DELTA, DELTA0, DLTFDC, F, FX, G, H, IERR, C 1 IPIVOT, IPIV0, KAGQT, KALM, LMAT, MODE, NFGCAL, QTR, C 2 RD, RSAVE, SAVEI, SWITCH, TOOBIG, W, XMSAVE C DATA HALF/0.5E0/, NEGPT5/-0.5E0/, ONE/1.0E0/, TWO/2.0E0/, + ZERO/0.0E0/ C DATA COVMAT/26/, COVREQ/15/, DELTA/50/, DELTA0/44/, + DLTFDC/40/, F/10/, FX/46/, G/28/, H/44/, IERR/32/, + IPIVOT/61/, IPIV0/60/, KAGQT/35/, KALM/36/, + LMAT/58/, MODE/38/, NFGCAL/7/, QTR/49/, + RD/51/, RSAVE/52/, SAVEI/54/, SWITCH/12/, + TOOBIG/2/, W/59/, XMSAVE/49/ C C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ C COV = IV(LMAT) C COVIRC = 4 KIND = IV(COVREQ) M = IV(MODE) IF (M .GT. 0) GO TO 10 IV(KAGQT) = -1 IF (IV(KALM) .GT. 0) IV(KALM) = 0 IF (ABS(KIND) .GE. 3) GO TO 300 V(FX) = V(F) K = IV(RSAVE) CALL VCOPY(N, V(K), R) 10 IF (M .GT. P) GO TO 200 IF (KIND .LT. 0) GO TO 100 C C *** COMPUTE FINITE-DIFFERENCE HESSIAN USING BOTH FUNCTION AND C *** GRADIENT VALUES. C GSAVE1 = IV(W) + P G1 = IV(G) IF (M .GT. 0) GO TO 15 C *** FIRST CALL ON COVCLC. SET GSAVE = G, TAKE FIRST STEP *** CALL VCOPY(P, V(GSAVE1), V(G1)) IV(SWITCH) = IV(NFGCAL) GO TO 80 C 15 DEL = V(DELTA) X(M) = V(XMSAVE) IF (IV(TOOBIG) .EQ. 0) GO TO 30 C C *** HANDLE OVERSIZE V(DELTA) *** C IF (DEL*X(M) .GT. ZERO) GO TO 20 C *** WE ALREADY TRIED SHRINKING V(DELTA), SO QUIT *** IV(COVMAT) = -2 GO TO 190 C C *** TRY SHRINKING V(DELTA) *** 20 DEL = NEGPT5 * DEL GO TO 90 C 30 COV = IV(LMAT) GP = G1 + P - 1 C C *** SET G = (G - GSAVE)/DEL *** C DO 40 I = G1, GP V(I) = (V(I) - V(GSAVE1)) / DEL GSAVE1 = GSAVE1 + 1 40 CONTINUE C C *** ADD G AS NEW COL. TO FINITE-DIFF. HESSIAN MATRIX *** C K = COV + M*(M-1)/2 L = K + M - 2 IF ( M .EQ. 1) GO TO 60 C C *** SET H(I,M) = 0.5 * (H(I,M) + G(I)) FOR I = 1 TO M-1 *** C DO 50 I = K, L V(I) = HALF * (V(I) + V(G1)) G1 = G1 + 1 50 CONTINUE C C *** ADD H(I,M) = G(I) FOR I = M TO P *** C 60 L = L + 1 DO 70 I = M, P V(L) = V(G1) L = L + I G1 = G1 + 1 70 CONTINUE C 80 M = M + 1 IV(MODE) = M IF (M .GT. P) GO TO 190 C C *** CHOOSE NEXT FINITE-DIFFERENCE STEP, RETURN TO GET G THERE *** C DEL = V(DELTA0) * MAX(ONE/D(M), ABS(X(M))) IF (X(M) .LT. ZERO) DEL = -DEL V(XMSAVE) = X(M) 90 X(M) = X(M) + DEL V(DELTA) = DEL COVIRC = 2 GO TO 999 C C *** COMPUTE FINITE-DIFFERENCE HESSIAN USING FUNCTION VALUES ONLY. C 100 STP0 = IV(W) + P - 1 MM1 = M - 1 MM1O2 = M*MM1/2 IF (M .GT. 0) GO TO 105 C *** FIRST CALL ON COVCLC. *** IV(SAVEI) = 0 GO TO 180 C 105 I = IV(SAVEI) IF (I .GT. 0) GO TO 160 IF (IV(TOOBIG) .EQ. 0) GO TO 120 C C *** HANDLE OVERSIZE STEP *** C STPM = STP0 + M DEL = V(STPM) IF (DEL*V(XMSAVE) .GT. ZERO) GO TO 110 C *** WE ALREADY TRIED SHRINKING THE STEP, SO QUIT *** IV(COVMAT) = -2 GO TO 999 C C *** TRY SHRINKING THE STEP *** 110 DEL = NEGPT5 * DEL X(M) = V(XMSAVE) + DEL V(STPM) = DEL COVIRC = 1 GO TO 999 C C *** SAVE F(X + STP(M)*E(M)) IN H(P,M) *** C 120 PP1O2 = P * (P-1) / 2 COV = IV(LMAT) HPM = COV + PP1O2 + MM1 V(HPM) = V(F) C C *** START COMPUTING ROW M OF THE FINITE-DIFFERENCE HESSIAN H. *** C HMI = COV + MM1O2 IF (MM1 .EQ. 0) GO TO 140 HPI = COV + PP1O2 DO 130 I = 1, MM1 V(HMI) = V(FX) - (V(F) + V(HPI)) HMI = HMI + 1 HPI = HPI + 1 130 CONTINUE 140 V(HMI) = V(F) - TWO*V(FX) C C *** COMPUTE FUNCTION VALUES NEEDED TO COMPLETE ROW M OF H. *** C I = 1 C 150 IV(SAVEI) = I STPI = STP0 + I V(DELTA) = X(I) X(I) = X(I) + V(STPI) IF (I .EQ. M) X(I) = V(XMSAVE) - V(STPI) COVIRC = 1 GO TO 999 C 160 X(I) = V(DELTA) IF (IV(TOOBIG) .EQ. 0) GO TO 170 C *** PUNT IN THE EVENT OF AN OVERSIZE STEP *** IV(COVMAT) = -2 GO TO 999 C C *** FINISH COMPUTING H(M,I) *** C 170 STPI = STP0 + I HMI = COV + MM1O2 + I - 1 STPM = STP0 + M V(HMI) = (V(HMI) + V(F)) / (V(STPI)*V(STPM)) I = I + 1 IF (I .LE. M) GO TO 150 IV(SAVEI) = 0 X(M) = V(XMSAVE) C 180 M = M + 1 IV(MODE) = M IF (M .GT. P) GO TO 190 C C *** PREPARE TO COMPUTE ROW M OF THE FINITE-DIFFERENCE HESSIAN H. C *** COMPUTE M-TH STEP SIZE STP(M), THEN RETURN TO OBTAIN C *** F(X + STP(M)*E(M)), WHERE E(M) = M-TH STD. UNIT VECTOR. C DEL = V(DLTFDC) * MAX(ONE/D(M), ABS(X(M))) IF (X(M) .LT. ZERO) DEL = -DEL V(XMSAVE) = X(M) X(M) = X(M) + DEL STPM = STP0 + M V(STPM) = DEL COVIRC = 1 GO TO 999 C C *** RESTORE R, V(F), ETC. *** C 190 K = IV(RSAVE) CALL VCOPY(N, R, V(K)) V(F) = V(FX) IF (KIND .LT. 0) GO TO 200 IV(NFGCAL) = IV(SWITCH) QTR1 = IV(QTR) CALL VCOPY(N, V(QTR1), R) IF (IV(COVMAT) .LT. 0) GO TO 999 COVIRC = 3 GO TO 999 C 200 COV = IV(LMAT) C C *** THE COMPLETE FINITE-DIFF. HESSIAN IS NOW STORED AT V(COV). *** C *** USE IT TO COMPUTE THE REQUESTED COVARIANCE MATRIX. *** C C *** COMPUTE CHOLESKY FACTOR C OF H = C*(C**T) *** C *** AND STORE IT AT V(HC). *** C HC = COV IF (ABS(KIND) .EQ. 2) GO TO 210 HC = ABS(IV(H)) IV(H) = -HC 210 CALL LSQRT(1, P, V(HC), V(COV), IRC) IV(COVMAT) = -1 IF (IRC .NE. 0) GO TO 999 C W1 = IV(W) + P IF (ABS(KIND) .GT. 1) GO TO 350 C C *** COVARIANCE = SCALE * H**-1 * (J**T * J) * H**-1 *** C CALL VSCOPY(P*(P+1)/2, V(COV), ZERO) HAVEJ = IV(KALM) .EQ. (-1) C *** HAVEJ = .TRUE. MEANS J IS IN ITS ORIGINAL FORM, WHILE C *** HAVEJ = .FALSE. MEANS QRFACT HAS BEEN APPLIED TO J. C M = P IF (HAVEJ) M = N W0 = W1 - 1 RD1 = IV(RD) DO 290 I = 1, M IF (HAVEJ) GO TO 240 C C *** SET W = IPIVOT * (ROW I OF R MATRIX FROM QRFACT). *** C CALL VSCOPY(P, V(W1), ZERO) IPIVI = IPIV0 + I L = W0 + IV(IPIVI) V(L) = V(RD1) RD1 = RD1 + 1 IF (I .EQ. P) GO TO 260 IP1 = I + 1 DO 230 K = IP1, P IPIVK = IPIV0 + K L = W0 + IV(IPIVK) V(L) = J(I,K) 230 CONTINUE GO TO 260 C C *** SET W = (ROW I OF J). *** C 240 L = W0 DO 250 K = 1, P L = L + 1 V(L) = J(I,K) 250 CONTINUE C C *** SET W = H**-1 * W. *** C 260 CALL LIVMUL(P, V(W1), V(HC), V(W1)) CALL LITVMU(P, V(W1), V(HC), V(W1)) C C *** ADD W * W**T TO COVARIANCE MATRIX. *** C KL = COV DO 280 K = 1, P L = W0 + K WK = V(L) DO 270 L = 1, K WL = W0 + L V(KL) = V(KL) + WK * V(WL) KL = KL + 1 270 CONTINUE 280 CONTINUE 290 CONTINUE GO TO 380 C C *** COVARIANCE = SCALE * (J**T * J)**-1. *** C 300 RD1 = IV(RD) IF (IV(KALM) .NE. (-1)) GO TO 310 C C *** APPLY QRFACT TO J *** C QTR1 = IV(QTR) CALL VCOPY(N, V(QTR1), R) W1 = IV(W) + P CALL QRFACT(NN, N, P, J, V(RD1), IV(IPIVOT), IV(IERR), 0, + V(W1)) IV(KALM) = -2 310 IV(COVMAT) = -1 IF (IV(IERR) .NE. 0) GO TO 999 COV = IV(LMAT) HC = ABS(IV(H)) IV(H) = -HC C C *** SET HC = (R MATRIX FROM QRFACT). *** C L = HC DO 340 I = 1, P IF (I .GT. 1) CALL VCOPY(I-1, V(L), J(1,I)) L = L + I - 1 V(L) = V(RD1) L = L + 1 RD1 = RD1 + 1 340 CONTINUE C C *** THE CHOLESKY FACTOR C OF THE UNSCALED INVERSE COVARIANCE MATRIX C *** (OR PERMUTATION THEREOF) IS STORED AT V(HC). C C *** SET C = C**-1. C 350 CALL LINVRT(P, V(HC), V(HC)) C C *** SET C = C**T * C. C CALL LTSQAR(P, V(HC), V(HC)) C IF (HC .EQ. COV) GO TO 380 C C *** C = PERMUTED, UNSCALED COVARIANCE. C *** SET COV = IPIVOT * C * IPIVOT**T. C DO 370 I = 1, P M = IPIV0 + I IPIVI = IV(M) KL = COV-1 + IPIVI*(IPIVI-1)/2 DO 360 K = 1, I M = IPIV0 + K IPIVK = IV(M) L = KL + IPIVK IF (IPIVK .GT. IPIVI) + L = L + (IPIVK-IPIVI)*(IPIVK+IPIVI-3)/2 V(L) = V(HC) HC = HC + 1 360 CONTINUE 370 CONTINUE C 380 IV(COVMAT) = COV C C *** APPLY SCALE FACTOR = (RESID. SUM OF SQUARES) / MAX(1,N-P). C T = V(F) / (HALF * MAX(1,N-P)) K = COV - 1 + P*(P+1)/2 DO 390 I = COV, K 390 V(I) = T * V(I) C 999 RETURN C *** LAST CARD OF COVCLC FOLLOWS *** END *DFAULT SUBROUTINE DFAULT(IV, V) C C LATEST REVISION - 03/15/90 (JRD) C C C C VARIABLE DECLARATIONS C C ARRAY ARGUMENTS REAL + V(45) INTEGER + IV(25) C C LOCAL SCALARS REAL + MACHEP,MEPCRT,ONE,SQTEPS,THREE INTEGER + AFCTOL,COSMIN,COVPRT,COVREQ,D0INIT,DECFAC,DELTA0,DFAC, + DINIT,DLTFDC,DLTFDJ,DTYPE,EPSLON,FUZZ,INCFAC,INITS,JTINIT, + LMAX0,MXFCAL,MXITER,OUTLEV,PARPRT,PHMNFC,PHMXFC,PRUNIT, + RDFCMN,RDFCMX,RFCTOL,RLIMIT,SOLPRT,STATPR,TUNER1,TUNER2, + TUNER3,TUNER4,TUNER5,X0PRT,XCTOL,XFTOL C C EXTERNAL FUNCTIONS REAL + RMDCON INTEGER + IMDCON EXTERNAL RMDCON,IMDCON C C INTRINSIC FUNCTIONS INTRINSIC MAX C C *** SUPPLY NL2SOL (VERSION 2.2) DEFAULT VALUES TO IV AND V *** C C INTEGER IV(25) C REAL V(45) C/+ C REAL MAX C/ C EXTERNAL IMDCON, RMDCON C INTEGER IMDCON C REAL RMDCON C C REAL MACHEP, MEPCRT, ONE, SQTEPS, THREE C C *** SUBSCRIPTS FOR IV AND V *** C C INTEGER AFCTOL, COSMIN, COVPRT, COVREQ, DECFAC, DELTA0, DFAC, C 1 DINIT, DLTFDC, DLTFDJ, DTYPE, D0INIT, EPSLON, FUZZ, C 2 INCFAC, INITS, JTINIT, LMAX0, MXFCAL, MXITER, OUTLEV, C 3 PARPRT, PHMNFC, PHMXFC, PRUNIT, RDFCMN, RDFCMX, C 4 RFCTOL, RLIMIT, SOLPRT, STATPR, TUNER1, TUNER2, TUNER3, C 5 TUNER4, TUNER5, XCTOL, XFTOL, X0PRT C DATA ONE/1.0E0/, THREE/3.0E0/ C C *** IV SUBSCRIPT VALUES *** C DATA COVPRT/14/, COVREQ/15/, DTYPE/16/, INITS/25/, + MXFCAL/17/, MXITER/18/, OUTLEV/19/, + PARPRT/20/, PRUNIT/21/, SOLPRT/22/, + STATPR/23/, X0PRT/24/ C C *** V SUBSCRIPT VALUES *** C DATA AFCTOL/31/, COSMIN/43/, DECFAC/22/, + DELTA0/44/, DFAC/41/, DINIT/38/, DLTFDC/40/, + DLTFDJ/36/, D0INIT/37/, EPSLON/19/, FUZZ/45/, + INCFAC/23/, JTINIT/39/, LMAX0/35/, PHMNFC/20/, + PHMXFC/21/, RDFCMN/24/, RDFCMX/25/, + RFCTOL/32/, RLIMIT/42/, TUNER1/26/, + TUNER2/27/, TUNER3/28/, TUNER4/29/, + TUNER5/30/, XCTOL/33/, XFTOL/34/ C C----------------------------------------------------------------------- C IV(1) = 12 IV(COVPRT) = 1 IV(COVREQ) = 1 IV(DTYPE) = 1 IV(INITS) = 0 IV(MXFCAL) = 200 IV(MXITER) = 150 IV(OUTLEV) = -1 IV(PARPRT) = 1 IV(PRUNIT) = IMDCON(1) IV(SOLPRT) = 1 IV(STATPR) = 1 IV(X0PRT) = 1 C MACHEP = RMDCON(3) V(AFCTOL) = 1.0E-20 IF (MACHEP .GT. 1.0E-10) V(AFCTOL) = MACHEP**2 V(COSMIN) = MAX(1.0E-6, 1.0E2 * MACHEP) V(DECFAC) = 0.5E0 SQTEPS = sqrt ( machep ) V(DELTA0) = SQTEPS V(DFAC) = 0.6E0 V(DINIT) = 0.0E0 MEPCRT = MACHEP ** (ONE/THREE) V(DLTFDC) = MEPCRT V(DLTFDJ) = SQTEPS V(D0INIT) = 1.0E0 V(EPSLON) = 0.1E0 V(FUZZ) = 1.5E0 V(INCFAC) = 2.0E0 V(JTINIT) = 1.0E-6 V(LMAX0) = 100.0E0 V(PHMNFC) = -0.1E0 V(PHMXFC) = 0.1E0 V(RDFCMN) = 0.1E0 V(RDFCMX) = 4.0E0 V(RFCTOL) = MAX(1.0E-10, MEPCRT**2) V(RLIMIT) = RMDCON(5) V(TUNER1) = 0.1E0 V(TUNER2) = 1.0E-4 V(TUNER3) = 0.75E0 V(TUNER4) = 0.5E0 V(TUNER5) = 0.75E0 V(XCTOL) = SQTEPS V(XFTOL) = 1.0E2 * MACHEP C RETURN C *** LAST CARD OF DFAULT FOLLOWS *** END *DOTPRD REAL FUNCTION DOTPRD(P, X, Y) C C LATEST REVISION - 03/15/90 (JRD) C C *** RETURN THE INNER PRODUCT OF THE P-VECTORS X AND Y. *** C C C VARIABLE DECLARATIONS C C SCALAR ARGUMENTS INTEGER + P C C ARRAY ARGUMENTS REAL + X(*),Y(*) C C LOCAL SCALARS REAL + ONE,SQTETA,T,ZERO INTEGER + I C C EXTERNAL FUNCTIONS REAL + RMDCON EXTERNAL RMDCON C C INTRINSIC FUNCTIONS INTRINSIC ABS,MAX C C INTEGER P C REAL X(*), Y(*) C C INTEGER I C REAL ONE, SQTETA, T, ZERO C/+ C REAL MAX, ABS C/ C EXTERNAL RMDCON C REAL RMDCON C C *** RMDCON(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 DATA ONE/1.0E0/, SQTETA/0.0E0/, ZERO/0.0E0/ C DOTPRD = ZERO IF (P .LE. 0) GO TO 999 IF (SQTETA .EQ. ZERO) SQTETA = RMDCON(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 DOTPRD = DOTPRD + X(I)*Y(I) 20 CONTINUE C 999 RETURN C *** LAST CARD OF DOTPRD FOLLOWS *** END *DUPDAT SUBROUTINE DUPDAT(D, IV, J, N, NN, P, V) C C LATEST REVISION - 03/15/90 (JRD) C C *** UPDATE SCALE VECTOR D FOR NL2ITR (NL2SOL VERSION 2.2) *** C C C VARIABLE DECLARATIONS C C SCALAR ARGUMENTS INTEGER + N,NN,P C C ARRAY ARGUMENTS REAL + D(P),J(NN,P),V(1) INTEGER + IV(1) C C LOCAL SCALARS REAL + SII,T,VDFAC,ZERO INTEGER + D0,DFAC,DTYPE,I,JTOL0,JTOLI,NITER,S,S1 C C EXTERNAL FUNCTIONS REAL + V2NORM EXTERNAL V2NORM C C INTRINSIC FUNCTIONS INTRINSIC MAX,SQRT C C *** PARAMETER DECLARATIONS *** C C INTEGER IV(1), N, NN, P C REAL D(P), J(NN,P), V(1) C DIMENSION IV(*), V(*) C C *** LOCAL VARIABLES *** C C INTEGER D0, I, JTOLI, S1 C REAL SII, T, VDFAC C C *** CONSTANTS *** C REAL ZERO C C/ C *** EXTERNAL FUNCTION *** C C EXTERNAL V2NORM C REAL V2NORM C C *** SUBSCRIPTS FOR IV AND V *** C C INTEGER DFAC, DTYPE, JTOL0, NITER, S DATA DFAC/41/, DTYPE/16/, JTOL0/86/, NITER/31/, S/53/ C DATA ZERO/0.0E0/ C C----------------------------------------------------------------------- C I = IV(DTYPE) IF (I .EQ. 1) GO TO 20 IF (IV(NITER) .GT. 0) GO TO 999 C 20 VDFAC = V(DFAC) D0 = JTOL0 + P S1 = IV(S) - 1 DO 30 I = 1, P S1 = S1 + I SII = V(S1) T = V2NORM(N, J(1,I)) IF (SII .GT. ZERO) T = SQRT(T*T + SII) JTOLI = JTOL0 + I D0 = D0 + 1 IF (T .LT. V(JTOLI)) T = MAX(V(D0), V(JTOLI)) D(I) = MAX(VDFAC*D(I), T) 30 CONTINUE C 999 RETURN C *** LAST CARD OF DUPDAT FOLLOWS *** END *GQTSTP SUBROUTINE GQTSTP(D, DIG, DIHDI, KA, L, P, STEP, V, W) C C LATEST REVISION - 03/15/90 (JRD) C C *** COMPUTE GOLDFELD-QUANDT-TROTTER STEP BY MORE-HEBDEN TECHNIQUE *** C *** (NL2SOL VERSION 2.2) *** C C C VARIABLE DECLARATIONS C C SCALAR ARGUMENTS INTEGER + KA,P C C ARRAY ARGUMENTS REAL + D(P),DIG(P),DIHDI(1),L(1),STEP(P),V(21),W(1) C C LOCAL SCALARS REAL + AKI,AKK,ALPHAK,DELTA,DGXFAC,DST,EPSFAC,EPSO6,FOUR,HALF,KAPPA, + LK,NEGONE,OLDPHI,ONE,P001,PHI,PHIMAX,PHIMIN,PSIFAC,RAD,ROOT, + SI,SIX,SK,SW,T,T1,THREE,TWO,TWOPSI,UK,WI,ZERO INTEGER + DGGDMX,DGNORM,DIAG,DIAG0,DST0,DSTNRM,DSTSAV,EMAX,EMIN, + EPSLON,GTSTEP,I,IM1,INC,IRC,J,K,K1,KALIM,LK0,NREDUC, + PHIPIN,PHMNFC,PHMXFC,PREDUC,Q,Q0,RAD0,RADIUS,STPPAR,UK0,X, + X0 LOGICAL + RESTRT C C EXTERNAL FUNCTIONS REAL + DOTPRD,LSVMIN,RMDCON,V2NORM EXTERNAL DOTPRD,LSVMIN,RMDCON,V2NORM C C EXTERNAL SUBROUTINES EXTERNAL LITVMU,LIVMUL,LSQRT C C INTRINSIC FUNCTIONS INTRINSIC ABS,MAX,MIN,SQRT C C *** PARAMETER DECLARATIONS *** C C INTEGER KA, P C REAL D(P), DIG(P), DIHDI(1), L(1), V(21), STEP(P), C 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 (GQTSTP ASSUMES DIG = D**-1 * G AND DIHDI = D**-1 * H * D**-1.) C IF G = 0, HOWEVER, STEP = 0 IS RETURNED (EVEN AT A SADDLE POINT). 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 INTIIAL 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. TO COMPUTE STEP FROM A SADDLE C POINT (WHERE THE TRUE GRADIENT VANISHES AND H HAS A NEGATIVE C EIGENVALUE), A NONZERO G WITH SMALL COMPONENTS SHOULD BE PASSED. C C *** APPLICATION AND USAGE RESTRICTIONS *** C C THIS ROUTINE IS CALLED AS PART OF THE NL2SOL (NONLINEAR LEAST- C SQUARES) PACKAGE (REF. 1), BUT IT COULD BE USED IN SOLVING ANY C UNCONSTRAINED MINIMIZATION PROBLEM. C C *** ALGORITHM NOTES *** C C THE DESIRED G-Q-T STEP (REF. 2, 3, 4) 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. 6. 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). C C *** FUNCTIONS AND SUBROUTINES CALLED *** C C DOTPRD - RETURNS INNER PRODUCT OF TWO VECTORS. C LITVMU - APPLIES INVERSE-TRANSPOSE OF COMPACT LOWER TRIANG. MATRIX. C LIVMUL - APPLIES INVERSE OF COMPACT LOWER TRIANG. MATRIX. C LSQRT - FINDS CHOLESKY FACTOR (OF COMPACTLY STORED LOWER TRIANG.). C LSVMIN - RETURNS APPROX. TO MIN. SING. VALUE OF LOWER TRIANG. MATRIX. C RMDCON - RETURNS MACHINE-DEPENDENT CONSTANTS. C V2NORM - RETURNS 2-NORM OF A VECTOR. C C *** REFERENCES *** C C 1. DENNIS, J.E., GAY, D.M., AND WELSCH, R.E. (1980), AN ADAPTIVE C NONLINEAR LEAST-SQUARES ALGORITHM, (SUBMITTED TO ACM C TRANS. MATH. SOFTWARE). C 2. GAY, D.M. (1979), COMPUTING OPTIMAL ELLIPTICALLY CONSTRAINED C STEPS, MRC TECH. SUMMARY REPORT NO. 2013, MATH. RESEARCH C CENTER, UNIV. OF WISCONSIN-MADISON. 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. 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 C LOGICAL RESTRT C INTEGER DGGDMX, DIAG, DIAG0, DSTSAV, EMAX, EMIN, I, IM1, INC, IRC, C 1 J, K, KALIM, K1, LK0, PHIPIN, Q, Q0, UK0, X, X0 C REAL ALPHAK, AKI, AKK, DELTA, DST, EPSO6, LK, C 1 OLDPHI, PHI, PHIMAX, PHIMIN, PSIFAC, RAD, C 2 ROOT, SI, SK, SW, T, TWOPSI, T1, UK, WI C C *** CONSTANTS *** C REAL DGXFAC, EPSFAC, FOUR, HALF, KAPPA, NEGONE, ONE, C 1 P001, SIX, THREE, TWO, ZERO C C/ C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C C EXTERNAL DOTPRD, LITVMU, LIVMUL, LSQRT, LSVMIN, RMDCON, V2NORM C REAL DOTPRD, LSVMIN, RMDCON, V2NORM C C *** SUBSCRIPTS FOR V *** C C INTEGER DGNORM, DSTNRM, DST0, EPSLON, GTSTEP, STPPAR, NREDUC, C 1 PHMNFC, PHMXFC, PREDUC, RADIUS, RAD0 DATA DGNORM/1/, DSTNRM/2/, DST0/3/, EPSLON/19/, + GTSTEP/4/, NREDUC/6/, PHMNFC/20/, + PHMXFC/21/, PREDUC/7/, RADIUS/8/, + RAD0/9/, STPPAR/5/ C DATA DGXFAC/0.0E0/, EPSFAC/50.0E0/, FOUR/4.0E0/, HALF/0.5E0/, + KAPPA/2.0E0/, NEGONE/-1.0E0/, ONE/1.0E0/, P001/1.0E-3/, + SIX/6.0E0/, THREE/3.0E0/, TWO/2.0E0/, ZERO/0.0E0/ C C *** BODY *** 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. UK = 0.0E0 PHI = 0.0E0 DST = 0.0E0 ALPHAK = 0.0E0 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 RAD = V(RADIUS) C *** PHITOL = MAX. ERROR ALLOWED IN DST = V(DSTNRM) = 2-NORM OF C *** D*STEP. PHIMAX = V(PHMXFC) * RAD PHIMIN = V(PHMNFC) * RAD C *** EPSO6 AND PSIFAC ARE USED IN CHECKING FOR THE SPECIAL CASE C *** OF (NEARLY) SINGULAR H + ALPHA*D**2 (SEE REF. 2). PSIFAC = TWO * V(EPSLON) / (THREE * (FOUR * (V(PHMNFC) + ONE) * + (KAPPA + ONE) + KAPPA + TWO) * RAD**2) 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 EPSO6 = V(EPSLON)/SIX IRC = 0 RESTRT = .FALSE. KALIM = KA + 50 C C *** START OR RESTART, DEPENDING ON KA *** C IF (KA .GE. 0) GO TO 310 C C *** FRESH START *** C K = 0 UK = NEGONE KA = 0 KALIM = 50 C C *** STORE DIAG(DIHDI) IN W(DIAG0+1),...,W(DIAG0+P) *** C J = 0 DO 20 I = 1, P J = J + I K1 = DIAG0 + I W(K1) = DIHDI(J) 20 CONTINUE C C *** DETERMINE W(DGGDMX), THE LARGEST ELEMENT OF DIHDI *** C T1 = ZERO J = P * (P + 1) / 2 DO 30 I = 1, J T = ABS(DIHDI(I)) IF (T1 .LT. T) T1 = T 30 CONTINUE W(DGGDMX) = T1 C C *** TRY ALPHA = 0 *** C 40 CALL LSQRT(1, P, L, DIHDI, IRC) IF (IRC .EQ. 0) GO TO 60 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 50 I = 1, IRC 50 W(I) = ZERO W(IRC) = ONE CALL LITVMU(IRC, W, L, W) T1 = V2NORM(IRC, W) LK = -T / T1 / T1 V(DST0) = -LK IF (RESTRT) GO TO 210 V(NREDUC) = ZERO GO TO 70 C C *** POSITIVE DEFINITE H -- COMPUTE UNMODIFIED NEWTON STEP. *** 60 LK = ZERO CALL LIVMUL(P, W(Q), L, DIG) V(NREDUC) = HALF * DOTPRD(P, W(Q), W(Q)) CALL LITVMU(P, W(Q), L, W(Q)) DST = V2NORM(P, W(Q)) V(DST0) = DST PHI = DST - RAD IF (PHI .LE. PHIMAX) GO TO 280 IF (RESTRT) GO TO 210 C C *** PREPARE TO COMPUTE GERSCHGORIN ESTIMATES OF LARGEST (AND C *** SMALLEST) EIGENVALUES. *** C 70 V(DGNORM) = V2NORM(P, DIG) IF (V(DGNORM) .EQ. ZERO) GO TO 450 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) 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 C IF (IRC .NE. 0) GO TO 210 C C *** COMPUTE L0 FOR POSITIVE DEFINITE H *** C CALL LIVMUL(P, W, L, W(Q)) T = V2NORM(P, W) W(PHIPIN) = DST / 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) + ALPHAK = UK * MAX(P001, SQRT(LK/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 LSQRT(1, P, L, DIHDI, IRC) IF (IRC .EQ. 0) GO TO 250 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 LITVMU(IRC, W, L, W) T1 = V2NORM(IRC, W) LK = ALPHAK - T/T1/T1 V(DST0) = -LK GO TO 210 C C *** ALPHAK MAKES (D**-1)*H*(D**-1) POSITIVE DEFINITE. C *** COMPUTE Q = -D*STEP, CHECK FOR CONVERGENCE. *** C 250 CALL LIVMUL(P, W(Q), L, DIG) CALL LITVMU(P, W(Q), L, W(Q)) DST = V2NORM(P, W(Q)) PHI = DST - RAD IF (PHI .LE. PHIMAX .AND. PHI .GE. PHIMIN) GO TO 290 IF (PHI .EQ. OLDPHI) GO TO 290 OLDPHI = PHI IF (PHI .GT. ZERO) GO TO 260 C *** CHECK FOR THE SPECIAL CASE OF H + ALPHA*D**2 (NEARLY) C *** SINGULAR. DELTA IS .GE. THE SMALLEST EIGENVALUE OF C *** (D**-1)*H*(D**-1) + ALPHAK*I. IF (V(DST0) .GT. ZERO) GO TO 260 DELTA = ALPHAK + V(DST0) TWOPSI = ALPHAK*DST*DST + DOTPRD(P, DIG, W(Q)) IF (DELTA .LT. PSIFAC*TWOPSI) GO TO 270 C C *** UNACCEPTABLE ALPHAK -- UPDATE LK, UK, ALPHAK *** C 260 IF (KA .GE. KALIM) GO TO 290 CALL LIVMUL(P, W, L, W(Q)) T1 = V2NORM(P, W) C *** THE FOLLOWING MIN IS NECESSARY BECAUSE OF RESTARTS *** IF (PHI .LT. ZERO) UK = MIN(UK, ALPHAK) ALPHAK = ALPHAK + (PHI/T1) * (DST/T1) * (DST/RAD) LK = MAX(LK, ALPHAK) GO TO 210 C C *** DECIDE HOW TO HANDLE (NEARLY) SINGULAR H + ALPHA*D**2 *** C C *** IF NOT YET AVAILABLE, OBTAIN MACHINE DEPENDENT VALUE DGXFAC. 270 IF (DGXFAC .EQ. ZERO) DGXFAC = EPSFAC * RMDCON(3) C C *** NOW DECIDE. *** IF (DELTA .GT. DGXFAC*W(DGGDMX)) GO TO 350 C *** DELTA IS SO SMALL WE CANNOT HANDLE THE SPECIAL CASE IN C *** THE AVAILABLE ARITHMETIC. ACCEPT STEP AS IT IS. GO TO 290 C C *** ACCEPTABLE STEP ON FIRST TRY *** C 280 ALPHAK = ZERO C C *** SUCCESSFUL STEP IN GENERAL. COMPUTE STEP = -(D**-1)*Q *** C 290 DO 300 I = 1, P J = Q0 + I STEP(I) = -W(J)/D(I) 300 CONTINUE V(GTSTEP) = -DOTPRD(P, DIG, W(Q)) V(PREDUC) = HALF * (ABS(ALPHAK)*DST*DST - V(GTSTEP)) GO TO 430 C C C *** RESTART WITH NEW RADIUS *** C 310 IF (V(DST0) .LE. ZERO .OR. V(DST0) - RAD .GT. PHIMAX) GO TO 330 C C *** PREPARE TO RETURN NEWTON STEP *** C RESTRT = .TRUE. KA = KA + 1 K = 0 DO 320 I = 1, P K = K + I J = DIAG0 + I DIHDI(K) = W(J) 320 CONTINUE UK = NEGONE GO TO 40 C 330 IF (KA .EQ. 0) GO TO 60 C DST = W(DSTSAV) ALPHAK = ABS(V(STPPAR)) PHI = DST - RAD T = V(DGNORM)/RAD IF (RAD .GT. V(RAD0)) GO TO 340 C C *** SMALLER RADIUS *** UK = T - W(EMIN) 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 260 C C *** BIGGER RADIUS *** 340 UK = T - W(EMIN) IF (ALPHAK .GT. ZERO) UK = MIN(UK, W(UK0)) LK = MAX(ZERO, -V(DST0), T - W(EMAX)) IF (V(DST0) .GT. ZERO) LK = MAX(LK, (V(DST0)-RAD)*W(PHIPIN)) GO TO 260 C C *** HANDLE (NEARLY) SINGULAR H + ALPHA*D**2 *** C C *** NEGATE ALPHAK TO INDICATE SPECIAL CASE *** 350 ALPHAK = -ALPHAK C *** ALLOCATE STORAGE FOR SCRATCH VECTOR X *** X0 = Q0 + P X = X0 + 1 C C *** USE INVERSE POWER METHOD WITH START FROM LSVMIN TO OBTAIN C *** APPROXIMATE EIGENVECTOR CORRESPONDING TO SMALLEST EIGENVALUE C *** OF (D**-1)*H*(D**-1). C DELTA = KAPPA*DELTA T = LSVMIN(P, L, W(X), W) C K = 0 C *** NORMALIZE W *** 360 DO 370 I = 1, P 370 W(I) = T*W(I) C *** COMPLETE CURRENT INV. POWER ITER. -- REPLACE W BY (L**-T)*W. CALL LITVMU(P, W, L, W) T1 = ONE/V2NORM(P, W) T = T1*T IF (T .LE. DELTA) GO TO 390 IF (K .GT. 30) GO TO 290 K = K + 1 C *** START NEXT INV. POWER ITER. BY STORING NORMALIZED W IN X. DO 380 I = 1, P J = X0 + I W(J) = T1*W(I) 380 CONTINUE C *** COMPUTE W = (L**-1)*X. CALL LIVMUL(P, W, L, W(X)) T = ONE/V2NORM(P, W) GO TO 360 C 390 DO 400 I = 1, P 400 W(I) = T1*W(I) C C *** NOW W IS THE DESIRED APPROXIMATE (UNIT) EIGENVECTOR AND C *** T*X = ((D**-1)*H*(D**-1) + ALPHAK*I)*W. C SW = DOTPRD(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 *** ACCEPT CURRENT STEP IF ADDING SI*W WOULD LEAD TO A C *** FURTHER RELATIVE REDUCTION IN PSI OF LESS THAN V(EPSLON)/3. V(PREDUC) = HALF*TWOPSI T1 = ZERO T = SI*(ALPHAK*SW - HALF*SI*(ALPHAK + T*DOTPRD(P,W(X),W))) IF (T .LT. EPSO6*TWOPSI) GO TO 410 V(PREDUC) = V(PREDUC) + T DST = RAD T1 = -SI 410 DO 420 I = 1, P J = Q0 + I W(J) = T1*W(I) - W(J) STEP(I) = W(J) / D(I) 420 CONTINUE V(GTSTEP) = DOTPRD(P, DIG, W(Q)) C C *** SAVE VALUES FOR USE IN A POSSIBLE RESTART *** C 430 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 440 I = 1, P J = J + I K = DIAG0 + I DIHDI(J) = W(K) 440 CONTINUE GO TO 999 C C *** SPECIAL CASE -- G = 0 *** C 450 V(STPPAR) = ZERO V(PREDUC) = ZERO V(DSTNRM) = ZERO V(GTSTEP) = ZERO DO 460 I = 1, P 460 STEP(I) = ZERO C 999 RETURN C C *** LAST CARD OF GQTSTP FOLLOWS *** END *IMDCON INTEGER FUNCTION IMDCON(K) C C LATEST REVISION - 03/15/90 (JRD) C C *** RETURN INTEGER MACHINE-DEPENDENT CONSTANTS *** C C *** K = 1 MEANS RETURN STANDARD OUTPUT UNIT NUMBER. *** C *** K = 2 MEANS RETURN ALTERNATE OUTPUT UNIT NUMBER. *** C *** K = 3 MEANS RETURN INPUT UNIT NUMBER. *** C (NOTE -- K = 2, 3 ARE USED ONLY BY TEST PROGRAMS.) C C C C VARIABLE DECLARATIONS C C SCALAR ARGUMENTS INTEGER + K C C LOCAL ARRAYS INTEGER + MDCON(3) C C EXTERNAL FUNCTIONS INTEGER + I1MACH EXTERNAL I1MACH C MDCON(1) = I1MACH(2) MDCON(2) = I1MACH(3) MDCON(3) = I1MACH(1) C IMDCON = MDCON(K) RETURN C *** LAST CARD OF IMDCON FOLLOWS *** END *ITSMRY SUBROUTINE ITSMRY(D, IV, P, V, X) C C LATEST REVISION - 03/15/90 (JRD) C C *** PRINT NL2SOL (VERSION 2.2) ITERATION SUMMARY *** C C C VARIABLE DECLARATIONS C C SCALAR ARGUMENTS INTEGER + P C C ARRAY ARGUMENTS REAL + D(P),V(*),X(P) INTEGER + IV(*) C C LOCAL SCALARS REAL + NRELDF,OLDF,PRELDF,RELDF,ZERO INTEGER + COV1,COVMAT,COVPRT,COVREQ,DSTNRM,F,F0,FDIF,G,G1,I,I1,ICH, + II,IV1,J,M,NEEDHD,NF,NFCALL,NFCOV,NG,NGCALL,NGCOV,NITER, + NREDUC,OL,OUTLEV,PREDUC,PRNTIT,PRUNIT,PU,RELDX,SIZE, + SOLPRT,STATPR,STPPAR,SUSED,X0PRT C C LOCAL ARRAYS CHARACTER + MODEL1(3,6)*1,MODEL2(4,6)*1 C C INTRINSIC FUNCTIONS INTRINSIC ABS C C *** PARAMETER DECLARATIONS *** C C INTEGER IV(1), P C REAL D(P), V(1), X(P) C DIMENSION IV(*), V(*) C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C *** LOCAL VARIABLES *** C C INTEGER COV1, G1, I, II, IV1, I1, J, M, NF, NG, OL, PU C CHARACTER*1 MODEL1(3, 6), MODEL2(4, 6) C REAL NRELDF, OLDF, PRELDF, RELDF, ZERO C C/ C *** NO EXTERNAL FUNCTIONS OR SUBROUTINES *** C C *** SUBSCRIPTS FOR IV AND V *** C C INTEGER COVMAT, COVPRT, COVREQ, DSTNRM, F, FDIF, F0, G, C 1 NEEDHD, NFCALL, NFCOV, NGCOV, NGCALL, NITER, NREDUC, C 2 OUTLEV, PREDUC, PRNTIT, PRUNIT, RELDX, SIZE, SOLPRT, C 3 STATPR, STPPAR, SUSED, X0PRT C C *** IV SUBSCRIPT VALUES *** C DATA COVMAT/26/, COVPRT/14/, G/28/, COVREQ/15/, + NEEDHD/39/, NFCALL/6/, NFCOV/40/, NGCOV/41/, + NGCALL/30/, NITER/31/, OUTLEV/19/, PRNTIT/48/, + PRUNIT/21/, SOLPRT/22/, STATPR/23/, SUSED/57/, + X0PRT/24/ C C *** V SUBSCRIPT VALUES *** C DATA DSTNRM/2/, F/10/, F0/13/, FDIF/11/, NREDUC/6/, + PREDUC/7/, RELDX/17/, SIZE/47/, STPPAR/5/ C DATA MODEL1(1, 1), MODEL1(2, 1), MODEL1(3, 1) + / ' ', ' ', ' ' / DATA MODEL1(1, 2), MODEL1(2, 2), MODEL1(3, 2) + / ' ', ' ', ' ' / DATA MODEL1(1, 3), MODEL1(2, 3), MODEL1(3, 3) + / ' ', ' ', ' ' / DATA MODEL1(1, 4), MODEL1(2, 4), MODEL1(3, 4) + / ' ', ' ', ' ' / DATA MODEL1(1, 5), MODEL1(2, 5), MODEL1(3, 5) + / ' ', 'G', ' ' / DATA MODEL1(1, 6), MODEL1(2, 6), MODEL1(3, 6) + / ' ', 'S', ' ' / DATA MODEL2(1, 1), MODEL2(2, 1), MODEL2(3, 1), MODEL2(4, 1) + / ' ', 'G', ' ', ' ' / DATA MODEL2(1, 2), MODEL2(2, 2), MODEL2(3, 2), MODEL2(4, 2) + / ' ', 'S', ' ', ' ' / DATA MODEL2(1, 3), MODEL2(2, 3), MODEL2(3, 3), MODEL2(4, 3) + / 'G', '-', 'S', ' ' / DATA MODEL2(1, 4), MODEL2(2, 4), MODEL2(3, 4), MODEL2(4, 4) + / 'S', '-', 'G', ' ' / DATA MODEL2(1, 5), MODEL2(2, 5), MODEL2(3, 5), MODEL2(4, 5) + / '-', 'S', '-', 'G' / DATA MODEL2(1, 6), MODEL2(2, 6), MODEL2(3, 6), MODEL2(4, 6) + / '-', 'G', '-', 'S' / DATA ZERO/0.0E0/ C C----------------------------------------------------------------------- C PU = IV(PRUNIT) IF (PU .EQ. 0) GO TO 999 IV1 = IV(1) OL = IV(OUTLEV) IF (IV1 .LT. 2 .OR. IV1 .GT. 15) GO TO 140 IF (OL .EQ. 0) GO TO 20 IF (IV1 .GE. 12) GO TO 20 IF (IV1 .GE. 10 .AND. IV(PRNTIT) .EQ. 0) GO TO 20 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 = V(F0) IF (OLDF .LE. ZERO) GO TO 12 RELDF = V(FDIF) / OLDF PRELDF = V(PREDUC) / OLDF 12 IF (OL .GT. 0) GO TO 15 C C *** PRINT SHORT SUMMARY LINE *** C IF (IV(NEEDHD) .EQ. 1) WRITE(PU, 1010) 1010 FORMAT(12H0 IT NF,6X,'F',8X,5HRELDF,6X,6HPRELDF,5X,5HRELDX) IV(NEEDHD) = 0 WRITE(PU,1017) IV(NITER), NF, V(F), RELDF, PRELDF, V(RELDX) GO TO 20 C C *** PRINT LONG SUMMARY LINE *** C 15 IF (IV(NEEDHD) .EQ. 1) WRITE(PU,1015) 1015 FORMAT(12H0 IT NF,6X,'F',8X,5HRELDF,6X,6HPRELDF,5X,5HRELDX, + 4X,15HMODEL STPPAR,6X,4HSIZE,6X,6HD*STEP,5X,7HNPRELDF) IV(NEEDHD) = 0 M = IV(SUSED) NRELDF = ZERO IF (OLDF .GT. ZERO) NRELDF = V(NREDUC) / OLDF WRITE(PU,1017) IV(NITER), NF, V(F), RELDF, PRELDF, V(RELDX), + (MODEL1(ICH, M), ICH = 1, 3), + (MODEL2(ICH, M), ICH = 1, 4), + V(STPPAR), V(SIZE), V(DSTNRM), NRELDF 1017 FORMAT(1X,I5,I6,4E11.3,7A1,4E11.3) C 20 GO TO (999,999,30,35,40,45,50,60,70,80,90,150,110,120,130), IV1 C 30 WRITE(PU,1030) 1030 FORMAT(26H0***** X-CONVERGENCE *****) GO TO 180 C 35 WRITE(PU,1035) 1035 FORMAT(42H0***** RELATIVE FUNCTION CONVERGENCE *****) GO TO 180 C 40 WRITE(PU,1040) 1040 FORMAT(49H0***** X- AND RELATIVE FUNCTION CONVERGENCE *****) GO TO 180 C 45 WRITE(PU,1045) 1045 FORMAT(42H0***** ABSOLUTE FUNCTION CONVERGENCE *****) GO TO 180 C 50 WRITE(PU,1050) 1050 FORMAT(33H0***** SINGULAR CONVERGENCE *****) GO TO 180 C 60 WRITE(PU,1060) 1060 FORMAT(30H0***** FALSE CONVERGENCE *****) GO TO 180 C 70 WRITE(PU,1070) 1070 FORMAT(38H0***** FUNCTION EVALUATION LIMIT *****) GO TO 180 C 80 WRITE(PU,1080) 1080 FORMAT(28H0***** ITERATION LIMIT *****) GO TO 180 C 90 WRITE(PU,1090) 1090 FORMAT(18H0***** STOPX *****) GO TO 180 C 110 WRITE(PU,1100) 1100 FORMAT(45H0***** INITIAL SUM OF SQUARES OVERFLOWS *****) C GO TO 150 C 120 WRITE(PU,1120) 1120 FORMAT(37H0***** BAD PARAMETERS TO ASSESS *****) GO TO 999 C 130 continue write ( pu, * ) ' ' WRITE(PU,1130) 1130 FORMAT('***** J COULD NOT BE COMPUTED *****') IF (IV(NITER) .GT. 0) GO TO 190 GO TO 150 C 140 write ( pu, * ) ' ' WRITE(PU,1140) IV1 1140 FORMAT('***** IV(1) =',I5,' *****') GO TO 999 C C *** INITIAL CALL ON ITSMRY *** C 150 IF (IV(X0PRT) .NE. 0) then write ( pu, * ) ' ' WRITE(PU,1150) (I, X(I), D(I), I = 1, P) end if 1150 FORMAT(' I INITIAL X(I) D(I)'//(1X,I5,E17.6,E14.3)) IF (IV1 .GE. 13) GO TO 999 IV(NEEDHD) = 0 IV(PRNTIT) = 0 IF (OL .EQ. 0) GO TO 999 IF (OL .LT. 0) WRITE(PU,1010) IF (OL .GT. 0) WRITE(PU,1015) write ( pu, * ) ' ' WRITE(PU,1160) V(F) 1160 FORMAT(' 0 1',E11.3,11X,E11.3) GO TO 999 C C *** PRINT VARIOUS INFORMATION REQUESTED ON SOLUTION *** C 180 IV(NEEDHD) = 1 IF (IV(STATPR) .EQ. 0) GO TO 190 OLDF = V(F0) PRELDF = ZERO NRELDF = ZERO IF (OLDF .LE. ZERO) GO TO 185 PRELDF = V(PREDUC) / OLDF NRELDF = V(NREDUC) / OLDF 185 NF = IV(NFCALL) - IV(NFCOV) NG = IV(NGCALL) - IV(NGCOV) WRITE(PU,1180) V(F), V(RELDX), NF, NG, PRELDF, NRELDF 1180 FORMAT(9H0FUNCTION,E17.6,8H RELDX,E20.6/12H FUNC. EVALS, + I8,9X,'GRAD. EVALS',I8/' PRELDF',E19.6,3X,'NPRELDF',E18.6) C IF (IV(NFCOV) .GT. 0) then write ( pu, * ) ' ' WRITE(PU,1185) IV(NFCOV) end if 1185 FORMAT(1x,I4,' EXTRA FUNC. EVALS FOR COVARIANCE.') IF (IV(NGCOV) .GT. 0) WRITE(PU,1186) IV(NGCOV) 1186 FORMAT(1X,I4,' EXTRA GRAD. EVALS FOR COVARIANCE.') 190 IF (IV(SOLPRT) .EQ. 0) GO TO 210 IV(NEEDHD) = 1 G1 = IV(G) write ( pu, * ) ' ' WRITE(PU,1190) 1190 FORMAT(' I FINAL X(I)',8X,'D(I)',10X,'G(I)'/) DO 200 I = 1, P WRITE(PU,1200) I, X(I), D(I), V(G1) G1 = G1 + 1 200 CONTINUE 1200 FORMAT(1X,I5,E17.6,2E14.3) C 210 IF (IV(COVPRT) .EQ. 0) GO TO 999 COV1 = IV(COVMAT) IV(NEEDHD) = 1 IF (COV1) 220, 230, 240 220 IF (-1 .EQ. COV1) then write ( pu, * ) ' ' WRITE(PU,1220) end if 1220 FORMAT('++++++ INDEFINITE COVARIANCE MATRIX ++++++') IF (-2 .EQ. COV1) then write ( pu, * ) ' ' WRITE(PU,1225) end if 1225 FORMAT('++++++ OVERSIZE STEPS IN COMPUTING COVARIANCE +++++') GO TO 999 230 continue write ( pu, * ) ' ' WRITE(PU,1230) 1230 FORMAT('++++++ COVARIANCE MATRIX NOT COMPUTED ++++++') GO TO 999 240 I = ABS(IV(COVREQ)) IF (I .LE. 1) then write ( pu, * ) ' ' WRITE(PU,1241) write ( pu, * ) ' ' end if 1241 FORMAT('COVARIANCE = SCALE * H**-1 * (J**T * J) * H**-1') IF (I .EQ. 2) then write ( pu, * ) ' ' WRITE(PU,1242) write ( pu, * ) ' ' end if 1242 FORMAT('COVARIANCE = SCALE * H**-1') IF (I .GE. 3) then write ( pu, * ) ' ' WRITE(PU,1243) write ( pu, * ) ' ' end if 1243 FORMAT('COVARIANCE = SCALE * (J**T * J)**-1') II = COV1 - 1 IF (OL .LE. 0) GO TO 260 DO 250 I = 1, P I1 = II + 1 II = II + I WRITE(PU,1250) I, (V(J), J = I1, II) 250 CONTINUE 1250 FORMAT(4H ROW,I3,2X,9E12.4/(9X,9E12.4)) GO TO 999 C 260 DO 270 I = 1, P I1 = II + 1 II = II + I WRITE(PU,1270) I, (V(J), J = I1, II) 270 CONTINUE 1270 FORMAT(4H ROW,I3,2X,5E12.4/(9X,5E12.4)) C 999 RETURN END SUBROUTINE LINVRT(N, LIN, L) C C *** COMPUTE LIN = L**-1, BOTH N X N LOWER TRIANG. STORED *** C *** COMPACTLY BY ROWS. LIN AND L MAY SHARE THE SAME STORAGE. *** C C C VARIABLE DECLARATIONS C C SCALAR ARGUMENTS INTEGER + N C C ARRAY ARGUMENTS REAL + L(*),LIN(*) C C LOCAL SCALARS REAL + ONE,T,ZERO INTEGER + I,II,IM1,J0,J1,JJ,K,K0,NP1 C C *** PARAMETERS *** C C INTEGER N C REAL L(*), LIN(*) C DIMENSION L(N*(N+1)/2), LIN(N*(N+1)/2) C C *** LOCAL VARIABLES *** C C INTEGER I, II, IM1, JJ, J0, J1, K, K0, NP1 C REAL ONE, T, ZERO DATA ONE/1.0E0/, ZERO/0.0E0/ C C *** BODY *** C NP1 = N + 1 J0 = N*(NP1)/2 DO 30 II = 1, N I = NP1 - II LIN(J0) = ONE/L(J0) IF (I .LE. 1) GO TO 999 J1 = J0 IM1 = I - 1 DO 20 JJ = 1, IM1 T = ZERO J0 = J1 K0 = J1 - JJ DO 10 K = 1, JJ T = T - L(K0)*LIN(J0) J0 = J0 - 1 K0 = K0 + K - I 10 CONTINUE LIN(J0) = T/L(K0) 20 CONTINUE J0 = J0 - 1 30 CONTINUE 999 RETURN C *** LAST CARD OF LINVRT FOLLOWS *** END *LITVMU SUBROUTINE LITVMU(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 C C VARIABLE DECLARATIONS C C SCALAR ARGUMENTS INTEGER + N C C ARRAY ARGUMENTS REAL + L(1),X(N),Y(N) C C LOCAL SCALARS REAL + XI,ZERO INTEGER + I,I0,II,IJ,IM1,J,NP1 C DATA ZERO/0.0E0/ 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 LITVMU FOLLOWS *** END *LIVMUL SUBROUTINE LIVMUL(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 C C VARIABLE DECLARATIONS C C SCALAR ARGUMENTS INTEGER + N C C ARRAY ARGUMENTS REAL + L(1),X(N),Y(N) C C LOCAL SCALARS REAL + T,ZERO INTEGER + I,J,K C C EXTERNAL FUNCTIONS REAL + DOTPRD EXTERNAL DOTPRD C DATA ZERO/0.0E0/ 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 = DOTPRD(I-1, L(J+1), X) J = J + I X(I) = (Y(I) - T)/L(J) 30 CONTINUE 999 RETURN C *** LAST CARD OF LIVMUL FOLLOWS *** END *LMSTEP SUBROUTINE LMSTEP(D, G, IERR, IPIVOT, KA, P, QTR, R, STEP, V, W) C C LATEST REVISION - 03/15/90 (JRD) C C C *** COMPUTE LEVENBERG-MARQUARDT STEP USING MORE-HEBDEN TECHNIQUE ** C *** NL2SOL VERSION 2.2. *** C C C VARIABLE DECLARATIONS C C SCALAR ARGUMENTS INTEGER + IERR,KA,P C C ARRAY ARGUMENTS REAL + D(P),G(P),QTR(P),R(1),STEP(P),V(21),W(1) INTEGER + IPIVOT(P) C C LOCAL SCALARS REAL + A,ADI,ALPHAK,B,D1,D2,DFAC,DFACSQ,DST,DTOL,EIGHT,HALF,LK, + NEGONE,OLDPHI,ONE,P001,PHI,PHIMAX,PHIMIN,PSIFAC,RAD,SI,SJ, + SQRTAK,T,THREE,TTOL,TWOPSI,UK,WL,ZERO INTEGER + DGNORM,DST0,DSTNRM,DSTSAV,EPSLON,GTSTEP,I,I1,IP1,J1,K, + KALIM,L,LK0,NREDUC,PHIPIN,PHMNFC,PHMXFC,PP1O2,PREDUC,RAD0, + RADIUS,RES,RES0,RMAT,RMAT0,STPPAR,UK0 C C EXTERNAL FUNCTIONS REAL + DOTPRD,V2NORM EXTERNAL DOTPRD,V2NORM C C EXTERNAL SUBROUTINES EXTERNAL LITVMU,LIVMUL,VCOPY C C INTRINSIC FUNCTIONS INTRINSIC ABS,MAX,MIN,SQRT C C *** PARAMETER DECLARATIONS *** C C INTEGER IERR, KA, P C INTEGER IPIVOT(P) C REAL D(P), G(P), QTR(P), R(1), STEP(P), V(21), W(1) C DIMENSION 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 QRFGS -- 0 MEANS R HAS C FULL RANK. C IPIVOT (I/O) = PERMUTATION ARRAY FROM QRFACT OR QRFGS, 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 LMSTEP 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 DOTPRD - RETURNS INNER PRODUCT OF TWO VECTORS. C LITVMU - APPLY INVERSE-TRANSPOSE OF COMPACT LOWER TRIANG. MATRIX. C LIVMUL - APPLY INVERSE OF COMPACT LOWER TRIANG. MATRIX. C VCOPY - COPIES ONE VECTOR TO ANOTHER. C V2NORM - RETURNS 2-NORM OF A VECTOR. C C *** REFERENCES *** C C 1. DENNIS, J.E., GAY, D.M., AND WELSCH, R.E. (1980), AN ADAPTIVE C NONLINEAR LEAST-SQUARES ALGORITHM, (SUBMITTED TO ACM C TRANS. MATH. SOFTWARE). C 2. GAY, D.M. (1979), COMPUTING OPTIMAL ELLIPTICALLY CONSTRAINED C STEPS, MRC TECH. SUMMARY REPORT NO. 2013, MATH RESEARCH C CENTER, UNIV. OF WISCONSIN-MADISON. 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 C INTEGER DSTSAV, I, IP1, I1, J1, K, KALIM, L, LK0, PHIPIN, C 1 PP1O2, RES, RES0, RMAT, RMAT0, UK0 C REAL A, ADI, ALPHAK, B, DFACSQ, DST, DTOL, D1, D2, C 1 LK, OLDPHI, PHI, PHIMAX, PHIMIN, PSIFAC, RAD, C 2 SI, SJ, SQRTAK, T, TWOPSI, UK, WL C C *** CONSTANTS *** C REAL DFAC, EIGHT, HALF, NEGONE, ONE, P001, THREE, C 1 TTOL, ZERO C C/ C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C C EXTERNAL DOTPRD, LITVMU, LIVMUL, VCOPY, V2NORM C REAL DOTPRD, V2NORM C C *** SUBSCRIPTS FOR V *** C C INTEGER DGNORM, DSTNRM, DST0, EPSLON, GTSTEP, NREDUC, PHMNFC, C 1 PHMXFC, PREDUC, RADIUS, RAD0, STPPAR DATA DGNORM/1/, DSTNRM/2/, DST0/3/, EPSLON/19/, + GTSTEP/4/, NREDUC/6/, PHMNFC/20/, + PHMXFC/21/, PREDUC/7/, RADIUS/8/, + RAD0/9/, STPPAR/5/ C DATA DFAC/256.0E0/, EIGHT/8.0E0/, HALF/0.5E0/, NEGONE/-1.0E0/, + ONE/1.0E0/, P001/1.0E-3/, THREE/3.0E0/, TTOL/2.5E0/, + ZERO/0.0E0/ 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. ALPHAK = 0.0E0 PSIFAC = 0.0E0 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) + PSIFAC = V(EPSLON)/((EIGHT*(V(PHMNFC) + ONE) + THREE) * RAD**2) 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*DOTPRD(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 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 LITVMU AND LIVMUL. CALL LITVMU(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 = V2NORM(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 LIVMUL(P, STEP, R, STEP) T = ONE / V2NORM(P, STEP) W(PHIPIN) = (T/DST)*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) = V2NORM(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. ALPHAK = ABS(V(STPPAR)) * V(RAD0)/RAD C C C *** TOP OF LOOP -- INCREMENT KA, COPY R TO RMAT, QTR TO RES *** C 110 KA = KA + 1 CALL VCOPY(PP1O2, W(RMAT), R) CALL VCOPY(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) + ALPHAK = UK * MAX(P001, SQRT(LK/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 LITVMU(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 = V2NORM(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 - DOTPRD(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 LIVMUL(P, STEP, W(RMAT), STEP) DO 340 I = 1, P 340 STEP(I) = STEP(I) / SQRT(W(I)) T = ONE / V2NORM(P, STEP) ALPHAK = ALPHAK + T*PHI*T/RAD LK = MAX(LK, ALPHAK) 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) = DOTPRD(P, STEP, G) 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 RETURN C C *** LAST CARD OF LMSTEP FOLLOWS *** END *LSQRT SUBROUTINE LSQRT(N1, N, L, A, IRC) C C LATEST REVISION - 03/15/90 (JRD) 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 C VARIABLE DECLARATIONS C C SCALAR ARGUMENTS INTEGER + IRC,N,N1 C C ARRAY ARGUMENTS REAL + A(1),L(1) C C LOCAL SCALARS REAL + T,TD,ZERO INTEGER + I,I0,IJ,IK,IM1,J,J0,JK,JM1,K C C INTRINSIC FUNCTIONS INTRINSIC SQRT C C *** PARAMETERS *** C C INTEGER N1, N, IRC C REAL L(1), A(1) C DIMENSION L(N*(N+1)/2), A(N*(N+1)/2) C C *** LOCAL VARIABLES *** C C INTEGER I, IJ, IK, IM1, I0, J, JK, JM1, J0, K C REAL T, TD, ZERO C C/ DATA ZERO/0.0E0/ 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) 50 CONTINUE C IRC = 0 GO TO 999 C 60 L(I0) = T IRC = I C 999 RETURN C C *** LAST CARD OF LSQRT *** END *LSVMIN REAL FUNCTION LSVMIN(P, L, X, Y) C C LATEST REVISION - 03/15/90 (JRD) C C *** ESTIMATE SMALLEST SING. VALUE OF PACKED LOWER TRIANG. MATRIX L C C C VARIABLE DECLARATIONS C C SCALAR ARGUMENTS INTEGER + P C C ARRAY ARGUMENTS REAL + L(1),X(P),Y(P) C C LOCAL SCALARS REAL + B,HALF,ONE,PSJ,R9973,SMINUS,SPLUS,T,XMINUS,XPLUS,ZERO INTEGER + I,II,IX,J,J0,JI,JJ,JJJ,JM1,PPLUS1 C C EXTERNAL FUNCTIONS REAL + V2NORM EXTERNAL V2NORM C C INTRINSIC FUNCTIONS INTRINSIC ABS,MOD C C *** PARAMETER DECLARATIONS *** C C INTEGER P C REAL L(1), 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 LSVMIN 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 LSVMIN RETURNS ZERO, THEN SOME COMPONENTS OF X C ARE ZERO AND THE REST RETAIN THEIR INPUT VALUES. C Y (OUT) IF LSVMIN 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 LSVMIN 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 LSVMIN RETURNS). C C *** APPLICATION AND USAGE RESTRICTIONS *** C C THERE ARE NO USAGE RESTRICTIONS. C C *** ALGORITHM NOTES *** C C THE ALGORITHM IS BASED ON (1), WITH THE ADDITIONAL PROVISION THAT C LSVMIN = 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 V2NORM - 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 C INTEGER I, II, IX, J, JI, JJ, JJJ, JM1, J0, PPLUS1 C REAL B, PSJ, SMINUS, SPLUS, T, XMINUS, XPLUS C C *** CONSTANTS *** C C REAL HALF, ONE, R9973, ZERO C C/ C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C C EXTERNAL V2NORM C REAL V2NORM C DATA IX/2/ DATA HALF/0.5E0/, ONE/1.0E0/, R9973/9973.0E0/, ZERO/0.0E0/ C C *** BODY *** C C *** FIRST CHECK WHETHER TO RETURN LSVMIN = 0 AND INITIALIZE X *** C II = 0 DO 10 I = 1, P X(I) = ZERO II = II + I IF (L(II) .EQ. ZERO) GO TO 300 10 CONTINUE IF (MOD(IX, 9973) .EQ. 0) IX = 2 PPLUS1 = P + 1 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 TO 1 BY -1... DO 100 JJJ = 1, P J = PPLUS1 - 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 + 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 .EQ. 0) GO TO 100 DO 40 I = 1, JM1 JI = J0 + I X(I) = X(I) + L(JI)*XPLUS 40 CONTINUE 100 CONTINUE C C *** NORMALIZE X *** C T = ONE/V2NORM(P, X) DO 110 I = 1, P 110 X(I) = T*X(I) C C *** SOLVE L*Y = X AND RETURN SVMIN = 1/TWONORM(Y) *** C DO 200 J = 1, P PSJ = ZERO JM1 = J - 1 J0 = J*JM1/2 IF (JM1 .EQ. 0) GO TO 130 DO 120 I = 1, JM1 JI = J0 + I PSJ = PSJ + L(JI)*Y(I) 120 CONTINUE 130 JJ = J0 + J Y(J) = (X(J) - PSJ)/L(JJ) 200 CONTINUE C LSVMIN = ONE/V2NORM(P, Y) GO TO 999 C 300 LSVMIN = ZERO 999 RETURN C *** LAST CARD OF LSVMIN FOLLOWS *** END *LTSQAR SUBROUTINE LTSQAR(N, A, L) C C LATEST REVISION - 03/15/90 (JRD) C C C *** SET A TO LOWER TRIANGLE OF (L**T) * L *** C C *** L = N X N LOWER TRIANG. MATRIX STORED ROWWISE. *** C *** A IS ALSO STORED ROWWISE AND MAY SHARE STORAGE WITH L. *** C C C VARIABLE DECLARATIONS C C SCALAR ARGUMENTS INTEGER + N C C ARRAY ARGUMENTS REAL + A(*),L(*) C C LOCAL SCALARS REAL + LII,LJ INTEGER + I,I1,II,IIM1,J,K,M C C INTEGER N C REAL A(1), L(1) C DIMENSION A(N*(N+1)/2), L(N*(N+1)/2) C C INTEGER I, II, IIM1, I1, J, K, M C REAL LII, LJ C II = 0 DO 50 I = 1, N I1 = II + 1 II = II + I M = 1 IF (I .EQ. 1) GO TO 30 IIM1 = II - 1 DO 20 J = I1, IIM1 LJ = L(J) DO 10 K = I1, J A(M) = A(M) + LJ*L(K) M = M + 1 10 CONTINUE 20 CONTINUE 30 LII = L(II) DO 40 J = I1, II 40 A(J) = LII * L(J) 50 CONTINUE C RETURN C *** LAST CARD OF LTSQAR FOLLOWS *** END *MADJ SUBROUTINE MADJ(N, P, X, NF, J, UIPARM, URPARM, UFPARM) C C VARIABLE DECLARATIONS C C SCALAR ARGUMENTS INTEGER + N,NF,P C C ARRAY ARGUMENTS REAL + J(N,P),URPARM(*),X(P) INTEGER + UIPARM(*) C C SUBROUTINE ARGUMENTS EXTERNAL UFPARM C C INTRINSIC FUNCTIONS INTRINSIC COS,SIN C J(1,1) = 2.0E0*X(1) + X(2) J(1,2) = 2.0E0*X(2) + X(1) J(2,1) = COS(X(1)) J(2,2) = 0.0E0 J(3,1) = 0.0E0 J(3,2) = -SIN(X(2)) RETURN END *MADR SUBROUTINE MADR(N, P, X, NF, R, UIPARM, URPARM, UFPARM) C C VARIABLE DECLARATIONS C C SCALAR ARGUMENTS INTEGER + N,NF,P C C ARRAY ARGUMENTS REAL + R(N),URPARM(*),X(P) INTEGER + UIPARM(*) C C SUBROUTINE ARGUMENTS EXTERNAL UFPARM C C INTRINSIC FUNCTIONS INTRINSIC COS,SIN C R(1) = X(1)**2 + X(2)**2 + X(1)*X(2) R(2) = SIN(X(1)) R(3) = COS(X(2)) RETURN END *NL2ITR SUBROUTINE NL2ITR (D, IV, J, N, NN, P, R, V, X) C C LATEST REVISION - 03/15/90 (JRD) C C C *** CARRY OUT NL2SOL (NONLINEAR LEAST-SQUARES) ITERATIONS *** C *** (NL2SOL VERSION 2.2) *** C C C VARIABLE DECLARATIONS C C SCALAR ARGUMENTS INTEGER + N,NN,P C C ARRAY ARGUMENTS REAL + D(P),J(NN,P),R(N),V(*),X(P) INTEGER + IV(*) C C LOCAL SCALARS REAL + E,HALF,NEGONE,ONE,RDOF1,STTSST,T,T1,ZERO INTEGER + CNVCOD,COSMIN,COVMAT,COVPRT,COVREQ,D0INIT,DGNORM,DIG,DIG1, + DINIT,DSTNRM,DTYPE,DUMMY,F,F0,FDIF,FUZZ,G,G01,G1,GTSTEP,H, + H0,H1,I,IERR,IM1,INCFAC,INITS,IPIV0,IPIV1,IPIVI,IPIVK, + IPIVOT,IPK,IRC,JTINIT,JTOL1,K,KAGQT,KALM,KM1,L,LKY,LKY1, + LMAT,LMAT1,LMAX0,LSTGST,M,MODE,MODEL,MXFCAL,MXITER,NFCALL, + NFCOV,NFGCAL,NGCALL,NGCOV,NITER,NVSAVE,PHMXFC,PP1O2, + PREDUC,QTR,QTR1,RAD0,RADFAC,RADINC,RADIUS,RD,RD0,RD1,RDK, + RESTOR,RLIMIT,RSAVE,RSAVE1,S,S1,SIZE,SMH,SSTEP,STEP,STEP1, + STGLIM,STLSTG,STPMOD,STPPAR,SUSED,SWITCH,TEMP1,TEMP2, + TOOBIG,TUNER4,TUNER5,VSAVE1,W,W1,WSCALE,X0,X01,XIRC C C EXTERNAL FUNCTIONS REAL + DOTPRD,R1MACH,V2NORM LOGICAL + STOPX EXTERNAL DOTPRD,R1MACH,V2NORM,STOPX C C EXTERNAL SUBROUTINES EXTERNAL ASSESS,COVCLC,DUPDAT,GQTSTP,ITSMRY,LMSTEP,PARCHK,QAPPLY, + QRFACT,RPTMUL,SLUPDT,SLVMUL,VAXPY,VCOPY,VSCOPY C C INTRINSIC FUNCTIONS INTRINSIC ABS,SQRT C C *** PARAMETER DECLARATIONS *** C C INTEGER IV(1), N, NN, P C REAL D(P), J(NN,P), R(N), V(1), X(P) C DIMENSION IV(60+P), V(93 + 2*N + P*(3*P+31)/2) C C C-------------------------- PARAMETER USAGE -------------------------- C C D.... SCALE VECTOR. C IV... INTEGER VALUE ARRAY. C J.... N BY P JACOBIAN MATRIX (LEAD DIMENSION NN). C N.... NUMBER OF OBSERVATIONS (COMPONENTS IN R). C NN... LEAD DIMENSION OF J. C P.... NUMBER OF PARAMETERS (COMPONENTS IN X). C R.... RESIDUAL VECTOR. C V.... FLOATING-POINT VALUE ARRAY. C X.... PARAMETER VECTOR. C C *** DISCUSSION *** C C PARAMETERS IV, N, P, V, AND X ARE THE SAME AS THE CORRESPOND- C ING ONES TO NL2SOL (WHICH SEE), EXCEPT THAT V CAN BE SHORTER C (SINCE THE PART OF V THAT NL2SOL USES FOR STORING D, J, AND R IS C NOT NEEDED). MOREOVER, COMPARED WITH NL2SOL, IV(1) MAY HAVE THE C TWO ADDITIONAL OUTPUT VALUES 1 AND 2, WHICH ARE EXPLAINED BELOW, C AS IS THE USE OF IV(TOOBIG) AND IV(NFGCAL). THE VALUES IV(D), C IV(J), AND IV(R), WHICH ARE OUTPUT VALUES FROM NL2SOL (AND C NL2SNO), ARE NOT REFERENCED BY NL2ITR OR THE SUBROUTINES IT CALLS. C ON A FRESH START, I.E., A CALL ON NL2ITR WITH IV(1) = 0 OR 12, C NL2ITR ASSUMES THAT R = R(X), THE RESIDUAL AT X, AND J = J(X), C THE CORRESPONDING JACOBIAN MATRIX OF R AT X. C C IV(1) = 1 MEANS THE CALLER SHOULD SET R TO R(X), THE RESIDUAL AT X, C AND CALL NL2ITR AGAIN, HAVING CHANGED NONE OF THE OTHER C PARAMETERS. AN EXCEPTION OCCURS IF R CANNOT BE EVALUATED C AT X (E.G. IF R WOULD OVERFLOW), WHICH MAY HAPPEN BECAUSE C OF AN OVERSIZED STEP. IN THIS CASE THE CALLER SHOULD SET C IV(TOOBIG) = IV(2) TO 1, WHICH WILL CAUSE NL2ITR TO IG- C NORE R AND TRY A SMALLER STEP. THE PARAMETER NF THAT C NL2SOL PASSES TO CALCR (FOR POSSIBLE USE BY CALCJ) IS A C COPY OF IV(NFCALL) = IV(6). C IV(1) = 2 MEANS THE CALLER SHOULD SET J TO J(X), THE JACOBIAN MATRIX C OF R AT X, AND CALL NL2ITR AGAIN. THE CALLER MAY CHANGE C D AT THIS TIME, BUT SHOULD NOT CHANGE ANY OF THE OTHER C PARAMETERS. THE PARAMETER NF THAT NL2SOL PASSES TO C CALCJ IS IV(NFGCAL) = IV(7). IF J CANNOT BE EVALUATED C AT X, THEN THE CALLER MAY SET IV(NFGCAL) TO 0, IN WHICH C CASE NL2ITR WILL RETURN WITH IV(1) = 15. 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 C MCS-7600324, DCR75-10143, 76-14311DSS, MCS76-11989, AND C MCS-7906671. C (SEE NL2SOL FOR REFERENCES.) C C+++++++++++++++++++++++++++ DECLARATIONS ++++++++++++++++++++++++++++ C C *** LOCAL VARIABLES *** C C INTEGER DUMMY, DIG1, G1, G01, H0, H1, I, IM1, IPIVI, IPIVK, IPIV1, C 1 IPK, K, KM1, L, LKY1, LMAT1, LSTGST, M, PP1O2, QTR1, C 2 RDK, RD0, RD1, RSAVE1, SMH, SSTEP, STEP1, STPMOD, S1, C 3 TEMP1, TEMP2, W1, X01 C REAL E, RDOF1, STTSST, T, T1 C C *** CONSTANTS *** C C REAL HALF, NEGONE, ONE, ZERO C C/ C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C C EXTERNAL ASSESS, COVCLC, DOTPRD, DUPDAT, GQTSTP, ITSMRY, LMSTEP, C 1 PARCHK, QAPPLY, QRFACT, RPTMUL, SLUPDT, SLVMUL, STOPX, C 2 VAXPY, VCOPY, VSCOPY, V2NORM C LOGICAL STOPX C REAL DOTPRD, R1MACH, V2NORM C C ASSESS... ASSESSES CANDIDATE STEP. C COVCLC... COMPUTES COVARIANCE MATRIX. C DOTPRD... RETURNS INNER PRODUCT OF TWO VECTORS. C DUPDAT... UPDATES SCALE VECTOR D. C GQTSTP... COMPUTES GOLDFELD-QUANDT-TROTTER STEP (AUGMENTED MODEL). C ITSMRY... PRINTS ITERATION SUMMARY AND INFO ABOUT INITIAL AND FINAL X. C LMSTEP... COMPUTES LEVENBERG-MARQUARDT STEP (GAUSS-NEWTON MODEL). C PARCHK... CHECKS VALIDITY OF INPUT IV AND V VALUES. C QAPPLY... APPLIES ORTHOGONAL MATRIX Q FROM QRFACT TO A VECTOR. C QRFACT... COMPUTES QR DECOMPOSITION OF A MATRIX VIA HOUSEHOLDER TRANS. C RPTMUL... MULTIPLIES VECTOR BY THE R MATRIX (AND/OR ITS TRANSPOSE) C STORED BY QRFACT. C SLUPDT... PERFORMS QUASI-NEWTON UPDATE ON COMPACTLY STORED LOWER TRI- C ANGLE OF A SYMMETRIC MATRIX. C STOPX.... RETURNS .TRUE. IF THE BREAK KEY HAS BEEN PRESSED. C VAXPY.... COMPUTES SCALAR TIMES ONE VECTOR PLUS ANOTHER. C VCOPY.... COPIES ONE VECTOR TO ANOTHER. C VSCOPY... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR. C V2NORM... RETURNS THE 2-NORM OF A VECTOR. C C *** SUBSCRIPTS FOR IV AND V *** C C INTEGER CNVCOD, COSMIN, COVMAT, COVPRT, COVREQ, DGNORM, DIG, C 1 DINIT, DSTNRM, DTYPE, D0INIT, F, FDIF, FUZZ, C 2 F0, G, GTSTEP, H, IERR, INCFAC, INITS, IPIVOT, IPIV0, IRC, C 3 JTINIT, JTOL1, KAGQT, KALM, LKY, LMAT, LMAX0, MODE, MODEL, C 4 MXFCAL, MXITER, NFCALL, NFGCAL, NFCOV, NGCOV, NGCALL, C 5 NITER, NVSAVE, PHMXFC, PREDUC, QTR, RADFAC, RADINC, C 6 RADIUS, RAD0, RD, RESTOR, RLIMIT, RSAVE, S, SIZE, STEP, C 7 STGLIM, STLSTG, STPPAR, SUSED, SWITCH, TOOBIG, TUNER4, C 8 TUNER5, VSAVE1, W, WSCALE, XIRC, X0 C C *** IV SUBSCRIPT VALUES *** C DATA CNVCOD/34/, COVMAT/26/, COVPRT/14/, + COVREQ/15/, DIG/43/, DTYPE/16/, G/28/, H/44/, + IERR/32/, INITS/25/, IPIVOT/61/, IPIV0/60/, + IRC/3/, KAGQT/35/, KALM/36/, LKY/37/, LMAT/58/, + MODE/38/, MODEL/5/, MXFCAL/17/, MXITER/18/, + NFCALL/6/, NFGCAL/7/, NFCOV/40/, NGCOV/41/, + NGCALL/30/, NITER/31/, QTR/49/, + RADINC/8/, RD/51/, RESTOR/9/, RSAVE/52/, S/53/, + STEP/55/, STGLIM/11/, STLSTG/56/, SUSED/57/, + SWITCH/12/, TOOBIG/2/, W/59/, XIRC/13/, X0/60/ C C *** V SUBSCRIPT VALUES *** C DATA COSMIN/43/, DGNORM/1/, DINIT/38/, DSTNRM/2/, + D0INIT/37/, F/10/, FDIF/11/, FUZZ/45/, + F0/13/, GTSTEP/4/, INCFAC/23/, + JTINIT/39/, JTOL1/87/, LMAX0/35/, + NVSAVE/9/, PHMXFC/21/, PREDUC/7/, + RADFAC/16/, RADIUS/8/, RAD0/9/, RLIMIT/42/, + SIZE/47/, STPPAR/5/, TUNER4/29/, TUNER5/30/, + VSAVE1/78/, WSCALE/48/ C C DATA HALF/0.5E0/, NEGONE/-1.0E0/, ONE/1.0E0/, ZERO/0.0E0/ C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C I = IV(1) IF (I .EQ. 1) GO TO 20 IF (I .EQ. 2) GO TO 50 C C *** CHECK VALIDITY OF IV AND V INPUT VALUES *** C C *** NOTE -- IF IV(1) = 0, THEN PARCHK CALLS DFAULT(IV, V) *** CALL PARCHK(IV, N, NN, P, V) I = IV(1) - 2 IF (I .GT. 10) GO TO 999 GO TO (350, 350, 350, 350, 350, 350, 195, 160, 195, 10), I C C *** INITIALIZATION AND STORAGE ALLOCATION *** C 10 IV(NITER) = 0 IV(NFCALL) = 1 IV(NGCALL) = 1 IV(NFGCAL) = 1 IV(MODE) = -1 IV(STGLIM) = 2 IV(TOOBIG) = 0 IV(CNVCOD) = 0 IV(COVMAT) = 0 IV(NFCOV) = 0 IV(NGCOV) = 0 IV(KALM) = -1 IV(RADINC) = 0 IV(S) = JTOL1 + 2*P PP1O2 = P * (P + 1) / 2 IV(X0) = IV(S) + PP1O2 IV(STEP) = IV(X0) + P IV(STLSTG) = IV(STEP) + P IV(DIG) = IV(STLSTG) + P IV(G) = IV(DIG) + P IV(LKY) = IV(G) + P IV(RD) = IV(LKY) + P IV(RSAVE) = IV(RD) + P IV(QTR) = IV(RSAVE) + N IV(H) = IV(QTR) + N IV(W) = IV(H) + PP1O2 IV(LMAT) = IV(W) + 4*P + 7 C +++ LENGTH OF W = P*(P+9)/2 + 7. LMAT IS CONTAINED IN W. IF (V(DINIT) .GE. ZERO) CALL VSCOPY(P, D, V(DINIT)) IF (V(JTINIT) .GT. ZERO) CALL VSCOPY(P, V(JTOL1), V(JTINIT)) I = JTOL1 + P IF (V(D0INIT) .GT. ZERO) CALL VSCOPY(P, V(I), V(D0INIT)) V(RAD0) = ZERO V(STPPAR) = ZERO V(RADIUS) = V(LMAX0) / (ONE + V(PHMXFC)) C C *** SET INITIAL MODEL AND S MATRIX *** C IV(MODEL) = 1 IF (IV(INITS) .EQ. 2) IV(MODEL) = 2 S1 = IV(S) IF (IV(INITS) .EQ. 0) CALL VSCOPY(PP1O2, V(S1), ZERO) C C *** COMPUTE FUNCTION VALUE (HALF THE SUM OF SQUARES) *** C 20 T = V2NORM(N, R) IF (T .GT. V(RLIMIT)) IV(TOOBIG) = 1 IF (IV(TOOBIG) .NE. 0) GO TO 30 V(F) = 0.0 IF (T.GT.SQRT(R1MACH(1))) V(F) = HALF * T**2 30 IF (IV(MODE)) 40, 350, 730 C 40 IF (IV(TOOBIG) .EQ. 0) GO TO 60 IV(1) = 13 GO TO 900 C C *** MAKE SURE JACOBIAN COULD BE COMPUTED *** C 50 IF (IV(NFGCAL) .NE. 0) GO TO 60 IV(1) = 15 GO TO 900 C C *** COMPUTE GRADIENT *** C 60 IV(KALM) = -1 G1 = IV(G) DO 70 I = 1, P V(G1) = DOTPRD(N, R, J(1,I)) G1 = G1 + 1 70 CONTINUE IF (IV(MODE) .GT. 0) GO TO 710 C C *** UPDATE D AND MAKE COPIES OF R FOR POSSIBLE USE LATER *** C IF (IV(DTYPE) .GT. 0) CALL DUPDAT(D, IV, J, N, NN, P, V) RSAVE1 = IV(RSAVE) CALL VCOPY(N, V(RSAVE1), R) QTR1 = IV(QTR) CALL VCOPY(N, V(QTR1), R) C C *** COMPUTE D**-1 * GRADIENT *** C G1 = IV(G) DIG1 = IV(DIG) K = DIG1 DO 80 I = 1, P V(K) = V(G1) / D(I) K = K + 1 G1 = G1 + 1 80 CONTINUE V(DGNORM) = V2NORM(P, V(DIG1)) C IF (IV(CNVCOD) .NE. 0) GO TO 700 IF (IV(MODE) .EQ. 0) GO TO 570 IV(MODE) = 0 C C C----------------------------- MAIN LOOP ----------------------------- C C C *** PRINT ITERATION SUMMARY, CHECK ITERATION LIMIT *** C 150 CALL ITSMRY(D, IV, P, V, X) 160 K = IV(NITER) IF (K .LT. IV(MXITER)) GO TO 170 IV(1) = 10 GO TO 900 170 IV(NITER) = K + 1 C C *** UPDATE RADIUS *** C IF (K .EQ. 0) GO TO 185 STEP1 = IV(STEP) DO 180 I = 1, P V(STEP1) = D(I) * V(STEP1) STEP1 = STEP1 + 1 180 CONTINUE STEP1 = IV(STEP) V(RADIUS) = V(RADFAC) * V2NORM(P, V(STEP1)) C C *** INITIALIZE FOR START OF NEXT ITERATION *** C 185 X01 = IV(X0) V(F0) = V(F) IV(KAGQT) = -1 IV(IRC) = 4 IV(H) = -ABS(IV(H)) IV(SUSED) = IV(MODEL) C C *** COPY X TO X0 *** C CALL VCOPY(P, V(X01), X) C C *** CHECK STOPX AND FUNCTION EVALUATION LIMIT *** C 190 IF (.NOT. STOPX(DUMMY)) GO TO 200 IV(1) = 11 GO TO 205 C C *** COME HERE WHEN RESTARTING AFTER FUNC. EVAL. LIMIT OR STOPX. C 195 IF (V(F) .GE. V(F0)) GO TO 200 V(RADFAC) = ONE K = IV(NITER) GO TO 170 C 200 IF (IV(NFCALL) .LT. IV(MXFCAL) + IV(NFCOV)) GO TO 210 IV(1) = 9 205 IF (V(F) .GE. V(F0)) GO TO 900 C C *** IN CASE OF STOPX OR FUNCTION EVALUATION LIMIT WITH C *** IMPROVED V(F), EVALUATE THE GRADIENT AT X. C IV(CNVCOD) = IV(1) GO TO 560 C C. . . . . . . . . . . . . COMPUTE CANDIDATE STEP . . . . . . . . . . C 210 STEP1 = IV(STEP) W1 = IV(W) IF (IV(MODEL) .EQ. 2) GO TO 240 C C *** COMPUTE LEVENBERG-MARQUARDT STEP *** C QTR1 = IV(QTR) IF (IV(KALM) .GE. 0) GO TO 215 RD1 = IV(RD) IF (-1 .EQ. IV(KALM)) CALL QRFACT(NN, N, P, J, V(RD1), + IV(IPIVOT), IV(IERR), 0, V(W1)) CALL QAPPLY(NN, N, P, J, V(QTR1), IV(IERR)) 215 H1 = IV(H) IF (H1 .GT. 0) GO TO 230 C C *** COPY R MATRIX TO H *** C H1 = -H1 IV(H) = H1 K = H1 RD1 = IV(RD) V(K) = V(RD1) IF (P .EQ. 1) GO TO 230 DO 220 I = 2, P CALL VCOPY(I-1, V(K+1), J(1,I)) K = K + I RD1 = RD1 + 1 V(K) = V(RD1) 220 CONTINUE C 230 G1 = IV(G) CALL LMSTEP(D, V(G1), IV(IERR), IV(IPIVOT), IV(KALM), P, + V(QTR1), V(H1), V(STEP1), V, V(W1)) GO TO 310 C C *** COMPUTE GOLDFELD-QUANDT-TROTTER STEP (AUGMENTED MODEL) *** C 240 IF (IV(H) .GT. 0) GO TO 300 C C *** SET H TO D**-1 * ( (J**T)*J + S) ) * D**-1. *** C H1 = -IV(H) IV(H) = H1 S1 = IV(S) IF (IV(KALM) .GE. 0) GO TO 270 C C *** J IS IN ITS ORIGINAL FORM *** C DO 260 I = 1, P T = ONE / D(I) DO 250 K = 1, I V(H1) = T*(DOTPRD(N,J(1,I),J(1,K))+V(S1)) / D(K) H1 = H1 + 1 S1 = S1 + 1 250 CONTINUE 260 CONTINUE GO TO 300 C C *** LMSTEP HAS APPLIED QRFACT TO J *** C 270 SMH = S1 - H1 H0 = H1 - 1 IPIV1 = IV(IPIVOT) T1 = ONE / D(IPIV1) RD0 = IV(RD) - 1 RDOF1 = V(RD0 + 1) DO 290 I = 1, P L = IPIV0 + I IPIVI = IV(L) H1 = H0 + IPIVI*(IPIVI-1)/2 L = H1 + IPIVI M = L + SMH C *** V(L) = H(IPIVOT(I), IPIVOT(I)) *** C *** V(M) = S(IPIVOT(I), IPIVOT(I)) *** T = ONE / D(IPIVI) RDK = RD0 + I E = V(RDK)**2 IF (I .GT. 1) E = E + DOTPRD(I-1, J(1,I), J(1,I)) V(L) = (E + V(M)) * T**2 IF (I .EQ. 1) GO TO 290 L = H1 + IPIV1 IF (IPIVI .LT. IPIV1) L = L + + ((IPIV1-IPIVI)*(IPIV1+IPIVI-3))/2 M = L + SMH C *** V(L) = H(IPIVOT(I), IPIVOT(1)) *** C *** V(M) = S(IPIVOT(I), IPIVOT(1)) *** V(L) = T * (RDOF1 * J(1,I) + V(M)) * T1 IF (I .EQ. 2) GO TO 290 IM1 = I - 1 DO 280 K = 2, IM1 IPK = IPIV0 + K IPIVK = IV(IPK) L = H1 + IPIVK IF (IPIVI .LT. IPIVK) L = L + + ((IPIVK-IPIVI)*(IPIVK+IPIVI-3))/2 M = L + SMH C *** V(L) = H(IPIVOT(I), IPIVOT(K)) *** C *** V(M) = S(IPIVOT(I), IPIVOT(K)) *** KM1 = K - 1 RDK = RD0 + K V(L) = T * (DOTPRD(KM1, J(1,I), J(1,K)) + + V(RDK)*J(K,I) + V(M)) / D(IPIVK) 280 CONTINUE 290 CONTINUE C C *** COMPUTE ACTUAL GOLDFELD-QUANDT-TROTTER STEP *** C 300 H1 = IV(H) DIG1 = IV(DIG) LMAT1 = IV(LMAT) CALL GQTSTP(D, V(DIG1), V(H1), IV(KAGQT), V(LMAT1), P, V(STEP1), + V, V(W1)) C C C *** COMPUTE R(X0 + STEP) *** C 310 IF (IV(IRC) .EQ. 6) GO TO 350 X01 = IV(X0) STEP1 = IV(STEP) CALL VAXPY(P, X, ONE, V(STEP1), V(X01)) IV(NFCALL) = IV(NFCALL) + 1 IV(1) = 1 IV(TOOBIG) = 0 GO TO 999 C C. . . . . . . . . . . . . ASSESS CANDIDATE STEP . . . . . . . . . . . C 350 STEP1 = IV(STEP) LSTGST = IV(STLSTG) X01 = IV(X0) CALL ASSESS(D, IV, P, V(STEP1), V(LSTGST), V, X, V(X01)) C C *** IF NECESSARY, SWITCH MODELS AND/OR RESTORE R *** C IF (IV(SWITCH) .EQ. 0) GO TO 360 IV(H) = -ABS(IV(H)) IV(SUSED) = IV(SUSED) + 2 CALL VCOPY(NVSAVE, V, V(VSAVE1)) 360 IF (IV(RESTOR) .EQ. 0) GO TO 390 RSAVE1 = IV(RSAVE) CALL VCOPY(N, R, V(RSAVE1)) 390 L = IV(IRC) - 4 STPMOD = IV(MODEL) IF (L .GT. 0) GO TO (410,440,450,450,450,450,450,450,640,570), L C C *** DECIDE WHETHER TO CHANGE MODELS *** C E = V(PREDUC) - V(FDIF) SSTEP = IV(LKY) S1 = IV(S) CALL SLVMUL(P, V(SSTEP), V(S1), V(STEP1)) STTSST = HALF * DOTPRD(P, V(STEP1), V(SSTEP)) IF (IV(MODEL) .EQ. 1) STTSST = -STTSST IF (ABS(E + STTSST) * V(FUZZ) .GE. ABS(E)) GO TO 400 C C *** SWITCH MODELS *** C IV(MODEL) = 3 - IV(MODEL) IF (IV(MODEL) .EQ. 1) IV(KAGQT) = -1 IF (IV(MODEL) .EQ. 2 .AND. IV(KALM) .GT. 0) IV(KALM) = 0 IF (-2 .LT. L) GO TO 480 IV(H) = -ABS(IV(H)) IV(SUSED) = IV(SUSED) + 2 CALL VCOPY(NVSAVE, V(VSAVE1), V) GO TO 420 C 400 IF (-3 .LT. L) GO TO 480 C C *** RECOMPUTE STEP WITH DECREASED RADIUS *** C V(RADIUS) = V(RADFAC) * V(DSTNRM) GO TO 190 C C *** RECOMPUTE STEP, SAVING V VALUES AND R IF NECESSARY *** C 410 V(RADIUS) = V(RADFAC) * V(DSTNRM) 420 IF (V(F) .GE. V(F0)) GO TO 190 RSAVE1 = IV(RSAVE) CALL VCOPY(N, V(RSAVE1), R) GO TO 190 C C *** COMPUTE STEP OF LENGTH V(LMAX0) FOR SINGULAR CONVERGENCE TEST C 440 V(RADIUS) = V(LMAX0) GO TO 210 C C *** CONVERGENCE OR FALSE CONVERGENCE *** C 450 IV(CNVCOD) = L IF (V(F) .GE. V(F0)) GO TO 700 IF (IV(XIRC) .EQ. 14) GO TO 700 IV(XIRC) = 14 C C. . . . . . . . . . . . PROCESS ACCEPTABLE STEP . . . . . . . . . . . C 480 IV(COVMAT) = 0 C C *** SET LKY = (J(X0)**T) * R(X) *** C LKY1 = IV(LKY) IF (IV(KALM) .GE. 0) GO TO 500 C C *** JACOBIAN HAS NOT BEEN MODIFIED *** C DO 490 I = 1, P V(LKY1) = DOTPRD(N, J(1,I), R) LKY1 = LKY1 + 1 490 CONTINUE GO TO 510 C C *** QRFACT HAS BEEN APPLIED TO J. STORE COPY OF R IN QTR AND *** C *** APPLY Q TO IT. *** C 500 QTR1 = IV(QTR) CALL VCOPY(N, V(QTR1), R) CALL QAPPLY(NN, N, P, J, V(QTR1), IV(IERR)) C C *** MULTIPLY TOP P-VECTOR IN QTR BY PERMUTED UPPER TRIANGLE *** C *** STORED BY QRFACT IN J AND RD. *** C RD1 = IV(RD) TEMP1 = IV(STLSTG) CALL RPTMUL(3, IV(IPIVOT), J, NN, P, V(RD1), V(QTR1), V(LKY1), + V(TEMP1)) C C *** SEE WHETHER TO SET V(RADFAC) BY GRADIENT TESTS *** C 510 IF (IV(IRC) .NE. 3) GO TO 560 STEP1 = IV(STEP) TEMP1 = IV(STLSTG) TEMP2 = IV(X0) C C *** SET TEMP1 = HESSIAN * STEP FOR USE IN GRADIENT TESTS *** C IF (STPMOD .EQ. 2) GO TO 530 C C *** STEP COMPUTED USING GAUSS-NEWTON MODEL *** C *** -- QRFACT HAS BEEN APPLIED TO J *** C RD1 = IV(RD) CALL RPTMUL(2, IV(IPIVOT), J, NN, P, V(RD1), + V(STEP1), V(TEMP1), V(TEMP2)) GO TO 560 C C *** STEP COMPUTED USING AUGMENTED MODEL *** C 530 H1 = IV(H) K = TEMP2 DO 540 I = 1, P V(K) = D(I) * V(STEP1) K = K + 1 STEP1 = STEP1 + 1 540 CONTINUE CALL SLVMUL(P, V(TEMP1), V(H1), V(TEMP2)) DO 550 I = 1, P V(TEMP1) = D(I) * V(TEMP1) TEMP1 = TEMP1 + 1 550 CONTINUE C C *** SAVE OLD GRADIENT AND COMPUTE NEW ONE *** C 560 IV(NGCALL) = IV(NGCALL) + 1 G1 = IV(G) G01 = IV(W) CALL VCOPY(P, V(G01), V(G1)) IV(1) = 2 GO TO 999 C C *** INITIALIZATIONS -- G0 = G - G0, ETC. *** C 570 G01 = IV(W) G1 = IV(G) CALL VAXPY(P, V(G01), NEGONE, V(G01), V(G1)) STEP1 = IV(STEP) TEMP1 = IV(STLSTG) TEMP2 = IV(X0) IF (IV(IRC) .NE. 3) GO TO 600 C C *** SET V(RADFAC) BY GRADIENT TESTS *** C C *** SET TEMP1 = D**-1 * (HESSIAN * STEP + (G(X0) - G(X))) *** C K = TEMP1 L = G01 DO 580 I = 1, P V(K) = (V(K) - V(L)) / D(I) K = K + 1 L = L + 1 580 CONTINUE C C *** DO GRADIENT TESTS *** C IF (V2NORM(P, V(TEMP1)) .LE. V(DGNORM) * V(TUNER4)) GO TO 590 IF (DOTPRD(P, V(G1), V(STEP1)) + .GE. V(GTSTEP) * V(TUNER5)) GO TO 600 590 V(RADFAC) = V(INCFAC) C C *** FINISH COMPUTING LKY = ((J(X) - J(X0))**T) * R *** C C *** CURRENTLY LKY = (J(X0)**T) * R *** C 600 LKY1 = IV(LKY) CALL VAXPY(P, V(LKY1), NEGONE, V(LKY1), V(G1)) C C *** DETERMINE SIZING FACTOR V(SIZE) *** C C *** SET TEMP1 = S * STEP *** S1 = IV(S) CALL SLVMUL(P, V(TEMP1), V(S1), V(STEP1)) C T1 = ABS(DOTPRD(P, V(STEP1), V(TEMP1))) T = ABS(DOTPRD(P, V(STEP1), V(LKY1))) V(SIZE) = ONE IF (T .LT. T1) V(SIZE) = T / T1 C C *** UPDATE S *** C CALL SLUPDT(V(S1), V(COSMIN), P, V(SIZE), V(STEP1), V(TEMP1), + V(TEMP2), V(G01), V(WSCALE), V(LKY1)) IV(1) = 2 GO TO 150 C C. . . . . . . . . . . . . . MISC. DETAILS . . . . . . . . . . . . . . C C *** BAD PARAMETERS TO ASSESS *** C 640 IV(1) = 14 GO TO 900 C C *** CONVERGENCE OBTAINED -- COMPUTE COVARIANCE MATRIX IF DESIRED *** C 700 IF (IV(COVREQ) .EQ. 0 .AND. IV(COVPRT) .EQ. 0) GO TO 760 IF (IV(COVMAT) .NE. 0) GO TO 760 IF (IV(CNVCOD) .GE. 7) GO TO 760 IV(MODE) = 0 710 CALL COVCLC(I, D, IV, J, N, NN, P, R, V, X) GO TO (720, 720, 740, 750), I 720 IV(NFCOV) = IV(NFCOV) + 1 IV(NFCALL) = IV(NFCALL) + 1 IV(RESTOR) = I IV(1) = 1 GO TO 999 C 730 IF (IV(RESTOR) .EQ. 1 .OR. IV(TOOBIG) .NE. 0) GO TO 710 IV(NFGCAL) = IV(NFCALL) 740 IV(NGCOV) = IV(NGCOV) + 1 IV(NGCALL) = IV(NGCALL) + 1 IV(1) = 2 GO TO 999 C 750 IV(MODE) = 0 IF (IV(NITER) .EQ. 0) IV(MODE) = -1 C 760 IV(1) = IV(CNVCOD) IV(CNVCOD) = 0 C C *** PRINT SUMMARY OF FINAL ITERATION AND OTHER REQUESTED ITEMS *** C 900 CALL ITSMRY(D, IV, P, V, X) C 999 RETURN C C *** LAST CARD OF NL2ITR FOLLOWS *** END *NL2SNO SUBROUTINE NL2SNO(N, P, X, CALCR, IV, V, UIPARM, URPARM, UFPARM) C C LATEST REVISION - 03/15/90 (JRD) C C *** LIKE NL2SOL, BUT WITHOUT CALCJ -- MINIMIZE NONLINEAR SUM OF *** C *** SQUARES USING FINITE-DIFFERENCE JACOBIAN APPROXIMATIONS *** C *** (NL2SOL VERSION 2.2) *** C C C VARIABLE DECLARATIONS C C SCALAR ARGUMENTS INTEGER + N,P C C ARRAY ARGUMENTS REAL + URPARM(1),V(1),X(P) INTEGER + IV(1),UIPARM(1) C C SUBROUTINE ARGUMENTS EXTERNAL CALCR,UFPARM C C LOCAL SCALARS REAL + H,HFAC,HLIM,NEGPT5,ONE,XK,ZERO INTEGER + COVPRT,COVREQ,D,D1,DK,DLTFDJ,DTYPE,I,J,J1,J1K,K,NF,NFCALL, + NFGCAL,R,R1,RN,TOOBIG LOGICAL + STRTED C C EXTERNAL FUNCTIONS REAL + RMDCON EXTERNAL RMDCON C C EXTERNAL SUBROUTINES EXTERNAL DFAULT,ITSMRY,NL2ITR,VSCOPY C C INTRINSIC FUNCTIONS INTRINSIC ABS,MAX C C INTEGER N, P, IV(1), UIPARM(1) C REAL X(P), V(1), URPARM(1) C DIMENSION IV(60+P), V(93 + N*P + 3*N + P*(3*P+33)/2) C EXTERNAL CALCR, UFPARM C C----------------------------- DISCUSSION ---------------------------- C C THE PARAMETERS FOR NL2SNO ARE THE SAME AS THOSE FOR NL2SOL C (WHICH SEE), EXCEPT THAT CALCJ IS OMITTED. INSTEAD OF CALLING C CALCJ TO OBTAIN THE JACOBIAN MATRIX OF R AT X, NL2SNO COMPUTES C AN APPROXIMATION TO IT BY FINITE (FORWARD) DIFFERENCES -- SEE C V(DLTFDJ) BELOW. NL2SNO USES FUNCTION VALUES ONLY WHEN COMPUT- C THE COVARIANCE MATRIX (RATHER THAN THE FUNCTIONS AND GRADIENTS C THAT NL2SOL MAY USE). TO DO SO, NL2SNO SETS IV(COVREQ) TO -1 IF C IV(COVPRT) = 1 WITH IV(COVREQ) = 0 AND TO MINUS ITS ABSOLUTE C VALUE OTHERWISE. THUS V(DELTA0) IS NEVER REFERENCED AND ONLY C V(DLTFDC) MATTERS -- SEE NL2SOL FOR A DESCRIPTION OF V(DLTFDC). C THE NUMBER OF EXTRA CALLS ON CALCR USED IN COMPUTING THE JACO- C BIAN APPROXIMATION ARE NOT INCLUDED IN THE FUNCTION EVALUATION C COUNT IV(NFCALL) AND ARE NOT OTHERWISE REPORTED. C C V(DLTFDJ)... V(36) HELPS CHOOSE THE STEP SIZE USED WHEN COMPUTING THE C FINITE-DIFFERENCE JACOBIAN MATRIX. FOR DIFFERENCES IN- C VOLVING X(I), THE STEP SIZE FIRST TRIED IS C V(DLTFDJ) * MAX(ABS(X(I)), 1/D(I)), C WHERE D IS THE CURRENT SCALE VECTOR (SEE REF. 1). (IF C THIS STEP IS TOO BIG, I.E., IF CALCR SETS NF TO 0, THEN C SMALLER STEPS ARE TRIED UNTIL THE STEP SIZE IS SHRUNK BE- C LOW 1000 * MACHEP, WHERE MACHEP IS THE UNIT ROUNDOFF. C DEFAULT = MACHEP**0.5. C C *** REFERENCES *** C C 1. DENNIS, J.E., GAY, D.M., AND WELSCH, R.E. (1980), AN ADAPTIVE C NONLINEAR LEAST-SQUARES ALGORITHM, SUBMITTED TO ACM TRANS. C MATH. SOFTWARE. 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+++++++++++++++++++++++++++ DECLARATIONS ++++++++++++++++++++++++++++ C C/ C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C C EXTERNAL DFAULT, ITSMRY, NL2ITR, RMDCON, VSCOPY C REAL RMDCON C C DFAULT... SUPPLIES DEFAULT PARAMETER VALUES. C ITSMRY... PRINTS ITERATION SUMMARY AND INFO ABOUT INITIAL AND FINAL X. C NL2ITR... REVERSE-COMMUNICATION ROUTINE THAT CARRIES OUT NL2SOL ALGO- C RITHM. C RMDCON... RETURNS MACHINE-DEPENDENT CONSTANTS. C VSCOPY... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR. C C LOGICAL STRTED C INTEGER DK, D1, I, J1, J1K, K, NF, RN, R1 C REAL H, HFAC, HLIM, NEGPT5, ONE, XK, ZERO C C *** SUBSCRIPTS FOR IV AND V *** C C INTEGER COVPRT, COVREQ, D, DLTFDJ, DTYPE, J, NFCALL, NFGCAL, R, C 1 TOOBIG C DATA HFAC/1.0E3/, HLIM/0.0E0/, NEGPT5/-0.5E0/, + ONE/1.0E0/, ZERO/0.0E0/ C C *** IV SUBSCRIPT VALUES *** C DATA COVPRT/14/, COVREQ/15/, D/27/, DTYPE/16/, J/33/, + NFCALL/6/, NFGCAL/7/, R/50/, TOOBIG/2/ C C *** V SUBSCRIPT VALUES *** C DATA DLTFDJ/36/ C C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ C D1 = 94 + 2*N + P*(3*P + 31)/2 IV(D) = D1 R1 = D1 + P IV(R) = R1 J1 = R1 + N IV(J) = J1 RN = J1 - 1 IF (IV(1) .EQ. 0) CALL DFAULT(IV, V) IV(COVREQ) = -ABS(IV(COVREQ)) IF (IV(COVPRT) .NE. 0 .AND. IV(COVREQ) .EQ. 0) IV(COVREQ) = -1 STRTED = .TRUE. IF (IV(1) .NE. 12) GO TO 80 STRTED = .FALSE. IV(NFCALL) = 1 IV(NFGCAL) = 1 C *** INITIALIZE SCALE VECTOR D TO ONES FOR COMPUTING C *** INITIAL JACOBIAN. IF (IV(DTYPE) .GT. 0) CALL VSCOPY(P, V(D1), ONE) C 10 NF = IV(NFCALL) CALL CALCR(N, P, X, NF, V(R1), UIPARM, URPARM, UFPARM) IF (STRTED) GO TO 20 IF (NF .GT. 0) GO TO 30 IV(1) = 13 GO TO 90 C 20 IF (NF .LE. 0) IV(TOOBIG) = 1 GO TO 80 C C *** COMPUTE FINITE-DIFFERENCE JACOBIAN *** C 30 J1K = J1 DK = D1 DO 70 K = 1, P XK = X(K) H = V(DLTFDJ) * MAX(ABS(XK), ONE/V(DK)) DK = DK + 1 40 X(K) = XK + H NF = IV(NFGCAL) CALL CALCR (N, P, X, NF, V(J1K), UIPARM, URPARM, UFPARM) IF (NF .GT. 0) GO TO 50 IF (HLIM .EQ. ZERO) HLIM = HFAC * RMDCON(3) C *** HLIM = HFAC TIMES THE UNIT ROUNDOFF *** H = NEGPT5 * H IF (ABS(H) .GE. HLIM) GO TO 40 IV(1) = 15 GO TO 90 50 X(K) = XK DO 60 I = R1, RN V(J1K) = (V(J1K) - V(I)) / H J1K = J1K + 1 60 CONTINUE 70 CONTINUE C STRTED = .TRUE. C 80 CALL NL2ITR(V(D1), IV, V(J1), N, N, P, V(R1), V, X) IF (IV(1) - 2) 10, 30, 999 C 90 CALL ITSMRY(V(D1), IV, P, V, X) C 999 RETURN C *** LAST CARD OF NL2SNO FOLLOWS *** END *NL2SOL SUBROUTINE NL2SOL(N, P, X, CALCR, CALCJ, IV, V, UIPARM, URPARM, + UFPARM) C C *** MINIMIZE NONLINEAR SUM OF SQUARES USING ANALYTIC JACOBIAN *** C *** (NL2SOL VERSION 2.2) *** C C VARIABLE DECLARATIONS C C SCALAR ARGUMENTS INTEGER + N,P C C ARRAY ARGUMENTS REAL + URPA