
* ======================================================================
* NIST Guide to Available Math Software.
* Fullsource for module DHUMSL from package CMLIB.
* Retrieved from CAMSUN on Sat Aug 30 09:32:07 1997.
* ======================================================================
      SUBROUTINE DHUMSL(N, D, X, CALCF, CALCGH, IV, LIV, LV, V,
     1                  UIPARM, URPARM, UFPARM)
C
C  ***  MINIMIZE GENERAL UNCONSTRAINED OBJECTIVE FUNCTION USING   ***
C  ***  (ANALYTIC) GRADIENT AND HESSIAN PROVIDED BY THE CALLER.   ***
C
      INTEGER LIV, LV, N
      INTEGER IV(LIV), UIPARM(1)
      DOUBLE PRECISION D(N), X(N), V(LV), URPARM(1)
C     DIMENSION V(78 + N*(N+12)), UIPARM(*), URPARM(*)
      EXTERNAL CALCF, CALCGH, UFPARM
C
C------------------------------  DISCUSSION  ---------------------------
C
C        THIS ROUTINE IS LIKE DSUMSL, EXCEPT THAT THE SUBROUTINE PARA-
C     METER CALCG OF DSUMSL (WHICH COMPUTES THE GRADIENT OF THE OBJEC-
C     TIVE FUNCTION) IS REPLACED BY THE SUBROUTINE PARAMETER CALCGH,
C     WHICH COMPUTES BOTH THE GRADIENT AND (LOWER TRIANGLE OF THE)
C     HESSIAN OF THE OBJECTIVE FUNCTION.  THE CALLING SEQUENCE IS...
C             CALL CALCGH(N, X, NF, G, H, UIPARM, URPARM, UFPARM)
C     PARAMETERS N, X, NF, G, UIPARM, URPARM, AND UFPARM ARE THE SAME
C     AS FOR DSUMSL, WHILE H IS AN ARRAY OF LENGTH N*(N+1)/2 IN WHICH
C     CALCGH MUST STORE THE LOWER TRIANGLE OF THE HESSIAN AT X.  START-
C     ING AT H(1), CALCGH MUST STORE THE HESSIAN ENTRIES IN THE ORDER
C     (1,1), (2,1), (2,2), (3,1), (3,2), (3,3), ...
C        THE VALUE PRINTED (BY DITSUM) IN THE COLUMN LABELLED STPPAR
C     IS THE LEVENBERG-MARQUARDT USED IN COMPUTING THE CURRENT STEP.
C     ZERO MEANS A FULL NEWTON STEP.  IF THE SPECIAL CASE DESCRIBED IN
C     REF. 1 IS DETECTED, THEN STPPAR IS NEGATED.  THE VALUE PRINTED
C     IN THE COLUMN LABELLED NPRELDF IS ZERO IF THE CURRENT HESSIAN
C     IS NOT POSITIVE DEFINITE.
C        IT SOMETIMES PROVES WORTHWHILE TO LET D BE DETERMINED FROM THE
C     DIAGONAL OF THE HESSIAN MATRIX BY SETTING IV(DTYPE) = 1 AND
C     V(DINIT) = 0.  THE FOLLOWING IV AND V COMPONENTS ARE RELEVANT...
C
C IV(DTOL)..... IV(59) GIVES THE STARTING SUBSCRIPT IN V OF THE DTOL
C             ARRAY USED WHEN D IS UPDATED.  (IV(DTOL) CAN BE
C             INITIALIZED BY CALLING DHUMSL WITH IV(1) = 13.)
C IV(DTYPE).... IV(16) TELLS HOW THE SCALE VECTOR D SHOULD BE CHOSEN.
C             IV(DTYPE) .LE. 0 MEANS THAT D SHOULD NOT BE UPDATED, AND
C             IV(DTYPE) .GE. 1 MEANS THAT D SHOULD BE UPDATED AS
C             DESCRIBED BELOW WITH V(DFAC).  DEFAULT = 0.
C V(DFAC)..... V(41) AND THE DTOL AND D0 ARRAYS (SEE V(DTINIT) AND
C             V(D0INIT)) ARE USED IN UPDATING THE SCALE VECTOR D WHEN
C             IV(DTYPE) .GT. 0.  (D IS INITIALIZED ACCORDING TO
C             V(DINIT), DESCRIBED IN DSUMSL.)  LET
C                  D1(I) = MAX(SQRT(ABS(H(I,I))), V(DFAC)*D(I)),
C             WHERE H(I,I) IS THE I-TH DIAGONAL ELEMENT OF THE CURRENT
C             HESSIAN.  IF IV(DTYPE) = 1, THEN D(I) IS SET TO D1(I)
C             UNLESS D1(I) .LT. DTOL(I), IN WHICH CASE D(I) IS SET TO
C                  MAX(D0(I), DTOL(I)).
C             IF IV(DTYPE) .GE. 2, THEN D IS UPDATED DURING THE FIRST
C             ITERATION AS FOR IV(DTYPE) = 1 (AFTER ANY INITIALIZATION
C             DUE TO V(DINIT)) AND IS LEFT UNCHANGED THEREAFTER.
C             DEFAULT = 0.6.
C V(DTINIT)... V(39), IF POSITIVE, IS THE VALUE TO WHICH ALL COMPONENTS
C             OF THE DTOL ARRAY (SEE V(DFAC)) ARE INITIALIZED.  IF
C             V(DTINIT) = 0, THEN IT IS ASSUMED THAT THE CALLER HAS
C             STORED DTOL IN V STARTING AT V(IV(DTOL)).
C             DEFAULT = 10**-6.
C V(D0INIT)... V(40), IF POSITIVE, IS THE VALUE TO WHICH ALL COMPONENTS
C             OF THE D0 VECTOR (SEE V(DFAC)) ARE INITIALIZED.  IF
C             V(DFAC) = 0, THEN IT IS ASSUMED THAT THE CALLER HAS
C             STORED D0 IN V STARTING AT V(IV(DTOL)+N).  DEFAULT = 1.0.
C
C  ***  REFERENCE  ***
C
C 1. GAY, D.M. (1981), COMPUTING OPTIMAL LOCALLY CONSTRAINED STEPS,
C         SIAM J. SCI. STATIST. COMPUT. 2, PP. 186-197.
C.
C  ***  GENERAL  ***
C
C     CODED BY DAVID M. GAY (WINTER 1980).  REVISED SEPT. 1982.
C     THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH SUPPORTED
C     IN PART BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS
C     MCS-7600324 AND MCS-7906671.
C
C----------------------------  DECLARATIONS  ---------------------------
C
      EXTERNAL DDEFLT, DHUMIT
C
C DDEFLT... PROVIDES DEFAULT INPUT VALUES FOR IV AND V.
C DHUMIT... REVERSE-COMMUNICATION ROUTINE THAT DOES DHUMSL ALGORITHM.
C
      INTEGER G1, H1, IV1, LH, NF
      DOUBLE PRECISION F
C
C  ***  SUBSCRIPTS FOR IV   ***
C
      INTEGER G, H, NEXTV, NFCALL, NFGCAL, TOOBIG, VNEED
C
C/6
      DATA NEXTV/47/, NFCALL/6/, NFGCAL/7/, G/28/, H/56/, TOOBIG/2/,
     1     VNEED/4/
C/7
C     PARAMETER (NEXTV=47, NFCALL=6, NFGCAL=7, G=28, H=56, TOOBIG=2,
C    1           VNEED=4)
C/
C
C+++++++++++++++++++++++++++++++  BODY  ++++++++++++++++++++++++++++++++
C
      LH = N * (N + 1) / 2
      IF (IV(1) .EQ. 0) CALL DDEFLT(2, IV, LIV, LV, V)
      IV(VNEED) = IV(VNEED) + N*(N+3)/2
      IV1 = IV(1)
      IF (IV1 .EQ. 14) GO TO 10
      IF (IV1 .GT. 2 .AND. IV1 .LT. 12) GO TO 10
      G1 = 1
      H1 = 1
      IF (IV1 .EQ. 12) IV(1) = 13
      GO TO 20
C
 10   G1 = IV(G)
      H1 = IV(H)
C
 20   CALL DHUMIT(D, F, V(G1), V(H1), IV, LH, LIV, LV, N, V, X)
      IF (IV(1) - 2) 30, 40, 50
C
 30   NF = IV(NFCALL)
      CALL CALCF(N, X, NF, F, UIPARM, URPARM, UFPARM)
      IF (NF .LE. 0) IV(TOOBIG) = 1
      GO TO 20
C
 40   CALL CALCGH(N, X, IV(NFGCAL), V(G1), V(H1), UIPARM, URPARM,
     1            UFPARM)
      GO TO 20
C
 50   IF (IV(1) .NE. 14) GO TO 999
C
C  ***  STORAGE ALLOCATION
C
      IV(G) = IV(NEXTV)
      IV(H) = IV(G) + N
      IV(NEXTV) = IV(H) + N*(N+1)/2
      IF (IV1 .NE. 13) GO TO 10
C
 999  RETURN
C  ***  LAST CARD OF DHUMSL FOLLOWS  ***
      END
      SUBROUTINE DDEFLT(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)
      DOUBLE PRECISION V(LV)
C
      EXTERNAL  DVDFLT
C DVDFLT.... PROVIDES DEFAULT VALUES TO V.
C
      INTEGER MIV, MV
      INTEGER MINIV(2), MINV(2)
C
C  ***  SUBSCRIPTS FOR IV  ***
C
      INTEGER ALGSAV, COVPRT, COVREQ, DTYPE, HC, IERR, INITH, INITS,
     1        IPIVOT, IVNEED, LASTIV, LASTV, LMAT, MXFCAL, MXITER,
     2        NFCOV, NGCOV, NVDFLT, OUTLEV, PARPRT, PARSAV, PERM,
     3        PRUNIT, QRTYP, RDREQ, RMAT, SOLPRT, STATPR, VNEED,
     4        VSAVE, X0PRT
C
C  ***  IV SUBSCRIPT VALUES  ***
C
C/6
      DATA ALGSAV/51/, COVPRT/14/, COVREQ/15/, DTYPE/16/, HC/71/,
     1     IERR/75/, INITH/25/, INITS/25/, IPIVOT/76/, IVNEED/3/,
     2     LASTIV/44/, LASTV/45/, LMAT/42/, MXFCAL/17/, MXITER/18/,
     3     NFCOV/52/, NGCOV/53/, NVDFLT/50/, OUTLEV/19/, PARPRT/20/,
     4     PARSAV/49/, PERM/58/, PRUNIT/21/, QRTYP/80/, RDREQ/57/,
     5     RMAT/78/, SOLPRT/22/, STATPR/23/, VNEED/4/, VSAVE/60/,
     6     X0PRT/24/
C/7
C     PARAMETER (ALGSAV=51, COVPRT=14, COVREQ=15, DTYPE=16, HC=71,
C    1           IERR=75, INITH=25, INITS=25, IPIVOT=76, IVNEED=3,
C    2           LASTIV=44, LASTV=45, LMAT=42, MXFCAL=17, MXITER=18,
C    3           NFCOV=52, NGCOV=53, NVDFLT=50, OUTLEV=19, PARPRT=20,
C    4           PARSAV=49, PERM=58, PRUNIT=21, QRTYP=80, RDREQ=57,
C    5           RMAT=78, SOLPRT=22, STATPR=23, VNEED=4, VSAVE=60,
C    6           X0PRT=24)
C/
      DATA MINIV(1)/80/, MINIV(2)/59/, MINV(1)/98/, MINV(2)/71/
C
C-------------------------------  BODY  --------------------------------
C
      IF (ALG .LT. 1 .OR. ALG .GT. 2) GO TO 40
      MIV = MINIV(ALG)
      IF (LIV .LT. MIV) GO TO 20
      MV = MINV(ALG)
      IF (LV .LT. MV) GO TO 30
      CALL DVDFLT(ALG, LV, V)
      IV(1) = 12
      IV(ALGSAV) = ALG
      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(PRUNIT) = I1MACH(2)
      IV(SOLPRT) = 1
      IV(STATPR) = 1
      IV(VNEED) = 0
      IV(X0PRT) = 1
C
      IF (ALG .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(PARSAV) = 67
      IV(QRTYP) = 1
      IV(RDREQ) = 3
      IV(RMAT) = 0
      IV(VSAVE) = 58
      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
      GO TO 999
C
 20   IV(1) = 15
      GO TO 999
C
 30   IV(1) = 16
      GO TO 999
C
 40   IV(1) = 67
C
 999  RETURN
C  ***  LAST CARD OF DDEFLT FOLLOWS  ***
      END
      SUBROUTINE DHUMIT(D, FX, G, H, IV, LH, LIV, LV, N, V, X)
C
C  ***  CARRY OUT DHUMSL (UNCONSTRAINED MINIMIZATION) ITERATIONS, USING
C  ***  HESSIAN MATRIX PROVIDED BY THE CALLER.
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER LH, LIV, LV, N
      INTEGER IV(LIV)
      DOUBLE PRECISION D(N), FX, G(N), H(LH), V(LV), X(N)
C
C--------------------------  PARAMETER USAGE  --------------------------
C
C D.... SCALE VECTOR.
C FX... FUNCTION VALUE.
C G.... GRADIENT VECTOR.
C H.... LOWER TRIANGLE OF THE HESSIAN, STORED ROWWISE.
C IV... INTEGER VALUE ARRAY.
C LH... LENGTH OF H = P*(P+1)/2.
C LIV.. LENGTH OF IV (AT LEAST 60).
C LV... LENGTH OF V (AT LEAST 78 + N*(N+21)/2).
C N.... NUMBER OF VARIABLES (COMPONENTS IN X AND G).
C V.... FLOATING-POINT VALUE ARRAY.
C X.... PARAMETER VECTOR.
C
C  ***  DISCUSSION  ***
C
C        PARAMETERS IV, N, V, AND X ARE THE SAME AS THE CORRESPONDING
C     ONES TO DHUMSL (WHICH SEE), EXCEPT THAT V CAN BE SHORTER (SINCE
C     THE PART OF V THAT DHUMSL USES FOR STORING G AND H IS NOT NEEDED).
C     MOREOVER, COMPARED WITH DHUMSL, IV(1) MAY HAVE THE TWO ADDITIONAL
C     OUTPUT VALUES 1 AND 2, WHICH ARE EXPLAINED BELOW, AS IS THE USE
C     OF IV(TOOBIG) AND IV(NFGCAL).  THE VALUE IV(G), WHICH IS AN
C     OUTPUT VALUE FROM DHUMSL, IS NOT REFERENCED BY DHUMIT OR THE
C     SUBROUTINES IT CALLS.
C
C IV(1) = 1 MEANS THE CALLER SHOULD SET FX TO F(X), THE FUNCTION VALUE
C             AT X, AND CALL DHUMIT AGAIN, HAVING CHANGED NONE OF THE
C             OTHER PARAMETERS.  AN EXCEPTION OCCURS IF F(X) CANNOT BE
C             COMPUTED (E.G. IF OVERFLOW WOULD OCCUR), WHICH MAY HAPPEN
C             BECAUSE OF AN OVERSIZED STEP.  IN THIS CASE THE CALLER
C             SHOULD SET IV(TOOBIG) = IV(2) TO 1, WHICH WILL CAUSE
C             DHUMIT TO IGNORE FX AND TRY A SMALLER STEP.  THE PARA-
C             METER NF THAT DHUMSL PASSES TO CALCF (FOR POSSIBLE USE BY
C             CALCGH) IS A COPY OF IV(NFCALL) = IV(6).
C IV(1) = 2 MEANS THE CALLER SHOULD SET G TO G(X), THE GRADIENT OF F AT
C             X, AND H TO THE LOWER TRIANGLE OF H(X), THE HESSIAN OF F
C             AT X, AND CALL DHUMIT AGAIN, HAVING CHANGED NONE OF THE
C             OTHER PARAMETERS EXCEPT PERHAPS THE SCALE VECTOR D.
C                  THE PARAMETER NF THAT DHUMSL PASSES TO CALCG IS
C             IV(NFGCAL) = IV(7).  IF G(X) AND H(X) CANNOT BE EVALUATED,
C             THEN THE CALLER MAY SET IV(NFGCAL) TO 0, IN WHICH CASE
C             DHUMIT WILL RETURN WITH IV(1) = 65.
C                  NOTE -- DHUMIT OVERWRITES H WITH THE LOWER TRIANGLE
C             OF  DIAG(D)**-1 * H(X) * DIAG(D)**-1.
C.
C  ***  GENERAL  ***
C
C     CODED BY DAVID M. GAY (WINTER 1980).  REVISED SEPT. 1982.
C     THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH SUPPORTED
C     IN PART BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS
C     MCS-7600324 AND MCS-7906671.
C
C        (SEE DSUMSL AND DHUMSL FOR REFERENCES.)
C
C+++++++++++++++++++++++++++  DECLARATIONS  ++++++++++++++++++++++++++++
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER DG1, DUMMY, I, J, K, L, LSTGST, NN1O2, STEP1,
     1        TEMP1, W1, X01
      DOUBLE PRECISION T
C
C     ***  CONSTANTS  ***
C
      DOUBLE PRECISION ONE, ZERO
C
C  ***  NO INTRINSIC FUNCTIONS  ***
C
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      EXTERNAL DASSST, DDEFLT, DDUPDU, DGQTST, DITSUM, DPARCK,
     1         DSLVMU, DSTOPX, DVAXPY, DVSCPY
      LOGICAL DSTOPX
      DOUBLE PRECISION DDOT, DNRM2
C
C DASSST.... ASSESSES CANDIDATE STEP.
C DDEFLT.... PROVIDES DEFAULT IV AND V INPUT VALUES.
C DDUPDU.... UPDATES SCALE VECTOR D.
C DGQTST.... COMPUTES OPTIMALLY LOCALLY CONSTRAINED STEP.
C DITSUM.... PRINTS ITERATION SUMMARY AND INFO ON INITIAL AND FINAL X.
C DPARCK.... CHECKS VALIDITY OF INPUT IV AND V VALUES.
C DSLVMU... MULTIPLIES SYMMETRIC MATRIX TIMES VECTOR, GIVEN THE LOWER
C             TRIANGLE OF THE MATRIX.
C DSTOPX.... RETURNS .TRUE. IF THE BREAK KEY HAS BEEN PRESSED.
C DVAXPY.... COMPUTES SCALAR TIMES ONE VECTOR PLUS ANOTHER.
C DVSCPY... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR.
C
C  ***  SUBSCRIPTS FOR IV AND V  ***
C
      INTEGER CNVCOD, DG, DGNORM, DINIT, DSTNRM, DTINIT, DTOL, DTYPE,
     1        D0INIT, F, F0, GTSTEP, INCFAC, IRC, KAGQT, LMAT, LMAX0,
     2        MODE, MODEL, MXFCAL, MXITER, NEXTV, NFCALL, NFGCAL,
     3        NGCALL, NITER, RADFAC, RADINC, RADIUS, RAD0, STEP, STGLIM,
     4        STLSTG, STPPAR, TOOBIG, TUNER4, TUNER5, VNEED, W, XIRC, X0
C
C  ***  IV SUBSCRIPT VALUES  ***
C
C/6
      DATA CNVCOD/55/, DG/37/, DTOL/59/, DTYPE/16/, IRC/29/, KAGQT/33/,
     1     LMAT/42/, MODE/35/, MODEL/5/, MXFCAL/17/, MXITER/18/,
     2     NEXTV/47/, NFCALL/6/, NFGCAL/7/, NGCALL/30/, NITER/31/,
     3     RADINC/8/, STEP/40/, STGLIM/11/, STLSTG/41/, TOOBIG/2/,
     4     VNEED/4/, W/34/, XIRC/13/, X0/43/
C/7
C     PARAMETER (CNVCOD=55, DG=37, DTOL=59, DTYPE=16, IRC=29, KAGQT=33,
C    1           LMAT=42, MODE=35, MODEL=5, MXFCAL=17, MXITER=18,
C    2           NEXTV=47, NFCALL=6, NFGCAL=7, NGCALL=30, NITER=31,
C    3           RADINC=8, STEP=40, STGLIM=11, STLSTG=41, TOOBIG=2,
C    4           VNEED=4, W=34, XIRC=13, X0=43)
C/
C
C  ***  V SUBSCRIPT VALUES  ***
C
C/6
      DATA DGNORM/1/, DINIT/38/, DSTNRM/2/, DTINIT/39/, D0INIT/40/,
     1     F/10/, F0/13/, GTSTEP/4/, INCFAC/23/, LMAX0/35/, RADFAC/16/,
     2     RADIUS/8/, RAD0/9/, STPPAR/5/, TUNER4/29/, TUNER5/30/
C/7
C     PARAMETER (DGNORM=1, DINIT=38, DSTNRM=2, DTINIT=39, D0INIT=40,
C    1           F=10, F0=13, GTSTEP=4, INCFAC=23, LMAX0=35, RADFAC=16,
C    2           RADIUS=8, RAD0=9, STPPAR=5, TUNER4=29, TUNER5=30)
C/
C
      DATA ONE/1.D+0/, ZERO/0.D+0/
C
C+++++++++++++++++++++++++++++++  BODY  ++++++++++++++++++++++++++++++++
C
      I = IV(1)
      IF (I .EQ. 1) GO TO 30
      IF (I .EQ. 2) GO TO 40
C
C  ***  CHECK VALIDITY OF IV AND V INPUT VALUES  ***
C
      IF (IV(1) .EQ. 0) CALL DDEFLT(2, IV, LIV, LV, V)
      IV(VNEED) = IV(VNEED) + N*(N+21)/2 + 7
      CALL DPARCK(2, D, IV, LIV, LV, N, V)
      I = IV(1) - 2
      IF (I .GT. 12) GO TO 999
      NN1O2 = N * (N + 1) / 2
      IF (LH .GE. NN1O2) GO TO (200,200,200,200,200,200,160,120,160,
     1                          10,10,20), I
         IV(1) = 66
         GO TO 330
C
C  ***  STORAGE ALLOCATION  ***
C
 10   IV(DTOL) = IV(LMAT) + NN1O2
      IV(X0) = IV(DTOL) + 2*N
      IV(STEP) = IV(X0) + N
      IV(STLSTG) = IV(STEP) + N
      IV(DG) = IV(STLSTG) + N
      IV(W) = IV(DG) + N
      IV(NEXTV) = IV(W) + 4*N + 7
      IF (IV(1) .NE. 13) GO TO 20
         IV(1) = 14
         GO TO 999
C
C  ***  INITIALIZATION  ***
C
 20   IV(NITER) = 0
      IV(NFCALL) = 1
      IV(NGCALL) = 1
      IV(NFGCAL) = 1
      IV(MODE) = -1
      IV(MODEL) = 1
      IV(STGLIM) = 1
      IV(TOOBIG) = 0
      IV(CNVCOD) = 0
      IV(RADINC) = 0
      V(RAD0) = ZERO
      V(STPPAR) = ZERO
      IF (V(DINIT) .GE. ZERO) CALL DVSCPY(N, D, V(DINIT))
      K = IV(DTOL)
      IF (V(DTINIT) .GT. ZERO) CALL DVSCPY(N, V(K), V(DTINIT))
      K = K + N
      IF (V(D0INIT) .GT. ZERO) CALL DVSCPY(N, V(K), V(D0INIT))
      IV(1) = 1
      GO TO 999
C
 30   V(F) = FX
      IF (IV(MODE) .GE. 0) GO TO 200
      IV(1) = 2
      IF (IV(TOOBIG) .EQ. 0) GO TO 999
         IV(1) = 63
         GO TO 330
C
C  ***  MAKE SURE GRADIENT COULD BE COMPUTED  ***
C
 40   IF (IV(NFGCAL) .NE. 0) GO TO 50
         IV(1) = 65
         GO TO 330
C
C  ***  UPDATE THE SCALE VECTOR D  ***
C
 50   DG1 = IV(DG)
      IF (IV(DTYPE) .LE. 0) GO TO 70
      K = DG1
      J = 0
      DO 60 I = 1, N
         J = J + I
         V(K) = H(J)
         K = K + 1
 60      CONTINUE
      CALL DDUPDU(D, V(DG1), IV, LIV, LV, N, V)
C
C  ***  COMPUTE SCALED GRADIENT AND ITS NORM  ***
C
 70   DG1 = IV(DG)
      K = DG1
      DO 80 I = 1, N
         V(K) = G(I) / D(I)
         K = K + 1
 80      CONTINUE
      V(DGNORM) = DNRM2(N, V(DG1),1)
C
C  ***  COMPUTE SCALED HESSIAN  ***
C
      K = 1
      DO 100 I = 1, N
         T = ONE / D(I)
         DO 90 J = 1, I
              H(K) = T * H(K) / D(J)
              K = K + 1
 90           CONTINUE
 100     CONTINUE
C
      IF (IV(CNVCOD) .NE. 0) GO TO 320
      IF (IV(MODE) .EQ. 0) GO TO 280
C
C  ***  ALLOW FIRST STEP TO HAVE SCALED 2-NORM AT MOST V(LMAX0)  ***
C
      V(RADIUS) = V(LMAX0)
C
      IV(MODE) = 0
C
C
C-----------------------------  MAIN LOOP  -----------------------------
C
C
C  ***  PRINT ITERATION SUMMARY, CHECK ITERATION LIMIT  ***
C
 110  CALL DITSUM(D, G, IV, LIV, LV, N, V, X)
 120  K = IV(NITER)
      IF (K .LT. IV(MXITER)) GO TO 130
         IV(1) = 10
         GO TO 330
C
 130  IV(NITER) = K + 1
C
C  ***  INITIALIZE FOR START OF NEXT ITERATION  ***
C
      DG1 = IV(DG)
      X01 = IV(X0)
      V(F0) = V(F)
      IV(IRC) = 4
      IV(KAGQT) = -1
C
C     ***  COPY X TO X0  ***
C
      CALL DCOPYGAMS(N, X,1,V(X01),1)
C
C  ***  UPDATE RADIUS  ***
C
      IF (K .EQ. 0) GO TO 150
      STEP1 = IV(STEP)
      K = STEP1
      DO 140 I = 1, N
         V(K) = D(I) * V(K)
         K = K + 1
 140     CONTINUE
      V(RADIUS) = V(RADFAC) * DNRM2(N, V(STEP1),1)
C
C  ***  CHECK DSTOPX AND FUNCTION EVALUATION LIMIT  ***
C
 150  IF (.NOT. DSTOPX(DUMMY)) GO TO 170
         IV(1) = 11
         GO TO 180
C
C     ***  COME HERE WHEN RESTARTING AFTER FUNC. EVAL. LIMIT OR DSTOPX.
C
 160  IF (V(F) .GE. V(F0)) GO TO 170
         V(RADFAC) = ONE
         K = IV(NITER)
         GO TO 130
C
 170  IF (IV(NFCALL) .LT. IV(MXFCAL)) GO TO 190
         IV(1) = 9
 180     IF (V(F) .GE. V(F0)) GO TO 330
C
C        ***  IN CASE OF DSTOPX OR FUNCTION EVALUATION LIMIT WITH
C        ***  IMPROVED V(F), EVALUATE THE GRADIENT AT X.
C
              IV(CNVCOD) = IV(1)
              GO TO 270
C
C. . . . . . . . . . . . .  COMPUTE CANDIDATE STEP  . . . . . . . . . .
C
 190  STEP1 = IV(STEP)
      DG1 = IV(DG)
      L = IV(LMAT)
      W1 = IV(W)
      CALL DGQTST(D, V(DG1), H, IV(KAGQT), V(L), N, V(STEP1), V, V(W1))
      IF (IV(IRC) .EQ. 6) GO TO 200
C
C  ***  COMPUTE F(X0 + STEP)  ***
C
      X01 = IV(X0)
      STEP1 = IV(STEP)
      CALL DVAXPY(N, X, ONE, V(STEP1), V(X01))
      IV(NFCALL) = IV(NFCALL) + 1
      IV(1) = 1
      IV(TOOBIG) = 0
      GO TO 999
C
C. . . . . . . . . . . . .  ASSESS CANDIDATE STEP  . . . . . . . . . . .
C
 200  STEP1 = IV(STEP)
      LSTGST = IV(STLSTG)
      X01 = IV(X0)
      CALL DASSST(D, IV, N, V(STEP1), V(LSTGST), V, X, V(X01))
C
      K = IV(IRC)
      GO TO (210,240,240,240,210,220,230,230,230,230,230,230,310,280), K
C
C     ***  RECOMPUTE STEP WITH NEW RADIUS  ***
C
 210     V(RADIUS) = V(RADFAC) * V(DSTNRM)
         GO TO 150
C
C  ***  COMPUTE STEP OF LENGTH V(LMAX0) FOR SINGULAR CONVERGENCE TEST.
C
 220  V(RADIUS) = V(LMAX0)
      GO TO 190
C
C  ***  CONVERGENCE OR FALSE CONVERGENCE  ***
C
 230  IV(CNVCOD) = K - 4
      IF (V(F) .GE. V(F0)) GO TO 320
         IF (IV(XIRC) .EQ. 14) GO TO 320
              IV(XIRC) = 14
C
C. . . . . . . . . . . .  PROCESS ACCEPTABLE STEP  . . . . . . . . . . .
C
 240  IF (IV(IRC) .NE. 3) GO TO 270
         TEMP1 = LSTGST
C
C     ***  PREPARE FOR GRADIENT TESTS  ***
C     ***  SET  TEMP1 = HESSIAN * STEP + G(X0)
C     ***             = DIAG(D) * (H * STEP + G(X0))
C
C        USE X0 VECTOR AS TEMPORARY.
         K = X01
         DO 250 I = 1, N
              V(K) = D(I) * V(STEP1)
              K = K + 1
              STEP1 = STEP1 + 1
 250          CONTINUE
         CALL DSLVMU(N, V(TEMP1), H, V(X01))
         DO 260 I = 1, N
              V(TEMP1) = D(I) * V(TEMP1) + G(I)
              TEMP1 = TEMP1 + 1
 260          CONTINUE
C
C  ***  COMPUTE GRADIENT AND HESSIAN  ***
C
 270  IV(NGCALL) = IV(NGCALL) + 1
      IV(1) = 2
      GO TO 999
C
 280  IV(1) = 2
      IF (IV(IRC) .NE. 3) GO TO 110
C
C  ***  SET V(RADFAC) BY GRADIENT TESTS  ***
C
      TEMP1 = IV(STLSTG)
      STEP1 = IV(STEP)
C
C     ***  SET  TEMP1 = DIAG(D)**-1 * (HESSIAN*STEP + (G(X0)-G(X)))  ***
C
      K = TEMP1
      DO 290 I = 1, N
         V(K) = (V(K) - G(I)) / D(I)
         K = K + 1
 290     CONTINUE
C
C     ***  DO GRADIENT TESTS  ***
C
      IF (DNRM2(N, V(TEMP1),1) .LE. V(DGNORM) * V(TUNER4)) GO TO 300
           IF (DDOT(N, G,1,V(STEP1),1)
     1               .GE. V(GTSTEP) * V(TUNER5))  GO TO 110
 300            V(RADFAC) = V(INCFAC)
                GO TO 110
C
C. . . . . . . . . . . . . .  MISC. DETAILS  . . . . . . . . . . . . . .
C
C  ***  BAD PARAMETERS TO ASSESS  ***
C
 310  IV(1) = 64
C
C  ***  PRINT SUMMARY OF FINAL ITERATION AND OTHER REQUESTED ITEMS  ***
C
 320  IV(1) = IV(CNVCOD)
      IV(CNVCOD) = 0
 330  CALL DITSUM(D, G, IV, LIV, LV, N, V, X)
C
 999  RETURN
C
C  ***  LAST CARD OF DHUMIT FOLLOWS  ***
      END
      SUBROUTINE DITSUM(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)
      DOUBLE PRECISION D(P), G(P), V(LV), X(P)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER ALG, I, IV1, M, NF, NG, OL, PU
C/6
      REAL MODEL1(6), MODEL2(6)
C/7
C     CHARACTER*4 MODEL1(6), MODEL2(6)
C/
      DOUBLE PRECISION NRELDF, OLDF, PRELDF, RELDF, ZERO
C
C  ***  INTRINSIC FUNCTIONS  ***
C/+
      INTEGER IABS
      DOUBLE PRECISION DABS, DMAX1
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
C/6
      DATA 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/7
C     PARAMETER (ALGSAV=51, NEEDHD=36, NFCALL=6, NFCOV=52, NGCALL=30,
C    1           NGCOV=53, NITER=31, OUTLEV=19, PRNTIT=39, PRUNIT=21,
C    2           SOLPRT=22, STATPR=23, SUSED=64, X0PRT=24)
C/
C
C  ***  V SUBSCRIPT VALUES  ***
C
C/6
      DATA DSTNRM/2/, F/10/, F0/13/, FDIF/11/, NREDUC/6/, PREDUC/7/,
     1     RELDX/17/, STPPAR/5/
C/7
C     PARAMETER (DSTNRM=2, F=10, F0=13, FDIF=11, NREDUC=6, PREDUC=7,
C    1           RELDX=17, STPPAR=5)
C/
C
C/6
      DATA ZERO/0.D+0/
C/7
C     PARAMETER (ZERO=0.D+0)
C/
C/6
      DATA MODEL1(1)/4H    /, MODEL1(2)/4H    /, MODEL1(3)/4H    /,
     1     MODEL1(4)/4H    /, MODEL1(5)/4H  G /, MODEL1(6)/4H  S /,
     2     MODEL2(1)/4H G  /, MODEL2(2)/4H S  /, MODEL2(3)/4HG-S /,
     3     MODEL2(4)/4HS-G /, MODEL2(5)/4H-S-G/, MODEL2(6)/4H-G-S/
C/7
C     DATA MODEL1/'    ','    ','    ','    ','  G ','  S '/,
C    1     MODEL2/' G  ',' S  ','G-S ','S-G ','-S-G','-G-S'/
C/
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 = IV(ALGSAV)
      IF (IV1 .LT. 2 .OR. IV1 .GT. 15) GO TO 370
      IF (OL .EQ. 0) GO TO 120
      IF (IV1 .GE. 12) GO TO 120
      IF (IV1 .EQ. 2 .AND. IV(NITER) .EQ. 0) GO TO 390
      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 = DMAX1(DABS(V(F0)), DABS(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,D10.3,2D9.2,D8.1,A3,A4,2D8.1,D9.2)
 110  FORMAT(I6,I5,D11.3,2D10.2,3D9.1,D10.2)
C
 120  IF (IV(STATPR) .LT. 0) GO TO 430
      GO TO (999, 999, 130, 150, 170, 190, 210, 230, 250, 270, 290, 310,
     1       330, 350, 520), 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 480
      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 DITSUM  ***
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,D17.6,D14.3))
      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) V(F)
      IF (ALG .EQ. 2) WRITE(PU,420) V(F)
 410  FORMAT(/11H     0    1,D10.3)
C365  FORMAT(/11H     0    1,E11.3)
 420  FORMAT(/11H     0    1,D11.3)
      GO TO 999
C
C  ***  PRINT VARIOUS INFORMATION REQUESTED ON SOLUTION  ***
C
 430  IV(NEEDHD) = 1
      IF (IV(STATPR) .EQ. 0) GO TO 480
         OLDF = DMAX1(DABS(V(F0)), DABS(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,D17.6,8H   RELDX,D17.3/12H FUNC. EVALS,
     1   I8,9X,11HGRAD. EVALS,I8/7H PRELDF,D16.3,6X,7HNPRELDF,D15.3)
C
         IF (IV(NFCOV) .GT. 0) WRITE(PU,460) IV(NFCOV)
 460     FORMAT(/1X,I4,50H EXTRA FUNC. EVALS FOR COVARIANCE AND DIAGNOST
     1ICS.)
         IF (IV(NGCOV) .GT. 0) WRITE(PU,470) IV(NGCOV)
 470     FORMAT(1X,I4,50H EXTRA GRAD. EVALS FOR COVARIANCE AND DIAGNOSTI
     1CS.)
C
 480  IF (IV(SOLPRT) .EQ. 0) GO TO 999
         IV(NEEDHD) = 1
         WRITE(PU,490)
 490  FORMAT(/22H     I      FINAL X(I),8X,4HD(I),10X,4HG(I)/)
         DO 500 I = 1, P
              WRITE(PU,510) I, X(I), D(I), G(I)
 500          CONTINUE
 510     FORMAT(1X,I5,D16.6,2D14.3)
      GO TO 999
C
 520  WRITE(PU,530)
 530  FORMAT(/24H INCONSISTENT DIMENSIONS)
 999  RETURN
C  ***  LAST CARD OF DITSUM FOLLOWS  ***
      END
      DOUBLE PRECISION FUNCTION DNRM2(N,DX,INCX)
C***BEGIN PROLOGUE  DNRM2
C***DATE WRITTEN   791001   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  D1A3B
C***KEYWORDS  BLAS,DOUBLE PRECISION,EUCLIDEAN,L2,LENGTH,LINEAR ALGEBRA,
C             NORM,VECTOR
C***AUTHOR  LAWSON, C. L., (JPL)
C           HANSON, R. J., (SNLA)
C           KINCAID, D. R., (U. OF TEXAS)
C           KROGH, F. T., (JPL)
C***PURPOSE  Euclidean length (L2 norm) of d.p. vector
C***DESCRIPTION
C
C                B L A S  Subprogram
C    Description of parameters
C
C     --Input--
C        N  number of elements in input vector(s)
C       DX  double precision vector with N elements
C     INCX  storage spacing between elements of DX
C
C     --Output--
C    DNRM2  double precision result (zero if N .LE. 0)
C
C     Euclidean norm of the N-vector stored in DX() with storage
C     increment INCX .
C     If    N .LE. 0 return with result = 0.
C     If N .GE. 1 then INCX must be .GE. 1
C
C           C.L. Lawson, 1978 Jan 08
C
C     Four phase method     using two built-in constants that are
C     hopefully applicable to all machines.
C         CUTLO = maximum of  DSQRT(U/EPS)  over all known machines.
C         CUTHI = minimum of  DSQRT(V)      over all known machines.
C     where
C         EPS = smallest no. such that EPS + 1. .GT. 1.
C         U   = smallest positive no.   (underflow limit)
C         V   = largest  no.            (overflow  limit)
C
C     Brief outline of algorithm..
C
C     Phase 1    scans zero components.
C     move to phase 2 when a component is nonzero and .LE. CUTLO
C     move to phase 3 when a component is .GT. CUTLO
C     move to phase 4 when a component is .GE. CUTHI/M
C     where M = N for X() real and M = 2*N for complex.
C
C     Values for CUTLO and CUTHI..
C     From the environmental parameters listed in the IMSL converter
C     document the limiting values are as followS..
C     CUTLO, S.P.   U/EPS = 2**(-102) for  Honeywell.  Close seconds are
C                   Univac and DEC at 2**(-103)
C                   Thus CUTLO = 2**(-51) = 4.44089E-16
C     CUTHI, S.P.   V = 2**127 for Univac, Honeywell, and DEC.
C                   Thus CUTHI = 2**(63.5) = 1.30438E19
C     CUTLO, D.P.   U/EPS = 2**(-67) for Honeywell and DEC.
C                   Thus CUTLO = 2**(-33.5) = 8.23181D-11
C     CUTHI, D.P.   same as S.P.  CUTHI = 1.30438D19
C     DATA CUTLO, CUTHI / 8.232D-11,  1.304D19 /
C     DATA CUTLO, CUTHI / 4.441E-16,  1.304E19 /
C***REFERENCES  LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T.,
C                 *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*,
C                 ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL
C                 SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323
C***ROUTINES CALLED  (NONE)
C***END PROLOGUE  DNRM2
      INTEGER          NEXT,N,INCX
      DOUBLE PRECISION   DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE
      DATA   ZERO, ONE /0.0D0, 1.0D0/
C
      DATA CUTLO, CUTHI / 8.232D-11,  1.304D19 /
C***FIRST EXECUTABLE STATEMENT  DNRM2
      IF(N .GT. 0) GO TO 10
         DNRM2  = ZERO
         GO TO 300
C
   10 ASSIGN 30 TO NEXT
      SUM = ZERO
      NN = N * INCX
C                                                 BEGIN MAIN LOOP
      I = 1
   20    GO TO NEXT,(30, 50, 70, 110)
   30 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85
      ASSIGN 50 TO NEXT
      XMAX = ZERO
C
C                        PHASE 1.  SUM IS ZERO
C
   50 IF( DX(I) .EQ. ZERO) GO TO 200
      IF( DABS(DX(I)) .GT. CUTLO) GO TO 85
C
C                                PREPARE FOR PHASE 2.
      ASSIGN 70 TO NEXT
      GO TO 105
C
C                                PREPARE FOR PHASE 4.
C
  100 I = J
      ASSIGN 110 TO NEXT
      SUM = (SUM / DX(I)) / DX(I)
  105 XMAX = DABS(DX(I))
      GO TO 115
C
C                   PHASE 2.  SUM IS SMALL.
C                             SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
C
   70 IF( DABS(DX(I)) .GT. CUTLO ) GO TO 75
C
C                     COMMON CODE FOR PHASES 2 AND 4.
C                     IN PHASE 4 SUM IS LARGE.  SCALE TO AVOID OVERFLOW.
C
  110 IF( DABS(DX(I)) .LE. XMAX ) GO TO 115
         SUM = ONE + SUM * (XMAX / DX(I))**2
         XMAX = DABS(DX(I))
         GO TO 200
C
  115 SUM = SUM + (DX(I)/XMAX)**2
      GO TO 200
C
C
C                  PREPARE FOR PHASE 3.
C
   75 SUM = (SUM * XMAX) * XMAX
C
C
C     FOR REAL OR D.P. SET HITEST = CUTHI/N
C     FOR COMPLEX      SET HITEST = CUTHI/(2*N)
C
   85 HITEST = CUTHI/FLOAT( N )
C
C                   PHASE 3.  SUM IS MID-RANGE.  NO SCALING.
C
      DO 95 J =I,NN,INCX
      IF(DABS(DX(J)) .GE. HITEST) GO TO 100
   95    SUM = SUM + DX(J)**2
      DNRM2 = DSQRT( SUM )
      GO TO 300
C
  200 CONTINUE
      I = I + INCX
      IF ( I .LE. NN ) GO TO 20
C
C              END OF MAIN LOOP.
C
C              COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
C
      DNRM2 = XMAX * DSQRT(SUM)
  300 CONTINUE
      RETURN
      END
      SUBROUTINE DPARCK(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)
      DOUBLE PRECISION D(N), V(LV)
C
      EXTERNAL  DVDFLT
      DOUBLE PRECISION D1MACH
C DVDFLT  -- SUPPLIES DEFAULT PARAMETER VALUES TO V ALONE.
C/+
      INTEGER MAX0
C/
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER I, II, IV1, J, K, L, M, MIV1, MIV2, NDFALT, PARSV1, PU
      INTEGER IJMP, JLIM(2), MINIV(2), NDFLT(2)
C/6
      INTEGER VARNM(2), SH(2)
      REAL CNGD(3), DFLT(3), VN(2,34), WHICH(3)
C/7
C     CHARACTER*1 VARNM(2), SH(2)
C     CHARACTER*4 CNGD(3), DFLT(3), VN(2,34), WHICH(3)
C/
      DOUBLE PRECISION 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
C/6
      DATA 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/
C/7
C     PARAMETER (ALGSAV=51, DINIT=38, DTYPE=16, DTYPE0=54, EPSLON=19,
C    1           INITS=25, IVNEED=3, LASTIV=44, LASTV=45, LMAT=42,
C    2           NEXTIV=46, NEXTV=47, NVDFLT=50, OLDN=38, PARPRT=20,
C    3           PARSAV=49, PERM=58, PRUNIT=21, VNEED=4)
C     SAVE BIG, MACHEP, TINY
C/
C
      DATA BIG/0.D+0/, MACHEP/-1.D+0/, TINY/1.D+0/, ZERO/0.D+0/
C/6
      DATA VN(1,1),VN(2,1)/4HEPSL,4HON../
      DATA VN(1,2),VN(2,2)/4HPHMN,4HFC../
      DATA VN(1,3),VN(2,3)/4HPHMX,4HFC../
      DATA VN(1,4),VN(2,4)/4HDECF,4HAC../
      DATA VN(1,5),VN(2,5)/4HINCF,4HAC../
      DATA VN(1,6),VN(2,6)/4HRDFC,4HMN../
      DATA VN(1,7),VN(2,7)/4HRDFC,4HMX../
      DATA VN(1,8),VN(2,8)/4HTUNE,4HR1../
      DATA VN(1,9),VN(2,9)/4HTUNE,4HR2../
      DATA VN(1,10),VN(2,10)/4HTUNE,4HR3../
      DATA VN(1,11),VN(2,11)/4HTUNE,4HR4../
      DATA VN(1,12),VN(2,12)/4HTUNE,4HR5../
      DATA VN(1,13),VN(2,13)/4HAFCT,4HOL../
      DATA VN(1,14),VN(2,14)/4HRFCT,4HOL../
      DATA VN(1,15),VN(2,15)/4HXCTO,4HL.../
      DATA VN(1,16),VN(2,16)/4HXFTO,4HL.../
      DATA VN(1,17),VN(2,17)/4HLMAX,4H0.../
      DATA VN(1,18),VN(2,18)/4HLMAX,4HS.../
      DATA VN(1,19),VN(2,19)/4HSCTO,4HL.../
      DATA VN(1,20),VN(2,20)/4HDINI,4HT.../
      DATA VN(1,21),VN(2,21)/4HDTIN,4HIT../
      DATA VN(1,22),VN(2,22)/4HD0IN,4HIT../
      DATA VN(1,23),VN(2,23)/4HDFAC,4H..../
      DATA VN(1,24),VN(2,24)/4HDLTF,4HDC../
      DATA VN(1,25),VN(2,25)/4HDLTF,4HDJ../
      DATA VN(1,26),VN(2,26)/4HDELT,4HA0../
      DATA VN(1,27),VN(2,27)/4HFUZZ,4H..../
      DATA VN(1,28),VN(2,28)/4HRLIM,4HIT../
      DATA VN(1,29),VN(2,29)/4HCOSM,4HIN../
      DATA VN(1,30),VN(2,30)/4HHUBE,4HRC../
      DATA VN(1,31),VN(2,31)/4HRSPT,4HOL../
      DATA VN(1,32),VN(2,32)/4HSIGM,4HIN../
      DATA VN(1,33),VN(2,33)/4HETA0,4H..../
      DATA VN(1,34),VN(2,34)/4HBIAS,4H..../
C/7
C     DATA VN(1,1),VN(2,1)/'EPSL','ON..'/
C     DATA VN(1,2),VN(2,2)/'PHMN','FC..'/
C     DATA VN(1,3),VN(2,3)/'PHMX','FC..'/
C     DATA VN(1,4),VN(2,4)/'DECF','AC..'/
C     DATA VN(1,5),VN(2,5)/'INCF','AC..'/
C     DATA VN(1,6),VN(2,6)/'RDFC','MN..'/
C     DATA VN(1,7),VN(2,7)/'RDFC','MX..'/
C     DATA VN(1,8),VN(2,8)/'TUNE','R1..'/
C     DATA VN(1,9),VN(2,9)/'TUNE','R2..'/
C     DATA VN(1,10),VN(2,10)/'TUNE','R3..'/
C     DATA VN(1,11),VN(2,11)/'TUNE','R4..'/
C     DATA VN(1,12),VN(2,12)/'TUNE','R5..'/
C     DATA VN(1,13),VN(2,13)/'AFCT','OL..'/
C     DATA VN(1,14),VN(2,14)/'RFCT','OL..'/
C     DATA VN(1,15),VN(2,15)/'XCTO','L...'/
C     DATA VN(1,16),VN(2,16)/'XFTO','L...'/
C     DATA VN(1,17),VN(2,17)/'LMAX','0...'/
C     DATA VN(1,18),VN(2,18)/'LMAX','S...'/
C     DATA VN(1,19),VN(2,19)/'SCTO','L...'/
C     DATA VN(1,20),VN(2,20)/'DINI','T...'/
C     DATA VN(1,21),VN(2,21)/'DTIN','IT..'/
C     DATA VN(1,22),VN(2,22)/'D0IN','IT..'/
C     DATA VN(1,23),VN(2,23)/'DFAC','....'/
C     DATA VN(1,24),VN(2,24)/'DLTF','DC..'/
C     DATA VN(1,25),VN(2,25)/'DLTF','DJ..'/
C     DATA VN(1,26),VN(2,26)/'DELT','A0..'/
C     DATA VN(1,27),VN(2,27)/'FUZZ','....'/
C     DATA VN(1,28),VN(2,28)/'RLIM','IT..'/
C     DATA VN(1,29),VN(2,29)/'COSM','IN..'/
C     DATA VN(1,30),VN(2,30)/'HUBE','RC..'/
C     DATA VN(1,31),VN(2,31)/'RSPT','OL..'/
C     DATA VN(1,32),VN(2,32)/'SIGM','IN..'/
C     DATA VN(1,33),VN(2,33)/'ETA0','....'/
C     DATA VN(1,34),VN(2,34)/'BIAS','....'/
C/
C
      DATA VM(1)/1.0D-3/, VM(2)/-0.99D+0/, VM(3)/1.0D-3/, VM(4)/1.0D-2/,
     1     VM(5)/1.2D+0/, VM(6)/1.D-2/, VM(7)/1.2D+0/, VM(8)/0.D+0/,
     2     VM(9)/0.D+0/, VM(10)/1.D-3/, VM(11)/-1.D+0/, VM(15)/0.D+0/,
     3     VM(16)/0.D+0/, VM(19)/0.D+0/, VM(20)/-10.D+0/, VM(21)/0.D+0/,
     4     VM(22)/0.D+0/, VM(23)/0.D+0/, VM(27)/1.01D+0/,
     5     VM(28)/1.D+10/, VM(30)/0.D+0/, VM(31)/0.D+0/, VM(32)/0.D+0/,
     6     VM(34)/0.D+0/
      DATA VX(1)/0.9D+0/, VX(2)/-1.D-3/, VX(3)/1.D+1/, VX(4)/0.8D+0/,
     1     VX(5)/1.D+2/, VX(6)/0.8D+0/, VX(7)/1.D+2/, VX(8)/0.5D+0/,
     2     VX(9)/0.5D+0/, VX(10)/1.D+0/, VX(11)/1.D+0/, VX(14)/0.1D+0/,
     3     VX(15)/1.D+0/, VX(16)/1.D+0/, VX(19)/1.D+0/, VX(23)/1.D+0/,
     4     VX(24)/1.D+0/, VX(25)/1.D+0/, VX(26)/1.D+0/, VX(27)/1.D+10/,
     5     VX(29)/1.D+0/, VX(31)/1.D+0/, VX(32)/1.D+0/, VX(33)/1.D+0/,
     6     VX(34)/1.D+0/
C
C/6
      DATA VARNM(1)/1HP/, VARNM(2)/1HN/, SH(1)/1HS/, SH(2)/1HH/
      DATA CNGD(1),CNGD(2),CNGD(3)/4H---C,4HHANG,4HED V/,
     1     DFLT(1),DFLT(2),DFLT(3)/4HNOND,4HEFAU,4HLT V/
C/7
C     DATA VARNM(1)/'P'/, VARNM(2)/'N'/, SH(1)/'S'/, SH(2)/'H'/
C     DATA CNGD(1),CNGD(2),CNGD(3)/'---C','HANG','ED V'/,
C    1     DFLT(1),DFLT(2),DFLT(3)/'NOND','EFAU','LT V'/
C/
      DATA IJMP/33/, JLIM(1)/0/, JLIM(2)/24/, NDFLT(1)/32/, NDFLT(2)/25/
      DATA MINIV(1)/80/, MINIV(2)/59/
C
C...............................  BODY  ................................
C
      IF (ALG .LT. 1 .OR. ALG .GT. 2) GO TO 330
      IF (IV(1) .EQ. 0) CALL DDEFLT(ALG, IV, LIV, LV, V)
      PU = IV(PRUNIT)
      MIV1 = MINIV(ALG)
      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 290
      IV(IVNEED) = 0
      IV(LASTV) = MAX0(IV(VNEED), 0) + IV(LMAT) - 1
      IF (LIV .LT. MIV2) GO TO 290
      IF (LV .LT. IV(LASTV)) GO TO 310
      IV(VNEED) = 0
      IF (ALG .EQ. IV(ALGSAV)) GO TO 20
         IF (PU .NE. 0) WRITE(PU,10) ALG, IV(ALGSAV)
 10      FORMAT(/' THE FIRST PARAMETER TO DDEFLT SHOULD BE',I3,
     1          12H RATHER THAN,I3)
         IV(1) = 82
         GO TO 999
 20   IV1 = IV(1)
      IF (IV1 .LT. 12 .OR. IV1 .GT. 14) GO TO 50
         IF (N .GE. 1) GO TO 40
              IV(1) = 81
              IF (PU .EQ. 0) GO TO 999
              WRITE(PU,30) VARNM(ALG), N
 30           FORMAT(/8H /// BAD,A1,2H =,I5)
              GO TO 999
 40      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 DVDFLT(ALG, LV-K, V(K+1))
         IV(DTYPE0) = 2 - ALG
         IV(OLDN) = N
         WHICH(1) = DFLT(1)
         WHICH(2) = DFLT(2)
         WHICH(3) = DFLT(3)
         GO TO 100
 50   IF (N .EQ. IV(OLDN)) GO TO 70
         IV(1) = 17
         IF (PU .EQ. 0) GO TO 999
         WRITE(PU,60) VARNM(ALG), IV(OLDN), N
 60      FORMAT(/5H /// ,1A1,14H CHANGED FROM ,I5,4H TO ,I5)
         GO TO 999
C
 70   IF (IV1 .LE. 11 .AND. IV1 .GE. 1) GO TO 90
         IV(1) = 80
         IF (PU .NE. 0) WRITE(PU,80) IV1
 80      FORMAT(/13H ///  IV(1) =,I5,28H SHOULD BE BETWEEN 0 AND 14.)
         GO TO 999
C
 90   WHICH(1) = CNGD(1)
      WHICH(2) = CNGD(2)
      WHICH(3) = CNGD(3)
C
 100  IF (IV1 .EQ. 14) IV1 = 12
      IF (BIG .GT. TINY) GO TO 110
         TINY = D1MACH(1)
         MACHEP = D1MACH(4)
         BIG = D1MACH(2)
         VM(12) = MACHEP
         VX(12) = BIG
         VM(13) = TINY
         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) = DSQRT(D1MACH(2))*16.
         VM(29) = MACHEP
         VX(30) = BIG
         VM(33) = MACHEP
 110  M = 0
      I = 1
      J = JLIM(ALG)
      K = EPSLON
      NDFALT = NDFLT(ALG)
      DO 140 L = 1, NDFALT
         VK = V(K)
         IF (VK .GE. VM(I) .AND. VK .LE. VX(I)) GO TO 130
              M = K
              IF (PU .NE. 0) WRITE(PU,120) VN(1,I), VN(2,I), K, VK,
     1                                    VM(I), VX(I)
 120          FORMAT(/6H ///  ,2A4,5H.. V(,I2,3H) =,D11.3,7H SHOULD,
     1               11H BE BETWEEN,D11.3,4H AND,D11.3)
 130     K = K + 1
         I = I + 1
         IF (I .EQ. J) I = IJMP
 140     CONTINUE
C
      IF (IV(NVDFLT) .EQ. NDFALT) GO TO 160
         IV(1) = 51
         IF (PU .EQ. 0) GO TO 999
         WRITE(PU,150) IV(NVDFLT), NDFALT
 150     FORMAT(/13H IV(NVDFLT) =,I5,13H RATHER THAN ,I5)
         GO TO 999
 160  IF ((IV(DTYPE) .GT. 0 .OR. V(DINIT) .GT. ZERO) .AND. IV1 .EQ. 12)
     1                  GO TO 190
      DO 180 I = 1, N
         IF (D(I) .GT. ZERO) GO TO 180
              M = 18
              IF (PU .NE. 0) WRITE(PU,170) I, D(I)
 170     FORMAT(/8H ///  D(,I3,3H) =,D11.3,19H SHOULD BE POSITIVE)
 180     CONTINUE
 190  IF (M .EQ. 0) GO TO 200
         IV(1) = M
         GO TO 999
C
 200  IF (PU .EQ. 0 .OR. IV(PARPRT) .EQ. 0) GO TO 999
      IF (IV1 .NE. 12 .OR. IV(INITS) .EQ. ALG-1) GO TO 220
         M = 1
         WRITE(PU,210) SH(ALG), IV(INITS)
 210     FORMAT(/22H NONDEFAULT VALUES..../5H INIT,A1,14H..... IV(25) =,
     1          I3)
 220  IF (IV(DTYPE) .EQ. IV(DTYPE0)) GO TO 240
         IF (M .EQ. 0) WRITE(PU,250) WHICH
         M = 1
         WRITE(PU,230) IV(DTYPE)
 230     FORMAT(20H DTYPE..... IV(16) =,I3)
 240  I = 1
      J = JLIM(ALG)
      K = EPSLON
      L = IV(PARSAV)
      NDFALT = NDFLT(ALG)
      DO 280 II = 1, NDFALT
         IF (V(K) .EQ. V(L)) GO TO 270
              IF (M .EQ. 0) WRITE(PU,250) WHICH
 250          FORMAT(/1H ,3A4,9HALUES..../)
              M = 1
              WRITE(PU,260) VN(1,I), VN(2,I), K, V(K)
 260          FORMAT(1X,2A4,5H.. V(,I2,3H) =,D15.7)
 270     K = K + 1
         L = L + 1
         I = I + 1
         IF (I .EQ. J) I = IJMP
 280     CONTINUE
C
      IV(DTYPE0) = IV(DTYPE)
      PARSV1 = IV(PARSAV)
      CALL DCOPYGAMS(IV(NVDFLT), V(EPSLON),1,V(PARSV1),1)
      GO TO 999
C
 290  IV(1) = 15
      IF (PU .EQ. 0) GO TO 999
      WRITE(PU,300) LIV, MIV2
 300  FORMAT(/10H /// LIV =,I5,17H MUST BE AT LEAST,I5)
      IF (LIV .LT. MIV1) GO TO 999
      IF (LV .LT. IV(LASTV)) GO TO 310
      GO TO 999
C
 310  IV(1) = 16
      IF (PU .EQ. 0) GO TO 999
      WRITE(PU,320) LV, IV(LASTV)
 320  FORMAT(/9H /// LV =,I5,17H MUST BE AT LEAST,I5)
      GO TO 999
C
 330  IV(1) = 67
C
 999  RETURN
C  ***  LAST CARD OF DPARCK FOLLOWS  ***
      END
      SUBROUTINE DSLVMU(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
      DOUBLE PRECISION S(1), X(P), Y(P)
C     DIMENSION S(P*(P+1)/2)
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER I, IM1, J, K
      DOUBLE PRECISION XI
C
C  ***  NO INTRINSIC FUNCTIONS  ***
C
C  ***  EXTERNAL FUNCTION  ***
C
      DOUBLE PRECISION DDOT
C
C-----------------------------------------------------------------------
C
      J = 1
      DO 10 I = 1, P
         Y(I) = DDOT(I, S(J),1,X,1)
         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 CARD OF DSLVMU FOLLOWS  ***
      END
      LOGICAL FUNCTION DSTOPX(IDUMMY)
C     *****PARAMETERS...
      INTEGER IDUMMY
C
C     ..................................................................
C
C     *****PURPOSE...
C     THIS FUNCTION MAY SERVE AS THE DSTOPX (ASYNCHRONOUS INTERRUPTION)
C     FUNCTION FOR THE NL2SOL (NONLINEAR LEAST-SQUARES) PACKAGE AT
C     THOSE INSTALLATIONS WHICH DO NOT WISH TO IMPLEMENT A
C     DYNAMIC DSTOPX.
C
C     *****ALGORITHM NOTES...
C     AT INSTALLATIONS WHERE THE NL2SOL SYSTEM IS USED
C     INTERACTIVELY, THIS DUMMY DSTOPX 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 DSTOPX.
C
C     ..................................................................
C
      DSTOPX = .FALSE.
      RETURN
      END
      SUBROUTINE DVAXPY(P, W, A, X, Y)
C
C  ***  SET W = A*X + Y  --  W, X, Y = P-VECTORS, A = SCALAR  ***
C
      INTEGER P
      DOUBLE PRECISION 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
      SUBROUTINE DVDFLT(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
      DOUBLE PRECISION V(LV)
C/+
      DOUBLE PRECISION DMAX1
C/
      DOUBLE PRECISION D1MACH
C
      DOUBLE PRECISION 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
C/6
      DATA ONE/1.D+0/, THREE/3.D+0/
C/7
C     PARAMETER (ONE=1.D+0, THREE=3.D+0)
C/
C
C  ***  V SUBSCRIPT VALUES  ***
C
C/6
      DATA 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/7
C     PARAMETER (AFCTOL=31, BIAS=43, COSMIN=47, DECFAC=22, DELTA0=44,
C    1           DFAC=41, DINIT=38, DLTFDC=42, DLTFDJ=43, DTINIT=39,
C    2           D0INIT=40, EPSLON=19, ETA0=42, FUZZ=45, HUBERC=48,
C    3           INCFAC=23, LMAX0=35, LMAXS=36, PHMNFC=20, PHMXFC=21,
C    4           RDFCMN=24, RDFCMX=25, RFCTOL=32, RLIMIT=46, RSPTOL=49,
C    5           SCTOL=37, SIGMIN=50, TUNER1=26, TUNER2=27, TUNER3=28,
C    6           TUNER4=29, TUNER5=30, XCTOL=33, XFTOL=34)
C/
C
C-------------------------------  BODY  --------------------------------
C
      MACHEP = D1MACH(4)
      V(AFCTOL) = 1.D-20
      IF (MACHEP .GT. 1.D-10) V(AFCTOL) = MACHEP**2
      V(DECFAC) = 0.5D+0
      SQTEPS = DSQRT(D1MACH(4))
      V(DFAC) = 0.6D+0
      V(DELTA0) = SQTEPS
      V(DTINIT) = 1.D-6
      MEPCRT = MACHEP ** (ONE/THREE)
      V(D0INIT) = 1.D+0
      V(EPSLON) = 0.1D+0
      V(INCFAC) = 2.D+0
      V(LMAX0) = 1.D+0
      V(LMAXS) = 1.D+0
      V(PHMNFC) = -0.1D+0
      V(PHMXFC) = 0.1D+0
      V(RDFCMN) = 0.1D+0
      V(RDFCMX) = 4.D+0
      V(RFCTOL) = DMAX1(1.D-10, MEPCRT**2)
      V(SCTOL) = V(RFCTOL)
      V(TUNER1) = 0.1D+0
      V(TUNER2) = 1.D-4
      V(TUNER3) = 0.75D+0
      V(TUNER4) = 0.5D+0
      V(TUNER5) = 0.75D+0
      V(XCTOL) = SQTEPS
      V(XFTOL) = 1.D+2 * MACHEP
C
      IF (ALG .GE. 2) GO TO 10
C
C  ***  REGRESSION  VALUES
C
      V(COSMIN) = DMAX1(1.D-6, 1.D+2 * MACHEP)
      V(DINIT) = 0.D+0
      V(DLTFDC) = MEPCRT
      V(DLTFDJ) = SQTEPS
      V(FUZZ) = 1.5D+0
      V(HUBERC) = 0.7D+0
      V(RLIMIT) = DSQRT(D1MACH(2))*16.
      V(RSPTOL) = 1.D-2
      V(SIGMIN) = 1.D-4
      GO TO 999
C
C  ***  GENERAL OPTIMIZATION VALUES
C
 10   V(BIAS) = 0.8D+0
      V(DINIT) = -1.0D+0
      V(ETA0) = 1.0D+3 * MACHEP
C
 999  RETURN
C  ***  LAST CARD OF DVDFLT FOLLOWS  ***
      END
      SUBROUTINE DVSCPY(P, Y, S)
C
C  ***  SET P-VECTOR Y TO SCALAR S  ***
C
      INTEGER P
      DOUBLE PRECISION S, Y(P)
C
      INTEGER I
C
      DO 10 I = 1, P
 10      Y(I) = S
      RETURN
      END
      SUBROUTINE DCOPYGAMS(N,DX,INCX,DY,INCY)
C***BEGIN PROLOGUE  DCOPYGAMS
C***DATE WRITTEN   791001   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  D1A5
C***KEYWORDS  BLAS,COPY,DOUBLE PRECISION,LINEAR ALGEBRA,VECTOR
C***AUTHOR  LAWSON, C. L., (JPL)
C           HANSON, R. J., (SNLA)
C           KINCAID, D. R., (U. OF TEXAS)
C           KROGH, F. T., (JPL)
C***PURPOSE  D.P. vector copy y = x
C***DESCRIPTION
C
C                B L A S  Subprogram
C    Description of Parameters
C
C     --Input--
C        N  number of elements in input vector(s)
C       DX  double precision vector with N elements
C     INCX  storage spacing between elements of DX
C       DY  double precision vector with N elements
C     INCY  storage spacing between elements of DY
C
C     --Output--
C       DY  copy of vector DX (unchanged if N .LE. 0)
C
C     Copy double precision DX to double precision DY.
C     For I = 0 to N-1, copy DX(LX+I*INCX) to DY(LY+I*INCY),
C     where LX = 1 if INCX .GE. 0, else LX = (-INCX)*N, and LY is
C     defined in a similar way using INCY.
C***REFERENCES  LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T.,
C                 *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*,
C                 ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL
C                 SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323
C***ROUTINES CALLED  (NONE)
C***END PROLOGUE  DCOPYGAMS
C
      DOUBLE PRECISION DX(1),DY(1)
C***FIRST EXECUTABLE STATEMENT  DCOPYGAMS
      IF(N.LE.0)RETURN
      IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60
    5 CONTINUE
C
C        CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS.
C
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        DY(IY) = DX(IX)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 7.
C
   20 M = MOD(N,7)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        DY(I) = DX(I)
   30 CONTINUE
      IF( N .LT. 7 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,7
        DY(I) = DX(I)
        DY(I + 1) = DX(I + 1)
        DY(I + 2) = DX(I + 2)
        DY(I + 3) = DX(I + 3)
        DY(I + 4) = DX(I + 4)
        DY(I + 5) = DX(I + 5)
        DY(I + 6) = DX(I + 6)
   50 CONTINUE
      RETURN
C
C        CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS.
C
   60 CONTINUE
      NS=N*INCX
          DO 70 I=1,NS,INCX
          DY(I) = DX(I)
   70     CONTINUE
      RETURN
      END
      DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)
C***BEGIN PROLOGUE  DDOT
C***DATE WRITTEN   791001   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  D1A4
C***KEYWORDS  BLAS,DOUBLE PRECISION,INNER PRODUCT,LINEAR ALGEBRA,VECTOR
C***AUTHOR  LAWSON, C. L., (JPL)
C           HANSON, R. J., (SNLA)
C           KINCAID, D. R., (U. OF TEXAS)
C           KROGH, F. T., (JPL)
C***PURPOSE  D.P. inner product of d.p. vectors
C***DESCRIPTION
C
C                B L A S  Subprogram
C    Description of Parameters
C
C     --Input--
C        N  number of elements in input vector(s)
C       DX  double precision vector with N elements
C     INCX  storage spacing between elements of DX
C       DY  double precision vector with N elements
C     INCY  storage spacing between elements of DY
C
C     --Output--
C     DDOT  double precision dot product (zero if N .LE. 0)
C
C     Returns the dot product of double precision DX and DY.
C     DDOT = sum for I = 0 to N-1 of  DX(LX+I*INCX) * DY(LY+I*INCY)
C     where LX = 1 if INCX .GE. 0, else LX = (-INCX)*N, and LY is
C     defined in a similar way using INCY.
C***REFERENCES  LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T.,
C                 *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*,
C                 ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL
C                 SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323
C***ROUTINES CALLED  (NONE)
C***END PROLOGUE  DDOT
C
      DOUBLE PRECISION DX(1),DY(1)
C***FIRST EXECUTABLE STATEMENT  DDOT
      DDOT = 0.D0
      IF(N.LE.0)RETURN
      IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60
    5 CONTINUE
C
C         CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS.
C
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
         DDOT = DDOT + DX(IX)*DY(IY)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1.
C
C
C        CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 5.
C
   20 M = MOD(N,5)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
         DDOT = DDOT + DX(I)*DY(I)
   30 CONTINUE
      IF( N .LT. 5 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,5
         DDOT = DDOT + DX(I)*DY(I) + DX(I+1)*DY(I+1) +
     1   DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4)
   50 CONTINUE
      RETURN
C
C         CODE FOR POSITIVE EQUAL INCREMENTS .NE.1.
C
   60 CONTINUE
      NS = N*INCX
          DO 70 I=1,NS,INCX
          DDOT = DDOT + DX(I)*DY(I)
   70     CONTINUE
      RETURN
      END
      SUBROUTINE DDUPDU(D, HDIAG, IV, LIV, LV, N, V)
C
C  ***  UPDATE SCALE VECTOR D FOR DHUMSL  ***
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER LIV, LV, N
      INTEGER IV(LIV)
      DOUBLE PRECISION D(N), HDIAG(N), V(LV)
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER DTOLI, D0I, I
      DOUBLE PRECISION T, VDFAC
C
C  ***  INTRINSIC FUNCTIONS  ***
C/+
      DOUBLE PRECISION DABS, DMAX1, DSQRT
C/
C  ***  SUBSCRIPTS FOR IV AND V  ***
C
      INTEGER DFAC, DTOL, DTYPE, NITER
C/6
      DATA DFAC/41/, DTOL/59/, DTYPE/16/, NITER/31/
C/7
C     PARAMETER (DFAC=41, DTOL=59, DTYPE=16, NITER=31)
C/
C
C-------------------------------  BODY  --------------------------------
C
      I = IV(DTYPE)
      IF (I .EQ. 1) GO TO 10
         IF (IV(NITER) .GT. 0) GO TO 999
C
 10   DTOLI = IV(DTOL)
      D0I = DTOLI + N
      VDFAC = V(DFAC)
      DO 20 I = 1, N
         T = DMAX1(DSQRT(DABS(HDIAG(I))), VDFAC*D(I))
         IF (T .LT. V(DTOLI)) T = DMAX1(V(DTOLI), V(D0I))
         D(I) = T
         DTOLI = DTOLI + 1
         D0I = D0I + 1
 20      CONTINUE
C
 999  RETURN
C  ***  LAST CARD OF DDUPDU FOLLOWS  ***
      END
      SUBROUTINE DGQTST(D, DIG, DIHDI, KA, L, P, STEP, V, W)
C
C  *** COMPUTE GOLDFELD-QUANDT-TROTTER STEP BY MORE-HEBDEN TECHNIQUE ***
C  ***  (NL2SOL VERSION 2.2), MODIFIED A LA MORE AND SORENSEN  ***
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER KA, P
      DOUBLE PRECISION D(P), DIG(P), DIHDI(1), L(1), V(21), STEP(P),
     1                 W(1)
C     DIMENSION DIHDI(P*(P+1)/2), L(P*(P+1)/2), W(4*P+7)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ***  PURPOSE  ***
C
C        GIVEN THE (COMPACTLY STORED) LOWER TRIANGLE OF A SCALED
C     HESSIAN (APPROXIMATION) AND A NONZERO SCALED GRADIENT VECTOR,
C     THIS SUBROUTINE COMPUTES A GOLDFELD-QUANDT-TROTTER STEP OF
C     APPROXIMATE LENGTH V(RADIUS) BY THE MORE-HEBDEN TECHNIQUE.  IN
C     OTHER WORDS, STEP IS COMPUTED TO (APPROXIMATELY) MINIMIZE
C     PSI(STEP) = (G**T)*STEP + 0.5*(STEP**T)*H*STEP  SUCH THAT THE
C     2-NORM OF D*STEP IS AT MOST (APPROXIMATELY) V(RADIUS), WHERE
C     G  IS THE GRADIENT,  H  IS THE HESSIAN, AND  D  IS A DIAGONAL
C     SCALE MATRIX WHOSE DIAGONAL IS STORED IN THE PARAMETER D.
C     (DGQTST 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 DLITVM - APPLIES INVERSE-TRANSPOSE OF COMPACT LOWER TRIANG. MATRIX.
C DLIVMU - APPLIES INVERSE OF COMPACT LOWER TRIANG. MATRIX.
C DLSQRT  - FINDS CHOLESKY FACTOR (OF COMPACTLY STORED LOWER TRIANG.).
C DLSVMN - RETURNS APPROX. TO MIN. SING. VALUE OF LOWER TRIANG. MATRIX.
C
C  ***  REFERENCES  ***
C
C 1.  DENNIS, J.E., GAY, D.M., AND WELSCH, R.E. (1981), AN ADAPTIVE
C             NONLINEAR LEAST-SQUARES ALGORITHM, ACM TRANS. MATH.
C             SOFTWARE, VOL. 7, NO. 3.
C 2.  GAY, D.M. (1981), COMPUTING OPTIMAL LOCALLY CONSTRAINED STEPS,
C             SIAM J. SCI. STATIST. COMPUTING, VOL. 2, NO. 2, PP.
C             186-197.
C 3.  GOLDFELD, S.M., QUANDT, R.E., AND TROTTER, H.F. (1966),
C             MAXIMIZATION BY QUADRATIC HILL-CLIMBING, ECONOMETRICA 34,
C             PP. 541-551.
C 4.  HEBDEN, M.D. (1973), AN ALGORITHM FOR MINIMIZATION USING EXACT
C             SECOND DERIVATIVES, REPORT T.P. 515, THEORETICAL PHYSICS
C             DIV., A.E.R.E. HARWELL, OXON., ENGLAND.
C 5.  MORE, J.J. (1978), THE LEVENBERG-MARQUARDT ALGORITHM, IMPLEMEN-
C             TATION AND THEORY, PP.105-116 OF SPRINGER LECTURE NOTES
C             IN MATHEMATICS NO. 630, EDITED BY G.A. WATSON, SPRINGER-
C             VERLAG, BERLIN AND NEW YORK.
C 6.  MORE, J.J., AND SORENSEN, D.C. (1981), COMPUTING A TRUST REGION
C             STEP, TECHNICAL REPORT ANL-81-83, ARGONNE NATIONAL LAB.
C 7.  VARGA, R.S. (1965), MINIMAL GERSCHGORIN SETS, PACIFIC J. MATH. 15,
C             PP. 719-729.
C
C  ***  GENERAL  ***
C
C     CODED BY DAVID M. GAY.
C     THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH
C     SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS
C     MCS-7600324, DCR75-10143, 76-14311DSS, MCS76-11989, AND
C     MCS-7906671.
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ***  LOCAL VARIABLES  ***
C
      LOGICAL RESTRT
      INTEGER DGGDMX, DIAG, DIAG0, DSTSAV, EMAX, EMIN, I, IM1, INC, IRC,
     1        J, K, KALIM, KAMIN, K1, LK0, PHIPIN, Q, Q0, UK0, X
      DOUBLE PRECISION ALPHAK, AKI, AKK, DELTA, DST, EPS, GTSTA, LK,
     1                 OLDPHI, PHI, PHIMAX, PHIMIN, PSIFAC, RAD, RADSQ,
     2                 ROOT, SI, SK, SW, T, TWOPSI, T1, T2, UK, WI
C
C     ***  CONSTANTS  ***
      DOUBLE PRECISION BIG, DGXFAC, EPSFAC, FOUR, HALF, KAPPA, NEGONE,
     1                 ONE, P001, SIX, THREE, TWO, ZERO
C
C  ***  INTRINSIC FUNCTIONS  ***
C/+
      DOUBLE PRECISION DABS, DMAX1, DMIN1, DSQRT
C/
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      EXTERNAL DLITVM, DLIVMU, DLSQRT, DLSVMN
      DOUBLE PRECISION DDOT, DLSVMN, D1MACH, DNRM2
C
C  ***  SUBSCRIPTS FOR V  ***
C
      INTEGER DGNORM, DSTNRM, DST0, EPSLON, GTSTEP, STPPAR, NREDUC,
     1        PHMNFC, PHMXFC, PREDUC, RADIUS, RAD0
C/6
      DATA 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/7
C     PARAMETER (DGNORM=1, DSTNRM=2, DST0=3, EPSLON=19, GTSTEP=4,
C    1           NREDUC=6, PHMNFC=20, PHMXFC=21, PREDUC=7, RADIUS=8,
C    2           RAD0=9, STPPAR=5)
C/
C
C/6
      DATA EPSFAC/50.0D+0/, FOUR/4.0D+0/, HALF/0.5D+0/,
     1     KAPPA/2.0D+0/, NEGONE/-1.0D+0/, ONE/1.0D+0/, P001/1.0D-3/,
     2     SIX/6.0D+0/, THREE/3.0D+0/, TWO/2.0D+0/, ZERO/0.0D+0/
C/7
C     PARAMETER (EPSFAC=50.0D+0, FOUR=4.0D+0, HALF=0.5D+0,
C    1     KAPPA=2.0D+0, NEGONE=-1.0D+0, ONE=1.0D+0, P001=1.0D-3,
C    2     SIX=6.0D+0, THREE=3.0D+0, TWO=2.0D+0, ZERO=0.0D+0)
C     SAVE DGXFAC
C/
      DATA BIG/0.D+0/, DGXFAC/0.D+0/
C
C  ***  BODY  ***
C
C     ***  STORE LARGEST ABS. ENTRY IN (D**-1)*H*(D**-1) AT W(DGGDMX).
      DGGDMX = P + 1
C     ***  STORE GERSCHGORIN OVER- AND UNDERESTIMATES OF THE LARGEST
C     ***  AND SMALLEST EIGENVALUES OF (D**-1)*H*(D**-1) AT W(EMAX)
C     ***  AND W(EMIN) RESPECTIVELY.
      EMAX = DGGDMX + 1
      EMIN = EMAX + 1
C     ***  FOR USE IN RECOMPUTING STEP, THE FINAL VALUES OF LK, UK, DST,
C     ***  AND THE INVERSE DERIVATIVE OF MORE*S PHI AT 0 (FOR POS. DEF.
C     ***  H) ARE STORED IN W(LK0), W(UK0), W(DSTSAV), AND W(PHIPIN)
C     ***  RESPECTIVELY.
      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 = TWO * V(EPSLON) / (THREE * (FOUR * (V(PHMNFC) + ONE) *
     1                       (KAPPA + ONE)  +  KAPPA  +  TWO) * RAD**2)
C     ***  OLDPHI IS USED TO DETECT LIMITS OF NUMERICAL ACCURACY.  IF
C     ***  WE RECOMPUTE STEP AND IT DOES NOT CHANGE, THEN WE ACCEPT IT.
      OLDPHI = ZERO
      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) = DNRM2(P, DIG,1)
      V(NREDUC) = ZERO
      V(DST0) = ZERO
      KAMIN = 3
      IF (V(DGNORM) .EQ. ZERO) KAMIN = 0
C
C     ***  STORE DIAG(DIHDI) IN W(DIAG0+1),...,W(DIAG0+P)  ***
C
      J = 0
      DO 10 I = 1, P
         J = J + I
         K1 = DIAG0 + I
         W(K1) = DIHDI(J)
 10      CONTINUE
C
C     ***  DETERMINE W(DGGDMX), THE LARGEST ELEMENT OF DIHDI  ***
C
      T1 = ZERO
      J = P * (P + 1) / 2
      DO 20 I = 1, J
         T = DABS(DIHDI(I))
         IF (T1 .LT. T) T1 = T
 20      CONTINUE
      W(DGGDMX) = T1
C
C  ***  TRY ALPHA = 0  ***
C
 30   CALL DLSQRT(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 DLITVM(IRC, W, L, W)
         T1 = DNRM2(IRC, W,1)
         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 = DLSVMN(P, L, W(Q), W(Q))
      IF (T .GE. ONE) GO TO 60
         IF (BIG .LE. ZERO) BIG = D1MACH(2)
         IF (V(DGNORM) .GE. T*T*BIG) GO TO 70
 60   CALL DLIVMU(P, W(Q), L, DIG)
      GTSTA = DDOT(P, W(Q),1,W(Q),1)
      V(NREDUC) = HALF * GTSTA
      CALL DLITVM(P, W(Q), L, W(Q))
      DST = DNRM2(P, W(Q),1)
      V(DST0) = DST
      PHI = DST - RAD
      IF (PHI .LE. PHIMAX) GO TO 260
      IF (RESTRT) GO TO 210
C
C  ***  PREPARE TO COMPUTE GERSCHGORIN ESTIMATES OF LARGEST (AND
C  ***  SMALLEST) EIGENVALUES.  ***
C
 70   K = 0
      DO 100 I = 1, P
         WI = ZERO
         IF (I .EQ. 1) GO TO 90
         IM1 = I - 1
         DO 80 J = 1, IM1
              K = K + 1
              T = DABS(DIHDI(K))
              WI = WI + T
              W(J) = W(J) + T
 80           CONTINUE
 90      W(I) = WI
         K = K + 1
 100     CONTINUE
C
C  ***  (UNDER-)ESTIMATE SMALLEST EIGENVALUE OF (D**-1)*H*(D**-1)  ***
C
      K = 1
      T1 = W(DIAG) - W(1)
      IF (P .LE. 1) GO TO 120
      DO 110 I = 2, P
         J = DIAG0 + I
         T = W(J) - W(I)
         IF (T .GE. T1) GO TO 110
              T1 = T
              K = I
 110     CONTINUE
C
 120  SK = W(K)
      J = DIAG0 + K
      AKK = W(J)
      K1 = K*(K-1)/2 + 1
      INC = 1
      T = ZERO
      DO 150 I = 1, P
         IF (I .EQ. K) GO TO 130
         AKI = DABS(DIHDI(K1))
         SI = W(I)
         J = DIAG0 + I
         T1 = HALF * (AKK - W(J) + SI - AKI)
         T1 = T1 + DSQRT(T1*T1 + SK*AKI)
         IF (T .LT. T1) T = T1
         IF (I .LT. K) GO TO 140
 130     INC = I
 140     K1 = K1 + INC
 150     CONTINUE
C
      W(EMIN) = AKK - T
      UK = V(DGNORM)/RAD - W(EMIN)
      IF (V(DGNORM) .EQ. ZERO) UK = UK + P001 + P001*UK
C
C  ***  COMPUTE GERSCHGORIN (OVER-)ESTIMATE OF LARGEST EIGENVALUE  ***
C
      K = 1
      T1 = W(DIAG) + W(1)
      IF (P .LE. 1) GO TO 170
      DO 160 I = 2, P
         J = DIAG0 + I
         T = W(J) + W(I)
         IF (T .LE. T1) GO TO 160
              T1 = T
              K = I
 160     CONTINUE
C
 170  SK = W(K)
      J = DIAG0 + K
      AKK = W(J)
      K1 = K*(K-1)/2 + 1
      INC = 1
      T = ZERO
      DO 200 I = 1, P
         IF (I .EQ. K) GO TO 180
         AKI = DABS(DIHDI(K1))
         SI = W(I)
         J = DIAG0 + I
         T1 = HALF * (W(J) + SI - AKI - AKK)
         T1 = T1 + DSQRT(T1*T1 + SK*AKI)
         IF (T .LT. T1) T = T1
         IF (I .LT. K) GO TO 190
 180     INC = I
 190     K1 = K1 + INC
 200     CONTINUE
C
      W(EMAX) = AKK + T
      LK = DMAX1(LK, V(DGNORM)/RAD - W(EMAX))
C
C     ***  ALPHAK = CURRENT VALUE OF ALPHA (SEE ALG. NOTES ABOVE).  WE
C     ***  USE MORE*S SCHEME FOR INITIALIZING IT.
      ALPHAK = DABS(V(STPPAR)) * V(RAD0)/RAD
C
      IF (IRC .NE. 0) GO TO 210
C
C  ***  COMPUTE L0 FOR POSITIVE DEFINITE H  ***
C
      CALL DLIVMU(P, W, L, W(Q))
      T = DNRM2(P, W,1)
      W(PHIPIN) = DST / T / T
      LK = DMAX1(LK, PHI*W(PHIPIN))
C
C  ***  SAFEGUARD ALPHAK AND ADD ALPHAK*I TO (D**-1)*H*(D**-1)  ***
C
 210  KA = KA + 1
      IF (-V(DST0) .GE. ALPHAK .OR. ALPHAK .LT. LK .OR. ALPHAK .GE. UK)
     1                      ALPHAK = UK * DMAX1(P001, DSQRT(LK/UK))
      IF(ALPHAK .LE. ZERO) ALPHAK = HALF * UK
      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 DLSQRT(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 DLITVM(IRC, W, L, W)
      T1 = DNRM2(IRC, W,1)
      LK = ALPHAK - T/T1/T1
      V(DST0) = -LK
      GO TO 210
C
C  ***  ALPHAK MAKES (D**-1)*H*(D**-1) POSITIVE DEFINITE.
C  ***  COMPUTE Q = -D*STEP, CHECK FOR CONVERGENCE.  ***
C
 240  CALL DLIVMU(P, W(Q), L, DIG)
      GTSTA = DDOT(P, W(Q),1,W(Q),1)
      CALL DLITVM(P, W(Q), L, W(Q))
      DST = DNRM2(P, W(Q),1)
      PHI = DST - RAD
      IF (PHI .LE. PHIMAX .AND. PHI .GE. PHIMIN) GO TO 270
      IF (PHI .EQ. OLDPHI) GO TO 270
      OLDPHI = PHI
      IF (PHI .LT. ZERO) GO TO 330
C
C  ***  UNACCEPTABLE ALPHAK -- UPDATE LK, UK, ALPHAK  ***
C
 250  IF (KA .GE. KALIM) GO TO 270
C     ***  THE FOLLOWING DMIN1 IS NECESSARY BECAUSE OF RESTARTS  ***
      IF (PHI .LT. ZERO) UK = DMIN1(UK, ALPHAK)
C     *** KAMIN = 0 ONLY IFF THE GRADIENT VANISHES  ***
      IF (KAMIN .EQ. 0) GO TO 210
      CALL DLIVMU(P, W, L, W(Q))
      T1 = DNRM2(P, W,1)
      ALPHAK = ALPHAK  +  (PHI/T1) * (DST/T1) * (DST/RAD)
      LK = DMAX1(LK, ALPHAK)
      GO TO 210
C
C  ***  ACCEPTABLE STEP ON FIRST TRY  ***
C
 260  ALPHAK = ZERO
C
C  ***  SUCCESSFUL STEP IN GENERAL.  COMPUTE STEP = -(D**-1)*Q  ***
C
 270  DO 280 I = 1, P
         J = Q0 + I
         STEP(I) = -W(J)/D(I)
 280     CONTINUE
      V(GTSTEP) = -GTSTA
      V(PREDUC) = HALF * (DABS(ALPHAK)*DST*DST + GTSTA)
      GO TO 410
C
C
C  ***  RESTART WITH NEW RADIUS  ***
C
 290  IF (V(DST0) .LE. ZERO .OR. V(DST0) - RAD .GT. PHIMAX) GO TO 310
C
C     ***  PREPARE TO RETURN NEWTON STEP  ***
C
         RESTRT = .TRUE.
         KA = KA + 1
         K = 0
         DO 300 I = 1, P
              K = K + I
              J = DIAG0 + I
              DIHDI(K) = W(J)
 300          CONTINUE
         UK = NEGONE
         GO TO 30
C
 310  KAMIN = KA + 3
      IF (V(DGNORM) .EQ. ZERO) KAMIN = 0
      IF (KA .EQ. 0) GO TO 50
C
      DST = W(DSTSAV)
      ALPHAK = DABS(V(STPPAR))
      PHI = DST - RAD
      T = V(DGNORM)/RAD
      UK = T - W(EMIN)
      IF (V(DGNORM) .EQ. ZERO) UK = UK + P001 + P001*UK
      IF (RAD .GT. V(RAD0)) GO TO 320
C
C        ***  SMALLER RADIUS  ***
         LK = ZERO
         IF (ALPHAK .GT. ZERO) LK = W(LK0)
         LK = DMAX1(LK, T - W(EMAX))
         IF (V(DST0) .GT. ZERO) LK = DMAX1(LK, (V(DST0)-RAD)*W(PHIPIN))
         GO TO 250
C
C     ***  BIGGER RADIUS  ***
 320  IF (ALPHAK .GT. ZERO) UK = DMIN1(UK, W(UK0))
      LK = DMAX1(ZERO, -V(DST0), T - W(EMAX))
      IF (V(DST0) .GT. ZERO) LK = DMAX1(LK, (V(DST0)-RAD)*W(PHIPIN))
      GO TO 250
C
C  ***  DECIDE WHETHER TO CHECK FOR SPECIAL CASE... IN PRACTICE (FROM
C  ***  THE STANDPOINT OF THE CALLING OPTIMIZATION CODE) IT SEEMS BEST
C  ***  NOT TO CHECK UNTIL A FEW ITERATIONS HAVE FAILED -- HENCE THE
C  ***  TEST ON KAMIN BELOW.
C
 330  DELTA = ALPHAK + DMIN1(ZERO, V(DST0))
      TWOPSI = ALPHAK*DST*DST + GTSTA
      IF (KA .GE. KAMIN) GO TO 340
C     *** IF THE TEST IN REF. 2 IS SATISFIED, FALL THROUGH TO HANDLE
C     *** THE SPECIAL CASE (AS SOON AS THE MORE-SORENSEN TEST DETECTS
C     *** IT).
      IF (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 DLSVMN TO OBTAIN APPROXIMATE EIGENVECTOR CORRESPONDING
C  ***  TO SMALLEST EIGENVALUE OF (D**-1)*H*(D**-1).  DLSVMN RETURNS
C  ***  X AND W WITH  L*W = X.
C
 340  T = DLSVMN(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 DLITVM(P, W, L, W)
      T2 = ONE/DNRM2(P, W,1)
      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 = DDOT(P, W(Q),1,W,1)
      T1 = (RAD + DST) * (RAD - DST)
      ROOT = DSQRT(SW*SW + T1)
      IF (SW .LT. ZERO) ROOT = -ROOT
      SI = T1 / (SW + ROOT)
C
C  ***  THE ACTUAL TEST FOR THE SPECIAL CASE...
C
      IF ((T2*SI)**2 .LE. EPS*(DST**2 + ALPHAK*RADSQ)) GO TO 380
C
C  ***  UPDATE UPPER BOUND ON SMALLEST EIGENVALUE (WHEN NOT POSITIVE)
C  ***  (AS RECOMMENDED BY MORE AND SORENSEN) AND CONTINUE...
C
      IF (V(DST0) .LE. ZERO) V(DST0) = DMIN1(V(DST0), T2**2 - ALPHAK)
      LK = DMAX1(LK, -V(DST0))
C
C  ***  CHECK WHETHER WE CAN HOPE TO DETECT THE SPECIAL CASE IN
C  ***  THE AVAILABLE ARITHMETIC.  ACCEPT STEP AS IT IS IF NOT.
C
C     ***  IF NOT YET AVAILABLE, OBTAIN MACHINE DEPENDENT VALUE DGXFAC.
 370  IF (DGXFAC .EQ. ZERO) DGXFAC = EPSFAC * D1MACH(4)
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*DDOT(P,W(X),1,W,1)))
      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) = DDOT(P, DIG,1,W(Q),1)
C
C  ***  SAVE VALUES FOR USE IN A POSSIBLE RESTART  ***
C
 410  V(DSTNRM) = DST
      V(STPPAR) = ALPHAK
      W(LK0) = LK
      W(UK0) = UK
      V(RAD0) = RAD
      W(DSTSAV) = DST
C
C     ***  RESTORE DIAGONAL OF DIHDI  ***
C
      J = 0
      DO 420 I = 1, P
         J = J + I
         K = DIAG0 + I
         DIHDI(J) = W(K)
 420     CONTINUE
C
 999  RETURN
C
C  ***  LAST CARD OF DGQTST FOLLOWS  ***
      END
      SUBROUTINE DASSST(D, IV, P, STEP, STLSTG, V, X, X0)
C
C  ***  ASSESS CANDIDATE STEP (***SOL VERSION 2.3)  ***
C
      INTEGER P, IV(32)
      DOUBLE PRECISION D(P), STEP(P), STLSTG(P), V(37), X(P), X0(P)
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      D (IN)  SCALE VECTOR USED IN COMPUTING V(RELDX) -- SEE BELOW.
C      P (IN)  NUMBER OF PARAMETERS BEING OPTIMIZED.
C   STEP (I/O) ON INPUT, STEP IS THE STEP TO BE ASSESSED.  IT IS UN-
C             CHANGED ON OUTPUT UNLESS A PREVIOUS STEP ACHIEVED A
C             BETTER OBJECTIVE FUNCTION REDUCTION, IN WHICH CASE STLSTG
C             WILL HAVE BEEN COPIED TO STEP.
C STLSTG (I/O) WHEN ASSESS RECOMMENDS RECOMPUTING STEP EVEN THOUGH THE
C             CURRENT (OR A PREVIOUS) STEP YIELDS AN OBJECTIVE FUNC-
C             TION DECREASE, IT SAVES IN STLSTG THE STEP THAT GAVE THE
C             BEST FUNCTION REDUCTION SEEN SO FAR (IN THE CURRENT ITERA-
C             TION).  IF THE RECOMPUTED STEP YIELDS A LARGER FUNCTION
C             VALUE, THEN STEP IS RESTORED FROM STLSTG AND
C             X = X0 + STEP IS RECOMPUTED.
C      V (I/O) REAL PARAMETER AND SCRATCH VECTOR -- SEE DESCRIPTION
C             BELOW OF V VALUES REFERENCED.
C      X (I/O) ON INPUT, X = X0 + STEP IS THE POINT AT WHICH THE OBJEC-
C             TIVE FUNCTION HAS JUST BEEN EVALUATED.  IF AN EARLIER
C             STEP YIELDED A BIGGER FUNCTION DECREASE, THEN X IS
C             RESTORED TO THE CORRESPONDING EARLIER VALUE.  OTHERWISE,
C             IF THE CURRENT STEP DOES NOT GIVE ANY FUNCTION DECREASE,
C             THEN X IS RESTORED TO X0.
C     X0 (IN)  INITIAL OBJECTIVE FUNCTION PARAMETER VECTOR (AT THE
C             START OF THE CURRENT ITERATION).
C
C  ***  IV VALUES REFERENCED  ***
C
C    IV(IRC) (I/O) ON INPUT FOR THE FIRST STEP TRIED IN A NEW ITERATION,
C             IV(IRC) SHOULD BE SET TO 3 OR 4 (THE VALUE TO WHICH IT IS
C             SET WHEN STEP IS DEFINITELY TO BE ACCEPTED).  ON INPUT
C             AFTER STEP HAS BEEN RECOMPUTED, IV(IRC) SHOULD BE
C             UNCHANGED SINCE THE PREVIOUS RETURN OF ASSESS.
C                ON OUTPUT, IV(IRC) IS A RETURN CODE HAVING ONE OF THE
C             FOLLOWING VALUES...
C                  1 = SWITCH MODELS OR TRY SMALLER STEP.
C                  2 = SWITCH MODELS OR ACCEPT STEP.
C                  3 = ACCEPT STEP AND DETERMINE V(RADFAC) BY GRADIENT
C                       TESTS.
C                  4 = ACCEPT STEP, V(RADFAC) HAS BEEN DETERMINED.
C                  5 = RECOMPUTE STEP (USING THE SAME MODEL).
C                  6 = RECOMPUTE STEP WITH RADIUS = V(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 0 UNLESS X AND V(F) HAVE BEEN RESTORED, IN
C             WHICH CASE ASSESS SETS IV(RESTOR) = 1.
C  IV(STAGE) (I/O) COUNT OF THE NUMBER OF MODELS TRIED SO FAR IN THE
C             CURRENT ITERATION.
C IV(STGLIM) (IN)  MAXIMUM NUMBER OF MODELS TO CONSIDER.
C IV(SWITCH) (OUT) SET TO 0 UNLESS A NEW MODEL IS BEING TRIED AND IT
C             GIVES A SMALLER FUNCTION VALUE THAN THE PREVIOUS MODEL,
C             IN WHICH CASE ASSESS SETS IV(SWITCH) = 1.
C IV(TOOBIG) (IN)  IS NONZERO IF STEP WAS TOO BIG (E.G. IF IT CAUSED
C             OVERFLOW).
C   IV(XIRC) (I/O) VALUE THAT IV(IRC) WOULD HAVE IN THE ABSENCE OF
C             CONVERGENCE, FALSE CONVERGENCE, AND OVERSIZED STEPS.
C
C  ***  V VALUES REFERENCED  ***
C
C V(AFCTOL) (IN)  ABSOLUTE FUNCTION CONVERGENCE TOLERANCE.  IF THE
C             ABSOLUTE VALUE OF THE CURRENT FUNCTION VALUE V(F) IS LESS
C             THAN V(AFCTOL), THEN ASSESS RETURNS WITH IV(IRC) = 10.
C V(DECFAC) (IN)  FACTOR BY WHICH TO DECREASE RADIUS WHEN IV(TOOBIG) IS
C             NONZERO.
C V(DSTNRM) (IN)  THE 2-NORM OF D*STEP.
C V(DSTSAV) (I/O) VALUE OF V(DSTNRM) ON SAVED STEP.
C   V(DST0) (IN)  THE 2-NORM OF D TIMES THE NEWTON STEP (WHEN DEFINED,
C             I.E., FOR V(NREDUC) .GE. 0).
C      V(F) (I/O) ON BOTH INPUT AND OUTPUT, V(F) IS THE OBJECTIVE FUNC-
C             TION VALUE AT X.  IF X IS RESTORED TO A PREVIOUS VALUE,
C             THEN V(F) IS RESTORED TO THE CORRESPONDING VALUE.
C   V(FDIF) (OUT) THE FUNCTION REDUCTION V(F0) - V(F) (FOR THE OUTPUT
C             VALUE OF V(F) IF AN EARLIER STEP GAVE A BIGGER FUNCTION
C             DECREASE, AND FOR THE INPUT VALUE OF V(F) OTHERWISE).
C V(FLSTGD) (I/O) SAVED VALUE OF V(F).
C     V(F0) (IN)  OBJECTIVE FUNCTION VALUE AT START OF ITERATION.
C V(GTSLST) (I/O) VALUE OF V(GTSTEP) ON SAVED STEP.
C V(GTSTEP) (IN)  INNER PRODUCT BETWEEN STEP AND GRADIENT.
C V(INCFAC) (IN)  MINIMUM FACTOR BY WHICH TO INCREASE RADIUS.
C  V(LMAXS) (IN)  MAXIMUM REASONABLE STEP SIZE (AND INITIAL STEP BOUND).
C             IF THE ACTUAL FUNCTION DECREASE IS NO MORE THAN TWICE
C             WHAT WAS PREDICTED, IF A RETURN WITH IV(IRC) = 7, 8, 9,
C             OR 10 DOES NOT OCCUR, IF V(DSTNRM) .GT. V(LMAXS), AND IF
C             V(PREDUC) .LE. V(SCTOL) * ABS(V(F0)), THEN ASSESS RE-
C             TURNS WITH IV(IRC) = 11.  IF SO DOING APPEARS WORTHWHILE,
C             THEN ASSESS REPEATS THIS TEST WITH V(PREDUC) COMPUTED FOR
C             A STEP OF LENGTH V(LMAXS) (BY A RETURN WITH IV(IRC) = 6).
C V(NREDUC) (I/O)  FUNCTION REDUCTION PREDICTED BY QUADRATIC MODEL FOR
C             NEWTON STEP.  IF ASSESS IS CALLED WITH IV(IRC) = 6, I.E.,
C             IF V(PREDUC) HAS BEEN COMPUTED WITH RADIUS = V(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) (OUT) SCALED RELATIVE CHANGE IN X CAUSED BY STEP, COMPUTED
C             BY FUNCTION  DRELST  AS
C                 MAX (D(I)*ABS(X(I)-X0(I)), 1 .LE. I .LE. P) /
C                    MAX (D(I)*(ABS(X(I))+ABS(X0(I))), 1 .LE. I .LE. P).
C             IF AN ACCEPTABLE STEP IS RETURNED, THEN V(RELDX) IS COM-
C             PUTED USING THE OUTPUT (POSSIBLY RESTORED) VALUES OF X
C             AND STEP.  OTHERWISE IT IS COMPUTED USING THE INPUT
C             VALUES.
C V(RFCTOL) (IN)  RELATIVE FUNCTION CONVERGENCE TOLERANCE.  IF THE
C             ACTUAL FUNCTION REDUCTION IS AT MOST TWICE WHAT WAS PRE-
C             DICTED AND  V(NREDUC) .LE. V(RFCTOL)*ABS(V(F0)),  THEN
C             ASSESS RETURNS WITH IV(IRC) = 8 OR 9.
C V(STPPAR) (IN)  MARQUARDT PARAMETER -- 0 MEANS FULL NEWTON STEP.
C V(TUNER1) (IN)  TUNING CONSTANT USED TO DECIDE IF THE FUNCTION
C             REDUCTION WAS MUCH LESS THAN EXPECTED.  SUGGESTED
C             VALUE = 0.1.
C V(TUNER2) (IN)  TUNING CONSTANT USED TO DECIDE IF THE FUNCTION
C             REDUCTION WAS LARGE ENOUGH TO ACCEPT STEP.  SUGGESTED
C             VALUE = 10**-4.
C V(TUNER3) (IN)  TUNING CONSTANT USED TO DECIDE IF THE RADIUS
C             SHOULD BE INCREASED.  SUGGESTED VALUE = 0.75.
C  V(XCTOL) (IN)  X-CONVERGENCE CRITERION.  IF STEP IS A NEWTON STEP
C             (V(STPPAR) = 0) HAVING V(RELDX) .LE. V(XCTOL) AND GIVING
C             AT MOST TWICE THE PREDICTED FUNCTION DECREASE, THEN
C             ASSESS RETURNS IV(IRC) = 7 OR 9.
C  V(XFTOL) (IN)  FALSE CONVERGENCE TOLERANCE.  IF STEP GAVE NO OR ONLY
C             A SMALL FUNCTION DECREASE AND V(RELDX) .LE. V(XFTOL),
C             THEN ASSESS RETURNS WITH IV(IRC) = 12.
C
C-------------------------------  NOTES  -------------------------------
C
C  ***  APPLICATION AND USAGE RESTRICTIONS  ***
C
C        THIS ROUTINE IS CALLED AS PART OF THE NL2SOL (NONLINEAR
C     LEAST-SQUARES) PACKAGE.  IT MAY BE USED IN ANY UNCONSTRAINED
C     MINIMIZATION SOLVER THAT USES DOGLEG, GOLDFELD-QUANDT-TROTTER,
C     OR LEVENBERG-MARQUARDT STEPS.
C
C  ***  ALGORITHM NOTES  ***
C
C        SEE (1) FOR FURTHER DISCUSSION OF THE ASSESSING AND MODEL
C     SWITCHING STRATEGIES.  WHILE NL2SOL CONSIDERS ONLY TWO MODELS,
C     ASSESS IS DESIGNED TO HANDLE ANY NUMBER OF MODELS.
C
C  ***  USAGE NOTES  ***
C
C        ON THE FIRST CALL OF AN ITERATION, ONLY THE I/O VARIABLES
C     STEP, X, IV(IRC), IV(MODEL), V(F), V(DSTNRM), V(GTSTEP), AND
C     V(PREDUC) NEED HAVE BEEN INITIALIZED.  BETWEEN CALLS, NO I/O
C     VALUES EXECPT STEP, X, IV(MODEL), V(F) AND THE STOPPING TOLER-
C     ANCES SHOULD BE CHANGED.
C        AFTER A RETURN FOR CONVERGENCE OR FALSE CONVERGENCE, ONE CAN
C     CHANGE THE STOPPING TOLERANCES AND CALL ASSESS AGAIN, IN WHICH
C     CASE THE STOPPING TESTS WILL BE REPEATED.
C
C  ***  REFERENCES  ***
C
C     (1) DENNIS, J.E., JR., GAY, D.M., AND WELSCH, R.E. (1981),
C        AN ADAPTIVE NONLINEAR LEAST-SQUARES ALGORITHM,
C        ACM TRANS. MATH. SOFTWARE, VOL. 7, NO. 3.
C
C     (2) POWELL, M.J.D. (1970)  A FORTRAN SUBROUTINE FOR SOLVING
C        SYSTEMS OF NONLINEAR ALGEBRAIC EQUATIONS, IN NUMERICAL
C        METHODS FOR NONLINEAR ALGEBRAIC EQUATIONS, EDITED BY
C        P. RABINOWITZ, GORDON AND BREACH, LONDON.
C
C  ***  HISTORY  ***
C
C        JOHN DENNIS DESIGNED MUCH OF THIS ROUTINE, STARTING WITH
C     IDEAS IN (2). ROY WELSCH SUGGESTED THE MODEL SWITCHING STRATEGY.
C        DAVID GAY AND STEPHEN PETERS CAST THIS SUBROUTINE INTO A MORE
C     PORTABLE FORM (WINTER 1977), AND DAVID GAY CAST IT INTO ITS
C     PRESENT FORM (FALL 1978).
C
C  ***  GENERAL  ***
C
C     THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH
C     SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS
C     MCS-7600324, DCR75-10143, 76-14311DSS, MCS76-11989, AND
C     MCS-7906671.
C
C------------------------  EXTERNAL QUANTITIES  ------------------------
C
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      EXTERNAL DRELST
      DOUBLE PRECISION DRELST
C
C DRELST... COMPUTES V(RELDX) = RELATIVE STEP SIZE.
C
C  ***  INTRINSIC FUNCTIONS  ***
C/+
      DOUBLE PRECISION DABS, DMAX1
C/
C  ***  NO COMMON BLOCKS  ***
C
C--------------------------  LOCAL VARIABLES  --------------------------
C
      LOGICAL GOODX
      INTEGER I, NFC
      DOUBLE PRECISION EMAX, EMAXS, GTS, HALF, ONE, RELDX1, RFAC1, TWO,
     1                 XMAX, 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
C/6
      DATA HALF/0.5D+0/, ONE/1.D+0/, TWO/2.D+0/, ZERO/0.D+0/
C/7
C     PARAMETER (HALF=0.5D+0, ONE=1.D+0, TWO=2.D+0, ZERO=0.D+0)
C/
C
C/6
      DATA 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/
C/7
C     PARAMETER (IRC=29, MLSTGD=32, MODEL=5, NFCALL=6, NFGCAL=7,
C    1           RADINC=8, RESTOR=9, STAGE=10, STGLIM=11, SWITCH=12,
C    2           TOOBIG=2, XIRC=13)
C/
C/6
      DATA 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/7
C     PARAMETER (AFCTOL=31, DECFAC=22, DSTNRM=2, DST0=3, DSTSAV=18,
C    1           F=10, FDIF=11, FLSTGD=12, F0=13, GTSLST=14, GTSTEP=4,
C    2           INCFAC=23, LMAXS=36, NREDUC=6, PLSTGD=15, PREDUC=7,
C    3           RADFAC=16, RDFCMN=24, RDFCMX=25, RELDX=17, RFCTOL=32,
C    4           SCTOL=37, STPPAR=5, TUNER1=26, TUNER2=27, TUNER3=28,
C    5           XCTOL=33, XFTOL=34)
C/
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 90
         IV(STAGE) = -1
         IV(XIRC) = I
         GO TO 60
C
C  ***  STEP WAS RECOMPUTED WITH NEW MODEL OR SMALLER RADIUS  ***
C  ***  FIRST DECIDE WHICH  ***
C
 20   IF (IV(MODEL) .NE. IV(MLSTGD)) GO TO 30
C        ***  OLD MODEL RETAINED, SMALLER RADIUS TRIED  ***
C        ***  DO NOT CONSIDER ANY MORE NEW MODELS THIS ITERATION  ***
         IV(STAGE) = IV(STGLIM)
         IV(RADINC) = -1
         GO TO 90
C
C  ***  A NEW MODEL IS BEING TRIED.  DECIDE WHETHER TO KEEP IT.  ***
C
 30   IV(STAGE) = IV(STAGE) + 1
C
C     ***  NOW WE ADD THE POSSIBILTIY THAT STEP WAS RECOMPUTED WITH  ***
C     ***  THE SAME MODEL, PERHAPS BECAUSE OF AN OVERSIZED STEP.     ***
C
 40   IF (IV(STAGE) .GT. 0) GO TO 50
C
C        ***  STEP WAS RECOMPUTED BECAUSE IT WAS TOO BIG.  ***
C
         IF (IV(TOOBIG) .NE. 0) GO TO 60
C
C        ***  RESTORE IV(STAGE) AND PICK UP WHERE WE LEFT OFF.  ***
C
         IV(STAGE) = -IV(STAGE)
         I = IV(XIRC)
         GO TO (20, 30, 90, 90, 70), I
C
 50   IF (IV(TOOBIG) .EQ. 0) GO TO 70
C
C  ***  HANDLE OVERSIZE STEP  ***
C
      IF (IV(RADINC) .GT. 0) GO TO 80
         IV(STAGE) = -IV(STAGE)
         IV(XIRC) = IV(IRC)
C
 60      V(RADFAC) = V(DECFAC)
         IV(RADINC) = IV(RADINC) - 1
         IV(IRC) = 5
         GO TO 999
C
 70   IF (V(F) .LT. V(FLSTGD)) GO TO 90
C
C     *** THE NEW STEP IS A LOSER.  RESTORE OLD MODEL.  ***
C
      IF (IV(MODEL) .EQ. IV(MLSTGD)) GO TO 80
         IV(MODEL) = IV(MLSTGD)
         IV(SWITCH) = 1
C
C     ***  RESTORE STEP, ETC. ONLY IF A PREVIOUS STEP DECREASED V(F).
C
 80   IF (V(FLSTGD) .GE. V(F0)) GO TO 90
         IV(RESTOR) = 1
         V(F) = V(FLSTGD)
         V(PREDUC) = V(PLSTGD)
         V(GTSTEP) = V(GTSLST)
         IF (IV(SWITCH) .EQ. 0) RFAC1 = V(DSTNRM) / V(DSTSAV)
         V(DSTNRM) = V(DSTSAV)
         NFC = IV(NFGCAL)
         GOODX = .FALSE.
C
C
C  ***  COMPUTE RELATIVE CHANGE IN X BY CURRENT STEP  ***
C
 90   RELDX1 = DRELST(P, D, X, X0)
C
C  ***  RESTORE X AND STEP IF NECESSARY  ***
C
      IF (GOODX) GO TO 110
      DO 100 I = 1, P
         STEP(I) = STLSTG(I)
         X(I) = X0(I) + STLSTG(I)
 100     CONTINUE
C
 110  V(FDIF) = V(F0) - V(F)
      IF (V(FDIF) .GT. V(TUNER2) * V(PREDUC)) GO TO 140
C
C        ***  NO (OR ONLY A TRIVIAL) FUNCTION DECREASE
C        ***  -- SO TRY NEW MODEL OR SMALLER RADIUS
C
         V(RELDX) = RELDX1
         IF (V(F) .LT. V(F0)) GO TO 120
              IV(MLSTGD) = IV(MODEL)
              V(FLSTGD) = V(F)
              V(F) = V(F0)
              CALL DCOPYGAMS(P,X0,1,X,1)
              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
      IF (GOODX) V(RELDX) = RELDX1
      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 * DMAX1(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) = DMAX1(V(INCFAC), HALF*GTS/(GTS + V(FDIF)))
         IV(IRC) = 4
         IF (V(STPPAR) .EQ. ZERO) 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)
      CALL DCOPYGAMS(P, STEP,1,STLSTG,1)
      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 (DABS(V(F)) .LT. V(AFCTOL)) IV(IRC) = 10
      IF (HALF * V(FDIF) .GT. V(PREDUC)) GO TO 999
      EMAX = V(RFCTOL) * DABS(V(F0))
      EMAXS = V(SCTOL) * DABS(V(F0))
      IF (V(DSTNRM) .GT. V(LMAXS) .AND. V(PREDUC) .LE. EMAXS)
     1                       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(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)
      IV(IRC) = 6
      CALL DCOPYGAMS(P, STEP,1,STLSTG,1)
      GO TO 999
C
C  ***  PERFORM SINGULAR CONVERGENCE TEST WITH RECOMPUTED V(PREDUC)  ***
C
 280  V(GTSTEP) = V(GTSLST)
      V(DSTNRM) = DABS(V(DSTSAV))
      CALL DCOPYGAMS(P, STLSTG,1,STEP,1)
      IV(IRC) = IV(XIRC)
      IF (V(DSTSAV) .LE. ZERO) IV(IRC) = 12
      V(NREDUC) = -V(PREDUC)
      V(PREDUC) = V(PLSTGD)
 290  IF (-V(NREDUC) .LE. V(RFCTOL) * DABS(V(F0))) IV(IRC) = 11
C
 999  RETURN
C
C  ***  LAST CARD OF ASSESS FOLLOWS  ***
      END
      SUBROUTINE DLITVM(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
      DOUBLE PRECISION X(N), L(1), Y(N)
      INTEGER I, II, IJ, IM1, I0, J, NP1
      DOUBLE PRECISION XI, ZERO
C/6
      DATA ZERO/0.D+0/
C/7
C     PARAMETER (ZERO=0.D+0)
C/
C
      DO 10 I = 1, N
 10      X(I) = Y(I)
      NP1 = N + 1
      I0 = N*(N+1)/2
      DO 30 II = 1, N
         I = NP1 - II
         XI = X(I)/L(I0)
         X(I) = XI
         IF (I .LE. 1) GO TO 999
         I0 = I0 - I
         IF (XI .EQ. ZERO) GO TO 30
         IM1 = I - 1
         DO 20 J = 1, IM1
              IJ = I0 + J
              X(J) = X(J) - XI*L(IJ)
 20           CONTINUE
 30      CONTINUE
 999  RETURN
C  ***  LAST CARD OF DLITVM FOLLOWS  ***
      END
      SUBROUTINE DLIVMU(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
      DOUBLE PRECISION X(N), L(1), Y(N)
      DOUBLE PRECISION DDOT
      INTEGER I, J, K
      DOUBLE PRECISION T, ZERO
C/6
      DATA ZERO/0.D+0/
C/7
C     PARAMETER (ZERO=0.D+0)
C/
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 = DDOT(I-1, L(J+1),1,X,1)
         J = J + I
         X(I) = (Y(I) - T)/L(J)
 30      CONTINUE
 999  RETURN
C  ***  LAST CARD OF DLIVMU FOLLOWS  ***
      END
      SUBROUTINE DLSQRT(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
      DOUBLE PRECISION 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
      DOUBLE PRECISION T, TD, ZERO
C
C  ***  INTRINSIC FUNCTIONS  ***
C/+
      DOUBLE PRECISION DSQRT
C/
C/6
      DATA ZERO/0.D+0/
C/7
C     PARAMETER (ZERO=0.D+0)
C/
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) = DSQRT(T)
 50      CONTINUE
C
      IRC = 0
      GO TO 999
C
 60   L(I0) = T
      IRC = I
C
 999  RETURN
C
C  ***  LAST CARD OF DLSQRT  ***
      END
      DOUBLE PRECISION FUNCTION DLSVMN(P, L, X, Y)
C
C  ***  ESTIMATE SMALLEST SING. VALUE OF PACKED LOWER TRIANG. MATRIX L
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER P
      DOUBLE PRECISION 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 DLSVMN 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 DLSVMN RETURNS ZERO, THEN SOME COMPONENTS OF X
C             ARE ZERO AND THE REST RETAIN THEIR INPUT VALUES.
C  Y (OUT) IF DLSVMN 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 DLSVMN 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 DLSVMN RETURNS).
C
C  ***  ALGORITHM NOTES  ***
C
C     THE ALGORITHM IS BASED ON (1), WITH THE ADDITIONAL PROVISION THAT
C     DLSVMN = 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
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
      DOUBLE PRECISION B, SMINUS, SPLUS, T, XMINUS, XPLUS
C
C  ***  CONSTANTS  ***
C
      DOUBLE PRECISION HALF, ONE, R9973, ZERO
C
C  ***  INTRINSIC FUNCTIONS  ***
C/+
      INTEGER MOD
      REAL FLOAT
      DOUBLE PRECISION DABS
C/
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      EXTERNAL   DVAXPY
      DOUBLE PRECISION DDOT, DNRM2
C
C/6
      DATA HALF/0.5D+0/, ONE/1.D+0/, R9973/9973.D+0/, ZERO/0.D+0/
C/7
C     PARAMETER (HALF=0.5D+0, ONE=1.D+0, R9973=9973.D+0, ZERO=0.D+0)
C/
C
C  ***  BODY  ***
C
      IX = 2
      PM1 = P - 1
C
C  ***  FIRST CHECK WHETHER TO RETURN DLSVMN = 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 = DABS(XPLUS)
         SMINUS = DABS(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 + DABS(X(I) + L(JI)*XPLUS)
              SMINUS = SMINUS + DABS(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 DVAXPY(JM1, X, XPLUS, L(J0+1), X)
 50      CONTINUE
C
C  ***  NORMALIZE X  ***
C
 60   T = ONE/DNRM2(P, X,1)
      DO 70 I = 1, P
 70      X(I) = T*X(I)
C
C  ***  SOLVE L*Y = X AND RETURN DLSVMN = 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 = DDOT(JM1,L(J0+1),1,Y,1)
      Y(J)=(X(J)-T)/L(JJ)
 100     CONTINUE
C
      DLSVMN = ONE/DNRM2(P, Y,1)
      GO TO 999
C
 110  DLSVMN = ZERO
 999  RETURN
C  ***  LAST CARD OF DLSVMN FOLLOWS  ***
      END
      DOUBLE PRECISION FUNCTION DRELST(P, D, X, X0)
C
C  ***  COMPUTE AND RETURN RELATIVE DIFFERENCE BETWEEN X AND X0  ***
C  ***  NL2SOL VERSION 2.2  ***
C
      INTEGER P
      DOUBLE PRECISION D(P), X(P), X0(P)
C/+
      DOUBLE PRECISION DABS
C/
      INTEGER I
      DOUBLE PRECISION EMAX, T, XMAX, ZERO
C/6
      DATA ZERO/0.D+0/
C/7
C     PARAMETER (ZERO=0.D+0)
C/
C
      EMAX = ZERO
      XMAX = ZERO
      DO 10 I = 1, P
         T = DABS(D(I) * (X(I) - X0(I)))
         IF (EMAX .LT. T) EMAX = T
         T = D(I) * (DABS(X(I)) + DABS(X0(I)))
         IF (XMAX .LT. T) XMAX = T
 10      CONTINUE
      DRELST = ZERO
      IF (XMAX .GT. ZERO) DRELST = EMAX / XMAX
 999  RETURN
C  ***  LAST CARD OF DRELST FOLLOWS  ***
      END
