      SUBROUTINE  A7SST(IV, LIV, LV, V)
C
C  ***  ASSESS CANDIDATE STEP (***SOL VERSION 2.3)  ***
C
      INTEGER LIV, LV
      INTEGER IV(LIV)
      REAL V(LV)
C
C  ***  PURPOSE  ***
C
C        THIS SUBROUTINE IS CALLED BY AN UNCONSTRAINED MINIMIZATION
C     ROUTINE TO ASSESS THE NEXT CANDIDATE STEP.  IT MAY RECOMMEND ONE
C     OF SEVERAL COURSES OF ACTION, SUCH AS ACCEPTING THE STEP, RECOM-
C     PUTING IT USING THE SAME OR A NEW QUADRATIC MODEL, OR HALTING DUE
C     TO CONVERGENCE OR FALSE CONVERGENCE.  SEE THE RETURN CODE LISTING
C     BELOW.
C
C--------------------------  PARAMETER USAGE  --------------------------
C
C  IV (I/O) INTEGER PARAMETER AND SCRATCH VECTOR -- SEE DESCRIPTION
C             BELOW OF IV VALUES REFERENCED.
C LIV (IN)  LENGTH OF IV ARRAY.
C  LV (IN)  LENGTH OF V ARRAY.
C   V (I/O) REAL PARAMETER AND SCRATCH VECTOR -- SEE DESCRIPTION
C             BELOW OF V VALUES REFERENCED.
C
C  ***  IV VALUES REFERENCED  ***
C
C    IV(IRC) (I/O) ON INPUT FOR THE FIRST STEP TRIED IN A NEW ITERATION,
C             IV(IRC) SHOULD BE SET TO 3 OR 4 (THE VALUE TO WHICH IT IS
C             SET WHEN STEP IS DEFINITELY TO BE ACCEPTED).  ON INPUT
C             AFTER STEP HAS BEEN RECOMPUTED, IV(IRC) SHOULD BE
C             UNCHANGED SINCE THE PREVIOUS RETURN OF  A7SST.
C                ON OUTPUT, IV(IRC) IS A RETURN CODE HAVING ONE OF THE
C             FOLLOWING VALUES...
C                  1 = SWITCH MODELS OR TRY SMALLER STEP.
C                  2 = SWITCH MODELS OR ACCEPT STEP.
C                  3 = ACCEPT STEP AND DETERMINE V(RADFAC) BY GRADIENT
C                       TESTS.
C                  4 = ACCEPT STEP, V(RADFAC) HAS BEEN DETERMINED.
C                  5 = RECOMPUTE STEP (USING THE SAME MODEL).
C                  6 = RECOMPUTE STEP WITH RADIUS = V(LMAXS) BUT DO NOT
C                       EVAULATE THE OBJECTIVE FUNCTION.
C                  7 = X-CONVERGENCE (SEE V(XCTOL)).
C                  8 = RELATIVE FUNCTION CONVERGENCE (SEE V(RFCTOL)).
C                  9 = BOTH X- AND RELATIVE FUNCTION CONVERGENCE.
C                 10 = ABSOLUTE FUNCTION CONVERGENCE (SEE V(AFCTOL)).
C                 11 = SINGULAR CONVERGENCE (SEE V(LMAXS)).
C                 12 = FALSE CONVERGENCE (SEE V(XFTOL)).
C                 13 = IV(IRC) WAS OUT OF RANGE ON INPUT.
C             RETURN CODE I HAS PRECDENCE OVER I+1 FOR I = 9, 10, 11.
C IV(MLSTGD) (I/O) SAVED VALUE OF IV(MODEL).
C  IV(MODEL) (I/O) ON INPUT, IV(MODEL) SHOULD BE AN INTEGER IDENTIFYING
C             THE CURRENT QUADRATIC MODEL OF THE OBJECTIVE FUNCTION.
C             IF A PREVIOUS STEP YIELDED A BETTER FUNCTION REDUCTION,
C             THEN IV(MODEL) WILL BE SET TO IV(MLSTGD) ON OUTPUT.
C IV(NFCALL) (IN)  INVOCATION COUNT FOR THE OBJECTIVE FUNCTION.
C IV(NFGCAL) (I/O) VALUE OF IV(NFCALL) AT STEP THAT GAVE THE BIGGEST
C             FUNCTION REDUCTION THIS ITERATION.  IV(NFGCAL) REMAINS
C             UNCHANGED UNTIL A FUNCTION REDUCTION IS OBTAINED.
C IV(RADINC) (I/O) THE NUMBER OF RADIUS INCREASES (OR MINUS THE NUMBER
C             OF DECREASES) SO FAR THIS ITERATION.
C IV(RESTOR) (OUT) SET TO 1 IF V(F) HAS BEEN RESTORED AND X SHOULD BE
C             RESTORED TO ITS INITIAL VALUE, TO 2 IF X SHOULD BE SAVED,
C             TO 3 IF X SHOULD BE RESTORED FROM THE SAVED VALUE, AND TO
C             0 OTHERWISE.
C  IV(STAGE) (I/O) COUNT OF THE NUMBER OF MODELS TRIED SO FAR IN THE
C             CURRENT ITERATION.
C IV(STGLIM) (IN)  MAXIMUM NUMBER OF MODELS TO CONSIDER.
C IV(SWITCH) (OUT) SET TO 0 UNLESS A NEW MODEL IS BEING TRIED AND IT
C             GIVES A SMALLER FUNCTION VALUE THAN THE PREVIOUS MODEL,
C             IN WHICH CASE  A7SST 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) AND  A7SST DOES NOT RETURN WITH
C             IV(IRC) = 11, THEN  A7SST RETURNS WITH IV(IRC) = 10.
C V(DECFAC) (IN)  FACTOR BY WHICH TO DECREASE RADIUS WHEN IV(TOOBIG) IS
C             NONZERO.
C V(DSTNRM) (IN)  THE 2-NORM OF D*STEP.
C V(DSTSAV) (I/O) VALUE OF V(DSTNRM) ON SAVED STEP.
C   V(DST0) (IN)  THE 2-NORM OF D TIMES THE NEWTON STEP (WHEN DEFINED,
C             I.E., FOR V(NREDUC) .GE. 0).
C      V(F) (I/O) ON BOTH INPUT AND OUTPUT, V(F) IS THE OBJECTIVE FUNC-
C             TION VALUE AT X.  IF X IS RESTORED TO A PREVIOUS VALUE,
C             THEN V(F) IS RESTORED TO THE CORRESPONDING VALUE.
C   V(FDIF) (OUT) THE FUNCTION REDUCTION V(F0) - V(F) (FOR THE OUTPUT
C             VALUE OF V(F) IF AN EARLIER STEP GAVE A BIGGER FUNCTION
C             DECREASE, AND FOR THE INPUT VALUE OF V(F) OTHERWISE).
C V(FLSTGD) (I/O) SAVED VALUE OF V(F).
C     V(F0) (IN)  OBJECTIVE FUNCTION VALUE AT START OF ITERATION.
C V(GTSLST) (I/O) VALUE OF V(GTSTEP) ON SAVED STEP.
C V(GTSTEP) (IN)  INNER PRODUCT BETWEEN STEP AND GRADIENT.
C V(INCFAC) (IN)  MINIMUM FACTOR BY WHICH TO INCREASE RADIUS.
C  V(LMAXS) (IN)  MAXIMUM REASONABLE STEP SIZE (AND INITIAL STEP BOUND).
C             IF THE ACTUAL FUNCTION DECREASE IS NO MORE THAN TWICE
C             WHAT WAS PREDICTED, IF A RETURN WITH IV(IRC) = 7, 8, OR 9
C             DOES NOT OCCUR, IF V(DSTNRM) .GT. V(LMAXS) OR THE CURRENT
C             STEP IS A NEWTON STEP, AND IF
C             V(PREDUC) .LE. V(SCTOL) * ABS(V(F0)), THEN  A7SST RETURNS
C             WITH IV(IRC) = 11.  IF SO DOING APPEARS WORTHWHILE, THEN
C             A7SST REPEATS THIS TEST (DISALLOWING A FULL NEWTON STEP)
C             WITH V(PREDUC) COMPUTED FOR A STEP OF LENGTH V(LMAXS)
C             (BY A RETURN WITH IV(IRC) = 6).
C V(NREDUC) (I/O)  FUNCTION REDUCTION PREDICTED BY QUADRATIC MODEL FOR
C             NEWTON STEP.  IF  A7SST IS CALLED WITH IV(IRC) = 6, I.E.,
C             IF V(PREDUC) HAS BEEN COMPUTED WITH RADIUS = V(LMAXS) FOR
C             USE IN THE SINGULAR CONVERVENCE TEST, THEN V(NREDUC) IS
C             SET TO -V(PREDUC) BEFORE THE LATTER IS RESTORED.
C V(PLSTGD) (I/O) VALUE OF V(PREDUC) ON SAVED STEP.
C V(PREDUC) (I/O) FUNCTION REDUCTION PREDICTED BY QUADRATIC MODEL FOR
C             CURRENT STEP.
C V(RADFAC) (OUT) FACTOR TO BE USED IN DETERMINING THE NEW RADIUS,
C             WHICH SHOULD BE V(RADFAC)*DST, WHERE  DST  IS EITHER THE
C             OUTPUT VALUE OF V(DSTNRM) OR THE 2-NORM OF
C             DIAG(NEWD)*STEP  FOR THE OUTPUT VALUE OF STEP AND THE
C             UPDATED VERSION, NEWD, OF THE SCALE VECTOR D.  FOR
C             IV(IRC) = 3, V(RADFAC) = 1.0 IS RETURNED.
C V(RDFCMN) (IN)  MINIMUM VALUE FOR V(RADFAC) IN TERMS OF THE INPUT
C             VALUE OF V(DSTNRM) -- SUGGESTED VALUE = 0.1.
C V(RDFCMX) (IN)  MAXIMUM VALUE FOR V(RADFAC) -- SUGGESTED VALUE = 4.0.
C  V(RELDX) (IN) SCALED RELATIVE CHANGE IN X CAUSED BY STEP, COMPUTED
C             (E.G.) BY FUNCTION   RLDST  AS
C                 MAX (D(I)*ABS(X(I)-X0(I)), 1 .LE. I .LE. P) /
C                    MAX (D(I)*(ABS(X(I))+ABS(X0(I))), 1 .LE. I .LE. P).
C V(RFCTOL) (IN)  RELATIVE FUNCTION CONVERGENCE TOLERANCE.  IF THE
C             ACTUAL FUNCTION REDUCTION IS AT MOST TWICE WHAT WAS PRE-
C             DICTED AND  V(NREDUC) .LE. V(RFCTOL)*ABS(V(F0)),  THEN
C             A7SST RETURNS WITH IV(IRC) = 8 OR 9.
C  V(SCTOL) (IN)  SINGULAR CONVERGENCE TOLERANCE -- SEE V(LMAXS).
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             A7SST 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  A7SST 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     A7SST 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  A7SST AGAIN, IN WHICH
C     CASE THE STOPPING TESTS WILL BE REPEATED.
C
C  ***  REFERENCES  ***
C
C     (1) DENNIS, J.E., JR., GAY, D.M., AND WELSCH, R.E. (1981),
C        AN ADAPTIVE NONLINEAR LEAST-SQUARES ALGORITHM,
C        ACM TRANS. MATH. SOFTWARE, VOL. 7, NO. 3.
C
C     (2) POWELL, M.J.D. (1970)  A FORTRAN SUBROUTINE FOR SOLVING
C        SYSTEMS OF NONLINEAR ALGEBRAIC EQUATIONS, IN NUMERICAL
C        METHODS FOR NONLINEAR ALGEBRAIC EQUATIONS, EDITED BY
C        P. RABINOWITZ, GORDON AND BREACH, LONDON.
C
C  ***  HISTORY  ***
C
C        JOHN DENNIS DESIGNED MUCH OF THIS ROUTINE, STARTING WITH
C     IDEAS IN (2). ROY WELSCH SUGGESTED THE MODEL SWITCHING STRATEGY.
C        DAVID GAY AND STEPHEN PETERS CAST THIS SUBROUTINE INTO A MORE
C     PORTABLE FORM (WINTER 1977), AND DAVID GAY CAST IT INTO ITS
C     PRESENT FORM (FALL 1978), WITH MINOR CHANGES TO THE SINGULAR
C     CONVERGENCE TEST IN MAY, 1984 (TO DEAL WITH FULL NEWTON STEPS).
C
C  ***  GENERAL  ***
C
C     THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH
C     SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS
C     MCS-7600324, DCR75-10143, 76-14311DSS, MCS76-11989, AND
C     MCS-7906671.
C
C------------------------  EXTERNAL QUANTITIES  ------------------------
C
C  ***  NO EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
C--------------------------  LOCAL VARIABLES  --------------------------
C
      LOGICAL GOODX
      INTEGER I, NFC
      REAL EMAX, EMAXS, GTS, RFAC1, XMAX
      REAL HALF, ONE, ONEP2, TWO, ZERO
C
C  ***  SUBSCRIPTS FOR IV AND V  ***
C
      INTEGER AFCTOL, DECFAC, DSTNRM, DSTSAV, DST0, F, FDIF, FLSTGD, F0,
     1        GTSLST, GTSTEP, INCFAC, IRC, LMAXS, MLSTGD, MODEL, NFCALL,
     2        NFGCAL, NREDUC, PLSTGD, PREDUC, RADFAC, RADINC, RDFCMN,
     3        RDFCMX, RELDX, RESTOR, RFCTOL, SCTOL, STAGE, STGLIM,
     4        STPPAR, SWITCH, TOOBIG, TUNER1, TUNER2, TUNER3, XCTOL,
     5        XFTOL, XIRC
C
C  ***  DATA INITIALIZATIONS  ***
C
      PARAMETER (HALF=0.5E+0, ONE=1.E+0, ONEP2=1.2E+0, TWO=2.E+0,
     1           ZERO=0.E+0)
C
      PARAMETER (IRC=29, MLSTGD=32, MODEL=5, NFCALL=6, NFGCAL=7,
     1           RADINC=8, RESTOR=9, STAGE=10, STGLIM=11, SWITCH=12,
     2           TOOBIG=2, XIRC=13)
      PARAMETER (AFCTOL=31, DECFAC=22, DSTNRM=2, DST0=3, DSTSAV=18,
     1           F=10, FDIF=11, FLSTGD=12, F0=13, GTSLST=14, GTSTEP=4,
     2           INCFAC=23, LMAXS=36, NREDUC=6, PLSTGD=15, PREDUC=7,
     3           RADFAC=16, RDFCMN=24, RDFCMX=25, RELDX=17, RFCTOL=32,
     4           SCTOL=37, STPPAR=5, TUNER1=26, TUNER2=27, TUNER3=28,
     5           XCTOL=33, XFTOL=34)
C
C+++++++++++++++++++++++++++++++  BODY  ++++++++++++++++++++++++++++++++
C
      NFC = IV(NFCALL)
      IV(SWITCH) = 0
      IV(RESTOR) = 0
      RFAC1 = ONE
      GOODX = .TRUE.
      I = IV(IRC)
      IF (I .GE. 1 .AND. I .LE. 12)
     1             GO TO (20,30,10,10,40,280,220,220,220,220,220,170), I
         IV(IRC) = 13
         GO TO 999
C
C  ***  INITIALIZE FOR NEW ITERATION  ***
C
 10   IV(STAGE) = 1
      IV(RADINC) = 0
      V(FLSTGD) = V(F0)
      IF (IV(TOOBIG) .EQ. 0) GO TO 110
         IV(STAGE) = -1
         IV(XIRC) = I
         GO TO 60
C
C  ***  STEP WAS RECOMPUTED WITH NEW MODEL OR SMALLER RADIUS  ***
C  ***  FIRST DECIDE WHICH  ***
C
 20   IF (IV(MODEL) .NE. IV(MLSTGD)) GO TO 30
C        ***  OLD MODEL RETAINED, SMALLER RADIUS TRIED  ***
C        ***  DO NOT CONSIDER ANY MORE NEW MODELS THIS ITERATION  ***
         IV(STAGE) = IV(STGLIM)
         IV(RADINC) = -1
         GO TO 110
C
C  ***  A NEW MODEL IS BEING TRIED.  DECIDE WHETHER TO KEEP IT.  ***
C
 30   IV(STAGE) = IV(STAGE) + 1
C
C     ***  NOW WE ADD THE POSSIBILTIY THAT STEP WAS RECOMPUTED WITH  ***
C     ***  THE SAME MODEL, PERHAPS BECAUSE OF AN OVERSIZED STEP.     ***
C
 40   IF (IV(STAGE) .GT. 0) GO TO 50
C
C        ***  STEP WAS RECOMPUTED BECAUSE IT WAS TOO BIG.  ***
C
         IF (IV(TOOBIG) .NE. 0) GO TO 60
C
C        ***  RESTORE IV(STAGE) AND PICK UP WHERE WE LEFT OFF.  ***
C
         IV(STAGE) = -IV(STAGE)
         I = IV(XIRC)
         GO TO (20, 30, 110, 110, 70), I
C
 50   IF (IV(TOOBIG) .EQ. 0) GO TO 70
C
C  ***  HANDLE OVERSIZE STEP  ***
C
      IF (IV(RADINC) .GT. 0) GO TO 80
         IV(STAGE) = -IV(STAGE)
         IV(XIRC) = IV(IRC)
C
 60      V(RADFAC) = V(DECFAC)
         IV(RADINC) = IV(RADINC) - 1
         IV(IRC) = 5
         IV(RESTOR) = 1
         GO TO 999
C
 70   IF (V(F) .LT. V(FLSTGD)) GO TO 110
C
C     *** THE NEW STEP IS A LOSER.  RESTORE OLD MODEL.  ***
C
      IF (IV(MODEL) .EQ. IV(MLSTGD)) GO TO 80
         IV(MODEL) = IV(MLSTGD)
         IV(SWITCH) = 1
C
C     ***  RESTORE STEP, ETC. ONLY IF A PREVIOUS STEP DECREASED V(F).
C
 80   IF (V(FLSTGD) .GE. V(F0)) GO TO 110
         IV(RESTOR) = 1
         V(F) = V(FLSTGD)
         V(PREDUC) = V(PLSTGD)
         V(GTSTEP) = V(GTSLST)
         IF (IV(SWITCH) .EQ. 0) RFAC1 = V(DSTNRM) / V(DSTSAV)
         V(DSTNRM) = V(DSTSAV)
         NFC = IV(NFGCAL)
         GOODX = .FALSE.
C
 110  V(FDIF) = V(F0) - V(F)
      IF (V(FDIF) .GT. V(TUNER2) * V(PREDUC)) GO TO 140
      IF (IV(RADINC) .GT. 0) GO TO 140
C
C        ***  NO (OR ONLY A TRIVIAL) FUNCTION DECREASE
C        ***  -- SO TRY NEW MODEL OR SMALLER RADIUS
C
         IF (V(F) .LT. V(F0)) GO TO 120
              IV(MLSTGD) = IV(MODEL)
              V(FLSTGD) = V(F)
              V(F) = V(F0)
              IV(RESTOR) = 1
              GO TO 130
 120     IV(NFGCAL) = NFC
 130     IV(IRC) = 1
         IF (IV(STAGE) .LT. IV(STGLIM)) GO TO 160
              IV(IRC) = 5
              IV(RADINC) = IV(RADINC) - 1
              GO TO 160
C
C  ***  NONTRIVIAL FUNCTION DECREASE ACHIEVED  ***
C
 140  IV(NFGCAL) = NFC
      RFAC1 = ONE
      V(DSTSAV) = V(DSTNRM)
      IF (V(FDIF) .GT. V(PREDUC)*V(TUNER1)) GO TO 190
C
C  ***  DECREASE WAS MUCH LESS THAN PREDICTED -- EITHER CHANGE MODELS
C  ***  OR ACCEPT STEP WITH DECREASED RADIUS.
C
      IF (IV(STAGE) .GE. IV(STGLIM)) GO TO 150
C        ***  CONSIDER SWITCHING MODELS  ***
         IV(IRC) = 2
         GO TO 160
C
C     ***  ACCEPT STEP WITH DECREASED RADIUS  ***
C
 150  IV(IRC) = 4
C
C  ***  SET V(RADFAC) TO FLETCHER*S DECREASE FACTOR  ***
C
 160  IV(XIRC) = IV(IRC)
      EMAX = V(GTSTEP) + V(FDIF)
      V(RADFAC) = HALF * RFAC1
      IF (EMAX .LT. V(GTSTEP)) V(RADFAC) = RFAC1 *   MAX(V(RDFCMN),
     1                                           HALF * V(GTSTEP)/EMAX)
C
C  ***  DO FALSE CONVERGENCE TEST  ***
C
 170  IF (V(RELDX) .LE. V(XFTOL)) GO TO 180
         IV(IRC) = IV(XIRC)
         IF (V(F) .LT. V(F0)) GO TO 200
              GO TO 230
C
 180  IV(IRC) = 12
      GO TO 240
C
C  ***  HANDLE GOOD FUNCTION DECREASE  ***
C
 190  IF (V(FDIF) .LT. (-V(TUNER3) * V(GTSTEP))) GO TO 210
C
C     ***  INCREASING RADIUS LOOKS WORTHWHILE.  SEE IF WE JUST
C     ***  RECOMPUTED STEP WITH A DECREASED RADIUS OR RESTORED STEP
C     ***  AFTER RECOMPUTING IT WITH A LARGER RADIUS.
C
      IF (IV(RADINC) .LT. 0) GO TO 210
      IF (IV(RESTOR) .EQ. 1) GO TO 210
C
C        ***  WE DID NOT.  TRY A LONGER STEP UNLESS THIS WAS A NEWTON
C        ***  STEP.
C
         V(RADFAC) = V(RDFCMX)
         GTS = V(GTSTEP)
         IF (V(FDIF) .LT. (HALF/V(RADFAC) - ONE) * GTS)
     1            V(RADFAC) =   MAX(V(INCFAC), HALF*GTS/(GTS + V(FDIF)))
         IV(IRC) = 4
         IF (V(STPPAR) .EQ. ZERO) GO TO 230
         IF (V(DST0) .GE. ZERO .AND. (V(DST0) .LT. TWO*V(DSTNRM)
     1             .OR. V(NREDUC) .LT. ONEP2*V(FDIF)))  GO TO 230
C             ***  STEP WAS NOT A NEWTON STEP.  RECOMPUTE IT WITH
C             ***  A LARGER RADIUS.
              IV(IRC) = 5
              IV(RADINC) = IV(RADINC) + 1
C
C  ***  SAVE VALUES CORRESPONDING TO GOOD STEP  ***
C
 200  V(FLSTGD) = V(F)
      IV(MLSTGD) = IV(MODEL)
      IF (IV(RESTOR) .NE. 1) IV(RESTOR) = 2
      V(DSTSAV) = V(DSTNRM)
      IV(NFGCAL) = NFC
      V(PLSTGD) = V(PREDUC)
      V(GTSLST) = V(GTSTEP)
      GO TO 230
C
C  ***  ACCEPT STEP WITH RADIUS UNCHANGED  ***
C
 210  V(RADFAC) = ONE
      IV(IRC) = 3
      GO TO 230
C
C  ***  COME HERE FOR A RESTART AFTER CONVERGENCE  ***
C
 220  IV(IRC) = IV(XIRC)
      IF (V(DSTSAV) .GE. ZERO) GO TO 240
         IV(IRC) = 12
         GO TO 240
C
C  ***  PERFORM CONVERGENCE TESTS  ***
C
 230  IV(XIRC) = IV(IRC)
 240  IF (IV(RESTOR) .EQ. 1 .AND. V(FLSTGD) .LT. V(F0)) IV(RESTOR) = 3
      IF ( ABS(V(F)) .LT. V(AFCTOL)) IV(IRC) = 10
      IF (HALF * V(FDIF) .GT. V(PREDUC)) GO TO 999
      EMAX = V(RFCTOL) *  ABS(V(F0))
      EMAXS = V(SCTOL) *  ABS(V(F0))
      IF (V(PREDUC) .LE. EMAXS .AND. (V(DSTNRM) .GT. V(LMAXS) .OR.
     1     V(STPPAR) .EQ. ZERO)) IV(IRC) = 11
      IF (V(DST0) .LT. ZERO) GO TO 250
      I = 0
      IF ((V(NREDUC) .GT. ZERO .AND. V(NREDUC) .LE. EMAX) .OR.
     1    (V(NREDUC) .EQ. ZERO. AND. V(PREDUC) .EQ. ZERO))  I = 2
      IF (V(STPPAR) .EQ. ZERO .AND. V(RELDX) .LE. V(XCTOL)
     1                        .AND. GOODX)                  I = I + 1
      IF (I .GT. 0) IV(IRC) = I + 6
C
C  ***  CONSIDER RECOMPUTING STEP OF LENGTH V(LMAXS) FOR SINGULAR
C  ***  CONVERGENCE TEST.
C
 250  IF (IV(IRC) .GT. 5 .AND. IV(IRC) .NE. 12) GO TO 999
      IF (V(STPPAR) .EQ. ZERO) GO TO 999
      IF (V(DSTNRM) .GT. V(LMAXS)) GO TO 260
         IF (V(PREDUC) .GE. EMAXS) GO TO 999
              IF (V(DST0) .LE. ZERO) GO TO 270
                   IF (HALF * V(DST0) .LE. V(LMAXS)) GO TO 999
                        GO TO 270
 260  IF (HALF * V(DSTNRM) .LE. V(LMAXS)) GO TO 999
      XMAX = V(LMAXS) / V(DSTNRM)
      IF (XMAX * (TWO - XMAX) * V(PREDUC) .GE. EMAXS) GO TO 999
 270  IF (V(NREDUC) .LT. ZERO) GO TO 290
C
C  ***  RECOMPUTE V(PREDUC) FOR USE IN SINGULAR CONVERGENCE TEST  ***
C
      V(GTSLST) = V(GTSTEP)
      V(DSTSAV) = V(DSTNRM)
      IF (IV(IRC) .EQ. 12) V(DSTSAV) = -V(DSTSAV)
      V(PLSTGD) = V(PREDUC)
      I = IV(RESTOR)
      IV(RESTOR) = 2
      IF (I .EQ. 3) IV(RESTOR) = 0
      IV(IRC) = 6
      GO TO 999
C
C  ***  PERFORM SINGULAR CONVERGENCE TEST WITH RECOMPUTED V(PREDUC)  ***
C
 280  V(GTSTEP) = V(GTSLST)
      V(DSTNRM) =  ABS(V(DSTSAV))
      IV(IRC) = IV(XIRC)
      IF (V(DSTSAV) .LE. ZERO) IV(IRC) = 12
      V(NREDUC) = -V(PREDUC)
      V(PREDUC) = V(PLSTGD)
      IV(RESTOR) = 3
 290  IF (-V(NREDUC) .LE. V(SCTOL) *  ABS(V(F0))) IV(IRC) = 11
C
 999  RETURN
C
C  ***  LAST LINE OF  A7SST FOLLOWS  ***
      END
      REAL FUNCTION  D7TPR(P, X, Y)
C
C  ***  RETURN THE INNER PRODUCT OF THE P-VECTORS X AND Y.  ***
C
      INTEGER P
      REAL X(P), Y(P)
C
      INTEGER I
C  ***  ACTIVATE THE *'ED COMMENT LINES BELOW IF UNDERFLOW IS A PROBLEM.
C  ***   R7MDC(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
*     REAL  R7MDC
*     EXTERNAL  R7MDC
      REAL ONE, ZERO
      PARAMETER (ONE=1.E+0, ZERO=0.E+0)
*     REAL SQTETA, T
*      DATA SQTETA/0.E+0/
C
       D7TPR = ZERO
*      IF (P .LE. 0) GO TO 999
*      IF (SQTETA .EQ. ZERO) SQTETA =  R7MDC(2)
      DO 20 I = 1, P
*         T = AMAX1( 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       D7TPR =  D7TPR + X(I)*Y(I)
 20   CONTINUE
C
 999  RETURN
C  ***  LAST LINE OF  D7TPR FOLLOWS  ***
      END
      SUBROUTINE  D7UP5(D, IV, LIV, LV, P, PS, V)
C
C  ***  UPDATE SCALE VECTOR D FOR  G7LIT  ***
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER LIV, LV, P, PS
      INTEGER IV(LIV)
      REAL D(P), V(LV)
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER D0, HII, I, JTOLI, JTOL0, R1I, S1
      REAL T, VDFAC
C
C     ***  CONSTANTS  ***
      REAL ZERO
C
C     ***  EXTERNAL FUNCTIONS  ***
C
      EXTERNAL  D7TPR
      REAL  D7TPR
C
C  ***  SUBSCRIPTS FOR IV AND V  ***
C
      INTEGER DFAC, DTYPE, HC, JTOL, NITER, RMAT, S
      PARAMETER (DFAC=41, DTYPE=16, HC=71, JTOL=59, NITER=31, RMAT=78,
     1           S=62)
C
      PARAMETER (ZERO=0.E+0)
C
C  ***  BODY  ***
C
      IF (IV(DTYPE) .NE. 1 .AND. IV(NITER) .GT. 0) GO TO 999
      R1I = IV(RMAT)
      HII = IV(HC) - 1
      VDFAC = V(DFAC)
      JTOL0 = IV(JTOL) - 1
      D0 = JTOL0 + P
      S1 = IV(S) - 1
      DO 30 I = 1, P
         IF (R1I .LE. 0) GO TO 10
             T =  D7TPR(I, V(R1I), V(R1I))
             R1I = R1I + I
             GO TO 20
 10      HII = HII + I
         T =  ABS(V(HII))
 20      S1 = S1 + I
         IF (I .LE. PS) T = T +   MAX(V(S1), ZERO)
         T =  SQRT(T)
         JTOLI = JTOL0 + I
         D0 = D0 + 1
         IF (T .LT. V(JTOLI)) T =   MAX(V(D0), V(JTOLI))
         D(I) =   MAX(VDFAC*D(I), T)
 30      CONTINUE
C
 999  RETURN
C  ***  LAST LINE OF  D7UP5 FOLLOWS  ***
      END
      SUBROUTINE  G7QTS(D, DIG, DIHDI, KA, L, P, STEP, V, W)
C
C  *** COMPUTE GOLDFELD-QUANDT-TROTTER STEP BY MORE-HEBDEN TECHNIQUE ***
C  ***  (NL2SOL VERSION 2.2), MODIFIED A LA MORE AND SORENSEN  ***
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER KA, P
      REAL D(P), DIG(P), DIHDI(1), L(1), V(21), STEP(P),
     1                 W(1)
C     DIMENSION DIHDI(P*(P+1)/2), L(P*(P+1)/2), W(4*P+7)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ***  PURPOSE  ***
C
C        GIVEN THE (COMPACTLY STORED) LOWER TRIANGLE OF A SCALED
C     HESSIAN (APPROXIMATION) AND A NONZERO SCALED GRADIENT VECTOR,
C     THIS SUBROUTINE COMPUTES A GOLDFELD-QUANDT-TROTTER STEP OF
C     APPROXIMATE LENGTH V(RADIUS) BY THE MORE-HEBDEN TECHNIQUE.  IN
C     OTHER WORDS, STEP IS COMPUTED TO (APPROXIMATELY) MINIMIZE
C     PSI(STEP) = (G**T)*STEP + 0.5*(STEP**T)*H*STEP  SUCH THAT THE
C     2-NORM OF D*STEP IS AT MOST (APPROXIMATELY) V(RADIUS), WHERE
C     G  IS THE GRADIENT,  H  IS THE HESSIAN, AND  D  IS A DIAGONAL
C     SCALE MATRIX WHOSE DIAGONAL IS STORED IN THE PARAMETER D.
C     ( G7QTS ASSUMES  DIG = D**-1 * G  AND  DIHDI = D**-1 * H * D**-1.)
C
C  ***  PARAMETER DESCRIPTION  ***
C
C     D (IN)  = THE SCALE VECTOR, I.E. THE DIAGONAL OF THE SCALE
C              MATRIX  D  MENTIONED ABOVE UNDER PURPOSE.
C   DIG (IN)  = THE SCALED GRADIENT VECTOR, D**-1 * G.  IF G = 0, THEN
C              STEP = 0  AND  V(STPPAR) = 0  ARE RETURNED.
C DIHDI (IN)  = LOWER TRIANGLE OF THE SCALED HESSIAN (APPROXIMATION),
C              I.E., D**-1 * H * D**-1, STORED COMPACTLY BY ROWS., I.E.,
C              IN THE ORDER (1,1), (2,1), (2,2), (3,1), (3,2), ETC.
C    KA (I/O) = THE NUMBER OF HEBDEN ITERATIONS (SO FAR) TAKEN TO DETER-
C              MINE STEP.  KA .LT. 0 ON INPUT MEANS THIS IS THE FIRST
C              ATTEMPT TO DETERMINE STEP (FOR THE PRESENT DIG AND DIHDI)
C              -- KA IS INITIALIZED TO 0 IN THIS CASE.  OUTPUT WITH
C              KA = 0  (OR V(STPPAR) = 0)  MEANS  STEP = -(H**-1)*G.
C     L (I/O) = WORKSPACE OF LENGTH P*(P+1)/2 FOR CHOLESKY FACTORS.
C     P (IN)  = NUMBER OF PARAMETERS -- THE HESSIAN IS A  P X P  MATRIX.
C  STEP (I/O) = THE STEP COMPUTED.
C     V (I/O) CONTAINS VARIOUS CONSTANTS AND VARIABLES DESCRIBED BELOW.
C     W (I/O) = WORKSPACE OF LENGTH 4*P + 6.
C
C  ***  ENTRIES IN V  ***
C
C V(DGNORM) (I/O) = 2-NORM OF (D**-1)*G.
C V(DSTNRM) (OUTPUT) = 2-NORM OF D*STEP.
C V(DST0)   (I/O) = 2-NORM OF D*(H**-1)*G (FOR POS. DEF. H ONLY), OR
C             OVERESTIMATE OF SMALLEST EIGENVALUE OF (D**-1)*H*(D**-1).
C V(EPSLON) (IN)  = MAX. REL. ERROR ALLOWED FOR PSI(STEP).  FOR THE
C             STEP RETURNED, PSI(STEP) WILL EXCEED ITS OPTIMAL VALUE
C             BY LESS THAN -V(EPSLON)*PSI(STEP).  SUGGESTED VALUE = 0.1.
C V(GTSTEP) (OUT) = INNER PRODUCT BETWEEN G AND STEP.
C V(NREDUC) (OUT) = PSI(-(H**-1)*G) = PSI(NEWTON STEP)  (FOR POS. DEF.
C             H ONLY -- V(NREDUC) IS SET TO ZERO OTHERWISE).
C V(PHMNFC) (IN)  = TOL. (TOGETHER WITH V(PHMXFC)) FOR ACCEPTING STEP
C             (MORE*S SIGMA).  THE ERROR V(DSTNRM) - V(RADIUS) MUST LIE
C             BETWEEN V(PHMNFC)*V(RADIUS) AND V(PHMXFC)*V(RADIUS).
C V(PHMXFC) (IN)  (SEE V(PHMNFC).)
C             SUGGESTED VALUES -- V(PHMNFC) = -0.25, V(PHMXFC) = 0.5.
C V(PREDUC) (OUT) = PSI(STEP) = PREDICTED OBJ. FUNC. REDUCTION FOR STEP.
C V(RADIUS) (IN)  = RADIUS OF CURRENT (SCALED) TRUST REGION.
C V(RAD0)   (I/O) = VALUE OF V(RADIUS) FROM PREVIOUS CALL.
C V(STPPAR) (I/O) IS NORMALLY THE MARQUARDT PARAMETER, I.E. THE ALPHA
C             DESCRIBED BELOW UNDER ALGORITHM NOTES.  IF H + ALPHA*D**2
C             (SEE ALGORITHM NOTES) IS (NEARLY) SINGULAR, HOWEVER,
C             THEN V(STPPAR) = -ALPHA.
C
C  ***  USAGE NOTES  ***
C
C     IF IT IS DESIRED TO RECOMPUTE STEP USING A DIFFERENT VALUE OF
C     V(RADIUS), THEN THIS ROUTINE MAY BE RESTARTED BY CALLING IT
C     WITH ALL PARAMETERS UNCHANGED EXCEPT V(RADIUS).  (THIS EXPLAINS
C     WHY STEP AND W ARE LISTED AS I/O).  ON AN INITIAL CALL (ONE WITH
C     KA .LT. 0), STEP AND W NEED NOT BE INITIALIZED AND ONLY COMPO-
C     NENTS V(EPSLON), V(STPPAR), V(PHMNFC), V(PHMXFC), V(RADIUS), AND
C     V(RAD0) OF V MUST BE INITIALIZED.
C
C  ***  ALGORITHM NOTES  ***
C
C        THE DESIRED G-Q-T STEP (REF. 2, 3, 4, 6) SATISFIES
C     (H + ALPHA*D**2)*STEP = -G  FOR SOME NONNEGATIVE ALPHA SUCH THAT
C     H + ALPHA*D**2 IS POSITIVE SEMIDEFINITE.  ALPHA AND STEP ARE
C     COMPUTED BY A SCHEME ANALOGOUS TO THE ONE DESCRIBED IN REF. 5.
C     ESTIMATES OF THE SMALLEST AND LARGEST EIGENVALUES OF THE HESSIAN
C     ARE OBTAINED FROM THE GERSCHGORIN CIRCLE THEOREM ENHANCED BY A
C     SIMPLE FORM OF THE SCALING DESCRIBED IN REF. 7.  CASES IN WHICH
C     H + ALPHA*D**2 IS NEARLY (OR EXACTLY) SINGULAR ARE HANDLED BY
C     THE TECHNIQUE DISCUSSED IN REF. 2.  IN THESE CASES, A STEP OF
C     (EXACT) LENGTH V(RADIUS) IS RETURNED FOR WHICH PSI(STEP) EXCEEDS
C     ITS OPTIMAL VALUE BY LESS THAN -V(EPSLON)*PSI(STEP).  THE TEST
C     SUGGESTED IN REF. 6 FOR DETECTING THE SPECIAL CASE IS PERFORMED
C     ONCE TWO MATRIX FACTORIZATIONS HAVE BEEN DONE -- DOING SO SOONER
C     SEEMS TO DEGRADE THE PERFORMANCE OF OPTIMIZATION ROUTINES THAT
C     CALL THIS ROUTINE.
C
C  ***  FUNCTIONS AND SUBROUTINES CALLED  ***
C
C  D7TPR - RETURNS INNER PRODUCT OF TWO VECTORS.
C  L7ITV - APPLIES INVERSE-TRANSPOSE OF COMPACT LOWER TRIANG. MATRIX.
C  L7IVM - APPLIES INVERSE OF COMPACT LOWER TRIANG. MATRIX.
C  L7SRT  - FINDS CHOLESKY FACTOR (OF COMPACTLY STORED LOWER TRIANG.).
C  L7SVN - RETURNS APPROX. TO MIN. SING. VALUE OF LOWER TRIANG. MATRIX.
C  R7MDC - RETURNS MACHINE-DEPENDENT CONSTANTS.
C  V2NRM - RETURNS 2-NORM OF A VECTOR.
C
C  ***  REFERENCES  ***
C
C 1.  DENNIS, J.E., GAY, D.M., AND WELSCH, R.E. (1981), AN ADAPTIVE
C             NONLINEAR LEAST-SQUARES ALGORITHM, ACM TRANS. MATH.
C             SOFTWARE, VOL. 7, NO. 3.
C 2.  GAY, D.M. (1981), COMPUTING OPTIMAL LOCALLY CONSTRAINED STEPS,
C             SIAM J. SCI. STATIST. COMPUTING, VOL. 2, NO. 2, PP.
C             186-197.
C 3.  GOLDFELD, S.M., QUANDT, R.E., AND TROTTER, H.F. (1966),
C             MAXIMIZATION BY QUADRATIC HILL-CLIMBING, ECONOMETRICA 34,
C             PP. 541-551.
C 4.  HEBDEN, M.D. (1973), AN ALGORITHM FOR MINIMIZATION USING EXACT
C             SECOND DERIVATIVES, REPORT T.P. 515, THEORETICAL PHYSICS
C             DIV., A.E.R.E. HARWELL, OXON., ENGLAND.
C 5.  MORE, J.J. (1978), THE LEVENBERG-MARQUARDT ALGORITHM, IMPLEMEN-
C             TATION AND THEORY, PP.105-116 OF SPRINGER LECTURE NOTES
C             IN MATHEMATICS NO. 630, EDITED BY G.A. WATSON, SPRINGER-
C             VERLAG, BERLIN AND NEW YORK.
C 6.  MORE, J.J., AND SORENSEN, D.C. (1981), COMPUTING A TRUST REGION
C             STEP, TECHNICAL REPORT ANL-81-83, ARGONNE NATIONAL LAB.
C 7.  VARGA, R.S. (1965), MINIMAL GERSCHGORIN SETS, PACIFIC J. MATH. 15,
C             PP. 719-729.
C
C  ***  GENERAL  ***
C
C     CODED BY DAVID M. GAY.
C     THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH
C     SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS
C     MCS-7600324, DCR75-10143, 76-14311DSS, MCS76-11989, AND
C     MCS-7906671.
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ***  LOCAL VARIABLES  ***
C
      LOGICAL RESTRT
      INTEGER DGGDMX, DIAG, DIAG0, DSTSAV, EMAX, EMIN, I, IM1, INC, IRC,
     1        J, K, KALIM, KAMIN, K1, LK0, PHIPIN, Q, Q0, UK0, X
      REAL ALPHAK, AKI, AKK, DELTA, DST, EPS, GTSTA, LK,
     1                 OLDPHI, PHI, PHIMAX, PHIMIN, PSIFAC, RAD, RADSQ,
     2                 ROOT, SI, SK, SW, T, TWOPSI, T1, T2, UK, WI
C
C     ***  CONSTANTS  ***
      REAL BIG, DGXFAC, EPSFAC, FOUR, HALF, KAPPA, NEGONE,
     1                 ONE, P001, SIX, THREE, TWO, ZERO
C
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      REAL  D7TPR,  L7SVN,  R7MDC,  V2NRM
      EXTERNAL  D7TPR,  L7ITV,  L7IVM, L7SRT,  L7SVN,  R7MDC,  V2NRM
C
C  ***  SUBSCRIPTS FOR V  ***
C
      INTEGER DGNORM, DSTNRM, DST0, EPSLON, GTSTEP, STPPAR, NREDUC,
     1        PHMNFC, PHMXFC, PREDUC, RADIUS, RAD0
      PARAMETER (DGNORM=1, DSTNRM=2, DST0=3, EPSLON=19, GTSTEP=4,
     1           NREDUC=6, PHMNFC=20, PHMXFC=21, PREDUC=7, RADIUS=8,
     2           RAD0=9, STPPAR=5)
C
      PARAMETER (EPSFAC=50.0E+0, FOUR=4.0E+0, HALF=0.5E+0,
     1     KAPPA=2.0E+0, NEGONE=-1.0E+0, ONE=1.0E+0, P001=1.0E-3,
     2     SIX=6.0E+0, THREE=3.0E+0, TWO=2.0E+0, ZERO=0.0E+0)
      SAVE DGXFAC
      DATA BIG/0.E+0/, DGXFAC/0.E+0/
C
C  ***  BODY  ***
C
      IF (BIG .LE. ZERO) BIG =  R7MDC(6)
C
C     ***  STORE LARGEST ABS. ENTRY IN (D**-1)*H*(D**-1) AT W(DGGDMX).
      DGGDMX = P + 1
C     ***  STORE GERSCHGORIN OVER- AND UNDERESTIMATES OF THE LARGEST
C     ***  AND SMALLEST EIGENVALUES OF (D**-1)*H*(D**-1) AT W(EMAX)
C     ***  AND W(EMIN) RESPECTIVELY.
      EMAX = DGGDMX + 1
      EMIN = EMAX + 1
C     ***  FOR USE IN RECOMPUTING STEP, THE FINAL VALUES OF LK, UK, DST,
C     ***  AND THE INVERSE DERIVATIVE OF MORE*S PHI AT 0 (FOR POS. DEF.
C     ***  H) ARE STORED IN W(LK0), W(UK0), W(DSTSAV), AND W(PHIPIN)
C     ***  RESPECTIVELY.
      LK0 = EMIN + 1
      PHIPIN = LK0 + 1
      UK0 = PHIPIN + 1
      DSTSAV = UK0 + 1
C     ***  STORE DIAG OF (D**-1)*H*(D**-1) IN W(DIAG),...,W(DIAG0+P).
      DIAG0 = DSTSAV
      DIAG = DIAG0 + 1
C     ***  STORE -D*STEP IN W(Q),...,W(Q0+P).
      Q0 = DIAG0 + P
      Q = Q0 + 1
C     ***  ALLOCATE STORAGE FOR SCRATCH VECTOR X  ***
      X = Q + P
      RAD = V(RADIUS)
      RADSQ = RAD**2
C     ***  PHITOL = MAX. ERROR ALLOWED IN DST = V(DSTNRM) = 2-NORM OF
C     ***  D*STEP.
      PHIMAX = V(PHMXFC) * RAD
      PHIMIN = V(PHMNFC) * RAD
      PSIFAC = BIG
      T1 = TWO * V(EPSLON) / (THREE * (FOUR * (V(PHMNFC) + ONE) *
     1                       (KAPPA + ONE)  +  KAPPA  +  TWO) * RAD)
      IF (T1 .LT. BIG*  MIN(RAD,ONE)) PSIFAC = T1 / RAD
C     ***  OLDPHI IS USED TO DETECT LIMITS OF NUMERICAL ACCURACY.  IF
C     ***  WE RECOMPUTE STEP AND IT DOES NOT CHANGE, THEN WE ACCEPT IT.
      OLDPHI = ZERO
      EPS = V(EPSLON)
      IRC = 0
      RESTRT = .FALSE.
      KALIM = KA + 50
C
C  ***  START OR RESTART, DEPENDING ON KA  ***
C
      IF (KA .GE. 0) GO TO 290
C
C  ***  FRESH START  ***
C
      K = 0
      UK = NEGONE
      KA = 0
      KALIM = 50
      V(DGNORM) =  V2NRM(P, DIG)
      V(NREDUC) = ZERO
      V(DST0) = ZERO
      KAMIN = 3
      IF (V(DGNORM) .EQ. ZERO) KAMIN = 0
C
C     ***  STORE DIAG(DIHDI) IN W(DIAG0+1),...,W(DIAG0+P)  ***
C
      J = 0
      DO 10 I = 1, P
         J = J + I
         K1 = DIAG0 + I
         W(K1) = DIHDI(J)
 10      CONTINUE
C
C     ***  DETERMINE W(DGGDMX), THE LARGEST ELEMENT OF DIHDI  ***
C
      T1 = ZERO
      J = P * (P + 1) / 2
      DO 20 I = 1, J
         T =  ABS(DIHDI(I))
         IF (T1 .LT. T) T1 = T
 20      CONTINUE
      W(DGGDMX) = T1
C
C  ***  TRY ALPHA = 0  ***
C
 30   CALL  L7SRT(1, P, L, DIHDI, IRC)
      IF (IRC .EQ. 0) GO TO 50
C        ***  INDEF. H -- UNDERESTIMATE SMALLEST EIGENVALUE, USE THIS
C        ***  ESTIMATE TO INITIALIZE LOWER BOUND LK ON ALPHA.
         J = IRC*(IRC+1)/2
         T = L(J)
         L(J) = ONE
         DO 40 I = 1, IRC
 40           W(I) = ZERO
         W(IRC) = ONE
         CALL  L7ITV(IRC, W, L, W)
         T1 =  V2NRM(IRC, W)
         LK = -T / T1 / T1
         V(DST0) = -LK
         IF (RESTRT) GO TO 210
         GO TO 70
C
C     ***  POSITIVE DEFINITE H -- COMPUTE UNMODIFIED NEWTON STEP.  ***
 50   LK = ZERO
      T =  L7SVN(P, L, W(Q), W(Q))
      IF (T .GE. ONE) GO TO 60
         IF (V(DGNORM) .GE. T*T*BIG) GO TO 70
 60   CALL  L7IVM(P, W(Q), L, DIG)
      GTSTA =  D7TPR(P, W(Q), W(Q))
      V(NREDUC) = HALF * GTSTA
      CALL  L7ITV(P, W(Q), L, W(Q))
      DST =  V2NRM(P, W(Q))
      V(DST0) = DST
      PHI = DST - RAD
      IF (PHI .LE. PHIMAX) GO TO 260
      IF (RESTRT) GO TO 210
C
C  ***  PREPARE TO COMPUTE GERSCHGORIN ESTIMATES OF LARGEST (AND
C  ***  SMALLEST) EIGENVALUES.  ***
C
 70   K = 0
      DO 100 I = 1, P
         WI = ZERO
         IF (I .EQ. 1) GO TO 90
         IM1 = I - 1
         DO 80 J = 1, IM1
              K = K + 1
              T =  ABS(DIHDI(K))
              WI = WI + T
              W(J) = W(J) + T
 80           CONTINUE
 90      W(I) = WI
         K = K + 1
 100     CONTINUE
C
C  ***  (UNDER-)ESTIMATE SMALLEST EIGENVALUE OF (D**-1)*H*(D**-1)  ***
C
      K = 1
      T1 = W(DIAG) - W(1)
      IF (P .LE. 1) GO TO 120
      DO 110 I = 2, P
         J = DIAG0 + I
         T = W(J) - W(I)
         IF (T .GE. T1) GO TO 110
              T1 = T
              K = I
 110     CONTINUE
C
 120  SK = W(K)
      J = DIAG0 + K
      AKK = W(J)
      K1 = K*(K-1)/2 + 1
      INC = 1
      T = ZERO
      DO 150 I = 1, P
         IF (I .EQ. K) GO TO 130
         AKI =  ABS(DIHDI(K1))
         SI = W(I)
         J = DIAG0 + I
         T1 = HALF * (AKK - W(J) + SI - AKI)
         T1 = T1 +  SQRT(T1*T1 + SK*AKI)
         IF (T .LT. T1) T = T1
         IF (I .LT. K) GO TO 140
 130     INC = I
 140     K1 = K1 + INC
 150     CONTINUE
C
      W(EMIN) = AKK - T
      UK = V(DGNORM)/RAD - W(EMIN)
      IF (V(DGNORM) .EQ. ZERO) UK = UK + P001 + P001*UK
      IF (UK .LE. ZERO) UK = P001
C
C  ***  COMPUTE GERSCHGORIN (OVER-)ESTIMATE OF LARGEST EIGENVALUE  ***
C
      K = 1
      T1 = W(DIAG) + W(1)
      IF (P .LE. 1) GO TO 170
      DO 160 I = 2, P
         J = DIAG0 + I
         T = W(J) + W(I)
         IF (T .LE. T1) GO TO 160
              T1 = T
              K = I
 160     CONTINUE
C
 170  SK = W(K)
      J = DIAG0 + K
      AKK = W(J)
      K1 = K*(K-1)/2 + 1
      INC = 1
      T = ZERO
      DO 200 I = 1, P
         IF (I .EQ. K) GO TO 180
         AKI =  ABS(DIHDI(K1))
         SI = W(I)
         J = DIAG0 + I
         T1 = HALF * (W(J) + SI - AKI - AKK)
         T1 = T1 +  SQRT(T1*T1 + SK*AKI)
         IF (T .LT. T1) T = T1
         IF (I .LT. K) GO TO 190
 180     INC = I
 190     K1 = K1 + INC
 200     CONTINUE
C
      W(EMAX) = AKK + T
      LK =   MAX(LK, V(DGNORM)/RAD - W(EMAX))
C
C     ***  ALPHAK = CURRENT VALUE OF ALPHA (SEE ALG. NOTES ABOVE).  WE
C     ***  USE MORE*S SCHEME FOR INITIALIZING IT.
      ALPHAK =  ABS(V(STPPAR)) * V(RAD0)/RAD
      ALPHAK =   MIN(UK,   MAX(ALPHAK, LK))
C
      IF (IRC .NE. 0) GO TO 210
C
C  ***  COMPUTE L0 FOR POSITIVE DEFINITE H  ***
C
      CALL  L7IVM(P, W, L, W(Q))
      T =  V2NRM(P, W)
      W(PHIPIN) = RAD / T / T
      LK =   MAX(LK, PHI*W(PHIPIN))
C
C  ***  SAFEGUARD ALPHAK AND ADD ALPHAK*I TO (D**-1)*H*(D**-1)  ***
C
 210  KA = KA + 1
      IF (-V(DST0) .GE. ALPHAK .OR. ALPHAK .LT. LK .OR. ALPHAK .GE. UK)
     1                      ALPHAK = UK *   MAX(P001,  SQRT(LK/UK))
      IF (ALPHAK .LE. ZERO) ALPHAK = HALF * UK
      IF (ALPHAK .LE. ZERO) ALPHAK = UK
      K = 0
      DO 220 I = 1, P
         K = K + I
         J = DIAG0 + I
         DIHDI(K) = W(J) + ALPHAK
 220     CONTINUE
C
C  ***  TRY COMPUTING CHOLESKY DECOMPOSITION  ***
C
      CALL  L7SRT(1, P, L, DIHDI, IRC)
      IF (IRC .EQ. 0) GO TO 240
C
C  ***  (D**-1)*H*(D**-1) + ALPHAK*I  IS INDEFINITE -- OVERESTIMATE
C  ***  SMALLEST EIGENVALUE FOR USE IN UPDATING LK  ***
C
      J = (IRC*(IRC+1))/2
      T = L(J)
      L(J) = ONE
      DO 230 I = 1, IRC
 230     W(I) = ZERO
      W(IRC) = ONE
      CALL  L7ITV(IRC, W, L, W)
      T1 =  V2NRM(IRC, W)
      LK = ALPHAK - T/T1/T1
      V(DST0) = -LK
      IF (UK .LT. LK) UK = LK
      IF (ALPHAK .LT. LK) GO TO 210
C
C  ***  NASTY CASE -- EXACT GERSCHGORIN BOUNDS.  FUDGE LK, UK...
C
      T = P001 * ALPHAK
      IF (T .LE. ZERO) T = P001
      LK = ALPHAK + T
      IF (UK .LE. LK) UK = LK + T
      GO TO 210
C
C  ***  ALPHAK MAKES (D**-1)*H*(D**-1) POSITIVE DEFINITE.
C  ***  COMPUTE Q = -D*STEP, CHECK FOR CONVERGENCE.  ***
C
 240  CALL  L7IVM(P, W(Q), L, DIG)
      GTSTA =  D7TPR(P, W(Q), W(Q))
      CALL  L7ITV(P, W(Q), L, W(Q))
      DST =  V2NRM(P, W(Q))
      PHI = DST - RAD
      IF (PHI .LE. PHIMAX .AND. PHI .GE. PHIMIN) GO TO 270
      IF (PHI .EQ. OLDPHI) GO TO 270
      OLDPHI = PHI
      IF (PHI .LT. ZERO) GO TO 330
C
C  ***  UNACCEPTABLE ALPHAK -- UPDATE LK, UK, ALPHAK  ***
C
 250  IF (KA .GE. KALIM) GO TO 270
C     ***  THE FOLLOWING   MIN IS NECESSARY BECAUSE OF RESTARTS  ***
      IF (PHI .LT. ZERO) UK =   MIN(UK, ALPHAK)
C     *** KAMIN = 0 ONLY IFF THE GRADIENT VANISHES  ***
      IF (KAMIN .EQ. 0) GO TO 210
      CALL  L7IVM(P, W, L, W(Q))
C     *** THE FOLLOWING, COMMENTED CALCULATION OF ALPHAK IS SOMETIMES
C     *** SAFER BUT WORSE IN PERFORMANCE...
C     T1 = DST /  V2NRM(P, W)
C     ALPHAK = ALPHAK  +  T1 * (PHI/RAD) * T1
      T1 =  V2NRM(P, W)
      ALPHAK = ALPHAK  +  (PHI/T1) * (DST/T1) * (DST/RAD)
      LK =   MAX(LK, ALPHAK)
      ALPHAK = LK
      GO TO 210
C
C  ***  ACCEPTABLE STEP ON FIRST TRY  ***
C
 260  ALPHAK = ZERO
C
C  ***  SUCCESSFUL STEP IN GENERAL.  COMPUTE STEP = -(D**-1)*Q  ***
C
 270  DO 280 I = 1, P
         J = Q0 + I
         STEP(I) = -W(J)/D(I)
 280     CONTINUE
      V(GTSTEP) = -GTSTA
      V(PREDUC) = HALF * ( ABS(ALPHAK)*DST*DST + GTSTA)
      GO TO 410
C
C
C  ***  RESTART WITH NEW RADIUS  ***
C
 290  IF (V(DST0) .LE. ZERO .OR. V(DST0) - RAD .GT. PHIMAX) GO TO 310
C
C     ***  PREPARE TO RETURN NEWTON STEP  ***
C
         RESTRT = .TRUE.
         KA = KA + 1
         K = 0
         DO 300 I = 1, P
              K = K + I
              J = DIAG0 + I
              DIHDI(K) = W(J)
 300          CONTINUE
         UK = NEGONE
         GO TO 30
C
 310  KAMIN = KA + 3
      IF (V(DGNORM) .EQ. ZERO) KAMIN = 0
      IF (KA .EQ. 0) GO TO 50
C
      DST = W(DSTSAV)
      ALPHAK =  ABS(V(STPPAR))
      PHI = DST - RAD
      T = V(DGNORM)/RAD
      UK = T - W(EMIN)
      IF (V(DGNORM) .EQ. ZERO) UK = UK + P001 + P001*UK
      IF (UK .LE. ZERO) UK = P001
      IF (RAD .GT. V(RAD0)) GO TO 320
C
C        ***  SMALLER RADIUS  ***
         LK = ZERO
         IF (ALPHAK .GT. ZERO) LK = W(LK0)
         LK =   MAX(LK, T - W(EMAX))
         IF (V(DST0) .GT. ZERO) LK =   MAX(LK, (V(DST0)-RAD)*W(PHIPIN))
         GO TO 250
C
C     ***  BIGGER RADIUS  ***
 320  IF (ALPHAK .GT. ZERO) UK =   MIN(UK, W(UK0))
      LK =   MAX(ZERO, -V(DST0), T - W(EMAX))
      IF (V(DST0) .GT. ZERO) LK =   MAX(LK, (V(DST0)-RAD)*W(PHIPIN))
      GO TO 250
C
C  ***  DECIDE WHETHER TO CHECK FOR SPECIAL CASE... IN PRACTICE (FROM
C  ***  THE STANDPOINT OF THE CALLING OPTIMIZATION CODE) IT SEEMS BEST
C  ***  NOT TO CHECK UNTIL A FEW ITERATIONS HAVE FAILED -- HENCE THE
C  ***  TEST ON KAMIN BELOW.
C
 330  DELTA = ALPHAK +   MIN(ZERO, V(DST0))
      TWOPSI = ALPHAK*DST*DST + GTSTA
      IF (KA .GE. KAMIN) GO TO 340
C     *** IF THE TEST IN REF. 2 IS SATISFIED, FALL THROUGH TO HANDLE
C     *** THE SPECIAL CASE (AS SOON AS THE MORE-SORENSEN TEST DETECTS
C     *** IT).
      IF (PSIFAC .GE. BIG) GO TO 340
      IF (DELTA .GE. PSIFAC*TWOPSI) GO TO 370
C
C  ***  CHECK FOR THE SPECIAL CASE OF  H + ALPHA*D**2  (NEARLY)
C  ***  SINGULAR.  USE ONE STEP OF INVERSE POWER METHOD WITH START
C  ***  FROM  L7SVN TO OBTAIN APPROXIMATE EIGENVECTOR CORRESPONDING
C  ***  TO SMALLEST EIGENVALUE OF (D**-1)*H*(D**-1).   L7SVN RETURNS
C  ***  X AND W WITH  L*W = X.
C
 340  T =  L7SVN(P, L, W(X), W)
C
C     ***  NORMALIZE W  ***
      DO 350 I = 1, P
 350     W(I) = T*W(I)
C     ***  COMPLETE CURRENT INV. POWER ITER. -- REPLACE W BY (L**-T)*W.
      CALL  L7ITV(P, W, L, W)
      T2 = ONE/ V2NRM(P, W)
      DO 360 I = 1, P
 360     W(I) = T2*W(I)
      T = T2 * T
C
C  ***  NOW W IS THE DESIRED APPROXIMATE (UNIT) EIGENVECTOR AND
C  ***  T*X = ((D**-1)*H*(D**-1) + ALPHAK*I)*W.
C
      SW =  D7TPR(P, W(Q), W)
      T1 = (RAD + DST) * (RAD - DST)
      ROOT =  SQRT(SW*SW + T1)
      IF (SW .LT. ZERO) ROOT = -ROOT
      SI = T1 / (SW + ROOT)
C
C  ***  THE ACTUAL TEST FOR THE SPECIAL CASE...
C
      IF ((T2*SI)**2 .LE. EPS*(DST**2 + ALPHAK*RADSQ)) GO TO 380
C
C  ***  UPDATE UPPER BOUND ON SMALLEST EIGENVALUE (WHEN NOT POSITIVE)
C  ***  (AS RECOMMENDED BY MORE AND SORENSEN) AND CONTINUE...
C
      IF (V(DST0) .LE. ZERO) V(DST0) =   MIN(V(DST0), T2**2 - ALPHAK)
      LK =   MAX(LK, -V(DST0))
C
C  ***  CHECK WHETHER WE CAN HOPE TO DETECT THE SPECIAL CASE IN
C  ***  THE AVAILABLE ARITHMETIC.  ACCEPT STEP AS IT IS IF NOT.
C
C     ***  IF NOT YET AVAILABLE, OBTAIN MACHINE DEPENDENT VALUE DGXFAC.
 370  IF (DGXFAC .EQ. ZERO) DGXFAC = EPSFAC *  R7MDC(3)
C
      IF (DELTA .GT. DGXFAC*W(DGGDMX)) GO TO 250
         GO TO 270
C
C  ***  SPECIAL CASE DETECTED... NEGATE ALPHAK TO INDICATE SPECIAL CASE
C
 380  ALPHAK = -ALPHAK
      V(PREDUC) = HALF * TWOPSI
C
C  ***  ACCEPT CURRENT STEP IF ADDING SI*W WOULD LEAD TO A
C  ***  FURTHER RELATIVE REDUCTION IN PSI OF LESS THAN V(EPSLON)/3.
C
      T1 = ZERO
      T = SI*(ALPHAK*SW - HALF*SI*(ALPHAK + T* D7TPR(P,W(X),W)))
      IF (T .LT. EPS*TWOPSI/SIX) GO TO 390
         V(PREDUC) = V(PREDUC) + T
         DST = RAD
         T1 = -SI
 390  DO 400 I = 1, P
         J = Q0 + I
         W(J) = T1*W(I) - W(J)
         STEP(I) = W(J) / D(I)
 400     CONTINUE
      V(GTSTEP) =  D7TPR(P, DIG, W(Q))
C
C  ***  SAVE VALUES FOR USE IN A POSSIBLE RESTART  ***
C
 410  V(DSTNRM) = DST
      V(STPPAR) = ALPHAK
      W(LK0) = LK
      W(UK0) = UK
      V(RAD0) = RAD
      W(DSTSAV) = DST
C
C     ***  RESTORE DIAGONAL OF DIHDI  ***
C
      J = 0
      DO 420 I = 1, P
         J = J + I
         K = DIAG0 + I
         DIHDI(J) = W(K)
 420     CONTINUE
C
 999  RETURN
C
C  ***  LAST LINE OF  G7QTS FOLLOWS  ***
      END
      SUBROUTINE  ITSUM(D, G, IV, LIV, LV, P, V, X)
C
C  ***  PRINT ITERATION SUMMARY FOR ***SOL (VERSION 2.3)  ***
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER LIV, LV, P
      INTEGER IV(LIV)
      REAL D(P), G(P), V(LV), X(P)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER ALG, I, IV1, M, NF, NG, OL, PU
      CHARACTER*4 MODEL1(6), MODEL2(6)
      REAL NRELDF, OLDF, PRELDF, RELDF, ZERO
C
C  ***  NO EXTERNAL FUNCTIONS OR SUBROUTINES  ***
C
C  ***  SUBSCRIPTS FOR IV AND V  ***
C
      INTEGER ALGSAV, DSTNRM, F, FDIF, F0, NEEDHD, NFCALL, NFCOV, NGCOV,
     1        NGCALL, NITER, NREDUC, OUTLEV, PREDUC, PRNTIT, PRUNIT,
     2        RELDX, SOLPRT, STATPR, STPPAR, SUSED, X0PRT
C
C  ***  IV SUBSCRIPT VALUES  ***
C
      PARAMETER (ALGSAV=51, NEEDHD=36, NFCALL=6, NFCOV=52, NGCALL=30,
     1           NGCOV=53, NITER=31, OUTLEV=19, PRNTIT=39, PRUNIT=21,
     2           SOLPRT=22, STATPR=23, SUSED=64, X0PRT=24)
C
C  ***  V SUBSCRIPT VALUES  ***
C
      PARAMETER (DSTNRM=2, F=10, F0=13, FDIF=11, NREDUC=6, PREDUC=7,
     1           RELDX=17, STPPAR=5)
C
      PARAMETER (ZERO=0.E+0)
      DATA MODEL1/'    ','    ','    ','    ','  G ','  S '/,
     1     MODEL2/' G  ',' S  ','G-S ','S-G ','-S-G','-G-S'/
C
C-------------------------------  BODY  --------------------------------
C
      PU = IV(PRUNIT)
      IF (PU .EQ. 0) GO TO 999
      IV1 = IV(1)
      IF (IV1 .GT. 62) IV1 = IV1 - 51
      OL = IV(OUTLEV)
      ALG = MOD(IV(ALGSAV)-1,2) + 1
      IF (IV1 .LT. 2 .OR. IV1 .GT. 15) GO TO 370
      IF (IV1 .GE. 12) GO TO 120
      IF (IV1 .EQ. 2 .AND. IV(NITER) .EQ. 0) GO TO 390
      IF (OL .EQ. 0) GO TO 120
      IF (IV1 .GE. 10 .AND. IV(PRNTIT) .EQ. 0) GO TO 120
      IF (IV1 .GT. 2) GO TO 10
         IV(PRNTIT) = IV(PRNTIT) + 1
         IF (IV(PRNTIT) .LT. IABS(OL)) GO TO 999
 10   NF = IV(NFCALL) - IABS(IV(NFCOV))
      IV(PRNTIT) = 0
      RELDF = ZERO
      PRELDF = ZERO
      OLDF =   MAX( ABS(V(F0)),  ABS(V(F)))
      IF (OLDF .LE. ZERO) GO TO 20
         RELDF = V(FDIF) / OLDF
         PRELDF = V(PREDUC) / OLDF
 20   IF (OL .GT. 0) GO TO 60
C
C        ***  PRINT SHORT SUMMARY LINE  ***
C
         IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 1) WRITE(PU,30)
 30   FORMAT(/10H   IT   NF,6X,1HF,7X,5HRELDF,3X,6HPRELDF,3X,5HRELDX,
     1       2X,13HMODEL  STPPAR)
         IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 2) WRITE(PU,40)
 40   FORMAT(/11H    IT   NF,7X,1HF,8X,5HRELDF,4X,6HPRELDF,4X,5HRELDX,
     1       3X,6HSTPPAR)
         IV(NEEDHD) = 0
         IF (ALG .EQ. 2) GO TO 50
         M = IV(SUSED)
         WRITE(PU,100) IV(NITER), NF, V(F), RELDF, PRELDF, V(RELDX),
     1                 MODEL1(M), MODEL2(M), V(STPPAR)
         GO TO 120
C
 50      WRITE(PU,110) IV(NITER), NF, V(F), RELDF, PRELDF, V(RELDX),
     1                 V(STPPAR)
         GO TO 120
C
C     ***  PRINT LONG SUMMARY LINE  ***
C
 60   IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 1) WRITE(PU,70)
 70   FORMAT(/11H    IT   NF,6X,1HF,7X,5HRELDF,3X,6HPRELDF,3X,5HRELDX,
     1       2X,13HMODEL  STPPAR,2X,6HD*STEP,2X,7HNPRELDF)
      IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 2) WRITE(PU,80)
 80   FORMAT(/11H    IT   NF,7X,1HF,8X,5HRELDF,4X,6HPRELDF,4X,5HRELDX,
     1       3X,6HSTPPAR,3X,6HD*STEP,3X,7HNPRELDF)
      IV(NEEDHD) = 0
      NRELDF = ZERO
      IF (OLDF .GT. ZERO) NRELDF = V(NREDUC) / OLDF
      IF (ALG .EQ. 2) GO TO 90
      M = IV(SUSED)
      WRITE(PU,100) IV(NITER), NF, V(F), RELDF, PRELDF, V(RELDX),
     1             MODEL1(M), MODEL2(M), V(STPPAR), V(DSTNRM), NRELDF
      GO TO 120
C
 90   WRITE(PU,110) IV(NITER), NF, V(F), RELDF, PRELDF,
     1             V(RELDX), V(STPPAR), V(DSTNRM), NRELDF
 100  FORMAT(I6,I5,E10.3,2E9.2,E8.1,A3,A4,2E8.1,E9.2)
 110  FORMAT(I6,I5,E11.3,2E10.2,3E9.1,E10.2)
C
 120  IF (IV1 .LE. 2) GO TO 999
      I = IV(STATPR)
      IF (I .EQ. (-1)) GO TO 460
      IF (I + IV1 .LT. 0) GO TO 460
      GO TO (999, 999, 130, 150, 170, 190, 210, 230, 250, 270, 290, 310,
     1       330, 350, 500),  IV1
C
 130  WRITE(PU,140)
 140  FORMAT(/26H ***** X-CONVERGENCE *****)
      GO TO 430
C
 150  WRITE(PU,160)
 160  FORMAT(/42H ***** RELATIVE FUNCTION CONVERGENCE *****)
      GO TO 430
C
 170  WRITE(PU,180)
 180  FORMAT(/49H ***** X- AND RELATIVE FUNCTION CONVERGENCE *****)
      GO TO 430
C
 190  WRITE(PU,200)
 200  FORMAT(/42H ***** ABSOLUTE FUNCTION CONVERGENCE *****)
      GO TO 430
C
 210  WRITE(PU,220)
 220  FORMAT(/33H ***** SINGULAR CONVERGENCE *****)
      GO TO 430
C
 230  WRITE(PU,240)
 240  FORMAT(/30H ***** FALSE CONVERGENCE *****)
      GO TO 430
C
 250  WRITE(PU,260)
 260  FORMAT(/38H ***** FUNCTION EVALUATION LIMIT *****)
      GO TO 430
C
 270  WRITE(PU,280)
 280  FORMAT(/28H ***** ITERATION LIMIT *****)
      GO TO 430
C
 290  WRITE(PU,300)
 300  FORMAT(/18H ***** STOPX *****)
      GO TO 430
C
 310  WRITE(PU,320)
 320  FORMAT(/44H ***** INITIAL F(X) CANNOT BE COMPUTED *****)
C
      GO TO 390
C
 330  WRITE(PU,340)
 340  FORMAT(/37H ***** BAD PARAMETERS TO ASSESS *****)
      GO TO 999
C
 350  WRITE(PU,360)
 360  FORMAT(/43H ***** GRADIENT COULD NOT BE COMPUTED *****)
      IF (IV(NITER) .GT. 0) GO TO 460
      GO TO 390
C
 370  WRITE(PU,380) IV(1)
 380  FORMAT(/14H ***** IV(1) =,I5,6H *****)
      GO TO 999
C
C  ***  INITIAL CALL ON  ITSUM  ***
C
 390  IF (IV(X0PRT) .NE. 0) WRITE(PU,400) (I, X(I), D(I), I = 1, P)
 400  FORMAT(/23H     I     INITIAL X(I),8X,4HD(I)//(1X,I5,E17.6,E14.3))
C     *** THE FOLLOWING ARE TO AVOID UNDEFINED VARIABLES WHEN THE
C     *** FUNCTION EVALUATION LIMIT IS 1...
      V(DSTNRM) = ZERO
      V(FDIF) = ZERO
      V(NREDUC) = ZERO
      V(PREDUC) = ZERO
      V(RELDX) = ZERO
      IF (IV1 .GE. 12) GO TO 999
      IV(NEEDHD) = 0
      IV(PRNTIT) = 0
      IF (OL .EQ. 0) GO TO 999
      IF (OL .LT. 0 .AND. ALG .EQ. 1) WRITE(PU,30)
      IF (OL .LT. 0 .AND. ALG .EQ. 2) WRITE(PU,40)
      IF (OL .GT. 0 .AND. ALG .EQ. 1) WRITE(PU,70)
      IF (OL .GT. 0 .AND. ALG .EQ. 2) WRITE(PU,80)
      IF (ALG .EQ. 1) WRITE(PU,410) IV(NFCALL), V(F)
      IF (ALG .EQ. 2) WRITE(PU,420) IV(NFCALL), V(F)
 410  FORMAT(/6H     0,I5,E10.3)
 420  FORMAT(/6H     0,I5,E11.3)
      GO TO 999
C
C  ***  PRINT VARIOUS INFORMATION REQUESTED ON SOLUTION  ***
C
 430  IV(NEEDHD) = 1
      IF (IV(STATPR) .LE. 0) GO TO 460
         OLDF =   MAX( ABS(V(F0)),  ABS(V(F)))
         PRELDF = ZERO
         NRELDF = ZERO
         IF (OLDF .LE. ZERO) GO TO 440
              PRELDF = V(PREDUC) / OLDF
              NRELDF = V(NREDUC) / OLDF
 440     NF = IV(NFCALL) - IV(NFCOV)
         NG = IV(NGCALL) - IV(NGCOV)
         WRITE(PU,450) V(F), V(RELDX), NF, NG, PRELDF, NRELDF
 450  FORMAT(/9H FUNCTION,E17.6,8H   RELDX,E17.3/12H FUNC. EVALS,
     1   I8,9X,11HGRAD. EVALS,I8/7H PRELDF,E16.3,6X,7HNPRELDF,E15.3)
C
 460  IF (IV(SOLPRT) .EQ. 0) GO TO 999
         IV(NEEDHD) = 1
         IF (IV(ALGSAV) .GT. 2) GO TO 999
         WRITE(PU,470)
 470  FORMAT(/22H     I      FINAL X(I),8X,4HD(I),10X,4HG(I)/)
         DO 480 I = 1, P
 480          WRITE(PU,490) I, X(I), D(I), G(I)
 490     FORMAT(1X,I5,E16.6,2E14.3)
      GO TO 999
C
 500  WRITE(PU,510)
 510  FORMAT(/24H INCONSISTENT DIMENSIONS)
 999  RETURN
C  ***  LAST LINE OF  ITSUM FOLLOWS  ***
      END
      SUBROUTINE  IVSET(ALG, IV, LIV, LV, V)
C
C  ***  SUPPLY ***SOL (VERSION 2.3) DEFAULT VALUES TO IV AND V  ***
C
C  ***  ALG = 1 MEANS REGRESSION CONSTANTS.
C  ***  ALG = 2 MEANS GENERAL UNCONSTRAINED OPTIMIZATION CONSTANTS.
C
      INTEGER LIV, LV
      INTEGER ALG, IV(LIV)
      REAL V(LV)
C
      INTEGER I7MDCN
      EXTERNAL I7MDCN, V7DFL
C I7MDCN... RETURNS MACHINE-DEPENDENT INTEGER CONSTANTS.
C  V7DFL.... PROVIDES DEFAULT VALUES TO V.
C
      INTEGER ALG1, MIV, MV
      INTEGER MINIV(4), MINV(4)
C
C  ***  SUBSCRIPTS FOR IV  ***
C
      INTEGER ALGSAV, COVPRT, COVREQ, DRADPR, DTYPE, HC, IERR, INITH,
     1        INITS, IPIVOT, IVNEED, LASTIV, LASTV, LMAT, MXFCAL,
     2        MXITER, NFCOV, NGCOV, NVDFLT, NVSAVE, OUTLEV, PARPRT,
     3        PARSAV, PERM, PRUNIT, QRTYP, RDREQ, RMAT, SOLPRT, STATPR,
     4        VNEED, VSAVE, X0PRT
C
C  ***  IV SUBSCRIPT VALUES  ***
C
      PARAMETER (ALGSAV=51, COVPRT=14, COVREQ=15, DRADPR=101, DTYPE=16,
     1           HC=71, IERR=75, INITH=25, INITS=25, IPIVOT=76,
     2           IVNEED=3, LASTIV=44, LASTV=45, LMAT=42, MXFCAL=17,
     3           MXITER=18, NFCOV=52, NGCOV=53, NVDFLT=50, NVSAVE=9,
     4           OUTLEV=19, PARPRT=20, PARSAV=49, PERM=58, PRUNIT=21,
     5           QRTYP=80, RDREQ=57, RMAT=78, SOLPRT=22, STATPR=23,
     6           VNEED=4, VSAVE=60, X0PRT=24)
      DATA MINIV(1)/82/, MINIV(2)/59/, MINIV(3)/103/, MINIV(4)/103/,
     1     MINV(1)/98/, MINV(2)/71/, MINV(3)/101/, MINV(4)/85/
C
C-------------------------------  BODY  --------------------------------
C
      IF (PRUNIT .LE. LIV) IV(PRUNIT) = I7MDCN(1)
      IF (ALGSAV .LE. LIV) IV(ALGSAV) = ALG
      IF (ALG .LT. 1 .OR. ALG .GT. 4) GO TO 40
      MIV = MINIV(ALG)
      IF (LIV .LT. MIV) GO TO 20
      MV = MINV(ALG)
      IF (LV .LT. MV) GO TO 30
      ALG1 = MOD(ALG-1,2) + 1
      CALL  V7DFL(ALG1, LV, V)
      IV(1) = 12
      IF (ALG .GT. 2) IV(DRADPR) = 1
      IV(IVNEED) = 0
      IV(LASTIV) = MIV
      IV(LASTV) = MV
      IV(LMAT) = MV + 1
      IV(MXFCAL) = 200
      IV(MXITER) = 150
      IV(OUTLEV) = 1
      IV(PARPRT) = 1
      IV(PERM) = MIV + 1
      IV(SOLPRT) = 1
      IV(STATPR) = 1
      IV(VNEED) = 0
      IV(X0PRT) = 1
C
      IF (ALG1 .GE. 2) GO TO 10
C
C  ***  REGRESSION  VALUES
C
      IV(COVPRT) = 3
      IV(COVREQ) = 1
      IV(DTYPE) = 1
      IV(HC) = 0
      IV(IERR) = 0
      IV(INITS) = 0
      IV(IPIVOT) = 0
      IV(NVDFLT) = 32
      IV(VSAVE) = 58
      IF (ALG .GT. 2) IV(VSAVE) = IV(VSAVE) + 3
      IV(PARSAV) = IV(VSAVE) + NVSAVE
      IV(QRTYP) = 1
      IV(RDREQ) = 3
      IV(RMAT) = 0
      GO TO 999
C
C  ***  GENERAL OPTIMIZATION VALUES
C
 10   IV(DTYPE) = 0
      IV(INITH) = 1
      IV(NFCOV) = 0
      IV(NGCOV) = 0
      IV(NVDFLT) = 25
      IV(PARSAV) = 47
      IF (ALG .GT. 2) IV(PARSAV) = 61
      GO TO 999
C
 20   IV(1) = 15
      GO TO 999
C
 30   IV(1) = 16
      GO TO 999
C
 40   IV(1) = 67
C
 999  RETURN
C  ***  LAST LINE OF  IVSET FOLLOWS  ***
      END
      SUBROUTINE  L7ITV(N, X, L, Y)
C
C  ***  SOLVE  (L**T)*X = Y,  WHERE  L  IS AN  N X N  LOWER TRIANGULAR
C  ***  MATRIX STORED COMPACTLY BY ROWS.  X AND Y MAY OCCUPY THE SAME
C  ***  STORAGE.  ***
C
      INTEGER N
      REAL X(N), L(1), Y(N)
      INTEGER I, II, IJ, IM1, I0, J, NP1
      REAL XI, ZERO
      PARAMETER (ZERO=0.E+0)
C
      DO 10 I = 1, N
 10      X(I) = Y(I)
      NP1 = N + 1
      I0 = N*(N+1)/2
      DO 30 II = 1, N
         I = NP1 - II
         XI = X(I)/L(I0)
         X(I) = XI
         IF (I .LE. 1) GO TO 999
         I0 = I0 - I
         IF (XI .EQ. ZERO) GO TO 30
         IM1 = I - 1
         DO 20 J = 1, IM1
              IJ = I0 + J
              X(J) = X(J) - XI*L(IJ)
 20           CONTINUE
 30      CONTINUE
 999  RETURN
C  ***  LAST LINE OF  L7ITV FOLLOWS  ***
      END
      SUBROUTINE  L7IVM(N, X, L, Y)
C
C  ***  SOLVE  L*X = Y, WHERE  L  IS AN  N X N  LOWER TRIANGULAR
C  ***  MATRIX STORED COMPACTLY BY ROWS.  X AND Y MAY OCCUPY THE SAME
C  ***  STORAGE.  ***
C
      INTEGER N
      REAL X(N), L(1), Y(N)
      REAL  D7TPR
      EXTERNAL  D7TPR
      INTEGER I, J, K
      REAL T, ZERO
      PARAMETER (ZERO=0.E+0)
C
      DO 10 K = 1, N
         IF (Y(K) .NE. ZERO) GO TO 20
         X(K) = ZERO
 10      CONTINUE
      GO TO 999
 20   J = K*(K+1)/2
      X(K) = Y(K) / L(J)
      IF (K .GE. N) GO TO 999
      K = K + 1
      DO 30 I = K, N
         T =  D7TPR(I-1, L(J+1), X)
         J = J + I
         X(I) = (Y(I) - T)/L(J)
 30      CONTINUE
 999  RETURN
C  ***  LAST LINE OF  L7IVM FOLLOWS  ***
      END
      SUBROUTINE  L7MST(D, G, IERR, IPIVOT, KA, P, QTR, R, STEP, V, W)
C
C  ***  COMPUTE LEVENBERG-MARQUARDT STEP USING MORE-HEBDEN TECHNIQUE  **
C  ***  NL2SOL VERSION 2.2.  ***
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER IERR, KA, P
      INTEGER IPIVOT(P)
      REAL D(P), G(P), QTR(P), R(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  Q7RGS -- 0 MEANS R HAS
C             FULL RANK.
C IPIVOT (I/O) = PERMUTATION ARRAY FROM QRFACT OR  Q7RGS, 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              L7MST 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  D7TPR - RETURNS INNER PRODUCT OF TWO VECTORS.
C  L7ITV - APPLY INVERSE-TRANSPOSE OF COMPACT LOWER TRIANG. MATRIX.
C  L7IVM - APPLY INVERSE OF COMPACT LOWER TRIANG. MATRIX.
C  V7CPY  - COPIES ONE VECTOR TO ANOTHER.
C  V2NRM - RETURNS 2-NORM OF A VECTOR.
C
C  ***  REFERENCES  ***
C
C 1.  DENNIS, J.E., GAY, D.M., AND WELSCH, R.E. (1981), AN ADAPTIVE
C             NONLINEAR LEAST-SQUARES ALGORITHM, ACM TRANS. MATH.
C             SOFTWARE, VOL. 7, NO. 3.
C 2.  GAY, D.M. (1981), COMPUTING OPTIMAL LOCALLY CONSTRAINED STEPS,
C             SIAM J. SCI. STATIST. COMPUTING, VOL. 2, NO. 2, PP.
C             186-197.
C 3.  LAWSON, C.L., AND HANSON, R.J. (1974), SOLVING LEAST SQUARES
C             PROBLEMS, PRENTICE-HALL, ENGLEWOOD CLIFFS, N.J.
C 4.  MORE, J.J. (1978), THE LEVENBERG-MARQUARDT ALGORITHM, IMPLEMEN-
C             TATION AND THEORY, PP.105-116 OF SPRINGER LECTURE NOTES
C             IN MATHEMATICS NO. 630, EDITED BY G.A. WATSON, SPRINGER-
C             VERLAG, BERLIN AND NEW YORK.
C
C  ***  GENERAL  ***
C
C     CODED BY DAVID M. GAY.
C     THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH
C     SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS
C     MCS-7600324, DCR75-10143, 76-14311DSS, MCS76-11989, AND
C     MCS-7906671.
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER DSTSAV, I, IP1, I1, J1, K, KALIM, L, LK0, PHIPIN,
     1        PP1O2, RES, RES0, RMAT, RMAT0, UK0
      REAL A, ADI, ALPHAK, B, DFACSQ, DST, DTOL, D1, D2,
     1                 LK, OLDPHI, PHI, PHIMAX, PHIMIN, PSIFAC, RAD,
     2                 SI, SJ, SQRTAK, T, TWOPSI, UK, WL
C
C     ***  CONSTANTS  ***
      REAL DFAC, EIGHT, HALF, NEGONE, ONE, P001, THREE,
     1                 TTOL, ZERO
      REAL BIG
C
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      REAL  D7TPR,  L7SVN,  R7MDC,  V2NRM
      EXTERNAL  D7TPR,  L7ITV,  L7IVM,  L7SVN,  R7MDC, V7CPY,  V2NRM
C
C  ***  SUBSCRIPTS FOR V  ***
C
      INTEGER DGNORM, DSTNRM, DST0, EPSLON, GTSTEP, NREDUC, PHMNFC,
     1        PHMXFC, PREDUC, RADIUS, RAD0, STPPAR
      PARAMETER (DGNORM=1, DSTNRM=2, DST0=3, EPSLON=19, GTSTEP=4,
     1           NREDUC=6, PHMNFC=20, PHMXFC=21, PREDUC=7, RADIUS=8,
     2           RAD0=9, STPPAR=5)
C
      PARAMETER (DFAC=256.E+0, EIGHT=8.E+0, HALF=0.5E+0, NEGONE=-1.E+0,
     1     ONE=1.E+0, P001=1.E-3, THREE=3.E+0, TTOL=2.5E+0,
     2     ZERO=0.E+0)
      SAVE BIG
      DATA BIG/0.E+0/
C
C  ***  BODY  ***
C
C     ***  FOR USE IN RECOMPUTING STEP, THE FINAL VALUES OF LK AND UK,
C     ***  THE INVERSE DERIVATIVE OF MORE*S PHI AT 0 (FOR NONSING. J)
C     ***  AND THE VALUE RETURNED AS V(DSTNRM) ARE STORED AT W(LK0),
C     ***  W(UK0), W(PHIPIN), AND W(DSTSAV) RESPECTIVELY.
      LK0 = P + 1
      PHIPIN = LK0 + 1
      UK0 = PHIPIN + 1
      DSTSAV = UK0 + 1
      RMAT0 = DSTSAV
C     ***  A COPY OF THE R-MATRIX FROM THE QR DECOMPOSITION OF J IS
C     ***  STORED IN W STARTING AT W(RMAT), AND A COPY OF THE RESIDUAL
C     ***  VECTOR IS STORED IN W STARTING AT W(RES).  THE LOOPS BELOW
C     ***  THAT UPDATE THE QR DECOMP. FOR A NONZERO MARQUARDT PARAMETER
C     ***  WORK ON THESE COPIES.
      RMAT = RMAT0 + 1
      PP1O2 = P * (P + 1) / 2
      RES0 = PP1O2 + RMAT0
      RES = RES0 + 1
      RAD = V(RADIUS)
      IF (RAD .GT. ZERO)
     1   PSIFAC = V(EPSLON)/((EIGHT*(V(PHMNFC) + ONE) + THREE) * RAD**2)
      IF (BIG .LE. ZERO) BIG =  R7MDC(6)
      PHIMAX = V(PHMXFC) * RAD
      PHIMIN = V(PHMNFC) * RAD
C     ***  DTOL, DFAC, AND DFACSQ ARE USED IN RESCALING THE FAST GIVENS
C     ***  REPRESENTATION OF THE UPDATED QR DECOMPOSITION.
      DTOL = ONE/DFAC
      DFACSQ = DFAC*DFAC
C     ***  OLDPHI IS USED TO DETECT LIMITS OF NUMERICAL ACCURACY.  IF
C     ***  WE RECOMPUTE STEP AND IT DOES NOT CHANGE, THEN WE ACCEPT IT.
      OLDPHI = ZERO
      LK = ZERO
      UK = ZERO
      KALIM = KA + 12
C
C  ***  START OR RESTART, DEPENDING ON KA  ***
C
      IF (KA) 10, 20, 370
C
C  ***  FRESH START -- COMPUTE V(NREDUC)  ***
C
 10   KA = 0
      KALIM = 12
      K = P
      IF (IERR .NE. 0) K = IABS(IERR) - 1
      V(NREDUC) = HALF* D7TPR(K, QTR, QTR)
C
C  ***  SET UP TO TRY INITIAL GAUSS-NEWTON STEP  ***
C
 20   V(DST0) = NEGONE
      IF (IERR .NE. 0) GO TO 90
      T =  L7SVN(P, R, STEP, W(RES))
      IF (T .GE. ONE) GO TO 30
         IF ( V2NRM(P, QTR) .GE. BIG*T) GO TO 90
C
C  ***  COMPUTE GAUSS-NEWTON STEP  ***
C
C     ***  NOTE -- THE R-MATRIX IS STORED COMPACTLY BY COLUMNS IN
C     ***  R(1), R(2), R(3), ...  IT IS THE TRANSPOSE OF A
C     ***  LOWER TRIANGULAR MATRIX STORED COMPACTLY BY ROWS, AND WE
C     ***  TREAT IT AS SUCH WHEN USING  L7ITV AND  L7IVM.
 30   CALL  L7ITV(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 =  V2NRM(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  L7IVM(P, STEP, R, STEP)
      T = ONE /  V2NRM(P, STEP)
      W(PHIPIN) = (T/RAD)*T
      LK = PHI*W(PHIPIN)
C
C  ***  COMPUTE U0  ***
C
 90   DO 100 I = 1, P
 100     W(I) = G(I)/D(I)
      V(DGNORM) =  V2NRM(P, W)
      UK = V(DGNORM)/RAD
      IF (UK .LE. ZERO) GO TO 390
C
C     ***  ALPHAK WILL BE USED AS THE CURRENT MARQUARDT PARAMETER.  WE
C     ***  USE MORE*S SCHEME FOR INITIALIZING IT.
C
      ALPHAK =  ABS(V(STPPAR)) * V(RAD0)/RAD
      ALPHAK =   MIN(UK,   MAX(ALPHAK, LK))
C
C
C  ***  TOP OF LOOP -- INCREMENT KA, COPY R TO RMAT, QTR TO RES  ***
C
 110  KA = KA + 1
      CALL  V7CPY(PP1O2, W(RMAT), R)
      CALL  V7CPY(P, W(RES), QTR)
C
C  ***  SAFEGUARD ALPHAK AND INITIALIZE FAST GIVENS SCALE VECTOR.  ***
C
      IF (ALPHAK .LE. ZERO .OR. ALPHAK .LT. LK .OR. ALPHAK .GE. UK)
     1             ALPHAK = UK *   MAX(P001,  SQRT(LK/UK))
      IF (ALPHAK .LE. ZERO) ALPHAK = HALF * UK
      SQRTAK =  SQRT(ALPHAK)
      DO 120 I = 1, P
 120     W(I) = ONE
C
C  ***  ADD ALPHAK*D AND UPDATE QR DECOMP. USING FAST GIVENS TRANS.  ***
C
      DO 270 I = 1, P
C        ***  GENERATE, APPLY 1ST GIVENS TRANS. FOR ROW I OF ALPHAK*D.
C        ***  (USE STEP TO STORE TEMPORARY ROW)  ***
         L = I*(I+1)/2 + RMAT0
         WL = W(L)
         D2 = ONE
         D1 = W(I)
         J1 = IPIVOT(I)
         ADI = SQRTAK*D(J1)
         IF (ADI .GE.  ABS(WL)) GO TO 150
 130     A = ADI/WL
         B = D2*A/D1
         T = A*B + ONE
         IF (T .GT. TTOL) GO TO 150
         W(I) = D1/T
         D2 = D2/T
         W(L) = T*WL
         A = -A
         DO 140 J1 = I, P
              L = L + J1
              STEP(J1) = A*W(L)
 140          CONTINUE
         GO TO 170
C
 150     B = WL/ADI
         A = D1*B/D2
         T = A*B + ONE
         IF (T .GT. TTOL) GO TO 130
         W(I) = D2/T
         D2 = D1/T
         W(L) = T*ADI
         DO 160 J1 = I, P
              L = L + J1
              WL = W(L)
              STEP(J1) = -WL
              W(L) = A*WL
 160          CONTINUE
C
 170     IF (I .EQ. P) GO TO 280
C
C        ***  NOW USE GIVENS TRANS. TO ZERO ELEMENTS OF TEMP. ROW  ***
C
         IP1 = I + 1
         DO 260 I1 = IP1, P
              L = I1*(I1+1)/2 + RMAT0
              WL = W(L)
              SI = STEP(I1-1)
              D1 = W(I1)
C
C             ***  RESCALE ROW I1 IF NECESSARY  ***
C
              IF (D1 .GE. DTOL) GO TO 190
                   D1 = D1*DFACSQ
                   WL = WL/DFAC
                   K = L
                   DO 180 J1 = I1, P
                        K = K + J1
                        W(K) = W(K)/DFAC
 180                    CONTINUE
C
C             ***  USE GIVENS TRANS. TO ZERO NEXT ELEMENT OF TEMP. ROW
C
 190          IF ( ABS(SI) .GT.  ABS(WL)) GO TO 220
              IF (SI .EQ. ZERO) GO TO 260
 200          A = SI/WL
              B = D2*A/D1
              T = A*B + ONE
              IF (T .GT. TTOL) GO TO 220
              W(L) = T*WL
              W(I1) = D1/T
              D2 = D2/T
              DO 210 J1 = I1, P
                   L = L + J1
                   WL = W(L)
                   SJ = STEP(J1)
                   W(L) = WL + B*SJ
                   STEP(J1) = SJ - A*WL
 210               CONTINUE
              GO TO 240
C
 220          B = WL/SI
              A = D1*B/D2
              T = A*B + ONE
              IF (T .GT. TTOL) GO TO 200
              W(I1) = D2/T
              D2 = D1/T
              W(L) = T*SI
              DO 230 J1 = I1, P
                   L = L + J1
                   WL = W(L)
                   SJ = STEP(J1)
                   W(L) = A*WL + SJ
                   STEP(J1) = B*SJ - WL
 230               CONTINUE
C
C             ***  RESCALE TEMP. ROW IF NECESSARY  ***
C
 240          IF (D2 .GE. DTOL) GO TO 260
                   D2 = D2*DFACSQ
                   DO 250 K = I1, P
 250                    STEP(K) = STEP(K)/DFAC
 260          CONTINUE
 270     CONTINUE
C
C  ***  COMPUTE STEP  ***
C
 280  CALL  L7ITV(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 =  V2NRM(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 -  D7TPR(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  L7IVM(P, STEP, W(RMAT), STEP)
      DO 340 I = 1, P
 340     STEP(I) = STEP(I) /  SQRT(W(I))
      T = ONE /  V2NRM(P, STEP)
      ALPHAK = ALPHAK + T*PHI*T/RAD
      LK =   MAX(LK, ALPHAK)
      ALPHAK = LK
      GO TO 110
C
C  ***  RESTART  ***
C
 370  LK = W(LK0)
      UK = W(UK0)
      IF (V(DST0) .GT. ZERO .AND. V(DST0) - RAD .LE. PHIMAX) GO TO 20
      ALPHAK =  ABS(V(STPPAR))
      DST = W(DSTSAV)
      PHI = DST - RAD
      T = V(DGNORM)/RAD
      IF (RAD .GT. V(RAD0)) GO TO 380
C
C        ***  SMALLER RADIUS  ***
         UK = T
         IF (ALPHAK .LE. ZERO) LK = ZERO
         IF (V(DST0) .GT. ZERO) LK =   MAX(LK, (V(DST0)-RAD)*W(PHIPIN))
         GO TO 300
C
C     ***  BIGGER RADIUS  ***
 380  IF (ALPHAK .LE. ZERO .OR. UK .GT. T) UK = T
      LK = ZERO
      IF (V(DST0) .GT. ZERO) LK =   MAX(LK, (V(DST0)-RAD)*W(PHIPIN))
      GO TO 300
C
C  ***  SPECIAL CASE -- RAD .LE. 0 OR (G = 0 AND J IS SINGULAR)  ***
C
 390  V(STPPAR) = ZERO
      DST = ZERO
      LK = ZERO
      UK = ZERO
      V(GTSTEP) = ZERO
      V(PREDUC) = ZERO
      DO 400 I = 1, P
 400     STEP(I) = ZERO
      GO TO 450
C
C  ***  ACCEPTABLE GAUSS-NEWTON STEP -- RECOVER STEP FROM W  ***
C
 410  ALPHAK = ZERO
      DO 420 I = 1, P
         J1 = IPIVOT(I)
         STEP(J1) = -W(I)
 420     CONTINUE
C
C  ***  SAVE VALUES FOR USE IN A POSSIBLE RESTART  ***
C
 430  V(STPPAR) = ALPHAK
 440  V(GTSTEP) =   MIN( D7TPR(P,STEP,G), ZERO)
      V(PREDUC) = HALF * (ALPHAK*DST*DST - V(GTSTEP))
 450  V(DSTNRM) = DST
      W(DSTSAV) = DST
      W(LK0) = LK
      W(UK0) = UK
      V(RAD0) = RAD
C
 999  RETURN
C
C  ***  LAST LINE OF  L7MST FOLLOWS  ***
      END
      SUBROUTINE  L7SQR(N, A, L)
C
C  ***  COMPUTE  A = LOWER TRIANGLE OF  L*(L**T),  WITH BOTH
C  ***  L  AND  A  STORED COMPACTLY BY ROWS.  (BOTH MAY OCCUPY THE
C  ***  SAME STORAGE.
C
C  ***  PARAMETERS  ***
C
      INTEGER N
      REAL A(1), L(1)
C     DIMENSION A(N*(N+1)/2), L(N*(N+1)/2)
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER I, II, IJ, IK, IP1, I0, J, JJ, JK, J0, K, NP1
      REAL T
C
      NP1 = N + 1
      I0 = N*(N+1)/2
      DO 30 II = 1, N
         I = NP1 - II
         IP1 = I + 1
         I0 = I0 - I
         J0 = I*(I+1)/2
         DO 20 JJ = 1, I
              J = IP1 - JJ
              J0 = J0 - J
              T = 0.0E0
              DO 10 K = 1, J
                   IK = I0 + K
                   JK = J0 + K
                   T = T + L(IK)*L(JK)
 10                CONTINUE
              IJ = I0 + J
              A(IJ) = T
 20           CONTINUE
 30      CONTINUE
 999  RETURN
      END
      SUBROUTINE  L7SRT(N1, N, L, A, IRC)
C
C  ***  COMPUTE ROWS N1 THROUGH N OF THE CHOLESKY FACTOR  L  OF
C  ***  A = L*(L**T),  WHERE  L  AND THE LOWER TRIANGLE OF  A  ARE BOTH
C  ***  STORED COMPACTLY BY ROWS (AND MAY OCCUPY THE SAME STORAGE).
C  ***  IRC = 0 MEANS ALL WENT WELL.  IRC = J MEANS THE LEADING
C  ***  PRINCIPAL  J X J  SUBMATRIX OF  A  IS NOT POSITIVE DEFINITE --
C  ***  AND  L(J*(J+1)/2)  CONTAINS THE (NONPOS.) REDUCED J-TH DIAGONAL.
C
C  ***  PARAMETERS  ***
C
      INTEGER N1, N, IRC
      REAL L(1), A(1)
C     DIMENSION L(N*(N+1)/2), A(N*(N+1)/2)
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER I, IJ, IK, IM1, I0, J, JK, JM1, J0, K
      REAL T, TD, ZERO
C
      PARAMETER (ZERO=0.E+0)
C
C  ***  BODY  ***
C
      I0 = N1 * (N1 - 1) / 2
      DO 50 I = N1, N
         TD = ZERO
         IF (I .EQ. 1) GO TO 40
         J0 = 0
         IM1 = I - 1
         DO 30 J = 1, IM1
              T = ZERO
              IF (J .EQ. 1) GO TO 20
              JM1 = J - 1
              DO 10 K = 1, JM1
                   IK = I0 + K
                   JK = J0 + K
                   T = T + L(IK)*L(JK)
 10                CONTINUE
 20           IJ = I0 + J
              J0 = J0 + J
              T = (A(IJ) - T) / L(J0)
              L(IJ) = T
              TD = TD + T*T
 30           CONTINUE
 40      I0 = I0 + I
         T = A(I0) - TD
         IF (T .LE. ZERO) GO TO 60
         L(I0) =  SQRT(T)
 50      CONTINUE
C
      IRC = 0
      GO TO 999
C
 60   L(I0) = T
      IRC = I
C
 999  RETURN
C
C  ***  LAST LINE OF  L7SRT  ***
      END
      REAL FUNCTION  L7SVN(P, L, X, Y)
C
C  ***  ESTIMATE SMALLEST SING. VALUE OF PACKED LOWER TRIANG. MATRIX L
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER P
      REAL L(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  L7SVN 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  L7SVN RETURNS ZERO, THEN SOME COMPONENTS OF X
C             ARE ZERO AND THE REST RETAIN THEIR INPUT VALUES.
C  Y (OUT) IF  L7SVN 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  L7SVN 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  L7SVN RETURNS).
C
C  ***  ALGORITHM NOTES  ***
C
C     THE ALGORITHM IS BASED ON (1), WITH THE ADDITIONAL PROVISION THAT
C      L7SVN = 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         V2NRM - FUNCTION, RETURNS THE 2-NORM OF A VECTOR.
C
C  ***  REFERENCES  ***
C
C     (1) CLINE, A., MOLER, C., STEWART, G., AND WILKINSON, J.H.(1977),
C         AN ESTIMATE FOR THE CONDITION NUMBER OF A MATRIX, REPORT
C         TM-310, APPLIED MATH. DIV., ARGONNE NATIONAL LABORATORY.
C
C     (2) HOAGLIN, D.C. (1976), THEORETICAL PROPERTIES OF CONGRUENTIAL
C         RANDOM-NUMBER GENERATORS --  AN EMPIRICAL VIEW,
C         MEMORANDUM NS-340, DEPT. OF STATISTICS, HARVARD UNIV.
C
C     (3) KNUTH, D.E. (1969), THE ART OF COMPUTER PROGRAMMING, VOL. 2
C         (SEMINUMERICAL ALGORITHMS), ADDISON-WESLEY, READING, MASS.
C
C     (4) SMITH, C.S. (1971), MULTIPLICATIVE PSEUDO-RANDOM NUMBER
C         GENERATORS WITH PRIME MODULUS, J. ASSOC. COMPUT. MACH. 18,
C         PP. 586-593.
C
C  ***  HISTORY  ***
C
C     DESIGNED AND CODED BY DAVID M. GAY (WINTER 1977/SUMMER 1978).
C
C  ***  GENERAL  ***
C
C     THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH
C     SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS
C     MCS-7600324, DCR75-10143, 76-14311DSS, AND MCS76-11989.
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER I, II, IX, J, JI, JJ, JJJ, JM1, J0, PM1
      REAL B, SMINUS, SPLUS, T, XMINUS, XPLUS
C
C  ***  CONSTANTS  ***
C
      REAL HALF, ONE, R9973, ZERO
C
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      REAL  D7TPR,  V2NRM
      EXTERNAL  D7TPR,  V2NRM, V2AXY
C
      PARAMETER (HALF=0.5E+0, ONE=1.E+0, R9973=9973.E+0, ZERO=0.E+0)
C
C  ***  BODY  ***
C
      IX = 2
      PM1 = P - 1
C
C  ***  FIRST CHECK WHETHER TO RETURN  L7SVN = 0 AND INITIALIZE X  ***
C
      II = 0
      J0 = P*PM1/2
      JJ = J0 + P
      IF (L(JJ) .EQ. ZERO) GO TO 110
      IX = MOD(3432*IX, 9973)
      B = HALF*(ONE + FLOAT(IX)/R9973)
      XPLUS = B / L(JJ)
      X(P) = XPLUS
      IF (P .LE. 1) GO TO 60
      DO 10 I = 1, PM1
         II = II + I
         IF (L(II) .EQ. ZERO) GO TO 110
         JI = J0 + I
         X(I) = XPLUS * L(JI)
 10      CONTINUE
C
C  ***  SOLVE (L**T)*X = B, WHERE THE COMPONENTS OF B HAVE RANDOMLY
C  ***  CHOSEN MAGNITUDES IN (.5,1) WITH SIGNS CHOSEN TO MAKE X LARGE.
C
C     DO J = P-1 TO 1 BY -1...
      DO 50 JJJ = 1, PM1
         J = P - JJJ
C       ***  DETERMINE X(J) IN THIS ITERATION. NOTE FOR I = 1,2,...,J
C       ***  THAT X(I) HOLDS THE CURRENT PARTIAL SUM FOR ROW I.
         IX = MOD(3432*IX, 9973)
         B = HALF*(ONE + FLOAT(IX)/R9973)
         XPLUS = (B - X(J))
         XMINUS = (-B - X(J))
         SPLUS =  ABS(XPLUS)
         SMINUS =  ABS(XMINUS)
         JM1 = J - 1
         J0 = J*JM1/2
         JJ = J0 + J
         XPLUS = XPLUS/L(JJ)
         XMINUS = XMINUS/L(JJ)
         IF (JM1 .EQ. 0) GO TO 30
         DO 20 I = 1, JM1
              JI = J0 + I
              SPLUS = SPLUS +  ABS(X(I) + L(JI)*XPLUS)
              SMINUS = SMINUS +  ABS(X(I) + L(JI)*XMINUS)
 20           CONTINUE
 30      IF (SMINUS .GT. SPLUS) XPLUS = XMINUS
         X(J) = XPLUS
C       ***  UPDATE PARTIAL SUMS  ***
         IF (JM1 .GT. 0) CALL  V2AXY(JM1, X, XPLUS, L(J0+1), X)
 50      CONTINUE
C
C  ***  NORMALIZE X  ***
C
 60   T = ONE/ V2NRM(P, X)
      DO 70 I = 1, P
 70      X(I) = T*X(I)
C
C  ***  SOLVE L*Y = X AND RETURN  L7SVN = 1/TWONORM(Y)  ***
C
      DO 100 J = 1, P
         JM1 = J - 1
         J0 = J*JM1/2
         JJ = J0 + J
         T = ZERO
         IF (JM1 .GT. 0) T =  D7TPR(JM1, L(J0+1), Y)
         Y(J) = (X(J) - T) / L(JJ)
 100     CONTINUE
C
       L7SVN = ONE/ V2NRM(P, Y)
      GO TO 999
C
 110   L7SVN = ZERO
 999  RETURN
C  ***  LAST LINE OF  L7SVN FOLLOWS  ***
      END
      REAL FUNCTION  L7SVX(P, L, X, Y)
C
C  ***  ESTIMATE LARGEST SING. VALUE OF PACKED LOWER TRIANG. MATRIX L
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER P
      REAL L(1), X(P), Y(P)
C     DIMENSION L(P*(P+1)/2)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ***  PURPOSE  ***
C
C     THIS FUNCTION RETURNS A GOOD UNDER-ESTIMATE OF THE LARGEST
C     SINGULAR VALUE OF THE PACKED LOWER TRIANGULAR MATRIX L.
C
C  ***  PARAMETER DESCRIPTION  ***
C
C  P (IN)  = THE ORDER OF L.  L IS A  P X P  LOWER TRIANGULAR MATRIX.
C  L (IN)  = ARRAY HOLDING THE ELEMENTS OF  L  IN ROW ORDER, I.E.
C             L(1,1), L(2,1), L(2,2), L(3,1), L(3,2), L(3,3), ETC.
C  X (OUT) IF  L7SVX RETURNS A POSITIVE VALUE, THEN X = (L**T)*Y IS AN
C             (UNNORMALIZED) APPROXIMATE RIGHT SINGULAR VECTOR
C             CORRESPONDING TO THE LARGEST SINGULAR VALUE.  THIS
C             APPROXIMATION MAY BE CRUDE.
C  Y (OUT) IF  L7SVX RETURNS A POSITIVE VALUE, THEN Y = L*X IS A
C             NORMALIZED APPROXIMATE LEFT SINGULAR VECTOR CORRESPOND-
C             ING TO THE LARGEST SINGULAR VALUE.  THIS APPROXIMATION
C             MAY BE VERY CRUDE.  THE CALLER MAY PASS THE SAME VECTOR
C             FOR X AND Y (NONSTANDARD FORTRAN USAGE), IN WHICH CASE X
C             OVER-WRITES Y.
C
C  ***  ALGORITHM NOTES  ***
C
C     THE ALGORITHM IS BASED ON ANALOGY WITH (1).  IT USES A
C     RANDOM NUMBER GENERATOR PROPOSED IN (4), WHICH PASSES THE
C     SPECTRAL TEST WITH FLYING COLORS -- SEE (2) AND (3).
C
C  ***  SUBROUTINES AND FUNCTIONS CALLED  ***
C
C         V2NRM - FUNCTION, RETURNS THE 2-NORM OF A VECTOR.
C
C  ***  REFERENCES  ***
C
C     (1) CLINE, A., MOLER, C., STEWART, G., AND WILKINSON, J.H.(1977),
C         AN ESTIMATE FOR THE CONDITION NUMBER OF A MATRIX, REPORT
C         TM-310, APPLIED MATH. DIV., ARGONNE NATIONAL LABORATORY.
C
C     (2) HOAGLIN, D.C. (1976), THEORETICAL PROPERTIES OF CONGRUENTIAL
C         RANDOM-NUMBER GENERATORS --  AN EMPIRICAL VIEW,
C         MEMORANDUM NS-340, DEPT. OF STATISTICS, HARVARD UNIV.
C
C     (3) KNUTH, D.E. (1969), THE ART OF COMPUTER PROGRAMMING, VOL. 2
C         (SEMINUMERICAL ALGORITHMS), ADDISON-WESLEY, READING, MASS.
C
C     (4) SMITH, C.S. (1971), MULTIPLICATIVE PSEUDO-RANDOM NUMBER
C         GENERATORS WITH PRIME MODULUS, J. ASSOC. COMPUT. MACH. 18,
C         PP. 586-593.
C
C  ***  HISTORY  ***
C
C     DESIGNED AND CODED BY DAVID M. GAY (WINTER 1977/SUMMER 1978).
C
C  ***  GENERAL  ***
C
C     THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH
C     SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS
C     MCS-7600324, DCR75-10143, 76-14311DSS, AND MCS76-11989.
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER I, IX, J, JI, JJ, JJJ, JM1, J0, PM1, PPLUS1
      REAL B, BLJI, SMINUS, SPLUS, T, YI
C
C  ***  CONSTANTS  ***
C
      REAL HALF, ONE, R9973, ZERO
C
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      REAL  D7TPR,  V2NRM
      EXTERNAL  D7TPR,  V2NRM, V2AXY
C
      PARAMETER (HALF=0.5E+0, ONE=1.E+0, R9973=9973.E+0, ZERO=0.E+0)
C
C  ***  BODY  ***
C
      IX = 2
      PPLUS1 = P + 1
      PM1 = P - 1
C
C  ***  FIRST INITIALIZE X TO PARTIAL SUMS  ***
C
      J0 = P*PM1/2
      JJ = J0 + P
      IX = MOD(3432*IX, 9973)
      B = HALF*(ONE + FLOAT(IX)/R9973)
      X(P) = B * L(JJ)
      IF (P .LE. 1) GO TO 40
      DO 10 I = 1, PM1
         JI = J0 + I
         X(I) = B * L(JI)
 10      CONTINUE
C
C  ***  COMPUTE X = (L**T)*B, WHERE THE COMPONENTS OF B HAVE RANDOMLY
C  ***  CHOSEN MAGNITUDES IN (.5,1) WITH SIGNS CHOSEN TO MAKE X LARGE.
C
C     DO J = P-1 TO 1 BY -1...
      DO 30 JJJ = 1, PM1
         J = P - JJJ
C       ***  DETERMINE X(J) IN THIS ITERATION. NOTE FOR I = 1,2,...,J
C       ***  THAT X(I) HOLDS THE CURRENT PARTIAL SUM FOR ROW I.
         IX = MOD(3432*IX, 9973)
         B = HALF*(ONE + FLOAT(IX)/R9973)
         JM1 = J - 1
         J0 = J*JM1/2
         SPLUS = ZERO
         SMINUS = ZERO
         DO 20 I = 1, J
              JI = J0 + I
              BLJI = B * L(JI)
              SPLUS = SPLUS +  ABS(BLJI + X(I))
              SMINUS = SMINUS +  ABS(BLJI - X(I))
 20           CONTINUE
         IF (SMINUS .GT. SPLUS) B = -B
         X(J) = ZERO
C        ***  UPDATE PARTIAL SUMS  ***
         CALL  V2AXY(J, X, B, L(J0+1), X)
 30      CONTINUE
C
C  ***  NORMALIZE X  ***
C
 40   T =  V2NRM(P, X)
      IF (T .LE. ZERO) GO TO 80
      T = ONE / T
      DO 50 I = 1, P
 50      X(I) = T*X(I)
C
C  ***  COMPUTE L*X = Y AND RETURN SVMAX = TWONORM(Y)  ***
C
      DO 60 JJJ = 1, P
         J = PPLUS1 - JJJ
         JI = J*(J-1)/2 + 1
         Y(J) =  D7TPR(J, L(JI), X)
 60      CONTINUE
C
C  ***  NORMALIZE Y AND SET X = (L**T)*Y  ***
C
      T = ONE /  V2NRM(P, Y)
      JI = 1
      DO 70 I = 1, P
         YI = T * Y(I)
         X(I) = ZERO
         CALL  V2AXY(I, X, YI, L(JI), X)
         JI = JI + I
 70      CONTINUE
       L7SVX =  V2NRM(P, X)
      GO TO 999
C
 80    L7SVX = ZERO
C
 999  RETURN
C  ***  LAST LINE OF  L7SVX FOLLOWS  ***
      END
      SUBROUTINE  L7TVM(N, X, L, Y)
C
C  ***  COMPUTE  X = (L**T)*Y, WHERE  L  IS AN  N X N  LOWER
C  ***  TRIANGULAR MATRIX STORED COMPACTLY BY ROWS.  X AND Y MAY
C  ***  OCCUPY THE SAME STORAGE.  ***
C
      INTEGER N
      REAL X(N), L(1), Y(N)
C     DIMENSION L(N*(N+1)/2)
      INTEGER I, IJ, I0, J
      REAL YI, ZERO
      PARAMETER (ZERO=0.E+0)
C
      I0 = 0
      DO 20 I = 1, N
         YI = Y(I)
         X(I) = ZERO
         DO 10 J = 1, I
              IJ = I0 + J
              X(J) = X(J) + YI*L(IJ)
 10           CONTINUE
         I0 = I0 + I
 20      CONTINUE
 999  RETURN
C  ***  LAST LINE OF  L7TVM FOLLOWS  ***
      END
      SUBROUTINE  L7VML(N, X, L, Y)
C
C  ***  COMPUTE  X = L*Y, WHERE  L  IS AN  N X N  LOWER TRIANGULAR
C  ***  MATRIX STORED COMPACTLY BY ROWS.  X AND Y MAY OCCUPY THE SAME
C  ***  STORAGE.  ***
C
      INTEGER N
      REAL X(N), L(1), Y(N)
C     DIMENSION L(N*(N+1)/2)
      INTEGER I, II, IJ, I0, J, NP1
      REAL T, ZERO
      PARAMETER (ZERO=0.E+0)
C
      NP1 = N + 1
      I0 = N*(N+1)/2
      DO 20 II = 1, N
         I = NP1 - II
         I0 = I0 - I
         T = ZERO
         DO 10 J = 1, I
              IJ = I0 + J
              T = T + L(IJ)*Y(J)
 10           CONTINUE
         X(I) = T
 20      CONTINUE
 999  RETURN
C  ***  LAST LINE OF  L7VML FOLLOWS  ***
      END
      SUBROUTINE  O7PRD(L, LS, P, S, W, Y, Z)
C
C  ***  FOR I = 1..L, SET S = S + W(I)*Y(.,I)*(Z(.,I)**T), I.E.,
C  ***        ADD W(I) TIMES THE OUTER PRODUCT OF Y(.,I) AND Z(.,I).
C
      INTEGER L, LS, P
      REAL S(LS), W(L), Y(P,L), Z(P,L)
C     DIMENSION S(P*(P+1)/2)
C
      INTEGER I, J, K, M
      REAL WK, YI, ZERO
      DATA ZERO/0.E+0/
C
      DO 30 K = 1, L
         WK = W(K)
         IF (WK .EQ. ZERO) GO TO 30
         M = 1
         DO 20 I = 1, P
              YI = WK * Y(I,K)
              DO 10 J = 1, I
                   S(M) = S(M) + YI*Z(J,K)
                   M = M + 1
 10                CONTINUE
 20           CONTINUE
 30      CONTINUE
C
 999  RETURN
C  ***  LAST LINE OF  O7PRD FOLLOWS  ***
      END
      SUBROUTINE  PARCK(ALG, D, IV, LIV, LV, N, V)
C
C  ***  CHECK ***SOL (VERSION 2.3) PARAMETERS, PRINT CHANGED VALUES  ***
C
C  ***  ALG = 1 FOR REGRESSION, ALG = 2 FOR GENERAL UNCONSTRAINED OPT.
C
      INTEGER ALG, LIV, LV, N
      INTEGER IV(LIV)
      REAL D(N), V(LV)
C
      REAL  R7MDC
      EXTERNAL  IVSET,  R7MDC, V7CPY, V7DFL
C  IVSET  -- SUPPLIES DEFAULT VALUES TO BOTH IV AND V.
C  R7MDC -- RETURNS MACHINE-DEPENDENT CONSTANTS.
C  V7CPY  -- COPIES ONE VECTOR TO ANOTHER.
C  V7DFL  -- SUPPLIES DEFAULT PARAMETER VALUES TO V ALONE.
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER ALG1, I, II, IV1, J, K, L, M, MIV1, MIV2, NDFALT, PARSV1,
     1        PU
      INTEGER IJMP, JLIM(4), MINIV(4), NDFLT(4)
      CHARACTER*1 VARNM(2), SH(2)
      CHARACTER*4 CNGD(3), DFLT(3), VN(2,34), WHICH(3)
      REAL BIG, MACHEP, TINY, VK, VM(34), VX(34), ZERO
C
C  ***  IV AND V SUBSCRIPTS  ***
C
      INTEGER ALGSAV, DINIT, DTYPE, DTYPE0, EPSLON, INITS, IVNEED,
     1        LASTIV, LASTV, LMAT, NEXTIV, NEXTV, NVDFLT, OLDN,
     2        PARPRT, PARSAV, PERM, PRUNIT, VNEED
C
C
      PARAMETER (ALGSAV=51, DINIT=38, DTYPE=16, DTYPE0=54, EPSLON=19,
     1           INITS=25, IVNEED=3, LASTIV=44, LASTV=45, LMAT=42,
     2           NEXTIV=46, NEXTV=47, NVDFLT=50, OLDN=38, PARPRT=20,
     3           PARSAV=49, PERM=58, PRUNIT=21, VNEED=4)
      SAVE BIG, MACHEP, TINY
C
      DATA BIG/0.E+0/, MACHEP/-1.E+0/, TINY/1.E+0/, ZERO/0.E+0/
      DATA VN(1,1),VN(2,1)/'EPSL','ON..'/
      DATA VN(1,2),VN(2,2)/'PHMN','FC..'/
      DATA VN(1,3),VN(2,3)/'PHMX','FC..'/
      DATA VN(1,4),VN(2,4)/'DECF','AC..'/
      DATA VN(1,5),VN(2,5)/'INCF','AC..'/
      DATA VN(1,6),VN(2,6)/'RDFC','MN..'/
      DATA VN(1,7),VN(2,7)/'RDFC','MX..'/
      DATA VN(1,8),VN(2,8)/'TUNE','R1..'/
      DATA VN(1,9),VN(2,9)/'TUNE','R2..'/
      DATA VN(1,10),VN(2,10)/'TUNE','R3..'/
      DATA VN(1,11),VN(2,11)/'TUNE','R4..'/
      DATA VN(1,12),VN(2,12)/'TUNE','R5..'/
      DATA VN(1,13),VN(2,13)/'AFCT','OL..'/
      DATA VN(1,14),VN(2,14)/'RFCT','OL..'/
      DATA VN(1,15),VN(2,15)/'XCTO','L...'/
      DATA VN(1,16),VN(2,16)/'XFTO','L...'/
      DATA VN(1,17),VN(2,17)/'LMAX','0...'/
      DATA VN(1,18),VN(2,18)/'LMAX','S...'/
      DATA VN(1,19),VN(2,19)/'SCTO','L...'/
      DATA VN(1,20),VN(2,20)/'DINI','T...'/
      DATA VN(1,21),VN(2,21)/'DTIN','IT..'/
      DATA VN(1,22),VN(2,22)/'D0IN','IT..'/
      DATA VN(1,23),VN(2,23)/'DFAC','....'/
      DATA VN(1,24),VN(2,24)/'DLTF','DC..'/
      DATA VN(1,25),VN(2,25)/'DLTF','DJ..'/
      DATA VN(1,26),VN(2,26)/'DELT','A0..'/
      DATA VN(1,27),VN(2,27)/'FUZZ','....'/
      DATA VN(1,28),VN(2,28)/'RLIM','IT..'/
      DATA VN(1,29),VN(2,29)/'COSM','IN..'/
      DATA VN(1,30),VN(2,30)/'HUBE','RC..'/
      DATA VN(1,31),VN(2,31)/'RSPT','OL..'/
      DATA VN(1,32),VN(2,32)/'SIGM','IN..'/
      DATA VN(1,33),VN(2,33)/'ETA0','....'/
      DATA VN(1,34),VN(2,34)/'BIAS','....'/
C
      DATA VM(1)/1.0E-3/, VM(2)/-0.99E+0/, VM(3)/1.0E-3/, VM(4)/1.0E-2/,
     1     VM(5)/1.2E+0/, VM(6)/1.E-2/, VM(7)/1.2E+0/, VM(8)/0.E+0/,
     2     VM(9)/0.E+0/, VM(10)/1.E-3/, VM(11)/-1.E+0/, VM(13)/0.E+0/,
     3     VM(15)/0.E+0/, VM(16)/0.E+0/, VM(19)/0.E+0/, VM(20)/-10.E+0/,
     4     VM(21)/0.E+0/, VM(22)/0.E+0/, VM(23)/0.E+0/, VM(27)/1.01E+0/,
     5     VM(28)/1.E+10/, VM(30)/0.E+0/, VM(31)/0.E+0/, VM(32)/0.E+0/,
     6     VM(34)/0.E+0/
      DATA VX(1)/0.9E+0/, VX(2)/-1.E-3/, VX(3)/1.E+1/, VX(4)/0.8E+0/,
     1     VX(5)/1.E+2/, VX(6)/0.8E+0/, VX(7)/1.E+2/, VX(8)/0.5E+0/,
     2     VX(9)/0.5E+0/, VX(10)/1.E+0/, VX(11)/1.E+0/, VX(14)/0.1E+0/,
     3     VX(15)/1.E+0/, VX(16)/1.E+0/, VX(19)/1.E+0/, VX(23)/1.E+0/,
     4     VX(24)/1.E+0/, VX(25)/1.E+0/, VX(26)/1.E+0/, VX(27)/1.E+10/,
     5     VX(29)/1.E+0/, VX(31)/1.E+0/, VX(32)/1.E+0/, VX(33)/1.E+0/,
     6     VX(34)/1.E+0/
C
      DATA VARNM(1)/'P'/, VARNM(2)/'P'/, SH(1)/'S'/, SH(2)/'H'/
      DATA CNGD(1),CNGD(2),CNGD(3)/'---C','HANG','ED V'/,
     1     DFLT(1),DFLT(2),DFLT(3)/'NOND','EFAU','LT V'/
      DATA IJMP/33/, JLIM(1)/0/, JLIM(2)/24/, JLIM(3)/0/, JLIM(4)/24/,
     1     NDFLT(1)/32/, NDFLT(2)/25/, NDFLT(3)/32/, NDFLT(4)/25/
      DATA MINIV(1)/82/, MINIV(2)/59/, MINIV(3)/103/, MINIV(4)/103/
C
C...............................  BODY  ................................
C
      PU = 0
      IF (PRUNIT .LE. LIV) PU = IV(PRUNIT)
      IF (ALGSAV .GT. LIV) GO TO 20
      IF (ALG .EQ. IV(ALGSAV)) GO TO 20
         IF (PU .NE. 0) WRITE(PU,10) ALG, IV(ALGSAV)
 10      FORMAT(/40H THE FIRST PARAMETER TO  IVSET SHOULD BE,I3,
     1          12H RATHER THAN,I3)
         IV(1) = 67
         GO TO 999
 20   IF (ALG .LT. 1 .OR. ALG .GT. 4) GO TO 340
      MIV1 = MINIV(ALG)
      IF (IV(1) .EQ. 15) GO TO 360
      ALG1 = MOD(ALG-1,2) + 1
      IF (IV(1) .EQ. 0) CALL  IVSET(ALG, IV, LIV, LV, V)
      IV1 = IV(1)
      IF (IV1 .NE. 13 .AND. IV1 .NE. 12) GO TO 30
      IF (PERM .LE. LIV) MIV1 = MAX0(MIV1, IV(PERM) - 1)
      IF (IVNEED .LE. LIV) MIV2 = MIV1 + MAX0(IV(IVNEED), 0)
      IF (LASTIV .LE. LIV) IV(LASTIV) = MIV2
      IF (LIV .LT. MIV1) GO TO 300
      IV(IVNEED) = 0
      IV(LASTV) = MAX0(IV(VNEED), 0) + IV(LMAT) - 1
      IV(VNEED) = 0
      IF (LIV .LT. MIV2) GO TO 300
      IF (LV .LT. IV(LASTV)) GO TO 320
 30   IF (IV1 .LT. 12 .OR. IV1 .GT. 14) GO TO 60
         IF (N .GE. 1) GO TO 50
              IV(1) = 81
              IF (PU .EQ. 0) GO TO 999
              WRITE(PU,40) VARNM(ALG1), N
 40           FORMAT(/8H /// BAD,A1,2H =,I5)
              GO TO 999
 50      IF (IV1 .NE. 14) IV(NEXTIV) = IV(PERM)
         IF (IV1 .NE. 14) IV(NEXTV) = IV(LMAT)
         IF (IV1 .EQ. 13) GO TO 999
         K = IV(PARSAV) - EPSLON
         CALL  V7DFL(ALG1, LV-K, V(K+1))
         IV(DTYPE0) = 2 - ALG1
         IV(OLDN) = N
         WHICH(1) = DFLT(1)
         WHICH(2) = DFLT(2)
         WHICH(3) = DFLT(3)
         GO TO 110
 60   IF (N .EQ. IV(OLDN)) GO TO 80
         IV(1) = 17
         IF (PU .EQ. 0) GO TO 999
         WRITE(PU,70) VARNM(ALG1), IV(OLDN), N
 70      FORMAT(/5H /// ,1A1,14H CHANGED FROM ,I5,4H TO ,I5)
         GO TO 999
C
 80   IF (IV1 .LE. 11 .AND. IV1 .GE. 1) GO TO 100
         IV(1) = 80
         IF (PU .NE. 0) WRITE(PU,90) IV1
 90      FORMAT(/13H ///  IV(1) =,I5,28H SHOULD BE BETWEEN 0 AND 14.)
         GO TO 999
C
 100  WHICH(1) = CNGD(1)
      WHICH(2) = CNGD(2)
      WHICH(3) = CNGD(3)
C
 110  IF (IV1 .EQ. 14) IV1 = 12
      IF (BIG .GT. TINY) GO TO 120
         TINY =  R7MDC(1)
         MACHEP =  R7MDC(3)
         BIG =  R7MDC(6)
         VM(12) = MACHEP
         VX(12) = BIG
         VX(13) = BIG
         VM(14) = MACHEP
         VM(17) = TINY
         VX(17) = BIG
         VM(18) = TINY
         VX(18) = BIG
         VX(20) = BIG
         VX(21) = BIG
         VX(22) = BIG
         VM(24) = MACHEP
         VM(25) = MACHEP
         VM(26) = MACHEP
         VX(28) =  R7MDC(5)
         VM(29) = MACHEP
         VX(30) = BIG
         VM(33) = MACHEP
 120  M = 0
      I = 1
      J = JLIM(ALG1)
      K = EPSLON
      NDFALT = NDFLT(ALG1)
      DO 150 L = 1, NDFALT
         VK = V(K)
         IF (VK .GE. VM(I) .AND. VK .LE. VX(I)) GO TO 140
              M = K
              IF (PU .NE. 0) WRITE(PU,130) VN(1,I), VN(2,I), K, VK,
     1                                    VM(I), VX(I)
 130          FORMAT(/6H ///  ,2A4,5H.. V(,I2,3H) =,E11.3,7H SHOULD,
     1               11H BE BETWEEN,E11.3,4H AND,E11.3)
 140     K = K + 1
         I = I + 1
         IF (I .EQ. J) I = IJMP
 150     CONTINUE
C
      IF (IV(NVDFLT) .EQ. NDFALT) GO TO 170
         IV(1) = 51
         IF (PU .EQ. 0) GO TO 999
         WRITE(PU,160) IV(NVDFLT), NDFALT
 160     FORMAT(/13H IV(NVDFLT) =,I5,13H RATHER THAN ,I5)
         GO TO 999
 170  IF ((IV(DTYPE) .GT. 0 .OR. V(DINIT) .GT. ZERO) .AND. IV1 .EQ. 12)
     1                  GO TO 200
      DO 190 I = 1, N
         IF (D(I) .GT. ZERO) GO TO 190
              M = 18
              IF (PU .NE. 0) WRITE(PU,180) I, D(I)
 180     FORMAT(/8H ///  D(,I3,3H) =,E11.3,19H SHOULD BE POSITIVE)
 190     CONTINUE
 200  IF (M .EQ. 0) GO TO 210
         IV(1) = M
         GO TO 999
C
 210  IF (PU .EQ. 0 .OR. IV(PARPRT) .EQ. 0) GO TO 999
      IF (IV1 .NE. 12 .OR. IV(INITS) .EQ. ALG1-1) GO TO 230
         M = 1
         WRITE(PU,220) SH(ALG1), IV(INITS)
 220     FORMAT(/22H NONDEFAULT VALUES..../5H INIT,A1,14H..... IV(25) =,
     1          I3)
 230  IF (IV(DTYPE) .EQ. IV(DTYPE0)) GO TO 250
         IF (M .EQ. 0) WRITE(PU,260) WHICH
         M = 1
         WRITE(PU,240) IV(DTYPE)
 240     FORMAT(20H DTYPE..... IV(16) =,I3)
 250  I = 1
      J = JLIM(ALG1)
      K = EPSLON
      L = IV(PARSAV)
      NDFALT = NDFLT(ALG1)
      DO 290 II = 1, NDFALT
         IF (V(K) .EQ. V(L)) GO TO 280
              IF (M .EQ. 0) WRITE(PU,260) WHICH
 260          FORMAT(/1H ,3A4,9HALUES..../)
              M = 1
              WRITE(PU,270) VN(1,I), VN(2,I), K, V(K)
 270          FORMAT(1X,2A4,5H.. V(,I2,3H) =,E15.7)
 280     K = K + 1
         L = L + 1
         I = I + 1
         IF (I .EQ. J) I = IJMP
 290     CONTINUE
C
      IV(DTYPE0) = IV(DTYPE)
      PARSV1 = IV(PARSAV)
      CALL  V7CPY(IV(NVDFLT), V(PARSV1), V(EPSLON))
      GO TO 999
C
 300  IV(1) = 15
      IF (PU .EQ. 0) GO TO 999
      WRITE(PU,310) LIV, MIV2
 310  FORMAT(/10H /// LIV =,I5,17H MUST BE AT LEAST,I5)
      IF (LIV .LT. MIV1) GO TO 999
      IF (LV .LT. IV(LASTV)) GO TO 320
      GO TO 999
C
 320  IV(1) = 16
      IF (PU .NE. 0) WRITE(PU,330) LV, IV(LASTV)
 330  FORMAT(/9H /// LV =,I5,17H MUST BE AT LEAST,I5)
      GO TO 999
C
 340  IV(1) = 67
      IF (PU .NE. 0) WRITE(PU,350) ALG
 350  FORMAT(/10H /// ALG =,I5,21H MUST BE 1 2, 3, OR 4)
      GO TO 999
 360  IF (PU .NE. 0) WRITE(PU,370) LIV, MIV1
 370  FORMAT(/10H /// LIV =,I5,17H MUST BE AT LEAST,I5,
     1       37H TO COMPUTE TRUE MIN. LIV AND MIN. LV)
      IF (LASTIV .LE. LIV) IV(LASTIV) = MIV1
      IF (LASTV .LE. LIV) IV(LASTV) = 0
C
 999  RETURN
C  ***  LAST LINE OF  PARCK FOLLOWS  ***
      END
      SUBROUTINE  Q7ADR(P, QTR, RMAT, W, Y)
C
C  ***  ADD ROW W TO QR FACTORIZATION WITH R MATRIX RMAT AND
C  ***  Q**T * RESIDUAL = QTR.  Y = NEW COMPONENT OF RESIDUAL
C  ***  CORRESPONDING TO W.
C
      INTEGER P
      REAL QTR(P), RMAT(1), W(P), Y
C     DIMENSION RMAT(P*(P+1)/2)
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER I, II, IJ, IP1, J
      REAL RI, RW, T, U1, U2, V, WI, WR
C
      REAL ONE, ZERO
      PARAMETER (ONE=1.E+0, ZERO=0.E+0)
C
C------------------------------ BODY -----------------------------------
C
      II = 0
      DO 60 I = 1, P
         II = II+I
         WI = W(I)
         IF (WI .EQ. ZERO) GOTO  60
         RI = RMAT(II)
         IF (RI .NE. ZERO) GOTO 20
            IJ = II
C           *** SWAP W AND ROW I OF RMAT ***
            DO 10 J = I, P
               T = RMAT(IJ)
               RMAT(IJ) = W(J)
               W(J) = T
               IJ = IJ+J
 10            CONTINUE
            T = QTR(I)
            QTR(I) = Y
            Y = T
            GO TO 60
 20      IP1 = I+1
         IJ = II+I
         IF ( ABS(WI) .LE.  ABS(RI)) GO TO 40
            RW = RI/WI
            T =  SQRT(ONE+RW**2)
            IF (RW .GT. ZERO) T = -T
            V = RW-T
            U1 = ONE/T
            U2 = ONE/(T*V)
            RMAT(II) = WI*T
            T = Y+V*QTR(I)
            QTR(I) = QTR(I)+T*U1
            Y = Y+T*U2
            IF (IP1 .GT. P) GO TO 60
            DO 30 J = IP1, P
               T = W(J)+V*RMAT(IJ)
               RMAT(IJ) = RMAT(IJ)+T*U1
               W(J) = W(J)+T*U2
               IJ = IJ+J
 30            CONTINUE
            GO TO 60
C
C        *** AT THIS POINT WE MUST HAVE ABS(WI) .LE. ABS(RI)...
C
 40      WR = WI/RI
         T = - SQRT(ONE+WR**2)
         V = WR/(ONE-T)
         U1 = ONE/T-ONE
         U2 = V*U1
         RMAT(II) = RI*T
         T = QTR(I)+V*Y
         QTR(I) = QTR(I)+T*U1
         Y = Y+T*U2
         IF (IP1 .GT. P) GO TO 60
         DO 50 J = IP1, P
            T = RMAT(IJ)+V*W(J)
            RMAT(IJ) = RMAT(IJ)+T*U1
            W(J) = W(J)+T*U2
            IJ = IJ+J
 50         CONTINUE
 60      CONTINUE
 999  RETURN
      END
      REAL FUNCTION  RLDST(P, D, X, X0)
C
C  ***  COMPUTE AND RETURN RELATIVE DIFFERENCE BETWEEN X AND X0  ***
C  ***  NL2SOL VERSION 2.2  ***
C
      INTEGER P
      REAL D(P), X(P), X0(P)
C
      INTEGER I
      REAL EMAX, T, XMAX, ZERO
      PARAMETER (ZERO=0.E+0)
C
C  ***  BODY  ***
C
      EMAX = ZERO
      XMAX = ZERO
      DO 10 I = 1, P
         T =  ABS(D(I) * (X(I) - X0(I)))
         IF (EMAX .LT. T) EMAX = T
         T = D(I) * ( ABS(X(I)) +  ABS(X0(I)))
         IF (XMAX .LT. T) XMAX = T
 10      CONTINUE
       RLDST = ZERO
      IF (XMAX .GT. ZERO)  RLDST = EMAX / XMAX
 999  RETURN
C  ***  LAST LINE OF  RLDST FOLLOWS  ***
      END
      SUBROUTINE  S7LUP(A, COSMIN, P, SIZE, STEP, U, W, WCHMTD, WSCALE,
     1                  Y)
C
C  ***  UPDATE SYMMETRIC  A  SO THAT  A * STEP = Y  ***
C  ***  (LOWER TRIANGLE OF  A  STORED ROWWISE       ***
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER P
      REAL A(1), COSMIN, SIZE, STEP(P), U(P), W(P),
     1                 WCHMTD(P), WSCALE, Y(P)
C     DIMENSION A(P*(P+1)/2)
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER I, J, K
      REAL DENMIN, SDOTWM, T, UI, WI
C
C     ***  CONSTANTS  ***
      REAL HALF, ONE, ZERO
C
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      REAL  D7TPR,  V2NRM
      EXTERNAL  D7TPR,  S7LVM,  V2NRM
C
      PARAMETER (HALF=0.5E+0, ONE=1.E+0, ZERO=0.E+0)
C
C-----------------------------------------------------------------------
C
      SDOTWM =  D7TPR(P, STEP, WCHMTD)
      DENMIN = COSMIN *  V2NRM(P,STEP) *  V2NRM(P,WCHMTD)
      WSCALE = ONE
      IF (DENMIN .NE. ZERO) WSCALE =   MIN(ONE,  ABS(SDOTWM/DENMIN))
      T = ZERO
      IF (SDOTWM .NE. ZERO) T = WSCALE / SDOTWM
      DO 10 I = 1, P
 10      W(I) = T * WCHMTD(I)
      CALL  S7LVM(P, U, A, STEP)
      T = HALF * (SIZE *  D7TPR(P, STEP, U)  -   D7TPR(P, STEP, Y))
      DO 20 I = 1, P
 20      U(I) = T*W(I) + Y(I) - SIZE*U(I)
C
C  ***  SET  A = A + U*(W**T) + W*(U**T)  ***
C
      K = 1
      DO 40 I = 1, P
         UI = U(I)
         WI = W(I)
         DO 30 J = 1, I
              A(K) = SIZE*A(K) + UI*W(J) + WI*U(J)
              K = K + 1
 30           CONTINUE
 40      CONTINUE
C
 999  RETURN
C  ***  LAST LINE OF  S7LUP FOLLOWS  ***
      END
      SUBROUTINE  S7LVM(P, Y, S, X)
C
C  ***  SET  Y = S * X,  S = P X P SYMMETRIC MATRIX.  ***
C  ***  LOWER TRIANGLE OF  S  STORED ROWWISE.         ***
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER P
      REAL S(1), X(P), Y(P)
C     DIMENSION S(P*(P+1)/2)
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER I, IM1, J, K
      REAL XI
C
C
C  ***  EXTERNAL FUNCTION  ***
C
      REAL  D7TPR
      EXTERNAL  D7TPR
C
C-----------------------------------------------------------------------
C
      J = 1
      DO 10 I = 1, P
         Y(I) =  D7TPR(I, S(J), X)
         J = J + I
 10      CONTINUE
C
      IF (P .LE. 1) GO TO 999
      J = 1
      DO 40 I = 2, P
         XI = X(I)
         IM1 = I - 1
         J = J + 1
         DO 30 K = 1, IM1
              Y(K) = Y(K) + S(J)*XI
              J = J + 1
 30           CONTINUE
 40      CONTINUE
C
 999  RETURN
C  ***  LAST LINE OF  S7LVM FOLLOWS  ***
      END
      SUBROUTINE  V2AXY(P, W, A, X, Y)
C
C  ***  SET W = A*X + Y  --  W, X, Y = P-VECTORS, A = SCALAR  ***
C
      INTEGER P
      REAL A, W(P), X(P), Y(P)
C
      INTEGER I
C
      DO 10 I = 1, P
 10      W(I) = A*X(I) + Y(I)
      RETURN
      END
      REAL FUNCTION  V2NRM(P, X)
C
C  ***  RETURN THE 2-NORM OF THE P-VECTOR X, TAKING  ***
C  ***  CARE TO AVOID THE MOST LIKELY UNDERFLOWS.    ***
C
      INTEGER P
      REAL X(P)
C
      INTEGER I, J
      REAL ONE, R, SCALE, SQTETA, T, XI, ZERO
      REAL  R7MDC
      EXTERNAL  R7MDC
C
      PARAMETER (ONE=1.E+0, ZERO=0.E+0)
      SAVE SQTETA
      DATA SQTETA/0.E+0/
C
      IF (P .GT. 0) GO TO 10
          V2NRM = ZERO
         GO TO 999
 10   DO 20 I = 1, P
         IF (X(I) .NE. ZERO) GO TO 30
 20      CONTINUE
       V2NRM = ZERO
      GO TO 999
C
 30   SCALE =  ABS(X(I))
      IF (I .LT. P) GO TO 40
          V2NRM = SCALE
         GO TO 999
 40   T = ONE
      IF (SQTETA .EQ. ZERO) SQTETA =  R7MDC(2)
C
C     ***  SQTETA IS (SLIGHTLY LARGER THAN) THE SQUARE ROOT OF THE
C     ***  SMALLEST POSITIVE FLOATING POINT NUMBER ON THE MACHINE.
C     ***  THE TESTS INVOLVING SQTETA ARE DONE TO PREVENT UNDERFLOWS.
C
      J = I + 1
      DO 60 I = J, P
         XI =  ABS(X(I))
         IF (XI .GT. SCALE) GO TO 50
              R = XI / SCALE
              IF (R .GT. SQTETA) T = T + R*R
              GO TO 60
 50           R = SCALE / XI
              IF (R .LE. SQTETA) R = ZERO
              T = ONE  +  T * R*R
              SCALE = XI
 60      CONTINUE
C
       V2NRM = SCALE *  SQRT(T)
 999  RETURN
C  ***  LAST LINE OF  V2NRM FOLLOWS  ***
      END
      SUBROUTINE  V7CPY(P, Y, X)
C
C  ***  SET Y = X, WHERE X AND Y ARE P-VECTORS  ***
C
      INTEGER P
      REAL X(P), Y(P)
C
      INTEGER I
C
      DO 10 I = 1, P
 10      Y(I) = X(I)
      RETURN
      END
      SUBROUTINE  V7DFL(ALG, LV, V)
C
C  ***  SUPPLY ***SOL (VERSION 2.3) DEFAULT VALUES TO V  ***
C
C  ***  ALG = 1 MEANS REGRESSION CONSTANTS.
C  ***  ALG = 2 MEANS GENERAL UNCONSTRAINED OPTIMIZATION CONSTANTS.
C
      INTEGER ALG, LV
      REAL V(LV)
C
      REAL  R7MDC
      EXTERNAL  R7MDC
C  R7MDC... RETURNS MACHINE-DEPENDENT CONSTANTS
C
      REAL MACHEP, MEPCRT, ONE, SQTEPS, THREE
C
C  ***  SUBSCRIPTS FOR V  ***
C
      INTEGER AFCTOL, BIAS, COSMIN, DECFAC, DELTA0, DFAC, DINIT, DLTFDC,
     1        DLTFDJ, DTINIT, D0INIT, EPSLON, ETA0, FUZZ, HUBERC,
     2        INCFAC, LMAX0, LMAXS, PHMNFC, PHMXFC, RDFCMN, RDFCMX,
     3        RFCTOL, RLIMIT, RSPTOL, SCTOL, SIGMIN, TUNER1, TUNER2,
     4        TUNER3, TUNER4, TUNER5, XCTOL, XFTOL
C
      PARAMETER (ONE=1.E+0, THREE=3.E+0)
C
C  ***  V SUBSCRIPT VALUES  ***
C
      PARAMETER (AFCTOL=31, BIAS=43, COSMIN=47, DECFAC=22, DELTA0=44,
     1           DFAC=41, DINIT=38, DLTFDC=42, DLTFDJ=43, DTINIT=39,
     2           D0INIT=40, EPSLON=19, ETA0=42, FUZZ=45, HUBERC=48,
     3           INCFAC=23, LMAX0=35, LMAXS=36, PHMNFC=20, PHMXFC=21,
     4           RDFCMN=24, RDFCMX=25, RFCTOL=32, RLIMIT=46, RSPTOL=49,
     5           SCTOL=37, SIGMIN=50, TUNER1=26, TUNER2=27, TUNER3=28,
     6           TUNER4=29, TUNER5=30, XCTOL=33, XFTOL=34)
C
C-------------------------------  BODY  --------------------------------
C
      MACHEP =  R7MDC(3)
      V(AFCTOL) = 1.E-20
      IF (MACHEP .GT. 1.E-10) V(AFCTOL) = MACHEP**2
      V(DECFAC) = 0.5E+0
      SQTEPS =  R7MDC(4)
      V(DFAC) = 0.6E+0
      V(DTINIT) = 1.E-6
      MEPCRT = MACHEP ** (ONE/THREE)
      V(D0INIT) = 1.E+0
      V(EPSLON) = 0.1E+0
      V(INCFAC) = 2.E+0
      V(LMAX0) = 1.E+0
      V(LMAXS) = 1.E+0
      V(PHMNFC) = -0.1E+0
      V(PHMXFC) = 0.1E+0
      V(RDFCMN) = 0.1E+0
      V(RDFCMX) = 4.E+0
      V(RFCTOL) =   MAX(1.E-10, MEPCRT**2)
      V(SCTOL) = V(RFCTOL)
      V(TUNER1) = 0.1E+0
      V(TUNER2) = 1.E-4
      V(TUNER3) = 0.75E+0
      V(TUNER4) = 0.5E+0
      V(TUNER5) = 0.75E+0
      V(XCTOL) = SQTEPS
      V(XFTOL) = 1.E+2 * MACHEP
C
      IF (ALG .GE. 2) GO TO 10
C
C  ***  REGRESSION  VALUES
C
      V(COSMIN) =   MAX(1.E-6, 1.E+2 * MACHEP)
      V(DINIT) = 0.E+0
      V(DELTA0) = SQTEPS
      V(DLTFDC) = MEPCRT
      V(DLTFDJ) = SQTEPS
      V(FUZZ) = 1.5E+0
      V(HUBERC) = 0.7E+0
      V(RLIMIT) =  R7MDC(5)
      V(RSPTOL) = 1.E-3
      V(SIGMIN) = 1.E-4
      GO TO 999
C
C  ***  GENERAL OPTIMIZATION VALUES
C
 10   V(BIAS) = 0.8E+0
      V(DINIT) = -1.0E+0
      V(ETA0) = 1.0E+3 * MACHEP
C
 999  RETURN
C  ***  LAST LINE OF  V7DFL FOLLOWS  ***
      END
      SUBROUTINE  V7SCL(N, X, A, Y)
C
C  ***  SET X(I) = A*Y(I), I = 1(1)N  ***
C
      INTEGER N
      REAL A, X(N), Y(N)
C
      INTEGER I
C
      DO 10 I = 1, N
 10       X(I) = A * Y(I)
 999    RETURN
C  ***  LAST LINE OF  V7SCL FOLLOWS  ***
      END
      SUBROUTINE  V7SCP(P, Y, S)
C
C  ***  SET P-VECTOR Y TO SCALAR S  ***
C
      INTEGER P
      REAL S, Y(P)
C
      INTEGER I
C
      DO 10 I = 1, P
 10      Y(I) = S
      RETURN
      END
      REAL FUNCTION  VSUM(N, X)
      INTEGER N
      REAL X(N)
      INTEGER I
C
       VSUM = 0.E+0
      DO 10 I = 1, N
 10       VSUM =  VSUM + X(I)
      END
      LOGICAL FUNCTION STOPX(IDUMMY)
C     *****PARAMETERS...
      INTEGER IDUMMY
C
C     ..................................................................
C
C     *****PURPOSE...
C     THIS FUNCTION MAY SERVE AS THE STOPX (ASYNCHRONOUS INTERRUPTION)
C     FUNCTION FOR THE NL2SOL (NONLINEAR LEAST-SQUARES) PACKAGE AT
C     THOSE INSTALLATIONS WHICH DO NOT WISH TO IMPLEMENT A
C     DYNAMIC STOPX.
C
C     *****ALGORITHM NOTES...
C     AT INSTALLATIONS WHERE THE NL2SOL SYSTEM IS USED
C     INTERACTIVELY, THIS DUMMY STOPX SHOULD BE REPLACED BY A
C     FUNCTION THAT RETURNS .TRUE. IF AND ONLY IF THE INTERRUPT
C     (BREAK) KEY HAS BEEN PRESSED SINCE THE LAST CALL ON STOPX.
C
C     ..................................................................
C
      STOPX = .FALSE.
      RETURN
      END
