      SUBROUTINE NL2SN(N, P, X, CALCR, IV, LIV, LV, V,
     1                  UIPARM, URPARM, UFPARM)
C
C  ***  MINIMIZE A NONLINEAR SUM OF SQUARES USING RESIDUAL VALUES ONLY..
C  ***  THIS AMOUNTS TO NL2S1 WITHOUT THE SUBROUTINE PARAMETER CALCJ.
C
C  ***  PARAMETERS  ***
C
      INTEGER N, P, LIV, LV
C/6
C     INTEGER IV(LIV), UIPARM(1)
C     REAL X(P), V(LV), URPARM(1)
C/7
      INTEGER IV(LIV), UIPARM(*)
      REAL X(P), V(LV), URPARM(*)
C/
      EXTERNAL CALCR, UFPARM
C
C-----------------------------  DISCUSSION  ----------------------------
C
C        THIS AMOUNTS TO SUBROUTINE NL2SNO (REF. 1) MODIFIED TO CALL
C     NL2IT.
C        THE PARAMETERS FOR NL2SN ARE THE SAME AS THOSE FOR NL2S1
C     (WHICH SEE), EXCEPT THAT CALCJ IS OMITTED.  INSTEAD OF CALLING
C     CALCJ TO OBTAIN THE JACOBIAN MATRIX OF R AT X, NL2SN COMPUTES
C     AN APPROXIMATION TO IT BY FINITE (FORWARD) DIFFERENCES -- SEE
C     V(DLTFDJ) BELOW.  NL2SN USES FUNCTION VALUES ONLY WHEN COMPUT-
C     THE COVARIANCE MATRIX (RATHER THAN THE FUNCTIONS AND GRADIENTS
C     THAT NL2S1 MAY USE).  TO DO SO, NL2SN SETS IV(COVREQ) TO -1 IF
C     IV(COVPRT) = 1 WITH IV(COVREQ) = 0 AND TO MINUS ITS ABSOLUTE
C     VALUE OTHERWISE.  THUS V(DELTA0) IS NEVER REFERENCED AND ONLY
C     V(DLTFDC) MATTERS -- SEE NL2SOL FOR A DESCRIPTION OF V(DLTFDC).
C        THE NUMBER OF EXTRA CALLS ON CALCR USED IN COMPUTING THE JACO-
C     BIAN APPROXIMATION ARE NOT INCLUDED IN THE FUNCTION EVALUATION
C     COUNT IV(NFCALL) AND ARE NOT OTHERWISE REPORTED.
C
C V(DLTFDJ)... V(43) HELPS CHOOSE THE STEP SIZE USED WHEN COMPUTING THE
C             FINITE-DIFFERENCE JACOBIAN MATRIX.  FOR DIFFERENCES IN-
C             VOLVING X(I), THE STEP SIZE FIRST TRIED IS
C                       V(DLTFDJ) * MAX(ABS(X(I)), 1/D(I)),
C             WHERE D IS THE CURRENT SCALE VECTOR (SEE REF. 1).  (IF
C             THIS STEP IS TOO BIG, I.E., IF CALCR SETS NF TO 0, THEN
C             SMALLER STEPS ARE TRIED UNTIL THE STEP SIZE IS SHRUNK BE-
C             LOW 1000 * MACHEP, WHERE MACHEP IS THE UNIT ROUNDOFF.
C             DEFAULT = MACHEP**0.5.
C
C  ***  REFERENCES  ***
C
C 1.  DENNIS, J.E., GAY, D.M., AND WELSCH, R.E. (1981), AN ADAPTIVE
C             NONLINEAR LEAST-SQUARES ALGORITHM, ACM TRANS. MATH.
C             SOFTWARE, VOL. 7, NO. 3.
C
C  ***  GENERAL  ***
C
C     CODED BY DAVID M. GAY.
C
C+++++++++++++++++++++++++++  DECLARATIONS  +++++++++++++++++++++++++++
C
C  ***  EXTERNAL SUBROUTINES  ***
C
      EXTERNAL DEFLT, NL2IT, MLPRDP
C
C DEFLT.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS.
C NL2IT.... CARRIES OUT OPTIMIZATION ITERATIONS.
C MLPRDP... PRINTS REGRESSION DIAGNOSTICS.
C
C  ***  INTRINSIC FUNCTIONS  ***
C/+
      INTEGER IABS
      REAL  ABS, AMAX1
C/
C
C  ***  LOCAL VARIABLES  ***
C
      LOGICAL ONERD, RSAVE, TWORD
      INTEGER D1, DK, DR1, I, IN, IV1, J1K, K, NF, RD1, RD2,
     1        RSAVE1, R1, RN
      REAL H, HFAC, HLIM, NEGPT5, ONE, XK, ZERO
C
C  ***  IV AND V COMPONENTS  ***
C
      INTEGER COVPRT, COVREQ, D, DINIT, DLTFDJ, J, LASTV, MODE, NEXTV,
     1        NF0, NF1, NFCALL, NFGCAL, R, RDREQ, REGD, TOOBIG, VNEED
C/6
C     DATA COVPRT/14/, COVREQ/15/, D/27/, DINIT/38/, DLTFDJ/43/, J/70/,
C    1     LASTV/45/, MODE/35/, NEXTV/47/, NF0/68/, NF1/69/, NFCALL/6/,
C    2     NFGCAL/7/, R/61/, RDREQ/57/, REGD/67/, TOOBIG/2/, VNEED/4/
C/7
      PARAMETER (COVPRT=14, COVREQ=15, D=27, DINIT=38, DLTFDJ=43, J=70,
     1           LASTV=45, MODE=35, NEXTV=47, NF0=68, NF1=69, NFCALL=6,
     2           NFGCAL=7, R=61, RDREQ=57, REGD=67, TOOBIG=2, VNEED=4)
C/
      DATA HLIM/0.E+0/
      DATA HFAC/1.E+3/, NEGPT5/-0.5E+0/, ONE/1.E+0/, ZERO/0.E+0/
C
C---------------------------------  BODY  ------------------------------
C
      IF (IV(1) .EQ. 0) CALL DEFLT(1, IV, LIV, LV, V)
      IV(COVREQ) = -IABS(IV(COVREQ))
      IF (IV(COVPRT) .NE. 0 .AND. IV(COVREQ) .EQ. 0) IV(COVREQ) = -1
      IV1 = IV(1)
      ONERD = IV(RDREQ) .GE. 1
      TWORD = IV(RDREQ) .GE. 3
      IF (IV1 .EQ. 14) GO TO 10
      IF (IV1 .GT. 2 .AND. IV1 .LT. 12) GO TO 10
      IF (IV1 .EQ. 12) IV(1) = 13
      R1 = 1
      D1 = 1
      DR1 = 1
      RD1 = 1
      RD2 = 1
      RSAVE = .FALSE.
      GO TO 20
C
 10   D1 = IV(D)
      DR1 = IV(J)
      R1 = IV(R)
      RN = R1 + N - 1
      RD1 = R1
      IF (ONERD) RD1 = RD1 + N
      RD2 = RD1
      IF (TWORD) RD2 = RD2 + N
      RSAVE1 = RD2 + N
      IF (ONERD) RSAVE1 = RSAVE1 + N
      RSAVE = RSAVE1 + N .LE. DR1
C
 20   IN = IV(VNEED) + P + N*(P+1)
      IF (ONERD) IN = IN + N
      IF (TWORD) IN = IN + N
      IF (RSAVE) IN = IN + N
      IV(VNEED) = IN
C
 30   CALL NL2IT(V(D1), V(DR1), IV, LIV, LV, N, 1, 1, N, P, V(R1),
     1           V(RD1), V(RD2), V, X)
      IF (IABS(IV(1))-2) 50, 40, 120
C
C   ***  NEW GRADIENT REQUESTED -- MAKE SURE R IS CURRENT  ***
C
 40   NF = IV(NFGCAL)
      IF (NF .EQ. IV(NF1)) GO TO 70
         IF (NF .NE. IV(NF0)) GO TO 60
              CALL SCOPY(N, V(R1),1, V(RSAVE1),1)
              IV(NF1) = IV(NF0)
              GO TO 70
C
C  ***  NEW FUNCTION VALUE (R VALUE) NEEDED -- SAVE OLD ONE FIRST  ***
C
 50   NF = IV(NFCALL)
      IF (IV(MODE) .NE. 0 .OR. .NOT. RSAVE) GO TO 60
         CALL SCOPY(N, V(RSAVE1),1, V(R1),1)
         IV(NF0) = IV(NF1)
C
 60   CALL CALCR(N, P, X, NF, V(R1), UIPARM, URPARM, UFPARM)
      IV(NF1) = NF
      IF (NF .LE. 0) IV(TOOBIG) = 1
      IF (IV(1) .EQ. 1) GO TO 30
C
C  ***  COMPUTE FINITE-DIFFERENCE APPROXIMATION TO DR = GRAD. OF R  ***
C
 70   IF (IV(1) .LT. 0) GO TO 30
C
C     *** INITIALIZE D IF NECESSARY ***
C
      IF (IV(MODE) .LT. 0 .AND. V(DINIT) .EQ. ZERO)
     1        CALL VSCOPY(P, V(D1), ONE)
C
      J1K = DR1
      DK = D1
      DO 110 K = 1, P
         XK = X(K)
         H = V(DLTFDJ) * AMAX1( ABS(XK), ONE/V(DK))
         DK = DK + 1
 80      X(K) = XK + H
         NF = IV(NFGCAL)
         CALL CALCR (N, P, X, NF, V(J1K), UIPARM, URPARM, UFPARM)
         IF (NF .GT. 0) GO TO 90
              IF (HLIM .EQ. ZERO) HLIM = HFAC * R1MACH(4)
C             ***  HLIM = HFAC TIMES THE UNIT ROUNDOFF  ***
              H = NEGPT5 * H
              IF ( ABS(H) .GE. HLIM) GO TO 80
                   IV(TOOBIG) = 1
                   GO TO 30
 90      X(K) = XK
         DO 100 I = R1, RN
              V(J1K) = (V(J1K) - V(I)) / H
              J1K = J1K + 1
 100          CONTINUE
 110     CONTINUE
      GO TO 30
C
C  ***  STORAGE ALLOCATION  ***
C
 120  IF (IV(1) .NE. 14) GO TO 140
      IV(D) = IV(NEXTV)
      IV(R) = IV(D) + P
      IN = IV(R) + N
      IF (ONERD) IN = IN + N
      IF (TWORD) IN = IN + N
      IF (IN + N*(P+1) .GT. LV + 1) GO TO 130
         IN = IN + N
         IV(LASTV) = IV(LASTV) + N
 130  IV(J) = IN
      IV(NEXTV) = IN + N*P
      IV(NF0) = 0
      IV(NF1) = 0
      IF (IV1 .EQ. 13) GO TO 999
      GO TO 10
C
 140  IF (IV(REGD) .GT. 0) IV(REGD) = RD1
      CALL MLPRDP(IV, LIV, LV, N, V(RD1), V(RD2), V)
C
 999  RETURN
C
C  ***  LAST CARD OF NL2SN FOLLOWS  ***
      END
      REAL FUNCTION R1MACH(I)
C***BEGIN PROLOGUE  R1MACH
C***DATE WRITTEN   790101   (YYMMDD)
C***REVISION DATE  860825   (YYMMDD)
C***CATEGORY NO.  R1
C***KEYWORDS  MACHINE CONSTANTS
C***AUTHOR  FOX, P. A., (BELL LABS)
C           HALL, A. D., (BELL LABS)
C           SCHRYER, N. L., (BELL LABS)
C***PURPOSE  Returns single precision machine dependent constants
C***DESCRIPTION
C
C     This is the CMLIB version of R1MACH, the real machine
C     constants subroutine originally developed for the PORT library.
C
C     R1MACH can be used to obtain machine-dependent parameters
C     for the local machine environment.  It is a function
C     subroutine with one (input) argument, and can be called
C     as follows, for example
C
C          A = R1MACH(I)
C
C     where I=1,...,5.  The (output) value of A above is
C     determined by the (input) value of I.  The results for
C     various values of I are discussed below.
C
C  Single-Precision Machine Constants
C  R1MACH(1) = B**(EMIN-1), the smallest positive magnitude.
C  R1MACH(2) = B**EMAX*(1 - B**(-T)), the largest magnitude.
C  R1MACH(3) = B**(-T), the smallest relative spacing.
C  R1MACH(4) = B**(1-T), the largest relative spacing.
C  R1MACH(5) = LOG10(B)
C***REFERENCES  FOX, P.A., HALL, A.D., SCHRYER, N.L, *FRAMEWORK FOR
C                 A PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHE-
C                 MATICAL SOFTWARE, VOL. 4, NO. 2, JUNE 1978,
C                 PP. 177-188.
C***ROUTINES CALLED  XERROR
C***END PROLOGUE  R1MACH
C
      INTEGER SMALL(2)
      INTEGER LARGE(2)
      INTEGER RIGHT(2)
      INTEGER DIVER(2)
      INTEGER LOG10(2)
C
      REAL RMACH(5)
C
      EQUIVALENCE (RMACH(1),SMALL(1))
      EQUIVALENCE (RMACH(2),LARGE(1))
      EQUIVALENCE (RMACH(3),RIGHT(1))
      EQUIVALENCE (RMACH(4),DIVER(1))
      EQUIVALENCE (RMACH(5),LOG10(1))
C
C
C     MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T
C     3B SERIES, MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T
C     PC 7300), SUN SPARCSTATIONS, SILICON GRAPHCS WORKSTATIONS, HP
C     9000 WORKSTATIONS, IBM RS/6000 WORKSTATIONS, AND 8087 BASED 
C     MICROS (E.G. IBM PC AND AT&T 6300).
C
C === MACHINE = ATT.3B
C === MACHINE = ATT.6300
C === MACHINE = ATT.7300
C === MACHINE = HP.9000
C === MACHINE = IBM.PC
C === MACHINE = IBM.RS6000
C === MACHINE = IEEE.MOST-SIG-BYTE-FIRST
C === MACHINE = IEEE.LEAST-SIG-BYTE-FIRST
C === MACHINE = SGI
C === MACHINE = SUN
C === MACHINE = 68000
C === MACHINE = 8087
C      DATA SMALL(1) /     8388608 /
C      DATA LARGE(1) /  2139095039 /
C      DATA RIGHT(1) /   864026624 /
C      DATA DIVER(1) /   872415232 /
C      DATA LOG10(1) /  1050288283 /
C
C     MACHINE CONSTANTS FOR SUN WORKSTATIONS.  f77 WITH -r8 OPTION.
C     MACHINE CONSTANTS FOR IBM RS/6000 WORKSTATIONS WITH -qautodbl=dblpad.
C
C === MACHINE = IBM.RS6000.XLF-WITH-AUTODBL-OPTION
C === MACHINE = SUN.F77-WITH-R8-OPTION
C === MACHINE = SGI.ORIGIN.F77-WITH-R8-D16-OPTIONS
C      DATA RMACH(1) / 2.2250738585072E-307 /
C      DATA RMACH(2) / 1.7976931348623E308  /
C      DATA RMACH(3) / 1.1102230246252E-16  /
C      DATA RMACH(4) / 2.2204460492503E-16  /
C      DATA RMACH(5) / 0.30102999566398     /
C
C     MACHINE CONSTANTS FOR AMDAHL MACHINES.
C
C === MACHINE = AMDAHL
C      DATA SMALL(1) /    1048576 /
C      DATA LARGE(1) / 2147483647 /
C      DATA RIGHT(1) /  990904320 /
C      DATA DIVER(1) / 1007681536 /
C      DATA LOG10(1) / 1091781651 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.
C
C === MACHINE = BURROUGHS.1700
C      DATA RMACH(1) / Z400800000 /
C      DATA RMACH(2) / Z5FFFFFFFF /
C      DATA RMACH(3) / Z4E9800000 /
C      DATA RMACH(4) / Z4EA800000 /
C      DATA RMACH(5) / Z500E730E8 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 5700/6700/7700 SYSTEMS.
C
C === MACHINE = BURROUGHS.5700
C === MACHINE = BURROUGHS.6700
C === MACHINE = BURROUGHS.7700
C      DATA RMACH(1) / O1771000000000000 /
C      DATA RMACH(2) / O0777777777777777 /
C      DATA RMACH(3) / O1311000000000000 /
C      DATA RMACH(4) / O1301000000000000 /
C      DATA RMACH(5) / O1157163034761675 /
C
C     MACHINE CONSTANTS FOR THE CONVEX C1, C2, C3 SERIES (NATIVE MODE)
C
C === MACHINE = CONVEX
C      DATA RMACH(1) / 2.9387360E-39 /
C      DATA RMACH(2) / 1.7014117E+38 /
C      DATA RMACH(3) / 5.9604645E-08 /
C      DATA RMACH(4) / 1.1920929E-07 /
C      DATA RMACH(5) / 3.0102999E-01 /
C
C     MACHINE CONSTANTS FOR THE CONVEX C1, C2, C3 (NATIVE MODE)
C     WITH -P8 OPTION
C
C === MACHINE = CONVEX.P8
C      DATA RMACH(1) / 5.562684646268007E-309 /
C      DATA RMACH(2) / 8.988465674311577E+307 /
C      DATA RMACH(3) / 1.110223024625157E-016 /
C      DATA RMACH(4) / 2.220446049250313E-016 /
C      DATA RMACH(5) / 3.010299956639812E-001 /
C
C     MACHINE CONSTANTS FOR THE CONVEX C1, C2, C3 (IEEE MODE)
C
C === MACHINE = CONVEX.IEEE
C      DATA RMACH(1) / 1.1754945E-38 /
C      DATA RMACH(2) / 3.4028234E+38 /
C      DATA RMACH(3) / 5.9604645E-08 /
C      DATA RMACH(4) / 1.1920929E-07 /
C      DATA RMACH(5) / 3.0102999E-01 /
C
C     MACHINE CONSTANTS FOR THE CONVEX C1, C2, C3 (IEEE MODE)
C     WITH -P8 OPTION
C
C === MACHINE = CONVEX.IEEE.P8
C      DATA RMACH(1) / 2.225073858507202E-308 /
C      DATA RMACH(2) / 1.797693134862315E+308 /
C      DATA RMACH(3) / 1.110223024625157E-016 /
C      DATA RMACH(4) / 2.220446049250313E-016 /
C      DATA RMACH(5) / 3.010299956639812E-001 /
C
C     MACHINE CONSTANTS FOR THE CYBER 170/180 SERIES USING NOS (FTN5).
C
C === MACHINE = CYBER.170.NOS
C === MACHINE = CYBER.180.NOS
C      DATA RMACH(1) / O"00014000000000000000" /
C      DATA RMACH(2) / O"37767777777777777777" /
C      DATA RMACH(3) / O"16404000000000000000" /
C      DATA RMACH(4) / O"16414000000000000000" /
C      DATA RMACH(5) / O"17164642023241175720" /
C
C     MACHINE CONSTANTS FOR THE CDC 180 SERIES USING NOS/VE
C
C === MACHINE = CYBER.180.NOS/VE
C      DATA RMACH(1) / Z"3001800000000000" /
C      DATA RMACH(2) / Z"4FFEFFFFFFFFFFFE" /
C      DATA RMACH(3) / Z"3FD2800000000000" /
C      DATA RMACH(4) / Z"3FD3800000000000" /
C      DATA RMACH(5) / Z"3FFF9A209A84FBCF" /
C
C     MACHINE CONSTANTS FOR THE CYBER 205
C
C === MACHINE = CYBER.205
C      DATA RMACH(1) / X'9000400000000000' /
C      DATA RMACH(2) / X'6FFF7FFFFFFFFFFF' /
C      DATA RMACH(3) / X'FFA3400000000000' /
C      DATA RMACH(4) / X'FFA4400000000000' /
C      DATA RMACH(5) / X'FFD04D104D427DE8' /
C
C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES.
C
C === MACHINE = CDC.6000
C === MACHINE = CDC.7000
C      DATA RMACH(1) / 00014000000000000000B /
C      DATA RMACH(2) / 37767777777777777777B /
C      DATA RMACH(3) / 16404000000000000000B /
C      DATA RMACH(4) / 16414000000000000000B /
C      DATA RMACH(5) / 17164642023241175720B /
C
C     MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3.
C
C === MACHINE = CRAY.46-BIT-INTEGER
C === MACHINE = CRAY.64-BIT-INTEGER
C      DATA RMACH(1) / 200034000000000000000B /
C      DATA RMACH(2) / 577767777777777777776B /
C      DATA RMACH(3) / 377224000000000000000B /
C      DATA RMACH(4) / 377234000000000000000B /
C      DATA RMACH(5) / 377774642023241175720B /
C
C     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200
C
C     NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING LINE -
C     STATIC RMACH(5)
C
C === MACHINE = DATA_GENERAL.ECLIPSE.S/200
C      DATA SMALL/20K,0/,LARGE/77777K,177777K/
C      DATA RIGHT/35420K,0/,DIVER/36020K,0/
C      DATA LOG10/40423K,42023K/
C
C     ELXSI 6400
C
C === MACHINE = ELSXI.6400
C      DATA SMALL(1) / '00800000'X /
C      DATA LARGE(1) / '7F7FFFFF'X /
C      DATA RIGHT(1) / '33800000'X /
C      DATA DIVER(1) / '34000000'X /
C      DATA LOG10(1) / '3E9A209B'X /
C
C     MACHINE CONSTANTS FOR THE HARRIS 220
C     MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7
C
C === MACHINE = HARRIS.220
C === MACHINE = HARRIS.SLASH6
C === MACHINE = HARRIS.SLASH7
C      DATA SMALL(1),SMALL(2) / '20000000, '00000201 /
C      DATA LARGE(1),LARGE(2) / '37777777, '00000177 /
C      DATA RIGHT(1),RIGHT(2) / '20000000, '00000352 /
C      DATA DIVER(1),DIVER(2) / '20000000, '00000353 /
C      DATA LOG10(1),LOG10(2) / '23210115, '00000377 /
C
C     MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES.
C     MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES.
C
C === MACHINE = HONEYWELL.600/6000
C === MACHINE = HONEYWELL.DPS.8/70
C      DATA RMACH(1) / O402400000000 /
C      DATA RMACH(2) / O376777777777 /
C      DATA RMACH(3) / O714400000000 /
C      DATA RMACH(4) / O716400000000 /
C      DATA RMACH(5) / O776464202324 /
C
C     MACHINE CONSTANTS FOR THE HP 2100
C     3 WORD DOUBLE PRECISION WITH FTN4
C
C === MACHINE = HP.2100.3_WORD_DP
C      DATA SMALL(1), SMALL(2) / 40000B,       1 /
C      DATA LARGE(1), LARGE(2) / 77777B, 177776B /
C      DATA RIGHT(1), RIGHT(2) / 40000B,    325B /
C      DATA DIVER(1), DIVER(2) / 40000B,    327B /
C      DATA LOG10(1), LOG10(2) / 46420B,  46777B /
C
C     MACHINE CONSTANTS FOR THE HP 2100
C     4 WORD DOUBLE PRECISION WITH FTN4
C
C === MACHINE = HP.2100.4_WORD_DP
C      DATA SMALL(1), SMALL(2) / 40000B,       1 /
C      DATA LARGE91), LARGE(2) / 77777B, 177776B /
C      DATA RIGHT(1), RIGHT(2) / 40000B,    325B /
C      DATA DIVER(1), DIVER(2) / 40000B,    327B /
C      DATA LOG10(1), LOG10(2) / 46420B,  46777B /
C
C     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
C     THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86 AND
C     THE INTERDATA 3230 AND INTERDATA 7/32.
C
C === MACHINE = IBM.360
C === MACHINE = IBM.370
C === MACHINE = XEROX.SIGMA.5
C === MACHINE = XEROX.SIGMA.7
C === MACHINE = XEROX.SIGMA.9
C === MACHINE = SEL.85
C === MACHINE = SEL.86
C === MACHINE = INTERDATA.3230
C === MACHINE = INTERDATA.7/32
C      DATA RMACH(1) / Z00100000 /
C      DATA RMACH(2) / Z7FFFFFFF /
C      DATA RMACH(3) / Z3B100000 /
C      DATA RMACH(4) / Z3C100000 /
C      DATA RMACH(5) / Z41134413 /
C
C     MACHINE CONSTANTS FOR THE INTERDATA 8/32
C     WITH THE UNIX SYSTEM FORTRAN 77 COMPILER.
C
C     FOR THE INTERDATA FORTRAN VII COMPILER REPLACE
C     THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S.
C
C === MACHINE = INTERDATA.8/32.UNIX
C      DATA RMACH(1) / Z'00100000' /
C      DATA RMACH(2) / Z'7EFFFFFF' /
C      DATA RMACH(3) / Z'3B100000' /
C      DATA RMACH(4) / Z'3C100000' /
C      DATA RMACH(5) / Z'41134413' /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KA OR KI PROCESSOR).
C
C === MACHINE = PDP-10.KA
C === MACHINE = PDP-10.KI
C      DATA RMACH(1) / "000400000000 /
C      DATA RMACH(2) / "377777777777 /
C      DATA RMACH(3) / "146400000000 /
C      DATA RMACH(4) / "147400000000 /
C      DATA RMACH(5) / "177464202324 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C     32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
C
C === MACHINE = PDP-11.32-BIT
C      DATA SMALL(1) /    8388608 /
C      DATA LARGE(1) / 2147483647 /
C      DATA RIGHT(1) /  880803840 /
C      DATA DIVER(1) /  889192448 /
C      DATA LOG10(1) / 1067065499 /
C
C      DATA RMACH(1) / O00040000000 /
C      DATA RMACH(2) / O17777777777 /
C      DATA RMACH(3) / O06440000000 /
C      DATA RMACH(4) / O06500000000 /
C      DATA RMACH(5) / O07746420233 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C     16-BIT INTEGERS  (EXPRESSED IN INTEGER AND OCTAL).
C
C === MACHINE = PDP-11.16-BIT
C      DATA SMALL(1),SMALL(2) /   128,     0 /
C      DATA LARGE(1),LARGE(2) / 32767,    -1 /
C      DATA RIGHT(1),RIGHT(2) / 13440,     0 /
C      DATA DIVER(1),DIVER(2) / 13568,     0 /
C      DATA LOG10(1),LOG10(2) / 16282,  8347 /
C
C      DATA SMALL(1),SMALL(2) / O000200, O000000 /
C      DATA LARGE(1),LARGE(2) / O077777, O177777 /
C      DATA RIGHT(1),RIGHT(2) / O032200, O000000 /
C      DATA DIVER(1),DIVER(2) / O032400, O000000 /
C      DATA LOG10(1),LOG10(2) / O037632, O020233 /
C
C     MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000.
C
C === MACHINE = SEQUENT.BALANCE.8000
C      DATA SMALL(1) / $00800000 /
C      DATA LARGE(1) / $7F7FFFFF /
C      DATA RIGHT(1) / $33800000 /
C      DATA DIVER(1) / $34000000 /
C      DATA LOG10(1) / $3E9A209B /
C
C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES.
C
C === MACHINE = UNIVAC.1100
C      DATA RMACH(1) / O000400000000 /
C      DATA RMACH(2) / O377777777777 /
C      DATA RMACH(3) / O146400000000 /
C      DATA RMACH(4) / O147400000000 /
C      DATA RMACH(5) / O177464202324 /
C
C     MACHINE CONSTANTS FOR THE VAX 11/780
C    (EXPRESSED IN INTEGER AND HEXADECIMAL)
C  *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***
C
C === MACHINE = VAX.11/780
C      DATA SMALL(1) /       128 /
C      DATA LARGE(1) /    -32769 /
C      DATA RIGHT(1) /     13440 /
C      DATA DIVER(1) /     13568 /
C      DATA LOG10(1) / 547045274 /
C
C  ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSTEMS***
C
C      DATA SMALL(1) / Z00000080 /
C      DATA LARGE(1) / ZFFFF7FFF /
C      DATA RIGHT(1) / Z00003480 /
C      DATA DIVER(1) / Z00003500 /
C      DATA LOG10(1) / Z209B3F9A /
C
C
C***FIRST EXECUTABLE STATEMENT  R1MACH
C
      R1MACH = RMACH(I)
      RETURN
C
      END
      SUBROUTINE SCOPY(N,SX,INCX,SY,INCY)
C***BEGIN PROLOGUE  SCOPY
C***DATE WRITTEN   791001   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***REVISION HISTORY  (YYMMDD)
C   000330  Modified array declarations.  (JEC)
C
C***CATEGORY NO.  D1A5
C***KEYWORDS  BLAS,COPY,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  Copy s.p. vector 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       SX  single precision vector with N elements
C     INCX  storage spacing between elements of SX
C       SY  single precision vector with N elements
C     INCY  storage spacing between elements of SY
C
C     --Output--
C       SY  copy of vector SX (unchanged if N .LE. 0)
C
C     Copy single precision SX to single precision SY.
C     For I = 0 to N-1, copy  SX(LX+I*INCX) to SY(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  SCOPY
C
      REAL SX(*),SY(*)
C***FIRST EXECUTABLE STATEMENT  SCOPY
      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
        SY(IY) = SX(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
        SY(I) = SX(I)
   30 CONTINUE
      IF( N .LT. 7 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,7
        SY(I) = SX(I)
        SY(I + 1) = SX(I + 1)
        SY(I + 2) = SX(I + 2)
        SY(I + 3) = SX(I + 3)
        SY(I + 4) = SX(I + 4)
        SY(I + 5) = SX(I + 5)
        SY(I + 6) = SX(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
          SY(I) = SX(I)
   70     CONTINUE
      RETURN
      END
      SUBROUTINE DEFLT(ALG, IV, LIV, LV, V)
C
C  ***  SUPPLY ***SOL (VERSION 2.3) DEFAULT VALUES TO IV AND V  ***
C
C  ***  ALG = 1 MEANS REGRESSION CONSTANTS.
C  ***  ALG = 2 MEANS GENERAL UNCONSTRAINED OPTIMIZATION CONSTANTS.
C
      INTEGER LIV, LV
      INTEGER ALG, IV(LIV)
      REAL V(LV)
C
      EXTERNAL  VDFLT
C VDFLT.... 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 VDFLT(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 DEFLT FOLLOWS  ***
      END
      SUBROUTINE VSCOPY(P, Y, S)
C
C  ***  SET P-VECTOR Y TO SCALAR S  ***
C
      INTEGER P
      REAL S, Y(P)
C
      INTEGER I
C
      DO 10 I = 1, P
 10      Y(I) = S
      RETURN
      END
      SUBROUTINE MLPRDP(IV, LIV, LV, N, RD1, RD2, V)
C
C  ***  PRINT REGRESSION DIAGNOSTICS FOR MLPSL AND NL2S1 ***
C
      INTEGER LIV, LV, N
      INTEGER IV(LIV)
      REAL RD1(N), RD2(N), V(LV)
C
C     ***  NOTE -- V IS PASSED FOR POSSIBLE USE BY REVISED VERSIONS OF
C     ***  THIS ROUTINE.
C
      INTEGER PU
C
C  ***  IV SUBSCRIPTS  ***
C
      INTEGER COVPRT, NEEDHD, OUTLEV, PRUNIT, RDREQ, REGD
C
C/6
C     DATA COVPRT/14/, NEEDHD/36/, OUTLEV/19/, PRUNIT/21/, RDREQ/57/,
C    1     REGD/67/
C/7
      PARAMETER (COVPRT=14, NEEDHD=36, OUTLEV=19, PRUNIT=21, RDREQ=57,
     1           REGD=67)
C/
C
C+++++++++++++++++++++++++++++++  BODY  ++++++++++++++++++++++++++++++++
C
      PU = IV(PRUNIT)
      IF (PU .EQ. 0) GO TO 999
      IF (IV(COVPRT) .LT. 2) GO TO 999
      IF (IV(REGD) .LE. 0) GO TO 999
      IV(NEEDHD) = 1
      IF (IV(OUTLEV) .LT. 0) GO TO 30
      IF (IV(RDREQ) .EQ. 1 .OR. IV(RDREQ) .GE. 3) WRITE(PU,10) RD1
 10   FORMAT(/61H REGRESSION DIAGNOSTIC 1 = SQRT(G(I)**T * H(I)**-1 * G(
     1I)).../(10E12.3))
      IF (IV(RDREQ) .GE. 2) WRITE(PU,20) RD2
 20   FORMAT(/76H REGRESSION DIAGNOSTIC 2 = SQRT(G(I)**T * H(I)**-1 * H
     1* H(I)**-1 * G(I)).../(10E12.3))
      GO TO 999
C
 30   IF (IV(RDREQ) .EQ. 1 .OR. IV(RDREQ) .GE. 3) WRITE(PU,40) RD1
 40   FORMAT(/61H REGRESSION DIAGNOSTIC 1 = SQRT(G(I)**T * H(I)**-1 * G(
     1I)).../(6E12.3))
      IF (IV(RDREQ) .GE. 2) WRITE(PU,50) RD2
 50   FORMAT(/76H REGRESSION DIAGNOSTIC 2 = SQRT(G(I)**T * H(I)**-1 * H
     1* H(I)**-1 * G(I)).../(6E12.3))
C
 999  RETURN
C  ***  LAST CARD OF MLPRDP FOLLOWS  ***
      END
      SUBROUTINE NL2IT(D, DR, IV, LIV, LV, N, NB, NN, N1, P, R,
     1                 RD1, RD2, V, X)
C
C  ***  REVISED ITERATION DRIVER FOR NL2SOL  ***
C
      INTEGER LIV, LV, N, NB, NN, N1, P
      INTEGER IV(LIV)
      REAL D(P), DR(N,P,NN), R(N,NN), RD1(N,NN),
     1                 RD2(N,NN), V(LV), X(P)
C
C--------------------------  PARAMETER USAGE  --------------------------
C
C D........ SCALE VECTOR.
C DR....... DERIVATIVES OF R AT X.
C IV....... INTEGER VALUES ARRAY.
C LIV...... LENGTH OF IV... LIV MUST BE AT LEAST P + 80.
C LV....... LENGTH OF V...  LV  MUST BE AT LEAST 105 + P*(2*P+16).
C N........ NUMBER OF RESIDUALS IN ONE BLOCK OF R.
C NB....... NUMBER OF BLOCKS OF RESIDUALS.
C NN....... NUMBER OF BLOCKS OF R AND DR SUPPLIED THIS TIME.
C N1....... HIGHEST BLOCK INDEX FOR WHICH ABOVE ARE SUPPLIED THIS TIME.
C P........ NUMBER OF PARAMETERS (COMPONENTS OF X) BEING ESTIMATED.
C R........ RESIDUALS.
C RD1...... RD1(I) = SQRT(G(I)**T * H(I)**-1 * G(I)) ON OUTPUT WHEN
C        IV(RDREQ) = 1 OR 3.
C RD2...... RD2(I) = SQRT(G(I)**T * H(I)**-1 * H * H(I)**-1 * G(I)) ON
C        OUTPUT IF IV(RDREQ) = 2 OR 3.  NL2IT SETS IV(REGD) = 1 IF RD1
C        AND RD2 ARE SUCCESSFULLY COMPUTED, TO 0 IF NO ATTEMPT WAS MADE
C        TO COMPUTE THEM, AND TO -1 IF H (THE FINITE-DIFFERENCE HESSIAN)
C        WAS INDEFINITE.
C V........ FLOATING-POINT VALUES ARRAY.
C X........ PARAMETER VECTOR BEING ESTIMATED (INPUT = INITIAL GUESS,
C             OUTPUT = BEST VALUE FOUND).
C
C  ***  DISCUSSION  ***
C
C        THIS ROUTINE CARRIES OUT ITERATIONS FOR SOLVING NONLINEAR
C     LEAST SQUARES PROBLEMS.  WHEN NB = 1, IT IS SIMILAR TO NL2ITR
C     (WITH J = DR), EXCEPT THAT R(X) AND DR(X) NEED NOT BE INITIALIZED
C     WHEN NL2IT IS CALLED WITH IV(1) = 0 OR 12.  NL2IT ALSO ALLOWS
C     R AND DR TO BE SUPPLIED ROW-WISE -- JUST SET N = 1 AND NB = NO.
C     OF OBSERVATIONS.
C        ANOTHER NEW FEATURE IS THAT CALLING NL2IT WITH IV(1) = 13
C     CAUSES STORAGE ALLOCATION ONLY TO BE PERFORMED -- ON RETURN, SUCH
C     COMPONENTS AS IV(G) (THE FIRST SUBSCRIPT IN G OF THE GRADIENT)
C     AND IV(S) (THE FIRST SUBSCRIPT IN V OF THE S LOWER TRIANGLE OF
C     THE S MATRIX) WILL HAVE BEEN SET (UNLESS LIV OR LV IS TOO SMALL),
C     AND IV(1) WILL HAVE BEEN SET TO 14. CALLING NL2IT WITH IV(1) = 14
C     CAUSES EXECUTION OF THE ALGORITHM TO BEGIN UNDER THE ASSUMPTION
C     THAT STORAGE HAS BEEN ALLOCATED.
C
C ***  SUPPLYING R AND DR  ***
C
C        NL2IT USES IV AND V IN THE SAME WAY AS NL2SOL, WITH A SMALL
C     NUMBER OF OBVIOUS CHANGES.  ONE DIFFERENCE BETWEEN NL2IT AND
C     NL2ITR IS THAT INITIAL FUNCTION AND GRADIENT INFORMATION NEED NOT
C     BE SUPPLIED IN THE VERY FIRST CALL ON NL2IT, THE ONE WITH
C     IV(1) = 0 OR 12.  ANOTHER DIFFERENCE IS THAT NL2IT RETURNS WITH
C     IV(1) = -2 WHEN IT WANTS ANOTHER LOOK AT THE OLD JACOBIAN MATRIX
C     AND THE CURRENT RESIDUAL -- THE ONE CORRESPONDING TO X AND
C     IV(NFGCAL).  NOTE THAT IV(NFGCAL) = IV(7) CONTAINS THE VALUE THAT
C     IV(NFCALL) = IV(6) HAD WHEN THE CURRENT RESIDUAL WAS EVALUATED.
C     ALSO NOTE THAT THE VALUE OF X CORRESPONDING TO THE OLD JACOBIAN
C     MATRIX IS STORED IN V, STARTING AT V(IV(X0)) = V(IV(43)).
C        FUNCTION INFORMATION MUST ONLY BE SUPPLIED WHEN NL2IT RETURNS
C     WITH IV(1) = 1.  THIS TAKES THE FORM OF VALUES OF R(K,I,X) PASSED
C     IN R. IT IS POSSIBLE TO PASS ALL THESE VALUES AT ONCE, BUT THIS
C     IS NOT REQUIRED.  IN PRACTICE, R(K,I) AND DR(K,L,I) MUST
C     ALL CORRESPOND TO THE SAME RESIDUALS (K = 1(1)Q-1, L = 1(1)P).
C     THE CALLER CHOOSES HOW MANY R BLOCKS TO PASS AT ONCE AND SETS NN
C     TO THIS NUMBER.  SO THAT NL2IT WILL KNOW WHEN ALL R BLOCKS HAVE
C     BEEN PASSED, THE CALLER SETS N1 TO THE HIGHEST VALUE OF I FOR THE
C     CURRENT CALL.  FOR EXAMPLE, SUPPOSE N = 80 AND THAT R IS TO BE
C     PASSED IN 8 BLOCKS OF SIZE 10.  THE FOLLOWING CODE WOULD DO THE
C     JOB.
C
C         N = 10
C         NB = 8
C         NN = 1
C         DO 10 N1 = 1, NB
C              ***  COMPUTE THE NEXT 10 VALUES OF R  ***
C              CALL NL2IT(...,N,NB,NN,N1,...)
C      10      CONTINUE
C
C        THE SITUATION IS SIMILAR WHEN GRADIENT INFORMATION IS
C     REQUIRED, I.E., WHEN NL2IT RETURNS WITH IV(1) = 2 OR -2.  IT MAY
C     BE HELPFUL TO KNOW THAT NL2IT DOES NOT OVERWRITE R.  IT
C     SHOULD BE CLEAR THAT THE PARTIAL DERIVATIVE OF R(K,I,X) WITH
C     RESPECT TO X(L) IS TO BE STORED IN DR(K,L,I), L = 1(1)P.
C
C  ***  COVARIANCE MATRIX  ***
C
C        IF IV(COVREQ) OR IV(COVPRT) IS NONZERO, THEN A COVARIANCE
C     MATRIX IS COMPUTED AS DESCRIBED IN THE NL2SOL DOCUMENTATION.
C
C  ***  GENERAL  ***
C
C     CODED BY DAVID M. GAY.
C
C+++++++++++++++++++++++++++++  DECLARATIONS  ++++++++++++++++++++++++++
C
C  ***  INTRINSIC FUNCTIONS  ***
C/+
      INTEGER IABS, MAX0
C/
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      EXTERNAL DEFLT,  DUPD2, GLRIT, ITSUM, LINVRT, LSQRT,
     1         LTSQAR, LVMUL, NL2CVP, NL2RD, QRADR, VAXPY, VCOPY,
     2         VSCOPY
C
C DEFLT.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS.
C DUPD2...  UPDATES SCALE VECTOR D.
C GLRIT.... PERFORMS BASIC MINIMIZATION ALGORITHM.
C ITSUM.... PRINTS ITERATION SUMMARY, INFO ABOUT INITIAL AND FINAL X.
C LINVRT... INVERTS LOWER TRIANGULAR MATRIX.
C LSQRT.... COMPUTES CHOLESKY FACTOR OF (LOWER TRIANG. OF) SYM. MATRIX.
C LTSQAR... GIVEN LOWER TRIANG. MATRIX L, COMPUTES (L**T)*L.
C LVMUL.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX.
C NL2CVP... PRINTS COVARIANCE MATRIX.
C NL2RD.... COMPUTES REGRESSION DIAGNOSTICS.
C QRADR.... ADDS A NEW ROW TO QR DECOMPOSITION.
C VAXPY.... ADDS A MULTIPLE OF ONE VECTOR TO ANOTHER.
C VCOPY.... COPIES ONE VECTOR TO ANOTHER.
C VSCOPY... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR.
C
      INTEGER COV, G1, GI, I, IV1, JTOL1, K, L, LH,
     1        QTR1, RMAT1, YI, Y1
      REAL RKL, T
C
      REAL HALF, ZERO
C
C  ***  SUBSCRIPTS FOR IV AND V  ***
C
      INTEGER CNVCOD, COVMAT, COVPRT, COVREQ, DINIT, DTYPE, DTINIT,
     1        D0INIT, F, FDH, G, INITS, IPIVOT, IVNEED, JCN, JTOL, LMAT,
     2        MODE, NEXTIV, NEXTV, NFCALL, NFGCAL, NGCALL, NGCOV, QTR,
     3        RDREQ, REGD, RLIMIT, RMAT, TOOBIG, VNEED, Y
C
C  ***  IV SUBSCRIPT VALUES  ***
C
C/6
C     DATA CNVCOD/55/, COVMAT/26/, COVPRT/14/, COVREQ/15/, DTYPE/16/,
C    1     FDH/74/, G/28/, INITS/25/, IPIVOT/76/, IVNEED/3/, JCN/66/,
C    2     JTOL/59/, LMAT/42/, MODE/35/, NEXTIV/46/, NEXTV/47/,
C    3     NFCALL/6/, NFGCAL/7/, NGCALL/30/, NGCOV/53/, QTR/77/,
C    4     RMAT/78/, RDREQ/57/, REGD/67/, TOOBIG/2/, VNEED/4/, Y/48/
C/7
      PARAMETER (CNVCOD=55, COVMAT=26, COVPRT=14, COVREQ=15, DTYPE=16,
     1           FDH=74, G=28, INITS=25, IPIVOT=76, IVNEED=3, JCN=66,
     2           JTOL=59, LMAT=42, MODE=35, NEXTIV=46, NEXTV=47,
     3           NFCALL=6, NFGCAL=7, NGCALL=30, NGCOV=53, QTR=77,
     4           RMAT=78, RDREQ=57, REGD=67, TOOBIG=2, VNEED=4, Y=48)
C/
C
C  ***  V SUBSCRIPT VALUES  ***
C
C/6
C     DATA DINIT/38/, DTINIT/39/, D0INIT/40/, F/10/, RLIMIT/46/
C/7
      PARAMETER (DINIT=38, DTINIT=39, D0INIT=40, F=10, RLIMIT=46)
C/
C/6
C     DATA HALF/0.5E+0/, ZERO/0.E+0/
C/7
      PARAMETER (HALF=0.5E+0, ZERO=0.E+0)
C/
C
C+++++++++++++++++++++++++++++++  BODY  ++++++++++++++++++++++++++++++++
C
      LH = P * (P+1) / 2
      IF (IV(1) .EQ. 0) CALL DEFLT(1, IV, LIV, LV, V)
      IV1 = IV(1)
      IF (IV1 .EQ. 1) GO TO 60
      IF (IABS(IV1) .EQ. 2) GO TO 70
C
C  ***  FRESH START OR RESTART -- CHECK INPUT INTEGERS  ***
C
      IF (NN .LE. 0 .OR. P .LE. 0 .OR. N .LE. 0 .OR. NB .LE. 0 .OR.
     1     NN .GT. NB) GO TO 190
      IV(IVNEED) = IV(IVNEED) + P
      IV(VNEED) = IV(VNEED) + P*(P+13)/2
      IF (IV1 .EQ. 14) GO TO 10
      IF (IV1 .LT. 12) GO TO 20
      IF (IV1 .EQ. 12) IV(1) = 13
      G1 = 1
      Y1 = 1
      GO TO 30
C
 10   JTOL1 = IV(JTOL)
      IF (V(DINIT) .GE. ZERO) CALL VSCOPY(P, D, V(DINIT))
      IF (V(DTINIT) .GT. ZERO) CALL VSCOPY(P, V(JTOL1), V(DTINIT))
      I = JTOL1 + P
      IF (V(D0INIT) .GT. ZERO) CALL VSCOPY(P, V(I), V(D0INIT))
C
 20   G1 = IV(G)
      Y1 = IV(Y)
 30   CALL GLRIT(D, V(G1), IV, LIV, LV, P, P, V, X, V(Y1))
      IF (IV(1) - 2) 40, 50, 210
C
 40   V(F) = ZERO
      GO TO 999
C
 50   CALL VSCOPY(P, V(G1), ZERO)
      IF (IV(MODE) .GT. 0) GO TO 999
      CALL VSCOPY(P, V(Y1), ZERO)
      RMAT1 = IV(RMAT)
      CALL VSCOPY(LH, V(RMAT1), ZERO)
      QTR1 = IV(QTR)
      CALL VSCOPY(P, V(QTR1), ZERO)
      IV(REGD) = 0
      IF (IV(MODE) .EQ. 0) IV(1) = -2
      GO TO 999
C
C  ***  COMPUTE F(X)  ***
C
 60   T = SNRM2(NN*N,R,1)
      IF (T .GT. V(RLIMIT)) GO TO 180
      V(F) = V(F)  +  HALF * T**2
      IF (N1 .LT. NB) GO TO 999
      GO TO 20
C
C  ***  PROCESS NEW GRADIENT INFORMATION  ***
C
 70   IF (IV(NFGCAL) .LE. 0) GO TO 200
      IF (IV(1) .EQ. 2) GO TO 100
C
C  ***  COMPUTE Y (UNLESS THIS IS A FRESH START)  ***
C
      Y1 = IV(Y)
      DO 90 I = 1, NN
         YI = Y1
         DO 80 L = 1, P
            V(YI) = V(YI) + SDOT(N, DR(1,L,I),1, R(1,I),1)
            YI = YI + 1
 80         CONTINUE
 90      CONTINUE
      IF (N1 .GE. NB) IV(1) = 2
      GO TO 999
C
C  ***  COMPUTE GRADIENT INFORMATION  ***
C
 100  IF (IV(MODE) .GT. P) GO TO 230
      G1 = IV(G)
      IF (IV(MODE) .LE. 0) GO TO 130
C
C  ***  COMPUTE GRADIENT ONLY (FOR USE IN COVARIANCE COMPUTATION)  ***
C
      DO 120 I = 1, NN
         GI = G1
         DO 110 L = 1, P
            V(GI) = V(GI) + SDOT(N, R(1,I),1, DR(1,L,I),1)
            GI = GI + 1
 110        CONTINUE
 120     CONTINUE
      GO TO 170
C
C  ***  COMPUTE RMAT AND QTR, USING G AS SCRATCH  ***
C
 130  QTR1 = IV(QTR)
      RMAT1 = IV(RMAT)
      DO 160 L = 1, NN
         DO 150 K = 1, N
              GI = G1
              DO 140 I = 1, P
                   V(GI) = DR(K,I,L)
                   GI = GI + 1
 140               CONTINUE
              RKL = R(K,L)
              CALL QRADR(P, V(QTR1), V(RMAT1), V(G1), RKL)
 150          CONTINUE
 160     CONTINUE
C
C  ***  UPDATE D IF DESIRED  ***
C
      IF (IV(DTYPE) .LE. 0) GO TO 170
      IF (IV(MODE) .GT. 0) GO TO 170
         CALL DUPD2(D, DR, IV, LIV, LV, N, NB, NN, N1, P, V)
C
 170  IF (N1 .LT. NB) GO TO 999
      IF (IV(MODE) .GT. 0) GO TO 20
C
C  ***  COMPUTE G FROM RMAT AND QTR  ***
C
      CALL LVMUL(P, V(G1), V(RMAT1), V(QTR1))
      GO TO 20
C
C  ***  MISC. DETAILS  ***
C
C     ***  X IS OUT OF RANGE (OVERSIZE STEP)  ***
C
 180  IV(TOOBIG) = 1
      GO TO 20
C
C     ***  BAD N, NB, NN, OR P  ***
C
 190  IV(1) = 66
      GO TO 280
C
C     ***  DR (I.E., THE GRADIENT) CANNOT BE COMPUTED AT X  ***
C
 200  IV(1) = 65
      GO TO 280
C
C  ***  STORAGE ALLOCATION  ***
C
 210  IF (IV(1) .NE. 14) GO TO 220
      IV(IPIVOT) = IV(NEXTIV)
      IV(NEXTIV) = IV(IPIVOT) + P
      IV(Y) = IV(NEXTV)
      IV(G) = IV(Y) + P
      IV(JCN) = IV(G) + P
      IV(RMAT) = IV(JCN) + P
      IV(QTR) = IV(RMAT) + LH
      JTOL1 = IV(QTR) + P
      IV(JTOL) = JTOL1
      IV(NEXTV) = JTOL1 + 2*P
      IF (IV1 .EQ. 13) GO TO 999
      GO TO 10
C
C  ***  CONVERGENCE OBTAINED -- SEE WHETHER TO COMPUTE COVARIANCE  ***
C
 220  IF (IV(FDH) .LE. 0) GO TO 280
      IF (IV(COVMAT) .NE. 0 .AND. IV(REGD) .NE. 0) GO TO 280
C
C     ***  COMPUTE CHOLESKY FACTOR OF HESSIAN  ***
C
      L = IV(LMAT)
      K = IV(FDH)
      CALL LSQRT(1, P, V(L), V(K), I)
      IF (I .NE. 0) GO TO 270
C
C     ***  COMPUTE REGRESSION DIAGNOSTICS AND DEFAULT COVARIANCE IF
C          DESIRED  ***
C
      I = 0
      IF (IV(RDREQ) .GT. 0 .AND. IV(REGD) .EQ. 0) I = 1
      IF ((IABS(IV(COVREQ)) .EQ. 1) .OR. (IV(COVPRT) .EQ. 1 .AND.
     1     (IABS(IV(COVREQ)) .LE. 1))) I = I + 2
      IF (I .EQ. 0) GO TO 240
      IV(MODE) = P + I
      IV(NGCALL) = IV(NGCALL) + 1
      IV(NGCOV) = IV(NGCOV) + 1
      IV(NFGCAL) = IV(NFCALL) + IV(NGCOV)
      IV(CNVCOD) = IV(1)
      IV(1) = 2
      IF (I .LT. 2) GO TO 999
         L = IV(FDH)
         IV(FDH) = -L
         CALL VSCOPY(LH, V(L), ZERO)
         GO TO 999
C
 230  L = IV(LMAT)
      CALL NL2RD(DR, IV, V(L), LIV, LV, N, NN, P, R, RD1, RD2, V)
      IF (N1 .LT. NB) GO TO 999
      I = IV(MODE) - P
      IF ((I-2)**2 .EQ. 1) IV(REGD) = 1
C
C     ***  FINISH COMPUTING COVARIANCE MATRIX = INVERSE OF F.D. HESSIAN.
C
      IV(1) = IV(CNVCOD)
      IV(MODE) = 0
      COV = L
      IF (I .LT. 2) GO TO 240
         COV = IABS(IV(FDH))
         IV(FDH) = 0
 240  IV(CNVCOD) = 0
C
      IF (IV(COVMAT) .NE. 0) GO TO 280
      IF (IV(COVPRT) .EQ. 0 .AND. IV(COVREQ) .EQ. 0) GO TO 280
      IF (I .GE. 2) GO TO 250
      CALL LINVRT(P, V(COV), V(COV))
      CALL LTSQAR(P, V(COV), V(COV))
 250  T = V(F) / (HALF * FLOAT(MAX0(1,N*NB-P)))
      L = COV + LH - 1
      DO 260 I = COV, L
 260     V(I) = T * V(I)
      IV(COVMAT) = COV
      GO TO 280
C
C  ***  COME HERE FOR INDEFINITE FINITE-DIFFERENCE HESSIAN  ***
C
 270  IV(COVMAT) = -1
      IV(REGD) = -1
C
C  ***  PRINT SUMMARY OF FINAL ITERATION AND OTHER REQUESTED ITEMS  ***
C
 280  CALL ITSUM(D, V(G1), IV, LIV, LV, P, V, X)
      CALL NL2CVP(IV, LIV, LV, P, V)
C
 999  RETURN
C  ***  LAST CARD OF NL2IT FOLLOWS  ***
      END
      SUBROUTINE LINVRT(N, LIN, L)
C
C  ***  COMPUTE  LIN = L**-1,  BOTH  N X N  LOWER TRIANG. STORED   ***
C  ***  COMPACTLY BY ROWS.  LIN AND L MAY SHARE THE SAME STORAGE.  ***
C
C***REVISION HISTORY  (YYMMDD)
C   000330  Modified array declarations.  (JEC)
C
C  ***  PARAMETERS  ***
C
      INTEGER N
      REAL L(*), LIN(*)
C     DIMENSION L(N*(N+1)/2), LIN(N*(N+1)/2)
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER I, II, IM1, JJ, J0, J1, K, K0, NP1
      REAL ONE, T, ZERO
C/6
C     DATA ONE/1.E+0/, ZERO/0.E+0/
C/7
      PARAMETER (ONE=1.E+0, ZERO=0.E+0)
C/
C
C  ***  BODY  ***
C
      NP1 = N + 1
      J0 = N*(NP1)/2
      DO 30 II = 1, N
         I = NP1 - II
         LIN(J0) = ONE/L(J0)
         IF (I .LE. 1) GO TO 999
         J1 = J0
         IM1 = I - 1
         DO 20 JJ = 1, IM1
              T = ZERO
              J0 = J1
              K0 = J1 - JJ
              DO 10 K = 1, JJ
                   T = T - L(K0)*LIN(J0)
                   J0 = J0 - 1
                   K0 = K0 + K - I
 10                CONTINUE
              LIN(J0) = T/L(K0)
 20           CONTINUE
         J0 = J0 - 1
 30      CONTINUE
 999  RETURN
C  ***  LAST CARD OF LINVRT FOLLOWS  ***
      END
      SUBROUTINE VDFLT(ALG, LV, V)
C
C  ***  SUPPLY ***SOL (VERSION 2.3) DEFAULT VALUES TO V  ***
C
C  ***  ALG = 1 MEANS REGRESSION CONSTANTS.
C  ***  ALG = 2 MEANS GENERAL UNCONSTRAINED OPTIMIZATION CONSTANTS.
C
      INTEGER ALG, LV
      REAL V(LV)
C/+
      REAL AMAX1
C/
C
      REAL MACHEP, MEPCRT, ONE, SQTEPS, THREE
C
C  ***  SUBSCRIPTS FOR V  ***
C
      INTEGER AFCTOL, BIAS, COSMIN, DECFAC, DELTA0, DFAC, DINIT, DLTFDC,
     1        DLTFDJ, DTINIT, D0INIT, EPSLON, ETA0, FUZZ, HUBERC,
     2        INCFAC, LMAX0, LMAXS, PHMNFC, PHMXFC, RDFCMN, RDFCMX,
     3        RFCTOL, RLIMIT, RSPTOL, SCTOL, SIGMIN, TUNER1, TUNER2,
     4        TUNER3, TUNER4, TUNER5, XCTOL, XFTOL
C
C/6
      DATA ONE/1.E+0/, THREE/3.E+0/
C/7
C     PARAMETER (ONE=1.E+0, THREE=3.E+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           E0INIT=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 = R1MACH(4)
      V(AFCTOL) = 1.E-20
      IF (MACHEP .GT. 1.E-10) V(AFCTOL) = MACHEP**2
      V(DECFAC) = 0.5E+0
      SQTEPS = SQRT(R1MACH(4))
      V(DFAC) = 0.6E+0
      V(DELTA0) = SQTEPS
      V(DTINIT) = 1.E-6
      MEPCRT = MACHEP ** (ONE/THREE)
      V(D0INIT) = 1.E+0
      V(EPSLON) = 0.1E+0
      V(INCFAC) = 2.E+0
      V(LMAX0) = 1.E+0
      V(LMAXS) = 1.E+0
      V(PHMNFC) = -0.1E+0
      V(PHMXFC) = 0.1E+0
      V(RDFCMN) = 0.1E+0
      V(RDFCMX) = 4.E+0
      V(RFCTOL) = AMAX1(1.E-10, MEPCRT**2)
      V(SCTOL) = V(RFCTOL)
      V(TUNER1) = 0.1E+0
      V(TUNER2) = 1.E-4
      V(TUNER3) = 0.75E+0
      V(TUNER4) = 0.5E+0
      V(TUNER5) = 0.75E+0
      V(XCTOL) = SQTEPS
      V(XFTOL) = 1.E+2 * MACHEP
C
      IF (ALG .GE. 2) GO TO 10
C
C  ***  REGRESSION  VALUES
C
      V(COSMIN) = AMAX1(1.E-6, 1.E+2 * MACHEP)
      V(DINIT) = 0.E+0
      V(DLTFDC) = MEPCRT
      V(DLTFDJ) = SQTEPS
      V(FUZZ) = 1.5E+0
      V(HUBERC) = 0.7E+0
      V(RLIMIT) = SQRT(R1MACH(2))*16.
      V(RSPTOL) = 1.E-2
      V(SIGMIN) = 1.E-4
      GO TO 999
C
C  ***  GENERAL OPTIMIZATION VALUES
C
 10   V(BIAS) = 0.8E+0
      V(DINIT) = -1.0E+0
      V(ETA0) = 1.0E+3 * MACHEP
C
 999  RETURN
C  ***  LAST CARD OF VDFLT FOLLOWS  ***
      END
      SUBROUTINE NL2RD(DR, IV, L, LIV, LV, N, NN, P, R, RD1, RD2, V)
C
C  ***  COMPUTE REGRESSION DIAGNOSTIC AND DEFAULT COVARIANCE MATRIX FOR
C       NL2IT  ***
C
C  ***  PARAMETERS  ***
C
      INTEGER LIV, LV, N, NN, P
      INTEGER IV(LIV)
C/6
C     REAL DR(N,P,NN), L(1), R(N,NN), RD1(N,NN),
C    1             RD2(N,NN), V(LV)
C     DIMENSI0N L(P*(P+1)/2)
C/7
      REAL DR(N,P,NN), L(P*(P+1)/2), R(N,NN), RD1(N,NN),
     1             RD2(N,NN), V(LV)
C/
C
C  ***  CODED BY DAVID M. GAY (WINTER 1982)  ***
C
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      EXTERNAL  LIVMUL, LITVMU, LTVMUL, OPRDS, VAXPY, VCOPY,
     1         VSCOPY
C
C  ***  LOCAL VARIABLES  ***
C
      LOGICAL D1, D2
      INTEGER COV, I, J, K, LH, M, STEP1
      REAL A, S, T
C
C  ***  CONSTANTS  ***
C
      REAL NEGONE, ONE, ONEV(1), ZERO
C
C  ***  INTRINSIC FUNCTIONS  ***
C/+
      INTEGER IABS
      REAL  SQRT
C/
C
C  ***  IV SUBSCRIPTS  ***
C
      INTEGER D, FDH, MODE, RDREQ, STEP
C/6
C     DATA D/27/, FDH/74/, MODE/35/, RDREQ/57/, STEP/40/
C/7
      PARAMETER (D=27, FDH=74, MODE=35, RDREQ=57, STEP=40)
C/
C/6
C     DATA NEGONE/-1.E+0/, ONE/1.E+0/, ZERO/0.E+0/
C/7
      PARAMETER (NEGONE=-1.E+0, ONE=1.E+0, ZERO=0.E+0)
C/
      DATA ONEV(1)/1.E+0/
C
C++++++++++++++++++++++++++++++++  BODY  +++++++++++++++++++++++++++++++
C
      STEP1 = IV(STEP)
      I = IV(RDREQ)
      IF (I .LE. 0) GO TO 40
      D1 = I .EQ. 1 .OR. I .GE. 3
      D2 = I .GE. 2
      IF (D1) CALL VSCOPY(NN*N, RD1, NEGONE)
      IF (D2) CALL VSCOPY(NN*N, RD2, NEGONE)
      DO 30 K = 1, NN
         DO 20 I = 1, N
            A = R(I,K)**2
            M = STEP1
            DO 10 J = 1, P
               V(M) = DR(I,J,K)
               M = M + 1
 10            CONTINUE
            CALL LIVMUL(P, V(STEP1), L, V(STEP1))
            S = SDOT(P,V(STEP1),1,V(STEP1),1)
            T = ONE - S
            IF (T .LE. ZERO) GO TO 20
            A = A * S / T
            IF (D1) RD1(I,K) =  SQRT(A)
            IF (D2) RD2(I,K) =  SQRT(A / T)
 20         CONTINUE
 30      CONTINUE
C
 40   IF (IV(MODE) - P .LT. 2) GO TO 999
C
C  ***  COMPUTE DEFAULT COVARIANCE MATRIX  ***
C
      LH = P*(P+1)/2
      COV = IABS(IV(FDH))
      DO 70 K = 1, NN
         DO 60 I = 1, N
            M = STEP1
            DO 50 J = 1, P
               V(M) = DR(I,J,K)
               M = M + 1
 50            CONTINUE
            CALL LIVMUL(P, V(STEP1), L, V(STEP1))
            CALL LITVMU(P, V(STEP1), L, V(STEP1))
            CALL OPRDS(1, LH, P, V(COV), ONEV, V(STEP1), V(STEP1))
 60         CONTINUE
 70      CONTINUE
C
 999  RETURN
C  ***  LAST CARD OF NL2RD FOLLOWS  ***
      END
      SUBROUTINE LTSQAR(N, A, L)
C
C  ***  SET A TO LOWER TRIANGLE OF (L**T) * L  ***
C
C  ***  L = N X N LOWER TRIANG. MATRIX STORED ROWWISE.  ***
C  ***  A IS ALSO STORED ROWWISE AND MAY SHARE STORAGE WITH L.  ***
C
C***REVISION HISTORY  (YYMMDD)
C   000330  Modified array declarations.  (JEC)
C
      INTEGER N
      REAL A(*), L(*)
C     DIMENSION A(N*(N+1)/2), L(N*(N+1)/2)
C
      INTEGER I, II, IIM1, I1, J, K, M
      REAL LII, LJ
C
      II = 0
      DO 50 I = 1, N
         I1 = II + 1
         II = II + I
         M = 1
         IF (I .EQ. 1) GO TO 30
         IIM1 = II - 1
         DO 20 J = I1, IIM1
              LJ = L(J)
              DO 10 K = I1, J
                   A(M) = A(M) + LJ*L(K)
                   M = M + 1
 10                CONTINUE
 20           CONTINUE
 30      LII = L(II)
         DO 40 J = I1, II
 40           A(J) = LII * L(J)
 50      CONTINUE
C
 999  RETURN
C  ***  LAST CARD OF LTSQAR FOLLOWS  ***
      END
      INTEGER FUNCTION I1MACH(I)
C***BEGIN PROLOGUE  I1MACH
C***DATE WRITTEN   750101   (YYMMDD)
C***REVISION DATE  840405   (YYMMDD)
C***CATEGORY NO.  R1
C***KEYWORDS  MACHINE CONSTANTS
C***AUTHOR  FOX, P. A., (BELL LABS)
C           HALL, A. D., (BELL LABS)
C           SCHRYER, N. L., (BELL LABS)
C***PURPOSE  Returns integer machine dependent constants
C***DESCRIPTION
C
C     This is the CMLIB version of I1MACH, the integer machine
C     constants subroutine originally developed for the PORT library.
C
C     I1MACH can be used to obtain machine-dependent parameters
C     for the local machine environment.  It is a function
C     subroutine with one (input) argument, and can be called
C     as follows, for example
C
C          K = I1MACH(I)
C
C     where I=1,...,16.  The (output) value of K above is
C     determined by the (input) value of I.  The results for
C     various values of I are discussed below.
C
C  I/O unit numbers.
C    I1MACH( 1) = the standard input unit.
C    I1MACH( 2) = the standard output unit.
C    I1MACH( 3) = the standard punch unit.
C    I1MACH( 4) = the standard error message unit.
C
C  Words.
C    I1MACH( 5) = the number of bits per integer storage unit.
C    I1MACH( 6) = the number of characters per integer storage unit.
C
C  Integers.
C    assume integers are represented in the S-digit, base-A form
C
C               sign ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) )
C
C               where 0 .LE. X(I) .LT. A for I=0,...,S-1.
C    I1MACH( 7) = A, the base.
C    I1MACH( 8) = S, the number of base-A digits.
C    I1MACH( 9) = A**S - 1, the largest magnitude.
C
C  Floating-Point Numbers.
C    Assume floating-point numbers are represented in the T-digit,
C    base-B form
C               sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
C
C               where 0 .LE. X(I) .LT. B for I=1,...,T,
C               0 .LT. X(1), and EMIN .LE. E .LE. EMAX.
C    I1MACH(10) = B, the base.
C
C  Single-Precision
C    I1MACH(11) = T, the number of base-B digits.
C    I1MACH(12) = EMIN, the smallest exponent E.
C    I1MACH(13) = EMAX, the largest exponent E.
C
C  Double-Precision
C    I1MACH(14) = T, the number of base-B digits.
C    I1MACH(15) = EMIN, the smallest exponent E.
C    I1MACH(16) = EMAX, the largest exponent E.
C
C  To alter this function for a particular environment,
C  the desired set of DATA statements should be activated by
C  removing the C from column 1.  Also, the values of
C  I1MACH(1) - I1MACH(4) should be checked for consistency
C  with the local operating system.
C***REFERENCES  FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A
C                 PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL
C                 SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188.
C***ROUTINES CALLED  (NONE)
C***END PROLOGUE  I1MACH
C
      INTEGER IMACH(16),OUTPUT
      EQUIVALENCE (IMACH(4),OUTPUT)
C
C     MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T
C     3B SERIES, MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T
C     PC 7300), SUN SPARCSTATIONS, SILICON GRAPHCS WORKSTATIONS, HP
C     9000 WORKSTATIONS, IBM RS/6000 WORKSTATIONS, AND 8087 BASED 
C     MICROS (E.G. IBM PC AND AT&T 6300).
C
C === MACHINE = ATT.3B
C === MACHINE = ATT.6300
C === MACHINE = ATT.7300
C === MACHINE = HP.9000
C === MACHINE = IBM.PC
C === MACHINE = IBM.RS6000
C === MACHINE = IEEE.MOST-SIG-BYTE-FIRST
C === MACHINE = IEEE.LEAST-SIG-BYTE-FIRST
C === MACHINE = SGI
C === MACHINE = SUN
C === MACHINE = 68000
C === MACHINE = 8087
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    7 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   32 /
C      DATA IMACH( 6) /    4 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   31 /
C      DATA IMACH( 9) / 2147483647 /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   24 /
C      DATA IMACH(12) / -125 /
C      DATA IMACH(13) /  128 /
C      DATA IMACH(14) /   53 /
C      DATA IMACH(15) / -1021 /
C      DATA IMACH(16) /  1024 /
C
C     MACHINE CONSTANTS FOR SUN WORKSTATIONS.  f77 WITH -r8 OPTION.
C
C === MACHINE = SUN.F77-WITH-R8-OPTION
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    6 /
C      DATA IMACH( 4) /    0 /
C      DATA IMACH( 5) /   32 /
C      DATA IMACH( 6) /    4 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   31 /
C      DATA IMACH( 9) / 2147483647 /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   53 /
C      DATA IMACH(12) / -1021 /
C      DATA IMACH(13) /  1024 /
C      DATA IMACH(14) /   113 /
C      DATA IMACH(15) / -16382 /
C      DATA IMACH(16) /  16384 /
C
C     MACHINE CONSTANTS FOR SGI Origin 2000 with -r8 -d16 options.
C
C === MACHINE = SGI.ORIGIN.F77-WITH-R8-D16-OPTIONS
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    6 /
C      DATA IMACH( 4) /    0 /
C      DATA IMACH( 5) /   32 /
C      DATA IMACH( 6) /    4 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   31 /
C      DATA IMACH( 9) / 2147483647 /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   53 /
C      DATA IMACH(12) / -1022 /
C      DATA IMACH(13) /  1024 /
C      DATA IMACH(14) /   107 /
C      DATA IMACH(15) /  -916 /
C      DATA IMACH(16) /  1025 /
C
C     MACHINE CONSTANTS FOR IBM RS/6000 WORKSTATIONS WITH -qautodbl=dblpad.
C
C === MACHINE = IBM.RS6000.XLF-WITH-AUTODBL-OPTION
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    6 /
C      DATA IMACH( 4) /    0 /
C      DATA IMACH( 5) /   32 /
C      DATA IMACH( 6) /    4 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   31 /
C      DATA IMACH( 9) / 2147483647 /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   53 /
C      DATA IMACH(12) / -1021 /
C      DATA IMACH(13) /  1024 /
C      DATA IMACH(14) /   114 /
C      DATA IMACH(15) / -1021 /
C      DATA IMACH(16) /  1024 /
C
C     MACHINE CONSTANTS FOR AMDAHL MACHINES.
C
C === MACHINE = AMDAHL
C      DATA IMACH( 1) /   5 /
C      DATA IMACH( 2) /   6 /
C      DATA IMACH( 3) /   7 /
C      DATA IMACH( 4) /   6 /
C      DATA IMACH( 5) /  32 /
C      DATA IMACH( 6) /   4 /
C      DATA IMACH( 7) /   2 /
C      DATA IMACH( 8) /  31 /
C      DATA IMACH( 9) / 2147483647 /
C      DATA IMACH(10) /  16 /
C      DATA IMACH(11) /   6 /
C      DATA IMACH(12) / -64 /
C      DATA IMACH(13) /  63 /
C      DATA IMACH(14) /  14 /
C      DATA IMACH(15) / -64 /
C      DATA IMACH(16) /  63 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.
C
C === MACHINE = BURROUGHS.1700
C      DATA IMACH( 1) /    7 /
C      DATA IMACH( 2) /    2 /
C      DATA IMACH( 3) /    2 /
C      DATA IMACH( 4) /    2 /
C      DATA IMACH( 5) /   36 /
C      DATA IMACH( 6) /    4 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   33 /
C      DATA IMACH( 9) / Z1FFFFFFFF /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   24 /
C      DATA IMACH(12) / -256 /
C      DATA IMACH(13) /  255 /
C      DATA IMACH(14) /   60 /
C      DATA IMACH(15) / -256 /
C      DATA IMACH(16) /  255 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM.
C
C === MACHINE = BURROUGHS.5700
C      DATA IMACH( 1) /   5 /
C      DATA IMACH( 2) /   6 /
C      DATA IMACH( 3) /   7 /
C      DATA IMACH( 4) /   6 /
C      DATA IMACH( 5) /  48 /
C      DATA IMACH( 6) /   6 /
C      DATA IMACH( 7) /   2 /
C      DATA IMACH( 8) /  39 /
C      DATA IMACH( 9) / O0007777777777777 /
C      DATA IMACH(10) /   8 /
C      DATA IMACH(11) /  13 /
C      DATA IMACH(12) / -50 /
C      DATA IMACH(13) /  76 /
C      DATA IMACH(14) /  26 /
C      DATA IMACH(15) / -50 /
C      DATA IMACH(16) /  76 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS.
C
C === MACHINE = BURROUGHS.6700
C === MACHINE = BURROUGHS.7700
C      DATA IMACH( 1) /   5 /
C      DATA IMACH( 2) /   6 /
C      DATA IMACH( 3) /   7 /
C      DATA IMACH( 4) /   6 /
C      DATA IMACH( 5) /  48 /
C      DATA IMACH( 6) /   6 /
C      DATA IMACH( 7) /   2 /
C      DATA IMACH( 8) /  39 /
C      DATA IMACH( 9) / O0007777777777777 /
C      DATA IMACH(10) /   8 /
C      DATA IMACH(11) /  13 /
C      DATA IMACH(12) / -50 /
C      DATA IMACH(13) /  76 /
C      DATA IMACH(14) /  26 /
C      DATA IMACH(15) / -32754 /
C      DATA IMACH(16) /  32780 /
C
C     MACHINE CONSTANTS FOR THE CONVEX C1, C2, C3 SERIES
C
C === MACHINE = CONVEX
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    0 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   32 /
C      DATA IMACH( 6) /    4 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   31 /
C      DATA IMACH( 9) / 2147483647 /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   24 /
C      DATA IMACH(12) / -127 /
C      DATA IMACH(13) /  127 /
C      DATA IMACH(14) /   53 /
C      DATA IMACH(15) / -1023 /
C      DATA IMACH(16) /  1023 /
C
C     MACHINE CONSTANTS FOR THE CONVEX C1, C2, C3 SERIES (NATIVE MODE)
C     WITH -P8 OPTION
C
C === MACHINE = CONVEX.P8
C      DATA IMACH( 1) /     5 /
C      DATA IMACH( 2) /     6 /
C      DATA IMACH( 3) /     0 /
C      DATA IMACH( 4) /     6 /
C      DATA IMACH( 5) /    64 /
C      DATA IMACH( 6) /     8 /
C      DATA IMACH( 7) /     2 /
C      DATA IMACH( 8) /    63 /
C      DATA IMACH( 9) / 9223372036854775807 /
C      DATA IMACH(10) /     2 /
C      DATA IMACH(11) /    53 /
C      DATA IMACH(12) / -1023 /
C      DATA IMACH(13) /  1023 /
C      DATA IMACH(14) /   112 /
C      DATA IMACH(15) / -16383 /
C      DATA IMACH(16) /  16383 /
C
C     MACHINE CONSTANTS FOR THE CONVEX C1, C2, C3 SERIES (IEEE MODE)
C
C === MACHINE = CONVEX.IEEE
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    0 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   32 /
C      DATA IMACH( 6) /    4 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   31 /
C      DATA IMACH( 9) / 2147483647 /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   24 /
C      DATA IMACH(12) / -125 /
C      DATA IMACH(13) /  128 /
C      DATA IMACH(14) /   53 /
C      DATA IMACH(15) / -1021 /
C      DATA IMACH(16) /  1024 /
C
C     MACHINE CONSTANTS FOR THE CONVEX C1, C2, C3 SERIES(IEEE MODE)
C     WITH -P8 OPTION
C
C === MACHINE = CONVEX.IEEE.P8
C      DATA IMACH( 1) /     5 /
C      DATA IMACH( 2) /     6 /
C      DATA IMACH( 3) /     0 /
C      DATA IMACH( 4) /     6 /
C      DATA IMACH( 5) /    32 /
C      DATA IMACH( 6) /     4 /
C      DATA IMACH( 7) /     2 /
C      DATA IMACH( 8) /    31 /
C      DATA IMACH( 9) / 2147483647 /
C      DATA IMACH(10) /     2 /
C      DATA IMACH(11) /    53 /
C      DATA IMACH(12) / -1021 /
C      DATA IMACH(13) /  1024 /
C      DATA IMACH(14) /    53 /
C      DATA IMACH(15) / -1021 /
C      DATA IMACH(16) /  1024 /
C
C     MACHINE CONSTANTS FOR THE CYBER 170/180 SERIES USING NOS (FTN5).
C
C === MACHINE = CYBER.170.NOS
C === MACHINE = CYBER.180.NOS
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    7 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   60 /
C      DATA IMACH( 6) /   10 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   48 /
C      DATA IMACH( 9) / O"00007777777777777777" /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   48 /
C      DATA IMACH(12) / -974 /
C      DATA IMACH(13) / 1070 /
C      DATA IMACH(14) /   96 /
C      DATA IMACH(15) / -927 /
C      DATA IMACH(16) / 1070 /
C
C     MACHINE CONSTANTS FOR THE CDC 180 SERIES USING NOS/VE
C
C === MACHINE = CYBER.180.NOS/VE
C      DATA IMACH( 1) /     5 /
C      DATA IMACH( 2) /     6 /
C      DATA IMACH( 3) /     7 /
C      DATA IMACH( 4) /     6 /
C      DATA IMACH( 5) /    64 /
C      DATA IMACH( 6) /     8 /
C      DATA IMACH( 7) /     2 /
C      DATA IMACH( 8) /    63 /
C      DATA IMACH( 9) / 9223372036854775807 /
C      DATA IMACH(10) /     2 /
C      DATA IMACH(11) /    47 /
C      DATA IMACH(12) / -4095 /
C      DATA IMACH(13) /  4094 /
C      DATA IMACH(14) /    94 /
C      DATA IMACH(15) / -4095 /
C      DATA IMACH(16) /  4094 /
C
C     MACHINE CONSTANTS FOR THE CYBER 205
C
C === MACHINE = CYBER.205
C      DATA IMACH( 1) /      5 /
C      DATA IMACH( 2) /      6 /
C      DATA IMACH( 3) /      7 /
C      DATA IMACH( 4) /      6 /
C      DATA IMACH( 5) /     64 /
C      DATA IMACH( 6) /      8 /
C      DATA IMACH( 7) /      2 /
C      DATA IMACH( 8) /     47 /
C      DATA IMACH( 9) / X'00007FFFFFFFFFFF' /
C      DATA IMACH(10) /      2 /
C      DATA IMACH(11) /     47 /
C      DATA IMACH(12) / -28625 /
C      DATA IMACH(13) /  28718 /
C      DATA IMACH(14) /     94 /
C      DATA IMACH(15) / -28625 /
C      DATA IMACH(16) /  28718 /
C
C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES.
C
C === MACHINE = CDC.6000
C === MACHINE = CDC.7000
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    7 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   60 /
C      DATA IMACH( 6) /   10 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   48 /
C      DATA IMACH( 9) / 00007777777777777777B /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   48 /
C      DATA IMACH(12) / -974 /
C      DATA IMACH(13) / 1070 /
C      DATA IMACH(14) /   96 /
C      DATA IMACH(15) / -927 /
C      DATA IMACH(16) / 1070 /
C
C     MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3.
C     USING THE 46 BIT INTEGER COMPILER OPTION
C
C === MACHINE = CRAY.46-BIT-INTEGER
C      DATA IMACH( 1) /     5 /
C      DATA IMACH( 2) /     6 /
C      DATA IMACH( 3) /   102 /
C      DATA IMACH( 4) /     6 /
C      DATA IMACH( 5) /    64 /
C      DATA IMACH( 6) /     8 /
C      DATA IMACH( 7) /     2 /
C      DATA IMACH( 8) /    46 /
C      DATA IMACH( 9) /  777777777777777777777B /
C      DATA IMACH(10) /     2 /
C      DATA IMACH(11) /    47 /
C      DATA IMACH(12) / -8189 /
C      DATA IMACH(13) /  8190 /
C      DATA IMACH(14) /    94 /
C      DATA IMACH(15) / -8099 /
C      DATA IMACH(16) /  8190 /
C
C     MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3.
C     USING THE 64 BIT INTEGER COMPILER OPTION
C
C === MACHINE = CRAY.64-BIT-INTEGER
C      DATA IMACH( 1) /     5 /
C      DATA IMACH( 2) /     6 /
C      DATA IMACH( 3) /   102 /
C      DATA IMACH( 4) /     6 /
C      DATA IMACH( 5) /    64 /
C      DATA IMACH( 6) /     8 /
C      DATA IMACH( 7) /     2 /
C      DATA IMACH( 8) /    63 /
C      DATA IMACH( 9) /  777777777777777777777B /
C      DATA IMACH(10) /     2 /
C      DATA IMACH(11) /    47 /
C      DATA IMACH(12) / -8189 /
C      DATA IMACH(13) /  8190 /
C      DATA IMACH(14) /    94 /
C      DATA IMACH(15) / -8099 /
C      DATA IMACH(16) /  8190 /C
C     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200
C
C === MACHINE = DATA_GENERAL.ECLIPSE.S/200
C      DATA IMACH( 1) /   11 /
C      DATA IMACH( 2) /   12 /
C      DATA IMACH( 3) /    8 /
C      DATA IMACH( 4) /   10 /
C      DATA IMACH( 5) /   16 /
C      DATA IMACH( 6) /    2 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   15 /
C      DATA IMACH( 9) /32767 /
C      DATA IMACH(10) /   16 /
C      DATA IMACH(11) /    6 /
C      DATA IMACH(12) /  -64 /
C      DATA IMACH(13) /   63 /
C      DATA IMACH(14) /   14 /
C      DATA IMACH(15) /  -64 /
C      DATA IMACH(16) /   63 /
C
C     ELXSI  6400
C
C === MACHINE = ELSXI.6400
C      DATA IMACH( 1) /     5 /
C      DATA IMACH( 2) /     6 /
C      DATA IMACH( 3) /     6 /
C      DATA IMACH( 4) /     6 /
C      DATA IMACH( 5) /    32 /
C      DATA IMACH( 6) /     4 /
C      DATA IMACH( 7) /     2 /
C      DATA IMACH( 8) /    32 /
C      DATA IMACH( 9) / 2147483647 /
C      DATA IMACH(10) /     2 /
C      DATA IMACH(11) /    24 /
C      DATA IMACH(12) /  -126 /
C      DATA IMACH(13) /   127 /
C      DATA IMACH(14) /    53 /
C      DATA IMACH(15) / -1022 /
C      DATA IMACH(16) /  1023 /
C
C     MACHINE CONSTANTS FOR THE HARRIS 220
C     MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7
C
C === MACHINE = HARRIS.220
C === MACHINE = HARRIS.SLASH6
C === MACHINE = HARRIS.SLASH7
C      DATA IMACH( 1) /       5 /
C      DATA IMACH( 2) /       6 /
C      DATA IMACH( 3) /       0 /
C      DATA IMACH( 4) /       6 /
C      DATA IMACH( 5) /      24 /
C      DATA IMACH( 6) /       3 /
C      DATA IMACH( 7) /       2 /
C      DATA IMACH( 8) /      23 /
C      DATA IMACH( 9) / 8388607 /
C      DATA IMACH(10) /       2 /
C      DATA IMACH(11) /      23 /
C      DATA IMACH(12) /    -127 /
C      DATA IMACH(13) /     127 /
C      DATA IMACH(14) /      38 /
C      DATA IMACH(15) /    -127 /
C      DATA IMACH(16) /     127 /
C
C     MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES.
C     MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES.
C
C === MACHINE = HONEYWELL.600/6000
C === MACHINE = HONEYWELL.DPS.8/70
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /   43 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   36 /
C      DATA IMACH( 6) /    4 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   35 /
C      DATA IMACH( 9) / O377777777777 /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   27 /
C      DATA IMACH(12) / -127 /
C      DATA IMACH(13) /  127 /
C      DATA IMACH(14) /   63 /
C      DATA IMACH(15) / -127 /
C      DATA IMACH(16) /  127 /
C
C     MACHINE CONSTANTS FOR THE HP 2100
C     3 WORD DOUBLE PRECISION OPTION WITH FTN4
C
C === MACHINE = HP.2100.3_WORD_DP
C      DATA IMACH(1) /      5/
C      DATA IMACH(2) /      6 /
C      DATA IMACH(3) /      4 /
C      DATA IMACH(4) /      1 /
C      DATA IMACH(5) /     16 /
C      DATA IMACH(6) /      2 /
C      DATA IMACH(7) /      2 /
C      DATA IMACH(8) /     15 /
C      DATA IMACH(9) /  32767 /
C      DATA IMACH(10)/      2 /
C      DATA IMACH(11)/     23 /
C      DATA IMACH(12)/   -128 /
C      DATA IMACH(13)/    127 /
C      DATA IMACH(14)/     39 /
C      DATA IMACH(15)/   -128 /
C      DATA IMACH(16)/    127 /
C
C     MACHINE CONSTANTS FOR THE HP 2100
C     4 WORD DOUBLE PRECISION OPTION WITH FTN4
C
C === MACHINE = HP.2100.4_WORD_DP
C      DATA IMACH(1) /      5 /
C      DATA IMACH(2) /      6 /
C      DATA IMACH(3) /      4 /
C      DATA IMACH(4) /      1 /
C      DATA IMACH(5) /     16 /
C      DATA IMACH(6) /      2 /
C      DATA IMACH(7) /      2 /
C      DATA IMACH(8) /     15 /
C      DATA IMACH(9) /  32767 /
C      DATA IMACH(10)/      2 /
C      DATA IMACH(11)/     23 /
C      DATA IMACH(12)/   -128 /
C      DATA IMACH(13)/    127 /
C      DATA IMACH(14)/     55 /
C      DATA IMACH(15)/   -128 /
C      DATA IMACH(16)/    127 /
C
C     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
C     THE XEROX SIGMA 5/7/9 AND THE SEL SYSTEMS 85/86 AND
C     THE INTERDATA 3230 AND INTERDATA 7/32.
C
C === MACHINE = IBM.360
C === MACHINE = IBM.370
C === MACHINE = XEROX.SIGMA.5
C === MACHINE = XEROX.SIGMA.7
C === MACHINE = XEROX.SIGMA.9
C === MACHINE = SEL.85
C === MACHINE = SEL.86
C === MACHINE = INTERDATA.3230
C === MACHINE = INTERDATA.7/32
C      DATA IMACH( 1) /   5 /
C      DATA IMACH( 2) /   6 /
C      DATA IMACH( 3) /   7 /
C      DATA IMACH( 4) /   6 /
C      DATA IMACH( 5) /  32 /
C      DATA IMACH( 6) /   4 /
C      DATA IMACH( 7) /   2 /
C      DATA IMACH( 8) /  31 /
C      DATA IMACH( 9) / Z7FFFFFFF /
C      DATA IMACH(10) /  16 /
C      DATA IMACH(11) /   6 /
C      DATA IMACH(12) / -64 /
C      DATA IMACH(13) /  63 /
C      DATA IMACH(14) /  14 /
C      DATA IMACH(15) / -64 /
C      DATA IMACH(16) /  63 /
C
C     MACHINE CONSTANTS FOR THE INTERDATA 8/32
C     WITH THE UNIX SYSTEM FORTRAN 77 COMPILER.
C
C     FOR THE INTERDATA FORTRAN VII COMPILER REPLACE
C     THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S.
C
C === MACHINE = INTERDATA.8/32.UNIX
C      DATA IMACH( 1) /   5 /
C      DATA IMACH( 2) /   6 /
C      DATA IMACH( 3) /   6 /
C      DATA IMACH( 4) /   6 /
C      DATA IMACH( 5) /  32 /
C      DATA IMACH( 6) /   4 /
C      DATA IMACH( 7) /   2 /
C      DATA IMACH( 8) /  31 /
C      DATA IMACH( 9) / Z'7FFFFFFF' /
C      DATA IMACH(10) /  16 /
C      DATA IMACH(11) /   6 /
C      DATA IMACH(12) / -64 /
C      DATA IMACH(13) /  62 /
C      DATA IMACH(14) /  14 /
C      DATA IMACH(15) / -64 /
C      DATA IMACH(16) /  62 /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR).
C
C === MACHINE = PDP-10.KA
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    7 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   36 /
C      DATA IMACH( 6) /    5 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   35 /
C      DATA IMACH( 9) / "377777777777 /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   27 /
C      DATA IMACH(12) / -128 /
C      DATA IMACH(13) /  127 /
C      DATA IMACH(14) /   54 /
C      DATA IMACH(15) / -101 /
C      DATA IMACH(16) /  127 /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR).
C
C === MACHINE = PDP-10.KI
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    7 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   36 /
C      DATA IMACH( 6) /    5 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   35 /
C      DATA IMACH( 9) / "377777777777 /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   27 /
C      DATA IMACH(12) / -128 /
C      DATA IMACH(13) /  127 /
C      DATA IMACH(14) /   62 /
C      DATA IMACH(15) / -128 /
C      DATA IMACH(16) /  127 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING
C     32-BIT INTEGER ARITHMETIC.
C
C === MACHINE = PDP-11.32-BIT
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    7 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   32 /
C      DATA IMACH( 6) /    4 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   31 /
C      DATA IMACH( 9) / 2147483647 /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   24 /
C      DATA IMACH(12) / -127 /
C      DATA IMACH(13) /  127 /
C      DATA IMACH(14) /   56 /
C      DATA IMACH(15) / -127 /
C      DATA IMACH(16) /  127 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING
C     16-BIT INTEGER ARITHMETIC. 
C
C === MACHINE = PDP-11.16-BIT
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    7 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   16 /
C      DATA IMACH( 6) /    2 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   15 /
C      DATA IMACH( 9) / 32767 /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   24 /
C      DATA IMACH(12) / -127 /
C      DATA IMACH(13) /  127 /
C      DATA IMACH(14) /   56 /
C      DATA IMACH(15) / -127 /
C      DATA IMACH(16) /  127 /
C
C     MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000.
C
C === MACHINE = SEQUENT.BALANCE.8000
C      DATA IMACH( 1) /     0 /
C      DATA IMACH( 2) /     0 /
C      DATA IMACH( 3) /     7 /
C      DATA IMACH( 4) /     0 /
C      DATA IMACH( 5) /    32 /
C      DATA IMACH( 6) /     1 /
C      DATA IMACH( 7) /     2 /
C      DATA IMACH( 8) /    31 /
C      DATA IMACH( 9) /  2147483647 /
C      DATA IMACH(10) /     2 /
C      DATA IMACH(11) /    24 /
C      DATA IMACH(12) /  -125 /
C      DATA IMACH(13) /   128 /
C      DATA IMACH(14) /    53 /
C      DATA IMACH(15) / -1021 /
C      DATA IMACH(16) /  1024 /
C
C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILER
C
C === MACHINE = UNIVAC.1100
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    1 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   36 /
C      DATA IMACH( 6) /    4 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   35 /
C      DATA IMACH( 9) / O377777777777 /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   27 /
C      DATA IMACH(12) / -128 /
C      DATA IMACH(13) /  127 /
C      DATA IMACH(14) /   60 /
C      DATA IMACH(15) /-1024 /
C      DATA IMACH(16) / 1023 /
C
C     MACHINE CONSTANTS FOR THE VAX 11/780
C
C === MACHINE = VAX.11/780
C      DATA IMACH(1) /    5 /
C      DATA IMACH(2) /    6 /
C      DATA IMACH(3) /    5 /
C      DATA IMACH(4) /    6 /
C      DATA IMACH(5) /   32 /
C      DATA IMACH(6) /    4 /
C      DATA IMACH(7) /    2 /
C      DATA IMACH(8) /   31 /
C      DATA IMACH(9) /2147483647 /
C      DATA IMACH(10)/    2 /
C      DATA IMACH(11)/   24 /
C      DATA IMACH(12)/ -127 /
C      DATA IMACH(13)/  127 /
C      DATA IMACH(14)/   56 /
C      DATA IMACH(15)/ -127 /
C      DATA IMACH(16)/  127 /
C
C
C***FIRST EXECUTABLE STATEMENT  I1MACH
C
      I1MACH=IMACH(I)
      RETURN
C
      END
      REAL FUNCTION SDOT(N,SX,INCX,SY,INCY)
C***BEGIN PROLOGUE  SDOT
C***DATE WRITTEN   791001   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***REVISION HISTORY  (YYMMDD)
C   000330  Modified array declarations.  (JEC)
C
C***CATEGORY NO.  D1A4
C***KEYWORDS  BLAS,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  S.P. inner product of s.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       SX  single precision vector with N elements
C     INCX  storage spacing between elements of SX
C       SY  single precision vector with N elements
C     INCY  storage spacing between elements of SY
C
C     --Output--
C     SDOT  single precision dot product (zero if N .LE. 0)
C
C     Returns the dot product of single precision SX and SY.
C     SDOT = sum for I = 0 to N-1 of  SX(LX+I*INCX) * SY(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  SDOT
C
      REAL SX(*),SY(*)
C***FIRST EXECUTABLE STATEMENT  SDOT
      SDOT = 0.0E0
      IF(N.LE.0)RETURN
      IF(INCX.EQ.INCY) IF(INCX-1)5,20,60
    5 CONTINUE
C
C        CODE FOR UNEQUAL INCREMENTS 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
        SDOT = SDOT + SX(IX)*SY(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
        SDOT = SDOT + SX(I)*SY(I)
   30 CONTINUE
      IF( N .LT. 5 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,5
        SDOT = SDOT + SX(I)*SY(I) + SX(I + 1)*SY(I + 1) +
     1   SX(I + 2)*SY(I + 2) + SX(I + 3)*SY(I + 3) + SX(I + 4)*SY(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
        SDOT = SDOT + SX(I)*SY(I)
   70   CONTINUE
      RETURN
      END
      REAL FUNCTION SNRM2(N,SX,INCX)
C***BEGIN PROLOGUE  SNRM2
C***DATE WRITTEN   791001   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***REVISION HISTORY  (YYMMDD)
C   000330  Modified array declarations.  (JEC)
C
C***CATEGORY NO.  D1A3B
C***KEYWORDS  BLAS,EUCLIDEAN,L2,LENGTH,LINEAR ALGEBRA,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 s.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       SX  single precision vector with N elements
C     INCX  storage spacing between elements of SX
C
C     --Output--
C    SNRM2  single precision result (zero if N .LE. 0)
C
C     Euclidean norm of the N-vector stored in SX() 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  SQRT(U/EPS)  over all known machines.
C         CUTHI = minimum of  SQRT(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  SNRM2
      INTEGER          NEXT
      REAL   SX(*),  CUTLO, CUTHI, HITEST, SUM, XMAX, ZERO, ONE
      DATA   ZERO, ONE /0.0E0, 1.0E0/
C
      DATA CUTLO, CUTHI / 4.441E-16,  1.304E19 /
C***FIRST EXECUTABLE STATEMENT  SNRM2
      IF(N .GT. 0) GO TO 10
         SNRM2  = 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( ABS(SX(I)) .GT. CUTLO) GO TO 85
      ASSIGN 50 TO NEXT
      XMAX = ZERO
C
C                        PHASE 1.  SUM IS ZERO
C
   50 IF( SX(I) .EQ. ZERO) GO TO 200
      IF( ABS(SX(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 / SX(I)) / SX(I)
  105 XMAX = ABS(SX(I))
      GO TO 115
C
C                   PHASE 2.  SUM IS SMALL.
C                             SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
C
   70 IF( ABS(SX(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( ABS(SX(I)) .LE. XMAX ) GO TO 115
         SUM = ONE + SUM * (XMAX / SX(I))**2
         XMAX = ABS(SX(I))
         GO TO 200
C
  115 SUM = SUM + (SX(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(ABS(SX(J)) .GE. HITEST) GO TO 100
   95    SUM = SUM + SX(J)**2
      SNRM2 = SQRT( 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
      SNRM2 = XMAX * SQRT(SUM)
  300 CONTINUE
      RETURN
      END
      SUBROUTINE NL2CVP(IV, LIV, LV, P, V)
C
C  ***  PRINT COVARIANCE MATRIX FOR NL2IT  ***
C
      INTEGER LIV, LV, P
      INTEGER IV(LIV)
      REAL V(LV)
C
C  ***  INTRINSIC FUNCTIONS  ***
C/+
      INTEGER IABS, MOD
C/
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER COV1, I, II, I1, J, PU
C
C     ***  IV SUBSCRIPTS  ***
C
      INTEGER COVMAT, COVPRT, COVREQ, NEEDHD, OUTLEV, PRUNIT
C
C/6
C     DATA COVMAT/26/, COVPRT/14/, COVREQ/15/, NEEDHD/36/, OUTLEV/19/,
C    1     PRUNIT/21/
C/7
      PARAMETER (COVMAT=26, COVPRT=14, COVREQ=15, NEEDHD=36, OUTLEV=19,
     1           PRUNIT=21)
C/
      IF (MOD(IV(COVPRT),2) .EQ. 0) GO TO 999
      IF (IV(1) .GT. 8) GO TO 999
      PU = IV(PRUNIT)
      IF (PU .EQ. 0) GO TO 999
      COV1 = IV(COVMAT)
      IV(NEEDHD) = 1
      IF (COV1) 10,40,60
 10   IF (-1 .EQ. COV1) WRITE(PU,20)
 20   FORMAT(/43H ++++++ INDEFINITE COVARIANCE MATRIX ++++++)
      IF (-2 .EQ. COV1) WRITE(PU,30)
 30   FORMAT(/52H ++++++ OVERSIZE STEPS IN COMPUTING COVARIANCE +++++)
      GO TO 999
C
 40   WRITE(PU,50)
 50   FORMAT(/45H ++++++ COVARIANCE MATRIX NOT COMPUTED ++++++)
      GO TO 999
C
 60   I = IABS(IV(COVREQ))
      IF (I .LE. 1) WRITE(PU,70)
 70   FORMAT(/48H COVARIANCE = SCALE * H**-1 * (J**T * J) * H**-1/
     1       23H WHERE H = F.D. HESSIAN/)
      IF (I .EQ. 2) WRITE(PU,80)
 80   FORMAT(/57H COVARIANCE = H**-1, WHERE H = FINITE-DIFFERENCE HESSIA
     1N/)
      IF (I .GT. 2) WRITE(PU,90)
 90   FORMAT(/30H COVARIANCE = SCALE * J**T * J/)
      II = COV1 - 1
      IF (IV(OUTLEV) .LE. 0) GO TO 120
      DO 100 I = 1, P
         I1 = II + 1
         II = II + I
         WRITE(PU,110) I, (V(J), J = I1, II)
 100     CONTINUE
 110  FORMAT(4H ROW,I3,2X,9E12.4/(9X,9E12.4))
      GO TO 999
C
 120  DO 130 I = 1, P
         I1 = II + 1
         II = II + I
         WRITE(PU,140) I, (V(J), J = I1, II)
 130     CONTINUE
 140  FORMAT(4H ROW,I3,2X,5E12.4/(9X,5E12.4))
C
 999  RETURN
C  ***  LAST CARD OF MLPCVP FOLLOWS  ***
      END
      SUBROUTINE GLRIT(D, G, IV, LIV, LV, P, PS, V, X, Y)
C
C  ***  CARRY OUT NL2SOL-LIKE ITERATIONS FOR GENERALIZED LINEAR   ***
C  ***  REGRESSION PROBLEMS (AND OTHERS OF SIMILAR STRUCTURE)     ***
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER LIV, LV, P, PS
      INTEGER IV(LIV)
      REAL D(P), G(P), V(LV), X(P), Y(P)
C
C--------------------------  PARAMETER USAGE  --------------------------
C
C D.... SCALE VECTOR.
C IV... INTEGER VALUE ARRAY.
C LIV.. LENGTH OF IV.  MUST BE AT LEAST 80.
C LH... LENGTH OF H = P*(P+1)/2.
C LV... LENGTH OF V.  MUST BE AT LEAST P*(3*P + 19)/2 + 7.
C G.... GRADIENT AT X (WHEN IV(1) = 2).
C HC... GAUSS-NEWTON HESSIAN AT X (WHEN IV(1) = 2).
C P.... NUMBER OF PARAMETERS (COMPONENTS IN X).
C PS... NUMBER OF NONZERO ROWS AND COLUMNS IN S.
C V.... FLOATING-POINT VALUE ARRAY.
C X.... PARAMETER VECTOR.
C Y.... PART OF YIELD VECTOR (WHEN IV(1)= 2, SCRATCH OTHERWISE).
C
C  ***  DISCUSSION  ***
C
C        GLRIT PERFORMS NL2SOL-LIKE ITERATIONS FOR A VARIETY OF
C     REGRESSION PROBLEMS THAT ARE SIMILAR TO NONLINEAR LEAST-SQUARES
C     IN THAT THE HESSIAN IS THE SUM OF TWO TERMS, A READILY-COMPUTED
C     FIRST-ORDER TERM AND A SECOND-ORDER TERM.  THE CALLER SUPPLIES
C     THE FIRST-ORDER TERM OF THE HESSIAN IN HC (LOWER TRIANGLE, STORED
C     COMPACTLY BY ROWS), AND GLRIT BUILDS AN APPROXIMATION, S, TO THE
C     SECOND-ORDER TERM.  THE CALLER ALSO PROVIDES THE FUNCTION VALUE,
C     GRADIENT, AND PART OF THE YIELD VECTOR USED IN UPDATING S.  GLRIT
C     DECIDES DYNAMICALLY WHETHER OR NOT TO USE S WHEN CHOOSING THE
C     NEXT STEP TO TRY...  THE HESSIAN APPROXIMATION USED IS EITHER HC
C     ALONE (GAUSS-NEWTON MODEL) OR HC + S (AUGMENTED MODEL).
C     IF PS .LT. P, THEN ROWS AND COLUMNS PS+1...P OF S ARE KEPT
C     CONSTANT.  THEY WILL BE ZERO UNLESS THE CALLER SETS IV(INITS) TO
C     1 OR 2 AND SUPPLIES NONZERO VALUES FOR THEM.
C
C        FOR UPDATING S, GLRIT ASSUMES THAT THE GRADIENT HAS THE FORM
C     OF A SUM OVER I OF RHO(I,X)*GRAD(R(I,X)), WHERE GRAD DENOTES THE
C     GRADIENT WITH RESPECT TO X.  THE TRUE SECOND-ORDER TERM THEN IS
C     THE SUM OVER I OF RHO(I,X)*HESSIAN(R(I,X)).  IF X = X0 + STEP,
C     THEN WE WISH TO UPDATE S SO THAT S*STEP IS THE SUM OVER I OF
C     RHO(I,X)*(GRAD(R(I,X)) - GRAD(R(I,X0))).  THE CALLER MUST SUPPLY
C     PART OF THIS IN Y, NAMELY THE SUM OVER I OF
C     RHO(I,X)*GRAD(R(I,X0)), WHEN CALLING GLRIT WITH IV(1) = 2 AND
C     IV(MODE) = 0 (WHERE MODE = 38).  G THEN CONTANS THE OTHER PART,
C     SO THAT THE DESIRED YIELD VECTOR IS G - Y.  IF PS .LT. P, THEN
C     THE ABOVE DISCUSSION APPLIES ONLY TO THE FIRST PS COMPONENTS OF
C     GRAD(R(I,X)), STEP, AND Y.
C
C        PARAMETERS IV, P, V, AND X ARE THE SAME AS THE CORRESPONDING
C     ONES TO NL2SOL (WHICH SEE), EXCEPT THAT V CAN BE SHORTER
C     (SINCE THE PART OF V THAT NL2SOL USES FOR STORING D, J, AND R IS
C     NOT NEEDED).  MOREOVER, COMPARED WITH NL2SOL, IV(1) MAY HAVE THE
C     TWO ADDITIONAL OUTPUT VALUES 1 AND 2, WHICH ARE EXPLAINED BELOW,
C     AS IS THE USE OF IV(TOOBIG) AND IV(NFGCAL).  THE VALUES IV(D),
C     IV(J), AND IV(R), WHICH ARE OUTPUT VALUES FROM NL2SOL (AND
C     NL2SNO), ARE NOT REFERENCED BY GLRIT OR THE SUBROUTINES IT CALLS.
C
C        WHEN GLRIT IS FIRST CALLED, I.E., WHEN GLRIT IS CALLED WITH
C     IV(1) = 0 OR 12, V(F), G, AND HC NEED NOT BE INITIALIZED.  TO
C     OBTAIN THESE STARTING VALUES, GLRIT RETURNS FIRST WITH IV(1) = 1,
C     THEN WITH IV(1) = 2, WITH IV(MODE) = -1 IN BOTH CASES.  ON
C     SUBSEQUENT RETURNS WITH IV(1) = 2, IV(MODE) = 0 IMPLIES THAT
C     Y MUST ALSO BE SUPPLIED.  (NOTE THAT Y IS USED FOR SCRATCH -- ITS
C     INPUT CONTENTS ARE LOST.  BY CONTRAST, HC IS NEVER CHANGED.)
C     ONCE CONVERGENCE HAS BEEN OBTAINED, IV(COVREQ) AND IV(COVPRT) MAY
C     IMPLY THAT A FINITE-DIFFERENCE HESSIAN SHOULD BE COMPUTED FOR USE
C     IN COMPUTING A COVARIANCE MATRIX.  IN THIS CASE GLRIT WILL MAKE A
C     NUMBER OF RETURNS WITH IV(1) = 1 OR 2 AND IV(MODE) POSITIVE.
C     WHEN IV(MODE) IS POSITIVE, Y SHOULD NOT BE CHANGED.
C
C IV(1) = 1 MEANS THE CALLER SHOULD SET V(F) (I.E., V(10)) TO F(X), THE
C             FUNCTION VALUE AT X, AND CALL GLRIT AGAIN, HAVING CHANGED
C             NONE OF THE OTHER PARAMETERS.  AN EXCEPTION OCCURS IF F(X)
C             CANNOT BE EVALUATED (E.G. IF OVERFLOW WOULD OCCUR), WHICH
C             MAY HAPPEN BECAUSE OF AN OVERSIZED STEP.  IN THIS CASE
C             THE CALLER SHOULD SET IV(TOOBIG) = IV(2) TO 1, WHICH WILL
C             CAUSE GLRIT TO IGNORE V(F) AND TRY A SMALLER STEP.  NOTE
C             THAT THE CURRENT FUNCTION EVALUATION COUNT IS AVAILABLE
C             IN IV(NFCALL) = IV(6).  THIS MAY BE USED TO IDENTIFY
C             WHICH COPY OF SAVED INFORMATION SHOULD BE USED IN COM-
C             PUTING G, HC, AND Y THE NEXT TIME GLRIT RETURNS WITH
C             IV(1) = 2.  SEE MLPIT FOR AN EXAMPLE OF THIS.
C IV(1) = 2 MEANS THE CALLER SHOULD SET G TO G(X), THE GRADIENT OF F AT
C             X.  THE CALLER SHOULD ALSO SET HC TO THE GAUSS-NEWTON
C             HESSIAN AT X.  IF IV(MODE) = 0, THEN THE CALLER SHOULD
C             ALSO COMPUTE THE PART OF THE YIELD VECTOR DESCRIBED ABOVE.
C             THE CALLER SHOULD THEN CALL GLRIT AGAIN (WITH IV(1) = 2).
C             THE CALLER MAY ALSO CHANGE D AT THIS TIME, BUT SHOULD NOT
C             CHANGE X.  NOTE THAT IV(NFGCAL) = IV(7) CONTAINS THE
C             VALUE THAT IV(NFCALL) HAD DURING THE RETURN WITH
C             IV(1) = 1 IN WHICH X HAD THE SAME VALUE AS IT NOW HAS.
C             IV(NFGCAL) IS EITHER IV(NFCALL) OR IV(NFCALL) - 1.  MLPIT
C             IS AN EXAMPLE WHERE THIS INFORMATION IS USED.  IF G OR HC
C             CANNOT BE EVALUATED AT X, THEN THE CALLER MAY SET
C             IV(NFGCAL) TO 0, IN WHICH CASE GLRIT WILL RETURN WITH
C             IV(1) = 15.
C
C  ***  GENERAL  ***
C
C     CODED BY DAVID M. GAY.
C     THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH
C     SUPPORTED IN PART BY D.O.E. GRANT EX-76-A-01-2295 TO MIT/CCREMS.
C
C        (SEE NL2SOL FOR REFERENCES.)
C
C+++++++++++++++++++++++++++  DECLARATIONS  ++++++++++++++++++++++++++++
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER DUMMY, DIG1, G01, H1, HC1, I, IPIV1, J, K, L, LMAT1,
     1        LSTGST, PP1O2, QTR1, RMAT1, STEP1, STPMOD, S1, TEMP1,
     2        TEMP2, W1, X01
      REAL E, STTSST, T, T1
C
C     ***  CONSTANTS  ***
C
      REAL HALF, NEGONE, ONE, ZERO
C
C  ***  INTRINSIC FUNCTIONS  ***
C/+
      INTEGER IABS
      REAL  ABS
C/
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      EXTERNAL ASSST,  FDHES, GQTST, ITSUM, LMSTEP, LSQUAR,
     1         LTVMUL, LVMUL, PARCK, SLUPDT, SLVMUL, STOPX, VAXPY,
     2          VSCOPY
      LOGICAL STOPX
C
C ASSST.... ASSESSES CANDIDATE STEP.
C FDHES.... COMPUTE FINITE-DIFFERENCE HESSIAN (FOR COVARIANCE).
C GQTST.... COMPUTES GOLDFELD-QUANDT-TROTTER STEP (AUGMENTED MODEL).
C ITSUM.... PRINTS ITERATION SUMMARY AND INFO ON INITIAL AND FINAL X.
C LMSTEP... COMPUTES LEVENBERG-MARQUARDT STEP (GAUSS-NEWTON MODEL).
C LSQUAR... COMPUTES L * L**T FROM LOWER TRIANGULAR MATRIX L.
C LTVMUL... COMPUTES L**T * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX.
C LVMUL.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX.
C PARCK.... CHECK VALIDITY OF IV AND V INPUT COMPONENTS.
C SLUPDT... PERFORMS QUASI-NEWTON UPDATE ON COMPACTLY STORED LOWER TRI-
C             ANGLE OF A SYMMETRIC MATRIX.
C STOPX.... RETURNS .TRUE. IF THE BREAK KEY HAS BEEN PRESSED.
C VAXPY.... COMPUTES SCALAR TIMES ONE VECTOR PLUS ANOTHER.
C VSCOPY... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR.
C
C  ***  SUBSCRIPTS FOR IV AND V  ***
C
      INTEGER CNVCOD, COSMIN, COVMAT, COVPRT, COVREQ, DGNORM, DIG,
     1        DSTNRM, F, FDH, FDIF, FUZZ, F0, GTSTEP, H, HC, IERR,
     2        INCFAC, INITS, IPIVOT, IRC, KAGQT, KALM, LMAT, LMAX0,
     3        LMAXS, MODE, MODEL, MXFCAL, MXITER, NEXTV, NFCALL, NFGCAL,
     4        NFCOV, NGCOV, NGCALL, NITER, NVSAVE, PHMXFC, PREDUC, QTR,
     5        RADFAC, RADINC, RADIUS, RAD0, RDREQ, REGD, RMAT, S, SIZE,
     6        STEP, STGLIM, STLSTG, STPPAR, SUSED, SWITCH, TOOBIG,
     7        TUNER4, TUNER5, VNEED, VSAVE, W, WSCALE, XIRC, X0
C
C  ***  IV SUBSCRIPT VALUES  ***
C
C/6
C     DATA CNVCOD/55/, COVMAT/26/, COVPRT/14/, COVREQ/15/, DIG/37/,
C    1     FDH/74/, H/56/, HC/71/, IERR/75/, INITS/25/, IPIVOT/76/,
C    2     IRC/29/, KAGQT/33/, KALM/34/, LMAT/42/, MODE/35/, MODEL/5/,
C    3     MXFCAL/17/, MXITER/18/, NEXTV/47/, NFCALL/6/, NFGCAL/7/,
C    4     NFCOV/52/, NGCOV/53/, NGCALL/30/, NITER/31/, QTR/77/,
C    5     RADINC/8/, RDREQ/57/, REGD/67/, RMAT/78/, S/62/, STEP/40/,
C    6     STGLIM/11/, STLSTG/41/, SUSED/64/, SWITCH/12/, TOOBIG/2/,
C    7     VNEED/4/, VSAVE/60/, W/65/, XIRC/13/, X0/43/
C/7
      PARAMETER (CNVCOD=55, COVMAT=26, COVPRT=14, COVREQ=15, DIG=37,
     1           FDH=74, H=56, HC=71, IERR=75, INITS=25, IPIVOT=76,
     2           IRC=29, KAGQT=33, KALM=34, LMAT=42, MODE=35, MODEL=5,
     3           MXFCAL=17, MXITER=18, NEXTV=47, NFCALL=6, NFGCAL=7,
     4           NFCOV=52, NGCOV=53, NGCALL=30, NITER=31, QTR=77,
     5           RADINC=8, RDREQ=57, REGD=67, RMAT=78, S=62, STEP=40,
     6           STGLIM=11, STLSTG=41, SUSED=64, SWITCH=12, TOOBIG=2,
     7           VNEED=4, VSAVE=60, W=65, XIRC=13, X0=43)
C/
C
C  ***  V SUBSCRIPT VALUES  ***
C
C/6
C     DATA COSMIN/47/, DGNORM/1/, DSTNRM/2/, F/10/, FDIF/11/, FUZZ/45/,
C    1     F0/13/, GTSTEP/4/, INCFAC/23/, LMAX0/35/, LMAXS/36/,
C    2     NVSAVE/9/, PHMXFC/21/, PREDUC/7/, RADFAC/16/, RADIUS/8/,
C    3     RAD0/9/, SIZE/55/, STPPAR/5/, TUNER4/29/, TUNER5/30/,
C    4     WSCALE/56/
C/7
      PARAMETER (COSMIN=47, DGNORM=1, DSTNRM=2, F=10, FDIF=11, FUZZ=45,
     1           F0=13, GTSTEP=4, INCFAC=23, LMAX0=35, LMAXS=36,
     2           NVSAVE=9, PHMXFC=21, PREDUC=7, RADFAC=16, RADIUS=8,
     3           RAD0=9, SIZE=55, STPPAR=5, TUNER4=29, TUNER5=30,
     4           WSCALE=56)
C/
C
C
C/6
C     DATA HALF/0.5E+0/, NEGONE/-1.E+0/, ONE/1.E+0/, ZERO/0.E+0/
C/7
      PARAMETER (HALF=0.5E+0, NEGONE=-1.E+0, ONE=1.E+0, ZERO=0.E+0)
C/
C
C+++++++++++++++++++++++++++++++  BODY  ++++++++++++++++++++++++++++++++
C
      I = IV(1)
      IF (I .EQ. 1) GO TO 40
      IF (I .EQ. 2) GO TO 50
C
      IV(VNEED) = IV(VNEED) + P*(3*P + 19)/2 + 7
      CALL PARCK(1, D, IV, LIV, LV, P, V)
      I = IV(1) - 2
      IF (I .GT. 12) GO TO 999
      GO TO (260, 260, 260, 260, 260, 260, 140, 90, 140, 10, 10, 20), I
C
C  ***  STORAGE ALLOCATION  ***
C
 10   PP1O2 = P * (P + 1) / 2
      IV(S) = IV(LMAT) + PP1O2
      IV(X0) = IV(S) + PP1O2
      IV(STEP) = IV(X0) + P
      IV(STLSTG) = IV(STEP) + P
      IV(DIG) = IV(STLSTG) + P
      IV(W) = IV(DIG) + P
      IV(H) = IV(W) + 4*P + 7
      IV(NEXTV) = IV(H) + PP1O2
      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(STGLIM) = 2
      IV(TOOBIG) = 0
      IV(CNVCOD) = 0
      IV(COVMAT) = 0
      IV(NFCOV) = 0
      IV(NGCOV) = 0
      IV(RADINC) = 0
      V(RAD0) = ZERO
      V(STPPAR) = ZERO
      V(RADIUS) = V(LMAX0) / (ONE + V(PHMXFC))
C
C  ***  SET INITIAL MODEL AND S MATRIX  ***
C
      IV(MODEL) = 1
      IF (IV(S) .LT. 0) GO TO 999
      IF (IV(INITS) .EQ. 2) IV(MODEL) = 2
      S1 = IV(S)
      IF (IV(INITS) .EQ. 0) CALL VSCOPY(PP1O2, V(S1), ZERO)
      IV(1) = 1
      J = IV(IPIVOT)
      IF (J .LE. 0) GO TO 999
      DO 30 I = 1, P
         IV(J) = I
         J = J + 1
 30      CONTINUE
      GO TO 999
C
C  ***  NEW FUNCTION VALUE  ***
C
 40   IF (IV(MODE) .EQ. 0) GO TO 260
      IF (IV(MODE) .GT. 0) GO TO 440
C
      IV(1) = 2
      IF (IV(TOOBIG) .EQ. 0) GO TO 999
         IV(1) = 63
         GO TO 999
C
C  ***  MAKE SURE GRADIENT COULD BE COMPUTED  ***
C
 50   IF (IV(NFGCAL) .NE. 0) GO TO 60
         IV(1) = 65
         GO TO 999
C
C  ***  NEW GRADIENT  ***
C
 60   IV(KALM) = -1
      IV(KAGQT) = -1
      IF (IV(MODE) .GT. 0) GO TO 440
      IF (IV(HC) .LE. 0 .AND. IV(RMAT) .LE. 0) GO TO 500
C
C  ***  COMPUTE  D**-1 * GRADIENT  ***
C
      DIG1 = IV(DIG)
      K = DIG1
      DO 70 I = 1, P
         V(K) = G(I) / D(I)
         K = K + 1
 70      CONTINUE
      V(DGNORM) = SNRM2(P,V(DIG1),1)
C
      IF (IV(CNVCOD) .NE. 0) GO TO 430
      IF (IV(MODE) .EQ. 0) GO TO 360
      IV(1) = 2
      IV(MODE) = 0
C
C
C-----------------------------  MAIN LOOP  -----------------------------
C
C
C  ***  PRINT ITERATION SUMMARY, CHECK ITERATION LIMIT  ***
C
 80   CALL ITSUM(D, G, IV, LIV, LV, P, V, X)
 90   K = IV(NITER)
      IF (K .LT. IV(MXITER)) GO TO 100
         IV(1) = 10
         GO TO 999
 100  IV(NITER) = K + 1
C
C  ***  UPDATE RADIUS  ***
C
      IF (K .EQ. 0) GO TO 120
      STEP1 = IV(STEP)
      DO 110 I = 1, P
         V(STEP1) = D(I) * V(STEP1)
         STEP1 = STEP1 + 1
 110     CONTINUE
      STEP1 = IV(STEP)
      V(RADIUS) = V(RADFAC) * SNRM2(P,V(STEP1),1)
C
C  ***  INITIALIZE FOR START OF NEXT ITERATION  ***
C
 120  X01 = IV(X0)
      V(F0) = V(F)
      IV(IRC) = 4
      IV(H) = -IABS(IV(H))
      IV(SUSED) = IV(MODEL)
C
C     ***  COPY X TO X0  ***
C
      CALL SCOPY(P,X,1,V(X01),1)
C
C  ***  CHECK STOPX AND FUNCTION EVALUATION LIMIT  ***
C
 130  IF (.NOT. STOPX(DUMMY)) GO TO 150
         IV(1) = 11
         GO TO 160
C
C     ***  COME HERE WHEN RESTARTING AFTER FUNC. EVAL. LIMIT OR STOPX.
C
 140  IF (V(F) .GE. V(F0)) GO TO 150
         V(RADFAC) = ONE
         K = IV(NITER)
         GO TO 100
C
 150  IF (IV(NFCALL) .LT. IV(MXFCAL) + IV(NFCOV)) GO TO 170
         IV(1) = 9
 160     IF (V(F) .GE. V(F0)) GO TO 999
C
C        ***  IN CASE OF STOPX OR FUNCTION EVALUATION LIMIT WITH
C        ***  IMPROVED V(F), EVALUATE THE GRADIENT AT X.
C
              IV(CNVCOD) = IV(1)
              GO TO 350
C
C. . . . . . . . . . . . .  COMPUTE CANDIDATE STEP  . . . . . . . . . .
C
 170  STEP1 = IV(STEP)
      W1 = IV(W)
      H1 = IV(H)
      T1 = ONE
      IF (IV(MODEL) .EQ. 2) GO TO 200
         T1 = ZERO
C
C        ***  COMPUTE LEVENBERG-MARQUARDT STEP IF POSSIBLE...
C
         RMAT1 = IV(RMAT)
         IF (RMAT1 .LE. 0) GO TO 200
         QTR1 = IV(QTR)
         IF (QTR1 .LE. 0) GO TO 200
         IPIV1 = IV(IPIVOT)
         CALL LMSTEP(D, G, IV(IERR), IV(IPIV1), IV(KALM), P, V(QTR1),
     1               V(RMAT1),V(STEP1), V, V(W1))
C        *** H IS STORED IN THE END OF W AND HAS JUST BEEN OVERWRITTEN,
C        *** SO WE MARK IT INVALID...
         IV(H) = -IABS(H1)
C        *** EVEN IF H WERE STORED ELSEWHERE, IT WOULD BE NECESSARY TO
C        *** MARK INVALID THE INFORMATION GQTST MAY HAVE STORED IN V...
         IV(KAGQT) = -1
         GO TO 250
C
 200  IF (H1 .GT. 0) GO TO 240
C
C     ***  SET H TO  D**-1 * (HC + T1*S) * D**-1.  ***
C
         H1 = -H1
         IV(H) = H1
         IV(FDH) = 0
         J = IV(HC)
         IF (J .GT. 0) GO TO 210
            J = H1
            RMAT1 = IV(RMAT)
            CALL LSQUAR(P, V(H1), V(RMAT1))
 210     S1 = IV(S)
         DO 230 I = 1, P
              T = ONE / D(I)
              DO 220 K = 1, I
                   V(H1) = T * (V(J) + T1*V(S1)) / D(K)
                   J = J + 1
                   H1 = H1 + 1
                   S1 = S1 + 1
 220               CONTINUE
 230          CONTINUE
         H1 = IV(H)
         IV(KAGQT) = -1
C
C  ***  COMPUTE ACTUAL GOLDFELD-QUANDT-TROTTER STEP  ***
C
 240  DIG1 = IV(DIG)
      LMAT1 = IV(LMAT)
      CALL GQTST(D, V(DIG1), V(H1), IV(KAGQT), V(LMAT1), P, V(STEP1),
     1            V, V(W1))
      IF (IV(KALM) .GT. 0) IV(KALM) = 0
C
C
C  ***  COMPUTE F(X0 + STEP)  ***
C
 250  IF (IV(IRC) .EQ. 6) GO TO 260
      X01 = IV(X0)
      STEP1 = IV(STEP)
      CALL VAXPY(P, X, ONE, V(STEP1), V(X01))
      IV(NFCALL) = IV(NFCALL) + 1
      IV(1) = 1
      IV(TOOBIG) = 0
      GO TO 999
C
C. . . . . . . . . . . . .  ASSESS CANDIDATE STEP  . . . . . . . . . . .
C
 260  STEP1 = IV(STEP)
      LSTGST = IV(STLSTG)
      X01 = IV(X0)
      CALL ASSST(D, IV, P, V(STEP1), V(LSTGST), V, X, V(X01))
C
C  ***  IF NECESSARY, SWITCH MODELS AND/OR RESTORE R  ***
C
      IF (IV(SWITCH) .EQ. 0) GO TO 270
         IV(H) = -IABS(IV(H))
         IV(SUSED) = IV(SUSED) + 2
         L = IV(VSAVE)
         CALL SCOPY(NVSAVE,V(L),1,V,1)
 270  L = IV(IRC) - 4
      STPMOD = IV(MODEL)
      IF (L .GT. 0) GO TO (290,300,310,310,310,310,310,310,420,360), L
C
C  ***  DECIDE WHETHER TO CHANGE MODELS  ***
C
      E = V(PREDUC) - V(FDIF)
      S1 = IV(S)
      CALL SLVMUL(P, Y, V(S1), V(STEP1))
      STTSST = HALF * SDOT(P, V(STEP1),1,Y,1)
      IF (IV(MODEL) .EQ. 1) STTSST = -STTSST
      IF ( ABS(E + STTSST) * V(FUZZ) .GE.  ABS(E)) GO TO 280
C
C     ***  SWITCH MODELS  ***
C
         IV(MODEL) = 3 - IV(MODEL)
         IF (-2 .LT. L) GO TO 320
              IV(H) = -IABS(IV(H))
              IV(SUSED) = IV(SUSED) + 2
              L = IV(VSAVE)
              CALL SCOPY(NVSAVE,V,1,V(L),1)
              GO TO 130
C
 280  IF (-3 .LT. L) GO TO 320
C
C     ***  RECOMPUTE STEP WITH DECREASED RADIUS  ***
C
         V(RADIUS) = V(RADFAC) * V(DSTNRM)
         GO TO 130
C
C  ***  RECOMPUTE STEP, SAVING V VALUES AND R IF NECESSARY  ***
C
 290  V(RADIUS) = V(RADFAC) * V(DSTNRM)
      GO TO 130
C
C  ***  COMPUTE STEP OF LENGTH V(LMAXS) FOR SINGULAR CONVERGENCE TEST
C
 300  V(RADIUS) = V(LMAXS)
      GO TO 170
C
C  ***  CONVERGENCE OR FALSE CONVERGENCE  ***
C
 310  IV(CNVCOD) = L
      IF (V(F) .GE. V(F0)) GO TO 430
         IF (IV(XIRC) .EQ. 14) GO TO 430
              IV(XIRC) = 14
C
C. . . . . . . . . . . .  PROCESS ACCEPTABLE STEP  . . . . . . . . . . .
C
 320  IV(COVMAT) = 0
      IV(REGD) = 0
C
C  ***  SEE WHETHER TO SET V(RADFAC) BY GRADIENT TESTS  ***
C
      IF (IV(IRC) .NE. 3) GO TO 350
         STEP1 = IV(STEP)
         TEMP1 = IV(STLSTG)
         TEMP2 = IV(X0)
C
C     ***  SET  TEMP1 = HESSIAN * STEP  FOR USE IN GRADIENT TESTS  ***
C
         HC1 = IV(HC)
         IF (HC1 .LE. 0) GO TO 330
              CALL SLVMUL(P, V(TEMP1), V(HC1), V(STEP1))
              GO TO 340
 330     RMAT1 = IV(RMAT)
         CALL LTVMUL(P, V(TEMP1), V(RMAT1), V(STEP1))
         CALL LVMUL(P, V(TEMP1), V(RMAT1), V(TEMP1))
C
 340     IF (STPMOD .EQ. 1) GO TO 350
              S1 = IV(S)
              CALL SLVMUL(P, V(TEMP2), V(S1), V(STEP1))
              CALL VAXPY(P, V(TEMP1), ONE, V(TEMP2), V(TEMP1))
C
C  ***  SAVE OLD GRADIENT AND COMPUTE NEW ONE  ***
C
 350  IV(NGCALL) = IV(NGCALL) + 1
      G01 = IV(W)
      CALL SCOPY(P,G,1,V(G01),1)
      IV(1) = 2
      GO TO 999
C
C  ***  INITIALIZATIONS -- G0 = G - G0, ETC.  ***
C
 360  G01 = IV(W)
      CALL VAXPY(P, V(G01), NEGONE, V(G01), G)
      STEP1 = IV(STEP)
      TEMP1 = IV(STLSTG)
      TEMP2 = IV(X0)
      IF (IV(IRC) .NE. 3) GO TO 390
C
C  ***  SET V(RADFAC) BY GRADIENT TESTS  ***
C
C     ***  SET  TEMP1 = D**-1 * (HESSIAN * STEP  +  (G(X0) - G(X)))  ***
C
         K = TEMP1
         L = G01
         DO 370 I = 1, P
              V(K) = (V(K) - V(L)) / D(I)
              K = K + 1
              L = L + 1
 370          CONTINUE
C
C        ***  DO GRADIENT TESTS  ***
C
         IF (SNRM2(P,V(TEMP1),1) .LE. V(DGNORM) * V(TUNER4))  GO TO 380
              IF (SDOT(P, G,1, V(STEP1),1)
     1                  .GE. V(GTSTEP) * V(TUNER5))  GO TO 390
 380               V(RADFAC) = V(INCFAC)
C
C  ***  COMPUTE Y VECTOR NEEDED FOR UPDATING S  ***
C
 390  CALL VAXPY(P, Y, NEGONE, Y, G)
C
C  ***  DETERMINE SIZING FACTOR V(SIZE)  ***
C
C     ***  SET TEMP1 = S * STEP  ***
      S1 = IV(S)
      CALL SLVMUL(P, V(TEMP1), V(S1), V(STEP1))
C
      T1 =  ABS(SDOT(P, V(STEP1),1,V(TEMP1),1))
      T =  ABS(SDOT(P, V(STEP1),1, Y,1))
      V(SIZE) = ONE
      IF (T .LT. T1) V(SIZE) = T / T1
C
C  ***  SET G0 TO WCHMTD CHOICE OF FLETCHER AND AL-BAALI  ***
C
      HC1 = IV(HC)
      IF (HC1 .LE. 0) GO TO 400
         CALL SLVMUL(P, V(G01), V(HC1), V(STEP1))
         GO TO 410
C
 400  RMAT1 = IV(RMAT)
      CALL LTVMUL(P, V(G01), V(RMAT1), V(STEP1))
      CALL LVMUL(P, V(G01), V(RMAT1), V(G01))
C
 410  CALL VAXPY(P, V(G01), ONE, Y, V(G01))
C
C  ***  UPDATE S  ***
C
      CALL SLUPDT(V(S1), V(COSMIN), PS, V(SIZE), V(STEP1), V(TEMP1),
     1            V(TEMP2), V(G01), V(WSCALE), Y)
      IV(1) = 2
      GO TO 80
C
C. . . . . . . . . . . . . .  MISC. DETAILS  . . . . . . . . . . . . . .
C
C  ***  BAD PARAMETERS TO ASSESS  ***
C
 420  IV(1) = 64
      GO TO 999
C
C
C  ***  CONVERGENCE OBTAINED -- SEE WHETHER TO COMPUTE COVARIANCE  ***
C
 430  IF (IV(COVREQ) .EQ. 0 .AND. IV(COVPRT) .EQ. 0 .AND.
     1                            IV(RDREQ) .LE. 0) GO TO 490
      IF (IV(FDH) .NE. 0) GO TO 490
      IF (IV(CNVCOD) .GE. 7) GO TO 490
      IF (IABS(IV(COVREQ)) .GE. 3) GO TO 470
C
C  ***  COMPUTE FINITE-DIFFERENCE HESSIAN FOR COMPUTING COVARIANCE  ***
C
 440  CALL FDHES(D, G, I, IV, LIV, LV, P, V, X)
      GO TO (450, 460, 490), I
 450  IV(NFCOV) = IV(NFCOV) + 1
      IV(NFCALL) = IV(NFCALL) + 1
      IV(1) = 1
      GO TO 999
C
 460  IV(NGCOV) = IV(NGCOV) + 1
      IV(NGCALL) = IV(NGCALL) + 1
      IV(NFGCAL) = IV(NFCALL) + IV(NGCOV)
      IV(1) = 2
      GO TO 999
C
 470  H1 = IABS(IV(H))
      IV(FDH) = H1
      IV(H) = -H1
      HC1 = IV(HC)
      IF (HC1 .LE. 0) GO TO 480
           CALL SCOPY(P*(P+1)/2,V(HC1),1,V(H1),1)
           GO TO 490
 480  RMAT1 = IV(RMAT)
      CALL LSQUAR(P, V(H1), V(RMAT1))
C
 490  IV(MODE) = 0
      IV(1) = IV(CNVCOD)
      IV(CNVCOD) = 0
      GO TO 999
C
C  ***  SPECIAL RETURN FOR MISSING HESSIAN INFORMATION -- BOTH
C  ***  IV(HC) .LE. 0 AND IV(RMAT) .LE. 0
C
 500  IV(1) = 1400
C
 999  RETURN
C
C  ***  LAST CARD OF GLRIT FOLLOWS  ***
      END
      SUBROUTINE QRADR(P, QTR, RMAT, W, Y)
C
C  ***  ADD ROW W TO QR FACTORIZATION WITH R MATRIX RMAT AND
C  ***  Q**T * RESIDUAL = QTR.  Y = NEW COMPONENT OF RESIDUAL
C  ***  CORRESPONDING TO W.
C
C***REVISION HISTORY  (YYMMDD)
C   000330  Modified array declarations.  (JEC)
C
      INTEGER P
      REAL QTR(P), RMAT(*), W(P), Y
C     DIMENSION RMAT(P*(P+1)/2)
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER I, J, II, IJ, IP1
      REAL ONE,  ABS, T, V, ZERO, U1
      REAL U2, RI, WI, RW, WR,  SQRT
      DATA ONE/1.E+0/
      DATA ZERO/0.E+0/
C
C------------------------------ BODY -----------------------------------
C
      II = 0
      DO 60 I = 1, P
         II = II+I
         WI = W(I)
         IF (WI .EQ. ZERO) GOTO  60
         RI = RMAT(II)
         IF (RI .NE. ZERO) GOTO 20
            IJ = II
C           *** SWAP W AND ROW I OF RMAT ***
            DO 10 J = I, P
               T = RMAT(IJ)
               RMAT(IJ) = W(J)
               W(J) = T
               IJ = IJ+J
 10            CONTINUE
            T = QTR(I)
            QTR(I) = Y
            Y = T
            GO TO 60
 20      IP1 = I+1
         IJ = II+I
         IF ( ABS(WI) .LE.  ABS(RI)) GO TO 40
            RW = RI/WI
            T =  SQRT(ONE+RW**2)
            IF (RW .GT. ZERO) T = -T
            V = RW-T
            U1 = ONE/T
            U2 = ONE/(T*V)
            RMAT(II) = WI*T
            T = Y+V*QTR(I)
            QTR(I) = QTR(I)+T*U1
            Y = Y+T*U2
            IF (IP1 .GT. P) GO TO 60
            DO 30 J = IP1, P
               T = W(J)+V*RMAT(IJ)
               RMAT(IJ) = RMAT(IJ)+T*U1
               W(J) = W(J)+T*U2
               IJ = IJ+J
 30            CONTINUE
            GO TO 60
C
C        *** AT THIS POINT WE MUST HAVE ABS(WI) .LE. ABS(RI)...
C
 40      WR = WI/RI
         T = - SQRT(ONE+WR**2)
         V = WR/(ONE-T)
         U1 = ONE/T-ONE
         U2 = V*U1
         RMAT(II) = RI*T
         T = QTR(I)+V*Y
         QTR(I) = QTR(I)+T*U1
         Y = Y+T*U2
         IF (IP1 .GT. P) GO TO 60
         DO 50 J = IP1, P
            T = RMAT(IJ)+V*W(J)
            RMAT(IJ) = RMAT(IJ)+T*U1
            W(J) = W(J)+T*U2
            IJ = IJ+J
 50         CONTINUE
 60      CONTINUE
 999  RETURN
      END
      SUBROUTINE LTVMUL(N, X, L, Y)
C
C  ***  COMPUTE  X = (L**T)*Y, WHERE  L  IS AN  N X N  LOWER
C  ***  TRIANGULAR MATRIX STORED COMPACTLY BY ROWS.  X AND Y MAY
C  ***  OCCUPY THE SAME STORAGE.  ***
C
C***REVISION HISTORY  (YYMMDD)
C   000330  Modified array declarations.  (JEC)
C
      INTEGER N
      REAL X(N), L(*), Y(N)
C     DIMENSION L(N*(N+1)/2)
      INTEGER I, IJ, I0, J
      REAL YI, ZERO
C/6
      DATA ZERO/0.E+0/
C/7
C     PARAMETER (ZERO=0.E+0)
C/
C
      I0 = 0
      DO 20 I = 1, N
         YI = Y(I)
         X(I) = ZERO
         DO 10 J = 1, I
              IJ = I0 + J
              X(J) = X(J) + YI*L(IJ)
 10           CONTINUE
         I0 = I0 + I
 20      CONTINUE
 999  RETURN
C  ***  LAST CARD OF LTVMUL FOLLOWS  ***
      END
      SUBROUTINE LSQUAR(N, A, L)
C
C  ***  COMPUTE  A = LOWER TRIANGLE OF  L*(L**T),  WITH BOTH
C  ***  L  AND  A  STORED COMPACTLY BY ROWS.  (BOTH MAY OCCUPY THE
C  ***  SAME STORAGE.
C
C***REVISION HISTORY  (YYMMDD)
C   000330  Modified array declarations.  (JEC)
C
C  ***  PARAMETERS  ***
C
      INTEGER N
      REAL A(*), L(*)
C     DIMENSION A(N*(N+1)/2), L(N*(N+1)/2)
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER I, II, IJ, IK, IP1, I0, J, JJ, JK, J0, K, NP1
      REAL T
C
      NP1 = N + 1
      I0 = N*(N+1)/2
      DO 30 II = 1, N
         I = NP1 - II
         IP1 = I + 1
         I0 = I0 - I
         J0 = I*(I+1)/2
         DO 20 JJ = 1, I
              J = IP1 - JJ
              J0 = J0 - J
              T = 0.0E0
              DO 10 K = 1, J
                   IK = I0 + K
                   JK = J0 + K
                   T = T + L(IK)*L(JK)
 10                CONTINUE
              IJ = I0 + J
              A(IJ) = T
 20           CONTINUE
 30      CONTINUE
 999  RETURN
      END
      SUBROUTINE OPRDS(L, LS, P, S, W, Y, Z)
C
C  ***  FOR I = 1..L, SET S = S + W(I)*Y(.,I)*(Z(.,I)**T), I.E.,
C  ***        ADD W(I) TIMES THE OUTER PRODUCT OF Y(.,I) AND Z(.,I).
C
      INTEGER L, LS, P
      REAL S(LS), W(L), Y(P,L), Z(P,L)
C     DIMENSION S(P*(P+1)/2)
C
      INTEGER I, J, K, M
      REAL WK, YI, ZERO
      DATA ZERO/0.E+0/
C
      DO 30 K = 1, L
         WK = W(K)
         IF (WK .EQ. ZERO) GO TO 30
         M = 1
         DO 20 I = 1, P
              YI = WK * Y(I,K)
              DO 10 J = 1, I
                   S(M) = S(M) + YI*Z(J,K)
                   M = M + 1
 10                CONTINUE
 20           CONTINUE
 30      CONTINUE
C
 999  RETURN
C  ***  LAST CARD OF OPRDS FOLLOWS  ***
      END
      SUBROUTINE GQTST(D, DIG, DIHDI, KA, L, P, STEP, V, W)
C
C  *** COMPUTE GOLDFELD-QUANDT-TROTTER STEP BY MORE-HEBDEN TECHNIQUE ***
C  ***  (NL2SOL VERSION 2.2), MODIFIED A LA MORE AND SORENSEN  ***
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER KA, P
      REAL D(P), DIG(P), DIHDI(*), L(*), V(21), STEP(P),
     1                 W(*)
C     DIMENSION DIHDI(P*(P+1)/2), L(P*(P+1)/2), W(4*P+7)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ***  PURPOSE  ***
C
C        GIVEN THE (COMPACTLY STORED) LOWER TRIANGLE OF A SCALED
C     HESSIAN (APPROXIMATION) AND A NONZERO SCALED GRADIENT VECTOR,
C     THIS SUBROUTINE COMPUTES A GOLDFELD-QUANDT-TROTTER STEP OF
C     APPROXIMATE LENGTH V(RADIUS) BY THE MORE-HEBDEN TECHNIQUE.  IN
C     OTHER WORDS, STEP IS COMPUTED TO (APPROXIMATELY) MINIMIZE
C     PSI(STEP) = (G**T)*STEP + 0.5*(STEP**T)*H*STEP  SUCH THAT THE
C     2-NORM OF D*STEP IS AT MOST (APPROXIMATELY) V(RADIUS), WHERE
C     G  IS THE GRADIENT,  H  IS THE HESSIAN, AND  D  IS A DIAGONAL
C     SCALE MATRIX WHOSE DIAGONAL IS STORED IN THE PARAMETER D.
C     (GQTST 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 LITVMU - APPLIES INVERSE-TRANSPOSE OF COMPACT LOWER TRIANG. MATRIX.
C LIVMUL - APPLIES INVERSE OF COMPACT LOWER TRIANG. MATRIX.
C LSQRT  - FINDS CHOLESKY FACTOR (OF COMPACTLY STORED LOWER TRIANG.).
C LSVMIN - RETURNS APPROX. TO MIN. SING. VALUE OF LOWER TRIANG. MATRIX.
C
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***REVISION HISTORY  (YYMMDD)
C   000330  Modified array declarations.  (JEC)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ***  LOCAL VARIABLES  ***
C
      LOGICAL RESTRT
      INTEGER DGGDMX, DIAG, DIAG0, DSTSAV, EMAX, EMIN, I, IM1, INC, IRC,
     1        J, K, KALIM, KAMIN, K1, LK0, PHIPIN, Q, Q0, UK0, X
      REAL ALPHAK, AKI, AKK, DELTA, DST, EPS, GTSTA, LK,
     1                 OLDPHI, PHI, PHIMAX, PHIMIN, PSIFAC, RAD, RADSQ,
     2                 ROOT, SI, SK, SW, T, TWOPSI, T1, T2, UK, WI
C
C     ***  CONSTANTS  ***
      REAL BIG, DGXFAC, EPSFAC, FOUR, HALF, KAPPA, NEGONE,
     1                 ONE, P001, SIX, THREE, TWO, ZERO
C
C  ***  INTRINSIC FUNCTIONS  ***
C/+
      REAL  ABS, AMAX1, AMIN1,  SQRT
C/
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      EXTERNAL LITVMU, LIVMUL, LSQRT, LSVMIN
      REAL LSVMIN
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.0E+0/, FOUR/4.0E+0/, HALF/0.5E+0/,
     1     KAPPA/2.0E+0/, NEGONE/-1.0E+0/, ONE/1.0E+0/, P001/1.0E-3/,
     2     SIX/6.0E+0/, THREE/3.0E+0/, TWO/2.0E+0/, ZERO/0.0E+0/
C/7
C     PARAMETER (EPSFAC=50.0E+0, FOUR=4.0E+0, HALF=0.5E+0,
C    1     KAPPA=2.0E+0, NEGONE=-1.0E+0, ONE=1.0E+0, P001=1.0E-3,
C    2     SIX=6.0E+0, THREE=3.0E+0, TWO=2.0E+0, ZERO=0.0E+0)
C     SAVE DGXFAC
C/
      DATA BIG/0.E+0/, DGXFAC/0.E+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) = SNRM2(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 =  ABS(DIHDI(I))
         IF (T1 .LT. T) T1 = T
 20      CONTINUE
      W(DGGDMX) = T1
C
C  ***  TRY ALPHA = 0  ***
C
 30   CALL LSQRT(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 LITVMU(IRC, W, L, W)
         T1 = SNRM2(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 = LSVMIN(P, L, W(Q), W(Q))
      IF (T .GE. ONE) GO TO 60
         IF (BIG .LE. ZERO) BIG = R1MACH(2)
         IF (V(DGNORM) .GE. T*T*BIG) GO TO 70
 60   CALL LIVMUL(P, W(Q), L, DIG)
      GTSTA = SDOT(P, W(Q),1,W(Q),1)
      V(NREDUC) = HALF * GTSTA
      CALL LITVMU(P, W(Q), L, W(Q))
      DST = SNRM2(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 =  ABS(DIHDI(K))
              WI = WI + T
              W(J) = W(J) + T
 80           CONTINUE
 90      W(I) = WI
         K = K + 1
 100     CONTINUE
C
C  ***  (UNDER-)ESTIMATE SMALLEST EIGENVALUE OF (D**-1)*H*(D**-1)  ***
C
      K = 1
      T1 = W(DIAG) - W(1)
      IF (P .LE. 1) GO TO 120
      DO 110 I = 2, P
         J = DIAG0 + I
         T = W(J) - W(I)
         IF (T .GE. T1) GO TO 110
              T1 = T
              K = I
 110     CONTINUE
C
 120  SK = W(K)
      J = DIAG0 + K
      AKK = W(J)
      K1 = K*(K-1)/2 + 1
      INC = 1
      T = ZERO
      DO 150 I = 1, P
         IF (I .EQ. K) GO TO 130
         AKI =  ABS(DIHDI(K1))
         SI = W(I)
         J = DIAG0 + I
         T1 = HALF * (AKK - W(J) + SI - AKI)
         T1 = T1 +  SQRT(T1*T1 + SK*AKI)
         IF (T .LT. T1) T = T1
         IF (I .LT. K) GO TO 140
 130     INC = I
 140     K1 = K1 + INC
 150     CONTINUE
C
      W(EMIN) = AKK - T
      UK = V(DGNORM)/RAD - W(EMIN)
      IF (V(DGNORM) .EQ. ZERO) UK = UK + P001 + P001*UK
C
C  ***  COMPUTE GERSCHGORIN (OVER-)ESTIMATE OF LARGEST EIGENVALUE  ***
C
      K = 1
      T1 = W(DIAG) + W(1)
      IF (P .LE. 1) GO TO 170
      DO 160 I = 2, P
         J = DIAG0 + I
         T = W(J) + W(I)
         IF (T .LE. T1) GO TO 160
              T1 = T
              K = I
 160     CONTINUE
C
 170  SK = W(K)
      J = DIAG0 + K
      AKK = W(J)
      K1 = K*(K-1)/2 + 1
      INC = 1
      T = ZERO
      DO 200 I = 1, P
         IF (I .EQ. K) GO TO 180
         AKI =  ABS(DIHDI(K1))
         SI = W(I)
         J = DIAG0 + I
         T1 = HALF * (W(J) + SI - AKI - AKK)
         T1 = T1 +  SQRT(T1*T1 + SK*AKI)
         IF (T .LT. T1) T = T1
         IF (I .LT. K) GO TO 190
 180     INC = I
 190     K1 = K1 + INC
 200     CONTINUE
C
      W(EMAX) = AKK + T
      LK = AMAX1(LK, V(DGNORM)/RAD - W(EMAX))
C
C     ***  ALPHAK = CURRENT VALUE OF ALPHA (SEE ALG. NOTES ABOVE).  WE
C     ***  USE MORE*S SCHEME FOR INITIALIZING IT.
      ALPHAK =  ABS(V(STPPAR)) * V(RAD0)/RAD
C
      IF (IRC .NE. 0) GO TO 210
C
C  ***  COMPUTE L0 FOR POSITIVE DEFINITE H  ***
C
      CALL LIVMUL(P, W, L, W(Q))
      T = SNRM2(P, W,1)
      W(PHIPIN) = DST / T / T
      LK = AMAX1(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 * AMAX1(P001,  SQRT(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 LSQRT(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 LITVMU(IRC, W, L, W)
      T1 = SNRM2(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 LIVMUL(P, W(Q), L, DIG)
      GTSTA = SDOT(P, W(Q),1,W(Q),1)
      CALL LITVMU(P, W(Q), L, W(Q))
      DST = SNRM2(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 AMIN1 IS NECESSARY BECAUSE OF RESTARTS  ***
      IF (PHI .LT. ZERO) UK = AMIN1(UK, ALPHAK)
C     *** KAMIN = 0 ONLY IFF THE GRADIENT VANISHES  ***
      IF (KAMIN .EQ. 0) GO TO 210
      CALL LIVMUL(P, W, L, W(Q))
      T1 = SNRM2(P, W,1)
      ALPHAK = ALPHAK  +  (PHI/T1) * (DST/T1) * (DST/RAD)
      LK = AMAX1(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 * ( ABS(ALPHAK)*DST*DST + GTSTA)
      GO TO 410
C
C
C  ***  RESTART WITH NEW RADIUS  ***
C
 290  IF (V(DST0) .LE. ZERO .OR. V(DST0) - RAD .GT. PHIMAX) GO TO 310
C
C     ***  PREPARE TO RETURN NEWTON STEP  ***
C
         RESTRT = .TRUE.
         KA = KA + 1
         K = 0
         DO 300 I = 1, P
              K = K + I
              J = DIAG0 + I
              DIHDI(K) = W(J)
 300          CONTINUE
         UK = NEGONE
         GO TO 30
C
 310  KAMIN = KA + 3
      IF (V(DGNORM) .EQ. ZERO) KAMIN = 0
      IF (KA .EQ. 0) GO TO 50
C
      DST = W(DSTSAV)
      ALPHAK =  ABS(V(STPPAR))
      PHI = DST - RAD
      T = V(DGNORM)/RAD
      UK = T - W(EMIN)
      IF (V(DGNORM) .EQ. ZERO) UK = UK + P001 + P001*UK
      IF (RAD .GT. V(RAD0)) GO TO 320
C
C        ***  SMALLER RADIUS  ***
         LK = ZERO
         IF (ALPHAK .GT. ZERO) LK = W(LK0)
         LK = AMAX1(LK, T - W(EMAX))
         IF (V(DST0) .GT. ZERO) LK = AMAX1(LK, (V(DST0)-RAD)*W(PHIPIN))
         GO TO 250
C
C     ***  BIGGER RADIUS  ***
 320  IF (ALPHAK .GT. ZERO) UK = AMIN1(UK, W(UK0))
      LK = AMAX1(ZERO, -V(DST0), T - W(EMAX))
      IF (V(DST0) .GT. ZERO) LK = AMAX1(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 + AMIN1(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 LSVMIN TO OBTAIN APPROXIMATE EIGENVECTOR CORRESPONDING
C  ***  TO SMALLEST EIGENVALUE OF (D**-1)*H*(D**-1).  LSVMIN RETURNS
C  ***  X AND W WITH  L*W = X.
C
 340  T = LSVMIN(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 LITVMU(P, W, L, W)
      T2 = ONE/SNRM2(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 = SDOT(P, W(Q),1,W,1)
      T1 = (RAD + DST) * (RAD - DST)
      ROOT =  SQRT(SW*SW + T1)
      IF (SW .LT. ZERO) ROOT = -ROOT
      SI = T1 / (SW + ROOT)
C
C  ***  THE ACTUAL TEST FOR THE SPECIAL CASE...
C
      IF ((T2*SI)**2 .LE. EPS*(DST**2 + ALPHAK*RADSQ)) GO TO 380
C
C  ***  UPDATE UPPER BOUND ON SMALLEST EIGENVALUE (WHEN NOT POSITIVE)
C  ***  (AS RECOMMENDED BY MORE AND SORENSEN) AND CONTINUE...
C
      IF (V(DST0) .LE. ZERO) V(DST0) = AMIN1(V(DST0), T2**2 - ALPHAK)
      LK = AMAX1(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 * R1MACH(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*SDOT(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) = SDOT(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 GQTST FOLLOWS  ***
      END
      SUBROUTINE LMSTEP(D, G, IERR, IPIVOT, KA, P, QTR, R, STEP, V, W)
C
C  ***  COMPUTE LEVENBERG-MARQUARDT STEP USING MORE-HEBDEN TECHNIQUE  **
C  ***  NL2SOL VERSION 2.2.  ***
C
C***REVISION HISTORY  (YYMMDD)
C   000330  Modified array declarations.  (JEC)
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER IERR, KA, P
      INTEGER IPIVOT(P)
      REAL D(P), G(P), QTR(P), R(*), STEP(P), V(21), W(*)
C     DIMENSION W(P*(P+5)/2 + 4)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ***  PURPOSE  ***
C
C        GIVEN THE R MATRIX FROM THE QR DECOMPOSITION OF A JACOBIAN
C     MATRIX, J, AS WELL AS Q-TRANSPOSE TIMES THE CORRESPONDING
C     RESIDUAL VECTOR, RESID, THIS SUBROUTINE COMPUTES A LEVENBERG-
C     MARQUARDT STEP OF APPROXIMATE LENGTH V(RADIUS) BY THE MORE-
C     TECHNIQUE.
C
C  ***  PARAMETER DESCRIPTION  ***
C
C      D (IN)  = THE SCALE VECTOR.
C      G (IN)  = THE GRADIENT VECTOR (J**T)*R.
C   IERR (I/O) = RETURN CODE FROM QRFACT OR QRFGS -- 0 MEANS R HAS
C             FULL RANK.
C IPIVOT (I/O) = PERMUTATION ARRAY FROM QRFACT OR QRFGS, WHICH COMPUTE
C             QR DECOMPOSITIONS WITH COLUMN PIVOTING.
C     KA (I/O).  KA .LT. 0 ON INPUT MEANS THIS IS THE FIRST CALL ON
C             LMSTEP FOR THE CURRENT R AND QTR.  ON OUTPUT KA CON-
C             TAINS THE NUMBER OF HEBDEN ITERATIONS NEEDED TO DETERMINE
C             STEP.  KA = 0 MEANS A GAUSS-NEWTON STEP.
C      P (IN)  = NUMBER OF PARAMETERS.
C    QTR (IN)  = (Q**T)*RESID = Q-TRANSPOSE TIMES THE RESIDUAL VECTOR.
C      R (IN)  = THE R MATRIX, STORED COMPACTLY BY COLUMNS.
C   STEP (OUT) = THE LEVENBERG-MARQUARDT STEP COMPUTED.
C      V (I/O) CONTAINS VARIOUS CONSTANTS AND VARIABLES DESCRIBED BELOW.
C      W (I/O) = WORKSPACE OF LENGTH P*(P+5)/2 + 4.
C
C  ***  ENTRIES IN V  ***
C
C V(DGNORM) (I/O) = 2-NORM OF (D**-1)*G.
C V(DSTNRM) (I/O) = 2-NORM OF D*STEP.
C V(DST0)   (I/O) = 2-NORM OF GAUSS-NEWTON STEP (FOR NONSING. J).
C V(EPSLON) (IN) = MAX. REL. ERROR ALLOWED IN TWONORM(R)**2 MINUS
C             TWONORM(R - J*STEP)**2.  (SEE ALGORITHM NOTES BELOW.)
C V(GTSTEP) (OUT) = INNER PRODUCT BETWEEN G AND STEP.
C V(NREDUC) (OUT) = HALF THE REDUCTION IN THE SUM OF SQUARES PREDICTED
C             FOR A GAUSS-NEWTON STEP.
C V(PHMNFC) (IN)  = TOL. (TOGETHER WITH V(PHMXFC)) FOR ACCEPTING STEP
C             (MORE*S SIGMA).  THE ERROR V(DSTNRM) - V(RADIUS) MUST LIE
C             BETWEEN V(PHMNFC)*V(RADIUS) AND V(PHMXFC)*V(RADIUS).
C V(PHMXFC) (IN)  (SEE V(PHMNFC).)
C V(PREDUC) (OUT) = HALF THE REDUCTION IN THE SUM OF SQUARES PREDICTED
C             BY THE STEP RETURNED.
C V(RADIUS) (IN)  = RADIUS OF CURRENT (SCALED) TRUST REGION.
C V(RAD0)   (I/O) = VALUE OF V(RADIUS) FROM PREVIOUS CALL.
C V(STPPAR) (I/O) = MARQUARDT PARAMETER (OR ITS NEGATIVE IF THE SPECIAL
C             CASE MENTIONED BELOW IN THE ALGORITHM NOTES OCCURS).
C
C NOTE -- SEE DATA STATEMENT BELOW FOR VALUES OF ABOVE SUBSCRIPTS.
C
C  ***  USAGE NOTES  ***
C
C     IF IT IS DESIRED TO RECOMPUTE STEP USING A DIFFERENT VALUE OF
C     V(RADIUS), THEN THIS ROUTINE MAY BE RESTARTED BY CALLING IT
C     WITH ALL PARAMETERS UNCHANGED EXCEPT V(RADIUS).  (THIS EXPLAINS
C     WHY MANY PARAMETERS ARE LISTED AS I/O).  ON AN INTIIAL CALL (ONE
C     WITH KA = -1), THE CALLER NEED ONLY HAVE INITIALIZED D, G, KA, P,
C     QTR, R, V(EPSLON), V(PHMNFC), V(PHMXFC), V(RADIUS), AND V(RAD0).
C
C  ***  APPLICATION AND USAGE RESTRICTIONS  ***
C
C     THIS ROUTINE IS CALLED AS PART OF THE NL2SOL (NONLINEAR LEAST-
C     SQUARES) PACKAGE (REF. 1).
C
C  ***  ALGORITHM NOTES  ***
C
C     THIS CODE IMPLEMENTS THE STEP COMPUTATION SCHEME DESCRIBED IN
C     REFS. 2 AND 4.  FAST GIVENS TRANSFORMATIONS (SEE REF. 3, PP. 60-
C     62) ARE USED TO COMPUTE STEP WITH A NONZERO MARQUARDT PARAMETER.
C        A SPECIAL CASE OCCURS IF J IS (NEARLY) SINGULAR AND V(RADIUS)
C     IS SUFFICIENTLY LARGE.  IN THIS CASE THE STEP RETURNED IS SUCH
C     THAT  TWONORM(R)**2 - TWONORM(R - J*STEP)**2  DIFFERS FROM ITS
C     OPTIMAL VALUE BY LESS THAN V(EPSLON) TIMES THIS OPTIMAL VALUE,
C     WHERE J AND R DENOTE THE ORIGINAL JACOBIAN AND RESIDUAL.  (SEE
C     REF. 2 FOR MORE DETAILS.)
C
C  ***  FUNCTIONS AND SUBROUTINES CALLED  ***
C
C LITVMU - APPLY INVERSE-TRANSPOSE OF COMPACT LOWER TRIANG. MATRIX.
C LIVMUL - APPLY INVERSE OF COMPACT 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.  LAWSON, C.L., AND HANSON, R.J. (1974), SOLVING LEAST SQUARES
C             PROBLEMS, PRENTICE-HALL, ENGLEWOOD CLIFFS, N.J.
C 4.  MORE, J.J. (1978), THE LEVENBERG-MARQUARDT ALGORITHM, IMPLEMEN-
C             TATION AND THEORY, PP.105-116 OF SPRINGER LECTURE NOTES
C             IN MATHEMATICS NO. 630, EDITED BY G.A. WATSON, SPRINGER-
C             VERLAG, BERLIN AND NEW YORK.
C
C  ***  GENERAL  ***
C
C     CODED BY DAVID M. GAY.
C     THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH
C     SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS
C     MCS-7600324, DCR75-10143, 76-14311DSS, MCS76-11989, AND
C     MCS-7906671.
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER DSTSAV, I, IP1, I1, J1, K, KALIM, L, LK0, PHIPIN,
     1        PP1O2, RES, RES0, RMAT, RMAT0, UK0
      REAL A, ADI, ALPHAK, B, DFACSQ, DST, DTOL, D1, D2,
     1                 LK, OLDPHI, PHI, PHIMAX, PHIMIN, PSIFAC, RAD,
     2                 SI, SJ, SQRTAK, T, TWOPSI, UK, WL
C
C     ***  CONSTANTS  ***
      REAL DFAC, EIGHT, HALF, NEGONE, ONE, P001, THREE,
     1                 TTOL, ZERO
      REAL BIG
C
C  ***  INTRINSIC FUNCTIONS  ***
C/+
      INTEGER IABS
      REAL  ABS, AMAX1, AMIN1,  SQRT
C/
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      EXTERNAL  LITVMU, LIVMUL, LSVMIN
      REAL  LSVMIN
C
C  ***  SUBSCRIPTS FOR V  ***
C
      INTEGER DGNORM, DSTNRM, DST0, EPSLON, GTSTEP, NREDUC, PHMNFC,
     1        PHMXFC, PREDUC, RADIUS, RAD0, STPPAR
C/6
C     DATA DGNORM/1/, DSTNRM/2/, DST0/3/, EPSLON/19/, GTSTEP/4/,
C    1     NREDUC/6/, PHMNFC/20/, PHMXFC/21/, PREDUC/7/, RADIUS/8/,
C    2     RAD0/9/, STPPAR/5/
C/7
      PARAMETER (DGNORM=1, DSTNRM=2, DST0=3, EPSLON=19, GTSTEP=4,
     1           NREDUC=6, PHMNFC=20, PHMXFC=21, PREDUC=7, RADIUS=8,
     2           RAD0=9, STPPAR=5)
C/
C
C/6
C     DATA DFAC/256.E+0/, EIGHT/8.E+0/, HALF/0.5E+0/, NEGONE/-1.E+0/,
C    1     ONE/1.E+0/, P001/1.E-3/, THREE/3.E+0/, TTOL/2.5E+0/,
C    2     ZERO/0.E+0/
C/7
      PARAMETER (DFAC=256.E+0, EIGHT=8.E+0, HALF=0.5E+0, NEGONE=-1.E+0,
     1     ONE=1.E+0, P001=1.E-3, THREE=3.E+0, TTOL=2.5E+0,
     2     ZERO=0.E+0)
      SAVE BIG
C/
      DATA BIG/0.E+0/
C
C  ***  BODY  ***
C
C     ***  FOR USE IN RECOMPUTING STEP, THE FINAL VALUES OF LK AND UK,
C     ***  THE INVERSE DERIVATIVE OF MORE*S PHI AT 0 (FOR NONSING. J)
C     ***  AND THE VALUE RETURNED AS V(DSTNRM) ARE STORED AT W(LK0),
C     ***  W(UK0), W(PHIPIN), AND W(DSTSAV) RESPECTIVELY.
      LK0 = P + 1
      PHIPIN = LK0 + 1
      UK0 = PHIPIN + 1
      DSTSAV = UK0 + 1
      RMAT0 = DSTSAV
C     ***  A COPY OF THE R-MATRIX FROM THE QR DECOMPOSITION OF J IS
C     ***  STORED IN W STARTING AT W(RMAT), AND A COPY OF THE RESIDUAL
C     ***  VECTOR IS STORED IN W STARTING AT W(RES).  THE LOOPS BELOW
C     ***  THAT UPDATE THE QR DECOMP. FOR A NONZERO MARQUARDT PARAMETER
C     ***  WORK ON THESE COPIES.
      RMAT = RMAT0 + 1
      PP1O2 = P * (P + 1) / 2
      RES0 = PP1O2 + RMAT0
      RES = RES0 + 1
      RAD = V(RADIUS)
      IF (RAD .GT. ZERO)
     1   PSIFAC = V(EPSLON)/((EIGHT*(V(PHMNFC) + ONE) + THREE) * RAD**2)
      IF (BIG .LE. ZERO) BIG = R1MACH(2)
      PHIMAX = V(PHMXFC) * RAD
      PHIMIN = V(PHMNFC) * RAD
C     ***  DTOL, DFAC, AND DFACSQ ARE USED IN RESCALING THE FAST GIVENS
C     ***  REPRESENTATION OF THE UPDATED QR DECOMPOSITION.
      DTOL = ONE/DFAC
      DFACSQ = DFAC*DFAC
C     ***  OLDPHI IS USED TO DETECT LIMITS OF NUMERICAL ACCURACY.  IF
C     ***  WE RECOMPUTE STEP AND IT DOES NOT CHANGE, THEN WE ACCEPT IT.
      OLDPHI = ZERO
      LK = ZERO
      UK = ZERO
      KALIM = KA + 12
C
C  ***  START OR RESTART, DEPENDING ON KA  ***
C
      IF (KA) 10, 20, 370
C
C  ***  FRESH START -- COMPUTE V(NREDUC)  ***
C
 10   KA = 0
      KALIM = 12
      K = P
      IF (IERR .NE. 0) K = IABS(IERR) - 1
      V(NREDUC) = HALF*SDOT(K, QTR,1,QTR,1)
C
C  ***  SET UP TO TRY INITIAL GAUSS-NEWTON STEP  ***
C
 20   V(DST0) = NEGONE
      IF (IERR .NE. 0) GO TO 90
      T = LSVMIN(P, R, STEP, W(RES))
      IF (T .GE. ONE) GO TO 30
         IF (SNRM2(P, QTR,1) .GE. BIG*T) GO TO 90
C
C  ***  COMPUTE GAUSS-NEWTON STEP  ***
C
C     ***  NOTE -- THE R-MATRIX IS STORED COMPACTLY BY COLUMNS IN
C     ***  R(1), R(2), R(3), ...  IT IS THE TRANSPOSE OF A
C     ***  LOWER TRIANGULAR MATRIX STORED COMPACTLY BY ROWS, AND WE
C     ***  TREAT IT AS SUCH WHEN USING LITVMU AND LIVMUL.
 30   CALL LITVMU(P, W, R, QTR)
C     ***  TEMPORARILY STORE PERMUTED -D*STEP IN STEP.
      DO 60 I = 1, P
         J1 = IPIVOT(I)
         STEP(I) = D(J1)*W(I)
 60      CONTINUE
      DST = SNRM2(P, STEP,1)
      V(DST0) = DST
      PHI = DST - RAD
      IF (PHI .LE. PHIMAX) GO TO 410
C     ***  IF THIS IS A RESTART, GO TO 110  ***
      IF (KA .GT. 0) GO TO 110
C
C  ***  GAUSS-NEWTON STEP WAS UNACCEPTABLE.  COMPUTE L0  ***
C
      DO 70 I = 1, P
         J1 = IPIVOT(I)
         STEP(I) = D(J1)*(STEP(I)/DST)
 70      CONTINUE
      CALL LIVMUL(P, STEP, R, STEP)
      T = ONE / SNRM2(P, STEP,1)
      W(PHIPIN) = (T/DST)*T
      LK = PHI*W(PHIPIN)
C
C  ***  COMPUTE U0  ***
C
 90   DO 100 I = 1, P
 100     W(I) = G(I)/D(I)
      V(DGNORM) = SNRM2(P, W,1)
      UK = V(DGNORM)/RAD
      IF (UK .LE. ZERO) GO TO 390
C
C     ***  ALPHAK WILL BE USED AS THE CURRENT MARQUARDT PARAMETER.  WE
C     ***  USE MORE*S SCHEME FOR INITIALIZING IT.
      ALPHAK =  ABS(V(STPPAR)) * V(RAD0)/RAD
C
C
C  ***  TOP OF LOOP -- INCREMENT KA, COPY R TO RMAT, QTR TO RES  ***
C
 110  KA = KA + 1
      CALL SCOPY(PP1O2, R,1,W(RMAT),1)
      CALL SCOPY(P, QTR,1,W(RES),1)
C
C  ***  SAFEGUARD ALPHAK AND INITIALIZE FAST GIVENS SCALE VECTOR.  ***
C
      IF (ALPHAK .LE. ZERO .OR. ALPHAK .LT. LK .OR. ALPHAK .GE. UK)
     1             ALPHAK = UK * AMAX1(P001,  SQRT(LK/UK))
      IF(ALPHAK .LE. ZERO) ALPHAK = HALF * UK
      SQRTAK =  SQRT(ALPHAK)
      DO 120 I = 1, P
 120     W(I) = ONE
C
C  ***  ADD ALPHAK*D AND UPDATE QR DECOMP. USING FAST GIVENS TRANS.  ***
C
      DO 270 I = 1, P
C        ***  GENERATE, APPLY 1ST GIVENS TRANS. FOR ROW I OF ALPHAK*D.
C        ***  (USE STEP TO STORE TEMPORARY ROW)  ***
         L = I*(I+1)/2 + RMAT0
         WL = W(L)
         D2 = ONE
         D1 = W(I)
         J1 = IPIVOT(I)
         ADI = SQRTAK*D(J1)
         IF (ADI .GE.  ABS(WL)) GO TO 150
 130     A = ADI/WL
         B = D2*A/D1
         T = A*B + ONE
         IF (T .GT. TTOL) GO TO 150
         W(I) = D1/T
         D2 = D2/T
         W(L) = T*WL
         A = -A
         DO 140 J1 = I, P
              L = L + J1
              STEP(J1) = A*W(L)
 140          CONTINUE
         GO TO 170
C
 150     B = WL/ADI
         A = D1*B/D2
         T = A*B + ONE
         IF (T .GT. TTOL) GO TO 130
         W(I) = D2/T
         D2 = D1/T
         W(L) = T*ADI
         DO 160 J1 = I, P
              L = L + J1
              WL = W(L)
              STEP(J1) = -WL
              W(L) = A*WL
 160          CONTINUE
C
 170     IF (I .EQ. P) GO TO 280
C
C        ***  NOW USE GIVENS TRANS. TO ZERO ELEMENTS OF TEMP. ROW  ***
C
         IP1 = I + 1
         DO 260 I1 = IP1, P
              L = I1*(I1+1)/2 + RMAT0
              WL = W(L)
              SI = STEP(I1-1)
              D1 = W(I1)
C
C             ***  RESCALE ROW I1 IF NECESSARY  ***
C
              IF (D1 .GE. DTOL) GO TO 190
                   D1 = D1*DFACSQ
                   WL = WL/DFAC
                   K = L
                   DO 180 J1 = I1, P
                        K = K + J1
                        W(K) = W(K)/DFAC
 180                    CONTINUE
C
C             ***  USE GIVENS TRANS. TO ZERO NEXT ELEMENT OF TEMP. ROW
C
 190          IF ( ABS(SI) .GT.  ABS(WL)) GO TO 220
              IF (SI .EQ. ZERO) GO TO 260
 200          A = SI/WL
              B = D2*A/D1
              T = A*B + ONE
              IF (T .GT. TTOL) GO TO 220
              W(L) = T*WL
              W(I1) = D1/T
              D2 = D2/T
              DO 210 J1 = I1, P
                   L = L + J1
                   WL = W(L)
                   SJ = STEP(J1)
                   W(L) = WL + B*SJ
                   STEP(J1) = SJ - A*WL
 210               CONTINUE
              GO TO 240
C
 220          B = WL/SI
              A = D1*B/D2
              T = A*B + ONE
              IF (T .GT. TTOL) GO TO 200
              W(I1) = D2/T
              D2 = D1/T
              W(L) = T*SI
              DO 230 J1 = I1, P
                   L = L + J1
                   WL = W(L)
                   SJ = STEP(J1)
                   W(L) = A*WL + SJ
                   STEP(J1) = B*SJ - WL
 230               CONTINUE
C
C             ***  RESCALE TEMP. ROW IF NECESSARY  ***
C
 240          IF (D2 .GE. DTOL) GO TO 260
                   D2 = D2*DFACSQ
                   DO 250 K = I1, P
 250                    STEP(K) = STEP(K)/DFAC
 260          CONTINUE
 270     CONTINUE
C
C  ***  COMPUTE STEP  ***
C
 280  CALL LITVMU(P, W(RES), W(RMAT), W(RES))
C     ***  RECOVER STEP AND STORE PERMUTED -D*STEP AT W(RES)  ***
      DO 290 I = 1, P
         J1 = IPIVOT(I)
         K = RES0 + I
         T = W(K)
         STEP(J1) = -T
         W(K) = T*D(J1)
 290     CONTINUE
      DST = SNRM2(P, W(RES),1)
      PHI = DST - RAD
      IF (PHI .LE. PHIMAX .AND. PHI .GE. PHIMIN) GO TO 430
      IF (OLDPHI .EQ. PHI) GO TO 430
      OLDPHI = PHI
C
C  ***  CHECK FOR (AND HANDLE) SPECIAL CASE  ***
C
      IF (PHI .GT. ZERO) GO TO 310
         IF (KA .GE. KALIM) GO TO 430
              TWOPSI = ALPHAK*DST*DST - SDOT(P, STEP,1,G,1)
              IF (ALPHAK .GE. TWOPSI*PSIFAC) GO TO 310
                   V(STPPAR) = -ALPHAK
                   GO TO 440
C
C  ***  UNACCEPTABLE STEP -- UPDATE LK, UK, ALPHAK, AND TRY AGAIN  ***
C
 300  IF (PHI .LT. ZERO) UK = AMIN1(UK, ALPHAK)
      GO TO 320
 310  IF (PHI .LT. ZERO) UK = ALPHAK
 320  DO 330 I = 1, P
         J1 = IPIVOT(I)
         K = RES0 + I
         STEP(I) = D(J1) * (W(K)/DST)
 330     CONTINUE
      CALL LIVMUL(P, STEP, W(RMAT), STEP)
      DO 340 I = 1, P
 340     STEP(I) = STEP(I) /  SQRT(W(I))
      T = ONE / SNRM2(P, STEP,1)
      ALPHAK = ALPHAK + T*PHI*T/RAD
      LK = AMAX1(LK, ALPHAK)
      GO TO 110
C
C  ***  RESTART  ***
C
 370  LK = W(LK0)
      UK = W(UK0)
      IF (V(DST0) .GT. ZERO .AND. V(DST0) - RAD .LE. PHIMAX) GO TO 20
      ALPHAK =  ABS(V(STPPAR))
      DST = W(DSTSAV)
      PHI = DST - RAD
      T = V(DGNORM)/RAD
      IF (RAD .GT. V(RAD0)) GO TO 380
C
C        ***  SMALLER RADIUS  ***
         UK = T
         IF (ALPHAK .LE. ZERO) LK = ZERO
         IF (V(DST0) .GT. ZERO) LK = AMAX1(LK, (V(DST0)-RAD)*W(PHIPIN))
         GO TO 300
C
C     ***  BIGGER RADIUS  ***
 380  IF (ALPHAK .LE. ZERO .OR. UK .GT. T) UK = T
      LK = ZERO
      IF (V(DST0) .GT. ZERO) LK = AMAX1(LK, (V(DST0)-RAD)*W(PHIPIN))
      GO TO 300
C
C  ***  SPECIAL CASE -- RAD .LE. 0 OR (G = 0 AND J IS SINGULAR)  ***
C
 390  V(STPPAR) = ZERO
      DST = ZERO
      LK = ZERO
      UK = ZERO
      V(GTSTEP) = ZERO
      V(PREDUC) = ZERO
      DO 400 I = 1, P
 400     STEP(I) = ZERO
      GO TO 450
C
C  ***  ACCEPTABLE GAUSS-NEWTON STEP -- RECOVER STEP FROM W  ***
C
 410  ALPHAK = ZERO
      DO 420 I = 1, P
         J1 = IPIVOT(I)
         STEP(J1) = -W(I)
 420     CONTINUE
C
C  ***  SAVE VALUES FOR USE IN A POSSIBLE RESTART  ***
C
 430  V(STPPAR) = ALPHAK
 440  V(GTSTEP) = AMIN1(SDOT(P, STEP,1,G,1),ZERO)
      V(PREDUC) = HALF * (ALPHAK*DST*DST - V(GTSTEP))
 450  V(DSTNRM) = DST
      W(DSTSAV) = DST
      W(LK0) = LK
      W(UK0) = UK
      V(RAD0) = RAD
C
 999  RETURN
C
C  ***  LAST CARD OF LMSTEP FOLLOWS  ***
      END
      SUBROUTINE LSQRT(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***REVISION HISTORY  (YYMMDD)
C   000330  Modified array declarations.  (JEC)
C
C  ***  PARAMETERS  ***
C
      INTEGER N1, N, IRC
      REAL L(*), A(*)
C     DIMENSION L(N*(N+1)/2), A(N*(N+1)/2)
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER I, IJ, IK, IM1, I0, J, JK, JM1, J0, K
      REAL T, TD, ZERO
C
C  ***  INTRINSIC FUNCTIONS  ***
C/+
      REAL  SQRT
C/
C/6
      DATA ZERO/0.E+0/
C/7
C     PARAMETER (ZERO=0.E+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) =  SQRT(T)
 50      CONTINUE
C
      IRC = 0
      GO TO 999
C
 60   L(I0) = T
      IRC = I
C
 999  RETURN
C
C  ***  LAST CARD OF LSQRT  ***
      END
      LOGICAL FUNCTION STOPX(IDUMMY)
C     *****PARAMETERS...
      INTEGER IDUMMY
C
C     ..................................................................
C
C     *****PURPOSE...
C     THIS FUNCTION MAY SERVE AS THE STOPX (ASYNCHRONOUS INTERRUPTION)
C     FUNCTION FOR THE NL2SOL (NONLINEAR LEAST-SQUARES) PACKAGE AT
C     THOSE INSTALLATIONS WHICH DO NOT WISH TO IMPLEMENT A
C     DYNAMIC STOPX.
C
C     *****ALGORITHM NOTES...
C     AT INSTALLATIONS WHERE THE NL2SOL SYSTEM IS USED
C     INTERACTIVELY, THIS DUMMY STOPX SHOULD BE REPLACED BY A
C     FUNCTION THAT RETURNS .TRUE. IF AND ONLY IF THE INTERRUPT
C     (BREAK) KEY HAS BEEN PRESSED SINCE THE LAST CALL ON STOPX.
C
C     ..................................................................
C
      STOPX = .FALSE.
      RETURN
      END
      SUBROUTINE PARCK(ALG, D, IV, LIV, LV, N, V)
C
C  ***  CHECK ***SOL (VERSION 2.3) PARAMETERS, PRINT CHANGED VALUES  ***
C
C  ***  ALG = 1 FOR REGRESSION, ALG = 2 FOR GENERAL UNCONSTRAINED OPT.
C
      INTEGER ALG, LIV, LV, N
      INTEGER IV(LIV)
      REAL D(N), V(LV)
C
      EXTERNAL VDFLT
C VDFLT  -- 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/
      REAL BIG, MACHEP, TINY, VK, VM(34), VX(34), ZERO
C
C  ***  IV AND V SUBSCRIPTS  ***
C
      INTEGER ALGSAV, DINIT, DTYPE, DTYPE0, EPSLON, INITS, IVNEED,
     1        LASTIV, LASTV, LMAT, NEXTIV, NEXTV, NVDFLT, OLDN,
     2        PARPRT, PARSAV, PERM, PRUNIT, VNEED
C
C
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.E+0/, MACHEP/-1.E+0/, TINY/1.E+0/, ZERO/0.E+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.0E-3/, VM(2)/-0.99E+0/, VM(3)/1.0E-3/, VM(4)/1.0E-2/,
     1     VM(5)/1.2E+0/, VM(6)/1.E-2/, VM(7)/1.2E+0/, VM(8)/0.E+0/,
     2     VM(9)/0.E+0/, VM(10)/1.E-3/, VM(11)/-1.E+0/, VM(15)/0.E+0/,
     3     VM(16)/0.E+0/, VM(19)/0.E+0/, VM(20)/-10.E+0/, VM(21)/0.E+0/,
     4     VM(22)/0.E+0/, VM(23)/0.E+0/, VM(27)/1.01E+0/,
     5     VM(28)/1.E+10/, VM(30)/0.E+0/, VM(31)/0.E+0/, VM(32)/0.E+0/,
     6     VM(34)/0.E+0/
      DATA VX(1)/0.9E+0/, VX(2)/-1.E-3/, VX(3)/1.E+1/, VX(4)/0.8E+0/,
     1     VX(5)/1.E+2/, VX(6)/0.8E+0/, VX(7)/1.E+2/, VX(8)/0.5E+0/,
     2     VX(9)/0.5E+0/, VX(10)/1.E+0/, VX(11)/1.E+0/, VX(14)/0.1E+0/,
     3     VX(15)/1.E+0/, VX(16)/1.E+0/, VX(19)/1.E+0/, VX(23)/1.E+0/,
     4     VX(24)/1.E+0/, VX(25)/1.E+0/, VX(26)/1.E+0/, VX(27)/1.E+10/,
     5     VX(29)/1.E+0/, VX(31)/1.E+0/, VX(32)/1.E+0/, VX(33)/1.E+0/,
     6     VX(34)/1.E+0/
C
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 DEFLT(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(/39H THE FIRST PARAMETER TO DEFLT 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 VDFLT(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 = R1MACH(1)
         MACHEP = R1MACH(4)
         BIG = R1MACH(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) = SQRT(R1MACH(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) =,E11.3,7H SHOULD,
     1               11H BE BETWEEN,E11.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) =,E11.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) =,E15.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 SCOPY(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 PARCK FOLLOWS  ***
      END
      SUBROUTINE LIVMUL(N, X, L, Y)
C
C  ***  SOLVE  L*X = Y, WHERE  L  IS AN  N X N  LOWER TRIANGULAR
C  ***  MATRIX STORED COMPACTLY BY ROWS.  X AND Y MAY OCCUPY THE SAME
C  ***  STORAGE.  ***
C
C***REVISION HISTORY  (YYMMDD)
C   000330  Modified array declarations.  (JEC)
C
      INTEGER N
      REAL X(N), L(*), Y(N)
      INTEGER I, J, K
      REAL T, ZERO
C/6
      DATA ZERO/0.E+0/
C/7
C     PARAMETER (ZERO=0.E+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 = SDOT(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 LIVMUL FOLLOWS  ***
      END
      SUBROUTINE LVMUL(N, X, L, Y)
C
C  ***  COMPUTE  X = L*Y, WHERE  L  IS AN  N X N  LOWER TRIANGULAR
C  ***  MATRIX STORED COMPACTLY BY ROWS.  X AND Y MAY OCCUPY THE SAME
C  ***  STORAGE.  ***
C
C***REVISION HISTORY  (YYMMDD)
C   000330  Modified array declarations.  (JEC)
C
      INTEGER N
      REAL X(N), L(*), Y(N)
C     DIMENSION L(N*(N+1)/2)
      INTEGER I, II, IJ, I0, J, NP1
      REAL T, ZERO
C/6
      DATA ZERO/0.E+0/
C/7
C     PARAMETER (ZERO=0.E+0)
C/
C
      NP1 = N + 1
      I0 = N*(N+1)/2
      DO 20 II = 1, N
         I = NP1 - II
         I0 = I0 - I
         T = ZERO
         DO 10 J = 1, I
              IJ = I0 + J
              T = T + L(IJ)*Y(J)
 10           CONTINUE
         X(I) = T
 20      CONTINUE
 999  RETURN
C  ***  LAST CARD OF LVMUL FOLLOWS  ***
      END
      SUBROUTINE ITSUM(D, G, IV, LIV, LV, P, V, X)
C
C  ***  PRINT ITERATION SUMMARY FOR ***SOL (VERSION 2.3)  ***
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER LIV, LV, P
      INTEGER IV(LIV)
      REAL D(P), G(P), V(LV), X(P)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER ALG, I, IV1, M, NF, NG, OL, PU
C/6
      REAL MODEL1(6), MODEL2(6)
C/7
C     CHARACTER*4 MODEL1(6), MODEL2(6)
C/
      REAL NRELDF, OLDF, PRELDF, RELDF, ZERO
C
C  ***  INTRINSIC FUNCTIONS  ***
C/+
      INTEGER IABS
      REAL  ABS, AMAX1
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.E+0/
C/7
C     PARAMETER (ZERO=0.E+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 = AMAX1( ABS(V(F0)),  ABS(V(F)))
      IF (OLDF .LE. ZERO) GO TO 20
         RELDF = V(FDIF) / OLDF
         PRELDF = V(PREDUC) / OLDF
 20   IF (OL .GT. 0) GO TO 60
C
C        ***  PRINT SHORT SUMMARY LINE  ***
C
         IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 1) WRITE(PU,30)
 30   FORMAT(/10H   IT   NF,6X,1HF,7X,5HRELDF,3X,6HPRELDF,3X,5HRELDX,
     1       2X,13HMODEL  STPPAR)
         IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 2) WRITE(PU,40)
 40   FORMAT(/11H    IT   NF,7X,1HF,8X,5HRELDF,4X,6HPRELDF,4X,5HRELDX,
     1       3X,6HSTPPAR)
         IV(NEEDHD) = 0
         IF (ALG .EQ. 2) GO TO 50
         M = IV(SUSED)
         WRITE(PU,100) IV(NITER), NF, V(F), RELDF, PRELDF, V(RELDX),
     1                 MODEL1(M), MODEL2(M), V(STPPAR)
         GO TO 120
C
 50      WRITE(PU,110) IV(NITER), NF, V(F), RELDF, PRELDF, V(RELDX),
     1                 V(STPPAR)
         GO TO 120
C
C     ***  PRINT LONG SUMMARY LINE  ***
C
 60   IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 1) WRITE(PU,70)
 70   FORMAT(/11H    IT   NF,6X,1HF,7X,5HRELDF,3X,6HPRELDF,3X,5HRELDX,
     1       2X,13HMODEL  STPPAR,2X,6HD*STEP,2X,7HNPRELDF)
      IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 2) WRITE(PU,80)
 80   FORMAT(/11H    IT   NF,7X,1HF,8X,5HRELDF,4X,6HPRELDF,4X,5HRELDX,
     1       3X,6HSTPPAR,3X,6HD*STEP,3X,7HNPRELDF)
      IV(NEEDHD) = 0
      NRELDF = ZERO
      IF (OLDF .GT. ZERO) NRELDF = V(NREDUC) / OLDF
      IF (ALG .EQ. 2) GO TO 90
      M = IV(SUSED)
      WRITE(PU,100) IV(NITER), NF, V(F), RELDF, PRELDF, V(RELDX),
     1             MODEL1(M), MODEL2(M), V(STPPAR), V(DSTNRM), NRELDF
      GO TO 120
C
 90   WRITE(PU,110) IV(NITER), NF, V(F), RELDF, PRELDF,
     1             V(RELDX), V(STPPAR), V(DSTNRM), NRELDF
 100  FORMAT(I6,I5,E10.3,2E9.2,E8.1,A3,A4,2E8.1,E9.2)
 110  FORMAT(I6,I5,E11.3,2E10.2,3E9.1,E10.2)
C
 120  IF (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 ITSUM  ***
C
 390  IF (IV(X0PRT) .NE. 0) WRITE(PU,400) (I, X(I), D(I), I = 1, P)
 400  FORMAT(/23H     I     INITIAL X(I),8X,4HD(I)//(1X,I5,E17.6,E14.3))
      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,E10.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 = AMAX1( ABS(V(F0)),  ABS(V(F)))
         PRELDF = ZERO
         NRELDF = ZERO
         IF (OLDF .LE. ZERO) GO TO 440
              PRELDF = V(PREDUC) / OLDF
              NRELDF = V(NREDUC) / OLDF
 440     NF = IV(NFCALL) - IV(NFCOV)
         NG = IV(NGCALL) - IV(NGCOV)
         WRITE(PU,450) V(F), V(RELDX), NF, NG, PRELDF, NRELDF
 450  FORMAT(/9H FUNCTION,E17.6,8H   RELDX,E17.3/12H FUNC. EVALS,
     1   I8,9X,11HGRAD. EVALS,I8/7H PRELDF,E16.3,6X,7HNPRELDF,E15.3)
C
         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,E16.6,2E14.3)
      GO TO 999
C
 520  WRITE(PU,530)
 530  FORMAT(/24H INCONSISTENT DIMENSIONS)
 999  RETURN
C  ***  LAST CARD OF ITSUM FOLLOWS  ***
      END
      SUBROUTINE SLUPDT(A, COSMIN, P, SIZE, STEP, U, W, WCHMTD, WSCALE,
     1                  Y)
C
C  ***  UPDATE SYMMETRIC  A  SO THAT  A * STEP = Y  ***
C  ***  (LOWER TRIANGLE OF  A  STORED ROWWISE       ***
C
C***REVISION HISTORY  (YYMMDD)
C   000330  Modified array declarations.  (JEC)
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER P
      REAL A(*), COSMIN, SIZE, STEP(P), U(P), W(P),
     1                 WCHMTD(P), WSCALE, Y(P)
C     DIMENSION A(P*(P+1)/2)
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER I, J, K
      REAL DENMIN, SDOTWM, T, UI, WI
C
C     ***  CONSTANTS  ***
      REAL HALF, ONE, ZERO
C
C  ***  INTRINSIC FUNCTIONS  ***
C/+
      REAL  ABS, AMIN1
C/
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      EXTERNAL  SLVMUL
C
C/6
C     DATA HALF/0.5E+0/, ONE/1.E+0/, ZERO/0.E+0/
C/7
      PARAMETER (HALF=0.5E+0, ONE=1.E+0, ZERO=0.E+0)
C/
C
C-----------------------------------------------------------------------
C
      SDOTWM = SDOT(P,STEP,1,WCHMTD,1)
      DENMIN = COSMIN * SNRM2(P,STEP,1) * SNRM2(P,WCHMTD,1)
      WSCALE = ONE
      IF (DENMIN .NE. ZERO) WSCALE = AMIN1(ONE,  ABS(SDOTWM/DENMIN))
      T = ZERO
      IF (SDOTWM .NE. ZERO) T = WSCALE / SDOTWM
      DO 10 I = 1, P
 10      W(I) = T * WCHMTD(I)
      CALL SLVMUL(P, U, A, STEP)
      T = HALF * (SIZE * SDOT(P,STEP,1,U,1)  -  SDOT(P,STEP,1,Y,1))
      DO 20 I = 1, P
 20      U(I) = T*W(I) + Y(I) - SIZE*U(I)
C
C  ***  SET  A = A + U*(W**T) + W*(U**T)  ***
C
      K = 1
      DO 40 I = 1, P
         UI = U(I)
         WI = W(I)
         DO 30 J = 1, I
              A(K) = SIZE*A(K) + UI*W(J) + WI*U(J)
              K = K + 1
 30           CONTINUE
 40      CONTINUE
C
 999  RETURN
C  ***  LAST CARD OF SLUPDT FOLLOWS  ***
      END
      SUBROUTINE DUPD2(D, DR, IV, LIV, LV, N, NB, NN, N1, P, V)
C
C  ***  UPDATE SCALE VECTOR D FOR NL2IT  ***
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER LIV, LV, N, NB, NN, N1, P
      INTEGER IV(LIV)
      REAL D(P), DR(N,P,NN), V(LV)
C     DIMENSION V(*)
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER D0, I, JCN0, JCN1, JCNI, JTOL0, JTOLI, K, L, SII
      REAL T, VDFAC
C
C     ***  CONSTANTS  ***
C
      REAL ZERO
C
C  ***  INTRINSIC FUNCTIONS  ***
C/+
      REAL  ABS, AMAX1,  SQRT
C/
C  ***  EXTERNAL SUBROUTINE  ***
C
      EXTERNAL VSCOPY
C
C VSCOPY... SETS ALL COMPONENTS OF A VECTOR TO A SCALAR.
C
C  ***  SUBSCRIPTS FOR IV AND V  ***
C
      INTEGER DFAC, DTYPE, JCN, JTOL, NITER, S
C/6
C     DATA DFAC/41/, DTYPE/16/, JCN/66/, JTOL/59/, NITER/31/, S/62/
C/7
      PARAMETER (DFAC=41, DTYPE=16, JCN=66, JTOL=59, NITER=31, S=62)
C/
C
C/6
C     DATA ZERO/0.E+0/
C/7
      PARAMETER (ZERO=0.E+0)
C/
C
C-------------------------------  BODY  --------------------------------
C
      IF (IV(DTYPE) .NE. 1 .AND. IV(NITER) .GT. 0) GO TO 999
      JCN1 = IV(JCN)
      JCN0 = IABS(JCN1) - 1
      IF (JCN1 .LT. 0) GO TO 10
         IV(JCN) = -JCN1
         CALL VSCOPY(P, V(JCN1), ZERO)
 10   DO 40 L = 1, NN
         DO 30 I = 1, P
              JCNI = JCN0 + I
              T  = V(JCNI)
              DO 20 K = 1, N
 20                T = AMAX1(T,  ABS(DR(K,I,L)))
              V(JCNI) = T
 30           CONTINUE
 40      CONTINUE
      IF (N1 .LT. NB) GO TO 999
      VDFAC = V(DFAC)
      JTOL0 = IV(JTOL) - 1
      D0 = JTOL0 + P
      SII = IV(S) - 1
      DO 50 I = 1, P
         SII = SII + I
         JCNI = JCN0 + I
         T = V(JCNI)
         IF (V(SII) .GT. ZERO) T = AMAX1( SQRT(V(SII)), T)
         JTOLI = JTOL0 + I
         D0 = D0 + 1
         IF (T .LT. V(JTOLI)) T = AMAX1(V(D0), V(JTOLI))
         D(I) = AMAX1(VDFAC*D(I), T)
 50      CONTINUE
C
 999  RETURN
C  ***  LAST CARD OF DUPD2 FOLLOWS  ***
      END
      SUBROUTINE VAXPY(P, W, A, X, Y)
C
C  ***  SET W = A*X + Y  --  W, X, Y = P-VECTORS, A = SCALAR  ***
C
      INTEGER P
      REAL A, W(P), X(P), Y(P)
C
      INTEGER I
C
      DO 10 I = 1, P
 10      W(I) = A*X(I) + Y(I)
      RETURN
      END
      SUBROUTINE LITVMU(N, X, L, Y)
C
C  ***  SOLVE  (L**T)*X = Y,  WHERE  L  IS AN  N X N  LOWER TRIANGULAR
C  ***  MATRIX STORED COMPACTLY BY ROWS.  X AND Y MAY OCCUPY THE SAME
C  ***  STORAGE.  ***
C
C***REVISION HISTORY  (YYMMDD)
C   000330  Modified array declarations.  (JEC)
C
      INTEGER N
      REAL X(N), L(*), Y(N)
      INTEGER I, II, IJ, IM1, I0, J, NP1
      REAL XI, ZERO
C/6
      DATA ZERO/0.E+0/
C/7
C     PARAMETER (ZERO=0.E+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 LITVMU FOLLOWS  ***
      END
      SUBROUTINE FDHES(D, G, IRT, IV, LIV, LV, P, V, X)
C
C  ***  COMPUTE FINITE-DIFFERENCE HESSIAN, STORE IT IN V STARTING
C  ***  AT V(IV(FDH)) = V(-IV(H)).
C
C  ***  IF IV(COVREQ) .GE. 0 THEN FDHES USES GRADIENT DIFFERENCES,
C  ***  OTHERWISE FUNCTION DIFFERENCES.  STORAGE IN V IS AS IN GLRIT.
C
C IRT VALUES...
C     1 = COMPUTE FUNCTION VALUE, I.E., V(F).
C     2 = COMPUTE G.
C     3 = DONE.
C
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER IRT, LIV, LV, P
      INTEGER IV(LIV)
      REAL D(P), G(P), V(LV), X(P)
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER GSAVE1, HES, HMI, HPI, HPM, I, K, KIND, L, M, MM1, MM1O2,
     1        PP1O2, STPI, STPM, STP0
      REAL DEL, HALF, NEGPT5, ONE, TWO, ZERO
C
C  ***  INTRINSIC FUNCTIONS  ***
C/+
      INTEGER IABS
      REAL ABS, AMAX1
C/
C  ***  EXTERNAL SUBROUTINES  ***
C
C
C
C  ***  SUBSCRIPTS FOR IV AND V  ***
C
      INTEGER COVMAT, COVREQ, DELTA, DELTA0, DLTFDC, F, FDH, FX, H,
     1        KAGQT, MODE, NFGCAL, SAVEI, SWITCH, TOOBIG, W, XMSAVE
C
C/6
C     DATA HALF/0.5E+0/, NEGPT5/-0.5E+0/, ONE/1.E+0/, TWO/2.E+0/,
C    1     ZERO/0.E+0/
C/7
      PARAMETER (HALF=0.5E+0, NEGPT5=-0.5E+0, ONE=1.E+0, TWO=2.E+0,
     1     ZERO=0.E+0)
C/
C
C/6
C     DATA COVMAT/26/, COVREQ/15/, DELTA/52/, DELTA0/44/, DLTFDC/42/,
C    1     F/10/, FDH/74/, FX/53/, H/56/, KAGQT/33/, MODE/35/,
C    2     NFGCAL/7/, SAVEI/63/, SWITCH/12/, TOOBIG/2/, W/65/,
C    3     XMSAVE/51/
C/7
      PARAMETER (COVMAT=26, COVREQ=15, DELTA=52, DELTA0=44, DLTFDC=42,
     1           F=10, FDH=74, FX=53, H=56, KAGQT=33, MODE=35, NFGCAL=7,
     2           SAVEI=63, SWITCH=12, TOOBIG=2, W=65, XMSAVE=51)
C/
C
C+++++++++++++++++++++++++++++++  BODY  ++++++++++++++++++++++++++++++++
C
      IRT = 4
      KIND = IV(COVREQ)
      M = IV(MODE)
      IF (M .GT. 0) GO TO 10
         IV(H) = -IABS(IV(H))
         IV(FDH) = 0
         IV(KAGQT) = -1
         V(FX) = V(F)
 10   IF (M .GT. P) GO TO 999
      IF (KIND .LT. 0) GO TO 110
C
C  ***  COMPUTE FINITE-DIFFERENCE HESSIAN USING BOTH FUNCTION AND
C  ***  GRADIENT VALUES.
C
      GSAVE1 = IV(W) + P
      IF (M .GT. 0) GO TO 20
C        ***  FIRST CALL ON FDHES.  SET GSAVE = G, TAKE FIRST STEP  ***
         CALL SCOPY(P,V(GSAVE1),1,G,1)
         IV(SWITCH) = IV(NFGCAL)
         GO TO 90
C
 20   DEL = V(DELTA)
      X(M) = V(XMSAVE)
      IF (IV(TOOBIG) .EQ. 0) GO TO 40
C
C     ***  HANDLE OVERSIZE V(DELTA)  ***
C
         IF (DEL*X(M) .GT. ZERO) GO TO 30
C             ***  WE ALREADY TRIED SHRINKING V(DELTA), SO QUIT  ***
              IV(COVMAT) = -2
              GO TO 220
C
C        ***  TRY SHRINKING V(DELTA)  ***
 30      DEL = NEGPT5 * DEL
         GO TO 100
C
 40   HES = -IV(H)
C
C  ***  SET  G = (G - GSAVE)/DEL  ***
C
      DO 50 I = 1, P
         G(I) = (G(I) - V(GSAVE1)) / DEL
         GSAVE1 = GSAVE1 + 1
 50      CONTINUE
C
C  ***  ADD G AS NEW COL. TO FINITE-DIFF. HESSIAN MATRIX  ***
C
      K = HES + M*(M-1)/2
      L = K + M - 2
      IF (M .EQ. 1) GO TO 70
C
C  ***  SET  H(I,M) = 0.5 * (H(I,M) + G(I))  FOR I = 1 TO M-1  ***
C
      MM1 = M - 1
      DO 60 I = 1, MM1
         V(K) = HALF * (V(K) + G(I))
         K = K + 1
 60      CONTINUE
C
C  ***  ADD  H(I,M) = G(I)  FOR I = M TO P  ***
C
 70   L = L + 1
      DO 80 I = M, P
         V(L) = G(I)
         L = L + I
 80      CONTINUE
C
 90   M = M + 1
      IV(MODE) = M
      IF (M .GT. P) GO TO 210
C
C  ***  CHOOSE NEXT FINITE-DIFFERENCE STEP, RETURN TO GET G THERE  ***
C
      DEL = V(DELTA0) * AMAX1(ONE/D(M),  ABS(X(M)))
      IF (X(M) .LT. ZERO) DEL = -DEL
      V(XMSAVE) = X(M)
 100  X(M) = X(M) + DEL
      V(DELTA) = DEL
      IRT = 2
      GO TO 999
C
C  ***  COMPUTE FINITE-DIFFERENCE HESSIAN USING FUNCTION VALUES ONLY.
C
 110  STP0 = IV(W) + P - 1
      MM1 = M - 1
      MM1O2 = M*MM1/2
      IF (M .GT. 0) GO TO 120
C        ***  FIRST CALL ON FDHES.  ***
         IV(SAVEI) = 0
         GO TO 200
C
 120  I = IV(SAVEI)
      IF (I .GT. 0) GO TO 180
      IF (IV(TOOBIG) .EQ. 0) GO TO 140
C
C     ***  HANDLE OVERSIZE STEP  ***
C
         STPM = STP0 + M
         DEL = V(STPM)
         IF (DEL*X(XMSAVE) .GT. ZERO) GO TO 130
C             ***  WE ALREADY TRIED SHRINKING THE STEP, SO QUIT  ***
              IV(COVMAT) = -2
              GO TO 220
C
C        ***  TRY SHRINKING THE STEP  ***
 130     DEL = NEGPT5 * DEL
         X(M) = X(XMSAVE) + DEL
         V(STPM) = DEL
         IRT = 1
         GO TO 999
C
C  ***  SAVE F(X + STP(M)*E(M)) IN H(P,M)  ***
C
 140  PP1O2 = P * (P-1) / 2
      HES = -IV(H)
      HPM = HES + PP1O2 + MM1
      V(HPM) = V(F)
C
C  ***  START COMPUTING ROW M OF THE FINITE-DIFFERENCE HESSIAN H.  ***
C
      HMI = HES + MM1O2
      IF (MM1 .EQ. 0) GO TO 160
      HPI = HES + PP1O2
      DO 150 I = 1, MM1
         V(HMI) = V(FX) - (V(F) + V(HPI))
         HMI = HMI + 1
         HPI = HPI + 1
 150     CONTINUE
 160  V(HMI) = V(F) - TWO*V(FX)
C
C  ***  COMPUTE FUNCTION VALUES NEEDED TO COMPLETE ROW M OF H.  ***
C
      I = 1
C
 170  IV(SAVEI) = I
      STPI = STP0 + I
      V(DELTA) = X(I)
      X(I) = X(I) + V(STPI)
      IF (I .EQ. M) X(I) = V(XMSAVE) - V(STPI)
      IRT = 1
      GO TO 999
C
 180  X(I) = V(DELTA)
      IF (IV(TOOBIG) .EQ. 0) GO TO 190
C        ***  PUNT IN THE EVENT OF AN OVERSIZE STEP  ***
         IV(COVMAT) = -2
         GO TO 220
C
C  ***  FINISH COMPUTING H(M,I)  ***
C
 190  STPI = STP0 + I
      HMI = HES + MM1O2 + I - 1
      STPM = STP0 + M
      V(HMI) = (V(HMI) + V(F)) / (V(STPI)*V(STPM))
      I = I + 1
      IF (I .LE. M) GO TO 170
      IV(SAVEI) = 0
      X(M) = V(XMSAVE)
C
 200  M = M + 1
      IV(MODE) = M
      IF (M .GT. P) GO TO 210
C
C  ***  PREPARE TO COMPUTE ROW M OF THE FINITE-DIFFERENCE HESSIAN H.
C  ***  COMPUTE M-TH STEP SIZE STP(M), THEN RETURN TO OBTAIN
C  ***  F(X + STP(M)*E(M)), WHERE E(M) = M-TH STD. UNIT VECTOR.
C
      DEL = V(DLTFDC) * AMAX1(ONE/D(M),  ABS(X(M)))
      IF (X(M) .LT. ZERO) DEL = -DEL
      V(XMSAVE) = X(M)
      X(M) = X(M) + DEL
      STPM = STP0 + M
      V(STPM) = DEL
      IRT = 1
      GO TO 999
C
C  ***  RESTORE V(F), ETC.  ***
C
 210  IV(FDH) = HES
 220  V(F) = V(FX)
      IRT = 3
      IF (KIND .LT. 0) GO TO 999
         IV(NFGCAL) = IV(SWITCH)
         GSAVE1 = IV(W) + P
         CALL SCOPY(P,G,1,V(GSAVE1),1)
         GO TO 999
C
 999  RETURN
C  ***  LAST CARD OF FDHES FOLLOWS  ***
      END
      SUBROUTINE SLVMUL(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***REVISION HISTORY  (YYMMDD)
C   000330  Modified array declarations.  (JEC)
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER P
      REAL S(*), X(P), Y(P)
C     DIMENSION S(P*(P+1)/2)
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER I, IM1, J, K
      REAL XI
C
C  ***  NO INTRINSIC FUNCTIONS  ***
C
C  ***  EXTERNAL FUNCTION  ***
C
C
C-----------------------------------------------------------------------
C
      J = 1
      DO 10 I = 1, P
         Y(I) = SDOT(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 SLVMUL FOLLOWS  ***
      END
      REAL FUNCTION LSVMIN(P, L, X, Y)
C
C  ***  ESTIMATE SMALLEST SING. VALUE OF PACKED LOWER TRIANG. MATRIX L
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER P
      REAL L(*), X(P), Y(P)
C     DIMENSION L(P*(P+1)/2)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ***  PURPOSE  ***
C
C     THIS FUNCTION RETURNS A GOOD OVER-ESTIMATE OF THE SMALLEST
C     SINGULAR VALUE OF THE PACKED LOWER TRIANGULAR MATRIX L.
C
C  ***  PARAMETER DESCRIPTION  ***
C
C  P (IN)  = THE ORDER OF L.  L IS A  P X P  LOWER TRIANGULAR MATRIX.
C  L (IN)  = ARRAY HOLDING THE ELEMENTS OF  L  IN ROW ORDER, I.E.
C             L(1,1), L(2,1), L(2,2), L(3,1), L(3,2), L(3,3), ETC.
C  X (OUT) IF LSVMIN RETURNS A POSITIVE VALUE, THEN X IS A NORMALIZED
C             APPROXIMATE LEFT SINGULAR VECTOR CORRESPONDING TO THE
C             SMALLEST SINGULAR VALUE.  THIS APPROXIMATION MAY BE VERY
C             CRUDE.  IF LSVMIN RETURNS ZERO, THEN SOME COMPONENTS OF X
C             ARE ZERO AND THE REST RETAIN THEIR INPUT VALUES.
C  Y (OUT) IF LSVMIN RETURNS A POSITIVE VALUE, THEN Y = (L**-1)*X IS AN
C             UNNORMALIZED APPROXIMATE RIGHT SINGULAR VECTOR CORRESPOND-
C             ING TO THE SMALLEST SINGULAR VALUE.  THIS APPROXIMATION
C             MAY BE CRUDE.  IF LSVMIN RETURNS ZERO, THEN Y RETAINS ITS
C             INPUT VALUE.  THE CALLER MAY PASS THE SAME VECTOR FOR X
C             AND Y (NONSTANDARD FORTRAN USAGE), IN WHICH CASE Y OVER-
C             WRITES X (FOR NONZERO LSVMIN RETURNS).
C
C  ***  ALGORITHM NOTES  ***
C
C     THE ALGORITHM IS BASED ON (1), WITH THE ADDITIONAL PROVISION THAT
C     LSVMIN = 0 IS RETURNED IF THE SMALLEST DIAGONAL ELEMENT OF L
C     (IN MAGNITUDE) IS NOT MORE THAN THE UNIT ROUNDOFF TIMES THE
C     LARGEST.  THE ALGORITHM USES A RANDOM NUMBER GENERATOR PROPOSED
C     IN (4), WHICH PASSES THE SPECTRAL TEST WITH FLYING COLORS -- SEE
C     (2) AND (3).
C
C  ***  SUBROUTINES AND FUNCTIONS CALLED  ***
C
C
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***REVISION HISTORY  (YYMMDD)
C   000330  Modified array declarations.  (JEC)
C
C  ***  GENERAL  ***
C
C     THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH
C     SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS
C     MCS-7600324, DCR75-10143, 76-14311DSS, AND MCS76-11989.
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER I, II, IX, J, JI, JJ, JJJ, JM1, J0, PM1
      REAL B, SMINUS, SPLUS, T, XMINUS, XPLUS
C
C  ***  CONSTANTS  ***
C
      REAL HALF, ONE, R9973, ZERO
C
C  ***  INTRINSIC FUNCTIONS  ***
C/+
      INTEGER MOD
      REAL ABS, FLOAT
C/
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      EXTERNAL VAXPY
C
C/6
      DATA HALF/0.5E+0/, ONE/1.E+0/, R9973/9973.E+0/, ZERO/0.E+0/
C/7
C     PARAMETER (HALF=0.5E+0, ONE=1.E+0, R9973=9973.E+0, ZERO=0.E+0)
C/
C
C  ***  BODY  ***
C
      IX = 2
      PM1 = P - 1
C
C  ***  FIRST CHECK WHETHER TO RETURN LSVMIN = 0 AND INITIALIZE X  ***
C
      II = 0
      J0 = P*PM1/2
      JJ = J0 + P
      IF (L(JJ) .EQ. ZERO) GO TO 110
      IX = MOD(3432*IX, 9973)
      B = HALF*(ONE + FLOAT(IX)/R9973)
      XPLUS = B / L(JJ)
      X(P) = XPLUS
      IF (P .LE. 1) GO TO 60
      DO 10 I = 1, PM1
         II = II + I
         IF (L(II) .EQ. ZERO) GO TO 110
         JI = J0 + I
         X(I) = XPLUS * L(JI)
 10      CONTINUE
C
C  ***  SOLVE (L**T)*X = B, WHERE THE COMPONENTS OF B HAVE RANDOMLY
C  ***  CHOSEN MAGNITUDES IN (.5,1) WITH SIGNS CHOSEN TO MAKE X LARGE.
C
C     DO J = P-1 TO 1 BY -1...
      DO 50 JJJ = 1, PM1
         J = P - JJJ
C       ***  DETERMINE X(J) IN THIS ITERATION. NOTE FOR I = 1,2,...,J
C       ***  THAT X(I) HOLDS THE CURRENT PARTIAL SUM FOR ROW I.
         IX = MOD(3432*IX, 9973)
         B = HALF*(ONE + FLOAT(IX)/R9973)
         XPLUS = (B - X(J))
         XMINUS = (-B - X(J))
         SPLUS =  ABS(XPLUS)
         SMINUS =  ABS(XMINUS)
         JM1 = J - 1
         J0 = J*JM1/2
         JJ = J0 + J
         XPLUS = XPLUS/L(JJ)
         XMINUS = XMINUS/L(JJ)
         IF (JM1 .EQ. 0) GO TO 30
         DO 20 I = 1, JM1
              JI = J0 + I
              SPLUS = SPLUS +  ABS(X(I) + L(JI)*XPLUS)
              SMINUS = SMINUS +  ABS(X(I) + L(JI)*XMINUS)
 20           CONTINUE
 30      IF (SMINUS .GT. SPLUS) XPLUS = XMINUS
         X(J) = XPLUS
C       ***  UPDATE PARTIAL SUMS  ***
         IF (JM1 .GT. 0) CALL VAXPY(JM1, X, XPLUS, L(J0+1), X)
 50      CONTINUE
C
C  ***  NORMALIZE X  ***
C
 60   T = ONE/SNRM2(P, X,1)
      DO 70 I = 1, P
 70      X(I) = T*X(I)
C
C  ***  SOLVE L*Y = X AND RETURN LSVMIN = 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=SDOT(JM1,L(J0+1),1,Y,1)
      Y(J)=(X(J)-T)/L(JJ)
 100     CONTINUE
C
      LSVMIN = ONE/SNRM2(P, Y,1)
      GO TO 999
C
 110  LSVMIN = ZERO
 999  RETURN
C  ***  LAST CARD OF LSVMIN FOLLOWS  ***
      END
      SUBROUTINE ASSST(D, IV, P, STEP, STLSTG, V, X, X0)
C
C  ***  ASSESS CANDIDATE STEP (***SOL VERSION 2.3)  ***
C
      INTEGER P, IV(32)
      REAL 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  RELDST  AS
C                 MAX (D(I)*ABS(X(I)-X0(I)), 1 .LE. I .LE. P) /
C                    MAX (D(I)*(ABS(X(I))+ABS(X0(I))), 1 .LE. I .LE. P).
C             IF AN ACCEPTABLE STEP IS RETURNED, THEN V(RELDX) IS COM-
C             PUTED USING THE OUTPUT (POSSIBLY RESTORED) VALUES OF X
C             AND STEP.  OTHERWISE IT IS COMPUTED USING THE INPUT
C             VALUES.
C V(RFCTOL) (IN)  RELATIVE FUNCTION CONVERGENCE TOLERANCE.  IF THE
C             ACTUAL FUNCTION REDUCTION IS AT MOST TWICE WHAT WAS PRE-
C             DICTED AND  V(NREDUC) .LE. V(RFCTOL)*ABS(V(F0)),  THEN
C             ASSESS RETURNS WITH IV(IRC) = 8 OR 9.
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 RELDST
      REAL RELDST
C
C RELDST... COMPUTES V(RELDX) = RELATIVE STEP SIZE.
C
C  ***  INTRINSIC FUNCTIONS  ***
C/+
      REAL  ABS, AMAX1
C/
C  ***  NO COMMON BLOCKS  ***
C
C--------------------------  LOCAL VARIABLES  --------------------------
C
      LOGICAL GOODX
      INTEGER I, NFC
      REAL 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.5E+0/, ONE/1.E+0/, TWO/2.E+0/, ZERO/0.E+0/
C/7
C     PARAMETER (HALF=0.5E+0, ONE=1.E+0, TWO=2.E+0, ZERO=0.E+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 = RELDST(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 SCOPY(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 * AMAX1(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) = AMAX1(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 SCOPY(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 ( ABS(V(F)) .LT. V(AFCTOL)) IV(IRC) = 10
      IF (HALF * V(FDIF) .GT. V(PREDUC)) GO TO 999
      EMAX = V(RFCTOL) *  ABS(V(F0))
      EMAXS = V(SCTOL) *  ABS(V(F0))
      IF (V(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 SCOPY(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) =  ABS(V(DSTSAV))
      CALL SCOPY(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) *  ABS(V(F0))) IV(IRC) = 11
C
 999  RETURN
C
C  ***  LAST CARD OF ASSESS FOLLOWS  ***
      END
      REAL FUNCTION RELDST(P, D, X, X0)
C
C  ***  COMPUTE AND RETURN RELATIVE DIFFERENCE BETWEEN X AND X0  ***
C  ***  NL2SOL VERSION 2.2  ***
C
      INTEGER P
      REAL D(P), X(P), X0(P)
C/+
      REAL  ABS
C/
      INTEGER I
      REAL EMAX, T, XMAX, ZERO
C/6
      DATA ZERO/0.E+0/
C/7
C     PARAMETER (ZERO=0.E+0)
C/
C
      EMAX = ZERO
      XMAX = ZERO
      DO 10 I = 1, P
         T =  ABS(D(I) * (X(I) - X0(I)))
         IF (EMAX .LT. T) EMAX = T
         T = D(I) * ( ABS(X(I)) +  ABS(X0(I)))
         IF (XMAX .LT. T) XMAX = T
 10      CONTINUE
      RELDST = ZERO
      IF (XMAX .GT. ZERO) RELDST = EMAX / XMAX
 999  RETURN
C  ***  LAST CARD OF RELDST FOLLOWS  ***
      END
