Logo Search packages:      
Sourcecode: tseries version File versions

dsumsl.f

C
C  Added to each function and subroutine a "save" and replaced
C  "stop" with a call to the external function "error". "error" 
C  has to be provided by the caller, A. Trapletti, 14.10.1999
C
C  *** These routines are from the NIST Core Math LIBrary CML ***
C
      SUBROUTINE DSUMSL(N, D, X, CALCF, CALCG, IV, LIV, LV, V,
     1                  UIPARM, URPARM, UFPARM)
      save
C
C  ***  MINIMIZE GENERAL UNCONSTRAINED OBJECTIVE FUNCTION USING   ***
C  ***  ANALYTIC GRADIENT AND HESSIAN APPROX. FROM SECANT UPDATE  ***
C
      INTEGER N, LIV, LV
      INTEGER IV(LIV), UIPARM(1)
      DOUBLE PRECISION D(N), X(N), V(LV), URPARM(1)
C     DIMENSION V(71 + N*(N+15)/2), UIPARM(*), URPARM(*)
      EXTERNAL CALCF, CALCG, UFPARM
C
C  ***  PURPOSE  ***
C
C        THIS ROUTINE INTERACTS WITH SUBROUTINE  DSUMIT  IN AN ATTEMPT
C     TO FIND AN N-VECTOR  X*  THAT MINIMIZES THE (UNCONSTRAINED)
C     OBJECTIVE FUNCTION COMPUTED BY  CALCF.  (OFTEN THE  X*  FOUND IS
C     A LOCAL MINIMIZER RATHER THAN A GLOBAL ONE.)
C
C--------------------------  PARAMETER USAGE  --------------------------
C
C N........ (INPUT) THE NUMBER OF VARIABLES ON WHICH  F  DEPENDS, I.E.,
C                  THE NUMBER OF COMPONENTS IN  X.
C D........ (INPUT/OUTPUT) A SCALE VECTOR SUCH THAT  D(I)*X(I),
C                  I = 1,2,...,N,  ARE ALL IN COMPARABLE UNITS.
C                  D CAN STRONGLY AFFECT THE BEHAVIOR OF DSUMSL.
C                  FINDING THE BEST CHOICE OF D IS GENERALLY A TRIAL-
C                  AND-ERROR PROCESS.  CHOOSING D SO THAT D(I)*X(I)
C                  HAS ABOUT THE SAME VALUE FOR ALL I OFTEN WORKS WELL.
C                  THE DEFAULTS PROVIDED BY SUBROUTINE DDEFLT (SEE IV
C                  BELOW) REQUIRE THE CALLER TO SUPPLY D.
C X........ (INPUT/OUTPUT) BEFORE (INITIALLY) CALLING DSUMSL, THE CALL-
C                  ER SHOULD SET  X  TO AN INITIAL GUESS AT  X*.  WHEN
C                  DSUMSL RETURNS,  X  CONTAINS THE BEST POINT SO FAR
C                  FOUND, I.E., THE ONE THAT GIVES THE LEAST VALUE SO
C                  FAR SEEN FOR  F(X).
C CALCF.... (INPUT) A SUBROUTINE THAT, GIVEN X, COMPUTES F(X).  CALCF
C                  MUST BE DECLARED EXTERNAL IN THE CALLING PROGRAM.
C                  IT IS INVOKED BY
C                       CALL CALCF(N, X, NF, F, UIPARM, URPARM, UFPARM)
C                  NF IS THE INVOCATION COUNT FOR CALCF.  IT IS INCLUD-
C                  ED FOR POSSIBLE USE WITH CALCG.  IF X IS OUT OF
C                  BOUNDS (E.G., IF IT WOULD CAUSE OVERFLOW IN COMPUT-
C                  ING F(X)), THEN CALCF SHOULD SET NF TO 0.  THIS WILL
C                  CAUSE A SHORTER STEP TO BE ATTEMPTED.  THE OTHER
C                  PARAMETERS ARE AS DESCRIBED ABOVE AND BELOW.  CALCF
C                  SHOULD NOT CHANGE N, P, OR X.
C CALCG.... (INPUT) A SUBROUTINE THAT, GIVEN X, COMPUTES G(X), THE GRA-
C                  DIENT OF F AT X.  CALCG MUST BE DECLARED EXTERNAL IN
C                  THE CALLING PROGRAM.  IT IS INVOKED BY
C                       CALL CALCG(N, X, NF, G, UIPARM, URPARM, UFAPRM)
C                  NF IS THE INVOCATION COUNT FOR CALCF AT THE TIME
C                  F(X) WAS EVALUATED.  THE X PASSED TO CALCG IS
C                  USUALLY THE ONE PASSED TO CALCF ON EITHER ITS MOST
C                  RECENT INVOCATION OR THE ONE PRIOR TO IT.  IF CALCF
C                  SAVES INTERMEDIATE RESULTS FOR USE BY CALCG, THEN IT
C                  IS POSSIBLE TO TELL FROM NF WHETHER THEY ARE VALID
C                  FOR THE CURRENT X (OR WHICH COPY IS VALID IF TWO
C                  COPIES ARE KEPT).  IF G CANNOT BE COMPUTED AT X,
C                  THEN CALCG SHOULD SET NF TO 0.  IN THIS CASE, DSUMSL
C                  WILL RETURN WITH IV(1) = 65.  THE OTHER PARAMETERS
C                  TO CALCG ARE AS DESCRIBED ABOVE AND BELOW.  CALCG
C                  SHOULD NOT CHANGE N OR X.
C IV....... (INPUT/OUTPUT) AN INTEGER VALUE ARRAY OF LENGTH LIV (SEE
C                  BELOW) THAT HELPS CONTROL THE DSUMSL ALGORITHM AND
C                  THAT IS USED TO STORE VARIOUS INTERMEDIATE QUANTI-
C                  TIES.  OF PARTICULAR INTEREST ARE THE INITIALIZATION/
C                  RETURN CODE IV(1) AND THE ENTRIES IN IV THAT CONTROL
C                  PRINTING AND LIMIT THE NUMBER OF ITERATIONS AND FUNC-
C                  TION EVALUATIONS.  SEE THE SECTION ON IV INPUT
C                  VALUES BELOW.
C V........ (INPUT/OUTPUT) A FLOATING-POINT VALUE ARRAY OF LENGTH LV
C                  (SEE BELOW) THAT HELPS CONTROL THE DSUMSL ALGORITHM
C                  AND THAT IS USED TO STORE VARIOUS INTERMEDIATE
C                  QUANTITIES.  OF PARTICULAR INTEREST ARE THE ENTRIES
C                  IN V THAT LIMIT THE LENGTH OF THE FIRST STEP
C                  ATTEMPTED (LMAX0) AND SPECIFY CONVERGENCE TOLERANCES
C                  (AFCTOL, LMAXS, RFCTOL, SCTOL, XCTOL, XFTOL).
C LIV...... (INPUT) LENGTH OF IV ARRAY.  MUST BE AT LEAST 60.  IF LIV
C                  IS TOO SMALL, THEN DSUMSL RETURNS WITH IV(1) = 15.
C                  IF LIV IS AT LEAST LASTIV (= 44), THEN THE MINIMUM
C                  ACCEPTABLE VALUE OF LIV IS STORED IN IV(LASTIV)
C                  WHEN DSUMSL RETURNS.  (THIS IS INTENDED FOR USE
C                  WITH EXTENSIONS OF DSUMSL THAT HANDLES CONSTRAINTS.)
C LV....... (INPUT) LENGTH OF V ARRAY.  MUST BE AT LEAST 71+N*(N+15)/2.
C                  (AT LEAST 77+N*(N+17)/2 FOR DSMSNO, AT LEAST
C                  78+N*(N+12) FOR DHUMSL).  IF LV IS TOO SMALL, THEN
C                  DSUMSL RETURNS WITH IV(1) = 16.  IF LIV IS AT LEAST
C                  LASTV (= 45), THEN THE MINIMUM ACCEPTABLE VALUE OF
C                  LV IS STORED IN IV(LASTV) WHEN DSUMSL RETURNS.
C UIPARM... (INPUT) USER INTEGER PARAMETER ARRAY PASSED WITHOUT CHANGE
C                  TO CALCF AND CALCG.
C URPARM... (INPUT) USER FLOATING-POINT PARAMETER ARRAY PASSED WITHOUT
C                  CHANGE TO CALCF AND CALCG.
C UFPARM... (INPUT) USER EXTERNAL SUBROUTINE OR FUNCTION PASSED WITHOUT
C                  CHANGE TO CALCF AND CALCG.
C
C  ***  IV INPUT VALUES (FROM SUBROUTINE DDEFLT)  ***
C
C IV(1)...  ON INPUT, IV(1) SHOULD HAVE A VALUE BETWEEN 0 AND 14......
C             0 AND 12 MEAN THIS IS A FRESH START.  0 MEANS THAT
C                  DDEFLT(2, IV, LIV, LV, V)
C             IS TO BE CALLED TO PROVIDE ALL DEFAULT VALUES TO IV AND
C             V.  12 (THE VALUE THAT DDEFLT ASSIGNS TO IV(1)) MEANS THE
C             CALLER HAS ALREADY CALLED DDEFLT AND HAS POSSIBLY CHANGED
C             SOME IV AND/OR V ENTRIES TO NON-DEFAULT VALUES.
C             13 MEANS DDEFLT HAS BEEN CALLED AND THAT DSUMSL (AND
C             DSUMIT) SHOULD ONLY ALLOCATE STORAGE IN IV AND V.
C             14 MEANS THAT A STORAGE HAS BEEN ALLOCATED (E.G. BY A
C             CALL WITH IV(1) = 13) AND THAT THE ALGORITHM SHOULD BE
C             STARTED.  WHEN CALLED WITH IV(1) = 13, DSUMSL RETURNS
C             IV(1) = 14 UNLESS LIV OR LV IS TOO SMALL (OR N IS NOT
C             POSITIVE).  DEFAULT = 12.
C IV(INITH).... IV(25) TELLS WHETHER THE HESSIAN APPROXIMATION H SHOULD
C             BE INITIALIZED.  1 (THE DEFAULT) MEANS DSUMIT SHOULD
C             INITIALIZE H TO THE DIAGONAL MATRIX WHOSE I-TH DIAGONAL
C             ELEMENT IS D(I)**2.  0 MEANS THE CALLER HAS SUPPLIED A
C             CHOLESKY FACTOR  L  OF THE INITIAL HESSIAN APPROXIMATION
C             H = L*(L**T)  IN V, STARTING AT V(IV(LMAT)) = V(IV(42))
C             (AND STORED COMPACTLY BY ROWS).  NOTE THAT IV(LMAT) MAY
C             BE INITIALIZED BY CALLING DSUMSL WITH IV(1) = 13 (SEE
C             THE IV(1) DISCUSSION ABOVE).  DEFAULT = 1.
C IV(MXFCAL)... IV(17) GIVES THE MAXIMUM NUMBER OF FUNCTION EVALUATIONS
C             (CALLS ON CALCF) ALLOWED.  IF THIS NUMBER DOES NOT SUF-
C             FICE, THEN DSUMSL RETURNS WITH IV(1) = 9.  DEFAULT = 200.
C IV(MXITER)... IV(18) GIVES THE MAXIMUM NUMBER OF ITERATIONS ALLOWED.
C             IT ALSO INDIRECTLY LIMITS THE NUMBER OF GRADIENT EVALUA-
C             TIONS (CALLS ON CALCG) TO IV(MXITER) + 1.  IF IV(MXITER)
C             ITERATIONS DO NOT SUFFICE, THEN DSUMSL RETURNS WITH
C             IV(1) = 10.  DEFAULT = 150.
C IV(OUTLEV)... IV(19) CONTROLS THE NUMBER AND LENGTH OF ITERATION SUM-
C             MARY LINES PRINTED (BY DITSUM).  IV(OUTLEV) = 0 MEANS DO
C             NOT PRINT ANY SUMMARY LINES.  OTHERWISE, PRINT A SUMMARY
C             LINE AFTER EACH ABS(IV(OUTLEV)) ITERATIONS.  IF IV(OUTLEV)
C             IS POSITIVE, THEN SUMMARY LINES OF LENGTH 78 (PLUS CARRI-
C             AGE CONTROL) ARE PRINTED, INCLUDING THE FOLLOWING...  THE
C             ITERATION AND FUNCTION EVALUATION COUNTS, F = THE CURRENT
C             FUNCTION VALUE, RELATIVE DIFFERENCE IN FUNCTION VALUES
C             ACHIEVED BY THE LATEST STEP (I.E., RELDF = (F0-V(F))/F01,
C             WHERE F01 IS THE MAXIMUM OF ABS(V(F)) AND ABS(V(F0)) AND
C             V(F0) IS THE FUNCTION VALUE FROM THE PREVIOUS ITERA-
C             TION), THE RELATIVE FUNCTION REDUCTION PREDICTED FOR THE
C             STEP JUST TAKEN (I.E., PRELDF = V(PREDUC) / F01, WHERE
C             V(PREDUC) IS DESCRIBED BELOW), THE SCALED RELATIVE CHANGE
C             IN X (SEE V(RELDX) BELOW), THE STEP PARAMETER FOR THE
C             STEP JUST TAKEN (STPPAR = 0 MEANS A FULL NEWTON STEP,
C             BETWEEN 0 AND 1 MEANS A RELAXED NEWTON STEP, BETWEEN 1
C             AND 2 MEANS A DOUBLE DOGLEG STEP, GREATER THAN 2 MEANS
C             A SCALED DOWN CAUCHY STEP -- SEE SUBROUTINE DBLDOG), THE
C             2-NORM OF THE SCALE VECTOR D TIMES THE STEP JUST TAKEN
C             (SEE V(DSTNRM) BELOW), AND NPRELDF, I.E.,
C             V(NREDUC)/F01, WHERE V(NREDUC) IS DESCRIBED BELOW -- IF
C             NPRELDF IS POSITIVE, THEN IT IS THE RELATIVE FUNCTION
C             REDUCTION PREDICTED FOR A NEWTON STEP (ONE WITH
C             STPPAR = 0).  IF NPRELDF IS NEGATIVE, THEN IT IS THE
C             NEGATIVE OF THE RELATIVE FUNCTION REDUCTION PREDICTED
C             FOR A STEP COMPUTED WITH STEP BOUND V(LMAXS) FOR USE IN
C             TESTING FOR SINGULAR CONVERGENCE.
C                  IF IV(OUTLEV) IS NEGATIVE, THEN LINES OF LENGTH 50
C             ARE PRINTED, INCLUDING ONLY THE FIRST 6 ITEMS LISTED
C             ABOVE (THROUGH RELDX).
C             DEFAULT = 1.
C IV(PARPRT)... IV(20) = 1 MEANS PRINT ANY NONDEFAULT V VALUES ON A
C             FRESH START OR ANY CHANGED V VALUES ON A RESTART.
C             IV(PARPRT) = 0 MEANS SKIP THIS PRINTING.  DEFAULT = 1.
C IV(PRUNIT)... IV(21) IS THE OUTPUT UNIT NUMBER ON WHICH ALL PRINTING
C             IS DONE.  IV(PRUNIT) = 0 MEANS SUPPRESS ALL PRINTING.
C             DEFAULT = STANDARD OUTPUT UNIT (UNIT 6 ON MOST SYSTEMS).
C IV(SOLPRT)... IV(22) = 1 MEANS PRINT OUT THE VALUE OF X RETURNED (AS
C             WELL AS THE GRADIENT AND THE SCALE VECTOR D).
C             IV(SOLPRT) = 0 MEANS SKIP THIS PRINTING.  DEFAULT = 1.
C IV(STATPR)... IV(23) = 1 MEANS PRINT SUMMARY STATISTICS UPON RETURN-
C             ING.  THESE CONSIST OF THE FUNCTION VALUE, THE SCALED
C             RELATIVE CHANGE IN X CAUSED BY THE MOST RECENT STEP (SEE
C             V(RELDX) BELOW), THE NUMBER OF FUNCTION AND GRADIENT
C             EVALUATIONS (CALLS ON CALCF AND CALCG), AND THE RELATIVE
C             FUNCTION REDUCTIONS PREDICTED FOR THE LAST STEP TAKEN AND
C             FOR A NEWTON STEP (OR PERHAPS A STEP BOUNDED BY V(LMAX0)
C             -- SEE THE DESCRIPTIONS OF PRELDF AND NPRELDF UNDER
C             IV(OUTLEV) ABOVE).
C             IV(STATPR) = 0 MEANS SKIP THIS PRINTING.
C             IV(STATPR) = -1 MEANS SKIP THIS PRINTING AS WELL AS THAT
C             OF THE ONE-LINE TERMINATION REASON MESSAGE.  DEFAULT = 1.
C IV(X0PRT).... IV(24) = 1 MEANS PRINT THE INITIAL X AND SCALE VECTOR D
C             (ON A FRESH START ONLY).  IV(X0PRT) = 0 MEANS SKIP THIS
C             PRINTING.  DEFAULT = 1.
C
C  ***  (SELECTED) IV OUTPUT VALUES  ***
C
C IV(1)........ ON OUTPUT, IV(1) IS A RETURN CODE....
C             3 = X-CONVERGENCE.  THE SCALED RELATIVE DIFFERENCE (SEE
C                  V(RELDX)) BETWEEN THE CURRENT PARAMETER VECTOR X AND
C                  A LOCALLY OPTIMAL PARAMETER VECTOR IS VERY LIKELY AT
C                  MOST V(XCTOL).
C             4 = RELATIVE FUNCTION CONVERGENCE.  THE RELATIVE DIFFER-
C                  ENCE BETWEEN THE CURRENT FUNCTION VALUE AND ITS LO-
C                  CALLY OPTIMAL VALUE IS VERY LIKELY AT MOST V(RFCTOL).
C             5 = BOTH X- AND RELATIVE FUNCTION CONVERGENCE (I.E., THE
C                  CONDITIONS FOR IV(1) = 3 AND IV(1) = 4 BOTH HOLD).
C             6 = ABSOLUTE FUNCTION CONVERGENCE.  THE CURRENT FUNCTION
C                  VALUE IS AT MOST V(AFCTOL) IN ABSOLUTE VALUE.
C             7 = SINGULAR CONVERGENCE.  THE HESSIAN NEAR THE CURRENT
C                  ITERATE APPEARS TO BE SINGULAR OR NEARLY SO, AND A
C                  STEP OF LENGTH AT MOST V(LMAX0) IS UNLIKELY TO YIELD
C                  A RELATIVE FUNCTION DECREASE OF MORE THAN V(SCTOL).
C             8 = FALSE CONVERGENCE.  THE ITERATES APPEAR TO BE CONVERG-
C                  ING TO A NONCRITICAL POINT.  THIS MAY MEAN THAT THE
C                  CONVERGENCE TOLERANCES (V(AFCTOL), V(RFCTOL),
C                  V(XCTOL)) ARE TOO SMALL FOR THE ACCURACY TO WHICH
C                  THE FUNCTION AND GRADIENT ARE BEING COMPUTED, THAT
C                  THERE IS AN ERROR IN COMPUTING THE GRADIENT, OR THAT
C                  THE FUNCTION OR GRADIENT IS DISCONTINUOUS NEAR X.
C             9 = FUNCTION EVALUATION LIMIT REACHED WITHOUT OTHER CON-
C                  VERGENCE (SEE IV(MXFCAL)).
C            10 = ITERATION LIMIT REACHED WITHOUT OTHER CONVERGENCE
C                  (SEE IV(MXITER)).
C            11 = DSTOPX RETURNED .TRUE. (EXTERNAL INTERRUPT).  SEE THE
C                  USAGE NOTES BELOW.
C            14 = STORAGE HAS BEEN ALLOCATED (AFTER A CALL WITH
C                  IV(1) = 13).
C            17 = RESTART ATTEMPTED WITH N CHANGED.
C            18 = D HAS A NEGATIVE COMPONENT AND IV(DTYPE) .LE. 0.
C            19...43 = V(IV(1)) IS OUT OF RANGE.
C            63 = F(X) CANNOT BE COMPUTED AT THE INITIAL X.
C            64 = BAD PARAMETERS PASSED TO ASSESS (WHICH SHOULD NOT
C                  OCCUR).
C            65 = THE GRADIENT COULD NOT BE COMPUTED AT X (SEE CALCG
C                  ABOVE).
C            67 = BAD FIRST PARAMETER TO DDEFLT.
C            80 = IV(1) WAS OUT OF RANGE.
C            81 = N IS NOT POSITIVE.
C IV(G)........ IV(28) IS THE STARTING SUBSCRIPT IN V OF THE CURRENT
C             GRADIENT VECTOR (THE ONE CORRESPONDING TO X).
C IV(NFCALL)... IV(6) IS THE NUMBER OF CALLS SO FAR MADE ON CALCF (I.E.,
C             FUNCTION EVALUATIONS).
C IV(NGCALL)... IV(30) IS THE NUMBER OF GRADIENT EVALUATIONS (CALLS ON
C             CALCG).
C IV(NITER).... IV(31) IS THE NUMBER OF ITERATIONS PERFORMED.
C
C  ***  (SELECTED) V INPUT VALUES (FROM SUBROUTINE DDEFLT)  ***
C
C V(BIAS)..... V(43) IS THE BIAS PARAMETER USED IN SUBROUTINE DBLDOG --
C             SEE THAT SUBROUTINE FOR DETAILS.  DEFAULT = 0.8.
C V(AFCTOL)... V(31) IS THE ABSOLUTE FUNCTION CONVERGENCE TOLERANCE.
C             IF DSUMSL FINDS A POINT WHERE THE FUNCTION VALUE IS LESS
C             THAN V(AFCTOL) IN ABSOLUTE VALUE, AND IF DSUMSL DOES NOT
C             RETURN WITH IV(1) = 3, 4, OR 5, THEN IT RETURNS WITH
C             IV(1) = 6.  DEFAULT = MAX(10**-20, MACHEP**2), WHERE
C             MACHEP IS THE UNIT ROUNDOFF.
C V(DINIT).... V(38), IF NONNEGATIVE, IS THE VALUE TO WHICH THE SCALE
C             VECTOR D IS INITIALIZED.  DEFAULT = -1.
C V(LMAX0).... V(35) GIVES THE MAXIMUM 2-NORM ALLOWED FOR D TIMES THE
C             VERY FIRST STEP THAT DSUMSL ATTEMPTS.  THIS PARAMETER CAN
C             MARKEDLY AFFECT THE PERFORMANCE OF DSUMSL.
C V(LMAXS).... V(36) IS USED IN TESTING FOR SINGULAR CONVERGENCE -- IF
C             THE FUNCTION REDUCTION PREDICTED FOR A STEP OF LENGTH
C             BOUNDED BY V(LMAXS) IS AT MOST V(SCTOL) * ABS(F0), WHERE
C             F0  IS THE FUNCTION VALUE AT THE START OF THE CURRENT
C             ITERATION, AND IF DSUMSL DOES NOT RETURN WITH IV(1) = 3,
C             4, 5, OR 6, THEN IT RETURNS WITH IV(1) = 7.  DEFAULT = 1.
C V(RFCTOL)... V(32) IS THE RELATIVE FUNCTION CONVERGENCE TOLERANCE.
C             IF THE CURRENT MODEL PREDICTS A MAXIMUM POSSIBLE FUNCTION
C             REDUCTION (SEE V(NREDUC)) OF AT MOST V(RFCTOL)*ABS(F0)
C             AT THE START OF THE CURRENT ITERATION, WHERE  F0  IS THE
C             THEN CURRENT FUNCTION VALUE, AND IF THE LAST STEP ATTEMPT-
C             ED ACHIEVED NO MORE THAN TWICE THE PREDICTED FUNCTION
C             DECREASE, THEN DSUMSL RETURNS WITH IV(1) = 4 (OR 5).
C             DEFAULT = MAX(10**-10, MACHEP**(2/3)), WHERE MACHEP IS
C             THE UNIT ROUNDOFF.
C V(SCTOL).... V(37) IS THE SINGULAR CONVERGENCE TOLERANCE -- SEE THE
C             DESCRIPTION OF V(LMAXS) ABOVE.
C V(TUNER1)... V(26) HELPS DECIDE WHEN TO CHECK FOR FALSE CONVERGENCE.
C             THIS IS DONE IF THE ACTUAL FUNCTION DECREASE FROM THE
C             CURRENT STEP IS NO MORE THAN V(TUNER1) TIMES ITS PREDICT-
C             ED VALUE.  DEFAULT = 0.1.
C V(XCTOL).... V(33) IS THE X-CONVERGENCE TOLERANCE.  IF A NEWTON STEP
C             (SEE V(NREDUC)) IS TRIED THAT HAS V(RELDX) .LE. V(XCTOL)
C             AND IF THIS STEP YIELDS AT MOST TWICE THE PREDICTED FUNC-
C             TION DECREASE, THEN DSUMSL RETURNS WITH IV(1) = 3 (OR 5).
C             (SEE THE DESCRIPTION OF V(RELDX) BELOW.)
C             DEFAULT = MACHEP**0.5, WHERE MACHEP IS THE UNIT ROUNDOFF.
C V(XFTOL).... V(34) IS THE FALSE CONVERGENCE TOLERANCE.  IF A STEP IS
C             TRIED THAT GIVES NO MORE THAN V(TUNER1) TIMES THE PREDICT-
C             ED FUNCTION DECREASE AND THAT HAS V(RELDX) .LE. V(XFTOL),
C             AND IF DSUMSL DOES NOT RETURN WITH IV(1) = 3, 4, 5, 6, OR
C             7, THEN IT RETURNS WITH IV(1) = 8.  (SEE THE DESCRIPTION
C             OF V(RELDX) BELOW.)  DEFAULT = 100*MACHEP, WHERE
C             MACHEP IS THE UNIT ROUNDOFF.
C V(*)........ DDEFLT SUPPLIES TO V A NUMBER OF TUNING CONSTANTS, WITH
C             WHICH IT SHOULD ORDINARILY BE UNNECESSARY TO TINKER.  SEE
C             SECTION 17 OF VERSION 2.2 OF THE NL2SOL USAGE SUMMARY
C             (I.E., THE APPENDIX TO REF. 1) FOR DETAILS ON V(I),
C             I = DECFAC, INCFAC, PHMNFC, PHMXFC, RDFCMN, RDFCMX,
C             TUNER2, TUNER3, TUNER4, TUNER5.
C
C  ***  (SELECTED) V OUTPUT VALUES  ***
C
C V(DGNORM)... V(1) IS THE 2-NORM OF (DIAG(D)**-1)*G, WHERE G IS THE
C             MOST RECENTLY COMPUTED GRADIENT.
C V(DSTNRM)... V(2) IS THE 2-NORM OF DIAG(D)*STEP, WHERE STEP IS THE
C             CURRENT STEP.
C V(F)........ V(10) IS THE CURRENT FUNCTION VALUE.
C V(F0)....... V(13) IS THE FUNCTION VALUE AT THE START OF THE CURRENT
C             ITERATION.
C V(NREDUC)... V(6), IF POSITIVE, IS THE MAXIMUM FUNCTION REDUCTION
C             POSSIBLE ACCORDING TO THE CURRENT MODEL, I.E., THE FUNC-
C             TION REDUCTION PREDICTED FOR A NEWTON STEP (I.E.,
C             STEP = -H**-1 * G,  WHERE  G  IS THE CURRENT GRADIENT AND
C             H IS THE CURRENT HESSIAN APPROXIMATION).
C                  IF V(NREDUC) IS NEGATIVE, THEN IT IS THE NEGATIVE OF
C             THE FUNCTION REDUCTION PREDICTED FOR A STEP COMPUTED WITH
C             A STEP BOUND OF V(LMAXS) FOR USE IN TESTING FOR SINGULAR
C             CONVERGENCE.
C V(PREDUC)... V(7) IS THE FUNCTION REDUCTION PREDICTED (BY THE CURRENT
C             QUADRATIC MODEL) FOR THE CURRENT STEP.  THIS (DIVIDED BY
C             V(F0)) IS USED IN TESTING FOR RELATIVE FUNCTION
C             CONVERGENCE.
C V(RELDX).... V(17) IS THE SCALED RELATIVE CHANGE IN X CAUSED BY THE
C             CURRENT STEP, COMPUTED AS
C                  MAX(ABS(D(I)*(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             WHERE X = X0 + STEP.
C
C-------------------------------  NOTES  -------------------------------
C
C  ***  ALGORITHM NOTES  ***
C
C        THIS ROUTINE USES A HESSIAN APPROXIMATION COMPUTED FROM THE
C     BFGS UPDATE (SEE REF 3).  ONLY A CHOLESKY FACTOR OF THE HESSIAN
C     APPROXIMATION IS STORED, AND THIS IS UPDATED USING IDEAS FROM
C     REF. 4.  STEPS ARE COMPUTED BY THE DOUBLE DOGLEG SCHEME DESCRIBED
C     IN REF. 2.  THE STEPS ARE ASSESSED AS IN REF. 1.
C
C  ***  USAGE NOTES  ***
C
C        AFTER A RETURN WITH IV(1) .LE. 11, IT IS POSSIBLE TO RESTART,
C     I.E., TO CHANGE SOME OF THE IV AND V INPUT VALUES DESCRIBED ABOVE
C     AND CONTINUE THE ALGORITHM FROM THE POINT WHERE IT WAS INTERRUPT-
C     ED.  IV(1) SHOULD NOT BE CHANGED, NOR SHOULD ANY ENTRIES OF IV
C     AND V OTHER THAN THE INPUT VALUES (THOSE SUPPLIED BY DDEFLT).
C        THOSE WHO DO NOT WISH TO WRITE A CALCG WHICH COMPUTES THE
C     GRADIENT ANALYTICALLY SHOULD CALL DSMSNO RATHER THAN DSUMSL.
C     DSMSNO USES FINITE DIFFERENCES TO COMPUTE AN APPROXIMATE GRADIENT.
C        THOSE WHO WOULD PREFER TO PROVIDE F AND G (THE FUNCTION AND
C     GRADIENT) BY REVERSE COMMUNICATION RATHER THAN BY WRITING SUBROU-
C     TINES CALCF AND CALCG MAY CALL ON DSUMIT DIRECTLY.  SEE THE COM-
C     MENTS AT THE BEGINNING OF DSUMIT.
C        THOSE WHO USE DSUMSL INTERACTIVELY MAY WISH TO SUPPLY THEIR
C     OWN DSTOPX FUNCTION, WHICH SHOULD RETURN .TRUE. IF THE BREAK KEY
C     HAS BEEN PRESSED SINCE DSTOPX WAS LAST INVOKED.  THIS MAKES IT
C     POSSIBLE TO EXTERNALLY INTERRUPT DSUMSL (WHICH WILL RETURN WITH
C     IV(1) = 11 IF DSTOPX RETURNS .TRUE.).
C        STORAGE FOR G IS ALLOCATED AT THE END OF V.  THUS THE CALLER
C     MAY MAKE V LONGER THAN SPECIFIED ABOVE AND MAY ALLOW CALCG TO USE
C     ELEMENTS OF G BEYOND THE FIRST N AS SCRATCH STORAGE.
C
C  ***  PORTABILITY NOTES  ***
C
C        THE DSUMSL DISTRIBUTION TAPE CONTAINS BOTH SINGLE- AND DOUBLE-
C     PRECISION VERSIONS OF THE DSUMSL SOURCE CODE, SO IT SHOULD BE UN-
C     NECESSARY TO CHANGE PRECISIONS.
C        INTRINSIC FUNCTIONS ARE EXPLICITLY DECLARED.  ON CERTAIN COM-
C     PUTERS (E.G. UNIVAC), IT MAY BE NECESSARY TO COMMENT OUT THESE
C     DECLARATIONS.  SO THAT THIS MAY BE DONE AUTOMATICALLY BY A SIMPLE
C     PROGRAM, SUCH DECLARATIONS ARE PRECEDED BY A COMMENT HAVING C/+
C     IN COLUMNS 1-3 AND BLANKS IN COLUMNS 4-72 AND ARE FOLLOWED BY
C     A COMMENT HAVING C/ IN COLUMNS 1 AND 2 AND BLANKS IN COLUMNS 3-72.
C        THE DSUMSL SOURCE CODE IS EXPRESSED IN 1966 ANSI STANDARD
C     FORTRAN.  IT MAY BE CONVERTED TO FORTRAN 77 BY COMMENTING OUT ALL
C     LINES THAT FALL BETWEEN A LINE HAVING C/6 IN COLUMNS 1-3 AND A
C     LINE HAVING C/7 IN COLUMNS 1-3 AND BY REMOVING (I.E., REPLACING
C     BY A BLANK) THE C IN COLUMN 1 OF THE LINES THAT FOLLOW THE C/7
C     LINE AND PRECEDE A LINE HAVING C/ IN COLUMNS 1-2 AND BLANKS IN
C     COLUMNS 3-72.  THESE CHANGES CONVERT SOME DATA STATEMENTS INTO
C     PARAMETER STATEMENTS, CONVERT SOME VARIABLES FROM REAL TO
C     CHARACTER*4, AND MAKE THE DATA STATEMENTS THAT INITIALIZE THESE
C     VARIABLES USE CHARACTER STRINGS DELIMITED BY PRIMES INSTEAD
C     OF HOLLERITH CONSTANTS.  (SUCH VARIABLES AND DATA STATEMENTS
C     APPEAR ONLY IN MODULES DITSUM AND DPARCK.  PARAMETER STATEMENTS
C     APPEAR NEARLY EVERYWHERE.)
C
C  ***  REFERENCES  ***
C
C 1.  DENNIS, J.E., GAY, D.M., AND WELSCH, R.E. (1981), ALGORITHM 573 --
C             AN ADAPTIVE NONLINEAR LEAST-SQUARES ALGORITHM, ACM TRANS.
C             MATH. SOFTWARE 7, PP. 369-383.
C
C 2.  DENNIS, J.E., AND MEI, H.H.W. (1979), TWO NEW UNCONSTRAINED OPTI-
C             MIZATION ALGORITHMS WHICH USE FUNCTION AND GRADIENT
C             VALUES, J. OPTIM. THEORY APPLIC. 28, PP. 453-482.
C
C 3.  DENNIS, J.E., AND MORE, J.J. (1977), QUASI-NEWTON METHODS, MOTIVA-
C             TION AND THEORY, SIAM REV. 19, PP. 46-89.
C
C 4.  GOLDFARB, D. (1976), FACTORIZED VARIABLE METRIC METHODS FOR UNCON-
C             STRAINED OPTIMIZATION, MATH. COMPUT. 30, PP. 796-811.
C
C  ***  GENERAL  ***
C
C     CODED BY DAVID M. GAY (WINTER 1980).  REVISED SUMMER 1982.
C     THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH
C     SUPPORTED IN PART BY THE NATIONAL SCIENCE FOUNDATION UNDER
C     GRANTS MCS-7600324, DCR75-10143, 76-14311DSS, MCS76-11989,
C     AND MCS-7906671.
C
C.
C----------------------------  DECLARATIONS  ---------------------------
C
      EXTERNAL DDEFLT, DSUMIT
C
C DDEFLT.... SUPPLIES DEFAULT IV AND V INPUT COMPONENTS.
C DSUMIT... REVERSE-COMMUNICATION ROUTINE THAT CARRIES OUT DSUMSL ALGO-
C             RITHM.
C
      INTEGER G1, IV1, NF
      DOUBLE PRECISION F
C
C  ***  SUBSCRIPTS FOR IV   ***
C
      INTEGER NEXTV, NFCALL, NFGCAL, G, TOOBIG, VNEED
C
C/6
      DATA NEXTV/47/, NFCALL/6/, NFGCAL/7/, G/28/, TOOBIG/2/, VNEED/4/
C/7
C     PARAMETER (NEXTV=47, NFCALL=6, NFGCAL=7, G=28, TOOBIG=2, VNEED=4)
C/
C
C+++++++++++++++++++++++++++++++  BODY  ++++++++++++++++++++++++++++++++
C
      IF (IV(1) .EQ. 0) CALL DDEFLT(2, IV, LIV, LV, V)
      IV(VNEED) = IV(VNEED) + N
      IV1 = IV(1)
      IF (IV1 .EQ. 14) GO TO 10
      IF (IV1 .GT. 2 .AND. IV1 .LT. 12) GO TO 10
      G1 = 1
      IF (IV1 .EQ. 12) IV(1) = 13
      GO TO 20
C
 10   G1 = IV(G)
C
 20   CALL DSUMIT(D, F, V(G1), IV, LIV, LV, N, V, X)
      IF (IV(1) - 2) 30, 40, 50
C
 30   NF = IV(NFCALL)
      CALL CALCF(N, X, NF, F, UIPARM, URPARM, UFPARM)
      IF (NF .LE. 0) IV(TOOBIG) = 1
      GO TO 20
C
 40   CALL CALCG(N, X, IV(NFGCAL), V(G1), UIPARM, URPARM, UFPARM)
      GO TO 20
C
 50   IF (IV(1) .NE. 14) GO TO 999
C
C  ***  STORAGE ALLOCATION
C
      IV(G) = IV(NEXTV)
      IV(NEXTV) = IV(G) + N
      IF (IV1 .NE. 13) GO TO 10
C
 999  RETURN
C  ***  LAST CARD OF DSUMSL FOLLOWS  ***
      END
      SUBROUTINE DDEFLT(ALG, IV, LIV, LV, V)
      save
C
C  ***  SUPPLY ***SOL (VERSION 2.3) DEFAULT VALUES TO IV AND V  ***
C
C  ***  ALG = 1 MEANS REGRESSION CONSTANTS.
C  ***  ALG = 2 MEANS GENERAL UNCONSTRAINED OPTIMIZATION CONSTANTS.
C
      INTEGER LIV, LV
      INTEGER ALG, IV(LIV)
      DOUBLE PRECISION V(LV)
C
      EXTERNAL  DVDFLT
C DVDFLT.... PROVIDES DEFAULT VALUES TO V.
C
      INTEGER MIV, MV
      INTEGER MINIV(2), MINV(2)
C
C  ***  SUBSCRIPTS FOR IV  ***
C
      INTEGER ALGSAV, COVPRT, COVREQ, DTYPE, HC, IERR, INITH, INITS,
     1        IPIVOT, IVNEED, LASTIV, LASTV, LMAT, MXFCAL, MXITER,
     2        NFCOV, NGCOV, NVDFLT, OUTLEV, PARPRT, PARSAV, PERM,
     3        PRUNIT, QRTYP, RDREQ, RMAT, SOLPRT, STATPR, VNEED,
     4        VSAVE, X0PRT
C
C  ***  IV SUBSCRIPT VALUES  ***
C
C/6
      DATA ALGSAV/51/, COVPRT/14/, COVREQ/15/, DTYPE/16/, HC/71/,
     1     IERR/75/, INITH/25/, INITS/25/, IPIVOT/76/, IVNEED/3/,
     2     LASTIV/44/, LASTV/45/, LMAT/42/, MXFCAL/17/, MXITER/18/,
     3     NFCOV/52/, NGCOV/53/, NVDFLT/50/, OUTLEV/19/, PARPRT/20/,
     4     PARSAV/49/, PERM/58/, PRUNIT/21/, QRTYP/80/, RDREQ/57/,
     5     RMAT/78/, SOLPRT/22/, STATPR/23/, VNEED/4/, VSAVE/60/,
     6     X0PRT/24/
C/7
C     PARAMETER (ALGSAV=51, COVPRT=14, COVREQ=15, DTYPE=16, HC=71,
C    1           IERR=75, INITH=25, INITS=25, IPIVOT=76, IVNEED=3,
C    2           LASTIV=44, LASTV=45, LMAT=42, MXFCAL=17, MXITER=18,
C    3           NFCOV=52, NGCOV=53, NVDFLT=50, OUTLEV=19, PARPRT=20,
C    4           PARSAV=49, PERM=58, PRUNIT=21, QRTYP=80, RDREQ=57,
C    5           RMAT=78, SOLPRT=22, STATPR=23, VNEED=4, VSAVE=60,
C    6           X0PRT=24)
C/
      DATA MINIV(1)/80/, MINIV(2)/59/, MINV(1)/98/, MINV(2)/71/
C
C-------------------------------  BODY  --------------------------------
C
      IF (ALG .LT. 1 .OR. ALG .GT. 2) GO TO 40
      MIV = MINIV(ALG)
      IF (LIV .LT. MIV) GO TO 20
      MV = MINV(ALG)
      IF (LV .LT. MV) GO TO 30
      CALL DVDFLT(ALG, LV, V)
      IV(1) = 12
      IV(ALGSAV) = ALG
      IV(IVNEED) = 0
      IV(LASTIV) = MIV
      IV(LASTV) = MV
      IV(LMAT) = MV + 1
      IV(MXFCAL) = 200
      IV(MXITER) = 150
      IV(OUTLEV) = 1
      IV(PARPRT) = 1
      IV(PERM) = MIV + 1
      IV(PRUNIT) = I1MACH(2)
      IV(SOLPRT) = 1
      IV(STATPR) = 1
      IV(VNEED) = 0
      IV(X0PRT) = 1
C
      IF (ALG .GE. 2) GO TO 10
C
C  ***  REGRESSION  VALUES
C
      IV(COVPRT) = 3
      IV(COVREQ) = 1
      IV(DTYPE) = 1
      IV(HC) = 0
      IV(IERR) = 0
      IV(INITS) = 0
      IV(IPIVOT) = 0
      IV(NVDFLT) = 32
      IV(PARSAV) = 67
      IV(QRTYP) = 1
      IV(RDREQ) = 3
      IV(RMAT) = 0
      IV(VSAVE) = 58
      GO TO 999
C
C  ***  GENERAL OPTIMIZATION VALUES
C
 10   IV(DTYPE) = 0
      IV(INITH) = 1
      IV(NFCOV) = 0
      IV(NGCOV) = 0
      IV(NVDFLT) = 25
      IV(PARSAV) = 47
      GO TO 999
C
 20   IV(1) = 15
      GO TO 999
C
 30   IV(1) = 16
      GO TO 999
C
 40   IV(1) = 67
C
 999  RETURN
C  ***  LAST CARD OF DDEFLT FOLLOWS  ***
      END
      SUBROUTINE DSUMIT(D, FX, G, IV, LIV, LV, N, V, X)
      save
C
C  ***  CARRY OUT DSUMSL (UNCONSTRAINED MINIMIZATION) ITERATIONS, USING
C  ***  DOUBLE-DOGLEG/BFGS STEPS.
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER LIV, LV, N
      INTEGER IV(LIV)
      DOUBLE PRECISION D(N), FX, G(N), V(LV), X(N)
C
C--------------------------  PARAMETER USAGE  --------------------------
C
C D.... SCALE VECTOR.
C FX... FUNCTION VALUE.
C G.... GRADIENT VECTOR.
C IV... INTEGER VALUE ARRAY.
C LIV.. LENGTH OF IV (AT LEAST 60).
C LV... LENGTH OF V (AT LEAST 71 + N*(N+13)/2).
C N.... NUMBER OF VARIABLES (COMPONENTS IN X AND G).
C V.... FLOATING-POINT VALUE ARRAY.
C X.... VECTOR OF PARAMETERS TO BE OPTIMIZED.
C
C  ***  DISCUSSION  ***
C
C        PARAMETERS IV, N, V, AND X ARE THE SAME AS THE CORRESPONDING
C     ONES TO DSUMSL (WHICH SEE), EXCEPT THAT V CAN BE SHORTER (SINCE
C     THE PART OF V THAT DSUMSL USES FOR STORING G IS NOT NEEDED).
C     MOREOVER, COMPARED WITH DSUMSL, IV(1) MAY HAVE THE TWO ADDITIONAL
C     OUTPUT VALUES 1 AND 2, WHICH ARE EXPLAINED BELOW, AS IS THE USE
C     OF IV(TOOBIG) AND IV(NFGCAL).  THE VALUE IV(G), WHICH IS AN
C     OUTPUT VALUE FROM DSUMSL (AND DSMSNO), IS NOT REFERENCED BY
C     DSUMIT OR THE SUBROUTINES IT CALLS.
C        FX AND G NEED NOT HAVE BEEN INITIALIZED WHEN DSUMIT IS CALLED
C     WITH IV(1) = 12, 13, OR 14.
C
C IV(1) = 1 MEANS THE CALLER SHOULD SET FX TO F(X), THE FUNCTION VALUE
C             AT X, AND CALL DSUMIT AGAIN, HAVING CHANGED NONE OF THE
C             OTHER PARAMETERS.  AN EXCEPTION OCCURS IF F(X) CANNOT BE
C             (E.G. IF OVERFLOW WOULD OCCUR), WHICH MAY HAPPEN BECAUSE
C             OF AN OVERSIZED STEP.  IN THIS CASE THE CALLER SHOULD SET
C             IV(TOOBIG) = IV(2) TO 1, WHICH WILL CAUSE DSUMIT TO IG-
C             NORE FX AND TRY A SMALLER STEP.  THE PARAMETER NF THAT
C             DSUMSL PASSES TO CALCF (FOR POSSIBLE USE BY CALCG) IS A
C             COPY OF IV(NFCALL) = IV(6).
C IV(1) = 2 MEANS THE CALLER SHOULD SET G TO G(X), THE GRADIENT VECTOR
C             OF F AT X, AND CALL DSUMIT AGAIN, HAVING CHANGED NONE OF
C             THE OTHER PARAMETERS EXCEPT POSSIBLY THE SCALE VECTOR D
C             WHEN IV(DTYPE) = 0.  THE PARAMETER NF THAT DSUMSL PASSES
C             TO CALCG IS IV(NFGCAL) = IV(7).  IF G(X) CANNOT BE
C             EVALUATED, THEN THE CALLER MAY SET IV(NFGCAL) TO 0, IN
C             WHICH CASE DSUMIT WILL RETURN WITH IV(1) = 65.
C.
C  ***  GENERAL  ***
C
C     CODED BY DAVID M. GAY (DECEMBER 1979).  REVISED SEPT. 1982.
C     THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH SUPPORTED
C     IN PART BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS
C     MCS-7600324 AND MCS-7906671.
C
C        (SEE DSUMSL FOR REFERENCES.)
C
C+++++++++++++++++++++++++++  DECLARATIONS  ++++++++++++++++++++++++++++
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER DG1, DUMMY, G01, I, K, L, LSTGST, NN1O2, NWTST1, STEP1,
     1        TEMP1, W, X01, Z
      DOUBLE PRECISION T
C
C     ***  CONSTANTS  ***
C
      DOUBLE PRECISION NEGONE, ONE, ZERO
C
C  ***  NO INTRINSIC FUNCTIONS  ***
C
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      EXTERNAL DASSST, DDBDOG, DDEFLT, DITSUM, DLITVM, DLIVMU,
     1         DLTVMU, DLUPDT, DLVMUL, DPARCK, DSTOPX, DVAXPY,
     2         DVSCPY, DVVMUP, DWZBFG
      LOGICAL DSTOPX
      DOUBLE PRECISION DDOT, DNRM2
C
C DASSST.... ASSESSES CANDIDATE STEP.
C DDBDOG.... COMPUTES DOUBLE-DOGLEG (CANDIDATE) STEP.
C DDEFLT.... SUPPLIES DEFAULT IV AND V INPUT COMPONENTS.
C DITSUM.... PRINTS ITERATION SUMMARY AND INFO ON INITIAL AND FINAL X.
C DLITVM... MULTIPLIES INVERSE TRANSPOSE OF LOWER TRIANGLE TIMES VECTOR.
C DLIVMU... MULTIPLIES INVERSE OF LOWER TRIANGLE TIMES VECTOR.
C DLTVMU... MULTIPLIES TRANSPOSE OF LOWER TRIANGLE TIMES VECTOR.
C LUPDT.... UPDATES CHOLESKY FACTOR OF HESSIAN APPROXIMATION.
C DLVMUL.... MULTIPLIES LOWER TRIANGLE TIMES VECTOR.
C DPARCK.... CHECKS VALIDITY OF INPUT IV AND V VALUES.
C DSTOPX.... RETURNS .TRUE. IF THE BREAK KEY HAS BEEN PRESSED.
C DVAXPY.... COMPUTES SCALAR TIMES ONE VECTOR PLUS ANOTHER.
C DVSCPY... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR.
C DVVMUP... MULTIPLIES VECTOR BY VECTOR RAISED TO POWER (COMPONENTWISE).
C DWZBFG... COMPUTES W AND Z FOR DLUPDT CORRESPONDING TO BFGS UPDATE.
C
C  ***  SUBSCRIPTS FOR IV AND V  ***
C
      INTEGER CNVCOD, DG, DGNORM, DINIT, DSTNRM, DST0, F, F0,
     1        GTHG, GTSTEP, G0, INCFAC, INITH, IRC, KAGQT, LMAT,
     2        LMAX0, MODE, MODEL, MXFCAL, MXITER, NEXTV, NFCALL, NFGCAL,
     3        NGCALL, NITER, NWTSTP, RADFAC, RADINC, RADIUS, RAD0, STEP,
     4        STGLIM, STLSTG, TOOBIG, TUNER4, TUNER5, VNEED, XIRC, X0
C
C  ***  IV SUBSCRIPT VALUES  ***
C
C/6
      DATA CNVCOD/55/, DG/37/, G0/48/, INITH/25/, IRC/29/, KAGQT/33/,
     1     MODE/35/, MODEL/5/, MXFCAL/17/, MXITER/18/, NFCALL/6/,
     2     NFGCAL/7/, NGCALL/30/, NITER/31/, NWTSTP/34/, RADINC/8/,
     3     STEP/40/, STGLIM/11/, STLSTG/41/, TOOBIG/2/, XIRC/13/, X0/43/
C/7
C     PARAMETER (CNVCOD=55, DG=37, G0=48, INITH=25, IRC=29, KAGQT=33,
C    1           MODE=35, MODEL=5, MXFCAL=17, MXITER=18, NFCALL=6,
C    2           NFGCAL=7, NGCALL=30, NITER=31, NWTSTP=34, RADINC=8,
C    3           STEP=40, STGLIM=11, STLSTG=41, TOOBIG=2, XIRC=13,
C    4           X0=43)
C/
C
C  ***  V SUBSCRIPT VALUES  ***
C
C/6
      DATA DGNORM/1/, DINIT/38/, DSTNRM/2/, DST0/3/, F/10/, F0/13/,
     1     GTHG/44/, GTSTEP/4/, INCFAC/23/, LMAT/42/, LMAX0/35/,
     2     NEXTV/47/, RADFAC/16/, RADIUS/8/, RAD0/9/, TUNER4/29/,
     3     TUNER5/30/, VNEED/4/
C/7
C     PARAMETER (DGNORM=1, DINIT=38, DSTNRM=2, DST0=3, F=10, F0=13,
C    1           GTHG=44, GTSTEP=4, INCFAC=23, LMAT=42, LMAX0=35,
C    2           NEXTV=47, RADFAC=16, RADIUS=8, RAD0=9, TUNER4=29,
C    3           TUNER5=30, VNEED=4)
C/
C
C/6
      DATA NEGONE/-1.D+0/, ONE/1.D+0/, ZERO/0.D+0/
C/7
C     PARAMETER (NEGONE=-1.D+0, ONE=1.D+0, ZERO=0.D+0)
C/
C
C+++++++++++++++++++++++++++++++  BODY  ++++++++++++++++++++++++++++++++
C
      I = IV(1)
      IF (I .EQ. 1) GO TO 40
      IF (I .EQ. 2) GO TO 50
C
C  ***  CHECK VALIDITY OF IV AND V INPUT VALUES  ***
C
      IF (IV(1) .EQ. 0) CALL DDEFLT(2, IV, LIV, LV, V)
      IV(VNEED) = IV(VNEED) + N*(N+13)/2
      CALL DPARCK(2, D, IV, LIV, LV, N, V)
      I = IV(1) - 2
      IF (I .GT. 12) GO TO 999
      GO TO (160, 160, 160, 160, 160, 160, 110, 80, 110, 10, 10, 20), I
C
C  ***  STORAGE ALLOCATION  ***
C
 10   NN1O2 = N * (N + 1) / 2
      L = IV(LMAT)
      IV(X0) = L + NN1O2
      IV(STEP) = IV(X0) + N
      IV(STLSTG) = IV(STEP) + N
      IV(G0) = IV(STLSTG) + N
      IV(NWTSTP) = IV(G0) + N
      IV(DG) = IV(NWTSTP) + N
      IV(NEXTV) = IV(DG) + N
      IF (IV(1) .NE. 13) GO TO 20
         IV(1) = 14
         GO TO 999
C
C  ***  INITIALIZATION  ***
C
 20   IV(NITER) = 0
      IV(NFCALL) = 1
      IV(NGCALL) = 1
      IV(NFGCAL) = 1
      IV(MODE) = -1
      IV(MODEL) = 1
      IV(STGLIM) = 1
      IV(TOOBIG) = 0
      IV(CNVCOD) = 0
      IV(RADINC) = 0
      V(RAD0) = ZERO
      IF (V(DINIT) .GE. ZERO) CALL DVSCPY(N, D, V(DINIT))
      IV(1) = 1
      IF (IV(INITH) .NE. 1) GO TO 999
C
C     ***  SET THE INITIAL HESSIAN APPROXIMATION TO DIAG(D)**-2  ***
C
         CALL DVSCPY(NN1O2, V(L), ZERO)
         K = L - 1
         DO 30 I = 1, N
              K = K + I
              T = D(I)
              IF (T .LE. ZERO) T = ONE
              V(K) = T
 30           CONTINUE
      GO TO 999
C
 40   V(F) = FX
      IF (IV(MODE) .GE. 0) GO TO 160
      IV(1) = 2
      IF (IV(TOOBIG) .EQ. 0) GO TO 999
         IV(1) = 63
         GO TO 270
C
C  ***  MAKE SURE GRADIENT COULD BE COMPUTED  ***
C
 50   IF (IV(NFGCAL) .NE. 0) GO TO 60
         IV(1) = 65
         GO TO 270
C
 60   DG1 = IV(DG)
      CALL DVVMUP(N, V(DG1), G, D, -1)
      V(DGNORM) = DNRM2(N, V(DG1),1)
C
      IF (IV(CNVCOD) .NE. 0) GO TO 260
      IF (IV(MODE) .EQ. 0) GO TO 220
C
C  ***  ALLOW FIRST STEP TO HAVE SCALED 2-NORM AT MOST V(LMAX0)  ***
C
      V(RADFAC) = V(LMAX0)
      V(DSTNRM) = ONE
C
      IV(MODE) = 0
C
C
C-----------------------------  MAIN LOOP  -----------------------------
C
C
C  ***  PRINT ITERATION SUMMARY, CHECK ITERATION LIMIT  ***
C
 70   CALL DITSUM(D, G, IV, LIV, LV, N, V, X)
 80   K = IV(NITER)
      IF (K .LT. IV(MXITER)) GO TO 90
         IV(1) = 10
         GO TO 270
C
C  ***  UPDATE RADIUS  ***
C
 90   IV(NITER) = K + 1
      V(RADIUS) = V(RADFAC) * V(DSTNRM)
C
C  ***  INITIALIZE FOR START OF NEXT ITERATION  ***
C
      G01 = IV(G0)
      X01 = IV(X0)
      V(F0) = V(F)
      IV(IRC) = 4
      IV(KAGQT) = -1
C
C     ***  COPY X TO X0, G TO G0  ***
C
      CALL DCOPY(N, X,1,V(X01),1)
      CALL DCOPY(N, G,1,V(G01),1)
C
C  ***  CHECK DSTOPX AND FUNCTION EVALUATION LIMIT  ***
C
 100  IF (.NOT. DSTOPX(DUMMY)) GO TO 120
         IV(1) = 11
         GO TO 130
C
C     ***  COME HERE WHEN RESTARTING AFTER FUNC. EVAL. LIMIT OR DSTOPX.
C
 110  IF (V(F) .GE. V(F0)) GO TO 120
         V(RADFAC) = ONE
         K = IV(NITER)
         GO TO 90
C
 120  IF (IV(NFCALL) .LT. IV(MXFCAL)) GO TO 140
         IV(1) = 9
 130     IF (V(F) .GE. V(F0)) GO TO 270
C
C        ***  IN CASE OF DSTOPX OR FUNCTION EVALUATION LIMIT WITH
C        ***  IMPROVED V(F), EVALUATE THE GRADIENT AT X.
C
              IV(CNVCOD) = IV(1)
              GO TO 210
C
C. . . . . . . . . . . . .  COMPUTE CANDIDATE STEP  . . . . . . . . . .
C
 140  STEP1 = IV(STEP)
      DG1 = IV(DG)
      NWTST1 = IV(NWTSTP)
      IF (IV(KAGQT) .GE. 0) GO TO 150
         L = IV(LMAT)
         CALL DLIVMU(N, V(NWTST1), V(L), G)
         CALL DLITVM(N, V(NWTST1), V(L), V(NWTST1))
         CALL DVVMUP(N, V(STEP1), V(NWTST1), D, 1)
         V(DST0) = DNRM2(N, V(STEP1),1)
         CALL DVVMUP(N, V(DG1), V(DG1), D, -1)
         CALL DLTVMU(N, V(STEP1), V(L), V(DG1))
         V(GTHG) = DNRM2(N, V(STEP1),1)
         IV(KAGQT) = 0
 150  CALL DDBDOG(V(DG1), G, LV, N, V(NWTST1), V(STEP1), V)
      IF (IV(IRC) .EQ. 6) GO TO 160
C
C  ***  COMPUTE F(X0 + STEP)  ***
C
      X01 = IV(X0)
      STEP1 = IV(STEP)
      CALL DVAXPY(N, X, ONE, V(STEP1), V(X01))
      IV(NFCALL) = IV(NFCALL) + 1
      IV(1) = 1
      IV(TOOBIG) = 0
      GO TO 999
C
C. . . . . . . . . . . . .  ASSESS CANDIDATE STEP  . . . . . . . . . . .
C
 160  STEP1 = IV(STEP)
      LSTGST = IV(STLSTG)
      X01 = IV(X0)
      CALL DASSST(D, IV, N, V(STEP1), V(LSTGST), V, X, V(X01))
C
      K = IV(IRC)
      GO TO (170,200,200,200,170,180,190,190,190,190,190,190,250,220), K
C
C     ***  RECOMPUTE STEP WITH CHANGED RADIUS  ***
C
 170     V(RADIUS) = V(RADFAC) * V(DSTNRM)
         GO TO 100
C
C  ***  COMPUTE STEP OF LENGTH V(LMAX0) FOR SINGULAR CONVERGENCE TEST.
C
 180  V(RADIUS) = V(LMAX0)
      GO TO 140
C
C  ***  CONVERGENCE OR FALSE CONVERGENCE  ***
C
 190  IV(CNVCOD) = K - 4
      IF (V(F) .GE. V(F0)) GO TO 260
         IF (IV(XIRC) .EQ. 14) GO TO 260
              IV(XIRC) = 14
C
C. . . . . . . . . . . .  PROCESS ACCEPTABLE STEP  . . . . . . . . . . .
C
 200  IF (IV(IRC) .NE. 3) GO TO 210
         STEP1 = IV(STEP)
         TEMP1 = IV(STLSTG)
C
C     ***  SET  TEMP1 = HESSIAN * STEP  FOR USE IN GRADIENT TESTS  ***
C
         L = IV(LMAT)
         CALL DLTVMU(N, V(TEMP1), V(L), V(STEP1))
         CALL DLVMUL(N, V(TEMP1), V(L), V(TEMP1))
C
C  ***  COMPUTE GRADIENT  ***
C
 210  IV(NGCALL) = IV(NGCALL) + 1
      IV(1) = 2
      GO TO 999
C
C  ***  INITIALIZATIONS -- G0 = G - G0, ETC.  ***
C
 220  G01 = IV(G0)
      CALL DVAXPY(N, V(G01), NEGONE, V(G01), G)
      STEP1 = IV(STEP)
      TEMP1 = IV(STLSTG)
      IF (IV(IRC) .NE. 3) GO TO 240
C
C  ***  SET V(RADFAC) BY GRADIENT TESTS  ***
C
C     ***  SET  TEMP1 = DIAG(D)**-1 * (HESSIAN*STEP + (G(X0)-G(X)))  ***
C
         CALL DVAXPY(N, V(TEMP1), NEGONE, V(G01), V(TEMP1))
         CALL DVVMUP(N, V(TEMP1), V(TEMP1), D, -1)
C
C        ***  DO GRADIENT TESTS  ***
C
         IF (DNRM2(N, V(TEMP1),1) .LE. V(DGNORM) * V(TUNER4))
     1                  GO TO 230
              IF (DDOT(N, G,1,V(STEP1),1)
     1                  .GE. V(GTSTEP) * V(TUNER5))  GO TO 240
 230               V(RADFAC) = V(INCFAC)
C
C  ***  UPDATE H, LOOP  ***
C
 240  W = IV(NWTSTP)
      Z = IV(X0)
      L = IV(LMAT)
      CALL DWZBFG(V(L), N, V(STEP1), V(W), V(G01), V(Z))
C
C     ** USE THE N-VECTORS STARTING AT V(STEP1) AND V(G01) FOR SCRATCH..
      CALL DLUPDT(V(TEMP1), V(STEP1), V(L), V(G01), V(L), N, V(W), V(Z))
      IV(1) = 2
      GO TO 70
C
C. . . . . . . . . . . . . .  MISC. DETAILS  . . . . . . . . . . . . . .
C
C  ***  BAD PARAMETERS TO ASSESS  ***
C
 250  IV(1) = 64
C
C  ***  PRINT SUMMARY OF FINAL ITERATION AND OTHER REQUESTED ITEMS  ***
C
 260  IV(1) = IV(CNVCOD)
      IV(CNVCOD) = 0
 270  CALL DITSUM(D, G, IV, LIV, LV, N, V, X)
C
 999  RETURN
C
C  ***  LAST CARD OF DSUMIT FOLLOWS  ***
      END
      SUBROUTINE DVAXPY(P, W, A, X, Y)
      save
C
C  ***  SET W = A*X + Y  --  W, X, Y = P-VECTORS, A = SCALAR  ***
C
      INTEGER P
      DOUBLE PRECISION A, W(P), X(P), Y(P)
C
      INTEGER I
C
      DO 10 I = 1, P
 10      W(I) = A*X(I) + Y(I)
      RETURN
      END
      SUBROUTINE DVDFLT(ALG, LV, V)
      save
C
C  ***  SUPPLY ***SOL (VERSION 2.3) DEFAULT VALUES TO V  ***
C
C  ***  ALG = 1 MEANS REGRESSION CONSTANTS.
C  ***  ALG = 2 MEANS GENERAL UNCONSTRAINED OPTIMIZATION CONSTANTS.
C
      INTEGER ALG, LV
      DOUBLE PRECISION V(LV)
C/+
      DOUBLE PRECISION DMAX1
C/
      DOUBLE PRECISION D1MACH
C
      DOUBLE PRECISION MACHEP, MEPCRT, ONE, SQTEPS, THREE
C
C  ***  SUBSCRIPTS FOR V  ***
C
      INTEGER AFCTOL, BIAS, COSMIN, DECFAC, DELTA0, DFAC, DINIT, DLTFDC,
     1        DLTFDJ, DTINIT, D0INIT, EPSLON, ETA0, FUZZ, HUBERC,
     2        INCFAC, LMAX0, LMAXS, PHMNFC, PHMXFC, RDFCMN, RDFCMX,
     3        RFCTOL, RLIMIT, RSPTOL, SCTOL, SIGMIN, TUNER1, TUNER2,
     4        TUNER3, TUNER4, TUNER5, XCTOL, XFTOL
C
C/6
      DATA ONE/1.D+0/, THREE/3.D+0/
C/7
C     PARAMETER (ONE=1.D+0, THREE=3.D+0)
C/
C
C  ***  V SUBSCRIPT VALUES  ***
C
C/6
      DATA AFCTOL/31/, BIAS/43/, COSMIN/47/, DECFAC/22/, DELTA0/44/,
     1     DFAC/41/, DINIT/38/, DLTFDC/42/, DLTFDJ/43/, DTINIT/39/,
     2     D0INIT/40/, EPSLON/19/, ETA0/42/, FUZZ/45/, HUBERC/48/,
     3     INCFAC/23/, LMAX0/35/, LMAXS/36/, PHMNFC/20/, PHMXFC/21/,
     4     RDFCMN/24/, RDFCMX/25/, RFCTOL/32/, RLIMIT/46/, RSPTOL/49/,
     5     SCTOL/37/, SIGMIN/50/, TUNER1/26/, TUNER2/27/, TUNER3/28/,
     6     TUNER4/29/, TUNER5/30/, XCTOL/33/, XFTOL/34/
C/7
C     PARAMETER (AFCTOL=31, BIAS=43, COSMIN=47, DECFAC=22, DELTA0=44,
C    1           DFAC=41, DINIT=38, DLTFDC=42, DLTFDJ=43, DTINIT=39,
C    2           D0INIT=40, EPSLON=19, ETA0=42, FUZZ=45, HUBERC=48,
C    3           INCFAC=23, LMAX0=35, LMAXS=36, PHMNFC=20, PHMXFC=21,
C    4           RDFCMN=24, RDFCMX=25, RFCTOL=32, RLIMIT=46, RSPTOL=49,
C    5           SCTOL=37, SIGMIN=50, TUNER1=26, TUNER2=27, TUNER3=28,
C    6           TUNER4=29, TUNER5=30, XCTOL=33, XFTOL=34)
C/
C
C-------------------------------  BODY  --------------------------------
C
      MACHEP = D1MACH(4)
      V(AFCTOL) = 1.D-20
      IF (MACHEP .GT. 1.D-10) V(AFCTOL) = MACHEP**2
      V(DECFAC) = 0.5D+0
      SQTEPS = DSQRT(D1MACH(4))
      V(DFAC) = 0.6D+0
      V(DELTA0) = SQTEPS
      V(DTINIT) = 1.D-6
      MEPCRT = MACHEP ** (ONE/THREE)
      V(D0INIT) = 1.D+0
      V(EPSLON) = 0.1D+0
      V(INCFAC) = 2.D+0
      V(LMAX0) = 1.D+0
      V(LMAXS) = 1.D+0
      V(PHMNFC) = -0.1D+0
      V(PHMXFC) = 0.1D+0
      V(RDFCMN) = 0.1D+0
      V(RDFCMX) = 4.D+0
      V(RFCTOL) = DMAX1(1.D-10, MEPCRT**2)
      V(SCTOL) = V(RFCTOL)
      V(TUNER1) = 0.1D+0
      V(TUNER2) = 1.D-4
      V(TUNER3) = 0.75D+0
      V(TUNER4) = 0.5D+0
      V(TUNER5) = 0.75D+0
      V(XCTOL) = SQTEPS
      V(XFTOL) = 1.D+2 * MACHEP
C
      IF (ALG .GE. 2) GO TO 10
C
C  ***  REGRESSION  VALUES
C
      V(COSMIN) = DMAX1(1.D-6, 1.D+2 * MACHEP)
      V(DINIT) = 0.D+0
      V(DLTFDC) = MEPCRT
      V(DLTFDJ) = SQTEPS
      V(FUZZ) = 1.5D+0
      V(HUBERC) = 0.7D+0
      V(RLIMIT) = DSQRT(D1MACH(2))*16.
      V(RSPTOL) = 1.D-2
      V(SIGMIN) = 1.D-4
      GO TO 999
C
C  ***  GENERAL OPTIMIZATION VALUES
C
 10   V(BIAS) = 0.8D+0
      V(DINIT) = -1.0D+0
      V(ETA0) = 1.0D+3 * MACHEP
C
 999  RETURN
C  ***  LAST CARD OF DVDFLT FOLLOWS  ***
      END
      SUBROUTINE DVSCPY(P, Y, S)
      save
C
C  ***  SET P-VECTOR Y TO SCALAR S  ***
C
      INTEGER P
      DOUBLE PRECISION S, Y(P)
C
      INTEGER I
C
      DO 10 I = 1, P
 10      Y(I) = S
      RETURN
      END
      SUBROUTINE DVVMUP(N, X, Y, Z, K)
      save
C
C ***  SET X(I) = Y(I) * Z(I)**K, 1 .LE. I .LE. N (FOR K = 1 OR -1)  ***
C
      INTEGER N, K
      DOUBLE PRECISION X(N), Y(N), Z(N)
      INTEGER I
C
      IF (K .GE. 0) GO TO 20
      DO 10 I = 1, N
 10      X(I) = Y(I) / Z(I)
      GO TO 999
C
 20   DO 30 I = 1, N
 30      X(I) = Y(I) * Z(I)
 999  RETURN
C  ***  LAST CARD OF DVVMUP FOLLOWS  ***
      END
      SUBROUTINE DWZBFG (L, N, S, W, Y, Z)
      save
C
C  ***  COMPUTE  Y  AND  Z  FOR  DLUPDT  CORRESPONDING TO BFGS UPDATE.
C
      INTEGER N
      DOUBLE PRECISION L(1), S(N), W(N), Y(N), Z(N)
C     DIMENSION L(N*(N+1)/2)
C
C--------------------------  PARAMETER USAGE  --------------------------
C
C L (I/O) CHOLESKY FACTOR OF HESSIAN, A LOWER TRIANG. MATRIX STORED
C             COMPACTLY BY ROWS.
C N (INPUT) ORDER OF  L  AND LENGTH OF  S,  W,  Y,  Z.
C S (INPUT) THE STEP JUST TAKEN.
C W (OUTPUT) RIGHT SINGULAR VECTOR OF RANK 1 CORRECTION TO L.
C Y (INPUT) CHANGE IN GRADIENTS CORRESPONDING TO S.
C Z (OUTPUT) LEFT SINGULAR VECTOR OF RANK 1 CORRECTION TO L.
C
C-------------------------------  NOTES  -------------------------------
C
C  ***  ALGORITHM NOTES  ***
C
C        WHEN  S  IS COMPUTED IN CERTAIN WAYS, E.G. BY  GQTSTP  OR
C     DBLDOG,  IT IS POSSIBLE TO SAVE N**2/2 OPERATIONS SINCE  (L**T)*S
C     OR  L*(L**T)*S IS THEN KNOWN.
C        IF THE BFGS UPDATE TO L*(L**T) WOULD REDUCE ITS DETERMINANT TO
C     LESS THAN EPS TIMES ITS OLD VALUE, THEN THIS ROUTINE IN EFFECT
C     REPLACES  Y  BY  THETA*Y + (1 - THETA)*L*(L**T)*S,  WHERE  THETA
C     (BETWEEN 0 AND 1) IS CHOSEN TO MAKE THE REDUCTION FACTOR = EPS.
C
C  ***  GENERAL  ***
C
C     CODED BY DAVID M. GAY (FALL 1979).
C     THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH SUPPORTED
C     BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS MCS-7600324 AND
C     MCS-7906671.
C
C------------------------  EXTERNAL QUANTITIES  ------------------------
C
C  ***  FUNCTIONS AND SUBROUTINES CALLED  ***
C
      EXTERNAL DLIVMU, DLTVMU
      DOUBLE PRECISION DDOT
C DLIVMU MULTIPLIES L**-1 TIMES A VECTOR.
C DLTVMU MULTIPLIES L**T TIMES A VECTOR.
C
C  ***  INTRINSIC FUNCTIONS  ***
C/+
      DOUBLE PRECISION DSQRT
C/
C--------------------------  LOCAL VARIABLES  --------------------------
C
      INTEGER I
      DOUBLE PRECISION CS, CY, EPS, EPSRT, ONE, SHS, YS, THETA
C
C  ***  DATA INITIALIZATIONS  ***
C
C/6
      DATA EPS/0.1D+0/, ONE/1.D+0/
C/7
C     PARAMETER (EPS=0.1D+0, ONE=1.D+0)
C/
C
C+++++++++++++++++++++++++++++++  BODY  ++++++++++++++++++++++++++++++++
C
      CALL DLTVMU(N, W, L, S)
      SHS = DDOT(N, W,1,W,1)
      YS = DDOT(N, Y,1,S,1)
      IF (YS .GE. EPS*SHS) GO TO 10
         THETA = (ONE - EPS) * SHS / (SHS - YS)
         EPSRT = DSQRT(EPS)
         CY = THETA / (SHS * EPSRT)
         CS = (ONE + (THETA-ONE)/EPSRT) / SHS
         GO TO 20
 10   CY = ONE / (DSQRT(YS) * DSQRT(SHS))
      CS = ONE / SHS
 20   CALL DLIVMU(N, Z, L, Y)
      DO 30 I = 1, N
 30      Z(I) = CY * Z(I)  -  CS * W(I)
C
 999  RETURN
C  ***  LAST CARD OF DWZBFG FOLLOWS  ***
      END
      INTEGER FUNCTION I1MACH(I)
      save
C***BEGIN PROLOGUE  I1MACH
C***DATE WRITTEN   750101   (YYMMDD)
C***REVISION DATE  910131   (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), AND 8087 BASED MICROS (E.G. IBM PC AND AT&T 6300).
C
C === MACHINE = IEEE.MOST-SIG-BYTE-FIRST
C === MACHINE = IEEE.LEAST-SIG-BYTE-FIRST
C === MACHINE = SUN
C === MACHINE = 68000
C === MACHINE = 8087
C === MACHINE = IBM.PC
C === MACHINE = ATT.3B
C === MACHINE = ATT.7300
C === MACHINE = ATT.6300
       DATA IMACH( 1) /    5 /
       DATA IMACH( 2) /    6 /
       DATA IMACH( 3) /    7 /
       DATA IMACH( 4) /    6 /
       DATA IMACH( 5) /   32 /
       DATA IMACH( 6) /    4 /
       DATA IMACH( 7) /    2 /
       DATA IMACH( 8) /   31 /
       DATA IMACH( 9) / 2147483647 /
       DATA IMACH(10) /    2 /
       DATA IMACH(11) /   24 /
       DATA IMACH(12) / -125 /
       DATA IMACH(13) /  128 /
       DATA IMACH(14) /   53 /
       DATA IMACH(15) / -1021 /
       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 C-120 (NATIVE MODE)
C
C === MACHINE = CONVEX.C1
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 (NATIVE MODE)
C     WITH -R8 OPTION
C
C === MACHINE = CONVEX.C1.R8
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) / -1023 /
C      DATA IMACH(13) /  1023 /
C      DATA IMACH(14) /    53 /
C      DATA IMACH(15) / -1023 /
C      DATA IMACH(16) /  1023 /
C
C     MACHINE CONSTANTS FOR THE CONVEX C-120 (IEEE MODE)
C
C === MACHINE = CONVEX.C1.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 (IEEE MODE)
C     WITH -R8 OPTION
C
C === MACHINE = CONVEX.C1.IEEE.R8
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     HP 9000
C
C === MACHINE = HP.9000
C      DATA IMACH( 1) /     5 /
C      DATA IMACH( 2) /     6 /
C      DATA IMACH( 3) /     6 /
C      DATA IMACH( 4) /     7 /
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) / -1015 /
C      DATA IMACH(16) /  1017 /
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 /CC     MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR).CC === MACHINE = PDP-10.KIC      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
      IF (I .LT. 1  .OR.  I .GT. 16)
     1   CALL XERROR ( 'I1MACH -- I OUT OF BOUNDS',25,1,2)
C
      I1MACH=IMACH(I)
      RETURN
C
      END
      SUBROUTINE XERROR(MESSG,NMESSG,NERR,LEVEL)
      save
C***BEGIN PROLOGUE  XERROR
C***DATE WRITTEN   790801   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  R3C
C***KEYWORDS  ERROR,XERROR PACKAGE
C***AUTHOR  JONES, R. E., (SNLA)
C***PURPOSE  Processes an error (diagnostic) message.
C***DESCRIPTION
C     Abstract
C        XERROR processes a diagnostic message, in a manner
C        determined by the value of LEVEL and the current value
C        of the library error control flag, KONTRL.
C        (See subroutine XSETF for details.)
C
C     Description of Parameters
C      --Input--
C        MESSG - the Hollerith message to be processed, containing
C                no more than 72 characters.
C        NMESSG- the actual number of characters in MESSG.
C        NERR  - the error number associated with this message.
C                NERR must not be zero.
C        LEVEL - error category.
C                =2 means this is an unconditionally fatal error.
C                =1 means this is a recoverable error.  (I.e., it is
C                   non-fatal if XSETF has been appropriately called.)
C                =0 means this is a warning message only.
C                =-1 means this is a warning message which is to be
C                   printed at most once, regardless of how many
C                   times this call is executed.
C
C     Examples
C        CALL XERROR('SMOOTH -- NUM WAS ZERO.',23,1,2)
C        CALL XERROR('INTEG  -- LESS THAN FULL ACCURACY ACHIEVED.',
C                    43,2,1)
C        CALL XERROR(
'ROOTER -- ACTUAL ZERO OF F FOUND BEFORE INTERVAL FC    1ULLY COLLAPSED.',65,3,0)
C        CALL XERROR('EXP    -- UNDERFLOWS BEING SET TO ZERO.',39,1,-1)
C
C     Latest revision ---  19 MAR 1980
C     Written by Ron Jones, with SLATEC Common Math Library Subcommittee
C***REFERENCES  JONES R.E., KAHANER D.K., 
"XERROR, THE SLATEC ERROR-C                 HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES,
C                 1982.
C***ROUTINES CALLED  XERRWV
C***END PROLOGUE  XERROR
      CHARACTER*(*) MESSG
C***FIRST EXECUTABLE STATEMENT  XERROR
      CALL XERRWV(MESSG,NMESSG,NERR,LEVEL,0,0,0,0,0.,0.)
      RETURN
      END
      SUBROUTINE XERRWV(MESSG,NMESSG,NERR,LEVEL,NI,I1,I2,NR,R1,R2)
      save
C***BEGIN PROLOGUE  XERRWV
C***DATE WRITTEN   800319   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  R3C
C***KEYWORDS  ERROR,XERROR PACKAGE
C***AUTHOR  JONES, R. E., (SNLA)
C***PURPOSE  Processes error message allowing 2 integer and two real
C            values to be included in the message.
C***DESCRIPTION
C     Abstract
C        XERRWV processes a diagnostic message, in a manner
C        determined by the value of LEVEL and the current value
C        of the library error control flag, KONTRL.
C        (See subroutine XSETF for details.)
C        In addition, up to two integer values and two real
C        values may be printed along with the message.
C
C     Description of Parameters
C      --Input--
C        MESSG - the Hollerith message to be processed.
C        NMESSG- the actual number of characters in MESSG.
C        NERR  - the error number associated with this message.
C                NERR must not be zero.
C        LEVEL - error category.
C                =2 means this is an unconditionally fatal error.
C                =1 means this is a recoverable error.  (I.e., it is
C                   non-fatal if XSETF has been appropriately called.)
C                =0 means this is a warning message only.
C                =-1 means this is a warning message which is to be
C                   printed at most once, regardless of how many
C                   times this call is executed.
C        NI    - number of integer values to be printed. (0 to 2)
C        I1    - first integer value.
C        I2    - second integer value.
C        NR    - number of real values to be printed. (0 to 2)
C        R1    - first real value.
C        R2    - second real value.
C
C     Examples
C        CALL XERRWV('SMOOTH -- NUM (=I1) WAS ZERO.',29,1,2,
C    1   1,NUM,0,0,0.,0.)
C        CALL XERRWV(





















'QUADXY -- REQUESTED ERROR (R1) LESS THAN MINIMUM (C    1R2).,54,77,1,0,0,0,2,ERRREQ,ERRMIN)CC     Latest revision ---  19 MAR 1980C     Written by Ron Jones, with SLATEC Common Math Library SubcommitteeC***REFERENCES  JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR-C                 HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES,C                 1982.C***ROUTINES CALLED  FDUMP,I1MACH,J4SAVE,XERABT,XERCTL,XERPRT,XERSAV,C                    XGETUAC***END PROLOGUE  XERRWV      CHARACTER*(*) MESSG      CHARACTER*20 LFIRST      CHARACTER*37 FORM      DIMENSION LUN(5)C     GET FLAGSC***FIRST EXECUTABLE STATEMENT  XERRWV      LKNTRL = J4SAVE(2,0,.FALSE.)      MAXMES = J4SAVE(4,0,.FALSE.)C     CHECK FOR VALID INPUT      IF ((NMESSG.GT.0).AND.(NERR.NE.0).AND.     1    (LEVEL.GE.(-1)).AND.(LEVEL.LE.2)) GO TO 10         IF (LKNTRL.GT.0) CALL XERPRT('FATAL ERROR IN...
',17)         CALL XERPRT('XERROR -- INVALID INPUT

',23)         IF (LKNTRL.GT.0) CALL FDUMP         IF (LKNTRL.GT.0) CALL XERPRT('JOB ABORT DUE TO FATAL ERROR.

',     1  29)         IF (LKNTRL.GT.0) CALL XERSAV(' 
',0,0,0,KDUMMY)         CALL XERABT('XERROR -- INVALID INPUT
























',23)         RETURN   10 CONTINUEC     RECORD MESSAGE      JUNK = J4SAVE(1,NERR,.TRUE.)      CALL XERSAV(MESSG,NMESSG,NERR,LEVEL,KOUNT)C     LET USER OVERRIDE      LFIRST = MESSG      LMESSG = NMESSG      LERR = NERR      LLEVEL = LEVEL      CALL XERCTL(LFIRST,LMESSG,LERR,LLEVEL,LKNTRL)C     RESET TO ORIGINAL VALUES      LMESSG = NMESSG      LERR = NERR      LLEVEL = LEVEL      LKNTRL = MAX0(-2,MIN0(2,LKNTRL))      MKNTRL = IABS(LKNTRL)C     DECIDE WHETHER TO PRINT MESSAGE      IF ((LLEVEL.LT.2).AND.(LKNTRL.EQ.0)) GO TO 100      IF (((LLEVEL.EQ.(-1)).AND.(KOUNT.GT.MIN0(1,MAXMES)))     1.OR.((LLEVEL.EQ.0)   .AND.(KOUNT.GT.MAXMES))     2.OR.((LLEVEL.EQ.1)   .AND.(KOUNT.GT.MAXMES).AND.(MKNTRL.EQ.1))     3.OR.((LLEVEL.EQ.2)   .AND.(KOUNT.GT.MAX0(1,MAXMES)))) GO TO 100         IF (LKNTRL.LE.0) GO TO 20            CALL XERPRT(' 


',1)C           INTRODUCTION            IF (LLEVEL.EQ.(-1)) CALL XERPRT     1('WARNING MESSAGE...THIS MESSAGE WILL ONLY BE PRINTED ONCE.
',57)            IF (LLEVEL.EQ.0) CALL XERPRT('WARNING IN...

',13)            IF (LLEVEL.EQ.1) CALL XERPRT     1      ('RECOVERABLE ERROR IN...
',23)            IF (LLEVEL.EQ.2) CALL XERPRT('FATAL ERROR IN...











',17)   20    CONTINUEC        MESSAGE         CALL XERPRT(MESSG,LMESSG)         CALL XGETUA(LUN,NUNIT)         ISIZEI = LOG10(FLOAT(I1MACH(9))) + 1.0         ISIZEF = LOG10(FLOAT(I1MACH(10))**I1MACH(11)) + 1.0         DO 50 KUNIT=1,NUNIT            IUNIT = LUN(KUNIT)            IF (IUNIT.EQ.0) IUNIT = I1MACH(4)            DO 22 I=1,MIN(NI,2)               WRITE (FORM,21) I,ISIZEI   21          FORMAT ('(11X,21HIN ABOVE MESSAGE, I',I1,'=,I',I2,')   





')               IF (I.EQ.1) WRITE (IUNIT,FORM) I1               IF (I.EQ.2) WRITE (IUNIT,FORM) I2   22       CONTINUE            DO 24 I=1,MIN(NR,2)               WRITE (FORM,23) I,ISIZEF+10,ISIZEF   23          FORMAT ('(11X,21HIN ABOVE MESSAGE, R',I1,'=,E
',     1         I2,'.',I2,')




















')               IF (I.EQ.1) WRITE (IUNIT,FORM) R1               IF (I.EQ.2) WRITE (IUNIT,FORM) R2   24       CONTINUE            IF (LKNTRL.LE.0) GO TO 40C              ERROR NUMBER               WRITE (IUNIT,30) LERR   30          FORMAT (15H ERROR NUMBER =,I10)   40       CONTINUE   50    CONTINUEC        TRACE-BACK         IF (LKNTRL.GT.0) CALL FDUMP  100 CONTINUE      IFATAL = 0      IF ((LLEVEL.EQ.2).OR.((LLEVEL.EQ.1).AND.(MKNTRL.EQ.2)))     1IFATAL = 1C     QUIT HERE IF MESSAGE IS NOT FATAL      IF (IFATAL.LE.0) RETURN      IF ((LKNTRL.LE.0).OR.(KOUNT.GT.MAX0(1,MAXMES))) GO TO 120C        PRINT REASON FOR ABORT         IF (LLEVEL.EQ.1) CALL XERPRT     1   ('JOB ABORT DUE TO UNRECOVERED ERROR.

',35)         IF (LLEVEL.EQ.2) CALL XERPRT     1   ('JOB ABORT DUE TO FATAL ERROR.

',29)C        PRINT ERROR SUMMARY         CALL XERSAV(' 


























































































































































































































































































































































',-1,0,0,KDUMMY)  120 CONTINUEC     ABORT      IF ((LLEVEL.EQ.2).AND.(KOUNT.GT.MAX0(1,MAXMES))) LMESSG = 0      CALL XERABT(MESSG,LMESSG)      RETURN      END      SUBROUTINE XERSAV(MESSG,NMESSG,NERR,LEVEL,ICOUNT)      saveC***BEGIN PROLOGUE  XERSAVC***DATE WRITTEN   800319   (YYMMDD)C***REVISION DATE  820801   (YYMMDD)C***CATEGORY NO.  ZC***KEYWORDS  ERROR,XERROR PACKAGEC***AUTHOR  JONES, R. E., (SNLA)C***PURPOSE  Records that an error occurred.C***DESCRIPTIONC     AbstractC        Record that this error occurred.CC     Description of ParametersC     --Input--C       MESSG, NMESSG, NERR, LEVEL are as in XERROR,C       except that when NMESSG=0 the tables will beC       dumped and cleared, and when NMESSG is less than zero theC       tables will be dumped and not cleared.C     --Output--C       ICOUNT will be the number of times this message hasC       been seen, or zero if the table has overflowed andC       does not contain this message specifically.C       When NMESSG=0, ICOUNT will not be altered.CC     Written by Ron Jones, with SLATEC Common Math Library SubcommitteeC     Latest revision ---  19 Mar 1980C***REFERENCES  JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR-C                 HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES,C                 1982.C***ROUTINES CALLED  I1MACH,S88FMT,XGETUAC***END PROLOGUE  XERSAV      INTEGER LUN(5)      CHARACTER*(*) MESSG      CHARACTER*20 MESTAB(10),MES      DIMENSION NERTAB(10),LEVTAB(10),KOUNT(10)      SAVE MESTAB,NERTAB,LEVTAB,KOUNT,KOUNTXC     NEXT TWO DATA STATEMENTS ARE NECESSARY TO PROVIDE A BLANKC     ERROR TABLE INITIALLY      DATA KOUNT(1),KOUNT(2),KOUNT(3),KOUNT(4),KOUNT(5),     1     KOUNT(6),KOUNT(7),KOUNT(8),KOUNT(9),KOUNT(10)     2     /0,0,0,0,0,0,0,0,0,0/      DATA KOUNTX/0/C***FIRST EXECUTABLE STATEMENT  XERSAV      IF (NMESSG.GT.0) GO TO 80C     DUMP THE TABLE         IF (KOUNT(1).EQ.0) RETURNC        PRINT TO EACH UNIT         CALL XGETUA(LUN,NUNIT)         DO 60 KUNIT=1,NUNIT            IUNIT = LUN(KUNIT)            IF (IUNIT.EQ.0) IUNIT = I1MACH(4)C           PRINT TABLE HEADER            WRITE (IUNIT,10)   10       FORMAT (32H0          ERROR MESSAGE SUMMARY/     1      51H MESSAGE START             NERR     LEVEL     COUNT)C           PRINT BODY OF TABLE            DO 20 I=1,10               IF (KOUNT(I).EQ.0) GO TO 30               WRITE (IUNIT,15) MESTAB(I),NERTAB(I),LEVTAB(I),KOUNT(I)   15          FORMAT (1X,A20,3I10)   20       CONTINUE   30       CONTINUEC           PRINT NUMBER OF OTHER ERRORS            IF (KOUNTX.NE.0) WRITE (IUNIT,40) KOUNTX   40       FORMAT (41H0OTHER ERRORS NOT INDIVIDUALLY TABULATED=,I10)            WRITE (IUNIT,50)   50       FORMAT (1X)   60    CONTINUE         IF (NMESSG.LT.0) RETURNC        CLEAR THE ERROR TABLES         DO 70 I=1,10   70       KOUNT(I) = 0         KOUNTX = 0         RETURN   80 CONTINUEC     PROCESS A MESSAGE...C     SEARCH FOR THIS MESSG, OR ELSE AN EMPTY SLOT FOR THIS MESSG,C     OR ELSE DETERMINE THAT THE ERROR TABLE IS FULL.      MES = MESSG      DO 90 I=1,10         II = I         IF (KOUNT(I).EQ.0) GO TO 110         IF (MES.NE.MESTAB(I)) GO TO 90         IF (NERR.NE.NERTAB(I)) GO TO 90         IF (LEVEL.NE.LEVTAB(I)) GO TO 90         GO TO 100   90 CONTINUEC     THREE POSSIBLE CASES...C     TABLE IS FULL         KOUNTX = KOUNTX+1         ICOUNT = 1         RETURNC     MESSAGE FOUND IN TABLE  100    KOUNT(II) = KOUNT(II) + 1         ICOUNT = KOUNT(II)         RETURNC     EMPTY SLOT FOUND FOR NEW MESSAGE  110    MESTAB(II) = MES         NERTAB(II) = NERR         LEVTAB(II) = LEVEL         KOUNT(II)  = 1         ICOUNT = 1         RETURN      END      SUBROUTINE XGETUA(IUNITA,N)      saveC***BEGIN PROLOGUE  XGETUAC***DATE WRITTEN   790801   (YYMMDD)C***REVISION DATE  820801   (YYMMDD)C***CATEGORY NO.  R3CC***KEYWORDS  ERROR,XERROR PACKAGEC***AUTHOR  JONES, R. E., (SNLA)C***PURPOSE  Returns unit number(s) to which error messages are beingC            sent.C***DESCRIPTIONC     AbstractC        XGETUA may be called to determine the unit number or numbersC        to which error messages are being sent.C        These unit numbers may have been set by a call to XSETUN,C        or a call to XSETUA, or may be a default value.CC     Description of ParametersC      --Output--C        IUNIT - an array of one to five unit numbers, dependingC                on the value of N.  A value of zero refers to theC                default unit, as defined by the I1MACH machineC                constant routine.  Only IUNIT(1),...,IUNIT(N) areC                defined by XGETUA.  The values of IUNIT(N+1),...,C                IUNIT(5) are not defined (for N .LT. 5) or alteredC                in any way by XGETUA.C        N     - the number of units to which copies of theC                error messages are being sent.  N will be in theC                range from 1 to 5.CC     Latest revision ---  19 MAR 1980C     Written by Ron Jones, with SLATEC Common Math Library SubcommitteeC***REFERENCES  JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR-C                 HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES,C                 1982.C***ROUTINES CALLED  J4SAVEC***END PROLOGUE  XGETUA      DIMENSION IUNITA(5)C***FIRST EXECUTABLE STATEMENT  XGETUA      N = J4SAVE(5,0,.FALSE.)      DO 30 I=1,N         INDEX = I+4         IF (I.EQ.1) INDEX = 3         IUNITA(I) = J4SAVE(INDEX,0,.FALSE.)   30 CONTINUE      RETURN      END      DOUBLE PRECISION FUNCTION D1MACH(I)      saveC***BEGIN PROLOGUE  D1MACHC***DATE WRITTEN   750101   (YYMMDD)C***REVISION DATE  910131   (YYMMDD)C***CATEGORY NO.  R1C***KEYWORDS  MACHINE CONSTANTSC***AUTHOR  FOX, P. A., (BELL LABS)C           HALL, A. D., (BELL LABS)C           SCHRYER, N. L., (BELL LABS)C***PURPOSE  Returns double precision machine dependent constantsC***DESCRIPTIONCC     This is the CMLIB version of D1MACH, the double precision machineC     constants subroutine originally developed for the PORT library.CC     D1MACH can be used to obtain machine-dependent parametersC     for the local machine environment.  It is a functionC     subprogram with one (input) argument, and can be calledC     as follows, for exampleCC          D = D1MACH(I)CC     where I=1,...,5.  The (output) value of D above isC     determined by the (input) value of I.  The results forC     various values of I are discussed below.CC  Double-precision machine constantsC  D1MACH( 1) = B**(EMIN-1), the smallest positive magnitude.C  D1MACH( 2) = B**EMAX*(1 - B**(-T)), the largest magnitude.C  D1MACH( 3) = B**(-T), the smallest relative spacing.C  D1MACH( 4) = B**(1-T), the largest relative spacing.C  D1MACH( 5) = LOG10(B)C***REFERENCES  FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR AC                 PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICALC                 SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188.C***ROUTINES CALLED  XERRORC***END PROLOGUE  D1MACHC      INTEGER SMALL(4)      INTEGER LARGE(4)      INTEGER RIGHT(4)      INTEGER DIVER(4)      INTEGER LOG10(4)C      DOUBLE PRECISION DMACH(5)C      EQUIVALENCE (DMACH(1),SMALL(1))      EQUIVALENCE (DMACH(2),LARGE(1))      EQUIVALENCE (DMACH(3),RIGHT(1))      EQUIVALENCE (DMACH(4),DIVER(1))      EQUIVALENCE (DMACH(5),LOG10(1))CC     MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&TC     3B SERIES AND MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&TC     PC 7300), IN WHICH THE MOST SIGNIFICANT BYTE IS STORED FIRST.CC === MACHINE = IEEE.MOST-SIG-BYTE-FIRSTC === MACHINE = SUNC === MACHINE = 68000C === MACHINE = ATT.3BC === MACHINE = ATT.7300C      DATA SMALL(1),SMALL(2) /    1048576,          0 /C      DATA LARGE(1),LARGE(2) / 2146435071,         -1 /C      DATA RIGHT(1),RIGHT(2) / 1017118720,          0 /C      DATA DIVER(1),DIVER(2) / 1018167296,          0 /C      DATA LOG10(1),LOG10(2) / 1070810131, 1352628735 /CC     MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES AND 8087-BASEDC     MICROS, SUCH AS THE IBM PC AND AT&T 6300, IN WHICH THE LEASTC     SIGNIFICANT BYTE IS STORED FIRST.CC === MACHINE = IEEE.LEAST-SIG-BYTE-FIRSTC === MACHINE = 8087C === MACHINE = IBM.PCC === MACHINE = ATT.6300       DATA SMALL(1),SMALL(2) /          0,    1048576 /       DATA LARGE(1),LARGE(2) /         -1, 2146435071 /       DATA RIGHT(1),RIGHT(2) /          0, 1017118720 /       DATA DIVER(1),DIVER(2) /          0, 1018167296 /       DATA LOG10(1),LOG10(2) / 1352628735, 1070810131 /CC     MACHINE CONSTANTS FOR AMDAHL MACHINES.CC === MACHINE = AMDAHLC      DATA SMALL(1),SMALL(2) /    1048576,          0 /C      DATA LARGE(1),LARGE(2) / 2147483647,         -1 /C      DATA RIGHT(1),RIGHT(2) /  856686592,          0 /C      DATA DIVER(1),DIVER(2) /  873463808,          0 /C      DATA LOG10(1),LOG10(2) / 1091781651, 1352628735 /CC     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.CC === MACHINE = BURROUGHS.1700C      DATA SMALL(1) / ZC00800000 /C      DATA SMALL(2) / Z000000000 /C      DATA LARGE(1) / ZDFFFFFFFF /C      DATA LARGE(2) / ZFFFFFFFFF /C      DATA RIGHT(1) / ZCC5800000 /C      DATA RIGHT(2) / Z000000000 /C      DATA DIVER(1) / ZCC6800000 /C      DATA DIVER(2) / Z000000000 /C      DATA LOG10(1) / ZD00E730E7 /C      DATA LOG10(2) / ZC77800DC0 /CC     MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM.CC === MACHINE = BURROUGHS.5700C      DATA SMALL(1) / O1771000000000000 /C      DATA SMALL(2) / O0000000000000000 /C      DATA LARGE(1) / O0777777777777777 /C      DATA LARGE(2) / O0007777777777777 /C      DATA RIGHT(1) / O1461000000000000 /C      DATA RIGHT(2) / O0000000000000000 /C      DATA DIVER(1) / O1451000000000000 /C      DATA DIVER(2) / O0000000000000000 /C      DATA LOG10(1) / O1157163034761674 /C      DATA LOG10(2) / O0006677466732724 /CC     MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS.CC === MACHINE = BURROUGHS.6700C === MACHINE = BURROUGHS.7700C      DATA SMALL(1) / O1771000000000000 /C      DATA SMALL(2) / O7770000000000000 /C      DATA LARGE(1) / O0777777777777777 /C      DATA LARGE(2) / O7777777777777777 /C      DATA RIGHT(1) / O1461000000000000 /C      DATA RIGHT(2) / O0000000000000000 /C      DATA DIVER(1) / O1451000000000000 /C      DATA DIVER(2) / O0000000000000000 /C      DATA LOG10(1) / O1157163034761674 /C      DATA LOG10(2) / O0006677466732724 /CC     MACHINE CONSTANTS FOR THE CONVEX C-120 (NATIVE MODE)C     WITH OR WITHOUT -R8 OPTIONCC === MACHINE = CONVEX.C1C === MACHINE = CONVEX.C1.R8C      DATA DMACH(1) / 5.562684646268007D-309 /C      DATA DMACH(2) / 8.988465674311577D+307 /C      DATA DMACH(3) / 1.110223024625157D-016 /C      DATA DMACH(4) / 2.220446049250313D-016 /C      DATA DMACH(5) / 3.010299956639812D-001 /CC     MACHINE CONSTANTS FOR THE CONVEX C-120 (IEEE MODE)C     WITH OR WITHOUT -R8 OPTIONCC === MACHINE = CONVEX.C1.IEEEC === MACHINE = CONVEX.C1.IEEE.R8C      DATA DMACH(1) / 2.225073858507202D-308 /C      DATA DMACH(2) / 1.797693134862315D+308 /C      DATA DMACH(3) / 1.110223024625157D-016 /C      DATA DMACH(4) / 2.220446049250313D-016 /C      DATA DMACH(5) / 3.010299956639812D-001 /CC     MACHINE CONSTANTS FOR THE CYBER 170/180 SERIES USING NOS (FTN5).CC === MACHINE = CYBER.170.NOSC === MACHINE = CYBER.180.NOSC      DATA SMALL(1) / O"00604000000000000000" /C      DATA SMALL(2) / O"00000000000000000000" /C      DATA LARGE(1) / O"37767777777777777777" /C      DATA LARGE(2) / O"37167777777777777777" /C      DATA RIGHT(1) / O"15604000000000000000" /C      DATA RIGHT(2) / O"15000000000000000000" /C      DATA DIVER(1) / O"15614000000000000000" /C      DATA DIVER(2) / O"15010000000000000000" /C      DATA LOG10(1) / O"17164642023241175717" /C      DATA LOG10(2) / O"16367571421742254654" /CC     MACHINE CONSTANTS FOR THE CDC 180 SERIES USING NOS/VECC === MACHINE = CYBER.180.NOS/VEC      DATA SMALL(1) / Z"3001800000000000" /C      DATA SMALL(2) / Z"3001000000000000" /C      DATA LARGE(1) / Z"4FFEFFFFFFFFFFFE" /C      DATA LARGE(2) / Z"4FFE000000000000" /C      DATA RIGHT(1) / Z"3FD2800000000000" /C      DATA RIGHT(2) / Z"3FD2000000000000" /C      DATA DIVER(1) / Z"3FD3800000000000" /C      DATA DIVER(2) / Z"3FD3000000000000" /C      DATA LOG10(1) / Z"3FFF9A209A84FBCF" /C      DATA LOG10(2) / Z"3FFFF7988F8959AC" /CC     MACHINE CONSTANTS FOR THE CYBER 205CC === MACHINE = CYBER.205C      DATA SMALL(1) / X'9000400000000000
' /C      DATA SMALL(2) / X'8FD1000000000000
' /C      DATA LARGE(1) / X'6FFF7FFFFFFFFFFF
' /C      DATA LARGE(2) / X'6FD07FFFFFFFFFFF
' /C      DATA RIGHT(1) / X'FF74400000000000
' /C      DATA RIGHT(2) / X'FF45000000000000
' /C      DATA DIVER(1) / X'FF75400000000000
' /C      DATA DIVER(2) / X'FF46000000000000
' /C      DATA LOG10(1) / X'FFD04D104D427DE7
' /C      DATA LOG10(2) / X'FFA17DE623E2566A












































' /CC     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES.CC === MACHINE = CDC.6000C === MACHINE = CDC.7000C      DATA SMALL(1) / 00604000000000000000B /C      DATA SMALL(2) / 00000000000000000000B /C      DATA LARGE(1) / 37767777777777777777B /C      DATA LARGE(2) / 37167777777777777777B /C      DATA RIGHT(1) / 15604000000000000000B /C      DATA RIGHT(2) / 15000000000000000000B /C      DATA DIVER(1) / 15614000000000000000B /C      DATA DIVER(2) / 15010000000000000000B /C      DATA LOG10(1) / 17164642023241175717B /C      DATA LOG10(2) / 16367571421742254654B /CC     MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3.CC === MACHINE = CRAY.46-BIT-INTEGERC === MACHINE = CRAY.64-BIT-INTEGERC      DATA SMALL(1) / 201354000000000000000B /C      DATA SMALL(2) / 000000000000000000000B /C      DATA LARGE(1) / 577767777777777777777B /C      DATA LARGE(2) / 000007777777777777776B /C      DATA RIGHT(1) / 376434000000000000000B /C      DATA RIGHT(2) / 000000000000000000000B /C      DATA DIVER(1) / 376444000000000000000B /C      DATA DIVER(2) / 000000000000000000000B /C      DATA LOG10(1) / 377774642023241175717B /C      DATA LOG10(2) / 000007571421742254654B /CC     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200CC     NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING LINE -C     STATIC DMACH(5)CC === MACHINE = DATA_GENERAL.ECLIPSE.S/200C      DATA SMALL/20K,3*0/,LARGE/77777K,3*177777K/C      DATA RIGHT/31420K,3*0/,DIVER/32020K,3*0/C      DATA LOG10/40423K,42023K,50237K,74776K/CC     ELXSI 6400CC === MACHINE = ELSXI.6400C      DATA SMALL(1), SMALL(2) / '00100000'X,'00000000
'X /C      DATA LARGE(1), LARGE(2) / '7FEFFFFF'X,'FFFFFFFF
'X /C      DATA RIGHT(1), RIGHT(2) / '3CB00000'X,'00000000
'X /C      DATA DIVER(1), DIVER(2) / '3CC00000'X,'00000000
'X /C      DATA LOG10(1), DIVER(2) / '3FD34413'X,'509F79FF







'X /CC     MACHINE CONSTANTS FOR THE HARRIS 220C     MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7CC === MACHINE = HARRIS.220C === MACHINE = HARRIS.SLASH6C === MACHINE = HARRIS.SLASH7C      DATA SMALL(1),SMALL(2) / '20000000, 
'00000201 /C      DATA LARGE(1),LARGE(2) / '37777777, 
'37777577 /C      DATA RIGHT(1),RIGHT(2) / '20000000, 
'00000333 /C      DATA DIVER(1),DIVER(2) / '20000000, 
'00000334 /C      DATA LOG10(1),LOG10(2) / '23210115, 











































































'10237777 /CC     MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES.C     MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES.CC === MACHINE = HONEYWELL.600/6000C === MACHINE = HONEYWELL.DPS.8/70C      DATA SMALL(1),SMALL(2) / O402400000000, O000000000000 /C      DATA LARGE(1),LARGE(2) / O376777777777, O777777777777 /C      DATA RIGHT(1),RIGHT(2) / O604400000000, O000000000000 /C      DATA DIVER(1),DIVER(2) / O606400000000, O000000000000 /C      DATA LOG10(1),LOG10(2) / O776464202324, O117571775714 /CC      MACHINE CONSTANTS FOR THE HP 2100C      3 WORD DOUBLE PRECISION OPTION WITH FTN4CC === MACHINE = HP.2100.3_WORD_DPC      DATA SMALL(1), SMALL(2), SMALL(3) / 40000B,       0,       1 /C      DATA LARGE(1), LARGE(2), LARGE(3) / 77777B, 177777B, 177776B /C      DATA RIGHT(1), RIGHT(2), RIGHT(3) / 40000B,       0,    265B /C      DATA DIVER(1), DIVER(2), DIVER(3) / 40000B,       0,    276B /C      DATA LOG10(1), LOG10(2), LOG10(3) / 46420B,  46502B,  77777B /CC      MACHINE CONSTANTS FOR THE HP 2100C      4 WORD DOUBLE PRECISION OPTION WITH FTN4CC === MACHINE = HP.2100.4_WORD_DPC      DATA SMALL(1), SMALL(2) /  40000B,       0 /C      DATA SMALL(3), SMALL(4) /       0,       1 /C      DATA LARGE(1), LARGE(2) /  77777B, 177777B /C      DATA LARGE(3), LARGE(4) / 177777B, 177776B /C      DATA RIGHT(1), RIGHT(2) /  40000B,       0 /C      DATA RIGHT(3), RIGHT(4) /       0,    225B /C      DATA DIVER(1), DIVER(2) /  40000B,       0 /C      DATA DIVER(3), DIVER(4) /       0,    227B /C      DATA LOG10(1), LOG10(2) /  46420B,  46502B /C      DATA LOG10(3), LOG10(4) /  76747B, 176377B /CC     HP 9000CC      D1MACH(1) = 2.8480954D-306C      D1MACH(2) = 1.40444776D+306C      D1MACH(3) = 2.22044605D-16C      D1MACH(4) = 4.44089210D-16C      D1MACH(5) = 3.01029996D-1CC === MACHINE = HP.9000C      DATA SMALL(1), SMALL(2) / 00040000000B, 00000000000B /C      DATA LARGE(1), LARGE(2) / 17737777777B, 37777777777B /C      DATA RIGHT(1), RIGHT(2) / 07454000000B, 00000000000B /C      DATA DIVER(1), DIVER(2) / 07460000000B, 00000000000B /C      DATA LOG10(1), LOG10(2) / 07764642023B, 12047674777B /CC     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,C     THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, ANDC     THE INTERDATA 3230 AND INTERDATA 7/32.CC === MACHINE = IBM.360C === MACHINE = IBM.370C === MACHINE = XEROX.SIGMA.5C === MACHINE = XEROX.SIGMA.7C === MACHINE = XEROX.SIGMA.9C === MACHINE = SEL.85C === MACHINE = SEL.86C === MACHINE = INTERDATA.3230C === MACHINE = INTERDATA.7/32C      DATA SMALL(1),SMALL(2) / Z00100000, Z00000000 /C      DATA LARGE(1),LARGE(2) / Z7FFFFFFF, ZFFFFFFFF /C      DATA RIGHT(1),RIGHT(2) / Z33100000, Z00000000 /C      DATA DIVER(1),DIVER(2) / Z34100000, Z00000000 /C      DATA LOG10(1),LOG10(2) / Z41134413, Z509F79FF /CC     MACHINE CONSTANTS FOR THE INTERDATA 8/32C     WITH THE UNIX SYSTEM FORTRAN 77 COMPILER.CC     FOR THE INTERDATA FORTRAN VII COMPILER REPLACEC     THE Z'S SPECIFYING HEX CONSTANTS WITH Y


'S.CC === MACHINE = INTERDATA.8/32.UNIXC      DATA SMALL(1),SMALL(2) / Z'00100000', Z'00000000
' /C      DATA LARGE(1),LARGE(2) / Z'7EFFFFFF', Z'FFFFFFFF
' /C      DATA RIGHT(1),RIGHT(2) / Z'33100000', Z'00000000
' /C      DATA DIVER(1),DIVER(2) / Z'34100000', Z'00000000
' /C      DATA LOG10(1),LOG10(2) / Z'41134413', Z'509F79FF





















































































































' /CC     MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR).CC === MACHINE = PDP-10.KAC      DATA SMALL(1),SMALL(2) / "033400000000, "000000000000 /C      DATA LARGE(1),LARGE(2) / "377777777777, "344777777777 /C      DATA RIGHT(1),RIGHT(2) / "113400000000, "000000000000 /C      DATA DIVER(1),DIVER(2) / "114400000000, "000000000000 /C      DATA LOG10(1),LOG10(2) / "177464202324, "144117571776 /CC     MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR).CC === MACHINE = PDP-10.KIC      DATA SMALL(1),SMALL(2) / "000400000000, "000000000000 /C      DATA LARGE(1),LARGE(2) / "377777777777, "377777777777 /C      DATA RIGHT(1),RIGHT(2) / "103400000000, "000000000000 /C      DATA DIVER(1),DIVER(2) / "104400000000, "000000000000 /C      DATA LOG10(1),LOG10(2) / "177464202324, "047674776746 /CC     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTINGC     32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).CC === MACHINE = PDP-11.32-BITC      DATA SMALL(1),SMALL(2) /    8388608,           0 /C      DATA LARGE(1),LARGE(2) / 2147483647,          -1 /C      DATA RIGHT(1),RIGHT(2) /  612368384,           0 /C      DATA DIVER(1),DIVER(2) /  620756992,           0 /C      DATA LOG10(1),LOG10(2) / 1067065498, -2063872008 /CC      DATA SMALL(1),SMALL(2) / O00040000000, O00000000000 /C      DATA LARGE(1),LARGE(2) / O17777777777, O37777777777 /C      DATA RIGHT(1),RIGHT(2) / O04440000000, O00000000000 /C      DATA DIVER(1),DIVER(2) / O04500000000, O00000000000 /C      DATA LOG10(1),LOG10(2) / O07746420232, O20476747770 /CC     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTINGC     16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).CC === MACHINE = PDP-11.16-BITC      DATA SMALL(1),SMALL(2) /    128,      0 /C      DATA SMALL(3),SMALL(4) /      0,      0 /C      DATA LARGE(1),LARGE(2) /  32767,     -1 /C      DATA LARGE(3),LARGE(4) /     -1,     -1 /C      DATA RIGHT(1),RIGHT(2) /   9344,      0 /C      DATA RIGHT(3),RIGHT(4) /      0,      0 /C      DATA DIVER(1),DIVER(2) /   9472,      0 /C      DATA DIVER(3),DIVER(4) /      0,      0 /C      DATA LOG10(1),LOG10(2) /  16282,   8346 /C      DATA LOG10(3),LOG10(4) / -31493, -12296 /CC      DATA SMALL(1),SMALL(2) / O000200, O000000 /C      DATA SMALL(3),SMALL(4) / O000000, O000000 /C      DATA LARGE(1),LARGE(2) / O077777, O177777 /C      DATA LARGE(3),LARGE(4) / O177777, O177777 /C      DATA RIGHT(1),RIGHT(2) / O022200, O000000 /C      DATA RIGHT(3),RIGHT(4) / O000000, O000000 /C      DATA DIVER(1),DIVER(2) / O022400, O000000 /C      DATA DIVER(3),DIVER(4) / O000000, O000000 /C      DATA LOG10(1),LOG10(2) / O037632, O020232 /C      DATA LOG10(3),LOG10(4) / O102373, O147770 /CC     MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000CC === MACHINE = SEQUENT.BALANCE.8000C      DATA SMALL(1),SMALL(2) / $00000000,  $00100000 /C      DATA LARGE(1),LARGE(2) / $FFFFFFFF,  $7FEFFFFF /C      DATA RIGHT(1),RIGHT(2) / $00000000,  $3CA00000 /C      DATA DIVER(1),DIVER(2) / $00000000,  $3CB00000 /C      DATA LOG10(1),LOG10(2) / $509F79FF,  $3FD34413 /CC     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILERCC === MACHINE = UNIVAC.1100C      DATA SMALL(1),SMALL(2) / O000040000000, O000000000000 /C      DATA LARGE(1),LARGE(2) / O377777777777, O777777777777 /C      DATA RIGHT(1),RIGHT(2) / O170540000000, O000000000000 /C      DATA DIVER(1),DIVER(2) / O170640000000, O000000000000 /C      DATA LOG10(1),LOG10(2) / O177746420232, O411757177572 /CC     MACHINE CONSTANTS FOR VAX 11/780C     (EXPRESSED IN INTEGER AND HEXADECIMAL)C    *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***CC === MACHINE = VAX.11/780C      DATA SMALL(1), SMALL(2) /        128,           0 /C      DATA LARGE(1), LARGE(2) /     -32769,          -1 /C      DATA RIGHT(1), RIGHT(2) /       9344,           0 /C      DATA DIVER(1), DIVER(2) /       9472,           0 /C      DATA LOG10(1), LOG10(2) /  546979738,  -805796613 /CC    ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS***C      DATA SMALL(1), SMALL(2) / Z00000080, Z00000000 /C      DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF /C      DATA RIGHT(1), RIGHT(2) / Z00002480, Z00000000 /C      DATA DIVER(1), DIVER(2) / Z00002500, Z00000000 /C      DATA LOG10(1), LOG10(2) / Z209A3F9A, ZCFF884FB /CC   MACHINE CONSTANTS FOR VAX 11/780 (G-FLOATING)C     (EXPRESSED IN INTEGER AND HEXADECIMAL)C    *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***CC      DATA SMALL(1), SMALL(2) /         16,           0 /C      DATA LARGE(1), LARGE(2) /     -32769,          -1 /C      DATA RIGHT(1), RIGHT(2) /      15552,           0 /C      DATA DIVER(1), DIVER(2) /      15568,           0 /C      DATA LOG10(1), LOG10(2) /  1142112243, 2046775455 /CC    ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS***C      DATA SMALL(1), SMALL(2) / Z00000010, Z00000000 /C      DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF /C      DATA RIGHT(1), RIGHT(2) / Z00003CC0, Z00000000 /C      DATA DIVER(1), DIVER(2) / Z00003CD0, Z00000000 /C      DATA LOG10(1), LOG10(2) / Z44133FF3, Z79FF509F /CCC***FIRST EXECUTABLE STATEMENT  D1MACH      IF (I .LT. 1  .OR.  I .GT. 5)     1   CALL XERROR( 'D1MACH -- I OUT OF BOUNDS













































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































',25,1,2)C      D1MACH = DMACH(I)      RETURNC      END      SUBROUTINE DASSST(D, IV, P, STEP, STLSTG, V, X, X0)      saveCC  ***  ASSESS CANDIDATE STEP (***SOL VERSION 2.3)  ***C      INTEGER P, IV(32)      DOUBLE PRECISION D(P), STEP(P), STLSTG(P), V(37), X(P), X0(P)CC  ***  PURPOSE  ***CC        THIS SUBROUTINE IS CALLED BY AN UNCONSTRAINED MINIMIZATIONC     ROUTINE TO ASSESS THE NEXT CANDIDATE STEP.  IT MAY RECOMMEND ONEC     OF SEVERAL COURSES OF ACTION, SUCH AS ACCEPTING THE STEP, RECOM-C     PUTING IT USING THE SAME OR A NEW QUADRATIC MODEL, OR HALTING DUEC     TO CONVERGENCE OR FALSE CONVERGENCE.  SEE THE RETURN CODE LISTINGC     BELOW.CC--------------------------  PARAMETER USAGE  --------------------------CC     IV (I/O) INTEGER PARAMETER AND SCRATCH VECTOR -- SEE DESCRIPTIONC             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 AC             BETTER OBJECTIVE FUNCTION REDUCTION, IN WHICH CASE STLSTGC             WILL HAVE BEEN COPIED TO STEP.C STLSTG (I/O) WHEN ASSESS RECOMMENDS RECOMPUTING STEP EVEN THOUGH THEC             CURRENT (OR A PREVIOUS) STEP YIELDS AN OBJECTIVE FUNC-C             TION DECREASE, IT SAVES IN STLSTG THE STEP THAT GAVE THEC             BEST FUNCTION REDUCTION SEEN SO FAR (IN THE CURRENT ITERA-C             TION).  IF THE RECOMPUTED STEP YIELDS A LARGER FUNCTIONC             VALUE, THEN STEP IS RESTORED FROM STLSTG ANDC             X = X0 + STEP IS RECOMPUTED.C      V (I/O) REAL PARAMETER AND SCRATCH VECTOR -- SEE DESCRIPTIONC             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 EARLIERC             STEP YIELDED A BIGGER FUNCTION DECREASE, THEN X ISC             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 THEC             START OF THE CURRENT ITERATION).CC  ***  IV VALUES REFERENCED  ***CC    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 ISC             SET WHEN STEP IS DEFINITELY TO BE ACCEPTED).  ON INPUTC             AFTER STEP HAS BEEN RECOMPUTED, IV(IRC) SHOULD BEC             UNCHANGED SINCE THE PREVIOUS RETURN OF ASSESS.C                ON OUTPUT, IV(IRC) IS A RETURN CODE HAVING ONE OF THEC             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 GRADIENTC                       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 NOTC                       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 IDENTIFYINGC             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 BIGGESTC             FUNCTION REDUCTION THIS ITERATION.  IV(NFGCAL) REMAINSC             UNCHANGED UNTIL A FUNCTION REDUCTION IS OBTAINED.C IV(RADINC) (I/O) THE NUMBER OF RADIUS INCREASES (OR MINUS THE NUMBERC             OF DECREASES) SO FAR THIS ITERATION.C IV(RESTOR) (OUT) SET TO 0 UNLESS X AND V(F) HAVE BEEN RESTORED, INC             WHICH CASE ASSESS SETS IV(RESTOR) = 1.C  IV(STAGE) (I/O) COUNT OF THE NUMBER OF MODELS TRIED SO FAR IN THEC             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 ITC             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 CAUSEDC             OVERFLOW).C   IV(XIRC) (I/O) VALUE THAT IV(IRC) WOULD HAVE IN THE ABSENCE OFC             CONVERGENCE, FALSE CONVERGENCE, AND OVERSIZED STEPS.CC  ***  V VALUES REFERENCED  ***CC V(AFCTOL) (IN)  ABSOLUTE FUNCTION CONVERGENCE TOLERANCE.  IF THEC             ABSOLUTE VALUE OF THE CURRENT FUNCTION VALUE V(F) IS LESSC             THAN V(AFCTOL), THEN ASSESS RETURNS WITH IV(IRC) = 10.C V(DECFAC) (IN)  FACTOR BY WHICH TO DECREASE RADIUS WHEN IV(TOOBIG) ISC             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 OUTPUTC             VALUE OF V(F) IF AN EARLIER STEP GAVE A BIGGER FUNCTIONC             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 TWICEC             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 IFC             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 FORC             A STEP OF LENGTH V(LMAXS) (BY A RETURN WITH IV(IRC) = 6).C V(NREDUC) (I/O)  FUNCTION REDUCTION PREDICTED BY QUADRATIC MODEL FORC             NEWTON STEP.  IF ASSESS IS CALLED WITH IV(IRC) = 6, I.E.,C             IF V(PREDUC) HAS BEEN COMPUTED WITH RADIUS = V(LMAXS) FORC             USE IN THE SINGULAR CONVERVENCE TEST, THEN V(NREDUC) ISC             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 FORC             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 THEC             OUTPUT VALUE OF V(DSTNRM) OR THE 2-NORM OFC             DIAG(NEWD)*STEP  FOR THE OUTPUT VALUE OF STEP AND THEC             UPDATED VERSION, NEWD, OF THE SCALE VECTOR D.  FORC             IV(IRC) = 3, V(RADFAC) = 1.0 IS RETURNED.C V(RDFCMN) (IN)  MINIMUM VALUE FOR V(RADFAC) IN TERMS OF THE INPUTC             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, COMPUTEDC             BY FUNCTION  DRELST  ASC                 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 XC             AND STEP.  OTHERWISE IT IS COMPUTED USING THE INPUTC             VALUES.C V(RFCTOL) (IN)  RELATIVE FUNCTION CONVERGENCE TOLERANCE.  IF THEC             ACTUAL FUNCTION REDUCTION IS AT MOST TWICE WHAT WAS PRE-C             DICTED AND  V(NREDUC) .LE. V(RFCTOL)*ABS(V(F0)),  THENC             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 FUNCTIONC             REDUCTION WAS MUCH LESS THAN EXPECTED.  SUGGESTEDC             VALUE = 0.1.C V(TUNER2) (IN)  TUNING CONSTANT USED TO DECIDE IF THE FUNCTIONC             REDUCTION WAS LARGE ENOUGH TO ACCEPT STEP.  SUGGESTEDC             VALUE = 10**-4.C V(TUNER3) (IN)  TUNING CONSTANT USED TO DECIDE IF THE RADIUSC             SHOULD BE INCREASED.  SUGGESTED VALUE = 0.75.C  V(XCTOL) (IN)  X-CONVERGENCE CRITERION.  IF STEP IS A NEWTON STEPC             (V(STPPAR) = 0) HAVING V(RELDX) .LE. V(XCTOL) AND GIVINGC             AT MOST TWICE THE PREDICTED FUNCTION DECREASE, THENC             ASSESS RETURNS IV(IRC) = 7 OR 9.C  V(XFTOL) (IN)  FALSE CONVERGENCE TOLERANCE.  IF STEP GAVE NO OR ONLYC             A SMALL FUNCTION DECREASE AND V(RELDX) .LE. V(XFTOL),C             THEN ASSESS RETURNS WITH IV(IRC) = 12.CC-------------------------------  NOTES  -------------------------------CC  ***  APPLICATION AND USAGE RESTRICTIONS  ***CC        THIS ROUTINE IS CALLED AS PART OF THE NL2SOL (NONLINEARC     LEAST-SQUARES) PACKAGE.  IT MAY BE USED IN ANY UNCONSTRAINEDC     MINIMIZATION SOLVER THAT USES DOGLEG, GOLDFELD-QUANDT-TROTTER,C     OR LEVENBERG-MARQUARDT STEPS.CC  ***  ALGORITHM NOTES  ***CC        SEE (1) FOR FURTHER DISCUSSION OF THE ASSESSING AND MODELC     SWITCHING STRATEGIES.  WHILE NL2SOL CONSIDERS ONLY TWO MODELS,C     ASSESS IS DESIGNED TO HANDLE ANY NUMBER OF MODELS.CC  ***  USAGE NOTES  ***CC        ON THE FIRST CALL OF AN ITERATION, ONLY THE I/O VARIABLESC     STEP, X, IV(IRC), IV(MODEL), V(F), V(DSTNRM), V(GTSTEP), ANDC     V(PREDUC) NEED HAVE BEEN INITIALIZED.  BETWEEN CALLS, NO I/OC     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 CANC     CHANGE THE STOPPING TOLERANCES AND CALL ASSESS AGAIN, IN WHICHC     CASE THE STOPPING TESTS WILL BE REPEATED.CC  ***  REFERENCES  ***CC     (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.CC     (2) POWELL, M.J.D. (1970)  A FORTRAN SUBROUTINE FOR SOLVINGC        SYSTEMS OF NONLINEAR ALGEBRAIC EQUATIONS, IN NUMERICALC        METHODS FOR NONLINEAR ALGEBRAIC EQUATIONS, EDITED BYC        P. RABINOWITZ, GORDON AND BREACH, LONDON.CC  ***  HISTORY  ***CC        JOHN DENNIS DESIGNED MUCH OF THIS ROUTINE, STARTING WITHC     IDEAS IN (2). ROY WELSCH SUGGESTED THE MODEL SWITCHING STRATEGY.C        DAVID GAY AND STEPHEN PETERS CAST THIS SUBROUTINE INTO A MOREC     PORTABLE FORM (WINTER 1977), AND DAVID GAY CAST IT INTO ITSC     PRESENT FORM (FALL 1978).CC  ***  GENERAL  ***CC     THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCHC     SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTSC     MCS-7600324, DCR75-10143, 76-14311DSS, MCS76-11989, ANDC     MCS-7906671.CC------------------------  EXTERNAL QUANTITIES  ------------------------CC  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***C      EXTERNAL DRELST      DOUBLE PRECISION DRELSTCC DRELST... COMPUTES V(RELDX) = RELATIVE STEP SIZE.CC  ***  INTRINSIC FUNCTIONS  ***C/+      DOUBLE PRECISION DABS, DMAX1C/C  ***  NO COMMON BLOCKS  ***CC--------------------------  LOCAL VARIABLES  --------------------------C      LOGICAL GOODX      INTEGER I, NFC      DOUBLE PRECISION EMAX, EMAXS, GTS, HALF, ONE, RELDX1, RFAC1, TWO,     1                 XMAX, ZEROCC  ***  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, XIRCCC  ***  DATA INITIALIZATIONS  ***CC/6      DATA HALF/0.5D+0/, ONE/1.D+0/, TWO/2.D+0/, ZERO/0.D+0/C/7C     PARAMETER (HALF=0.5D+0, ONE=1.D+0, TWO=2.D+0, ZERO=0.D+0)C/CC/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/7C     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/7C     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/CC+++++++++++++++++++++++++++++++  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 999CC  ***  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 60CC  ***  STEP WAS RECOMPUTED WITH NEW MODEL OR SMALLER RADIUS  ***C  ***  FIRST DECIDE WHICH  ***C 20   IF (IV(MODEL) .NE. IV(MLSTGD)) GO TO 30C        ***  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 90CC  ***  A NEW MODEL IS BEING TRIED.  DECIDE WHETHER TO KEEP IT.  ***C 30   IV(STAGE) = IV(STAGE) + 1CC     ***  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 50CC        ***  STEP WAS RECOMPUTED BECAUSE IT WAS TOO BIG.  ***C         IF (IV(TOOBIG) .NE. 0) GO TO 60CC        ***  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), IC 50   IF (IV(TOOBIG) .EQ. 0) GO TO 70CC  ***  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 999C 70   IF (V(F) .LT. V(FLSTGD)) GO TO 90CC     *** 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) = 1CC     ***  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.CCC  ***  COMPUTE RELATIVE CHANGE IN X BY CURRENT STEP  ***C 90   RELDX1 = DRELST(P, D, X, X0)CC  ***  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     CONTINUEC 110  V(FDIF) = V(F0) - V(F)      IF (V(FDIF) .GT. V(TUNER2) * V(PREDUC)) GO TO 140CC        ***  NO (OR ONLY A TRIVIAL) FUNCTION DECREASEC        ***  -- SO TRY NEW MODEL OR SMALLER RADIUSC         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 DCOPY(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 160CC  ***  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 190CC  ***  DECREASE WAS MUCH LESS THAN PREDICTED -- EITHER CHANGE MODELSC  ***  OR ACCEPT STEP WITH DECREASED RADIUS.C      IF (IV(STAGE) .GE. IV(STGLIM)) GO TO 150C        ***  CONSIDER SWITCHING MODELS  ***         IV(IRC) = 2         GO TO 160CC     ***  ACCEPT STEP WITH DECREASED RADIUS  ***C 150  IV(IRC) = 4CC  ***  SET V(RADFAC) TO FLETCHER*S DECREASE FACTOR  ***C 160  IV(XIRC) = IV(IRC)      EMAX = V(GTSTEP) + V(FDIF)      V(RADFAC) = HALF * RFAC1      IF (EMAX .LT. V(GTSTEP)) V(RADFAC) = RFAC1 * DMAX1(V(RDFCMN),     1                                           HALF * V(GTSTEP)/EMAX)CC  ***  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 230C 180  IV(IRC) = 12      GO TO 240CC  ***  HANDLE GOOD FUNCTION DECREASE  ***C 190  IF (V(FDIF) .LT. (-V(TUNER3) * V(GTSTEP))) GO TO 210CC     ***  INCREASING RADIUS LOOKS WORTHWHILE.  SEE IF WE JUSTC     ***  RECOMPUTED STEP WITH A DECREASED RADIUS OR RESTORED STEPC     ***  AFTER RECOMPUTING IT WITH A LARGER RADIUS.C      IF (IV(RADINC) .LT. 0) GO TO 210      IF (IV(RESTOR) .EQ. 1) GO TO 210CC        ***  WE DID NOT.  TRY A LONGER STEP UNLESS THIS WAS A NEWTONC        ***  STEP.C         V(RADFAC) = V(RDFCMX)         GTS = V(GTSTEP)         IF (V(FDIF) .LT. (HALF/V(RADFAC) - ONE) * GTS)     1            V(RADFAC) = DMAX1(V(INCFAC), HALF*GTS/(GTS + V(FDIF)))         IV(IRC) = 4         IF (V(STPPAR) .EQ. ZERO) GO TO 230C             ***  STEP WAS NOT A NEWTON STEP.  RECOMPUTE IT WITHC             ***  A LARGER RADIUS.              IV(IRC) = 5              IV(RADINC) = IV(RADINC) + 1CC  ***  SAVE VALUES CORRESPONDING TO GOOD STEP  ***C 200  V(FLSTGD) = V(F)      IV(MLSTGD) = IV(MODEL)      CALL DCOPY(P, STEP,1,STLSTG,1)      V(DSTSAV) = V(DSTNRM)      IV(NFGCAL) = NFC      V(PLSTGD) = V(PREDUC)      V(GTSLST) = V(GTSTEP)      GO TO 230CC  ***  ACCEPT STEP WITH RADIUS UNCHANGED  ***C 210  V(RADFAC) = ONE      IV(IRC) = 3      GO TO 230CC  ***  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 240CC  ***  PERFORM CONVERGENCE TESTS  ***C 230  IV(XIRC) = IV(IRC) 240  IF (DABS(V(F)) .LT. V(AFCTOL)) IV(IRC) = 10      IF (HALF * V(FDIF) .GT. V(PREDUC)) GO TO 999      EMAX = V(RFCTOL) * DABS(V(F0))      EMAXS = V(SCTOL) * DABS(V(F0))      IF (V(DSTNRM) .GT. V(LMAXS) .AND. V(PREDUC) .LE. EMAXS)     1                       IV(IRC) = 11      IF (V(DST0) .LT. ZERO) GO TO 250      I = 0      IF ((V(NREDUC) .GT. ZERO .AND. V(NREDUC) .LE. EMAX) .OR.     1    (V(NREDUC) .EQ. ZERO. AND. V(PREDUC) .EQ. ZERO))  I = 2      IF (V(STPPAR) .EQ. ZERO .AND. V(RELDX) .LE. V(XCTOL)     1                        .AND. GOODX)                  I = I + 1      IF (I .GT. 0) IV(IRC) = I + 6CC  ***  CONSIDER RECOMPUTING STEP OF LENGTH V(LMAXS) FOR SINGULARC  ***  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 290CC  ***  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 DCOPY(P, STEP,1,STLSTG,1)      GO TO 999CC  ***  PERFORM SINGULAR CONVERGENCE TEST WITH RECOMPUTED V(PREDUC)  ***C 280  V(GTSTEP) = V(GTSLST)      V(DSTNRM) = DABS(V(DSTSAV))      CALL DCOPY(P, STLSTG,1,STEP,1)      IV(IRC) = IV(XIRC)      IF (V(DSTSAV) .LE. ZERO) IV(IRC) = 12      V(NREDUC) = -V(PREDUC)      V(PREDUC) = V(PLSTGD) 290  IF (-V(NREDUC) .LE. V(RFCTOL) * DABS(V(F0))) IV(IRC) = 11C 999  RETURNCC  ***  LAST CARD OF ASSESS FOLLOWS  ***      END      SUBROUTINE DCOPY(N,DX,INCX,DY,INCY)      saveC***BEGIN PROLOGUE  DCOPYC***DATE WRITTEN   791001   (YYMMDD)C***REVISION DATE  820801   (YYMMDD)C***CATEGORY NO.  D1A5C***KEYWORDS  BLAS,COPY,DOUBLE PRECISION,LINEAR ALGEBRA,VECTORC***AUTHOR  LAWSON, C. L., (JPL)C           HANSON, R. J., (SNLA)C           KINCAID, D. R., (U. OF TEXAS)C           KROGH, F. T., (JPL)C***PURPOSE  D.P. vector copy y = xC***DESCRIPTIONCC                B L A S  SubprogramC    Description of ParametersCC     --Input--C        N  number of elements in input vector(s)C       DX  double precision vector with N elementsC     INCX  storage spacing between elements of DXC       DY  double precision vector with N elementsC     INCY  storage spacing between elements of DYCC     --Output--C       DY  copy of vector DX (unchanged if N .LE. 0)CC     Copy double precision DX to double precision DY.C     For I = 0 to N-1, copy DX(LX+I*INCX) to DY(LY+I*INCY),C     where LX = 1 if INCX .GE. 0, else LX = (-INCX)*N, and LY isC     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 MATHEMATICALC                 SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323C***ROUTINES CALLED  (NONE)C***END PROLOGUE  DCOPYC      DOUBLE PRECISION DX(1),DY(1)C***FIRST EXECUTABLE STATEMENT  DCOPY      IF(N.LE.0)RETURN      IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60    5 CONTINUECC        CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS.C      IX = 1      IY = 1      IF(INCX.LT.0)IX = (-N+1)*INCX + 1      IF(INCY.LT.0)IY = (-N+1)*INCY + 1      DO 10 I = 1,N        DY(IY) = DX(IX)        IX = IX + INCX        IY = IY + INCY   10 CONTINUE      RETURNCC        CODE FOR BOTH INCREMENTS EQUAL TO 1CCC        CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 7.C   20 M = MOD(N,7)      IF( M .EQ. 0 ) GO TO 40      DO 30 I = 1,M        DY(I) = DX(I)   30 CONTINUE      IF( N .LT. 7 ) RETURN   40 MP1 = M + 1      DO 50 I = MP1,N,7        DY(I) = DX(I)        DY(I + 1) = DX(I + 1)        DY(I + 2) = DX(I + 2)        DY(I + 3) = DX(I + 3)        DY(I + 4) = DX(I + 4)        DY(I + 5) = DX(I + 5)        DY(I + 6) = DX(I + 6)   50 CONTINUE      RETURNCC        CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS.C   60 CONTINUE      NS=N*INCX          DO 70 I=1,NS,INCX          DY(I) = DX(I)   70     CONTINUE      RETURN      END      SUBROUTINE DDBDOG(DIG, G, LV, N, NWTSTP, STEP, V)      saveCC  ***  COMPUTE DOUBLE DOGLEG STEP  ***CC  ***  PARAMETER DECLARATIONS  ***C      INTEGER LV, N      DOUBLE PRECISION DIG(N), G(N), NWTSTP(N), STEP(N), V(LV)CC  ***  PURPOSE  ***CC        THIS SUBROUTINE COMPUTES A CANDIDATE STEP (FOR USE IN AN UNCON-C     STRAINED MINIMIZATION CODE) BY THE DOUBLE DOGLEG ALGORITHM OFC     DENNIS AND MEI (REF. 1), WHICH IS A VARIATION ON POWELL*S DOGLEGC     SCHEME (REF. 2, P. 95).CC--------------------------  PARAMETER USAGE  --------------------------CC    DIG (INPUT) DIAG(D)**-2 * G -- SEE ALGORITHM NOTES.C      G (INPUT) THE CURRENT GRADIENT VECTOR.C     LV (INPUT) LENGTH OF V.C      N (INPUT) NUMBER OF COMPONENTS IN  DIG, G, NWTSTP,  AND  STEP.C NWTSTP (INPUT) NEGATIVE NEWTON STEP -- SEE ALGORITHM NOTES.C   STEP (OUTPUT) THE COMPUTED STEP.C      V (I/O) VALUES ARRAY, THE FOLLOWING COMPONENTS OF WHICH AREC             USED HERE...C V(BIAS)   (INPUT) BIAS FOR RELAXED NEWTON STEP, WHICH IS V(BIAS) OFC             THE WAY FROM THE FULL NEWTON TO THE FULLY RELAXED NEWTONC             STEP.  RECOMMENDED VALUE = 0.8 .C V(DGNORM) (INPUT) 2-NORM OF DIAG(D)**-1 * G -- SEE ALGORITHM NOTES.C V(DSTNRM) (OUTPUT) 2-NORM OF DIAG(D) * STEP, WHICH IS V(RADIUS)C             UNLESS V(STPPAR) = 0 -- SEE ALGORITHM NOTES.C V(DST0) (INPUT) 2-NORM OF DIAG(D) * NWTSTP -- SEE ALGORITHM NOTES.C V(GRDFAC) (OUTPUT) THE COEFFICIENT OF  DIG  IN THE STEP RETURNED --C             STEP(I) = V(GRDFAC)*DIG(I) + V(NWTFAC)*NWTSTP(I).C V(GTHG)   (INPUT) SQUARE-ROOT OF (DIG**T) * (HESSIAN) * DIG -- SEEC             ALGORITHM NOTES.C V(GTSTEP) (OUTPUT) INNER PRODUCT BETWEEN G AND STEP.C V(NREDUC) (OUTPUT) FUNCTION REDUCTION PREDICTED FOR THE FULL NEWTONC             STEP.C V(NWTFAC) (OUTPUT) THE COEFFICIENT OF  NWTSTP  IN THE STEP RETURNED --C             SEE V(GRDFAC) ABOVE.C V(PREDUC) (OUTPUT) FUNCTION REDUCTION PREDICTED FOR THE STEP RETURNED.C V(RADIUS) (INPUT) THE TRUST REGION RADIUS.  D TIMES THE STEP RETURNEDC             HAS 2-NORM V(RADIUS) UNLESS V(STPPAR) = 0.C V(STPPAR) (OUTPUT) CODE TELLING HOW STEP WAS COMPUTED... 0 MEANS AC             FULL NEWTON STEP.  BETWEEN 0 AND 1 MEANS V(STPPAR) OF THEC             WAY FROM THE NEWTON TO THE RELAXED NEWTON STEP.  BETWEENC             1 AND 2 MEANS A TRUE DOUBLE DOGLEG STEP, V(STPPAR) - 1 OFC             THE WAY FROM THE RELAXED NEWTON TO THE CAUCHY STEP.C             GREATER THAN 2 MEANS 1 / (V(STPPAR) - 1) TIMES THE CAUCHYC             STEP.CC-------------------------------  NOTES  -------------------------------CC  ***  ALGORITHM NOTES  ***CC        LET  G  AND  H  BE THE CURRENT GRADIENT AND HESSIAN APPROXIMA-C     TION RESPECTIVELY AND LET D BE THE CURRENT SCALE VECTOR.  THISC     ROUTINE ASSUMES DIG = DIAG(D)**-2 * G  AND  NWTSTP = H**-1 * G.C     THE STEP COMPUTED IS THE SAME ONE WOULD GET BY REPLACING G AND HC     BY  DIAG(D)**-1 * G  AND  DIAG(D)**-1 * H * DIAG(D)**-1,C     COMPUTING STEP, AND TRANSLATING STEP BACK TO THE ORIGINALC     VARIABLES, I.E., PREMULTIPLYING IT BY DIAG(D)**-1.CC  ***  REFERENCES  ***CC 1.  DENNIS, J.E., AND MEI, H.H.W. (1979), TWO NEW UNCONSTRAINED OPTI-C             MIZATION ALGORITHMS WHICH USE FUNCTION AND GRADIENTC             VALUES, J. OPTIM. THEORY APPLIC. 28, PP. 453-482.C 2. POWELL, M.J.D. (1970), A HYBRID METHOD FOR NON-LINEAR EQUATIONS,C             IN NUMERICAL METHODS FOR NON-LINEAR EQUATIONS, EDITED BYC             P. RABINOWITZ, GORDON AND BREACH, LONDON.CC  ***  GENERAL  ***CC     CODED BY DAVID M. GAY.C     THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH SUPPORTEDC     BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS MCS-7600324 ANDC     MCS-7906671.CC------------------------  EXTERNAL QUANTITIES  ------------------------CC  ***  FUNCTIONS AND SUBROUTINES CALLED  ***C      DOUBLE PRECISION DDOTCCC  ***  INTRINSIC FUNCTIONS  ***C/+      DOUBLE PRECISION DSQRTC/C--------------------------  LOCAL VARIABLES  --------------------------C      INTEGER I      DOUBLE PRECISION CFACT, CNORM, CTRNWT, GHINVG, FEMNSQ, GNORM,     1                 NWTNRM, RELAX, RLAMBD, T, T1, T2      DOUBLE PRECISION HALF, ONE, TWO, ZEROCC  ***  V SUBSCRIPTS  ***C      INTEGER BIAS, DGNORM, DSTNRM, DST0, GRDFAC, GTHG, GTSTEP,     1        NREDUC, NWTFAC, PREDUC, RADIUS, STPPARCC  ***  DATA INITIALIZATIONS  ***CC/6      DATA HALF/0.5D+0/, ONE/1.D+0/, TWO/2.D+0/, ZERO/0.D+0/C/7C     PARAMETER (HALF=0.5D+0, ONE=1.D+0, TWO=2.D+0, ZERO=0.D+0)C/CC/6      DATA BIAS/43/, DGNORM/1/, DSTNRM/2/, DST0/3/, GRDFAC/45/,     1     GTHG/44/, GTSTEP/4/, NREDUC/6/, NWTFAC/46/, PREDUC/7/,     2     RADIUS/8/, STPPAR/5/C/7C     PARAMETER (BIAS=43, DGNORM=1, DSTNRM=2, DST0=3, GRDFAC=45,C    1           GTHG=44, GTSTEP=4, NREDUC=6, NWTFAC=46, PREDUC=7,C    2           RADIUS=8, STPPAR=5)C/CC+++++++++++++++++++++++++++++++  BODY  ++++++++++++++++++++++++++++++++C      NWTNRM = V(DST0)      RLAMBD = ONE      IF (NWTNRM .GT. ZERO) RLAMBD = V(RADIUS) / NWTNRM      GNORM = V(DGNORM)      DO 10 I = 1, N 10      STEP(I) = G(I) / GNORM      GHINVG = DDOT(N, STEP,1,NWTSTP,1)      V(NREDUC) = HALF * GHINVG * GNORM      V(GRDFAC) = ZERO      V(NWTFAC) = ZERO      IF (RLAMBD .LT. ONE) GO TO 30CC        ***  THE NEWTON STEP IS INSIDE THE TRUST REGION  ***C         V(STPPAR) = ZERO         V(DSTNRM) = NWTNRM         V(GTSTEP) = -GHINVG * GNORM         V(PREDUC) = V(NREDUC)         V(NWTFAC) = -ONE         DO 20 I = 1, N 20           STEP(I) = -NWTSTP(I)         GO TO 999C 30   V(DSTNRM) = V(RADIUS)      CFACT = (GNORM / V(GTHG))**2C     ***  CAUCHY STEP = -CFACT * G.      CNORM = GNORM * CFACT      RELAX = ONE - V(BIAS) * (ONE - CNORM/GHINVG)      IF (RLAMBD .LT. RELAX) GO TO 50CC        ***  STEP IS BETWEEN RELAXED NEWTON AND FULL NEWTON STEPS  ***C         V(STPPAR)  =  ONE  -  (RLAMBD - RELAX) / (ONE - RELAX)         T = -RLAMBD         V(GTSTEP) = T * GHINVG * GNORM         V(PREDUC) = RLAMBD * (ONE - HALF*RLAMBD) * GHINVG * GNORM         V(NWTFAC) = T         DO 40 I = 1, N 40           STEP(I) = T * NWTSTP(I)         GO TO 999C 50   IF (CNORM .LT. V(RADIUS)) GO TO 70CC        ***  THE CAUCHY STEP LIES OUTSIDE THE TRUST REGION --C        ***  STEP = SCALED CAUCHY STEP  ***C         T = -V(RADIUS) / GNORM         V(GRDFAC) = T         V(STPPAR) = ONE  +  CNORM / V(RADIUS)         V(GTSTEP) = -V(RADIUS) * GNORM      V(PREDUC) = V(RADIUS)*(GNORM - HALF*V(RADIUS)*(V(GTHG)/GNORM)**2)         DO 60 I = 1, N 60           STEP(I) = T * DIG(I)         GO TO 999CC     ***  COMPUTE DOGLEG STEP BETWEEN CAUCHY AND RELAXED NEWTON  ***C     ***  FEMUR = RELAXED NEWTON STEP MINUS CAUCHY STEP  ***C 70   CTRNWT = CFACT * RELAX * GHINVG / GNORMC     *** CTRNWT = INNER PROD. OF CAUCHY AND RELAXED NEWTON STEPS,C     *** SCALED BY GNORM**-2.      T1 = CTRNWT - CFACT**2C     ***  T1 = INNER PROD. OF FEMUR AND CAUCHY STEP, SCALED BYC     ***  GNORM**-2.      T2 = (V(RADIUS)/GNORM)**2 - CFACT**2      FEMNSQ = (RELAX*NWTNRM/GNORM)**2 - CTRNWT - T1C     ***  FEMNSQ = SQUARE OF 2-NORM OF FEMUR, SCALED BY GNORM**-2.      T = T2 / (T1 + DSQRT(T1**2 + FEMNSQ*T2))C     ***  DOGLEG STEP  =  CAUCHY STEP  +  T * FEMUR.      T1 = (T - ONE) * CFACT      V(GRDFAC) = T1      T2 = -T * RELAX      V(NWTFAC) = T2      V(STPPAR) = TWO - T      V(GTSTEP) = GNORM * (T1*GNORM + T2*GHINVG)      V(PREDUC) = -(T1*GNORM) * ((T2 + ONE)*GNORM)     1                  - (T2*GNORM) * (ONE + HALF*T2)*GHINVG     2                  - HALF * (V(GTHG)*T1)**2      DO 80 I = 1, N 80      STEP(I) = T1*DIG(I) + T2*NWTSTP(I)C 999  RETURNC  ***  LAST CARD OF DDBDOG FOLLOWS  ***      END      DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)      saveC***BEGIN PROLOGUE  DDOTC***DATE WRITTEN   791001   (YYMMDD)C***REVISION DATE  820801   (YYMMDD)C***CATEGORY NO.  D1A4C***KEYWORDS  BLAS,DOUBLE PRECISION,INNER PRODUCT,LINEAR ALGEBRA,VECTORC***AUTHOR  LAWSON, C. L., (JPL)C           HANSON, R. J., (SNLA)C           KINCAID, D. R., (U. OF TEXAS)C           KROGH, F. T., (JPL)C***PURPOSE  D.P. inner product of d.p. vectorsC***DESCRIPTIONCC                B L A S  SubprogramC    Description of ParametersCC     --Input--C        N  number of elements in input vector(s)C       DX  double precision vector with N elementsC     INCX  storage spacing between elements of DXC       DY  double precision vector with N elementsC     INCY  storage spacing between elements of DYCC     --Output--C     DDOT  double precision dot product (zero if N .LE. 0)CC     Returns the dot product of double precision DX and DY.C     DDOT = sum for I = 0 to N-1 of  DX(LX+I*INCX) * DY(LY+I*INCY)C     where LX = 1 if INCX .GE. 0, else LX = (-INCX)*N, and LY isC     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 MATHEMATICALC                 SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323C***ROUTINES CALLED  (NONE)C***END PROLOGUE  DDOTC      DOUBLE PRECISION DX(1),DY(1)C***FIRST EXECUTABLE STATEMENT  DDOT      DDOT = 0.D0      IF(N.LE.0)RETURN      IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60    5 CONTINUECC         CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS.C      IX = 1      IY = 1      IF(INCX.LT.0)IX = (-N+1)*INCX + 1      IF(INCY.LT.0)IY = (-N+1)*INCY + 1      DO 10 I = 1,N         DDOT = DDOT + DX(IX)*DY(IY)        IX = IX + INCX        IY = IY + INCY   10 CONTINUE      RETURNCC        CODE FOR BOTH INCREMENTS EQUAL TO 1.CCC        CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 5.C   20 M = MOD(N,5)      IF( M .EQ. 0 ) GO TO 40      DO 30 I = 1,M         DDOT = DDOT + DX(I)*DY(I)   30 CONTINUE      IF( N .LT. 5 ) RETURN   40 MP1 = M + 1      DO 50 I = MP1,N,5         DDOT = DDOT + DX(I)*DY(I) + DX(I+1)*DY(I+1) +     1   DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4)   50 CONTINUE      RETURNCC         CODE FOR POSITIVE EQUAL INCREMENTS .NE.1.C   60 CONTINUE      NS = N*INCX          DO 70 I=1,NS,INCX          DDOT = DDOT + DX(I)*DY(I)   70     CONTINUE      RETURN      END      SUBROUTINE DITSUM(D, G, IV, LIV, LV, P, V, X)      saveCC  ***  PRINT ITERATION SUMMARY FOR ***SOL (VERSION 2.3)  ***CC  ***  PARAMETER DECLARATIONS  ***C      INTEGER LIV, LV, P      INTEGER IV(LIV)      DOUBLE PRECISION D(P), G(P), V(LV), X(P)CC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++CC  ***  LOCAL VARIABLES  ***C      INTEGER ALG, I, IV1, M, NF, NG, OL, PUC/6      REAL MODEL1(6), MODEL2(6)C/7C     CHARACTER*4 MODEL1(6), MODEL2(6)C/      DOUBLE PRECISION NRELDF, OLDF, PRELDF, RELDF, ZEROCC  ***  INTRINSIC FUNCTIONS  ***C/+      INTEGER IABS      DOUBLE PRECISION DABS, DMAX1C/C  ***  NO EXTERNAL FUNCTIONS OR SUBROUTINES  ***CC  ***  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, X0PRTCC  ***  IV SUBSCRIPT VALUES  ***CC/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/7C     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/CC  ***  V SUBSCRIPT VALUES  ***CC/6      DATA DSTNRM/2/, F/10/, F0/13/, FDIF/11/, NREDUC/6/, PREDUC/7/,     1     RELDX/17/, STPPAR/5/C/7C     PARAMETER (DSTNRM=2, F=10, F0=13, FDIF=11, NREDUC=6, PREDUC=7,C    1           RELDX=17, STPPAR=5)C/CC/6      DATA ZERO/0.D+0/C/7C     PARAMETER (ZERO=0.D+0)C/C/6      DATA MODEL1(1)/4H    /, MODEL1(2)/4H    /, MODEL1(3)/4H    /,     1     MODEL1(4)/4H    /, MODEL1(5)/4H  G /, MODEL1(6)/4H  S /,     2     MODEL2(1)/4H G  /, MODEL2(2)/4H S  /, MODEL2(3)/4HG-S /,     3     MODEL2(4)/4HS-G /, MODEL2(5)/4H-S-G/, MODEL2(6)/4H-G-S/C/7C     DATA MODEL1/'    ','    ','    ','    ','  G ','  S 
'/,C    1     MODEL2/' G  ',' S  ','G-S ','S-G ','-S-G','-G-S
































































































































































































































































































































































































































































































































































































































































































































'/C/CC-------------------------------  BODY  --------------------------------C      PU = IV(PRUNIT)      IF (PU .EQ. 0) GO TO 999      IV1 = IV(1)      IF (IV1 .GT. 62) IV1 = IV1 - 51      OL = IV(OUTLEV)      ALG = IV(ALGSAV)      IF (IV1 .LT. 2 .OR. IV1 .GT. 15) GO TO 370      IF (OL .EQ. 0) GO TO 120      IF (IV1 .GE. 12) GO TO 120      IF (IV1 .EQ. 2 .AND. IV(NITER) .EQ. 0) GO TO 390      IF (IV1 .GE. 10 .AND. IV(PRNTIT) .EQ. 0) GO TO 120      IF (IV1 .GT. 2) GO TO 10         IV(PRNTIT) = IV(PRNTIT) + 1         IF (IV(PRNTIT) .LT. IABS(OL)) GO TO 999 10   NF = IV(NFCALL) - IABS(IV(NFCOV))      IV(PRNTIT) = 0      RELDF = ZERO      PRELDF = ZERO      OLDF = DMAX1(DABS(V(F0)), DABS(V(F)))      IF (OLDF .LE. ZERO) GO TO 20         RELDF = V(FDIF) / OLDF         PRELDF = V(PREDUC) / OLDF 20   IF (OL .GT. 0) GO TO 60CC        ***  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 120C 50      WRITE(PU,110) IV(NITER), NF, V(F), RELDF, PRELDF, V(RELDX),     1                 V(STPPAR)         GO TO 120CC     ***  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 120C 90   WRITE(PU,110) IV(NITER), NF, V(F), RELDF, PRELDF,     1             V(RELDX), V(STPPAR), V(DSTNRM), NRELDF 100  FORMAT(I6,I5,D10.3,2D9.2,D8.1,A3,A4,2D8.1,D9.2) 110  FORMAT(I6,I5,D11.3,2D10.2,3D9.1,D10.2)C 120  IF (IV(STATPR) .LT. 0) GO TO 430      GO TO (999, 999, 130, 150, 170, 190, 210, 230, 250, 270, 290, 310,     1       330, 350, 520), IV1C 130  WRITE(PU,140) 140  FORMAT(/26H ***** X-CONVERGENCE *****)      GO TO 430C 150  WRITE(PU,160) 160  FORMAT(/42H ***** RELATIVE FUNCTION CONVERGENCE *****)      GO TO 430C 170  WRITE(PU,180) 180  FORMAT(/49H ***** X- AND RELATIVE FUNCTION CONVERGENCE *****)      GO TO 430C 190  WRITE(PU,200) 200  FORMAT(/42H ***** ABSOLUTE FUNCTION CONVERGENCE *****)      GO TO 430C 210  WRITE(PU,220) 220  FORMAT(/33H ***** SINGULAR CONVERGENCE *****)      GO TO 430C 230  WRITE(PU,240) 240  FORMAT(/30H ***** FALSE CONVERGENCE *****)      GO TO 430C 250  WRITE(PU,260) 260  FORMAT(/38H ***** FUNCTION EVALUATION LIMIT *****)      GO TO 430C 270  WRITE(PU,280) 280  FORMAT(/28H ***** ITERATION LIMIT *****)      GO TO 430C 290  WRITE(PU,300) 300  FORMAT(/18H ***** STOPX *****)      GO TO 430C 310  WRITE(PU,320) 320  FORMAT(/44H ***** INITIAL F(X) CANNOT BE COMPUTED *****)C      GO TO 390C 330  WRITE(PU,340) 340  FORMAT(/37H ***** BAD PARAMETERS TO ASSESS *****)      GO TO 999C 350  WRITE(PU,360) 360  FORMAT(/43H ***** GRADIENT COULD NOT BE COMPUTED *****)      IF (IV(NITER) .GT. 0) GO TO 480      GO TO 390C 370  WRITE(PU,380) IV(1) 380  FORMAT(/14H ***** IV(1) =,I5,6H *****)      GO TO 999CC  ***  INITIAL CALL ON DITSUM  ***C 390  IF (IV(X0PRT) .NE. 0) WRITE(PU,400) (I, X(I), D(I), I = 1, P) 400  FORMAT(/23H     I     INITIAL X(I),8X,4HD(I)//(1X,I5,D17.6,D14.3))      IF (IV1 .GE. 12) GO TO 999      IV(NEEDHD) = 0      IV(PRNTIT) = 0      IF (OL .EQ. 0) GO TO 999      IF (OL .LT. 0 .AND. ALG .EQ. 1) WRITE(PU,30)      IF (OL .LT. 0 .AND. ALG .EQ. 2) WRITE(PU,40)      IF (OL .GT. 0 .AND. ALG .EQ. 1) WRITE(PU,70)      IF (OL .GT. 0 .AND. ALG .EQ. 2) WRITE(PU,80)      IF (ALG .EQ. 1) WRITE(PU,410) V(F)      IF (ALG .EQ. 2) WRITE(PU,420) V(F) 410  FORMAT(/11H     0    1,D10.3)C365  FORMAT(/11H     0    1,E11.3) 420  FORMAT(/11H     0    1,D11.3)      GO TO 999CC  ***  PRINT VARIOUS INFORMATION REQUESTED ON SOLUTION  ***C 430  IV(NEEDHD) = 1      IF (IV(STATPR) .EQ. 0) GO TO 480         OLDF = DMAX1(DABS(V(F0)), DABS(V(F)))         PRELDF = ZERO         NRELDF = ZERO         IF (OLDF .LE. ZERO) GO TO 440              PRELDF = V(PREDUC) / OLDF              NRELDF = V(NREDUC) / OLDF 440     NF = IV(NFCALL) - IV(NFCOV)         NG = IV(NGCALL) - IV(NGCOV)         WRITE(PU,450) V(F), V(RELDX), NF, NG, PRELDF, NRELDF 450  FORMAT(/9H FUNCTION,D17.6,8H   RELDX,D17.3/12H FUNC. EVALS,     1   I8,9X,11HGRAD. EVALS,I8/7H PRELDF,D16.3,6X,7HNPRELDF,D15.3)C         IF (IV(NFCOV) .GT. 0) WRITE(PU,460) IV(NFCOV) 460     FORMAT(/1X,I4,50H EXTRA FUNC. EVALS FOR COVARIANCE AND DIAGNOST     1ICS.)         IF (IV(NGCOV) .GT. 0) WRITE(PU,470) IV(NGCOV) 470     FORMAT(1X,I4,50H EXTRA GRAD. EVALS FOR COVARIANCE AND DIAGNOSTI     1CS.)C 480  IF (IV(SOLPRT) .EQ. 0) GO TO 999         IV(NEEDHD) = 1         WRITE(PU,490) 490  FORMAT(/22H     I      FINAL X(I),8X,4HD(I),10X,4HG(I)/)         DO 500 I = 1, P              WRITE(PU,510) I, X(I), D(I), G(I) 500          CONTINUE 510     FORMAT(1X,I5,D16.6,2D14.3)      GO TO 999C 520  WRITE(PU,530) 530  FORMAT(/24H INCONSISTENT DIMENSIONS) 999  RETURNC  ***  LAST CARD OF DITSUM FOLLOWS  ***      END      SUBROUTINE DLITVM(N, X, L, Y)      saveCC  ***  SOLVE  (L**T)*X = Y,  WHERE  L  IS AN  N X N  LOWER TRIANGULARC  ***  MATRIX STORED COMPACTLY BY ROWS.  X AND Y MAY OCCUPY THE SAMEC  ***  STORAGE.  ***C      INTEGER N      DOUBLE PRECISION X(N), L(1), Y(N)      INTEGER I, II, IJ, IM1, I0, J, NP1      DOUBLE PRECISION XI, ZEROC/6      DATA ZERO/0.D+0/C/7C     PARAMETER (ZERO=0.D+0)C/C      DO 10 I = 1, N 10      X(I) = Y(I)      NP1 = N + 1      I0 = N*(N+1)/2      DO 30 II = 1, N         I = NP1 - II         XI = X(I)/L(I0)         X(I) = XI         IF (I .LE. 1) GO TO 999         I0 = I0 - I         IF (XI .EQ. ZERO) GO TO 30         IM1 = I - 1         DO 20 J = 1, IM1              IJ = I0 + J              X(J) = X(J) - XI*L(IJ) 20           CONTINUE 30      CONTINUE 999  RETURNC  ***  LAST CARD OF DLITVM FOLLOWS  ***      END      SUBROUTINE DLIVMU(N, X, L, Y)      saveCC  ***  SOLVE  L*X = Y, WHERE  L  IS AN  N X N  LOWER TRIANGULARC  ***  MATRIX STORED COMPACTLY BY ROWS.  X AND Y MAY OCCUPY THE SAMEC  ***  STORAGE.  ***C      INTEGER N      DOUBLE PRECISION X(N), L(1), Y(N)      DOUBLE PRECISION DDOT      INTEGER I, J, K      DOUBLE PRECISION T, ZEROC/6      DATA ZERO/0.D+0/C/7C     PARAMETER (ZERO=0.D+0)C/C      DO 10 K = 1, N         IF (Y(K) .NE. ZERO) GO TO 20         X(K) = ZERO 10      CONTINUE      GO TO 999 20   J = K*(K+1)/2      X(K) = Y(K) / L(J)      IF (K .GE. N) GO TO 999      K = K + 1      DO 30 I = K, N         T = DDOT(I-1, L(J+1),1,X,1)         J = J + I         X(I) = (Y(I) - T)/L(J) 30      CONTINUE 999  RETURNC  ***  LAST CARD OF DLIVMU FOLLOWS  ***      END      SUBROUTINE DLTVMU(N, X, L, Y)      saveCC  ***  COMPUTE  X = (L**T)*Y, WHERE  L  IS AN  N X N  LOWERC  ***  TRIANGULAR MATRIX STORED COMPACTLY BY ROWS.  X AND Y MAYC  ***  OCCUPY THE SAME STORAGE.  ***C      INTEGER N      DOUBLE PRECISION X(N), L(1), Y(N)C     DIMENSION L(N*(N+1)/2)      INTEGER I, IJ, I0, J      DOUBLE PRECISION YI, ZEROC/6      DATA ZERO/0.D+0/C/7C     PARAMETER (ZERO=0.D+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  RETURNC  ***  LAST CARD OF DLTVMU FOLLOWS  ***      END      SUBROUTINE DLUPDT(BETA, GAMMA, L, LAMBDA, LPLUS, N, W, Z)      saveCC  ***  COMPUTE LPLUS = SECANT UPDATE OF L  ***CC  ***  PARAMETER DECLARATIONS  ***C      INTEGER N      DOUBLE PRECISION BETA(N), GAMMA(N), L(1), LAMBDA(N), LPLUS(1),     1                 W(N), Z(N)C     DIMENSION L(N*(N+1)/2), LPLUS(N*(N+1)/2)CC--------------------------  PARAMETER USAGE  --------------------------CC   BETA = SCRATCH VECTOR.C  GAMMA = SCRATCH VECTOR.C      L (INPUT) LOWER TRIANGULAR MATRIX, STORED ROWWISE.C LAMBDA = SCRATCH VECTOR.C  LPLUS (OUTPUT) LOWER TRIANGULAR MATRIX, STORED ROWWISE, WHICH MAYC             OCCUPY THE SAME STORAGE AS  L.C      N (INPUT) LENGTH OF VECTOR PARAMETERS AND ORDER OF MATRICES.C      W (INPUT, DESTROYED ON OUTPUT) RIGHT SINGULAR VECTOR OF RANK 1C             CORRECTION TO  L.C      Z (INPUT, DESTROYED ON OUTPUT) LEFT SINGULAR VECTOR OF RANK 1C             CORRECTION TO  L.CC-------------------------------  NOTES  -------------------------------CC  ***  APPLICATION AND USAGE RESTRICTIONS  ***CC        THIS ROUTINE UPDATES THE CHOLESKY FACTOR  L  OF A SYMMETRICC     POSITIVE DEFINITE MATRIX TO WHICH A SECANT UPDATE IS BEINGC     APPLIED -- IT COMPUTES A CHOLESKY FACTOR  LPLUS  OFC     L * (I + Z*W**T) * (I + W*Z**T) * L**T.  IT IS ASSUMED THAT  WC     AND  Z  HAVE BEEN CHOSEN SO THAT THE UPDATED MATRIX IS STRICTLYC     POSITIVE DEFINITE.CC  ***  ALGORITHM NOTES  ***CC        THIS CODE USES RECURRENCE 3 OF REF. 1 (WITH D(J) = 1 FOR ALL J)C     TO COMPUTE  LPLUS  OF THE FORM  L * (I + Z*W**T) * Q,  WHERE  QC     IS AN ORTHOGONAL MATRIX THAT MAKES THE RESULT LOWER TRIANGULAR.C        LPLUS MAY HAVE SOME NEGATIVE DIAGONAL ELEMENTS.CC  ***  REFERENCES  ***CC 1.  GOLDFARB, D. (1976), FACTORIZED VARIABLE METRIC METHODS FOR UNCON-C             STRAINED OPTIMIZATION, MATH. COMPUT. 30, PP. 796-811.CC  ***  GENERAL  ***CC     CODED BY DAVID M. GAY (FALL 1979).C     THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH SUPPORTEDC     BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS MCS-7600324 ANDC     MCS-7906671.CC------------------------  EXTERNAL QUANTITIES  ------------------------CC  ***  INTRINSIC FUNCTIONS  ***C/+      DOUBLE PRECISION DSQRTC/C--------------------------  LOCAL VARIABLES  --------------------------C      INTEGER I, IJ, J, JJ, JP1, K, NM1, NP1      DOUBLE PRECISION A, B, BJ, ETA, GJ, LJ, LIJ, LJJ, NU, S, THETA,     1                 WJ, ZJ      DOUBLE PRECISION ONE, ZEROCC  ***  DATA INITIALIZATIONS  ***CC/6      DATA ONE/1.D+0/, ZERO/0.D+0/C/7C     PARAMETER (ONE=1.D+0, ZERO=0.D+0)C/CC+++++++++++++++++++++++++++++++  BODY  ++++++++++++++++++++++++++++++++C      NU = ONE      ETA = ZERO      IF (N .LE. 1) GO TO 30      NM1 = N - 1CC  ***  TEMPORARILY STORE S(J) = SUM OVER K = J+1 TO N OF W(K)**2 INC  ***  LAMBDA(J).C      S = ZERO      DO 10 I = 1, NM1         J = N - I         S = S + W(J+1)**2         LAMBDA(J) = S 10      CONTINUECC  ***  COMPUTE LAMBDA, GAMMA, AND BETA BY GOLDFARB*S RECURRENCE 3.C      DO 20 J = 1, NM1         WJ = W(J)         A = NU*Z(J) - ETA*WJ         THETA = ONE + A*WJ         S = A*LAMBDA(J)         LJ = DSQRT(THETA**2 + A*S)         IF (THETA .GT. ZERO) LJ = -LJ         LAMBDA(J) = LJ         B = THETA*WJ + S         GAMMA(J) = B * NU / LJ         BETA(J) = (A - B*ETA) / LJ         NU = -NU / LJ         ETA = -(ETA + (A**2)/(THETA - LJ)) / LJ 20      CONTINUE 30   LAMBDA(N) = ONE + (NU*Z(N) - ETA*W(N))*W(N)CC  ***  UPDATE L, GRADUALLY OVERWRITING  W  AND  Z  WITH  L*W  AND  L*Z.C      NP1 = N + 1      JJ = N * (N + 1) / 2      DO 60 K = 1, N         J = NP1 - K         LJ = LAMBDA(J)         LJJ = L(JJ)         LPLUS(JJ) = LJ * LJJ         WJ = W(J)         W(J) = LJJ * WJ         ZJ = Z(J)         Z(J) = LJJ * ZJ         IF (K .EQ. 1) GO TO 50         BJ = BETA(J)         GJ = GAMMA(J)         IJ = JJ + J         JP1 = J + 1         DO 40 I = JP1, N              LIJ = L(IJ)              LPLUS(IJ) = LJ*LIJ + BJ*W(I) + GJ*Z(I)              W(I) = W(I) + LIJ*WJ              Z(I) = Z(I) + LIJ*ZJ              IJ = IJ + I 40           CONTINUE 50      JJ = JJ - J 60      CONTINUEC 999  RETURNC  ***  LAST CARD OF DLUPDT FOLLOWS  ***      END      SUBROUTINE DLVMUL(N, X, L, Y)      saveCC  ***  COMPUTE  X = L*Y, WHERE  L  IS AN  N X N  LOWER TRIANGULARC  ***  MATRIX STORED COMPACTLY BY ROWS.  X AND Y MAY OCCUPY THE SAMEC  ***  STORAGE.  ***C      INTEGER N      DOUBLE PRECISION X(N), L(1), Y(N)C     DIMENSION L(N*(N+1)/2)      INTEGER I, II, IJ, I0, J, NP1      DOUBLE PRECISION T, ZEROC/6      DATA ZERO/0.D+0/C/7C     PARAMETER (ZERO=0.D+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  RETURNC  ***  LAST CARD OF DLVMUL FOLLOWS  ***      END      DOUBLE PRECISION FUNCTION DNRM2(N,DX,INCX)      saveC***BEGIN PROLOGUE  DNRM2C***DATE WRITTEN   791001   (YYMMDD)C***REVISION DATE  820801   (YYMMDD)C***CATEGORY NO.  D1A3BC***KEYWORDS  BLAS,DOUBLE PRECISION,EUCLIDEAN,L2,LENGTH,LINEAR ALGEBRA,C             NORM,VECTORC***AUTHOR  LAWSON, C. L., (JPL)C           HANSON, R. J., (SNLA)C           KINCAID, D. R., (U. OF TEXAS)C           KROGH, F. T., (JPL)C***PURPOSE  Euclidean length (L2 norm) of d.p. vectorC***DESCRIPTIONCC                B L A S  SubprogramC    Description of parametersCC     --Input--C        N  number of elements in input vector(s)C       DX  double precision vector with N elementsC     INCX  storage spacing between elements of DXCC     --Output--C    DNRM2  double precision result (zero if N .LE. 0)CC     Euclidean norm of the N-vector stored in DX() with storageC     increment INCX .C     If    N .LE. 0 return with result = 0.C     If N .GE. 1 then INCX must be .GE. 1CC           C.L. Lawson, 1978 Jan 08CC     Four phase method     using two built-in constants that areC     hopefully applicable to all machines.C         CUTLO = maximum of  DSQRT(U/EPS)  over all known machines.C         CUTHI = minimum of  DSQRT(V)      over all known machines.C     whereC         EPS = smallest no. such that EPS + 1. .GT. 1.C         U   = smallest positive no.   (underflow limit)C         V   = largest  no.            (overflow  limit)CC     Brief outline of algorithm..CC     Phase 1    scans zero components.C     move to phase 2 when a component is nonzero and .LE. CUTLOC     move to phase 3 when a component is .GT. CUTLOC     move to phase 4 when a component is .GE. CUTHI/MC     where M = N for X() real and M = 2*N for complex.CC     Values for CUTLO and CUTHI..C     From the environmental parameters listed in the IMSL converterC     document the limiting values are as followS..C     CUTLO, S.P.   U/EPS = 2**(-102) for  Honeywell.  Close seconds areC                   Univac and DEC at 2**(-103)C                   Thus CUTLO = 2**(-51) = 4.44089E-16C     CUTHI, S.P.   V = 2**127 for Univac, Honeywell, and DEC.C                   Thus CUTHI = 2**(63.5) = 1.30438E19C     CUTLO, D.P.   U/EPS = 2**(-67) for Honeywell and DEC.C                   Thus CUTLO = 2**(-33.5) = 8.23181D-11C     CUTHI, D.P.   same as S.P.  CUTHI = 1.30438D19C     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 MATHEMATICALC                 SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323C***ROUTINES CALLED  (NONE)C***END PROLOGUE  DNRM2      INTEGER          NEXT      DOUBLE PRECISION   DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE      DATA   ZERO, ONE /0.0D0, 1.0D0/C      DATA CUTLO, CUTHI / 8.232D-11,  1.304D19 /C***FIRST EXECUTABLE STATEMENT  DNRM2      IF(N .GT. 0) GO TO 10         DNRM2  = ZERO         GO TO 300C   10 ASSIGN 30 TO NEXT      SUM = ZERO      NN = N * INCXC                                                 BEGIN MAIN LOOP      I = 1   20    GO TO NEXT,(30, 50, 70, 110)   30 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85      ASSIGN 50 TO NEXT      XMAX = ZEROCC                        PHASE 1.  SUM IS ZEROC   50 IF( DX(I) .EQ. ZERO) GO TO 200      IF( DABS(DX(I)) .GT. CUTLO) GO TO 85CC                                PREPARE FOR PHASE 2.      ASSIGN 70 TO NEXT      GO TO 105CC                                PREPARE FOR PHASE 4.C  100 I = J      ASSIGN 110 TO NEXT      SUM = (SUM / DX(I)) / DX(I)  105 XMAX = DABS(DX(I))      GO TO 115CC                   PHASE 2.  SUM IS SMALL.C                             SCALE TO AVOID DESTRUCTIVE UNDERFLOW.C   70 IF( DABS(DX(I)) .GT. CUTLO ) GO TO 75CC                     COMMON CODE FOR PHASES 2 AND 4.C                     IN PHASE 4 SUM IS LARGE.  SCALE TO AVOID OVERFLOW.C  110 IF( DABS(DX(I)) .LE. XMAX ) GO TO 115         SUM = ONE + SUM * (XMAX / DX(I))**2         XMAX = DABS(DX(I))         GO TO 200C  115 SUM = SUM + (DX(I)/XMAX)**2      GO TO 200CCC                  PREPARE FOR PHASE 3.C   75 SUM = (SUM * XMAX) * XMAXCCC     FOR REAL OR D.P. SET HITEST = CUTHI/NC     FOR COMPLEX      SET HITEST = CUTHI/(2*N)C   85 HITEST = CUTHI/FLOAT( N )CC                   PHASE 3.  SUM IS MID-RANGE.  NO SCALING.C      DO 95 J =I,NN,INCX      IF(DABS(DX(J)) .GE. HITEST) GO TO 100   95    SUM = SUM + DX(J)**2      DNRM2 = DSQRT( SUM )      GO TO 300C  200 CONTINUE      I = I + INCX      IF ( I .LE. NN ) GO TO 20CC              END OF MAIN LOOP.CC              COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.C      DNRM2 = XMAX * DSQRT(SUM)  300 CONTINUE      RETURN      END      SUBROUTINE DPARCK(ALG, D, IV, LIV, LV, N, V)      saveCC  ***  CHECK ***SOL (VERSION 2.3) PARAMETERS, PRINT CHANGED VALUES  ***CC  ***  ALG = 1 FOR REGRESSION, ALG = 2 FOR GENERAL UNCONSTRAINED OPT.C      INTEGER ALG, LIV, LV, N      INTEGER IV(LIV)      DOUBLE PRECISION D(N), V(LV)C      EXTERNAL  DVDFLT      DOUBLE PRECISION D1MACHC DVDFLT  -- SUPPLIES DEFAULT PARAMETER VALUES TO V ALONE.C/+      INTEGER MAX0C/CC  ***  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/7C     CHARACTER*1 VARNM(2), SH(2)C     CHARACTER*4 CNGD(3), DFLT(3), VN(2,34), WHICH(3)C/      DOUBLE PRECISION BIG, MACHEP, TINY, VK, VM(34), VX(34), ZEROCC  ***  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, VNEEDCCC/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/7C     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, TINYC/C      DATA BIG/0.D+0/, MACHEP/-1.D+0/, TINY/1.D+0/, ZERO/0.D+0/C/6      DATA VN(1,1),VN(2,1)/4HEPSL,4HON../      DATA VN(1,2),VN(2,2)/4HPHMN,4HFC../      DATA VN(1,3),VN(2,3)/4HPHMX,4HFC../      DATA VN(1,4),VN(2,4)/4HDECF,4HAC../      DATA VN(1,5),VN(2,5)/4HINCF,4HAC../      DATA VN(1,6),VN(2,6)/4HRDFC,4HMN../      DATA VN(1,7),VN(2,7)/4HRDFC,4HMX../      DATA VN(1,8),VN(2,8)/4HTUNE,4HR1../      DATA VN(1,9),VN(2,9)/4HTUNE,4HR2../      DATA VN(1,10),VN(2,10)/4HTUNE,4HR3../      DATA VN(1,11),VN(2,11)/4HTUNE,4HR4../      DATA VN(1,12),VN(2,12)/4HTUNE,4HR5../      DATA VN(1,13),VN(2,13)/4HAFCT,4HOL../      DATA VN(1,14),VN(2,14)/4HRFCT,4HOL../      DATA VN(1,15),VN(2,15)/4HXCTO,4HL.../      DATA VN(1,16),VN(2,16)/4HXFTO,4HL.../      DATA VN(1,17),VN(2,17)/4HLMAX,4H0.../      DATA VN(1,18),VN(2,18)/4HLMAX,4HS.../      DATA VN(1,19),VN(2,19)/4HSCTO,4HL.../      DATA VN(1,20),VN(2,20)/4HDINI,4HT.../      DATA VN(1,21),VN(2,21)/4HDTIN,4HIT../      DATA VN(1,22),VN(2,22)/4HD0IN,4HIT../      DATA VN(1,23),VN(2,23)/4HDFAC,4H..../      DATA VN(1,24),VN(2,24)/4HDLTF,4HDC../      DATA VN(1,25),VN(2,25)/4HDLTF,4HDJ../      DATA VN(1,26),VN(2,26)/4HDELT,4HA0../      DATA VN(1,27),VN(2,27)/4HFUZZ,4H..../      DATA VN(1,28),VN(2,28)/4HRLIM,4HIT../      DATA VN(1,29),VN(2,29)/4HCOSM,4HIN../      DATA VN(1,30),VN(2,30)/4HHUBE,4HRC../      DATA VN(1,31),VN(2,31)/4HRSPT,4HOL../      DATA VN(1,32),VN(2,32)/4HSIGM,4HIN../      DATA VN(1,33),VN(2,33)/4HETA0,4H..../      DATA VN(1,34),VN(2,34)/4HBIAS,4H..../C/7C     DATA VN(1,1),VN(2,1)/'EPSL','ON..
'/C     DATA VN(1,2),VN(2,2)/'PHMN','FC..
'/C     DATA VN(1,3),VN(2,3)/'PHMX','FC..
'/C     DATA VN(1,4),VN(2,4)/'DECF','AC..
'/C     DATA VN(1,5),VN(2,5)/'INCF','AC..
'/C     DATA VN(1,6),VN(2,6)/'RDFC','MN..
'/C     DATA VN(1,7),VN(2,7)/'RDFC','MX..
'/C     DATA VN(1,8),VN(2,8)/'TUNE','R1..
'/C     DATA VN(1,9),VN(2,9)/'TUNE','R2..
'/C     DATA VN(1,10),VN(2,10)/'TUNE','R3..
'/C     DATA VN(1,11),VN(2,11)/'TUNE','R4..
'/C     DATA VN(1,12),VN(2,12)/'TUNE','R5..
'/C     DATA VN(1,13),VN(2,13)/'AFCT','OL..
'/C     DATA VN(1,14),VN(2,14)/'RFCT','OL..
'/C     DATA VN(1,15),VN(2,15)/'XCTO','L...
'/C     DATA VN(1,16),VN(2,16)/'XFTO','L...
'/C     DATA VN(1,17),VN(2,17)/'LMAX','0...
'/C     DATA VN(1,18),VN(2,18)/'LMAX','S...
'/C     DATA VN(1,19),VN(2,19)/'SCTO','L...
'/C     DATA VN(1,20),VN(2,20)/'DINI','T...
'/C     DATA VN(1,21),VN(2,21)/'DTIN','IT..
'/C     DATA VN(1,22),VN(2,22)/'D0IN','IT..
'/C     DATA VN(1,23),VN(2,23)/'DFAC','....
'/C     DATA VN(1,24),VN(2,24)/'DLTF','DC..
'/C     DATA VN(1,25),VN(2,25)/'DLTF','DJ..
'/C     DATA VN(1,26),VN(2,26)/'DELT','A0..
'/C     DATA VN(1,27),VN(2,27)/'FUZZ','....
'/C     DATA VN(1,28),VN(2,28)/'RLIM','IT..
'/C     DATA VN(1,29),VN(2,29)/'COSM','IN..
'/C     DATA VN(1,30),VN(2,30)/'HUBE','RC..
'/C     DATA VN(1,31),VN(2,31)/'RSPT','OL..
'/C     DATA VN(1,32),VN(2,32)/'SIGM','IN..
'/C     DATA VN(1,33),VN(2,33)/'ETA0','....
'/C     DATA VN(1,34),VN(2,34)/'BIAS','....






















'/C/C      DATA VM(1)/1.0D-3/, VM(2)/-0.99D+0/, VM(3)/1.0D-3/, VM(4)/1.0D-2/,     1     VM(5)/1.2D+0/, VM(6)/1.D-2/, VM(7)/1.2D+0/, VM(8)/0.D+0/,     2     VM(9)/0.D+0/, VM(10)/1.D-3/, VM(11)/-1.D+0/, VM(15)/0.D+0/,     3     VM(16)/0.D+0/, VM(19)/0.D+0/, VM(20)/-10.D+0/, VM(21)/0.D+0/,     4     VM(22)/0.D+0/, VM(23)/0.D+0/, VM(27)/1.01D+0/,     5     VM(28)/1.D+10/, VM(30)/0.D+0/, VM(31)/0.D+0/, VM(32)/0.D+0/,     6     VM(34)/0.D+0/      DATA VX(1)/0.9D+0/, VX(2)/-1.D-3/, VX(3)/1.D+1/, VX(4)/0.8D+0/,     1     VX(5)/1.D+2/, VX(6)/0.8D+0/, VX(7)/1.D+2/, VX(8)/0.5D+0/,     2     VX(9)/0.5D+0/, VX(10)/1.D+0/, VX(11)/1.D+0/, VX(14)/0.1D+0/,     3     VX(15)/1.D+0/, VX(16)/1.D+0/, VX(19)/1.D+0/, VX(23)/1.D+0/,     4     VX(24)/1.D+0/, VX(25)/1.D+0/, VX(26)/1.D+0/, VX(27)/1.D+10/,     5     VX(29)/1.D+0/, VX(31)/1.D+0/, VX(32)/1.D+0/, VX(33)/1.D+0/,     6     VX(34)/1.D+0/CC/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/7C     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/CC...............................  BODY  ................................C      IF (ALG .LT. 1 .OR. ALG .GT. 2) GO TO 330      IF (IV(1) .EQ. 0) CALL DDEFLT(ALG, IV, LIV, LV, V)      PU = IV(PRUNIT)      MIV1 = MINIV(ALG)      IF (PERM .LE. LIV) MIV1 = MAX0(MIV1, IV(PERM) - 1)      IF (IVNEED .LE. LIV) MIV2 = MIV1 + MAX0(IV(IVNEED), 0)      IF (LASTIV .LE. LIV) IV(LASTIV) = MIV2      IF (LIV .LT. MIV1) GO TO 290      IV(IVNEED) = 0      IV(LASTV) = MAX0(IV(VNEED), 0) + IV(LMAT) - 1      IF (LIV .LT. MIV2) GO TO 290      IF (LV .LT. IV(LASTV)) GO TO 310      IV(VNEED) = 0      IF (ALG .EQ. IV(ALGSAV)) GO TO 20         IF (PU .NE. 0) WRITE(PU,10) ALG, IV(ALGSAV) 10      FORMAT(/' THE FIRST PARAMETER TO DDEFLT SHOULD BE




























































































































































































































































































































',I3,     1          12H RATHER THAN,I3)         IV(1) = 82         GO TO 999 20   IV1 = IV(1)      IF (IV1 .LT. 12 .OR. IV1 .GT. 14) GO TO 50         IF (N .GE. 1) GO TO 40              IV(1) = 81              IF (PU .EQ. 0) GO TO 999              WRITE(PU,30) VARNM(ALG), N 30           FORMAT(/8H /// BAD,A1,2H =,I5)              GO TO 999 40      IF (IV1 .NE. 14) IV(NEXTIV) = IV(PERM)         IF (IV1 .NE. 14) IV(NEXTV) = IV(LMAT)         IF (IV1 .EQ. 13) GO TO 999         K = IV(PARSAV) - EPSLON         CALL DVDFLT(ALG, LV-K, V(K+1))         IV(DTYPE0) = 2 - ALG         IV(OLDN) = N         WHICH(1) = DFLT(1)         WHICH(2) = DFLT(2)         WHICH(3) = DFLT(3)         GO TO 100 50   IF (N .EQ. IV(OLDN)) GO TO 70         IV(1) = 17         IF (PU .EQ. 0) GO TO 999         WRITE(PU,60) VARNM(ALG), IV(OLDN), N 60      FORMAT(/5H /// ,1A1,14H CHANGED FROM ,I5,4H TO ,I5)         GO TO 999C 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 999C 90   WHICH(1) = CNGD(1)      WHICH(2) = CNGD(2)      WHICH(3) = CNGD(3)C 100  IF (IV1 .EQ. 14) IV1 = 12      IF (BIG .GT. TINY) GO TO 110         TINY = D1MACH(1)         MACHEP = D1MACH(4)         BIG = D1MACH(2)         VM(12) = MACHEP         VX(12) = BIG         VM(13) = TINY         VX(13) = BIG         VM(14) = MACHEP         VM(17) = TINY         VX(17) = BIG         VM(18) = TINY         VX(18) = BIG         VX(20) = BIG         VX(21) = BIG         VX(22) = BIG         VM(24) = MACHEP         VM(25) = MACHEP         VM(26) = MACHEP         VX(28) = DSQRT(D1MACH(2))*16.         VM(29) = MACHEP         VX(30) = BIG         VM(33) = MACHEP 110  M = 0      I = 1      J = JLIM(ALG)      K = EPSLON      NDFALT = NDFLT(ALG)      DO 140 L = 1, NDFALT         VK = V(K)         IF (VK .GE. VM(I) .AND. VK .LE. VX(I)) GO TO 130              M = K              IF (PU .NE. 0) WRITE(PU,120) VN(1,I), VN(2,I), K, VK,     1                                    VM(I), VX(I) 120          FORMAT(/6H ///  ,2A4,5H.. V(,I2,3H) =,D11.3,7H SHOULD,     1               11H BE BETWEEN,D11.3,4H AND,D11.3) 130     K = K + 1         I = I + 1         IF (I .EQ. J) I = IJMP 140     CONTINUEC      IF (IV(NVDFLT) .EQ. NDFALT) GO TO 160         IV(1) = 51         IF (PU .EQ. 0) GO TO 999         WRITE(PU,150) IV(NVDFLT), NDFALT 150     FORMAT(/13H IV(NVDFLT) =,I5,13H RATHER THAN ,I5)         GO TO 999 160  IF ((IV(DTYPE) .GT. 0 .OR. V(DINIT) .GT. ZERO) .AND. IV1 .EQ. 12)     1                  GO TO 190      DO 180 I = 1, N         IF (D(I) .GT. ZERO) GO TO 180              M = 18              IF (PU .NE. 0) WRITE(PU,170) I, D(I) 170     FORMAT(/8H ///  D(,I3,3H) =,D11.3,19H SHOULD BE POSITIVE) 180     CONTINUE 190  IF (M .EQ. 0) GO TO 200         IV(1) = M         GO TO 999C 200  IF (PU .EQ. 0 .OR. IV(PARPRT) .EQ. 0) GO TO 999      IF (IV1 .NE. 12 .OR. IV(INITS) .EQ. ALG-1) GO TO 220         M = 1         WRITE(PU,210) SH(ALG), IV(INITS) 210     FORMAT(/22H NONDEFAULT VALUES..../5H INIT,A1,14H..... IV(25) =,     1          I3) 220  IF (IV(DTYPE) .EQ. IV(DTYPE0)) GO TO 240         IF (M .EQ. 0) WRITE(PU,250) WHICH         M = 1         WRITE(PU,230) IV(DTYPE) 230     FORMAT(20H DTYPE..... IV(16) =,I3) 240  I = 1      J = JLIM(ALG)      K = EPSLON      L = IV(PARSAV)      NDFALT = NDFLT(ALG)      DO 280 II = 1, NDFALT         IF (V(K) .EQ. V(L)) GO TO 270              IF (M .EQ. 0) WRITE(PU,250) WHICH 250          FORMAT(/1H ,3A4,9HALUES..../)              M = 1              WRITE(PU,260) VN(1,I), VN(2,I), K, V(K) 260          FORMAT(1X,2A4,5H.. V(,I2,3H) =,D15.7) 270     K = K + 1         L = L + 1         I = I + 1         IF (I .EQ. J) I = IJMP 280     CONTINUEC      IV(DTYPE0) = IV(DTYPE)      PARSV1 = IV(PARSAV)      CALL DCOPY(IV(NVDFLT), V(EPSLON),1,V(PARSV1),1)      GO TO 999C 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 999C 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 999C 330  IV(1) = 67C 999  RETURNC  ***  LAST CARD OF DPARCK FOLLOWS  ***      END      DOUBLE PRECISION FUNCTION DRELST(P, D, X, X0)      saveCC  ***  COMPUTE AND RETURN RELATIVE DIFFERENCE BETWEEN X AND X0  ***C  ***  NL2SOL VERSION 2.2  ***C      INTEGER P      DOUBLE PRECISION D(P), X(P), X0(P)C/+      DOUBLE PRECISION DABSC/      INTEGER I      DOUBLE PRECISION EMAX, T, XMAX, ZEROC/6      DATA ZERO/0.D+0/C/7C     PARAMETER (ZERO=0.D+0)C/C      EMAX = ZERO      XMAX = ZERO      DO 10 I = 1, P         T = DABS(D(I) * (X(I) - X0(I)))         IF (EMAX .LT. T) EMAX = T         T = D(I) * (DABS(X(I)) + DABS(X0(I)))         IF (XMAX .LT. T) XMAX = T 10      CONTINUE      DRELST = ZERO      IF (XMAX .GT. ZERO) DRELST = EMAX / XMAX 999  RETURNC  ***  LAST CARD OF DRELST FOLLOWS  ***      END      LOGICAL FUNCTION DSTOPX(IDUMMY)      saveC     *****PARAMETERS...      INTEGER IDUMMYCC     ..................................................................CC     *****PURPOSE...C     THIS FUNCTION MAY SERVE AS THE DSTOPX (ASYNCHRONOUS INTERRUPTION)C     FUNCTION FOR THE NL2SOL (NONLINEAR LEAST-SQUARES) PACKAGE ATC     THOSE INSTALLATIONS WHICH DO NOT WISH TO IMPLEMENT AC     DYNAMIC DSTOPX.CC     *****ALGORITHM NOTES...C     AT INSTALLATIONS WHERE THE NL2SOL SYSTEM IS USEDC     INTERACTIVELY, THIS DUMMY DSTOPX SHOULD BE REPLACED BY AC     FUNCTION THAT RETURNS .TRUE. IF AND ONLY IF THE INTERRUPTC     (BREAK) KEY HAS BEEN PRESSED SINCE THE LAST CALL ON DSTOPX.CC     ..................................................................C      DSTOPX = .FALSE.      RETURN      END      SUBROUTINE FDUMP      saveC***BEGIN PROLOGUE  FDUMPC***DATE WRITTEN   790801   (YYMMDD)C***REVISION DATE  820801   (YYMMDD)C***CATEGORY NO.  ZC***KEYWORDS  ERROR,XERROR PACKAGEC***AUTHOR  JONES, R. E., (SNLA)C***PURPOSE  Symbolic dump (should be locally written).C***DESCRIPTIONC        ***Note*** Machine Dependent RoutineC        FDUMP is intended to be replaced by a locally writtenC        version which produces a symbolic dump.  Failing this,C        it should be replaced by a version which prints theC        subprogram nesting list.  Note that this dump must beC        printed on each of up to five files, as indicated by theC        XGETUA routine.  See XSETUA and XGETUA for details.CC     Written by Ron Jones, with SLATEC Common Math Library SubcommitteeC     Latest revision ---  23 May 1979C***ROUTINES CALLED  (NONE)C***END PROLOGUE  FDUMPC***FIRST EXECUTABLE STATEMENT  FDUMP      RETURN      END      FUNCTION J4SAVE(IWHICH,IVALUE,ISET)      saveC***BEGIN PROLOGUE  J4SAVEC***REFER TO  XERRORC     AbstractC        J4SAVE saves and recalls several global variables neededC        by the library error handling routines.CC     Description of ParametersC      --Input--C        IWHICH - Index of item desired.C                = 1 Refers to current error number.C                = 2 Refers to current error control flag.C                 = 3 Refers to current unit number to which errorC                    messages are to be sent.  (0 means use standard.)C                 = 4 Refers to the maximum number of times anyC                     message is to be printed (as set by XERMAX).C                 = 5 Refers to the total number of units to whichC                     each error message is to be written.C                 = 6 Refers to the 2nd unit for error messagesC                 = 7 Refers to the 3rd unit for error messagesC                 = 8 Refers to the 4th unit for error messagesC                 = 9 Refers to the 5th unit for error messagesC        IVALUE - The value to be set for the IWHICH-th parameter,C                 if ISET is .TRUE. .C        ISET   - If ISET=.TRUE., the IWHICH-th parameter will BEC                 given the value, IVALUE.  If ISET=.FALSE., theC                 IWHICH-th parameter will be unchanged, and IVALUEC                 is a dummy parameter.C      --Output--C        The (old) value of the IWHICH-th parameter will be returnedC        in the function value, J4SAVE.CC     Written by Ron Jones, with SLATEC Common Math Library SubcommitteeC    Adapted from Bell Laboratories PORT Library Error HandlerC     Latest revision ---  23 MAY 1979C***REFERENCES  JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR-C                 HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES,C                 1982.C***ROUTINES CALLED  (NONE)C***END PROLOGUE  J4SAVE      LOGICAL ISET      INTEGER IPARAM(9)      SAVE IPARAM      DATA IPARAM(1),IPARAM(2),IPARAM(3),IPARAM(4)/0,2,0,10/      DATA IPARAM(5)/1/      DATA IPARAM(6),IPARAM(7),IPARAM(8),IPARAM(9)/0,0,0,0/C***FIRST EXECUTABLE STATEMENT  J4SAVE      J4SAVE = IPARAM(IWHICH)      IF (ISET) IPARAM(IWHICH) = IVALUE      RETURN      END      SUBROUTINE XERABT(MESSG,NMESSG)      saveC***BEGIN PROLOGUE  XERABTC***DATE WRITTEN   790801   (YYMMDD)C***REVISION DATE  820801   (YYMMDD)C***CATEGORY NO.  R3CC***KEYWORDS  ERROR,XERROR PACKAGEC***AUTHOR  JONES, R. E., (SNLA)C***PURPOSE  Aborts program execution and prints error message.C***DESCRIPTIONC     AbstractC        ***Note*** machine dependent routineC        XERABT aborts the execution of the program.C        The error message causing the abort is given in the callingC        sequence, in case one needs it for printing on a dayfile,C        for example.CC     Description of ParametersC        MESSG and NMESSG are as in XERROR, except that NMESSG mayC        be zero, in which case no message is being supplied.CC     Written by Ron Jones, with SLATEC Common Math Library SubcommitteeC     Latest revision ---  19 MAR 1980C***REFERENCES  JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR-C                 HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES,C                 1982.C***ROUTINES CALLED  (NONE)C***END PROLOGUE  XERABT      CHARACTER*(*) MESSGC      EXTERNAL ERRORC***FIRST EXECUTABLE STATEMENT  XERABTC      CALL ERROR('dsumsl(), dsmsno() aborted\n















































































')      END      SUBROUTINE XERCTL(MESSG1,NMESSG,NERR,LEVEL,KONTRL)      saveC***BEGIN PROLOGUE  XERCTLC***DATE WRITTEN   790801   (YYMMDD)C***REVISION DATE  820801   (YYMMDD)C***CATEGORY NO.  R3CC***KEYWORDS  ERROR,XERROR PACKAGEC***AUTHOR  JONES, R. E., (SNLA)C***PURPOSE  Allows user control over handling of individual errors.C***DESCRIPTIONC     AbstractC        Allows user control over handling of individual errors.C        Just after each message is recorded, but before it isC        processed any further (i.e., before it is printed orC        a decision to abort is made), a call is made to XERCTL.C        If the user has provided his own version of XERCTL, heC        can then override the value of KONTROL used in processingC        this message by redefining its value.C        KONTRL may be set to any value from -2 to 2.C        The meanings for KONTRL are the same as in XSETF, exceptC        that the value of KONTRL changes only for this message.C        If KONTRL is set to a value outside the range from -2 to 2,C        it will be moved back into that range.CC     Description of ParametersCC      --Input--C        MESSG1 - the first word (only) of the error message.C        NMESSG - same as in the call to XERROR or XERRWV.C        NERR   - same as in the call to XERROR or XERRWV.C        LEVEL  - same as in the call to XERROR or XERRWV.C        KONTRL - the current value of the control flag as setC                 by a call to XSETF.CC      --Output--C        KONTRL - the new value of KONTRL.  If KONTRL is notC                 defined, it will remain at its original value.C                 This changed value of control affects onlyC                 the current occurrence of the current message.C***REFERENCES  JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR-C                 HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES,C                 1982.C***ROUTINES CALLED  (NONE)C***END PROLOGUE  XERCTL      CHARACTER*20 MESSG1C***FIRST EXECUTABLE STATEMENT  XERCTL      RETURN      END      SUBROUTINE XERPRT(MESSG,NMESSG)      saveC***BEGIN PROLOGUE  XERPRTC***DATE WRITTEN   790801   (YYMMDD)C***REVISION DATE  820801   (YYMMDD)C***CATEGORY NO.  ZC***KEYWORDS  ERROR,XERROR PACKAGEC***AUTHOR  JONES, R. E., (SNLA)C***PURPOSE  Prints error messages.C***DESCRIPTIONC     AbstractC        Print the Hollerith message in MESSG, of length NMESSG,C        on each file indicated by XGETUA.C     Latest revision ---  19 MAR 1980C***REFERENCES  JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR-C                 HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES,C                 1982.C***ROUTINES CALLED  I1MACH,S88FMT,XGETUAC***END PROLOGUE  XERPRT      INTEGER LUN(5)      CHARACTER*(*) MESSGC     OBTAIN UNIT NUMBERS AND WRITE LINE TO EACH UNITC***FIRST EXECUTABLE STATEMENT  XERPRT      CALL XGETUA(LUN,NUNIT)      LENMES = LEN(MESSG)      DO 20 KUNIT=1,NUNIT         IUNIT = LUN(KUNIT)         IF (IUNIT.EQ.0) IUNIT = I1MACH(4)         DO 10 ICHAR=1,LENMES,72            LAST = MIN0(ICHAR+71 , LENMES)            WRITE (IUNIT,'(1X,A)































































































































































































































































































































































































































































































Generated by  Doxygen 1.6.0   Back to index