https://github.com/cran/tseries
Tip revision: d37971a68eb391ea02b41bdc4740dfd1c5bd2acf authored by Kurt Hornik on 13 May 2013, 00:00:00 UTC
version 0.10-32
version 0.10-32
Tip revision: d37971a
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(*)
DOUBLE PRECISION D(N), X(N), V(LV), URPARM(*)
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
C DATA NEXTV/47/, NFCALL/6/, NFGCAL/7/, G/28/, TOOBIG/2/, VNEED/4/
C/7
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)
c IF (IV(1) - 2) 30, 40, 50
IF (IV(1) .EQ. 2) GO TO 40
IF (IV(1) .GT. 2) GO TO 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
C DATA 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/7
PARAMETER (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/
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
c standard output unit: unused
IV(PRUNIT) = 6
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
C DATA 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/, X0/43/
C/7
PARAMETER (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,
4 X0=43)
C/
C
C *** V SUBSCRIPT VALUES ***
C
C/6
C DATA 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/7
PARAMETER (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/
C
C/6
C DATA NEGONE/-1.D+0/, ONE/1.D+0/, ZERO/0.D+0/
C/7
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
C DATA ONE/1.D+0/, THREE/3.D+0/
C/7
PARAMETER (ONE=1.D+0, THREE=3.D+0)
C/
C
C *** V SUBSCRIPT VALUES ***
C
C/6
C DATA 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/7
PARAMETER (AFCTOL=31, BIAS=43, COSMIN=47, DECFAC=22, DELTA0=44,
1 DFAC=41, DINIT=38, DLTFDC=42, DLTFDJ=43, DTINIT=39,
2 D0INIT=40, EPSLON=19, ETA0=42, FUZZ=45, HUBERC=48,
3 INCFAC=23, LMAX0=35, LMAXS=36, PHMNFC=20, PHMXFC=21,
4 RDFCMN=24, RDFCMX=25, RFCTOL=32, RLIMIT=46, RSPTOL=49,
5 SCTOL=37, SIGMIN=50, TUNER1=26, TUNER2=27, TUNER3=28,
6 TUNER4=29, TUNER5=30, XCTOL=33, XFTOL=34)
C/
C
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, S(N), W(N), Y(N), Z(N)
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
C DATA EPS/0.1D+0/, ONE/1.D+0/
C/7
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
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
C DATA HALF/0.5D+0/, ONE/1.D+0/, TWO/2.D+0/, ZERO/0.D+0/
C/7
PARAMETER (HALF=0.5D+0, ONE=1.D+0, TWO=2.D+0, ZERO=0.D+0)
C/
C
C/6
C DATA 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/7
PARAMETER (IRC=29, MLSTGD=32, MODEL=5, NFCALL=6, NFGCAL=7,
1 RADINC=8, RESTOR=9, STAGE=10, STGLIM=11, SWITCH=12,
2 TOOBIG=2, XIRC=13)
C/
C/6
C DATA 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/7
PARAMETER (AFCTOL=31, DECFAC=22, DSTNRM=2, DST0=3, DSTSAV=18,
1 F=10, FDIF=11, FLSTGD=12, F0=13, GTSLST=14, GTSTEP=4,
2 INCFAC=23, LMAXS=36, NREDUC=6, PLSTGD=15, PREDUC=7,
3 RADFAC=16, RDFCMN=24, RDFCMX=25, RELDX=17, RFCTOL=32,
4 SCTOL=37, STPPAR=5, TUNER1=26, TUNER2=27, TUNER3=28,
5 XCTOL=33, XFTOL=34)
C/
C
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 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
C DATA HALF/0.5D+0/, ONE/1.D+0/, TWO/2.D+0/, ZERO/0.D+0/
C/7
PARAMETER (HALF=0.5D+0, ONE=1.D+0, TWO=2.D+0, ZERO=0.D+0)
C/
C
C/6
C DATA 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/7
PARAMETER (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/
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
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
C REAL MODEL1(6), MODEL2(6)
C/7
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
C DATA 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/7
PARAMETER (ALGSAV=51, NEEDHD=36, NFCALL=6, NFCOV=52, NGCALL=30,
1 NGCOV=53, NITER=31, OUTLEV=19, PRNTIT=39, PRUNIT=21,
2 SOLPRT=22, STATPR=23, SUSED=64, X0PRT=24)
C/
C
C *** V SUBSCRIPT VALUES ***
C
C/6
C DATA DSTNRM/2/, F/10/, F0/13/, FDIF/11/, NREDUC/6/, PREDUC/7/,
C 1 RELDX/17/, STPPAR/5/
C/7
PARAMETER (DSTNRM=2, F=10, F0=13, FDIF=11, NREDUC=6, PREDUC=7,
1 RELDX=17, STPPAR=5)
C/
C
C/6
C DATA ZERO/0.D+0/
C/7
PARAMETER (ZERO=0.D+0)
C/
C/6
C DATA MODEL1(1)/4H /, MODEL1(2)/4H /, MODEL1(3)/4H /,
C 1 MODEL1(4)/4H /, MODEL1(5)/4H G /, MODEL1(6)/4H S /,
C 2 MODEL2(1)/4H G /, MODEL2(2)/4H S /, MODEL2(3)/4HG-S /,
C 3 MODEL2(4)/4HS-G /, MODEL2(5)/4H-S-G/, MODEL2(6)/4H-G-S/
C/7
DATA MODEL1/' ',' ',' ',' ',' G ',' S '/,
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) call h30()
IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 2) call h40()
IV(NEEDHD) = 0
IF (ALG .EQ. 2) GO TO 50
M = IV(SUSED)
call h100(IV(NITER), NF, V(F), RELDF, PRELDF, V(RELDX),
1 MODEL1(M), MODEL2(M), V(STPPAR))
GO TO 120
C
50 call h110(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) call h70()
IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 2) call h80()
IV(NEEDHD) = 0
NRELDF = ZERO
IF (OLDF .GT. ZERO) NRELDF = V(NREDUC) / OLDF
IF (ALG .EQ. 2) GO TO 90
M = IV(SUSED)
call h100(IV(NITER), NF, V(F), RELDF, PRELDF, V(RELDX),
1 MODEL1(M), MODEL2(M), V(STPPAR), V(DSTNRM), NRELDF)
GO TO 120
C
90 call h110(IV(NITER), NF, V(F), RELDF, PRELDF,
1 V(RELDX), V(STPPAR), V(DSTNRM), NRELDF)
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 call cnlprt(' ***** X-CONVERGENCE *****', 26)
GO TO 430
C
150 call cnlprt(' ***** RELATIVE FUNCTION CONVERGENCE *****', 42)
GO TO 430
C
170 call cnlprt
1(' ***** X- AND RELATIVE FUNCTION CONVERGENCE *****', 49)
GO TO 430
C
190 call cnlprt(' ***** ABSOLUTE FUNCTION CONVERGENCE *****', 42)
GO TO 430
C
210 call cnlprt(' ***** SINGULAR CONVERGENCE *****', 33)
GO TO 430
C
230 call cnlprt(' ***** FALSE CONVERGENCE *****', 30)
GO TO 430
C
250 call cnlprt(' ***** FUNCTION EVALUATION LIMIT *****', 38)
GO TO 430
C
270 call cnlprt(' ***** ITERATION LIMIT *****', 28)
GO TO 430
C
290 call cnlprt(' ***** STOPX *****', 18)
GO TO 430
C
310 call cnlprt(' ***** INITIAL F(X) CANNOT BE COMPUTED *****', 44)
GO TO 390
C
330 call cnlprt(' ***** BAD PARAMETERS TO ASSESS *****', 37)
GO TO 999
C
350 call cnlprt(' ***** GRADIENT COULD NOT BE COMPUTED *****', 43)
IF (IV(NITER) .GT. 0) GO TO 480
GO TO 390
C
370 call h380(IV(1))
GO TO 999
C
C *** INITIAL CALL ON DITSUM ***
C
390 call h400(P, X, D)
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) call h30()
IF (OL .LT. 0 .AND. ALG .EQ. 2) call h40()
IF (OL .GT. 0 .AND. ALG .EQ. 1) call h70()
IF (OL .GT. 0 .AND. ALG .EQ. 2) call h80()
IF (ALG .EQ. 1) call h410(V(F))
IF (ALG .EQ. 2) call h420(V(F))
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)
call h450(V(F), V(RELDX), NF, NG, PRELDF, NRELDF)
C
IF (IV(NFCOV) .GT. 0) call h460(IV(NFCOV))
IF (IV(NGCOV) .GT. 0) call h470(IV(NGCOV))
C
480 IF (IV(SOLPRT) .EQ. 0) GO TO 999
IV(NEEDHD) = 1
call cnlprt(' I FINAL X(I) D(I) G(I)',
1 48)
call h500(P, X, D, G)
GO TO 999
C
520 call cnlprt(' INCONSISTENT DIMENSIONS', 24)
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, Y(N)
DIMENSION L(N*(N+1)/2)
INTEGER I, II, IJ, IM1, I0, J, NP1
DOUBLE PRECISION XI, ZERO
C/6
C DATA ZERO/0.D+0/
C/7
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, Y(N)
DIMENSION L(N*(N+1)/2)
DOUBLE PRECISION DDOT
INTEGER I, J, K
DOUBLE PRECISION T, ZERO
C/6
C DATA ZERO/0.D+0/
C/7
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, Y(N)
DIMENSION L(N*(N+1)/2)
INTEGER I, IJ, I0, J
DOUBLE PRECISION YI, ZERO
C/6
C DATA ZERO/0.D+0/
C/7
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, LAMBDA(N), LPLUS,
1 W(N), Z(N)
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
C DATA ONE/1.D+0/, ZERO/0.D+0/
C/7
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, Y(N)
DIMENSION L(N*(N+1)/2)
INTEGER I, II, IJ, I0, J, NP1
DOUBLE PRECISION T, ZERO
C/6
C DATA ZERO/0.D+0/
C/7
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
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
C INTEGER VARNM(2), SH(2)
C REAL CNGD(3), DFLT(3), VN(2,34), WHICH(3)
C/7
CHARACTER*1 VARNM(2), SH(2)
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
C DATA 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/7
PARAMETER (ALGSAV=51, DINIT=38, DTYPE=16, DTYPE0=54, EPSLON=19,
1 INITS=25, IVNEED=3, LASTIV=44, LASTV=45, LMAT=42,
2 NEXTIV=46, NEXTV=47, NVDFLT=50, OLDN=38, PARPRT=20,
3 PARSAV=49, PERM=58, PRUNIT=21, VNEED=4)
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
C DATA VN(1,1),VN(2,1)/4HEPSL,4HON../
C DATA VN(1,2),VN(2,2)/4HPHMN,4HFC../
C DATA VN(1,3),VN(2,3)/4HPHMX,4HFC../
C DATA VN(1,4),VN(2,4)/4HDECF,4HAC../
C DATA VN(1,5),VN(2,5)/4HINCF,4HAC../
C DATA VN(1,6),VN(2,6)/4HRDFC,4HMN../
C DATA VN(1,7),VN(2,7)/4HRDFC,4HMX../
C DATA VN(1,8),VN(2,8)/4HTUNE,4HR1../
C DATA VN(1,9),VN(2,9)/4HTUNE,4HR2../
C DATA VN(1,10),VN(2,10)/4HTUNE,4HR3../
C DATA VN(1,11),VN(2,11)/4HTUNE,4HR4../
C DATA VN(1,12),VN(2,12)/4HTUNE,4HR5../
C DATA VN(1,13),VN(2,13)/4HAFCT,4HOL../
C DATA VN(1,14),VN(2,14)/4HRFCT,4HOL../
C DATA VN(1,15),VN(2,15)/4HXCTO,4HL.../
C DATA VN(1,16),VN(2,16)/4HXFTO,4HL.../
C DATA VN(1,17),VN(2,17)/4HLMAX,4H0.../
C DATA VN(1,18),VN(2,18)/4HLMAX,4HS.../
C DATA VN(1,19),VN(2,19)/4HSCTO,4HL.../
C DATA VN(1,20),VN(2,20)/4HDINI,4HT.../
C DATA VN(1,21),VN(2,21)/4HDTIN,4HIT../
C DATA VN(1,22),VN(2,22)/4HD0IN,4HIT../
C DATA VN(1,23),VN(2,23)/4HDFAC,4H..../
C DATA VN(1,24),VN(2,24)/4HDLTF,4HDC../
C DATA VN(1,25),VN(2,25)/4HDLTF,4HDJ../
C DATA VN(1,26),VN(2,26)/4HDELT,4HA0../
C DATA VN(1,27),VN(2,27)/4HFUZZ,4H..../
C DATA VN(1,28),VN(2,28)/4HRLIM,4HIT../
C DATA VN(1,29),VN(2,29)/4HCOSM,4HIN../
C DATA VN(1,30),VN(2,30)/4HHUBE,4HRC../
C DATA VN(1,31),VN(2,31)/4HRSPT,4HOL../
C DATA VN(1,32),VN(2,32)/4HSIGM,4HIN../
C DATA VN(1,33),VN(2,33)/4HETA0,4H..../
C DATA VN(1,34),VN(2,34)/4HBIAS,4H..../
C/7
DATA VN(1,1),VN(2,1)/'EPSL','ON..'/
DATA VN(1,2),VN(2,2)/'PHMN','FC..'/
DATA VN(1,3),VN(2,3)/'PHMX','FC..'/
DATA VN(1,4),VN(2,4)/'DECF','AC..'/
DATA VN(1,5),VN(2,5)/'INCF','AC..'/
DATA VN(1,6),VN(2,6)/'RDFC','MN..'/
DATA VN(1,7),VN(2,7)/'RDFC','MX..'/
DATA VN(1,8),VN(2,8)/'TUNE','R1..'/
DATA VN(1,9),VN(2,9)/'TUNE','R2..'/
DATA VN(1,10),VN(2,10)/'TUNE','R3..'/
DATA VN(1,11),VN(2,11)/'TUNE','R4..'/
DATA VN(1,12),VN(2,12)/'TUNE','R5..'/
DATA VN(1,13),VN(2,13)/'AFCT','OL..'/
DATA VN(1,14),VN(2,14)/'RFCT','OL..'/
DATA VN(1,15),VN(2,15)/'XCTO','L...'/
DATA VN(1,16),VN(2,16)/'XFTO','L...'/
DATA VN(1,17),VN(2,17)/'LMAX','0...'/
DATA VN(1,18),VN(2,18)/'LMAX','S...'/
DATA VN(1,19),VN(2,19)/'SCTO','L...'/
DATA VN(1,20),VN(2,20)/'DINI','T...'/
DATA VN(1,21),VN(2,21)/'DTIN','IT..'/
DATA VN(1,22),VN(2,22)/'D0IN','IT..'/
DATA VN(1,23),VN(2,23)/'DFAC','....'/
DATA VN(1,24),VN(2,24)/'DLTF','DC..'/
DATA VN(1,25),VN(2,25)/'DLTF','DJ..'/
DATA VN(1,26),VN(2,26)/'DELT','A0..'/
DATA VN(1,27),VN(2,27)/'FUZZ','....'/
DATA VN(1,28),VN(2,28)/'RLIM','IT..'/
DATA VN(1,29),VN(2,29)/'COSM','IN..'/
DATA VN(1,30),VN(2,30)/'HUBE','RC..'/
DATA VN(1,31),VN(2,31)/'RSPT','OL..'/
DATA VN(1,32),VN(2,32)/'SIGM','IN..'/
DATA VN(1,33),VN(2,33)/'ETA0','....'/
DATA VN(1,34),VN(2,34)/'BIAS','....'/
C/
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
C DATA VARNM(1)/1HP/, VARNM(2)/1HN/, SH(1)/1HS/, SH(2)/1HH/
C DATA CNGD(1),CNGD(2),CNGD(3)/4H---C,4HHANG,4HED V/,
C 1 DFLT(1),DFLT(2),DFLT(3)/4HNOND,4HEFAU,4HLT V/
C/7
DATA VARNM(1)/'P'/, VARNM(2)/'N'/, SH(1)/'S'/, SH(2)/'H'/
DATA CNGD(1),CNGD(2),CNGD(3)/'---C','HANG','ED V'/,
1 DFLT(1),DFLT(2),DFLT(3)/'NOND','EFAU','LT V'/
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
c IF (PU .NE. 0) WRITE(PU,10) ALG, IV(ALGSAV)
c 10 FORMAT(/' THE FIRST PARAMETER TO DDEFLT SHOULD BE',I3,
c 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
c WRITE(PU,30) VARNM(ALG), N
c 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
c WRITE(PU,60) VARNM(ALG), IV(OLDN), N
c 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
c IF (PU .NE. 0) WRITE(PU,80) IV1
c 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
c IF (PU .NE. 0) WRITE(PU,120) VN(1,I), VN(2,I), K, VK,
c 1 VM(I), VX(I)
c 120 FORMAT(/6H /// ,2A4,5H.. V(,I2,3H) =,D11.3,7H SHOULD,
c 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
c WRITE(PU,150) IV(NVDFLT), NDFALT
c 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
c IF (PU .NE. 0) WRITE(PU,170) I, D(I)
c 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
c WRITE(PU,210) SH(ALG), IV(INITS)
c 210 FORMAT(/22H NONDEFAULT VALUES..../5H INIT,A1,14H..... IV(25) =,
c 1 I3)
220 IF (IV(DTYPE) .EQ. IV(DTYPE0)) GO TO 240
c IF (M .EQ. 0) WRITE(PU,250) WHICH
M = 1
c WRITE(PU,230) IV(DTYPE)
c 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
c IF (M .EQ. 0) WRITE(PU,250) WHICH
c 250 FORMAT(/1H ,3A4,9HALUES..../)
M = 1
c WRITE(PU,260) VN(1,I), VN(2,I), K, V(K)
c 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
c WRITE(PU,300) LIV, MIV2
c 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
c WRITE(PU,320) LV, IV(LASTV)
c 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
C DATA ZERO/0.D+0/
C/7
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 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(*)
DOUBLE PRECISION D(N), X(N), V(LV), URPARM(*)
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
C DATA NFCALL/6/, TOOBIG/2/
C/7
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
C DATA ETA0/42/, F/10/, G/28/, LMAT/42/, NEXTV/47/,
C 1 NFGCAL/7/, NGCALL/30/, NITER/31/, SGIRC/57/,
C 2 TOOBIG/2/, VNEED/4/
C/7
PARAMETER (DTYPE=16, ETA0=42, F=10, G=28, LMAT=42, NEXTV=47,
1 NFCALL=6, NFGCAL=7, NGCALL=30, NITER=31, SGIRC=57,
2 TOOBIG=2, VNEED=4)
C/
C/6
C DATA ZERO/0.D+0/
C/7
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)
c IF (IV(1) - 2) 999, 30, 70
IF (IV(1) .LT. 2) GO TO 999
IF (IV(1) .GT. 2) GO TO 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
C DATA 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/7
PARAMETER (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/
C/6
C DATA FH/3/, FX0/4/, HSAVE/5/, XISAVE/6/
C/7
PARAMETER (FH=3, FX0=4, HSAVE=5, XISAVE=6)
C/
C
C--------------------------------- BODY ------------------------------
C
c IF (IRC) 140, 100, 210
IF (IRC .LT. 0) GO TO 140
IF (IRC .GT. 0) GO TO 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
