https://github.com/cran/tseries
Tip revision: 54980448e6652a75784de03888c21e6b3ae0a3df authored by Compiled by Adrian Trapletti on 02 March 2000, 00:00:00 UTC
version 0.5-2
version 0.5-2
Tip revision: 5498044
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 /
C
C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR).
C
C === MACHINE = PDP-10.KI
C DATA IMACH( 1) / 5 /
C DATA IMACH( 2) / 6 /
C DATA IMACH( 3) / 7 /
C DATA IMACH( 4) / 6 /
C DATA IMACH( 5) / 36 /
C DATA IMACH( 6) / 5 /
C DATA IMACH( 7) / 2 /
C DATA IMACH( 8) / 35 /
C DATA IMACH( 9) / "377777777777 /
C DATA IMACH(10) / 2 /
C DATA IMACH(11) / 27 /
C DATA IMACH(12) / -128 /
C DATA IMACH(13) / 127 /
C DATA IMACH(14) / 62 /
C DATA IMACH(15) / -128 /
C DATA IMACH(16) / 127 /
C
C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING
C 32-BIT INTEGER ARITHMETIC.
C
C === MACHINE = PDP-11.32-BIT
C DATA IMACH( 1) / 5 /
C DATA IMACH( 2) / 6 /
C DATA IMACH( 3) / 7 /
C DATA IMACH( 4) / 6 /
C DATA IMACH( 5) / 32 /
C DATA IMACH( 6) / 4 /
C DATA IMACH( 7) / 2 /
C DATA IMACH( 8) / 31 /
C DATA IMACH( 9) / 2147483647 /
C DATA IMACH(10) / 2 /
C DATA IMACH(11) / 24 /
C DATA IMACH(12) / -127 /
C DATA IMACH(13) / 127 /
C DATA IMACH(14) / 56 /
C DATA IMACH(15) / -127 /
C DATA IMACH(16) / 127 /
C
C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING
C 16-BIT INTEGER ARITHMETIC.
C
C === MACHINE = PDP-11.16-BIT
C DATA IMACH( 1) / 5 /
C DATA IMACH( 2) / 6 /
C DATA IMACH( 3) / 7 /
C DATA IMACH( 4) / 6 /
C DATA IMACH( 5) / 16 /
C DATA IMACH( 6) / 2 /
C DATA IMACH( 7) / 2 /
C DATA IMACH( 8) / 15 /
C DATA IMACH( 9) / 32767 /
C DATA IMACH(10) / 2 /
C DATA IMACH(11) / 24 /
C DATA IMACH(12) / -127 /
C DATA IMACH(13) / 127 /
C DATA IMACH(14) / 56 /
C DATA IMACH(15) / -127 /
C DATA IMACH(16) / 127 /
C
C MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000.
C
C === MACHINE = SEQUENT.BALANCE.8000
C DATA IMACH( 1) / 0 /
C DATA IMACH( 2) / 0 /
C DATA IMACH( 3) / 7 /
C DATA IMACH( 4) / 0 /
C DATA IMACH( 5) / 32 /
C DATA IMACH( 6) / 1 /
C DATA IMACH( 7) / 2 /
C DATA IMACH( 8) / 31 /
C DATA IMACH( 9) / 2147483647 /
C DATA IMACH(10) / 2 /
C DATA IMACH(11) / 24 /
C DATA IMACH(12) / -125 /
C DATA IMACH(13) / 128 /
C DATA IMACH(14) / 53 /
C DATA IMACH(15) / -1021 /
C DATA IMACH(16) / 1024 /
C
C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILER
C
C === MACHINE = UNIVAC.1100
C DATA IMACH( 1) / 5 /
C DATA IMACH( 2) / 6 /
C DATA IMACH( 3) / 1 /
C DATA IMACH( 4) / 6 /
C DATA IMACH( 5) / 36 /
C DATA IMACH( 6) / 4 /
C DATA IMACH( 7) / 2 /
C DATA IMACH( 8) / 35 /
C DATA IMACH( 9) / O377777777777 /
C DATA IMACH(10) / 2 /
C DATA IMACH(11) / 27 /
C DATA IMACH(12) / -128 /
C DATA IMACH(13) / 127 /
C DATA IMACH(14) / 60 /
C DATA IMACH(15) /-1024 /
C DATA IMACH(16) / 1023 /
C
C MACHINE CONSTANTS FOR THE VAX 11/780
C
C === MACHINE = VAX.11/780
C DATA IMACH(1) / 5 /
C DATA IMACH(2) / 6 /
C DATA IMACH(3) / 5 /
C DATA IMACH(4) / 6 /
C DATA IMACH(5) / 32 /
C DATA IMACH(6) / 4 /
C DATA IMACH(7) / 2 /
C DATA IMACH(8) / 31 /
C DATA IMACH(9) /2147483647 /
C DATA IMACH(10)/ 2 /
C DATA IMACH(11)/ 24 /
C DATA IMACH(12)/ -127 /
C DATA IMACH(13)/ 127 /
C DATA IMACH(14)/ 56 /
C DATA IMACH(15)/ -127 /
C DATA IMACH(16)/ 127 /
C
C
C***FIRST EXECUTABLE STATEMENT I1MACH
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 F
C 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)
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 FDUMP,I1MACH,J4SAVE,XERABT,XERCTL,XERPRT,XERSAV,
C XGETUA
C***END PROLOGUE XERRWV
CHARACTER*(*) MESSG
CHARACTER*20 LFIRST
CHARACTER*37 FORM
DIMENSION LUN(5)
C GET FLAGS
C***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 CONTINUE
C 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 CONTINUE
C 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 40
C ERROR NUMBER
WRITE (IUNIT,30) LERR
30 FORMAT (15H ERROR NUMBER =,I10)
40 CONTINUE
50 CONTINUE
C 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 = 1
C QUIT HERE IF MESSAGE IS NOT FATAL
IF (IFATAL.LE.0) RETURN
IF ((LKNTRL.LE.0).OR.(KOUNT.GT.MAX0(1,MAXMES))) GO TO 120
C 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 CONTINUE
C 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)
save
C***BEGIN PROLOGUE XERSAV
C***DATE WRITTEN 800319 (YYMMDD)
C***REVISION DATE 820801 (YYMMDD)
C***CATEGORY NO. Z
C***KEYWORDS ERROR,XERROR PACKAGE
C***AUTHOR JONES, R. E., (SNLA)
C***PURPOSE Records that an error occurred.
C***DESCRIPTION
C Abstract
C Record that this error occurred.
C
C Description of Parameters
C --Input--
C MESSG, NMESSG, NERR, LEVEL are as in XERROR,
C except that when NMESSG=0 the tables will be
C dumped and cleared, and when NMESSG is less than zero the
C tables will be dumped and not cleared.
C --Output--
C ICOUNT will be the number of times this message has
C been seen, or zero if the table has overflowed and
C does not contain this message specifically.
C When NMESSG=0, ICOUNT will not be altered.
C
C Written by Ron Jones, with SLATEC Common Math Library Subcommittee
C Latest revision --- 19 Mar 1980
C***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,XGETUA
C***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,KOUNTX
C NEXT TWO DATA STATEMENTS ARE NECESSARY TO PROVIDE A BLANK
C 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 80
C DUMP THE TABLE
IF (KOUNT(1).EQ.0) RETURN
C 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 CONTINUE
C 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) RETURN
C CLEAR THE ERROR TABLES
DO 70 I=1,10
70 KOUNT(I) = 0
KOUNTX = 0
RETURN
80 CONTINUE
C 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 CONTINUE
C THREE POSSIBLE CASES...
C TABLE IS FULL
KOUNTX = KOUNTX+1
ICOUNT = 1
RETURN
C MESSAGE FOUND IN TABLE
100 KOUNT(II) = KOUNT(II) + 1
ICOUNT = KOUNT(II)
RETURN
C 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)
save
C***BEGIN PROLOGUE XGETUA
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 Returns unit number(s) to which error messages are being
C sent.
C***DESCRIPTION
C Abstract
C XGETUA may be called to determine the unit number or numbers
C 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.
C
C Description of Parameters
C --Output--
C IUNIT - an array of one to five unit numbers, depending
C on the value of N. A value of zero refers to the
C default unit, as defined by the I1MACH machine
C constant routine. Only IUNIT(1),...,IUNIT(N) are
C defined by XGETUA. The values of IUNIT(N+1),...,
C IUNIT(5) are not defined (for N .LT. 5) or altered
C in any way by XGETUA.
C N - the number of units to which copies of the
C error messages are being sent. N will be in the
C range from 1 to 5.
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 J4SAVE
C***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)
save
C***BEGIN PROLOGUE D1MACH
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 double precision machine dependent constants
C***DESCRIPTION
C
C This is the CMLIB version of D1MACH, the double precision machine
C constants subroutine originally developed for the PORT library.
C
C D1MACH can be used to obtain machine-dependent parameters
C for the local machine environment. It is a function
C subprogram with one (input) argument, and can be called
C as follows, for example
C
C D = D1MACH(I)
C
C where I=1,...,5. The (output) value of D above is
C determined by the (input) value of I. The results for
C various values of I are discussed below.
C
C Double-precision machine constants
C 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 A
C PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL
C SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188.
C***ROUTINES CALLED XERROR
C***END PROLOGUE D1MACH
C
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))
C
C MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T
C 3B SERIES AND MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T
C PC 7300), IN WHICH THE MOST SIGNIFICANT BYTE IS STORED FIRST.
C
C === MACHINE = IEEE.MOST-SIG-BYTE-FIRST
C === MACHINE = SUN
C === MACHINE = 68000
C === MACHINE = ATT.3B
C === MACHINE = ATT.7300
C 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 /
C
C MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES AND 8087-BASED
C MICROS, SUCH AS THE IBM PC AND AT&T 6300, IN WHICH THE LEAST
C SIGNIFICANT BYTE IS STORED FIRST.
C
C === MACHINE = IEEE.LEAST-SIG-BYTE-FIRST
C === MACHINE = 8087
C === MACHINE = IBM.PC
C === 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 /
C
C MACHINE CONSTANTS FOR AMDAHL MACHINES.
C
C === MACHINE = AMDAHL
C 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 /
C
C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.
C
C === MACHINE = BURROUGHS.1700
C 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 /
C
C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM.
C
C === MACHINE = BURROUGHS.5700
C 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 /
C
C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS.
C
C === MACHINE = BURROUGHS.6700
C === MACHINE = BURROUGHS.7700
C 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 /
C
C MACHINE CONSTANTS FOR THE CONVEX C-120 (NATIVE MODE)
C WITH OR WITHOUT -R8 OPTION
C
C === MACHINE = CONVEX.C1
C === MACHINE = CONVEX.C1.R8
C 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 /
C
C MACHINE CONSTANTS FOR THE CONVEX C-120 (IEEE MODE)
C WITH OR WITHOUT -R8 OPTION
C
C === MACHINE = CONVEX.C1.IEEE
C === MACHINE = CONVEX.C1.IEEE.R8
C 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 /
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 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" /
C
C MACHINE CONSTANTS FOR THE CDC 180 SERIES USING NOS/VE
C
C === MACHINE = CYBER.180.NOS/VE
C 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" /
C
C MACHINE CONSTANTS FOR THE CYBER 205
C
C === MACHINE = CYBER.205
C 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' /
C
C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES.
C
C === MACHINE = CDC.6000
C === MACHINE = CDC.7000
C 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 /
C
C MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3.
C
C === MACHINE = CRAY.46-BIT-INTEGER
C === MACHINE = CRAY.64-BIT-INTEGER
C DATA 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 /
C
C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200
C
C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING LINE -
C STATIC DMACH(5)
C
C === MACHINE = DATA_GENERAL.ECLIPSE.S/200
C 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/
C
C ELXSI 6400
C
C === MACHINE = ELSXI.6400
C 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 /
C
C MACHINE CONSTANTS FOR THE HARRIS 220
C MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7
C
C === MACHINE = HARRIS.220
C === MACHINE = HARRIS.SLASH6
C === MACHINE = HARRIS.SLASH7
C DATA SMALL(1),SMALL(2) / '20000000, '00000201 /
C DATA LARGE(1),LARGE(2) / '37777777, '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 /
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 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 /
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 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 /
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 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 /
C
C HP 9000
C
C D1MACH(1) = 2.8480954D-306
C D1MACH(2) = 1.40444776D+306
C D1MACH(3) = 2.22044605D-16
C D1MACH(4) = 4.44089210D-16
C D1MACH(5) = 3.01029996D-1
C
C === MACHINE = HP.9000
C 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 /
C
C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
C THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND
C THE INTERDATA 3230 AND INTERDATA 7/32.
C
C === MACHINE = IBM.360
C === MACHINE = IBM.370
C === MACHINE = XEROX.SIGMA.5
C === MACHINE = XEROX.SIGMA.7
C === MACHINE = XEROX.SIGMA.9
C === MACHINE = SEL.85
C === MACHINE = SEL.86
C === MACHINE = INTERDATA.3230
C === MACHINE = INTERDATA.7/32
C DATA 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 /
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 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' /
C
C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR).
C
C === MACHINE = PDP-10.KA
C 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 /
C
C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR).
C
C === MACHINE = PDP-10.KI
C 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 /
C
C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
C
C === MACHINE = PDP-11.32-BIT
C DATA SMALL(1),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 /
C
C 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 /
C
C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
C
C === MACHINE = PDP-11.16-BIT
C DATA SMALL(1),SMALL(2) / 128, 0 /
C DATA 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 /
C
C 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 /
C
C MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000
C
C === MACHINE = SEQUENT.BALANCE.8000
C 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 /
C
C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILER
C
C === MACHINE = UNIVAC.1100
C 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 /
C
C MACHINE CONSTANTS FOR VAX 11/780
C (EXPRESSED IN INTEGER AND HEXADECIMAL)
C *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***
C
C === MACHINE = VAX.11/780
C DATA SMALL(1), 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 /
C
C ***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 /
C
C MACHINE CONSTANTS FOR VAX 11/780 (G-FLOATING)
C (EXPRESSED IN INTEGER AND HEXADECIMAL)
C *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***
C
C 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 /
C
C ***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 /
C
C
C***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)
RETURN
C
END
SUBROUTINE DASSST(D, IV, P, STEP, STLSTG, V, X, X0)
save
C
C *** ASSESS CANDIDATE STEP (***SOL VERSION 2.3) ***
C
INTEGER P, IV(32)
DOUBLE PRECISION D(P), STEP(P), STLSTG(P), V(37), X(P), X0(P)
C
C *** PURPOSE ***
C
C THIS SUBROUTINE IS CALLED BY AN UNCONSTRAINED MINIMIZATION
C ROUTINE TO ASSESS THE NEXT CANDIDATE STEP. IT MAY RECOMMEND ONE
C OF SEVERAL COURSES OF ACTION, SUCH AS ACCEPTING THE STEP, RECOM-
C PUTING IT USING THE SAME OR A NEW QUADRATIC MODEL, OR HALTING DUE
C TO CONVERGENCE OR FALSE CONVERGENCE. SEE THE RETURN CODE LISTING
C BELOW.
C
C-------------------------- PARAMETER USAGE --------------------------
C
C IV (I/O) INTEGER PARAMETER AND SCRATCH VECTOR -- SEE DESCRIPTION
C BELOW OF IV VALUES REFERENCED.
C D (IN) SCALE VECTOR USED IN COMPUTING V(RELDX) -- SEE BELOW.
C P (IN) NUMBER OF PARAMETERS BEING OPTIMIZED.
C STEP (I/O) ON INPUT, STEP IS THE STEP TO BE ASSESSED. IT IS UN-
C CHANGED ON OUTPUT UNLESS A PREVIOUS STEP ACHIEVED A
C BETTER OBJECTIVE FUNCTION REDUCTION, IN WHICH CASE STLSTG
C WILL HAVE BEEN COPIED TO STEP.
C STLSTG (I/O) WHEN ASSESS RECOMMENDS RECOMPUTING STEP EVEN THOUGH THE
C CURRENT (OR A PREVIOUS) STEP YIELDS AN OBJECTIVE FUNC-
C TION DECREASE, IT SAVES IN STLSTG THE STEP THAT GAVE THE
C BEST FUNCTION REDUCTION SEEN SO FAR (IN THE CURRENT ITERA-
C TION). IF THE RECOMPUTED STEP YIELDS A LARGER FUNCTION
C VALUE, THEN STEP IS RESTORED FROM STLSTG AND
C X = X0 + STEP IS RECOMPUTED.
C V (I/O) REAL PARAMETER AND SCRATCH VECTOR -- SEE DESCRIPTION
C BELOW OF V VALUES REFERENCED.
C X (I/O) ON INPUT, X = X0 + STEP IS THE POINT AT WHICH THE OBJEC-
C TIVE FUNCTION HAS JUST BEEN EVALUATED. IF AN EARLIER
C STEP YIELDED A BIGGER FUNCTION DECREASE, THEN X IS
C RESTORED TO THE CORRESPONDING EARLIER VALUE. OTHERWISE,
C IF THE CURRENT STEP DOES NOT GIVE ANY FUNCTION DECREASE,
C THEN X IS RESTORED TO X0.
C X0 (IN) INITIAL OBJECTIVE FUNCTION PARAMETER VECTOR (AT THE
C START OF THE CURRENT ITERATION).
C
C *** IV VALUES REFERENCED ***
C
C IV(IRC) (I/O) ON INPUT FOR THE FIRST STEP TRIED IN A NEW ITERATION,
C IV(IRC) SHOULD BE SET TO 3 OR 4 (THE VALUE TO WHICH IT IS
C SET WHEN STEP IS DEFINITELY TO BE ACCEPTED). ON INPUT
C AFTER STEP HAS BEEN RECOMPUTED, IV(IRC) SHOULD BE
C UNCHANGED SINCE THE PREVIOUS RETURN OF ASSESS.
C ON OUTPUT, IV(IRC) IS A RETURN CODE HAVING ONE OF THE
C FOLLOWING VALUES...
C 1 = SWITCH MODELS OR TRY SMALLER STEP.
C 2 = SWITCH MODELS OR ACCEPT STEP.
C 3 = ACCEPT STEP AND DETERMINE V(RADFAC) BY GRADIENT
C TESTS.
C 4 = ACCEPT STEP, V(RADFAC) HAS BEEN DETERMINED.
C 5 = RECOMPUTE STEP (USING THE SAME MODEL).
C 6 = RECOMPUTE STEP WITH RADIUS = V(LMAXS) BUT DO NOT
C EVAULATE THE OBJECTIVE FUNCTION.
C 7 = X-CONVERGENCE (SEE V(XCTOL)).
C 8 = RELATIVE FUNCTION CONVERGENCE (SEE V(RFCTOL)).
C 9 = BOTH X- AND RELATIVE FUNCTION CONVERGENCE.
C 10 = ABSOLUTE FUNCTION CONVERGENCE (SEE V(AFCTOL)).
C 11 = SINGULAR CONVERGENCE (SEE V(LMAXS)).
C 12 = FALSE CONVERGENCE (SEE V(XFTOL)).
C 13 = IV(IRC) WAS OUT OF RANGE ON INPUT.
C RETURN CODE I HAS PRECDENCE OVER I+1 FOR I = 9, 10, 11.
C IV(MLSTGD) (I/O) SAVED VALUE OF IV(MODEL).
C IV(MODEL) (I/O) ON INPUT, IV(MODEL) SHOULD BE AN INTEGER IDENTIFYING
C THE CURRENT QUADRATIC MODEL OF THE OBJECTIVE FUNCTION.
C IF A PREVIOUS STEP YIELDED A BETTER FUNCTION REDUCTION,
C THEN IV(MODEL) WILL BE SET TO IV(MLSTGD) ON OUTPUT.
C IV(NFCALL) (IN) INVOCATION COUNT FOR THE OBJECTIVE FUNCTION.
C IV(NFGCAL) (I/O) VALUE OF IV(NFCALL) AT STEP THAT GAVE THE BIGGEST
C FUNCTION REDUCTION THIS ITERATION. IV(NFGCAL) REMAINS
C UNCHANGED UNTIL A FUNCTION REDUCTION IS OBTAINED.
C IV(RADINC) (I/O) THE NUMBER OF RADIUS INCREASES (OR MINUS THE NUMBER
C OF DECREASES) SO FAR THIS ITERATION.
C IV(RESTOR) (OUT) SET TO 0 UNLESS X AND V(F) HAVE BEEN RESTORED, IN
C WHICH CASE ASSESS SETS IV(RESTOR) = 1.
C IV(STAGE) (I/O) COUNT OF THE NUMBER OF MODELS TRIED SO FAR IN THE
C CURRENT ITERATION.
C IV(STGLIM) (IN) MAXIMUM NUMBER OF MODELS TO CONSIDER.
C IV(SWITCH) (OUT) SET TO 0 UNLESS A NEW MODEL IS BEING TRIED AND IT
C GIVES A SMALLER FUNCTION VALUE THAN THE PREVIOUS MODEL,
C IN WHICH CASE ASSESS SETS IV(SWITCH) = 1.
C IV(TOOBIG) (IN) IS NONZERO IF STEP WAS TOO BIG (E.G. IF IT CAUSED
C OVERFLOW).
C IV(XIRC) (I/O) VALUE THAT IV(IRC) WOULD HAVE IN THE ABSENCE OF
C CONVERGENCE, FALSE CONVERGENCE, AND OVERSIZED STEPS.
C
C *** V VALUES REFERENCED ***
C
C V(AFCTOL) (IN) ABSOLUTE FUNCTION CONVERGENCE TOLERANCE. IF THE
C ABSOLUTE VALUE OF THE CURRENT FUNCTION VALUE V(F) IS LESS
C THAN V(AFCTOL), THEN ASSESS RETURNS WITH IV(IRC) = 10.
C V(DECFAC) (IN) FACTOR BY WHICH TO DECREASE RADIUS WHEN IV(TOOBIG) IS
C NONZERO.
C V(DSTNRM) (IN) THE 2-NORM OF D*STEP.
C V(DSTSAV) (I/O) VALUE OF V(DSTNRM) ON SAVED STEP.
C V(DST0) (IN) THE 2-NORM OF D TIMES THE NEWTON STEP (WHEN DEFINED,
C I.E., FOR V(NREDUC) .GE. 0).
C V(F) (I/O) ON BOTH INPUT AND OUTPUT, V(F) IS THE OBJECTIVE FUNC-
C TION VALUE AT X. IF X IS RESTORED TO A PREVIOUS VALUE,
C THEN V(F) IS RESTORED TO THE CORRESPONDING VALUE.
C V(FDIF) (OUT) THE FUNCTION REDUCTION V(F0) - V(F) (FOR THE OUTPUT
C VALUE OF V(F) IF AN EARLIER STEP GAVE A BIGGER FUNCTION
C DECREASE, AND FOR THE INPUT VALUE OF V(F) OTHERWISE).
C V(FLSTGD) (I/O) SAVED VALUE OF V(F).
C V(F0) (IN) OBJECTIVE FUNCTION VALUE AT START OF ITERATION.
C V(GTSLST) (I/O) VALUE OF V(GTSTEP) ON SAVED STEP.
C V(GTSTEP) (IN) INNER PRODUCT BETWEEN STEP AND GRADIENT.
C V(INCFAC) (IN) MINIMUM FACTOR BY WHICH TO INCREASE RADIUS.
C V(LMAXS) (IN) MAXIMUM REASONABLE STEP SIZE (AND INITIAL STEP BOUND).
C IF THE ACTUAL FUNCTION DECREASE IS NO MORE THAN TWICE
C WHAT WAS PREDICTED, IF A RETURN WITH IV(IRC) = 7, 8, 9,
C OR 10 DOES NOT OCCUR, IF V(DSTNRM) .GT. V(LMAXS), AND IF
C V(PREDUC) .LE. V(SCTOL) * ABS(V(F0)), THEN ASSESS RE-
C TURNS WITH IV(IRC) = 11. IF SO DOING APPEARS WORTHWHILE,
C THEN ASSESS REPEATS THIS TEST WITH V(PREDUC) COMPUTED FOR
C A STEP OF LENGTH V(LMAXS) (BY A RETURN WITH IV(IRC) = 6).
C V(NREDUC) (I/O) FUNCTION REDUCTION PREDICTED BY QUADRATIC MODEL FOR
C NEWTON STEP. IF ASSESS IS CALLED WITH IV(IRC) = 6, I.E.,
C IF V(PREDUC) HAS BEEN COMPUTED WITH RADIUS = V(LMAXS) FOR
C USE IN THE SINGULAR CONVERVENCE TEST, THEN V(NREDUC) IS
C SET TO -V(PREDUC) BEFORE THE LATTER IS RESTORED.
C V(PLSTGD) (I/O) VALUE OF V(PREDUC) ON SAVED STEP.
C V(PREDUC) (I/O) FUNCTION REDUCTION PREDICTED BY QUADRATIC MODEL FOR
C CURRENT STEP.
C V(RADFAC) (OUT) FACTOR TO BE USED IN DETERMINING THE NEW RADIUS,
C WHICH SHOULD BE V(RADFAC)*DST, WHERE DST IS EITHER THE
C OUTPUT VALUE OF V(DSTNRM) OR THE 2-NORM OF
C DIAG(NEWD)*STEP FOR THE OUTPUT VALUE OF STEP AND THE
C UPDATED VERSION, NEWD, OF THE SCALE VECTOR D. FOR
C IV(IRC) = 3, V(RADFAC) = 1.0 IS RETURNED.
C V(RDFCMN) (IN) MINIMUM VALUE FOR V(RADFAC) IN TERMS OF THE INPUT
C VALUE OF V(DSTNRM) -- SUGGESTED VALUE = 0.1.
C V(RDFCMX) (IN) MAXIMUM VALUE FOR V(RADFAC) -- SUGGESTED VALUE = 4.0.
C V(RELDX) (OUT) SCALED RELATIVE CHANGE IN X CAUSED BY STEP, COMPUTED
C BY FUNCTION DRELST AS
C MAX (D(I)*ABS(X(I)-X0(I)), 1 .LE. I .LE. P) /
C MAX (D(I)*(ABS(X(I))+ABS(X0(I))), 1 .LE. I .LE. P).
C IF AN ACCEPTABLE STEP IS RETURNED, THEN V(RELDX) IS COM-
C PUTED USING THE OUTPUT (POSSIBLY RESTORED) VALUES OF X
C AND STEP. OTHERWISE IT IS COMPUTED USING THE INPUT
C VALUES.
C V(RFCTOL) (IN) RELATIVE FUNCTION CONVERGENCE TOLERANCE. IF THE
C ACTUAL FUNCTION REDUCTION IS AT MOST TWICE WHAT WAS PRE-
C DICTED AND V(NREDUC) .LE. V(RFCTOL)*ABS(V(F0)), THEN
C ASSESS RETURNS WITH IV(IRC) = 8 OR 9.
C V(STPPAR) (IN) MARQUARDT PARAMETER -- 0 MEANS FULL NEWTON STEP.
C V(TUNER1) (IN) TUNING CONSTANT USED TO DECIDE IF THE FUNCTION
C REDUCTION WAS MUCH LESS THAN EXPECTED. SUGGESTED
C VALUE = 0.1.
C V(TUNER2) (IN) TUNING CONSTANT USED TO DECIDE IF THE FUNCTION
C REDUCTION WAS LARGE ENOUGH TO ACCEPT STEP. SUGGESTED
C VALUE = 10**-4.
C V(TUNER3) (IN) TUNING CONSTANT USED TO DECIDE IF THE RADIUS
C SHOULD BE INCREASED. SUGGESTED VALUE = 0.75.
C V(XCTOL) (IN) X-CONVERGENCE CRITERION. IF STEP IS A NEWTON STEP
C (V(STPPAR) = 0) HAVING V(RELDX) .LE. V(XCTOL) AND GIVING
C AT MOST TWICE THE PREDICTED FUNCTION DECREASE, THEN
C ASSESS RETURNS IV(IRC) = 7 OR 9.
C V(XFTOL) (IN) FALSE CONVERGENCE TOLERANCE. IF STEP GAVE NO OR ONLY
C A SMALL FUNCTION DECREASE AND V(RELDX) .LE. V(XFTOL),
C THEN ASSESS RETURNS WITH IV(IRC) = 12.
C
C------------------------------- NOTES -------------------------------
C
C *** APPLICATION AND USAGE RESTRICTIONS ***
C
C THIS ROUTINE IS CALLED AS PART OF THE NL2SOL (NONLINEAR
C LEAST-SQUARES) PACKAGE. IT MAY BE USED IN ANY UNCONSTRAINED
C MINIMIZATION SOLVER THAT USES DOGLEG, GOLDFELD-QUANDT-TROTTER,
C OR LEVENBERG-MARQUARDT STEPS.
C
C *** ALGORITHM NOTES ***
C
C SEE (1) FOR FURTHER DISCUSSION OF THE ASSESSING AND MODEL
C SWITCHING STRATEGIES. WHILE NL2SOL CONSIDERS ONLY TWO MODELS,
C ASSESS IS DESIGNED TO HANDLE ANY NUMBER OF MODELS.
C
C *** USAGE NOTES ***
C
C ON THE FIRST CALL OF AN ITERATION, ONLY THE I/O VARIABLES
C STEP, X, IV(IRC), IV(MODEL), V(F), V(DSTNRM), V(GTSTEP), AND
C V(PREDUC) NEED HAVE BEEN INITIALIZED. BETWEEN CALLS, NO I/O
C VALUES EXECPT STEP, X, IV(MODEL), V(F) AND THE STOPPING TOLER-
C ANCES SHOULD BE CHANGED.
C AFTER A RETURN FOR CONVERGENCE OR FALSE CONVERGENCE, ONE CAN
C CHANGE THE STOPPING TOLERANCES AND CALL ASSESS AGAIN, IN WHICH
C CASE THE STOPPING TESTS WILL BE REPEATED.
C
C *** REFERENCES ***
C
C (1) DENNIS, J.E., JR., GAY, D.M., AND WELSCH, R.E. (1981),
C AN ADAPTIVE NONLINEAR LEAST-SQUARES ALGORITHM,
C ACM TRANS. MATH. SOFTWARE, VOL. 7, NO. 3.
C
C (2) POWELL, M.J.D. (1970) A FORTRAN SUBROUTINE FOR SOLVING
C SYSTEMS OF NONLINEAR ALGEBRAIC EQUATIONS, IN NUMERICAL
C METHODS FOR NONLINEAR ALGEBRAIC EQUATIONS, EDITED BY
C P. RABINOWITZ, GORDON AND BREACH, LONDON.
C
C *** HISTORY ***
C
C JOHN DENNIS DESIGNED MUCH OF THIS ROUTINE, STARTING WITH
C IDEAS IN (2). ROY WELSCH SUGGESTED THE MODEL SWITCHING STRATEGY.
C DAVID GAY AND STEPHEN PETERS CAST THIS SUBROUTINE INTO A MORE
C PORTABLE FORM (WINTER 1977), AND DAVID GAY CAST IT INTO ITS
C PRESENT FORM (FALL 1978).
C
C *** GENERAL ***
C
C THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH
C SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS
C MCS-7600324, DCR75-10143, 76-14311DSS, MCS76-11989, AND
C MCS-7906671.
C
C------------------------ EXTERNAL QUANTITIES ------------------------
C
C *** EXTERNAL FUNCTIONS AND SUBROUTINES ***
C
EXTERNAL DRELST
DOUBLE PRECISION DRELST
C
C DRELST... COMPUTES V(RELDX) = RELATIVE STEP SIZE.
C
C *** INTRINSIC FUNCTIONS ***
C/+
DOUBLE PRECISION DABS, DMAX1
C/
C *** NO COMMON BLOCKS ***
C
C-------------------------- LOCAL VARIABLES --------------------------
C
LOGICAL GOODX
INTEGER I, NFC
DOUBLE PRECISION EMAX, EMAXS, GTS, HALF, ONE, RELDX1, RFAC1, TWO,
1 XMAX, ZERO
C
C *** SUBSCRIPTS FOR IV AND V ***
C
INTEGER AFCTOL, DECFAC, DSTNRM, DSTSAV, DST0, F, FDIF, FLSTGD, F0,
1 GTSLST, GTSTEP, INCFAC, IRC, LMAXS, MLSTGD, MODEL, NFCALL,
2 NFGCAL, NREDUC, PLSTGD, PREDUC, RADFAC, RADINC, RDFCMN,
3 RDFCMX, RELDX, RESTOR, RFCTOL, SCTOL, STAGE, STGLIM,
4 STPPAR, SWITCH, TOOBIG, TUNER1, TUNER2, TUNER3, XCTOL,
5 XFTOL, XIRC
C
C *** DATA INITIALIZATIONS ***
C
C/6
DATA HALF/0.5D+0/, ONE/1.D+0/, TWO/2.D+0/, ZERO/0.D+0/
C/7
C PARAMETER (HALF=0.5D+0, ONE=1.D+0, TWO=2.D+0, ZERO=0.D+0)
C/
C
C/6
DATA IRC/29/, MLSTGD/32/, MODEL/5/, NFCALL/6/, NFGCAL/7/,
1 RADINC/8/, RESTOR/9/, STAGE/10/, STGLIM/11/, SWITCH/12/,
2 TOOBIG/2/, XIRC/13/
C/7
C PARAMETER (IRC=29, MLSTGD=32, MODEL=5, NFCALL=6, NFGCAL=7,
C 1 RADINC=8, RESTOR=9, STAGE=10, STGLIM=11, SWITCH=12,
C 2 TOOBIG=2, XIRC=13)
C/
C/6
DATA AFCTOL/31/, DECFAC/22/, DSTNRM/2/, DST0/3/, DSTSAV/18/,
1 F/10/, FDIF/11/, FLSTGD/12/, F0/13/, GTSLST/14/, GTSTEP/4/,
2 INCFAC/23/, LMAXS/36/, NREDUC/6/, PLSTGD/15/, PREDUC/7/,
3 RADFAC/16/, RDFCMN/24/, RDFCMX/25/, RELDX/17/, RFCTOL/32/,
4 SCTOL/37/, STPPAR/5/, TUNER1/26/, TUNER2/27/, TUNER3/28/,
5 XCTOL/33/, XFTOL/34/
C/7
C PARAMETER (AFCTOL=31, DECFAC=22, DSTNRM=2, DST0=3, DSTSAV=18,
C 1 F=10, FDIF=11, FLSTGD=12, F0=13, GTSLST=14, GTSTEP=4,
C 2 INCFAC=23, LMAXS=36, NREDUC=6, PLSTGD=15, PREDUC=7,
C 3 RADFAC=16, RDFCMN=24, RDFCMX=25, RELDX=17, RFCTOL=32,
C 4 SCTOL=37, STPPAR=5, TUNER1=26, TUNER2=27, TUNER3=28,
C 5 XCTOL=33, XFTOL=34)
C/
C
C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++
C
NFC = IV(NFCALL)
IV(SWITCH) = 0
IV(RESTOR) = 0
RFAC1 = ONE
GOODX = .TRUE.
I = IV(IRC)
IF (I .GE. 1 .AND. I .LE. 12)
1 GO TO (20,30,10,10,40,280,220,220,220,220,220,170), I
IV(IRC) = 13
GO TO 999
C
C *** INITIALIZE FOR NEW ITERATION ***
C
10 IV(STAGE) = 1
IV(RADINC) = 0
V(FLSTGD) = V(F0)
IF (IV(TOOBIG) .EQ. 0) GO TO 90
IV(STAGE) = -1
IV(XIRC) = I
GO TO 60
C
C *** STEP WAS RECOMPUTED WITH NEW MODEL OR SMALLER RADIUS ***
C *** FIRST DECIDE WHICH ***
C
20 IF (IV(MODEL) .NE. IV(MLSTGD)) GO TO 30
C *** OLD MODEL RETAINED, SMALLER RADIUS TRIED ***
C *** DO NOT CONSIDER ANY MORE NEW MODELS THIS ITERATION ***
IV(STAGE) = IV(STGLIM)
IV(RADINC) = -1
GO TO 90
C
C *** A NEW MODEL IS BEING TRIED. DECIDE WHETHER TO KEEP IT. ***
C
30 IV(STAGE) = IV(STAGE) + 1
C
C *** NOW WE ADD THE POSSIBILTIY THAT STEP WAS RECOMPUTED WITH ***
C *** THE SAME MODEL, PERHAPS BECAUSE OF AN OVERSIZED STEP. ***
C
40 IF (IV(STAGE) .GT. 0) GO TO 50
C
C *** STEP WAS RECOMPUTED BECAUSE IT WAS TOO BIG. ***
C
IF (IV(TOOBIG) .NE. 0) GO TO 60
C
C *** RESTORE IV(STAGE) AND PICK UP WHERE WE LEFT OFF. ***
C
IV(STAGE) = -IV(STAGE)
I = IV(XIRC)
GO TO (20, 30, 90, 90, 70), I
C
50 IF (IV(TOOBIG) .EQ. 0) GO TO 70
C
C *** HANDLE OVERSIZE STEP ***
C
IF (IV(RADINC) .GT. 0) GO TO 80
IV(STAGE) = -IV(STAGE)
IV(XIRC) = IV(IRC)
C
60 V(RADFAC) = V(DECFAC)
IV(RADINC) = IV(RADINC) - 1
IV(IRC) = 5
GO TO 999
C
70 IF (V(F) .LT. V(FLSTGD)) GO TO 90
C
C *** THE NEW STEP IS A LOSER. RESTORE OLD MODEL. ***
C
IF (IV(MODEL) .EQ. IV(MLSTGD)) GO TO 80
IV(MODEL) = IV(MLSTGD)
IV(SWITCH) = 1
C
C *** RESTORE STEP, ETC. ONLY IF A PREVIOUS STEP DECREASED V(F).
C
80 IF (V(FLSTGD) .GE. V(F0)) GO TO 90
IV(RESTOR) = 1
V(F) = V(FLSTGD)
V(PREDUC) = V(PLSTGD)
V(GTSTEP) = V(GTSLST)
IF (IV(SWITCH) .EQ. 0) RFAC1 = V(DSTNRM) / V(DSTSAV)
V(DSTNRM) = V(DSTSAV)
NFC = IV(NFGCAL)
GOODX = .FALSE.
C
C
C *** COMPUTE RELATIVE CHANGE IN X BY CURRENT STEP ***
C
90 RELDX1 = DRELST(P, D, X, X0)
C
C *** RESTORE X AND STEP IF NECESSARY ***
C
IF (GOODX) GO TO 110
DO 100 I = 1, P
STEP(I) = STLSTG(I)
X(I) = X0(I) + STLSTG(I)
100 CONTINUE
C
110 V(FDIF) = V(F0) - V(F)
IF (V(FDIF) .GT. V(TUNER2) * V(PREDUC)) GO TO 140
C
C *** NO (OR ONLY A TRIVIAL) FUNCTION DECREASE
C *** -- SO TRY NEW MODEL OR SMALLER RADIUS
C
V(RELDX) = RELDX1
IF (V(F) .LT. V(F0)) GO TO 120
IV(MLSTGD) = IV(MODEL)
V(FLSTGD) = V(F)
V(F) = V(F0)
CALL 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 160
C
C *** NONTRIVIAL FUNCTION DECREASE ACHIEVED ***
C
140 IV(NFGCAL) = NFC
RFAC1 = ONE
IF (GOODX) V(RELDX) = RELDX1
V(DSTSAV) = V(DSTNRM)
IF (V(FDIF) .GT. V(PREDUC)*V(TUNER1)) GO TO 190
C
C *** DECREASE WAS MUCH LESS THAN PREDICTED -- EITHER CHANGE MODELS
C *** OR ACCEPT STEP WITH DECREASED RADIUS.
C
IF (IV(STAGE) .GE. IV(STGLIM)) GO TO 150
C *** CONSIDER SWITCHING MODELS ***
IV(IRC) = 2
GO TO 160
C
C *** ACCEPT STEP WITH DECREASED RADIUS ***
C
150 IV(IRC) = 4
C
C *** SET V(RADFAC) TO FLETCHER*S DECREASE FACTOR ***
C
160 IV(XIRC) = IV(IRC)
EMAX = V(GTSTEP) + V(FDIF)
V(RADFAC) = HALF * RFAC1
IF (EMAX .LT. V(GTSTEP)) V(RADFAC) = RFAC1 * DMAX1(V(RDFCMN),
1 HALF * V(GTSTEP)/EMAX)
C
C *** DO FALSE CONVERGENCE TEST ***
C
170 IF (V(RELDX) .LE. V(XFTOL)) GO TO 180
IV(IRC) = IV(XIRC)
IF (V(F) .LT. V(F0)) GO TO 200
GO TO 230
C
180 IV(IRC) = 12
GO TO 240
C
C *** HANDLE GOOD FUNCTION DECREASE ***
C
190 IF (V(FDIF) .LT. (-V(TUNER3) * V(GTSTEP))) GO TO 210
C
C *** INCREASING RADIUS LOOKS WORTHWHILE. SEE IF WE JUST
C *** RECOMPUTED STEP WITH A DECREASED RADIUS OR RESTORED STEP
C *** AFTER RECOMPUTING IT WITH A LARGER RADIUS.
C
IF (IV(RADINC) .LT. 0) GO TO 210
IF (IV(RESTOR) .EQ. 1) GO TO 210
C
C *** WE DID NOT. TRY A LONGER STEP UNLESS THIS WAS A NEWTON
C *** STEP.
C
V(RADFAC) = V(RDFCMX)
GTS = V(GTSTEP)
IF (V(FDIF) .LT. (HALF/V(RADFAC) - ONE) * GTS)
1 V(RADFAC) = DMAX1(V(INCFAC), HALF*GTS/(GTS + V(FDIF)))
IV(IRC) = 4
IF (V(STPPAR) .EQ. ZERO) GO TO 230
C *** STEP WAS NOT A NEWTON STEP. RECOMPUTE IT WITH
C *** A LARGER RADIUS.
IV(IRC) = 5
IV(RADINC) = IV(RADINC) + 1
C
C *** SAVE VALUES CORRESPONDING TO GOOD STEP ***
C
200 V(FLSTGD) = V(F)
IV(MLSTGD) = IV(MODEL)
CALL DCOPY(P, STEP,1,STLSTG,1)
V(DSTSAV) = V(DSTNRM)
IV(NFGCAL) = NFC
V(PLSTGD) = V(PREDUC)
V(GTSLST) = V(GTSTEP)
GO TO 230
C
C *** ACCEPT STEP WITH RADIUS UNCHANGED ***
C
210 V(RADFAC) = ONE
IV(IRC) = 3
GO TO 230
C
C *** COME HERE FOR A RESTART AFTER CONVERGENCE ***
C
220 IV(IRC) = IV(XIRC)
IF (V(DSTSAV) .GE. ZERO) GO TO 240
IV(IRC) = 12
GO TO 240
C
C *** PERFORM CONVERGENCE TESTS ***
C
230 IV(XIRC) = IV(IRC)
240 IF (DABS(V(F)) .LT. V(AFCTOL)) IV(IRC) = 10
IF (HALF * V(FDIF) .GT. V(PREDUC)) GO TO 999
EMAX = V(RFCTOL) * DABS(V(F0))
EMAXS = V(SCTOL) * DABS(V(F0))
IF (V(DSTNRM) .GT. V(LMAXS) .AND. V(PREDUC) .LE. EMAXS)
1 IV(IRC) = 11
IF (V(DST0) .LT. ZERO) GO TO 250
I = 0
IF ((V(NREDUC) .GT. ZERO .AND. V(NREDUC) .LE. EMAX) .OR.
1 (V(NREDUC) .EQ. ZERO. AND. V(PREDUC) .EQ. ZERO)) I = 2
IF (V(STPPAR) .EQ. ZERO .AND. V(RELDX) .LE. V(XCTOL)
1 .AND. GOODX) I = I + 1
IF (I .GT. 0) IV(IRC) = I + 6
C
C *** CONSIDER RECOMPUTING STEP OF LENGTH V(LMAXS) FOR SINGULAR
C *** CONVERGENCE TEST.
C
250 IF (IV(IRC) .GT. 5 .AND. IV(IRC) .NE. 12) GO TO 999
IF (V(DSTNRM) .GT. V(LMAXS)) GO TO 260
IF (V(PREDUC) .GE. EMAXS) GO TO 999
IF (V(DST0) .LE. ZERO) GO TO 270
IF (HALF * V(DST0) .LE. V(LMAXS)) GO TO 999
GO TO 270
260 IF (HALF * V(DSTNRM) .LE. V(LMAXS)) GO TO 999
XMAX = V(LMAXS) / V(DSTNRM)
IF (XMAX * (TWO - XMAX) * V(PREDUC) .GE. EMAXS) GO TO 999
270 IF (V(NREDUC) .LT. ZERO) GO TO 290
C
C *** RECOMPUTE V(PREDUC) FOR USE IN SINGULAR CONVERGENCE TEST ***
C
V(GTSLST) = V(GTSTEP)
V(DSTSAV) = V(DSTNRM)
IF (IV(IRC) .EQ. 12) V(DSTSAV) = -V(DSTSAV)
V(PLSTGD) = V(PREDUC)
IV(IRC) = 6
CALL DCOPY(P, STEP,1,STLSTG,1)
GO TO 999
C
C *** PERFORM SINGULAR CONVERGENCE TEST WITH RECOMPUTED V(PREDUC) ***
C
280 V(GTSTEP) = V(GTSLST)
V(DSTNRM) = DABS(V(DSTSAV))
CALL 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) = 11
C
999 RETURN
C
C *** LAST CARD OF ASSESS FOLLOWS ***
END
SUBROUTINE DCOPY(N,DX,INCX,DY,INCY)
save
C***BEGIN PROLOGUE DCOPY
C***DATE WRITTEN 791001 (YYMMDD)
C***REVISION DATE 820801 (YYMMDD)
C***CATEGORY NO. D1A5
C***KEYWORDS BLAS,COPY,DOUBLE PRECISION,LINEAR ALGEBRA,VECTOR
C***AUTHOR LAWSON, C. L., (JPL)
C HANSON, R. J., (SNLA)
C KINCAID, D. R., (U. OF TEXAS)
C KROGH, F. T., (JPL)
C***PURPOSE D.P. vector copy y = x
C***DESCRIPTION
C
C B L A S Subprogram
C Description of Parameters
C
C --Input--
C N number of elements in input vector(s)
C DX double precision vector with N elements
C INCX storage spacing between elements of DX
C DY double precision vector with N elements
C INCY storage spacing between elements of DY
C
C --Output--
C DY copy of vector DX (unchanged if N .LE. 0)
C
C Copy double precision DX to double precision DY.
C For I = 0 to N-1, copy DX(LX+I*INCX) to DY(LY+I*INCY),
C where LX = 1 if INCX .GE. 0, else LX = (-INCX)*N, and LY is
C defined in a similar way using INCY.
C***REFERENCES LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T.,
C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*,
C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL
C SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323
C***ROUTINES CALLED (NONE)
C***END PROLOGUE DCOPY
C
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 CONTINUE
C
C CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS.
C
IX = 1
IY = 1
IF(INCX.LT.0)IX = (-N+1)*INCX + 1
IF(INCY.LT.0)IY = (-N+1)*INCY + 1
DO 10 I = 1,N
DY(IY) = DX(IX)
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
RETURN
C
C CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 7.
C
20 M = MOD(N,7)
IF( M .EQ. 0 ) GO TO 40
DO 30 I = 1,M
DY(I) = DX(I)
30 CONTINUE
IF( N .LT. 7 ) RETURN
40 MP1 = M + 1
DO 50 I = MP1,N,7
DY(I) = DX(I)
DY(I + 1) = DX(I + 1)
DY(I + 2) = DX(I + 2)
DY(I + 3) = DX(I + 3)
DY(I + 4) = DX(I + 4)
DY(I + 5) = DX(I + 5)
DY(I + 6) = DX(I + 6)
50 CONTINUE
RETURN
C
C CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS.
C
60 CONTINUE
NS=N*INCX
DO 70 I=1,NS,INCX
DY(I) = DX(I)
70 CONTINUE
RETURN
END
SUBROUTINE DDBDOG(DIG, G, LV, N, NWTSTP, STEP, V)
save
C
C *** COMPUTE DOUBLE DOGLEG STEP ***
C
C *** PARAMETER DECLARATIONS ***
C
INTEGER LV, N
DOUBLE PRECISION DIG(N), G(N), NWTSTP(N), STEP(N), V(LV)
C
C *** PURPOSE ***
C
C THIS SUBROUTINE COMPUTES A CANDIDATE STEP (FOR USE IN AN UNCON-
C STRAINED MINIMIZATION CODE) BY THE DOUBLE DOGLEG ALGORITHM OF
C DENNIS AND MEI (REF. 1), WHICH IS A VARIATION ON POWELL*S DOGLEG
C SCHEME (REF. 2, P. 95).
C
C-------------------------- PARAMETER USAGE --------------------------
C
C 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 ARE
C USED HERE...
C V(BIAS) (INPUT) BIAS FOR RELAXED NEWTON STEP, WHICH IS V(BIAS) OF
C THE WAY FROM THE FULL NEWTON TO THE FULLY RELAXED NEWTON
C 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 -- SEE
C ALGORITHM NOTES.
C V(GTSTEP) (OUTPUT) INNER PRODUCT BETWEEN G AND STEP.
C V(NREDUC) (OUTPUT) FUNCTION REDUCTION PREDICTED FOR THE FULL NEWTON
C 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 RETURNED
C HAS 2-NORM V(RADIUS) UNLESS V(STPPAR) = 0.
C V(STPPAR) (OUTPUT) CODE TELLING HOW STEP WAS COMPUTED... 0 MEANS A
C FULL NEWTON STEP. BETWEEN 0 AND 1 MEANS V(STPPAR) OF THE
C WAY FROM THE NEWTON TO THE RELAXED NEWTON STEP. BETWEEN
C 1 AND 2 MEANS A TRUE DOUBLE DOGLEG STEP, V(STPPAR) - 1 OF
C THE WAY FROM THE RELAXED NEWTON TO THE CAUCHY STEP.
C GREATER THAN 2 MEANS 1 / (V(STPPAR) - 1) TIMES THE CAUCHY
C STEP.
C
C------------------------------- NOTES -------------------------------
C
C *** ALGORITHM NOTES ***
C
C LET G AND H BE THE CURRENT GRADIENT AND HESSIAN APPROXIMA-
C TION RESPECTIVELY AND LET D BE THE CURRENT SCALE VECTOR. THIS
C 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 H
C BY DIAG(D)**-1 * G AND DIAG(D)**-1 * H * DIAG(D)**-1,
C COMPUTING STEP, AND TRANSLATING STEP BACK TO THE ORIGINAL
C VARIABLES, I.E., PREMULTIPLYING IT BY DIAG(D)**-1.
C
C *** REFERENCES ***
C
C 1. 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 2. POWELL, M.J.D. (1970), A HYBRID METHOD FOR NON-LINEAR EQUATIONS,
C IN NUMERICAL METHODS FOR NON-LINEAR EQUATIONS, EDITED BY
C P. RABINOWITZ, GORDON AND BREACH, LONDON.
C
C *** GENERAL ***
C
C CODED BY DAVID M. GAY.
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
DOUBLE PRECISION DDOT
C
C
C *** INTRINSIC FUNCTIONS ***
C/+
DOUBLE PRECISION DSQRT
C/
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, ZERO
C
C *** V SUBSCRIPTS ***
C
INTEGER BIAS, DGNORM, DSTNRM, DST0, GRDFAC, GTHG, GTSTEP,
1 NREDUC, NWTFAC, PREDUC, RADIUS, STPPAR
C
C *** DATA INITIALIZATIONS ***
C
C/6
DATA HALF/0.5D+0/, ONE/1.D+0/, TWO/2.D+0/, ZERO/0.D+0/
C/7
C PARAMETER (HALF=0.5D+0, ONE=1.D+0, TWO=2.D+0, ZERO=0.D+0)
C/
C
C/6
DATA 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/7
C 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/
C
C+++++++++++++++++++++++++++++++ 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 30
C
C *** 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 999
C
30 V(DSTNRM) = V(RADIUS)
CFACT = (GNORM / V(GTHG))**2
C *** CAUCHY STEP = -CFACT * G.
CNORM = GNORM * CFACT
RELAX = ONE - V(BIAS) * (ONE - CNORM/GHINVG)
IF (RLAMBD .LT. RELAX) GO TO 50
C
C *** 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 999
C
50 IF (CNORM .LT. V(RADIUS)) GO TO 70
C
C *** 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 999
C
C *** COMPUTE DOGLEG STEP BETWEEN CAUCHY AND RELAXED NEWTON ***
C *** FEMUR = RELAXED NEWTON STEP MINUS CAUCHY STEP ***
C
70 CTRNWT = CFACT * RELAX * GHINVG / GNORM
C *** CTRNWT = INNER PROD. OF CAUCHY AND RELAXED NEWTON STEPS,
C *** SCALED BY GNORM**-2.
T1 = CTRNWT - CFACT**2
C *** T1 = INNER PROD. OF FEMUR AND CAUCHY STEP, SCALED BY
C *** GNORM**-2.
T2 = (V(RADIUS)/GNORM)**2 - CFACT**2
FEMNSQ = (RELAX*NWTNRM/GNORM)**2 - CTRNWT - T1
C *** 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 RETURN
C *** LAST CARD OF DDBDOG FOLLOWS ***
END
DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)
save
C***BEGIN PROLOGUE DDOT
C***DATE WRITTEN 791001 (YYMMDD)
C***REVISION DATE 820801 (YYMMDD)
C***CATEGORY NO. D1A4
C***KEYWORDS BLAS,DOUBLE PRECISION,INNER PRODUCT,LINEAR ALGEBRA,VECTOR
C***AUTHOR LAWSON, C. L., (JPL)
C HANSON, R. J., (SNLA)
C KINCAID, D. R., (U. OF TEXAS)
C KROGH, F. T., (JPL)
C***PURPOSE D.P. inner product of d.p. vectors
C***DESCRIPTION
C
C B L A S Subprogram
C Description of Parameters
C
C --Input--
C N number of elements in input vector(s)
C DX double precision vector with N elements
C INCX storage spacing between elements of DX
C DY double precision vector with N elements
C INCY storage spacing between elements of DY
C
C --Output--
C DDOT double precision dot product (zero if N .LE. 0)
C
C Returns the dot product of double precision DX and DY.
C DDOT = sum for I = 0 to N-1 of DX(LX+I*INCX) * DY(LY+I*INCY)
C where LX = 1 if INCX .GE. 0, else LX = (-INCX)*N, and LY is
C defined in a similar way using INCY.
C***REFERENCES LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T.,
C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*,
C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL
C SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323
C***ROUTINES CALLED (NONE)
C***END PROLOGUE DDOT
C
DOUBLE PRECISION DX(1),DY(1)
C***FIRST EXECUTABLE STATEMENT DDOT
DDOT = 0.D0
IF(N.LE.0)RETURN
IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60
5 CONTINUE
C
C CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS.
C
IX = 1
IY = 1
IF(INCX.LT.0)IX = (-N+1)*INCX + 1
IF(INCY.LT.0)IY = (-N+1)*INCY + 1
DO 10 I = 1,N
DDOT = DDOT + DX(IX)*DY(IY)
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
RETURN
C
C CODE FOR BOTH INCREMENTS EQUAL TO 1.
C
C
C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 5.
C
20 M = MOD(N,5)
IF( M .EQ. 0 ) GO TO 40
DO 30 I = 1,M
DDOT = DDOT + DX(I)*DY(I)
30 CONTINUE
IF( N .LT. 5 ) RETURN
40 MP1 = M + 1
DO 50 I = MP1,N,5
DDOT = DDOT + DX(I)*DY(I) + DX(I+1)*DY(I+1) +
1 DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4)
50 CONTINUE
RETURN
C
C CODE FOR POSITIVE EQUAL INCREMENTS .NE.1.
C
60 CONTINUE
NS = N*INCX
DO 70 I=1,NS,INCX
DDOT = DDOT + DX(I)*DY(I)
70 CONTINUE
RETURN
END
SUBROUTINE DITSUM(D, G, IV, LIV, LV, P, V, X)
save
C
C *** PRINT ITERATION SUMMARY FOR ***SOL (VERSION 2.3) ***
C
C *** PARAMETER DECLARATIONS ***
C
INTEGER LIV, LV, P
INTEGER IV(LIV)
DOUBLE PRECISION D(P), G(P), V(LV), X(P)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C *** LOCAL VARIABLES ***
C
INTEGER ALG, I, IV1, M, NF, NG, OL, PU
C/6
REAL MODEL1(6), MODEL2(6)
C/7
C CHARACTER*4 MODEL1(6), MODEL2(6)
C/
DOUBLE PRECISION NRELDF, OLDF, PRELDF, RELDF, ZERO
C
C *** INTRINSIC FUNCTIONS ***
C/+
INTEGER IABS
DOUBLE PRECISION DABS, DMAX1
C/
C *** NO EXTERNAL FUNCTIONS OR SUBROUTINES ***
C
C *** SUBSCRIPTS FOR IV AND V ***
C
INTEGER ALGSAV, DSTNRM, F, FDIF, F0, NEEDHD, NFCALL, NFCOV, NGCOV,
1 NGCALL, NITER, NREDUC, OUTLEV, PREDUC, PRNTIT, PRUNIT,
2 RELDX, SOLPRT, STATPR, STPPAR, SUSED, X0PRT
C
C *** IV SUBSCRIPT VALUES ***
C
C/6
DATA ALGSAV/51/, NEEDHD/36/, NFCALL/6/, NFCOV/52/, NGCALL/30/,
1 NGCOV/53/, NITER/31/, OUTLEV/19/, PRNTIT/39/, PRUNIT/21/,
2 SOLPRT/22/, STATPR/23/, SUSED/64/, X0PRT/24/
C/7
C PARAMETER (ALGSAV=51, NEEDHD=36, NFCALL=6, NFCOV=52, NGCALL=30,
C 1 NGCOV=53, NITER=31, OUTLEV=19, PRNTIT=39, PRUNIT=21,
C 2 SOLPRT=22, STATPR=23, SUSED=64, X0PRT=24)
C/
C
C *** V SUBSCRIPT VALUES ***
C
C/6
DATA DSTNRM/2/, F/10/, F0/13/, FDIF/11/, NREDUC/6/, PREDUC/7/,
1 RELDX/17/, STPPAR/5/
C/7
C PARAMETER (DSTNRM=2, F=10, F0=13, FDIF=11, NREDUC=6, PREDUC=7,
C 1 RELDX=17, STPPAR=5)
C/
C
C/6
DATA ZERO/0.D+0/
C/7
C PARAMETER (ZERO=0.D+0)
C/
C/6
DATA MODEL1(1)/4H /, MODEL1(2)/4H /, MODEL1(3)/4H /,
1 MODEL1(4)/4H /, MODEL1(5)/4H G /, MODEL1(6)/4H S /,
2 MODEL2(1)/4H G /, MODEL2(2)/4H S /, MODEL2(3)/4HG-S /,
3 MODEL2(4)/4HS-G /, MODEL2(5)/4H-S-G/, MODEL2(6)/4H-G-S/
C/7
C DATA MODEL1/' ',' ',' ',' ',' G ',' S '/,
C 1 MODEL2/' G ',' S ','G-S ','S-G ','-S-G','-G-S'/
C/
C
C------------------------------- BODY --------------------------------
C
PU = IV(PRUNIT)
IF (PU .EQ. 0) GO TO 999
IV1 = IV(1)
IF (IV1 .GT. 62) IV1 = IV1 - 51
OL = IV(OUTLEV)
ALG = IV(ALGSAV)
IF (IV1 .LT. 2 .OR. IV1 .GT. 15) GO TO 370
IF (OL .EQ. 0) GO TO 120
IF (IV1 .GE. 12) GO TO 120
IF (IV1 .EQ. 2 .AND. IV(NITER) .EQ. 0) GO TO 390
IF (IV1 .GE. 10 .AND. IV(PRNTIT) .EQ. 0) GO TO 120
IF (IV1 .GT. 2) GO TO 10
IV(PRNTIT) = IV(PRNTIT) + 1
IF (IV(PRNTIT) .LT. IABS(OL)) GO TO 999
10 NF = IV(NFCALL) - IABS(IV(NFCOV))
IV(PRNTIT) = 0
RELDF = ZERO
PRELDF = ZERO
OLDF = DMAX1(DABS(V(F0)), DABS(V(F)))
IF (OLDF .LE. ZERO) GO TO 20
RELDF = V(FDIF) / OLDF
PRELDF = V(PREDUC) / OLDF
20 IF (OL .GT. 0) GO TO 60
C
C *** PRINT SHORT SUMMARY LINE ***
C
IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 1) WRITE(PU,30)
30 FORMAT(/10H IT NF,6X,1HF,7X,5HRELDF,3X,6HPRELDF,3X,5HRELDX,
1 2X,13HMODEL STPPAR)
IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 2) WRITE(PU,40)
40 FORMAT(/11H IT NF,7X,1HF,8X,5HRELDF,4X,6HPRELDF,4X,5HRELDX,
1 3X,6HSTPPAR)
IV(NEEDHD) = 0
IF (ALG .EQ. 2) GO TO 50
M = IV(SUSED)
WRITE(PU,100) IV(NITER), NF, V(F), RELDF, PRELDF, V(RELDX),
1 MODEL1(M), MODEL2(M), V(STPPAR)
GO TO 120
C
50 WRITE(PU,110) IV(NITER), NF, V(F), RELDF, PRELDF, V(RELDX),
1 V(STPPAR)
GO TO 120
C
C *** PRINT LONG SUMMARY LINE ***
C
60 IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 1) WRITE(PU,70)
70 FORMAT(/11H IT NF,6X,1HF,7X,5HRELDF,3X,6HPRELDF,3X,5HRELDX,
1 2X,13HMODEL STPPAR,2X,6HD*STEP,2X,7HNPRELDF)
IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 2) WRITE(PU,80)
80 FORMAT(/11H IT NF,7X,1HF,8X,5HRELDF,4X,6HPRELDF,4X,5HRELDX,
1 3X,6HSTPPAR,3X,6HD*STEP,3X,7HNPRELDF)
IV(NEEDHD) = 0
NRELDF = ZERO
IF (OLDF .GT. ZERO) NRELDF = V(NREDUC) / OLDF
IF (ALG .EQ. 2) GO TO 90
M = IV(SUSED)
WRITE(PU,100) IV(NITER), NF, V(F), RELDF, PRELDF, V(RELDX),
1 MODEL1(M), MODEL2(M), V(STPPAR), V(DSTNRM), NRELDF
GO TO 120
C
90 WRITE(PU,110) IV(NITER), NF, V(F), RELDF, PRELDF,
1 V(RELDX), V(STPPAR), V(DSTNRM), NRELDF
100 FORMAT(I6,I5,D10.3,2D9.2,D8.1,A3,A4,2D8.1,D9.2)
110 FORMAT(I6,I5,D11.3,2D10.2,3D9.1,D10.2)
C
120 IF (IV(STATPR) .LT. 0) GO TO 430
GO TO (999, 999, 130, 150, 170, 190, 210, 230, 250, 270, 290, 310,
1 330, 350, 520), IV1
C
130 WRITE(PU,140)
140 FORMAT(/26H ***** X-CONVERGENCE *****)
GO TO 430
C
150 WRITE(PU,160)
160 FORMAT(/42H ***** RELATIVE FUNCTION CONVERGENCE *****)
GO TO 430
C
170 WRITE(PU,180)
180 FORMAT(/49H ***** X- AND RELATIVE FUNCTION CONVERGENCE *****)
GO TO 430
C
190 WRITE(PU,200)
200 FORMAT(/42H ***** ABSOLUTE FUNCTION CONVERGENCE *****)
GO TO 430
C
210 WRITE(PU,220)
220 FORMAT(/33H ***** SINGULAR CONVERGENCE *****)
GO TO 430
C
230 WRITE(PU,240)
240 FORMAT(/30H ***** FALSE CONVERGENCE *****)
GO TO 430
C
250 WRITE(PU,260)
260 FORMAT(/38H ***** FUNCTION EVALUATION LIMIT *****)
GO TO 430
C
270 WRITE(PU,280)
280 FORMAT(/28H ***** ITERATION LIMIT *****)
GO TO 430
C
290 WRITE(PU,300)
300 FORMAT(/18H ***** STOPX *****)
GO TO 430
C
310 WRITE(PU,320)
320 FORMAT(/44H ***** INITIAL F(X) CANNOT BE COMPUTED *****)
C
GO TO 390
C
330 WRITE(PU,340)
340 FORMAT(/37H ***** BAD PARAMETERS TO ASSESS *****)
GO TO 999
C
350 WRITE(PU,360)
360 FORMAT(/43H ***** GRADIENT COULD NOT BE COMPUTED *****)
IF (IV(NITER) .GT. 0) GO TO 480
GO TO 390
C
370 WRITE(PU,380) IV(1)
380 FORMAT(/14H ***** IV(1) =,I5,6H *****)
GO TO 999
C
C *** INITIAL CALL ON DITSUM ***
C
390 IF (IV(X0PRT) .NE. 0) WRITE(PU,400) (I, X(I), D(I), I = 1, P)
400 FORMAT(/23H I INITIAL X(I),8X,4HD(I)//(1X,I5,D17.6,D14.3))
IF (IV1 .GE. 12) GO TO 999
IV(NEEDHD) = 0
IV(PRNTIT) = 0
IF (OL .EQ. 0) GO TO 999
IF (OL .LT. 0 .AND. ALG .EQ. 1) WRITE(PU,30)
IF (OL .LT. 0 .AND. ALG .EQ. 2) WRITE(PU,40)
IF (OL .GT. 0 .AND. ALG .EQ. 1) WRITE(PU,70)
IF (OL .GT. 0 .AND. ALG .EQ. 2) WRITE(PU,80)
IF (ALG .EQ. 1) WRITE(PU,410) V(F)
IF (ALG .EQ. 2) WRITE(PU,420) V(F)
410 FORMAT(/11H 0 1,D10.3)
C365 FORMAT(/11H 0 1,E11.3)
420 FORMAT(/11H 0 1,D11.3)
GO TO 999
C
C *** PRINT VARIOUS INFORMATION REQUESTED ON SOLUTION ***
C
430 IV(NEEDHD) = 1
IF (IV(STATPR) .EQ. 0) GO TO 480
OLDF = DMAX1(DABS(V(F0)), DABS(V(F)))
PRELDF = ZERO
NRELDF = ZERO
IF (OLDF .LE. ZERO) GO TO 440
PRELDF = V(PREDUC) / OLDF
NRELDF = V(NREDUC) / OLDF
440 NF = IV(NFCALL) - IV(NFCOV)
NG = IV(NGCALL) - IV(NGCOV)
WRITE(PU,450) V(F), V(RELDX), NF, NG, PRELDF, NRELDF
450 FORMAT(/9H FUNCTION,D17.6,8H RELDX,D17.3/12H FUNC. EVALS,
1 I8,9X,11HGRAD. EVALS,I8/7H PRELDF,D16.3,6X,7HNPRELDF,D15.3)
C
IF (IV(NFCOV) .GT. 0) WRITE(PU,460) IV(NFCOV)
460 FORMAT(/1X,I4,50H EXTRA FUNC. EVALS FOR COVARIANCE AND DIAGNOST
1ICS.)
IF (IV(NGCOV) .GT. 0) WRITE(PU,470) IV(NGCOV)
470 FORMAT(1X,I4,50H EXTRA GRAD. EVALS FOR COVARIANCE AND DIAGNOSTI
1CS.)
C
480 IF (IV(SOLPRT) .EQ. 0) GO TO 999
IV(NEEDHD) = 1
WRITE(PU,490)
490 FORMAT(/22H I FINAL X(I),8X,4HD(I),10X,4HG(I)/)
DO 500 I = 1, P
WRITE(PU,510) I, X(I), D(I), G(I)
500 CONTINUE
510 FORMAT(1X,I5,D16.6,2D14.3)
GO TO 999
C
520 WRITE(PU,530)
530 FORMAT(/24H INCONSISTENT DIMENSIONS)
999 RETURN
C *** LAST CARD OF DITSUM FOLLOWS ***
END
SUBROUTINE DLITVM(N, X, L, Y)
save
C
C *** SOLVE (L**T)*X = Y, WHERE L IS AN N X N LOWER TRIANGULAR
C *** MATRIX STORED COMPACTLY BY ROWS. X AND Y MAY OCCUPY THE SAME
C *** STORAGE. ***
C
INTEGER N
DOUBLE PRECISION X(N), L(1), Y(N)
INTEGER I, II, IJ, IM1, I0, J, NP1
DOUBLE PRECISION XI, ZERO
C/6
DATA ZERO/0.D+0/
C/7
C PARAMETER (ZERO=0.D+0)
C/
C
DO 10 I = 1, N
10 X(I) = Y(I)
NP1 = N + 1
I0 = N*(N+1)/2
DO 30 II = 1, N
I = NP1 - II
XI = X(I)/L(I0)
X(I) = XI
IF (I .LE. 1) GO TO 999
I0 = I0 - I
IF (XI .EQ. ZERO) GO TO 30
IM1 = I - 1
DO 20 J = 1, IM1
IJ = I0 + J
X(J) = X(J) - XI*L(IJ)
20 CONTINUE
30 CONTINUE
999 RETURN
C *** LAST CARD OF DLITVM FOLLOWS ***
END
SUBROUTINE DLIVMU(N, X, L, Y)
save
C
C *** SOLVE L*X = Y, WHERE L IS AN N X N LOWER TRIANGULAR
C *** MATRIX STORED COMPACTLY BY ROWS. X AND Y MAY OCCUPY THE SAME
C *** STORAGE. ***
C
INTEGER N
DOUBLE PRECISION X(N), L(1), Y(N)
DOUBLE PRECISION DDOT
INTEGER I, J, K
DOUBLE PRECISION T, ZERO
C/6
DATA ZERO/0.D+0/
C/7
C PARAMETER (ZERO=0.D+0)
C/
C
DO 10 K = 1, N
IF (Y(K) .NE. ZERO) GO TO 20
X(K) = ZERO
10 CONTINUE
GO TO 999
20 J = K*(K+1)/2
X(K) = Y(K) / L(J)
IF (K .GE. N) GO TO 999
K = K + 1
DO 30 I = K, N
T = DDOT(I-1, L(J+1),1,X,1)
J = J + I
X(I) = (Y(I) - T)/L(J)
30 CONTINUE
999 RETURN
C *** LAST CARD OF DLIVMU FOLLOWS ***
END
SUBROUTINE DLTVMU(N, X, L, Y)
save
C
C *** COMPUTE X = (L**T)*Y, WHERE L IS AN N X N LOWER
C *** TRIANGULAR MATRIX STORED COMPACTLY BY ROWS. X AND Y MAY
C *** OCCUPY THE SAME STORAGE. ***
C
INTEGER N
DOUBLE PRECISION X(N), L(1), Y(N)
C DIMENSION L(N*(N+1)/2)
INTEGER I, IJ, I0, J
DOUBLE PRECISION YI, ZERO
C/6
DATA ZERO/0.D+0/
C/7
C 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 RETURN
C *** LAST CARD OF DLTVMU FOLLOWS ***
END
SUBROUTINE DLUPDT(BETA, GAMMA, L, LAMBDA, LPLUS, N, W, Z)
save
C
C *** COMPUTE LPLUS = SECANT UPDATE OF L ***
C
C *** 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)
C
C-------------------------- PARAMETER USAGE --------------------------
C
C 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 MAY
C 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 1
C CORRECTION TO L.
C Z (INPUT, DESTROYED ON OUTPUT) LEFT SINGULAR VECTOR OF RANK 1
C CORRECTION TO L.
C
C------------------------------- NOTES -------------------------------
C
C *** APPLICATION AND USAGE RESTRICTIONS ***
C
C THIS ROUTINE UPDATES THE CHOLESKY FACTOR L OF A SYMMETRIC
C POSITIVE DEFINITE MATRIX TO WHICH A SECANT UPDATE IS BEING
C APPLIED -- IT COMPUTES A CHOLESKY FACTOR LPLUS OF
C L * (I + Z*W**T) * (I + W*Z**T) * L**T. IT IS ASSUMED THAT W
C AND Z HAVE BEEN CHOSEN SO THAT THE UPDATED MATRIX IS STRICTLY
C POSITIVE DEFINITE.
C
C *** ALGORITHM NOTES ***
C
C 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 Q
C IS AN ORTHOGONAL MATRIX THAT MAKES THE RESULT LOWER TRIANGULAR.
C LPLUS MAY HAVE SOME NEGATIVE DIAGONAL ELEMENTS.
C
C *** REFERENCES ***
C
C 1. 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 (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 *** INTRINSIC FUNCTIONS ***
C/+
DOUBLE PRECISION DSQRT
C/
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, ZERO
C
C *** DATA INITIALIZATIONS ***
C
C/6
DATA ONE/1.D+0/, ZERO/0.D+0/
C/7
C PARAMETER (ONE=1.D+0, ZERO=0.D+0)
C/
C
C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++
C
NU = ONE
ETA = ZERO
IF (N .LE. 1) GO TO 30
NM1 = N - 1
C
C *** TEMPORARILY STORE S(J) = SUM OVER K = J+1 TO N OF W(K)**2 IN
C *** LAMBDA(J).
C
S = ZERO
DO 10 I = 1, NM1
J = N - I
S = S + W(J+1)**2
LAMBDA(J) = S
10 CONTINUE
C
C *** 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)
C
C *** 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 CONTINUE
C
999 RETURN
C *** LAST CARD OF DLUPDT FOLLOWS ***
END
SUBROUTINE DLVMUL(N, X, L, Y)
save
C
C *** COMPUTE X = L*Y, WHERE L IS AN N X N LOWER TRIANGULAR
C *** MATRIX STORED COMPACTLY BY ROWS. X AND Y MAY OCCUPY THE SAME
C *** STORAGE. ***
C
INTEGER N
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, ZERO
C/6
DATA ZERO/0.D+0/
C/7
C 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 RETURN
C *** LAST CARD OF DLVMUL FOLLOWS ***
END
DOUBLE PRECISION FUNCTION DNRM2(N,DX,INCX)
save
C***BEGIN PROLOGUE DNRM2
C***DATE WRITTEN 791001 (YYMMDD)
C***REVISION DATE 820801 (YYMMDD)
C***CATEGORY NO. D1A3B
C***KEYWORDS BLAS,DOUBLE PRECISION,EUCLIDEAN,L2,LENGTH,LINEAR ALGEBRA,
C NORM,VECTOR
C***AUTHOR LAWSON, C. L., (JPL)
C HANSON, R. J., (SNLA)
C KINCAID, D. R., (U. OF TEXAS)
C KROGH, F. T., (JPL)
C***PURPOSE Euclidean length (L2 norm) of d.p. vector
C***DESCRIPTION
C
C B L A S Subprogram
C Description of parameters
C
C --Input--
C N number of elements in input vector(s)
C DX double precision vector with N elements
C INCX storage spacing between elements of DX
C
C --Output--
C DNRM2 double precision result (zero if N .LE. 0)
C
C Euclidean norm of the N-vector stored in DX() with storage
C increment INCX .
C If N .LE. 0 return with result = 0.
C If N .GE. 1 then INCX must be .GE. 1
C
C C.L. Lawson, 1978 Jan 08
C
C Four phase method using two built-in constants that are
C hopefully applicable to all machines.
C CUTLO = maximum of DSQRT(U/EPS) over all known machines.
C CUTHI = minimum of DSQRT(V) over all known machines.
C where
C EPS = smallest no. such that EPS + 1. .GT. 1.
C U = smallest positive no. (underflow limit)
C V = largest no. (overflow limit)
C
C Brief outline of algorithm..
C
C Phase 1 scans zero components.
C move to phase 2 when a component is nonzero and .LE. CUTLO
C move to phase 3 when a component is .GT. CUTLO
C move to phase 4 when a component is .GE. CUTHI/M
C where M = N for X() real and M = 2*N for complex.
C
C Values for CUTLO and CUTHI..
C From the environmental parameters listed in the IMSL converter
C document the limiting values are as followS..
C CUTLO, S.P. U/EPS = 2**(-102) for Honeywell. Close seconds are
C Univac and DEC at 2**(-103)
C Thus CUTLO = 2**(-51) = 4.44089E-16
C CUTHI, S.P. V = 2**127 for Univac, Honeywell, and DEC.
C Thus CUTHI = 2**(63.5) = 1.30438E19
C CUTLO, D.P. U/EPS = 2**(-67) for Honeywell and DEC.
C Thus CUTLO = 2**(-33.5) = 8.23181D-11
C CUTHI, D.P. same as S.P. CUTHI = 1.30438D19
C DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 /
C DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 /
C***REFERENCES LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T.,
C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*,
C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL
C SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323
C***ROUTINES CALLED (NONE)
C***END PROLOGUE DNRM2
INTEGER NEXT
DOUBLE PRECISION DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE
DATA ZERO, ONE /0.0D0, 1.0D0/
C
DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 /
C***FIRST EXECUTABLE STATEMENT DNRM2
IF(N .GT. 0) GO TO 10
DNRM2 = ZERO
GO TO 300
C
10 ASSIGN 30 TO NEXT
SUM = ZERO
NN = N * INCX
C BEGIN MAIN LOOP
I = 1
20 GO TO NEXT,(30, 50, 70, 110)
30 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85
ASSIGN 50 TO NEXT
XMAX = ZERO
C
C PHASE 1. SUM IS ZERO
C
50 IF( DX(I) .EQ. ZERO) GO TO 200
IF( DABS(DX(I)) .GT. CUTLO) GO TO 85
C
C PREPARE FOR PHASE 2.
ASSIGN 70 TO NEXT
GO TO 105
C
C PREPARE FOR PHASE 4.
C
100 I = J
ASSIGN 110 TO NEXT
SUM = (SUM / DX(I)) / DX(I)
105 XMAX = DABS(DX(I))
GO TO 115
C
C PHASE 2. SUM IS SMALL.
C SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
C
70 IF( DABS(DX(I)) .GT. CUTLO ) GO TO 75
C
C COMMON CODE FOR PHASES 2 AND 4.
C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW.
C
110 IF( DABS(DX(I)) .LE. XMAX ) GO TO 115
SUM = ONE + SUM * (XMAX / DX(I))**2
XMAX = DABS(DX(I))
GO TO 200
C
115 SUM = SUM + (DX(I)/XMAX)**2
GO TO 200
C
C
C PREPARE FOR PHASE 3.
C
75 SUM = (SUM * XMAX) * XMAX
C
C
C FOR REAL OR D.P. SET HITEST = CUTHI/N
C FOR COMPLEX SET HITEST = CUTHI/(2*N)
C
85 HITEST = CUTHI/FLOAT( N )
C
C PHASE 3. SUM IS MID-RANGE. NO SCALING.
C
DO 95 J =I,NN,INCX
IF(DABS(DX(J)) .GE. HITEST) GO TO 100
95 SUM = SUM + DX(J)**2
DNRM2 = DSQRT( SUM )
GO TO 300
C
200 CONTINUE
I = I + INCX
IF ( I .LE. NN ) GO TO 20
C
C END OF MAIN LOOP.
C
C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
C
DNRM2 = XMAX * DSQRT(SUM)
300 CONTINUE
RETURN
END
SUBROUTINE DPARCK(ALG, D, IV, LIV, LV, N, V)
save
C
C *** CHECK ***SOL (VERSION 2.3) PARAMETERS, PRINT CHANGED VALUES ***
C
C *** ALG = 1 FOR REGRESSION, ALG = 2 FOR GENERAL UNCONSTRAINED OPT.
C
INTEGER ALG, LIV, LV, N
INTEGER IV(LIV)
DOUBLE PRECISION D(N), V(LV)
C
EXTERNAL DVDFLT
DOUBLE PRECISION D1MACH
C DVDFLT -- SUPPLIES DEFAULT PARAMETER VALUES TO V ALONE.
C/+
INTEGER MAX0
C/
C
C *** LOCAL VARIABLES ***
C
INTEGER I, II, IV1, J, K, L, M, MIV1, MIV2, NDFALT, PARSV1, PU
INTEGER IJMP, JLIM(2), MINIV(2), NDFLT(2)
C/6
INTEGER VARNM(2), SH(2)
REAL CNGD(3), DFLT(3), VN(2,34), WHICH(3)
C/7
C CHARACTER*1 VARNM(2), SH(2)
C CHARACTER*4 CNGD(3), DFLT(3), VN(2,34), WHICH(3)
C/
DOUBLE PRECISION BIG, MACHEP, TINY, VK, VM(34), VX(34), ZERO
C
C *** IV AND V SUBSCRIPTS ***
C
INTEGER ALGSAV, DINIT, DTYPE, DTYPE0, EPSLON, INITS, IVNEED,
1 LASTIV, LASTV, LMAT, NEXTIV, NEXTV, NVDFLT, OLDN,
2 PARPRT, PARSAV, PERM, PRUNIT, VNEED
C
C
C/6
DATA ALGSAV/51/, DINIT/38/, DTYPE/16/, DTYPE0/54/, EPSLON/19/,
1 INITS/25/, IVNEED/3/, LASTIV/44/, LASTV/45/, LMAT/42/,
2 NEXTIV/46/, NEXTV/47/, NVDFLT/50/, OLDN/38/, PARPRT/20/,
3 PARSAV/49/, PERM/58/, PRUNIT/21/, VNEED/4/
C/7
C PARAMETER (ALGSAV=51, DINIT=38, DTYPE=16, DTYPE0=54, EPSLON=19,
C 1 INITS=25, IVNEED=3, LASTIV=44, LASTV=45, LMAT=42,
C 2 NEXTIV=46, NEXTV=47, NVDFLT=50, OLDN=38, PARPRT=20,
C 3 PARSAV=49, PERM=58, PRUNIT=21, VNEED=4)
C SAVE BIG, MACHEP, TINY
C/
C
DATA BIG/0.D+0/, MACHEP/-1.D+0/, TINY/1.D+0/, ZERO/0.D+0/
C/6
DATA VN(1,1),VN(2,1)/4HEPSL,4HON../
DATA VN(1,2),VN(2,2)/4HPHMN,4HFC../
DATA VN(1,3),VN(2,3)/4HPHMX,4HFC../
DATA VN(1,4),VN(2,4)/4HDECF,4HAC../
DATA VN(1,5),VN(2,5)/4HINCF,4HAC../
DATA VN(1,6),VN(2,6)/4HRDFC,4HMN../
DATA VN(1,7),VN(2,7)/4HRDFC,4HMX../
DATA VN(1,8),VN(2,8)/4HTUNE,4HR1../
DATA VN(1,9),VN(2,9)/4HTUNE,4HR2../
DATA VN(1,10),VN(2,10)/4HTUNE,4HR3../
DATA VN(1,11),VN(2,11)/4HTUNE,4HR4../
DATA VN(1,12),VN(2,12)/4HTUNE,4HR5../
DATA VN(1,13),VN(2,13)/4HAFCT,4HOL../
DATA VN(1,14),VN(2,14)/4HRFCT,4HOL../
DATA VN(1,15),VN(2,15)/4HXCTO,4HL.../
DATA VN(1,16),VN(2,16)/4HXFTO,4HL.../
DATA VN(1,17),VN(2,17)/4HLMAX,4H0.../
DATA VN(1,18),VN(2,18)/4HLMAX,4HS.../
DATA VN(1,19),VN(2,19)/4HSCTO,4HL.../
DATA VN(1,20),VN(2,20)/4HDINI,4HT.../
DATA VN(1,21),VN(2,21)/4HDTIN,4HIT../
DATA VN(1,22),VN(2,22)/4HD0IN,4HIT../
DATA VN(1,23),VN(2,23)/4HDFAC,4H..../
DATA VN(1,24),VN(2,24)/4HDLTF,4HDC../
DATA VN(1,25),VN(2,25)/4HDLTF,4HDJ../
DATA VN(1,26),VN(2,26)/4HDELT,4HA0../
DATA VN(1,27),VN(2,27)/4HFUZZ,4H..../
DATA VN(1,28),VN(2,28)/4HRLIM,4HIT../
DATA VN(1,29),VN(2,29)/4HCOSM,4HIN../
DATA VN(1,30),VN(2,30)/4HHUBE,4HRC../
DATA VN(1,31),VN(2,31)/4HRSPT,4HOL../
DATA VN(1,32),VN(2,32)/4HSIGM,4HIN../
DATA VN(1,33),VN(2,33)/4HETA0,4H..../
DATA VN(1,34),VN(2,34)/4HBIAS,4H..../
C/7
C DATA VN(1,1),VN(2,1)/'EPSL','ON..'/
C DATA VN(1,2),VN(2,2)/'PHMN','FC..'/
C DATA VN(1,3),VN(2,3)/'PHMX','FC..'/
C DATA VN(1,4),VN(2,4)/'DECF','AC..'/
C DATA VN(1,5),VN(2,5)/'INCF','AC..'/
C DATA VN(1,6),VN(2,6)/'RDFC','MN..'/
C DATA VN(1,7),VN(2,7)/'RDFC','MX..'/
C DATA VN(1,8),VN(2,8)/'TUNE','R1..'/
C DATA VN(1,9),VN(2,9)/'TUNE','R2..'/
C DATA VN(1,10),VN(2,10)/'TUNE','R3..'/
C DATA VN(1,11),VN(2,11)/'TUNE','R4..'/
C DATA VN(1,12),VN(2,12)/'TUNE','R5..'/
C DATA VN(1,13),VN(2,13)/'AFCT','OL..'/
C DATA VN(1,14),VN(2,14)/'RFCT','OL..'/
C DATA VN(1,15),VN(2,15)/'XCTO','L...'/
C DATA VN(1,16),VN(2,16)/'XFTO','L...'/
C DATA VN(1,17),VN(2,17)/'LMAX','0...'/
C DATA VN(1,18),VN(2,18)/'LMAX','S...'/
C DATA VN(1,19),VN(2,19)/'SCTO','L...'/
C DATA VN(1,20),VN(2,20)/'DINI','T...'/
C DATA VN(1,21),VN(2,21)/'DTIN','IT..'/
C DATA VN(1,22),VN(2,22)/'D0IN','IT..'/
C DATA VN(1,23),VN(2,23)/'DFAC','....'/
C DATA VN(1,24),VN(2,24)/'DLTF','DC..'/
C DATA VN(1,25),VN(2,25)/'DLTF','DJ..'/
C DATA VN(1,26),VN(2,26)/'DELT','A0..'/
C DATA VN(1,27),VN(2,27)/'FUZZ','....'/
C DATA VN(1,28),VN(2,28)/'RLIM','IT..'/
C DATA VN(1,29),VN(2,29)/'COSM','IN..'/
C DATA VN(1,30),VN(2,30)/'HUBE','RC..'/
C DATA VN(1,31),VN(2,31)/'RSPT','OL..'/
C DATA VN(1,32),VN(2,32)/'SIGM','IN..'/
C DATA VN(1,33),VN(2,33)/'ETA0','....'/
C DATA VN(1,34),VN(2,34)/'BIAS','....'/
C/
C
DATA VM(1)/1.0D-3/, VM(2)/-0.99D+0/, VM(3)/1.0D-3/, VM(4)/1.0D-2/,
1 VM(5)/1.2D+0/, VM(6)/1.D-2/, VM(7)/1.2D+0/, VM(8)/0.D+0/,
2 VM(9)/0.D+0/, VM(10)/1.D-3/, VM(11)/-1.D+0/, VM(15)/0.D+0/,
3 VM(16)/0.D+0/, VM(19)/0.D+0/, VM(20)/-10.D+0/, VM(21)/0.D+0/,
4 VM(22)/0.D+0/, VM(23)/0.D+0/, VM(27)/1.01D+0/,
5 VM(28)/1.D+10/, VM(30)/0.D+0/, VM(31)/0.D+0/, VM(32)/0.D+0/,
6 VM(34)/0.D+0/
DATA VX(1)/0.9D+0/, VX(2)/-1.D-3/, VX(3)/1.D+1/, VX(4)/0.8D+0/,
1 VX(5)/1.D+2/, VX(6)/0.8D+0/, VX(7)/1.D+2/, VX(8)/0.5D+0/,
2 VX(9)/0.5D+0/, VX(10)/1.D+0/, VX(11)/1.D+0/, VX(14)/0.1D+0/,
3 VX(15)/1.D+0/, VX(16)/1.D+0/, VX(19)/1.D+0/, VX(23)/1.D+0/,
4 VX(24)/1.D+0/, VX(25)/1.D+0/, VX(26)/1.D+0/, VX(27)/1.D+10/,
5 VX(29)/1.D+0/, VX(31)/1.D+0/, VX(32)/1.D+0/, VX(33)/1.D+0/,
6 VX(34)/1.D+0/
C
C/6
DATA VARNM(1)/1HP/, VARNM(2)/1HN/, SH(1)/1HS/, SH(2)/1HH/
DATA CNGD(1),CNGD(2),CNGD(3)/4H---C,4HHANG,4HED V/,
1 DFLT(1),DFLT(2),DFLT(3)/4HNOND,4HEFAU,4HLT V/
C/7
C DATA VARNM(1)/'P'/, VARNM(2)/'N'/, SH(1)/'S'/, SH(2)/'H'/
C DATA CNGD(1),CNGD(2),CNGD(3)/'---C','HANG','ED V'/,
C 1 DFLT(1),DFLT(2),DFLT(3)/'NOND','EFAU','LT V'/
C/
DATA IJMP/33/, JLIM(1)/0/, JLIM(2)/24/, NDFLT(1)/32/, NDFLT(2)/25/
DATA MINIV(1)/80/, MINIV(2)/59/
C
C............................... BODY ................................
C
IF (ALG .LT. 1 .OR. ALG .GT. 2) GO TO 330
IF (IV(1) .EQ. 0) CALL DDEFLT(ALG, IV, LIV, LV, V)
PU = IV(PRUNIT)
MIV1 = MINIV(ALG)
IF (PERM .LE. LIV) MIV1 = MAX0(MIV1, IV(PERM) - 1)
IF (IVNEED .LE. LIV) MIV2 = MIV1 + MAX0(IV(IVNEED), 0)
IF (LASTIV .LE. LIV) IV(LASTIV) = MIV2
IF (LIV .LT. MIV1) GO TO 290
IV(IVNEED) = 0
IV(LASTV) = MAX0(IV(VNEED), 0) + IV(LMAT) - 1
IF (LIV .LT. MIV2) GO TO 290
IF (LV .LT. IV(LASTV)) GO TO 310
IV(VNEED) = 0
IF (ALG .EQ. IV(ALGSAV)) GO TO 20
IF (PU .NE. 0) WRITE(PU,10) ALG, IV(ALGSAV)
10 FORMAT(/' THE FIRST PARAMETER TO DDEFLT SHOULD BE',I3,
1 12H RATHER THAN,I3)
IV(1) = 82
GO TO 999
20 IV1 = IV(1)
IF (IV1 .LT. 12 .OR. IV1 .GT. 14) GO TO 50
IF (N .GE. 1) GO TO 40
IV(1) = 81
IF (PU .EQ. 0) GO TO 999
WRITE(PU,30) VARNM(ALG), N
30 FORMAT(/8H /// BAD,A1,2H =,I5)
GO TO 999
40 IF (IV1 .NE. 14) IV(NEXTIV) = IV(PERM)
IF (IV1 .NE. 14) IV(NEXTV) = IV(LMAT)
IF (IV1 .EQ. 13) GO TO 999
K = IV(PARSAV) - EPSLON
CALL DVDFLT(ALG, LV-K, V(K+1))
IV(DTYPE0) = 2 - ALG
IV(OLDN) = N
WHICH(1) = DFLT(1)
WHICH(2) = DFLT(2)
WHICH(3) = DFLT(3)
GO TO 100
50 IF (N .EQ. IV(OLDN)) GO TO 70
IV(1) = 17
IF (PU .EQ. 0) GO TO 999
WRITE(PU,60) VARNM(ALG), IV(OLDN), N
60 FORMAT(/5H /// ,1A1,14H CHANGED FROM ,I5,4H TO ,I5)
GO TO 999
C
70 IF (IV1 .LE. 11 .AND. IV1 .GE. 1) GO TO 90
IV(1) = 80
IF (PU .NE. 0) WRITE(PU,80) IV1
80 FORMAT(/13H /// IV(1) =,I5,28H SHOULD BE BETWEEN 0 AND 14.)
GO TO 999
C
90 WHICH(1) = CNGD(1)
WHICH(2) = CNGD(2)
WHICH(3) = CNGD(3)
C
100 IF (IV1 .EQ. 14) IV1 = 12
IF (BIG .GT. TINY) GO TO 110
TINY = D1MACH(1)
MACHEP = D1MACH(4)
BIG = D1MACH(2)
VM(12) = MACHEP
VX(12) = BIG
VM(13) = TINY
VX(13) = BIG
VM(14) = MACHEP
VM(17) = TINY
VX(17) = BIG
VM(18) = TINY
VX(18) = BIG
VX(20) = BIG
VX(21) = BIG
VX(22) = BIG
VM(24) = MACHEP
VM(25) = MACHEP
VM(26) = MACHEP
VX(28) = DSQRT(D1MACH(2))*16.
VM(29) = MACHEP
VX(30) = BIG
VM(33) = MACHEP
110 M = 0
I = 1
J = JLIM(ALG)
K = EPSLON
NDFALT = NDFLT(ALG)
DO 140 L = 1, NDFALT
VK = V(K)
IF (VK .GE. VM(I) .AND. VK .LE. VX(I)) GO TO 130
M = K
IF (PU .NE. 0) WRITE(PU,120) VN(1,I), VN(2,I), K, VK,
1 VM(I), VX(I)
120 FORMAT(/6H /// ,2A4,5H.. V(,I2,3H) =,D11.3,7H SHOULD,
1 11H BE BETWEEN,D11.3,4H AND,D11.3)
130 K = K + 1
I = I + 1
IF (I .EQ. J) I = IJMP
140 CONTINUE
C
IF (IV(NVDFLT) .EQ. NDFALT) GO TO 160
IV(1) = 51
IF (PU .EQ. 0) GO TO 999
WRITE(PU,150) IV(NVDFLT), NDFALT
150 FORMAT(/13H IV(NVDFLT) =,I5,13H RATHER THAN ,I5)
GO TO 999
160 IF ((IV(DTYPE) .GT. 0 .OR. V(DINIT) .GT. ZERO) .AND. IV1 .EQ. 12)
1 GO TO 190
DO 180 I = 1, N
IF (D(I) .GT. ZERO) GO TO 180
M = 18
IF (PU .NE. 0) WRITE(PU,170) I, D(I)
170 FORMAT(/8H /// D(,I3,3H) =,D11.3,19H SHOULD BE POSITIVE)
180 CONTINUE
190 IF (M .EQ. 0) GO TO 200
IV(1) = M
GO TO 999
C
200 IF (PU .EQ. 0 .OR. IV(PARPRT) .EQ. 0) GO TO 999
IF (IV1 .NE. 12 .OR. IV(INITS) .EQ. ALG-1) GO TO 220
M = 1
WRITE(PU,210) SH(ALG), IV(INITS)
210 FORMAT(/22H NONDEFAULT VALUES..../5H INIT,A1,14H..... IV(25) =,
1 I3)
220 IF (IV(DTYPE) .EQ. IV(DTYPE0)) GO TO 240
IF (M .EQ. 0) WRITE(PU,250) WHICH
M = 1
WRITE(PU,230) IV(DTYPE)
230 FORMAT(20H DTYPE..... IV(16) =,I3)
240 I = 1
J = JLIM(ALG)
K = EPSLON
L = IV(PARSAV)
NDFALT = NDFLT(ALG)
DO 280 II = 1, NDFALT
IF (V(K) .EQ. V(L)) GO TO 270
IF (M .EQ. 0) WRITE(PU,250) WHICH
250 FORMAT(/1H ,3A4,9HALUES..../)
M = 1
WRITE(PU,260) VN(1,I), VN(2,I), K, V(K)
260 FORMAT(1X,2A4,5H.. V(,I2,3H) =,D15.7)
270 K = K + 1
L = L + 1
I = I + 1
IF (I .EQ. J) I = IJMP
280 CONTINUE
C
IV(DTYPE0) = IV(DTYPE)
PARSV1 = IV(PARSAV)
CALL DCOPY(IV(NVDFLT), V(EPSLON),1,V(PARSV1),1)
GO TO 999
C
290 IV(1) = 15
IF (PU .EQ. 0) GO TO 999
WRITE(PU,300) LIV, MIV2
300 FORMAT(/10H /// LIV =,I5,17H MUST BE AT LEAST,I5)
IF (LIV .LT. MIV1) GO TO 999
IF (LV .LT. IV(LASTV)) GO TO 310
GO TO 999
C
310 IV(1) = 16
IF (PU .EQ. 0) GO TO 999
WRITE(PU,320) LV, IV(LASTV)
320 FORMAT(/9H /// LV =,I5,17H MUST BE AT LEAST,I5)
GO TO 999
C
330 IV(1) = 67
C
999 RETURN
C *** LAST CARD OF DPARCK FOLLOWS ***
END
DOUBLE PRECISION FUNCTION DRELST(P, D, X, X0)
save
C
C *** COMPUTE AND RETURN RELATIVE DIFFERENCE BETWEEN X AND X0 ***
C *** NL2SOL VERSION 2.2 ***
C
INTEGER P
DOUBLE PRECISION D(P), X(P), X0(P)
C/+
DOUBLE PRECISION DABS
C/
INTEGER I
DOUBLE PRECISION EMAX, T, XMAX, ZERO
C/6
DATA ZERO/0.D+0/
C/7
C PARAMETER (ZERO=0.D+0)
C/
C
EMAX = ZERO
XMAX = ZERO
DO 10 I = 1, P
T = DABS(D(I) * (X(I) - X0(I)))
IF (EMAX .LT. T) EMAX = T
T = D(I) * (DABS(X(I)) + DABS(X0(I)))
IF (XMAX .LT. T) XMAX = T
10 CONTINUE
DRELST = ZERO
IF (XMAX .GT. ZERO) DRELST = EMAX / XMAX
999 RETURN
C *** LAST CARD OF DRELST FOLLOWS ***
END
LOGICAL FUNCTION DSTOPX(IDUMMY)
save
C *****PARAMETERS...
INTEGER IDUMMY
C
C ..................................................................
C
C *****PURPOSE...
C THIS FUNCTION MAY SERVE AS THE DSTOPX (ASYNCHRONOUS INTERRUPTION)
C FUNCTION FOR THE NL2SOL (NONLINEAR LEAST-SQUARES) PACKAGE AT
C THOSE INSTALLATIONS WHICH DO NOT WISH TO IMPLEMENT A
C DYNAMIC DSTOPX.
C
C *****ALGORITHM NOTES...
C AT INSTALLATIONS WHERE THE NL2SOL SYSTEM IS USED
C INTERACTIVELY, THIS DUMMY DSTOPX SHOULD BE REPLACED BY A
C FUNCTION THAT RETURNS .TRUE. IF AND ONLY IF THE INTERRUPT
C (BREAK) KEY HAS BEEN PRESSED SINCE THE LAST CALL ON DSTOPX.
C
C ..................................................................
C
DSTOPX = .FALSE.
RETURN
END
SUBROUTINE FDUMP
save
C***BEGIN PROLOGUE FDUMP
C***DATE WRITTEN 790801 (YYMMDD)
C***REVISION DATE 820801 (YYMMDD)
C***CATEGORY NO. Z
C***KEYWORDS ERROR,XERROR PACKAGE
C***AUTHOR JONES, R. E., (SNLA)
C***PURPOSE Symbolic dump (should be locally written).
C***DESCRIPTION
C ***Note*** Machine Dependent Routine
C FDUMP is intended to be replaced by a locally written
C version which produces a symbolic dump. Failing this,
C it should be replaced by a version which prints the
C subprogram nesting list. Note that this dump must be
C printed on each of up to five files, as indicated by the
C XGETUA routine. See XSETUA and XGETUA for details.
C
C Written by Ron Jones, with SLATEC Common Math Library Subcommittee
C Latest revision --- 23 May 1979
C***ROUTINES CALLED (NONE)
C***END PROLOGUE FDUMP
C***FIRST EXECUTABLE STATEMENT FDUMP
RETURN
END
FUNCTION J4SAVE(IWHICH,IVALUE,ISET)
save
C***BEGIN PROLOGUE J4SAVE
C***REFER TO XERROR
C Abstract
C J4SAVE saves and recalls several global variables needed
C by the library error handling routines.
C
C Description of Parameters
C --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 error
C messages are to be sent. (0 means use standard.)
C = 4 Refers to the maximum number of times any
C message is to be printed (as set by XERMAX).
C = 5 Refers to the total number of units to which
C each error message is to be written.
C = 6 Refers to the 2nd unit for error messages
C = 7 Refers to the 3rd unit for error messages
C = 8 Refers to the 4th unit for error messages
C = 9 Refers to the 5th unit for error messages
C 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 BE
C given the value, IVALUE. If ISET=.FALSE., the
C IWHICH-th parameter will be unchanged, and IVALUE
C is a dummy parameter.
C --Output--
C The (old) value of the IWHICH-th parameter will be returned
C in the function value, J4SAVE.
C
C Written by Ron Jones, with SLATEC Common Math Library Subcommittee
C Adapted from Bell Laboratories PORT Library Error Handler
C Latest revision --- 23 MAY 1979
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 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)
save
C***BEGIN PROLOGUE XERABT
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 Aborts program execution and prints error message.
C***DESCRIPTION
C Abstract
C ***Note*** machine dependent routine
C XERABT aborts the execution of the program.
C The error message causing the abort is given in the calling
C sequence, in case one needs it for printing on a dayfile,
C for example.
C
C Description of Parameters
C MESSG and NMESSG are as in XERROR, except that NMESSG may
C be zero, in which case no message is being supplied.
C
C Written by Ron Jones, with SLATEC Common Math Library Subcommittee
C Latest revision --- 19 MAR 1980
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 XERABT
CHARACTER*(*) MESSG
EXTERNAL ERROR
C***FIRST EXECUTABLE STATEMENT XERABT
CALL ERROR('dsumsl(), dsmsno() aborted\n')
END
SUBROUTINE XERCTL(MESSG1,NMESSG,NERR,LEVEL,KONTRL)
save
C***BEGIN PROLOGUE XERCTL
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 Allows user control over handling of individual errors.
C***DESCRIPTION
C Abstract
C Allows user control over handling of individual errors.
C Just after each message is recorded, but before it is
C processed any further (i.e., before it is printed or
C a decision to abort is made), a call is made to XERCTL.
C If the user has provided his own version of XERCTL, he
C can then override the value of KONTROL used in processing
C 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, except
C 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.
C
C Description of Parameters
C
C --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 set
C by a call to XSETF.
C
C --Output--
C KONTRL - the new value of KONTRL. If KONTRL is not
C defined, it will remain at its original value.
C This changed value of control affects only
C 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 MESSG1
C***FIRST EXECUTABLE STATEMENT XERCTL
RETURN
END
SUBROUTINE XERPRT(MESSG,NMESSG)
save
C***BEGIN PROLOGUE XERPRT
C***DATE WRITTEN 790801 (YYMMDD)
C***REVISION DATE 820801 (YYMMDD)
C***CATEGORY NO. Z
C***KEYWORDS ERROR,XERROR PACKAGE
C***AUTHOR JONES, R. E., (SNLA)
C***PURPOSE Prints error messages.
C***DESCRIPTION
C Abstract
C Print the Hollerith message in MESSG, of length NMESSG,
C on each file indicated by XGETUA.
C Latest revision --- 19 MAR 1980
C***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,XGETUA
C***END PROLOGUE XERPRT
INTEGER LUN(5)
CHARACTER*(*) MESSG
C OBTAIN UNIT NUMBERS AND WRITE LINE TO EACH UNIT
C***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)') MESSG(ICHAR:LAST)
10 CONTINUE
20 CONTINUE
RETURN
END
SUBROUTINE DSMSNO(N, D, X, CALCF, IV, LIV, LV, V,
1 UIPARM, URPARM, UFPARM)
save
C
C *** MINIMIZE GENERAL UNCONSTRAINED OBJECTIVE FUNCTION USING
C *** FINITE-DIFFERENCE GRADIENTS AND SECANT HESSIAN APPROXIMATIONS.
C
INTEGER N, LIV, LV
INTEGER IV(LIV), UIPARM(1)
DOUBLE PRECISION D(N), X(N), V(LV), URPARM(1)
C DIMENSION V(77 + N*(N+17)/2), UIPARM(*), URPARM(*)
EXTERNAL CALCF, UFPARM
C
C *** PURPOSE ***
C
C THIS ROUTINE INTERACTS WITH SUBROUTINE DSNOIT 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 *** PARAMETERS ***
C
C THE PARAMETERS FOR DSMSNO ARE THE SAME AS THOSE FOR DSUMSL
C (WHICH SEE), EXCEPT THAT CALCG IS OMITTED. INSTEAD OF CALLING
C CALCG TO OBTAIN THE GRADIENT OF THE OBJECTIVE FUNCTION AT X,
C DSMSNO CALLS DSGRD2, WHICH COMPUTES AN APPROXIMATION TO THE
C GRADIENT BY FINITE (FORWARD AND CENTRAL) DIFFERENCES USING THE
C METHOD OF REF. 1. THE FOLLOWING INPUT COMPONENT IS OF INTEREST
C IN THIS REGARD (AND IS NOT DESCRIBED IN DSUMSL).
C
C V(ETA0)..... V(42) IS AN ESTIMATED BOUND ON THE RELATIVE ERROR IN THE
C OBJECTIVE FUNCTION VALUE COMPUTED BY CALCF...
C (TRUE VALUE) = (COMPUTED VALUE) * (1 + E),
C WHERE ABS(E) .LE. V(ETA0). DEFAULT = MACHEP * 10**3,
C WHERE MACHEP IS THE UNIT ROUNDOFF.
C
C THE OUTPUT VALUES IV(NFCALL) AND IV(NGCALL) HAVE DIFFERENT
C MEANINGS FOR DSMSNO THAN FOR DSUMSL...
C
C IV(NFCALL)... IV(6) IS THE NUMBER OF CALLS SO FAR MADE ON CALCF (I.E.,
C FUNCTION EVALUATIONS) EXCLUDING THOSE MADE ONLY FOR
C COMPUTING GRADIENTS. THE INPUT VALUE IV(MXFCAL) IS A
C LIMIT ON IV(NFCALL).
C IV(NGCALL)... IV(30) IS THE NUMBER OF FUNCTION EVALUATIONS MADE ONLY
C FOR COMPUTING GRADIENTS. THE TOTAL NUMBER OF FUNCTION
C EVALUATIONS IS THUS IV(NFCALL) + IV(NGCALL).
C
C *** REFERENCES ***
C
C 1. STEWART, G.W. (1967), A MODIFICATION OF DAVIDON*S MINIMIZATION
C METHOD TO ACCEPT DIFFERENCE APPROXIMATIONS OF DERIVATIVES,
C J. ASSOC. COMPUT. MACH. 14, PP. 72-83.
C.
C *** GENERAL ***
C
C CODED BY DAVID M. GAY (WINTER 1980). REVISED SEPT. 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 DSNOIT
C
C DSNOIT.... OVERSEES COMPUTATION OF FINITE-DIFFERENCE GRADIENT AND
C CALLS DSUMIT TO CARRY OUT DSUMSL ALGORITHM.
C
INTEGER NF
DOUBLE PRECISION FX
C
C *** SUBSCRIPTS FOR IV ***
C
INTEGER NFCALL, TOOBIG
C
C/6
DATA NFCALL/6/, TOOBIG/2/
C/7
C PARAMETER (NFCALL=6, TOOBIG=2)
C/
C
C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++
C
10 CALL DSNOIT(D, FX, IV, LIV, LV, N, V, X)
IF (IV(1) .GT. 2) GO TO 999
C
C *** COMPUTE FUNCTION ***
C
NF = IV(NFCALL)
CALL CALCF(N, X, NF, FX, UIPARM, URPARM, UFPARM)
IF (NF .LE. 0) IV(TOOBIG) = 1
GO TO 10
C
C
999 RETURN
C *** LAST CARD OF DSMSNO FOLLOWS ***
END
SUBROUTINE DSNOIT(D, FX, IV, LIV, LV, N, V, X)
save
C
C *** ITERATION DRIVER FOR DSMSNO...
C *** MINIMIZE GENERAL UNCONSTRAINED OBJECTIVE FUNCTION USING
C *** FINITE-DIFFERENCE GRADIENTS AND SECANT HESSIAN APPROXIMATIONS.
C
INTEGER LIV, LV, N
INTEGER IV(LIV)
DOUBLE PRECISION D(N), FX, X(N), V(LV)
C DIMENSION V(77 + N*(N+17)/2)
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 FX = F(X) COMPUTED BY THE CALLER. (OFTEN
C THE X* FOUND IS A LOCAL MINIMIZER RATHER THAN A GLOBAL ONE.)
C
C *** PARAMETERS ***
C
C THE PARAMETERS FOR DSNOIT ARE THE SAME AS THOSE FOR DSUMSL
C (WHICH SEE), EXCEPT THAT CALCF, CALCG, UIPARM, URPARM, AND UFPARM
C ARE OMITTED, AND A PARAMETER FX FOR THE OBJECTIVE FUNCTION
C VALUE AT X IS ADDED. INSTEAD OF CALLING CALCG TO OBTAIN THE
C GRADIENT OF THE OBJECTIVE FUNCTION AT X, DSNOIT CALLS DSGRD2,
C WHICH COMPUTES AN APPROXIMATION TO THE GRADIENT BY FINITE
C (FORWARD AND CENTRAL) DIFFERENCES USING THE METHOD OF REF. 1.
C THE FOLLOWING INPUT COMPONENT IS OF INTEREST IN THIS REGARD
C (AND IS NOT DESCRIBED IN DSUMSL).
C
C V(ETA0)..... V(42) IS AN ESTIMATED BOUND ON THE RELATIVE ERROR IN THE
C OBJECTIVE FUNCTION VALUE COMPUTED BY CALCF...
C (TRUE VALUE) = (COMPUTED VALUE) * (1 + E),
C WHERE ABS(E) .LE. V(ETA0). DEFAULT = MACHEP * 10**3,
C WHERE MACHEP IS THE UNIT ROUNDOFF.
C
C THE OUTPUT VALUES IV(NFCALL) AND IV(NGCALL) HAVE DIFFERENT
C MEANINGS FOR DSMSNO THAN FOR DSUMSL...
C
C IV(NFCALL)... IV(6) IS THE NUMBER OF CALLS SO FAR MADE ON CALCF (I.E.,
C FUNCTION EVALUATIONS) EXCLUDING THOSE MADE ONLY FOR
C COMPUTING GRADIENTS. THE INPUT VALUE IV(MXFCAL) IS A
C LIMIT ON IV(NFCALL).
C IV(NGCALL)... IV(30) IS THE NUMBER OF FUNCTION EVALUATIONS MADE ONLY
C FOR COMPUTING GRADIENTS. THE TOTAL NUMBER OF FUNCTION
C EVALUATIONS IS THUS IV(NFCALL) + IV(NGCALL).
C
C *** REFERENCES ***
C
C 1. STEWART, G.W. (1967), A MODIFICATION OF DAVIDON*S MINIMIZATION
C METHOD TO ACCEPT DIFFERENCE APPROXIMATIONS OF DERIVATIVES,
C J. ASSOC. COMPUT. MACH. 14, PP. 72-83.
C.
C *** GENERAL ***
C
C CODED BY DAVID M. GAY (AUGUST 1982).
C
C---------------------------- DECLARATIONS ---------------------------
C
EXTERNAL DDEFLT, DSGRD2, DSUMIT, DVSCPY
DOUBLE PRECISION DDOT
C
C DDEFLT.... SUPPLIES DEFAULT PARAMETER VALUES.
C DSGRD2... COMPUTES FINITE-DIFFERENCE GRADIENT APPROXIMATION.
C DSUMIT.... REVERSE-COMMUNICATION ROUTINE THAT DOES DSUMSL ALGORITHM.
C DVSCPY... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR.
C
INTEGER ALPHA, G1, I, IV1, J, K, W
DOUBLE PRECISION ZERO
C
C *** SUBSCRIPTS FOR IV ***
C
INTEGER ETA0, F, G, LMAT, NEXTV, NFGCAL, NGCALL,
1 NITER, SGIRC, TOOBIG, VNEED
C
C/6
DATA ETA0/42/, F/10/, G/28/, LMAT/42/, NEXTV/47/,
1 NFGCAL/7/, NGCALL/30/, NITER/31/, SGIRC/57/,
2 TOOBIG/2/, VNEED/4/
C/7
C PARAMETER (DTYPE=16, ETA0=42, F=10, G=28, LMAT=42, NEXTV=47,
C 1 NFCALL=6, NFGCAL=7, NGCALL=30, NITER=31, SGIRC=57,
C 2 TOOBIG=2, VNEED=4)
C/
C/6
DATA ZERO/0.D+0/
C/7
C PARAMETER (ONE=1.D+0, ZERO=0.D+0)
C/
C
C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++
C
IV1 = IV(1)
IF (IV1 .EQ. 1) GO TO 10
IF (IV1 .EQ. 2) GO TO 50
IF (IV(1) .EQ. 0) CALL DDEFLT(2, IV, LIV, LV, V)
IV(VNEED) = IV(VNEED) + 2*N + 6
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, FX, V(G1), IV, LIV, LV, N, V, X)
IF (IV(1) - 2) 999, 30, 70
C
C *** COMPUTE GRADIENT ***
C
30 IF (IV(NITER) .EQ. 0) CALL DVSCPY(N, V(G1), ZERO)
J = IV(LMAT)
K = G1 - N
DO 40 I = 1, N
V(K) = DDOT(I, V(J),1,V(J),1)
K = K + 1
J = J + I
40 CONTINUE
C *** UNDO INCREMENT OF IV(NGCALL) DONE BY DSUMIT ***
IV(NGCALL) = IV(NGCALL) - 1
C *** STORE RETURN CODE FROM DSGRD2 IN IV(SGIRC) ***
IV(SGIRC) = 0
C *** X MAY HAVE BEEN RESTORED, SO COPY BACK FX... ***
FX = V(F)
GO TO 60
C
C *** GRADIENT LOOP ***
C
50 IF (IV(TOOBIG) .EQ. 0) GO TO 60
IV(NFGCAL) = 0
GO TO 10
C
60 G1 = IV(G)
ALPHA = G1 - N
W = ALPHA - 6
CALL DSGRD2(V(ALPHA), D, V(ETA0), FX, V(G1), IV(SGIRC), N, V(W),X)
IF (IV(SGIRC) .EQ. 0) GO TO 10
IV(NGCALL) = IV(NGCALL) + 1
GO TO 999
C
70 IF (IV(1) .NE. 14) GO TO 999
C
C *** STORAGE ALLOCATION ***
C
IV(G) = IV(NEXTV) + N + 6
IV(NEXTV) = IV(G) + N
IF (IV1 .NE. 13) GO TO 10
C
999 RETURN
C *** LAST CARD OF DSNOIT FOLLOWS ***
END
SUBROUTINE DSGRD2 (ALPHA, D, ETA0, FX, G, IRC, N, W, X)
save
C
C *** COMPUTE FINITE DIFFERENCE GRADIENT BY STWEART*S SCHEME ***
C
C *** PARAMETERS ***
C
INTEGER IRC, N
DOUBLE PRECISION ALPHA(N), D(N), ETA0, FX, G(N), W(6), X(N)
C
C.......................................................................
C
C *** PURPOSE ***
C
C THIS SUBROUTINE USES AN EMBELLISHED FORM OF THE FINITE-DIFFER-
C ENCE SCHEME PROPOSED BY STEWART (REF. 1) TO APPROXIMATE THE
C GRADIENT OF THE FUNCTION F(X), WHOSE VALUES ARE SUPPLIED BY
C REVERSE COMMUNICATION.
C
C *** PARAMETER DESCRIPTION ***
C
C ALPHA IN (APPROXIMATE) DIAGONAL ELEMENTS OF THE HESSIAN OF F(X).
C D IN SCALE VECTOR SUCH THAT D(I)*X(I), I = 1,...,N, ARE IN
C COMPARABLE UNITS.
C ETA0 IN ESTIMATED BOUND ON RELATIVE ERROR IN THE FUNCTION VALUE...
C (TRUE VALUE) = (COMPUTED VALUE)*(1+E), WHERE
C ABS(E) .LE. ETA0.
C FX I/O ON INPUT, FX MUST BE THE COMPUTED VALUE OF F(X). ON
C OUTPUT WITH IRC = 0, FX HAS BEEN RESTORED TO ITS ORIGINAL
C VALUE, THE ONE IT HAD WHEN DSGRD2 WAS LAST CALLED WITH
C IRC = 0.
C G I/O ON INPUT WITH IRC = 0, G SHOULD CONTAIN AN APPROXIMATION
C TO THE GRADIENT OF F NEAR X, E.G., THE GRADIENT AT THE
C PREVIOUS ITERATE. WHEN DSGRD2 RETURNS WITH IRC = 0, G IS
C THE DESIRED FINITE-DIFFERENCE APPROXIMATION TO THE
C GRADIENT AT X.
C IRC I/O INPUT/RETURN CODE... BEFORE THE VERY FIRST CALL ON DSGRD2,
C THE CALLER MUST SET IRC TO 0. WHENEVER DSGRD2 RETURNS A
C NONZERO VALUE FOR IRC, IT HAS PERTURBED SOME COMPONENT OF
C X... THE CALLER SHOULD EVALUATE F(X) AND CALL DSGRD2
C AGAIN WITH FX = F(X).
C N IN THE NUMBER OF VARIABLES (COMPONENTS OF X) ON WHICH F
C DEPENDS.
C X I/O ON INPUT WITH IRC = 0, X IS THE POINT AT WHICH THE
C GRADIENT OF F IS DESIRED. ON OUTPUT WITH IRC NONZERO, X
C IS THE POINT AT WHICH F SHOULD BE EVALUATED. ON OUTPUT
C WITH IRC = 0, X HAS BEEN RESTORED TO ITS ORIGINAL VALUE
C (THE ONE IT HAD WHEN DSGRD2 WAS LAST CALLED WITH IRC = 0)
C AND G CONTAINS THE DESIRED GRADIENT APPROXIMATION.
C W I/O WORK VECTOR OF LENGTH 6 IN WHICH DSGRD2 SAVES CERTAIN
C QUANTITIES WHILE THE CALLER IS EVALUATING F(X) AT A
C PERTURBED X.
C
C *** APPLICATION AND USAGE RESTRICTIONS ***
C
C THIS ROUTINE IS INTENDED FOR USE WITH QUASI-NEWTON ROUTINES
C FOR UNCONSTRAINED MINIMIZATION (IN WHICH CASE ALPHA COMES FROM
C THE DIAGONAL OF THE QUASI-NEWTON HESSIAN APPROXIMATION).
C
C *** ALGORITHM NOTES ***
C
C THIS CODE DEPARTS FROM THE SCHEME PROPOSED BY STEWART (REF. 1)
C IN ITS GUARDING AGAINST OVERLY LARGE OR SMALL STEP SIZES AND ITS
C HANDLING OF SPECIAL CASES (SUCH AS ZERO COMPONENTS OF ALPHA OR G).
C
C *** REFERENCES ***
C
C 1. STEWART, G.W. (1967), A MODIFICATION OF DAVIDON*S MINIMIZATION
C METHOD TO ACCEPT DIFFERENCE APPROXIMATIONS OF DERIVATIVES,
C J. ASSOC. COMPUT. MACH. 14, PP. 72-83.
C
C *** HISTORY ***
C
C DESIGNED AND CODED BY DAVID M. GAY (SUMMER 1977/SUMMER 1980).
C
C *** GENERAL ***
C
C THIS ROUTINE WAS PREPARED IN CONNECTION WITH WORK SUPPORTED BY
C THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS MCS76-00324 AND
C MCS-7906671.
C
C.......................................................................
C
C ***** EXTERNAL FUNCTION *****
C
DOUBLE PRECISION D1MACH
C
C ***** INTRINSIC FUNCTIONS *****
C/+
INTEGER IABS
DOUBLE PRECISION DABS, DMAX1, DSQRT
C/
C ***** LOCAL VARIABLES *****
C
INTEGER FH, FX0, HSAVE, I, XISAVE
DOUBLE PRECISION AAI, AFX, AFXETA, AGI, ALPHAI, AXI, AXIBAR,
1 DISCON, ETA, GI, H, HMIN
DOUBLE PRECISION C2000, FOUR, HMAX0, HMIN0, H0, MACHEP, ONE, P002,
1 THREE, TWO, ZERO
C
C/6
DATA C2000/2.0D+3/, FOUR/4.0D+0/, HMAX0/0.02D+0/, HMIN0/5.0D+1/,
1 ONE/1.0D+0/, P002/0.002D+0/, THREE/3.0D+0/,
2 TWO/2.0D+0/, ZERO/0.0D+0/
C/7
C PARAMETER (C2000=2.0D+3, FOUR=4.0D+0, HMAX0=0.02D+0, HMIN0=5.0D+1,
C 1 ONE=1.0D+0, P002=0.002D+0, THREE=3.0D+0,
C 2 TWO=2.0D+0, ZERO=0.0D+0)
C/
C/6
DATA FH/3/, FX0/4/, HSAVE/5/, XISAVE/6/
C/7
C PARAMETER (FH=3, FX0=4, HSAVE=5, XISAVE=6)
C/
C
C--------------------------------- BODY ------------------------------
C
IF (IRC) 140, 100, 210
C
C *** FRESH START -- GET MACHINE-DEPENDENT CONSTANTS ***
C
C STORE MACHEP IN W(1) AND H0 IN W(2), WHERE MACHEP IS THE UNIT
C ROUNDOFF (THE SMALLEST POSITIVE NUMBER SUCH THAT
C 1 + MACHEP .GT. 1 AND 1 - MACHEP .LT. 1), AND H0 IS THE
C SQUARE-ROOT OF MACHEP.
C
100 W(1) = D1MACH(4)
W(2) = DSQRT(W(1))
C
W(FX0) = FX
C
C *** INCREMENT I AND START COMPUTING G(I) ***
C
110 I = IABS(IRC) + 1
IF (I .GT. N) GO TO 300
IRC = I
AFX = DABS(W(FX0))
MACHEP = W(1)
H0 = W(2)
HMIN = HMIN0 * MACHEP
W(XISAVE) = X(I)
AXI = DABS(X(I))
AXIBAR = DMAX1(AXI, ONE/D(I))
GI = G(I)
AGI = DABS(GI)
ETA = DABS(ETA0)
IF (AFX .GT. ZERO) ETA = DMAX1(ETA, AGI*AXI*MACHEP/AFX)
ALPHAI = ALPHA(I)
IF (ALPHAI .EQ. ZERO) GO TO 170
IF (GI .EQ. ZERO .OR. FX .EQ. ZERO) GO TO 180
AFXETA = AFX*ETA
AAI = DABS(ALPHAI)
C
C *** COMPUTE H = STEWART*S FORWARD-DIFFERENCE STEP SIZE.
C
IF (GI**2 .LE. AFXETA*AAI) GO TO 120
H = TWO*DSQRT(AFXETA/AAI)
H = H*(ONE - AAI*H/(THREE*AAI*H + FOUR*AGI))
GO TO 130
120 H = TWO*(AFXETA*AGI/(AAI**2))**(ONE/THREE)
H = H*(ONE - TWO*AGI/(THREE*AAI*H + FOUR*AGI))
C
C *** ENSURE THAT H IS NOT INSIGNIFICANTLY SMALL ***
C
130 H = DMAX1(H, HMIN*AXIBAR)
C
C *** USE FORWARD DIFFERENCE IF BOUND ON TRUNCATION ERROR IS AT
C *** MOST 10**-3.
C
IF (AAI*H .LE. P002*AGI) GO TO 160
C
C *** COMPUTE H = STEWART*S STEP FOR CENTRAL DIFFERENCE.
C
DISCON = C2000*AFXETA
H = DISCON/(AGI + DSQRT(GI**2 + AAI*DISCON))
C
C *** ENSURE THAT H IS NEITHER TOO SMALL NOR TOO BIG ***
C
H = DMAX1(H, HMIN*AXIBAR)
IF (H .GE. HMAX0*AXIBAR) H = AXIBAR * H0**(TWO/THREE)
C
C *** COMPUTE CENTRAL DIFFERENCE ***
C
IRC = -I
GO TO 200
C
140 H = -W(HSAVE)
I = IABS(IRC)
IF (H .GT. ZERO) GO TO 150
W(FH) = FX
GO TO 200
C
150 G(I) = (W(FH) - FX) / (TWO * H)
X(I) = W(XISAVE)
GO TO 110
C
C *** COMPUTE FORWARD DIFFERENCES IN VARIOUS CASES ***
C
160 IF (H .GE. HMAX0*AXIBAR) H = H0 * AXIBAR
IF (ALPHAI*GI .LT. ZERO) H = -H
GO TO 200
170 H = AXIBAR
GO TO 200
180 H = H0 * AXIBAR
C
200 X(I) = W(XISAVE) + H
W(HSAVE) = H
GO TO 999
C
C *** COMPUTE ACTUAL FORWARD DIFFERENCE ***
C
210 G(IRC) = (FX - W(FX0)) / W(HSAVE)
X(IRC) = W(XISAVE)
GO TO 110
C
C *** RESTORE FX AND INDICATE THAT G HAS BEEN COMPUTED ***
C
300 FX = W(FX0)
IRC = 0
C
999 RETURN
C *** LAST CARD OF DSGRD2 FOLLOWS ***
END