1C 2C Added to each function and subroutine a "save" and replaced 3C "stop" with a call to the external function "error". "error" 4C has to be provided by the caller, A. Trapletti, 14.10.1999 5C 6C *** These routines are from the NIST Core Math LIBrary CML *** 7C 8 SUBROUTINE DSUMSL(N, D, X, CALCF, CALCG, IV, LIV, LV, V, 9 1 UIPARM, URPARM, UFPARM) 10 save 11C 12C *** MINIMIZE GENERAL UNCONSTRAINED OBJECTIVE FUNCTION USING *** 13C *** ANALYTIC GRADIENT AND HESSIAN APPROX. FROM SECANT UPDATE *** 14C 15 INTEGER N, LIV, LV 16 INTEGER IV(LIV), UIPARM(*) 17 DOUBLE PRECISION D(N), X(N), V(LV), URPARM(*) 18C DIMENSION V(71 + N*(N+15)/2), UIPARM(*), URPARM(*) 19 EXTERNAL CALCF, CALCG, UFPARM 20C 21C *** PURPOSE *** 22C 23C THIS ROUTINE INTERACTS WITH SUBROUTINE DSUMIT IN AN ATTEMPT 24C TO FIND AN N-VECTOR X* THAT MINIMIZES THE (UNCONSTRAINED) 25C OBJECTIVE FUNCTION COMPUTED BY CALCF. (OFTEN THE X* FOUND IS 26C A LOCAL MINIMIZER RATHER THAN A GLOBAL ONE.) 27C 28C-------------------------- PARAMETER USAGE -------------------------- 29C 30C N........ (INPUT) THE NUMBER OF VARIABLES ON WHICH F DEPENDS, I.E., 31C THE NUMBER OF COMPONENTS IN X. 32C D........ (INPUT/OUTPUT) A SCALE VECTOR SUCH THAT D(I)*X(I), 33C I = 1,2,...,N, ARE ALL IN COMPARABLE UNITS. 34C D CAN STRONGLY AFFECT THE BEHAVIOR OF DSUMSL. 35C FINDING THE BEST CHOICE OF D IS GENERALLY A TRIAL- 36C AND-ERROR PROCESS. CHOOSING D SO THAT D(I)*X(I) 37C HAS ABOUT THE SAME VALUE FOR ALL I OFTEN WORKS WELL. 38C THE DEFAULTS PROVIDED BY SUBROUTINE DDEFLT (SEE IV 39C BELOW) REQUIRE THE CALLER TO SUPPLY D. 40C X........ (INPUT/OUTPUT) BEFORE (INITIALLY) CALLING DSUMSL, THE CALL- 41C ER SHOULD SET X TO AN INITIAL GUESS AT X*. WHEN 42C DSUMSL RETURNS, X CONTAINS THE BEST POINT SO FAR 43C FOUND, I.E., THE ONE THAT GIVES THE LEAST VALUE SO 44C FAR SEEN FOR F(X). 45C CALCF.... (INPUT) A SUBROUTINE THAT, GIVEN X, COMPUTES F(X). CALCF 46C MUST BE DECLARED EXTERNAL IN THE CALLING PROGRAM. 47C IT IS INVOKED BY 48C CALL CALCF(N, X, NF, F, UIPARM, URPARM, UFPARM) 49C NF IS THE INVOCATION COUNT FOR CALCF. IT IS INCLUD- 50C ED FOR POSSIBLE USE WITH CALCG. IF X IS OUT OF 51C BOUNDS (E.G., IF IT WOULD CAUSE OVERFLOW IN COMPUT- 52C ING F(X)), THEN CALCF SHOULD SET NF TO 0. THIS WILL 53C CAUSE A SHORTER STEP TO BE ATTEMPTED. THE OTHER 54C PARAMETERS ARE AS DESCRIBED ABOVE AND BELOW. CALCF 55C SHOULD NOT CHANGE N, P, OR X. 56C CALCG.... (INPUT) A SUBROUTINE THAT, GIVEN X, COMPUTES G(X), THE GRA- 57C DIENT OF F AT X. CALCG MUST BE DECLARED EXTERNAL IN 58C THE CALLING PROGRAM. IT IS INVOKED BY 59C CALL CALCG(N, X, NF, G, UIPARM, URPARM, UFAPRM) 60C NF IS THE INVOCATION COUNT FOR CALCF AT THE TIME 61C F(X) WAS EVALUATED. THE X PASSED TO CALCG IS 62C USUALLY THE ONE PASSED TO CALCF ON EITHER ITS MOST 63C RECENT INVOCATION OR THE ONE PRIOR TO IT. IF CALCF 64C SAVES INTERMEDIATE RESULTS FOR USE BY CALCG, THEN IT 65C IS POSSIBLE TO TELL FROM NF WHETHER THEY ARE VALID 66C FOR THE CURRENT X (OR WHICH COPY IS VALID IF TWO 67C COPIES ARE KEPT). IF G CANNOT BE COMPUTED AT X, 68C THEN CALCG SHOULD SET NF TO 0. IN THIS CASE, DSUMSL 69C WILL RETURN WITH IV(1) = 65. THE OTHER PARAMETERS 70C TO CALCG ARE AS DESCRIBED ABOVE AND BELOW. CALCG 71C SHOULD NOT CHANGE N OR X. 72C IV....... (INPUT/OUTPUT) AN INTEGER VALUE ARRAY OF LENGTH LIV (SEE 73C BELOW) THAT HELPS CONTROL THE DSUMSL ALGORITHM AND 74C THAT IS USED TO STORE VARIOUS INTERMEDIATE QUANTI- 75C TIES. OF PARTICULAR INTEREST ARE THE INITIALIZATION/ 76C RETURN CODE IV(1) AND THE ENTRIES IN IV THAT CONTROL 77C PRINTING AND LIMIT THE NUMBER OF ITERATIONS AND FUNC- 78C TION EVALUATIONS. SEE THE SECTION ON IV INPUT 79C VALUES BELOW. 80C V........ (INPUT/OUTPUT) A FLOATING-POINT VALUE ARRAY OF LENGTH LV 81C (SEE BELOW) THAT HELPS CONTROL THE DSUMSL ALGORITHM 82C AND THAT IS USED TO STORE VARIOUS INTERMEDIATE 83C QUANTITIES. OF PARTICULAR INTEREST ARE THE ENTRIES 84C IN V THAT LIMIT THE LENGTH OF THE FIRST STEP 85C ATTEMPTED (LMAX0) AND SPECIFY CONVERGENCE TOLERANCES 86C (AFCTOL, LMAXS, RFCTOL, SCTOL, XCTOL, XFTOL). 87C LIV...... (INPUT) LENGTH OF IV ARRAY. MUST BE AT LEAST 60. IF LIV 88C IS TOO SMALL, THEN DSUMSL RETURNS WITH IV(1) = 15. 89C IF LIV IS AT LEAST LASTIV (= 44), THEN THE MINIMUM 90C ACCEPTABLE VALUE OF LIV IS STORED IN IV(LASTIV) 91C WHEN DSUMSL RETURNS. (THIS IS INTENDED FOR USE 92C WITH EXTENSIONS OF DSUMSL THAT HANDLES CONSTRAINTS.) 93C LV....... (INPUT) LENGTH OF V ARRAY. MUST BE AT LEAST 71+N*(N+15)/2. 94C (AT LEAST 77+N*(N+17)/2 FOR DSMSNO, AT LEAST 95C 78+N*(N+12) FOR DHUMSL). IF LV IS TOO SMALL, THEN 96C DSUMSL RETURNS WITH IV(1) = 16. IF LIV IS AT LEAST 97C LASTV (= 45), THEN THE MINIMUM ACCEPTABLE VALUE OF 98C LV IS STORED IN IV(LASTV) WHEN DSUMSL RETURNS. 99C UIPARM... (INPUT) USER INTEGER PARAMETER ARRAY PASSED WITHOUT CHANGE 100C TO CALCF AND CALCG. 101C URPARM... (INPUT) USER FLOATING-POINT PARAMETER ARRAY PASSED WITHOUT 102C CHANGE TO CALCF AND CALCG. 103C UFPARM... (INPUT) USER EXTERNAL SUBROUTINE OR FUNCTION PASSED WITHOUT 104C CHANGE TO CALCF AND CALCG. 105C 106C *** IV INPUT VALUES (FROM SUBROUTINE DDEFLT) *** 107C 108C IV(1)... ON INPUT, IV(1) SHOULD HAVE A VALUE BETWEEN 0 AND 14...... 109C 0 AND 12 MEAN THIS IS A FRESH START. 0 MEANS THAT 110C DDEFLT(2, IV, LIV, LV, V) 111C IS TO BE CALLED TO PROVIDE ALL DEFAULT VALUES TO IV AND 112C V. 12 (THE VALUE THAT DDEFLT ASSIGNS TO IV(1)) MEANS THE 113C CALLER HAS ALREADY CALLED DDEFLT AND HAS POSSIBLY CHANGED 114C SOME IV AND/OR V ENTRIES TO NON-DEFAULT VALUES. 115C 13 MEANS DDEFLT HAS BEEN CALLED AND THAT DSUMSL (AND 116C DSUMIT) SHOULD ONLY ALLOCATE STORAGE IN IV AND V. 117C 14 MEANS THAT A STORAGE HAS BEEN ALLOCATED (E.G. BY A 118C CALL WITH IV(1) = 13) AND THAT THE ALGORITHM SHOULD BE 119C STARTED. WHEN CALLED WITH IV(1) = 13, DSUMSL RETURNS 120C IV(1) = 14 UNLESS LIV OR LV IS TOO SMALL (OR N IS NOT 121C POSITIVE). DEFAULT = 12. 122C IV(INITH).... IV(25) TELLS WHETHER THE HESSIAN APPROXIMATION H SHOULD 123C BE INITIALIZED. 1 (THE DEFAULT) MEANS DSUMIT SHOULD 124C INITIALIZE H TO THE DIAGONAL MATRIX WHOSE I-TH DIAGONAL 125C ELEMENT IS D(I)**2. 0 MEANS THE CALLER HAS SUPPLIED A 126C CHOLESKY FACTOR L OF THE INITIAL HESSIAN APPROXIMATION 127C H = L*(L**T) IN V, STARTING AT V(IV(LMAT)) = V(IV(42)) 128C (AND STORED COMPACTLY BY ROWS). NOTE THAT IV(LMAT) MAY 129C BE INITIALIZED BY CALLING DSUMSL WITH IV(1) = 13 (SEE 130C THE IV(1) DISCUSSION ABOVE). DEFAULT = 1. 131C IV(MXFCAL)... IV(17) GIVES THE MAXIMUM NUMBER OF FUNCTION EVALUATIONS 132C (CALLS ON CALCF) ALLOWED. IF THIS NUMBER DOES NOT SUF- 133C FICE, THEN DSUMSL RETURNS WITH IV(1) = 9. DEFAULT = 200. 134C IV(MXITER)... IV(18) GIVES THE MAXIMUM NUMBER OF ITERATIONS ALLOWED. 135C IT ALSO INDIRECTLY LIMITS THE NUMBER OF GRADIENT EVALUA- 136C TIONS (CALLS ON CALCG) TO IV(MXITER) + 1. IF IV(MXITER) 137C ITERATIONS DO NOT SUFFICE, THEN DSUMSL RETURNS WITH 138C IV(1) = 10. DEFAULT = 150. 139C IV(OUTLEV)... IV(19) CONTROLS THE NUMBER AND LENGTH OF ITERATION SUM- 140C MARY LINES PRINTED (BY DITSUM). IV(OUTLEV) = 0 MEANS DO 141C NOT PRINT ANY SUMMARY LINES. OTHERWISE, PRINT A SUMMARY 142C LINE AFTER EACH ABS(IV(OUTLEV)) ITERATIONS. IF IV(OUTLEV) 143C IS POSITIVE, THEN SUMMARY LINES OF LENGTH 78 (PLUS CARRI- 144C AGE CONTROL) ARE PRINTED, INCLUDING THE FOLLOWING... THE 145C ITERATION AND FUNCTION EVALUATION COUNTS, F = THE CURRENT 146C FUNCTION VALUE, RELATIVE DIFFERENCE IN FUNCTION VALUES 147C ACHIEVED BY THE LATEST STEP (I.E., RELDF = (F0-V(F))/F01, 148C WHERE F01 IS THE MAXIMUM OF ABS(V(F)) AND ABS(V(F0)) AND 149C V(F0) IS THE FUNCTION VALUE FROM THE PREVIOUS ITERA- 150C TION), THE RELATIVE FUNCTION REDUCTION PREDICTED FOR THE 151C STEP JUST TAKEN (I.E., PRELDF = V(PREDUC) / F01, WHERE 152C V(PREDUC) IS DESCRIBED BELOW), THE SCALED RELATIVE CHANGE 153C IN X (SEE V(RELDX) BELOW), THE STEP PARAMETER FOR THE 154C STEP JUST TAKEN (STPPAR = 0 MEANS A FULL NEWTON STEP, 155C BETWEEN 0 AND 1 MEANS A RELAXED NEWTON STEP, BETWEEN 1 156C AND 2 MEANS A DOUBLE DOGLEG STEP, GREATER THAN 2 MEANS 157C A SCALED DOWN CAUCHY STEP -- SEE SUBROUTINE DBLDOG), THE 158C 2-NORM OF THE SCALE VECTOR D TIMES THE STEP JUST TAKEN 159C (SEE V(DSTNRM) BELOW), AND NPRELDF, I.E., 160C V(NREDUC)/F01, WHERE V(NREDUC) IS DESCRIBED BELOW -- IF 161C NPRELDF IS POSITIVE, THEN IT IS THE RELATIVE FUNCTION 162C REDUCTION PREDICTED FOR A NEWTON STEP (ONE WITH 163C STPPAR = 0). IF NPRELDF IS NEGATIVE, THEN IT IS THE 164C NEGATIVE OF THE RELATIVE FUNCTION REDUCTION PREDICTED 165C FOR A STEP COMPUTED WITH STEP BOUND V(LMAXS) FOR USE IN 166C TESTING FOR SINGULAR CONVERGENCE. 167C IF IV(OUTLEV) IS NEGATIVE, THEN LINES OF LENGTH 50 168C ARE PRINTED, INCLUDING ONLY THE FIRST 6 ITEMS LISTED 169C ABOVE (THROUGH RELDX). 170C DEFAULT = 1. 171C IV(PARPRT)... IV(20) = 1 MEANS PRINT ANY NONDEFAULT V VALUES ON A 172C FRESH START OR ANY CHANGED V VALUES ON A RESTART. 173C IV(PARPRT) = 0 MEANS SKIP THIS PRINTING. DEFAULT = 1. 174C IV(PRUNIT)... IV(21) IS THE OUTPUT UNIT NUMBER ON WHICH ALL PRINTING 175C IS DONE. IV(PRUNIT) = 0 MEANS SUPPRESS ALL PRINTING. 176C DEFAULT = STANDARD OUTPUT UNIT (UNIT 6 ON MOST SYSTEMS). 177C IV(SOLPRT)... IV(22) = 1 MEANS PRINT OUT THE VALUE OF X RETURNED (AS 178C WELL AS THE GRADIENT AND THE SCALE VECTOR D). 179C IV(SOLPRT) = 0 MEANS SKIP THIS PRINTING. DEFAULT = 1. 180C IV(STATPR)... IV(23) = 1 MEANS PRINT SUMMARY STATISTICS UPON RETURN- 181C ING. THESE CONSIST OF THE FUNCTION VALUE, THE SCALED 182C RELATIVE CHANGE IN X CAUSED BY THE MOST RECENT STEP (SEE 183C V(RELDX) BELOW), THE NUMBER OF FUNCTION AND GRADIENT 184C EVALUATIONS (CALLS ON CALCF AND CALCG), AND THE RELATIVE 185C FUNCTION REDUCTIONS PREDICTED FOR THE LAST STEP TAKEN AND 186C FOR A NEWTON STEP (OR PERHAPS A STEP BOUNDED BY V(LMAX0) 187C -- SEE THE DESCRIPTIONS OF PRELDF AND NPRELDF UNDER 188C IV(OUTLEV) ABOVE). 189C IV(STATPR) = 0 MEANS SKIP THIS PRINTING. 190C IV(STATPR) = -1 MEANS SKIP THIS PRINTING AS WELL AS THAT 191C OF THE ONE-LINE TERMINATION REASON MESSAGE. DEFAULT = 1. 192C IV(X0PRT).... IV(24) = 1 MEANS PRINT THE INITIAL X AND SCALE VECTOR D 193C (ON A FRESH START ONLY). IV(X0PRT) = 0 MEANS SKIP THIS 194C PRINTING. DEFAULT = 1. 195C 196C *** (SELECTED) IV OUTPUT VALUES *** 197C 198C IV(1)........ ON OUTPUT, IV(1) IS A RETURN CODE.... 199C 3 = X-CONVERGENCE. THE SCALED RELATIVE DIFFERENCE (SEE 200C V(RELDX)) BETWEEN THE CURRENT PARAMETER VECTOR X AND 201C A LOCALLY OPTIMAL PARAMETER VECTOR IS VERY LIKELY AT 202C MOST V(XCTOL). 203C 4 = RELATIVE FUNCTION CONVERGENCE. THE RELATIVE DIFFER- 204C ENCE BETWEEN THE CURRENT FUNCTION VALUE AND ITS LO- 205C CALLY OPTIMAL VALUE IS VERY LIKELY AT MOST V(RFCTOL). 206C 5 = BOTH X- AND RELATIVE FUNCTION CONVERGENCE (I.E., THE 207C CONDITIONS FOR IV(1) = 3 AND IV(1) = 4 BOTH HOLD). 208C 6 = ABSOLUTE FUNCTION CONVERGENCE. THE CURRENT FUNCTION 209C VALUE IS AT MOST V(AFCTOL) IN ABSOLUTE VALUE. 210C 7 = SINGULAR CONVERGENCE. THE HESSIAN NEAR THE CURRENT 211C ITERATE APPEARS TO BE SINGULAR OR NEARLY SO, AND A 212C STEP OF LENGTH AT MOST V(LMAX0) IS UNLIKELY TO YIELD 213C A RELATIVE FUNCTION DECREASE OF MORE THAN V(SCTOL). 214C 8 = FALSE CONVERGENCE. THE ITERATES APPEAR TO BE CONVERG- 215C ING TO A NONCRITICAL POINT. THIS MAY MEAN THAT THE 216C CONVERGENCE TOLERANCES (V(AFCTOL), V(RFCTOL), 217C V(XCTOL)) ARE TOO SMALL FOR THE ACCURACY TO WHICH 218C THE FUNCTION AND GRADIENT ARE BEING COMPUTED, THAT 219C THERE IS AN ERROR IN COMPUTING THE GRADIENT, OR THAT 220C THE FUNCTION OR GRADIENT IS DISCONTINUOUS NEAR X. 221C 9 = FUNCTION EVALUATION LIMIT REACHED WITHOUT OTHER CON- 222C VERGENCE (SEE IV(MXFCAL)). 223C 10 = ITERATION LIMIT REACHED WITHOUT OTHER CONVERGENCE 224C (SEE IV(MXITER)). 225C 11 = DSTOPX RETURNED .TRUE. (EXTERNAL INTERRUPT). SEE THE 226C USAGE NOTES BELOW. 227C 14 = STORAGE HAS BEEN ALLOCATED (AFTER A CALL WITH 228C IV(1) = 13). 229C 17 = RESTART ATTEMPTED WITH N CHANGED. 230C 18 = D HAS A NEGATIVE COMPONENT AND IV(DTYPE) .LE. 0. 231C 19...43 = V(IV(1)) IS OUT OF RANGE. 232C 63 = F(X) CANNOT BE COMPUTED AT THE INITIAL X. 233C 64 = BAD PARAMETERS PASSED TO ASSESS (WHICH SHOULD NOT 234C OCCUR). 235C 65 = THE GRADIENT COULD NOT BE COMPUTED AT X (SEE CALCG 236C ABOVE). 237C 67 = BAD FIRST PARAMETER TO DDEFLT. 238C 80 = IV(1) WAS OUT OF RANGE. 239C 81 = N IS NOT POSITIVE. 240C IV(G)........ IV(28) IS THE STARTING SUBSCRIPT IN V OF THE CURRENT 241C GRADIENT VECTOR (THE ONE CORRESPONDING TO X). 242C IV(NFCALL)... IV(6) IS THE NUMBER OF CALLS SO FAR MADE ON CALCF (I.E., 243C FUNCTION EVALUATIONS). 244C IV(NGCALL)... IV(30) IS THE NUMBER OF GRADIENT EVALUATIONS (CALLS ON 245C CALCG). 246C IV(NITER).... IV(31) IS THE NUMBER OF ITERATIONS PERFORMED. 247C 248C *** (SELECTED) V INPUT VALUES (FROM SUBROUTINE DDEFLT) *** 249C 250C V(BIAS)..... V(43) IS THE BIAS PARAMETER USED IN SUBROUTINE DBLDOG -- 251C SEE THAT SUBROUTINE FOR DETAILS. DEFAULT = 0.8. 252C V(AFCTOL)... V(31) IS THE ABSOLUTE FUNCTION CONVERGENCE TOLERANCE. 253C IF DSUMSL FINDS A POINT WHERE THE FUNCTION VALUE IS LESS 254C THAN V(AFCTOL) IN ABSOLUTE VALUE, AND IF DSUMSL DOES NOT 255C RETURN WITH IV(1) = 3, 4, OR 5, THEN IT RETURNS WITH 256C IV(1) = 6. DEFAULT = MAX(10**-20, MACHEP**2), WHERE 257C MACHEP IS THE UNIT ROUNDOFF. 258C V(DINIT).... V(38), IF NONNEGATIVE, IS THE VALUE TO WHICH THE SCALE 259C VECTOR D IS INITIALIZED. DEFAULT = -1. 260C V(LMAX0).... V(35) GIVES THE MAXIMUM 2-NORM ALLOWED FOR D TIMES THE 261C VERY FIRST STEP THAT DSUMSL ATTEMPTS. THIS PARAMETER CAN 262C MARKEDLY AFFECT THE PERFORMANCE OF DSUMSL. 263C V(LMAXS).... V(36) IS USED IN TESTING FOR SINGULAR CONVERGENCE -- IF 264C THE FUNCTION REDUCTION PREDICTED FOR A STEP OF LENGTH 265C BOUNDED BY V(LMAXS) IS AT MOST V(SCTOL) * ABS(F0), WHERE 266C F0 IS THE FUNCTION VALUE AT THE START OF THE CURRENT 267C ITERATION, AND IF DSUMSL DOES NOT RETURN WITH IV(1) = 3, 268C 4, 5, OR 6, THEN IT RETURNS WITH IV(1) = 7. DEFAULT = 1. 269C V(RFCTOL)... V(32) IS THE RELATIVE FUNCTION CONVERGENCE TOLERANCE. 270C IF THE CURRENT MODEL PREDICTS A MAXIMUM POSSIBLE FUNCTION 271C REDUCTION (SEE V(NREDUC)) OF AT MOST V(RFCTOL)*ABS(F0) 272C AT THE START OF THE CURRENT ITERATION, WHERE F0 IS THE 273C THEN CURRENT FUNCTION VALUE, AND IF THE LAST STEP ATTEMPT- 274C ED ACHIEVED NO MORE THAN TWICE THE PREDICTED FUNCTION 275C DECREASE, THEN DSUMSL RETURNS WITH IV(1) = 4 (OR 5). 276C DEFAULT = MAX(10**-10, MACHEP**(2/3)), WHERE MACHEP IS 277C THE UNIT ROUNDOFF. 278C V(SCTOL).... V(37) IS THE SINGULAR CONVERGENCE TOLERANCE -- SEE THE 279C DESCRIPTION OF V(LMAXS) ABOVE. 280C V(TUNER1)... V(26) HELPS DECIDE WHEN TO CHECK FOR FALSE CONVERGENCE. 281C THIS IS DONE IF THE ACTUAL FUNCTION DECREASE FROM THE 282C CURRENT STEP IS NO MORE THAN V(TUNER1) TIMES ITS PREDICT- 283C ED VALUE. DEFAULT = 0.1. 284C V(XCTOL).... V(33) IS THE X-CONVERGENCE TOLERANCE. IF A NEWTON STEP 285C (SEE V(NREDUC)) IS TRIED THAT HAS V(RELDX) .LE. V(XCTOL) 286C AND IF THIS STEP YIELDS AT MOST TWICE THE PREDICTED FUNC- 287C TION DECREASE, THEN DSUMSL RETURNS WITH IV(1) = 3 (OR 5). 288C (SEE THE DESCRIPTION OF V(RELDX) BELOW.) 289C DEFAULT = MACHEP**0.5, WHERE MACHEP IS THE UNIT ROUNDOFF. 290C V(XFTOL).... V(34) IS THE FALSE CONVERGENCE TOLERANCE. IF A STEP IS 291C TRIED THAT GIVES NO MORE THAN V(TUNER1) TIMES THE PREDICT- 292C ED FUNCTION DECREASE AND THAT HAS V(RELDX) .LE. V(XFTOL), 293C AND IF DSUMSL DOES NOT RETURN WITH IV(1) = 3, 4, 5, 6, OR 294C 7, THEN IT RETURNS WITH IV(1) = 8. (SEE THE DESCRIPTION 295C OF V(RELDX) BELOW.) DEFAULT = 100*MACHEP, WHERE 296C MACHEP IS THE UNIT ROUNDOFF. 297C V(*)........ DDEFLT SUPPLIES TO V A NUMBER OF TUNING CONSTANTS, WITH 298C WHICH IT SHOULD ORDINARILY BE UNNECESSARY TO TINKER. SEE 299C SECTION 17 OF VERSION 2.2 OF THE NL2SOL USAGE SUMMARY 300C (I.E., THE APPENDIX TO REF. 1) FOR DETAILS ON V(I), 301C I = DECFAC, INCFAC, PHMNFC, PHMXFC, RDFCMN, RDFCMX, 302C TUNER2, TUNER3, TUNER4, TUNER5. 303C 304C *** (SELECTED) V OUTPUT VALUES *** 305C 306C V(DGNORM)... V(1) IS THE 2-NORM OF (DIAG(D)**-1)*G, WHERE G IS THE 307C MOST RECENTLY COMPUTED GRADIENT. 308C V(DSTNRM)... V(2) IS THE 2-NORM OF DIAG(D)*STEP, WHERE STEP IS THE 309C CURRENT STEP. 310C V(F)........ V(10) IS THE CURRENT FUNCTION VALUE. 311C V(F0)....... V(13) IS THE FUNCTION VALUE AT THE START OF THE CURRENT 312C ITERATION. 313C V(NREDUC)... V(6), IF POSITIVE, IS THE MAXIMUM FUNCTION REDUCTION 314C POSSIBLE ACCORDING TO THE CURRENT MODEL, I.E., THE FUNC- 315C TION REDUCTION PREDICTED FOR A NEWTON STEP (I.E., 316C STEP = -H**-1 * G, WHERE G IS THE CURRENT GRADIENT AND 317C H IS THE CURRENT HESSIAN APPROXIMATION). 318C IF V(NREDUC) IS NEGATIVE, THEN IT IS THE NEGATIVE OF 319C THE FUNCTION REDUCTION PREDICTED FOR A STEP COMPUTED WITH 320C A STEP BOUND OF V(LMAXS) FOR USE IN TESTING FOR SINGULAR 321C CONVERGENCE. 322C V(PREDUC)... V(7) IS THE FUNCTION REDUCTION PREDICTED (BY THE CURRENT 323C QUADRATIC MODEL) FOR THE CURRENT STEP. THIS (DIVIDED BY 324C V(F0)) IS USED IN TESTING FOR RELATIVE FUNCTION 325C CONVERGENCE. 326C V(RELDX).... V(17) IS THE SCALED RELATIVE CHANGE IN X CAUSED BY THE 327C CURRENT STEP, COMPUTED AS 328C MAX(ABS(D(I)*(X(I)-X0(I)), 1 .LE. I .LE. P) / 329C MAX(D(I)*(ABS(X(I))+ABS(X0(I))), 1 .LE. I .LE. P), 330C WHERE X = X0 + STEP. 331C 332C------------------------------- NOTES ------------------------------- 333C 334C *** ALGORITHM NOTES *** 335C 336C THIS ROUTINE USES A HESSIAN APPROXIMATION COMPUTED FROM THE 337C BFGS UPDATE (SEE REF 3). ONLY A CHOLESKY FACTOR OF THE HESSIAN 338C APPROXIMATION IS STORED, AND THIS IS UPDATED USING IDEAS FROM 339C REF. 4. STEPS ARE COMPUTED BY THE DOUBLE DOGLEG SCHEME DESCRIBED 340C IN REF. 2. THE STEPS ARE ASSESSED AS IN REF. 1. 341C 342C *** USAGE NOTES *** 343C 344C AFTER A RETURN WITH IV(1) .LE. 11, IT IS POSSIBLE TO RESTART, 345C I.E., TO CHANGE SOME OF THE IV AND V INPUT VALUES DESCRIBED ABOVE 346C AND CONTINUE THE ALGORITHM FROM THE POINT WHERE IT WAS INTERRUPT- 347C ED. IV(1) SHOULD NOT BE CHANGED, NOR SHOULD ANY ENTRIES OF IV 348C AND V OTHER THAN THE INPUT VALUES (THOSE SUPPLIED BY DDEFLT). 349C THOSE WHO DO NOT WISH TO WRITE A CALCG WHICH COMPUTES THE 350C GRADIENT ANALYTICALLY SHOULD CALL DSMSNO RATHER THAN DSUMSL. 351C DSMSNO USES FINITE DIFFERENCES TO COMPUTE AN APPROXIMATE GRADIENT. 352C THOSE WHO WOULD PREFER TO PROVIDE F AND G (THE FUNCTION AND 353C GRADIENT) BY REVERSE COMMUNICATION RATHER THAN BY WRITING SUBROU- 354C TINES CALCF AND CALCG MAY CALL ON DSUMIT DIRECTLY. SEE THE COM- 355C MENTS AT THE BEGINNING OF DSUMIT. 356C THOSE WHO USE DSUMSL INTERACTIVELY MAY WISH TO SUPPLY THEIR 357C OWN DSTOPX FUNCTION, WHICH SHOULD RETURN .TRUE. IF THE BREAK KEY 358C HAS BEEN PRESSED SINCE DSTOPX WAS LAST INVOKED. THIS MAKES IT 359C POSSIBLE TO EXTERNALLY INTERRUPT DSUMSL (WHICH WILL RETURN WITH 360C IV(1) = 11 IF DSTOPX RETURNS .TRUE.). 361C STORAGE FOR G IS ALLOCATED AT THE END OF V. THUS THE CALLER 362C MAY MAKE V LONGER THAN SPECIFIED ABOVE AND MAY ALLOW CALCG TO USE 363C ELEMENTS OF G BEYOND THE FIRST N AS SCRATCH STORAGE. 364C 365C *** PORTABILITY NOTES *** 366C 367C THE DSUMSL DISTRIBUTION TAPE CONTAINS BOTH SINGLE- AND DOUBLE- 368C PRECISION VERSIONS OF THE DSUMSL SOURCE CODE, SO IT SHOULD BE UN- 369C NECESSARY TO CHANGE PRECISIONS. 370C INTRINSIC FUNCTIONS ARE EXPLICITLY DECLARED. ON CERTAIN COM- 371C PUTERS (E.G. UNIVAC), IT MAY BE NECESSARY TO COMMENT OUT THESE 372C DECLARATIONS. SO THAT THIS MAY BE DONE AUTOMATICALLY BY A SIMPLE 373C PROGRAM, SUCH DECLARATIONS ARE PRECEDED BY A COMMENT HAVING C/+ 374C IN COLUMNS 1-3 AND BLANKS IN COLUMNS 4-72 AND ARE FOLLOWED BY 375C A COMMENT HAVING C/ IN COLUMNS 1 AND 2 AND BLANKS IN COLUMNS 3-72. 376C THE DSUMSL SOURCE CODE IS EXPRESSED IN 1966 ANSI STANDARD 377C FORTRAN. IT MAY BE CONVERTED TO FORTRAN 77 BY COMMENTING OUT ALL 378C LINES THAT FALL BETWEEN A LINE HAVING C/6 IN COLUMNS 1-3 AND A 379C LINE HAVING C/7 IN COLUMNS 1-3 AND BY REMOVING (I.E., REPLACING 380C BY A BLANK) THE C IN COLUMN 1 OF THE LINES THAT FOLLOW THE C/7 381C LINE AND PRECEDE A LINE HAVING C/ IN COLUMNS 1-2 AND BLANKS IN 382C COLUMNS 3-72. THESE CHANGES CONVERT SOME DATA STATEMENTS INTO 383C PARAMETER STATEMENTS, CONVERT SOME VARIABLES FROM REAL TO 384C CHARACTER*4, AND MAKE THE DATA STATEMENTS THAT INITIALIZE THESE 385C VARIABLES USE CHARACTER STRINGS DELIMITED BY PRIMES INSTEAD 386C OF HOLLERITH CONSTANTS. (SUCH VARIABLES AND DATA STATEMENTS 387C APPEAR ONLY IN MODULES DITSUM AND DPARCK. PARAMETER STATEMENTS 388C APPEAR NEARLY EVERYWHERE.) 389C 390C *** REFERENCES *** 391C 392C 1. DENNIS, J.E., GAY, D.M., AND WELSCH, R.E. (1981), ALGORITHM 573 -- 393C AN ADAPTIVE NONLINEAR LEAST-SQUARES ALGORITHM, ACM TRANS. 394C MATH. SOFTWARE 7, PP. 369-383. 395C 396C 2. DENNIS, J.E., AND MEI, H.H.W. (1979), TWO NEW UNCONSTRAINED OPTI- 397C MIZATION ALGORITHMS WHICH USE FUNCTION AND GRADIENT 398C VALUES, J. OPTIM. THEORY APPLIC. 28, PP. 453-482. 399C 400C 3. DENNIS, J.E., AND MORE, J.J. (1977), QUASI-NEWTON METHODS, MOTIVA- 401C TION AND THEORY, SIAM REV. 19, PP. 46-89. 402C 403C 4. GOLDFARB, D. (1976), FACTORIZED VARIABLE METRIC METHODS FOR UNCON- 404C STRAINED OPTIMIZATION, MATH. COMPUT. 30, PP. 796-811. 405C 406C *** GENERAL *** 407C 408C CODED BY DAVID M. GAY (WINTER 1980). REVISED SUMMER 1982. 409C THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH 410C SUPPORTED IN PART BY THE NATIONAL SCIENCE FOUNDATION UNDER 411C GRANTS MCS-7600324, DCR75-10143, 76-14311DSS, MCS76-11989, 412C AND MCS-7906671. 413C 414C. 415C---------------------------- DECLARATIONS --------------------------- 416C 417 EXTERNAL DDEFLT, DSUMIT 418C 419C DDEFLT.... SUPPLIES DEFAULT IV AND V INPUT COMPONENTS. 420C DSUMIT... REVERSE-COMMUNICATION ROUTINE THAT CARRIES OUT DSUMSL ALGO- 421C RITHM. 422C 423 INTEGER G1, IV1, NF 424 DOUBLE PRECISION F 425C 426C *** SUBSCRIPTS FOR IV *** 427C 428 INTEGER NEXTV, NFCALL, NFGCAL, G, TOOBIG, VNEED 429C 430C/6 431C DATA NEXTV/47/, NFCALL/6/, NFGCAL/7/, G/28/, TOOBIG/2/, VNEED/4/ 432C/7 433 PARAMETER (NEXTV=47, NFCALL=6, NFGCAL=7, G=28, TOOBIG=2, VNEED=4) 434C/ 435C 436C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ 437C 438 IF (IV(1) .EQ. 0) CALL DDEFLT(2, IV, LIV, LV, V) 439 IV(VNEED) = IV(VNEED) + N 440 IV1 = IV(1) 441 IF (IV1 .EQ. 14) GO TO 10 442 IF (IV1 .GT. 2 .AND. IV1 .LT. 12) GO TO 10 443 G1 = 1 444 IF (IV1 .EQ. 12) IV(1) = 13 445 GO TO 20 446C 447 10 G1 = IV(G) 448C 449 20 CALL DSUMIT(D, F, V(G1), IV, LIV, LV, N, V, X) 450c IF (IV(1) - 2) 30, 40, 50 451 IF (IV(1) .EQ. 2) GO TO 40 452 IF (IV(1) .GT. 2) GO TO 50 453C 454 30 NF = IV(NFCALL) 455 CALL CALCF(N, X, NF, F, UIPARM, URPARM, UFPARM) 456 IF (NF .LE. 0) IV(TOOBIG) = 1 457 GO TO 20 458C 459 40 CALL CALCG(N, X, IV(NFGCAL), V(G1), UIPARM, URPARM, UFPARM) 460 GO TO 20 461C 462 50 IF (IV(1) .NE. 14) GO TO 999 463C 464C *** STORAGE ALLOCATION 465C 466 IV(G) = IV(NEXTV) 467 IV(NEXTV) = IV(G) + N 468 IF (IV1 .NE. 13) GO TO 10 469C 470 999 RETURN 471C *** LAST CARD OF DSUMSL FOLLOWS *** 472 END 473 SUBROUTINE DDEFLT(ALG, IV, LIV, LV, V) 474 save 475C 476C *** SUPPLY ***SOL (VERSION 2.3) DEFAULT VALUES TO IV AND V *** 477C 478C *** ALG = 1 MEANS REGRESSION CONSTANTS. 479C *** ALG = 2 MEANS GENERAL UNCONSTRAINED OPTIMIZATION CONSTANTS. 480C 481 INTEGER LIV, LV 482 INTEGER ALG, IV(LIV) 483 DOUBLE PRECISION V(LV) 484C 485 EXTERNAL DVDFLT 486C DVDFLT.... PROVIDES DEFAULT VALUES TO V. 487C 488 INTEGER MIV, MV 489 INTEGER MINIV(2), MINV(2) 490C 491C *** SUBSCRIPTS FOR IV *** 492C 493 INTEGER ALGSAV, COVPRT, COVREQ, DTYPE, HC, IERR, INITH, INITS, 494 1 IPIVOT, IVNEED, LASTIV, LASTV, LMAT, MXFCAL, MXITER, 495 2 NFCOV, NGCOV, NVDFLT, OUTLEV, PARPRT, PARSAV, PERM, 496 3 PRUNIT, QRTYP, RDREQ, RMAT, SOLPRT, STATPR, VNEED, 497 4 VSAVE, X0PRT 498C 499C *** IV SUBSCRIPT VALUES *** 500C 501C/6 502C DATA ALGSAV/51/, COVPRT/14/, COVREQ/15/, DTYPE/16/, HC/71/, 503C 1 IERR/75/, INITH/25/, INITS/25/, IPIVOT/76/, IVNEED/3/, 504C 2 LASTIV/44/, LASTV/45/, LMAT/42/, MXFCAL/17/, MXITER/18/, 505C 3 NFCOV/52/, NGCOV/53/, NVDFLT/50/, OUTLEV/19/, PARPRT/20/, 506C 4 PARSAV/49/, PERM/58/, PRUNIT/21/, QRTYP/80/, RDREQ/57/, 507C 5 RMAT/78/, SOLPRT/22/, STATPR/23/, VNEED/4/, VSAVE/60/, 508C 6 X0PRT/24/ 509C/7 510 PARAMETER (ALGSAV=51, COVPRT=14, COVREQ=15, DTYPE=16, HC=71, 511 1 IERR=75, INITH=25, INITS=25, IPIVOT=76, IVNEED=3, 512 2 LASTIV=44, LASTV=45, LMAT=42, MXFCAL=17, MXITER=18, 513 3 NFCOV=52, NGCOV=53, NVDFLT=50, OUTLEV=19, PARPRT=20, 514 4 PARSAV=49, PERM=58, PRUNIT=21, QRTYP=80, RDREQ=57, 515 5 RMAT=78, SOLPRT=22, STATPR=23, VNEED=4, VSAVE=60, 516 6 X0PRT=24) 517C/ 518 DATA MINIV(1)/80/, MINIV(2)/59/, MINV(1)/98/, MINV(2)/71/ 519C 520C------------------------------- BODY -------------------------------- 521C 522 IF (ALG .LT. 1 .OR. ALG .GT. 2) GO TO 40 523 MIV = MINIV(ALG) 524 IF (LIV .LT. MIV) GO TO 20 525 MV = MINV(ALG) 526 IF (LV .LT. MV) GO TO 30 527 CALL DVDFLT(ALG, LV, V) 528 IV(1) = 12 529 IV(ALGSAV) = ALG 530 IV(IVNEED) = 0 531 IV(LASTIV) = MIV 532 IV(LASTV) = MV 533 IV(LMAT) = MV + 1 534 IV(MXFCAL) = 200 535 IV(MXITER) = 150 536 IV(OUTLEV) = 1 537 IV(PARPRT) = 1 538 IV(PERM) = MIV + 1 539c standard output unit: unused 540 IV(PRUNIT) = 6 541 IV(SOLPRT) = 1 542 IV(STATPR) = 1 543 IV(VNEED) = 0 544 IV(X0PRT) = 1 545C 546 IF (ALG .GE. 2) GO TO 10 547C 548C *** REGRESSION VALUES 549C 550 IV(COVPRT) = 3 551 IV(COVREQ) = 1 552 IV(DTYPE) = 1 553 IV(HC) = 0 554 IV(IERR) = 0 555 IV(INITS) = 0 556 IV(IPIVOT) = 0 557 IV(NVDFLT) = 32 558 IV(PARSAV) = 67 559 IV(QRTYP) = 1 560 IV(RDREQ) = 3 561 IV(RMAT) = 0 562 IV(VSAVE) = 58 563 GO TO 999 564C 565C *** GENERAL OPTIMIZATION VALUES 566C 567 10 IV(DTYPE) = 0 568 IV(INITH) = 1 569 IV(NFCOV) = 0 570 IV(NGCOV) = 0 571 IV(NVDFLT) = 25 572 IV(PARSAV) = 47 573 GO TO 999 574C 575 20 IV(1) = 15 576 GO TO 999 577C 578 30 IV(1) = 16 579 GO TO 999 580C 581 40 IV(1) = 67 582C 583 999 RETURN 584C *** LAST CARD OF DDEFLT FOLLOWS *** 585 END 586 SUBROUTINE DSUMIT(D, FX, G, IV, LIV, LV, N, V, X) 587 save 588C 589C *** CARRY OUT DSUMSL (UNCONSTRAINED MINIMIZATION) ITERATIONS, USING 590C *** DOUBLE-DOGLEG/BFGS STEPS. 591C 592C *** PARAMETER DECLARATIONS *** 593C 594 INTEGER LIV, LV, N 595 INTEGER IV(LIV) 596 DOUBLE PRECISION D(N), FX, G(N), V(LV), X(N) 597C 598C-------------------------- PARAMETER USAGE -------------------------- 599C 600C D.... SCALE VECTOR. 601C FX... FUNCTION VALUE. 602C G.... GRADIENT VECTOR. 603C IV... INTEGER VALUE ARRAY. 604C LIV.. LENGTH OF IV (AT LEAST 60). 605C LV... LENGTH OF V (AT LEAST 71 + N*(N+13)/2). 606C N.... NUMBER OF VARIABLES (COMPONENTS IN X AND G). 607C V.... FLOATING-POINT VALUE ARRAY. 608C X.... VECTOR OF PARAMETERS TO BE OPTIMIZED. 609C 610C *** DISCUSSION *** 611C 612C PARAMETERS IV, N, V, AND X ARE THE SAME AS THE CORRESPONDING 613C ONES TO DSUMSL (WHICH SEE), EXCEPT THAT V CAN BE SHORTER (SINCE 614C THE PART OF V THAT DSUMSL USES FOR STORING G IS NOT NEEDED). 615C MOREOVER, COMPARED WITH DSUMSL, IV(1) MAY HAVE THE TWO ADDITIONAL 616C OUTPUT VALUES 1 AND 2, WHICH ARE EXPLAINED BELOW, AS IS THE USE 617C OF IV(TOOBIG) AND IV(NFGCAL). THE VALUE IV(G), WHICH IS AN 618C OUTPUT VALUE FROM DSUMSL (AND DSMSNO), IS NOT REFERENCED BY 619C DSUMIT OR THE SUBROUTINES IT CALLS. 620C FX AND G NEED NOT HAVE BEEN INITIALIZED WHEN DSUMIT IS CALLED 621C WITH IV(1) = 12, 13, OR 14. 622C 623C IV(1) = 1 MEANS THE CALLER SHOULD SET FX TO F(X), THE FUNCTION VALUE 624C AT X, AND CALL DSUMIT AGAIN, HAVING CHANGED NONE OF THE 625C OTHER PARAMETERS. AN EXCEPTION OCCURS IF F(X) CANNOT BE 626C (E.G. IF OVERFLOW WOULD OCCUR), WHICH MAY HAPPEN BECAUSE 627C OF AN OVERSIZED STEP. IN THIS CASE THE CALLER SHOULD SET 628C IV(TOOBIG) = IV(2) TO 1, WHICH WILL CAUSE DSUMIT TO IG- 629C NORE FX AND TRY A SMALLER STEP. THE PARAMETER NF THAT 630C DSUMSL PASSES TO CALCF (FOR POSSIBLE USE BY CALCG) IS A 631C COPY OF IV(NFCALL) = IV(6). 632C IV(1) = 2 MEANS THE CALLER SHOULD SET G TO G(X), THE GRADIENT VECTOR 633C OF F AT X, AND CALL DSUMIT AGAIN, HAVING CHANGED NONE OF 634C THE OTHER PARAMETERS EXCEPT POSSIBLY THE SCALE VECTOR D 635C WHEN IV(DTYPE) = 0. THE PARAMETER NF THAT DSUMSL PASSES 636C TO CALCG IS IV(NFGCAL) = IV(7). IF G(X) CANNOT BE 637C EVALUATED, THEN THE CALLER MAY SET IV(NFGCAL) TO 0, IN 638C WHICH CASE DSUMIT WILL RETURN WITH IV(1) = 65. 639C. 640C *** GENERAL *** 641C 642C CODED BY DAVID M. GAY (DECEMBER 1979). REVISED SEPT. 1982. 643C THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH SUPPORTED 644C IN PART BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS 645C MCS-7600324 AND MCS-7906671. 646C 647C (SEE DSUMSL FOR REFERENCES.) 648C 649C+++++++++++++++++++++++++++ DECLARATIONS ++++++++++++++++++++++++++++ 650C 651C *** LOCAL VARIABLES *** 652C 653 INTEGER DG1, DUMMY, G01, I, K, L, LSTGST, NN1O2, NWTST1, STEP1, 654 1 TEMP1, W, X01, Z 655 DOUBLE PRECISION T 656C 657C *** CONSTANTS *** 658C 659 DOUBLE PRECISION NEGONE, ONE, ZERO 660C 661C *** NO INTRINSIC FUNCTIONS *** 662C 663C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** 664C 665 EXTERNAL DASSST, DDBDOG, DDEFLT, DITSUM, DLITVM, DLIVMU, 666 1 DLTVMU, DLUPDT, DLVMUL, DPARCK, DSTOPX, DVAXPY, 667 2 DVSCPY, DVVMUP, DWZBFG 668 LOGICAL DSTOPX 669 DOUBLE PRECISION DDOT, DNRM2 670C 671C DASSST.... ASSESSES CANDIDATE STEP. 672C DDBDOG.... COMPUTES DOUBLE-DOGLEG (CANDIDATE) STEP. 673C DDEFLT.... SUPPLIES DEFAULT IV AND V INPUT COMPONENTS. 674C DITSUM.... PRINTS ITERATION SUMMARY AND INFO ON INITIAL AND FINAL X. 675C DLITVM... MULTIPLIES INVERSE TRANSPOSE OF LOWER TRIANGLE TIMES VECTOR. 676C DLIVMU... MULTIPLIES INVERSE OF LOWER TRIANGLE TIMES VECTOR. 677C DLTVMU... MULTIPLIES TRANSPOSE OF LOWER TRIANGLE TIMES VECTOR. 678C LUPDT.... UPDATES CHOLESKY FACTOR OF HESSIAN APPROXIMATION. 679C DLVMUL.... MULTIPLIES LOWER TRIANGLE TIMES VECTOR. 680C DPARCK.... CHECKS VALIDITY OF INPUT IV AND V VALUES. 681C DSTOPX.... RETURNS .TRUE. IF THE BREAK KEY HAS BEEN PRESSED. 682C DVAXPY.... COMPUTES SCALAR TIMES ONE VECTOR PLUS ANOTHER. 683C DVSCPY... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR. 684C DVVMUP... MULTIPLIES VECTOR BY VECTOR RAISED TO POWER (COMPONENTWISE). 685C DWZBFG... COMPUTES W AND Z FOR DLUPDT CORRESPONDING TO BFGS UPDATE. 686C 687C *** SUBSCRIPTS FOR IV AND V *** 688C 689 INTEGER CNVCOD, DG, DGNORM, DINIT, DSTNRM, DST0, F, F0, 690 1 GTHG, GTSTEP, G0, INCFAC, INITH, IRC, KAGQT, LMAT, 691 2 LMAX0, MODE, MODEL, MXFCAL, MXITER, NEXTV, NFCALL, NFGCAL, 692 3 NGCALL, NITER, NWTSTP, RADFAC, RADINC, RADIUS, RAD0, STEP, 693 4 STGLIM, STLSTG, TOOBIG, TUNER4, TUNER5, VNEED, XIRC, X0 694C 695C *** IV SUBSCRIPT VALUES *** 696C 697C/6 698C DATA CNVCOD/55/, DG/37/, G0/48/, INITH/25/, IRC/29/, KAGQT/33/, 699C 1 MODE/35/, MODEL/5/, MXFCAL/17/, MXITER/18/, NFCALL/6/, 700C 2 NFGCAL/7/, NGCALL/30/, NITER/31/, NWTSTP/34/, RADINC/8/, 701C 3 STEP/40/, STGLIM/11/, STLSTG/41/, TOOBIG/2/, XIRC/13/, X0/43/ 702C/7 703 PARAMETER (CNVCOD=55, DG=37, G0=48, INITH=25, IRC=29, KAGQT=33, 704 1 MODE=35, MODEL=5, MXFCAL=17, MXITER=18, NFCALL=6, 705 2 NFGCAL=7, NGCALL=30, NITER=31, NWTSTP=34, RADINC=8, 706 3 STEP=40, STGLIM=11, STLSTG=41, TOOBIG=2, XIRC=13, 707 4 X0=43) 708C/ 709C 710C *** V SUBSCRIPT VALUES *** 711C 712C/6 713C DATA DGNORM/1/, DINIT/38/, DSTNRM/2/, DST0/3/, F/10/, F0/13/, 714C 1 GTHG/44/, GTSTEP/4/, INCFAC/23/, LMAT/42/, LMAX0/35/, 715C 2 NEXTV/47/, RADFAC/16/, RADIUS/8/, RAD0/9/, TUNER4/29/, 716C 3 TUNER5/30/, VNEED/4/ 717C/7 718 PARAMETER (DGNORM=1, DINIT=38, DSTNRM=2, DST0=3, F=10, F0=13, 719 1 GTHG=44, GTSTEP=4, INCFAC=23, LMAT=42, LMAX0=35, 720 2 NEXTV=47, RADFAC=16, RADIUS=8, RAD0=9, TUNER4=29, 721 3 TUNER5=30, VNEED=4) 722C/ 723C 724C/6 725C DATA NEGONE/-1.D+0/, ONE/1.D+0/, ZERO/0.D+0/ 726C/7 727 PARAMETER (NEGONE=-1.D+0, ONE=1.D+0, ZERO=0.D+0) 728C/ 729C 730C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ 731C 732 I = IV(1) 733 IF (I .EQ. 1) GO TO 40 734 IF (I .EQ. 2) GO TO 50 735C 736C *** CHECK VALIDITY OF IV AND V INPUT VALUES *** 737C 738 IF (IV(1) .EQ. 0) CALL DDEFLT(2, IV, LIV, LV, V) 739 IV(VNEED) = IV(VNEED) + N*(N+13)/2 740 CALL DPARCK(2, D, IV, LIV, LV, N, V) 741 I = IV(1) - 2 742 IF (I .GT. 12) GO TO 999 743 GO TO (160, 160, 160, 160, 160, 160, 110, 80, 110, 10, 10, 20), I 744C 745C *** STORAGE ALLOCATION *** 746C 747 10 NN1O2 = N * (N + 1) / 2 748 L = IV(LMAT) 749 IV(X0) = L + NN1O2 750 IV(STEP) = IV(X0) + N 751 IV(STLSTG) = IV(STEP) + N 752 IV(G0) = IV(STLSTG) + N 753 IV(NWTSTP) = IV(G0) + N 754 IV(DG) = IV(NWTSTP) + N 755 IV(NEXTV) = IV(DG) + N 756 IF (IV(1) .NE. 13) GO TO 20 757 IV(1) = 14 758 GO TO 999 759C 760C *** INITIALIZATION *** 761C 762 20 IV(NITER) = 0 763 IV(NFCALL) = 1 764 IV(NGCALL) = 1 765 IV(NFGCAL) = 1 766 IV(MODE) = -1 767 IV(MODEL) = 1 768 IV(STGLIM) = 1 769 IV(TOOBIG) = 0 770 IV(CNVCOD) = 0 771 IV(RADINC) = 0 772 V(RAD0) = ZERO 773 IF (V(DINIT) .GE. ZERO) CALL DVSCPY(N, D, V(DINIT)) 774 IV(1) = 1 775 IF (IV(INITH) .NE. 1) GO TO 999 776C 777C *** SET THE INITIAL HESSIAN APPROXIMATION TO DIAG(D)**-2 *** 778C 779 CALL DVSCPY(NN1O2, V(L), ZERO) 780 K = L - 1 781 DO 30 I = 1, N 782 K = K + I 783 T = D(I) 784 IF (T .LE. ZERO) T = ONE 785 V(K) = T 786 30 CONTINUE 787 GO TO 999 788C 789 40 V(F) = FX 790 IF (IV(MODE) .GE. 0) GO TO 160 791 IV(1) = 2 792 IF (IV(TOOBIG) .EQ. 0) GO TO 999 793 IV(1) = 63 794 GO TO 270 795C 796C *** MAKE SURE GRADIENT COULD BE COMPUTED *** 797C 798 50 IF (IV(NFGCAL) .NE. 0) GO TO 60 799 IV(1) = 65 800 GO TO 270 801C 802 60 DG1 = IV(DG) 803 CALL DVVMUP(N, V(DG1), G, D, -1) 804 V(DGNORM) = DNRM2(N, V(DG1),1) 805C 806 IF (IV(CNVCOD) .NE. 0) GO TO 260 807 IF (IV(MODE) .EQ. 0) GO TO 220 808C 809C *** ALLOW FIRST STEP TO HAVE SCALED 2-NORM AT MOST V(LMAX0) *** 810C 811 V(RADFAC) = V(LMAX0) 812 V(DSTNRM) = ONE 813C 814 IV(MODE) = 0 815C 816C 817C----------------------------- MAIN LOOP ----------------------------- 818C 819C 820C *** PRINT ITERATION SUMMARY, CHECK ITERATION LIMIT *** 821C 822 70 CALL DITSUM(D, G, IV, LIV, LV, N, V, X) 823 80 K = IV(NITER) 824 IF (K .LT. IV(MXITER)) GO TO 90 825 IV(1) = 10 826 GO TO 270 827C 828C *** UPDATE RADIUS *** 829C 830 90 IV(NITER) = K + 1 831 V(RADIUS) = V(RADFAC) * V(DSTNRM) 832C 833C *** INITIALIZE FOR START OF NEXT ITERATION *** 834C 835 G01 = IV(G0) 836 X01 = IV(X0) 837 V(F0) = V(F) 838 IV(IRC) = 4 839 IV(KAGQT) = -1 840C 841C *** COPY X TO X0, G TO G0 *** 842C 843 CALL DCOPY(N, X,1,V(X01),1) 844 CALL DCOPY(N, G,1,V(G01),1) 845C 846C *** CHECK DSTOPX AND FUNCTION EVALUATION LIMIT *** 847C 848 100 IF (.NOT. DSTOPX(DUMMY)) GO TO 120 849 IV(1) = 11 850 GO TO 130 851C 852C *** COME HERE WHEN RESTARTING AFTER FUNC. EVAL. LIMIT OR DSTOPX. 853C 854 110 IF (V(F) .GE. V(F0)) GO TO 120 855 V(RADFAC) = ONE 856 K = IV(NITER) 857 GO TO 90 858C 859 120 IF (IV(NFCALL) .LT. IV(MXFCAL)) GO TO 140 860 IV(1) = 9 861 130 IF (V(F) .GE. V(F0)) GO TO 270 862C 863C *** IN CASE OF DSTOPX OR FUNCTION EVALUATION LIMIT WITH 864C *** IMPROVED V(F), EVALUATE THE GRADIENT AT X. 865C 866 IV(CNVCOD) = IV(1) 867 GO TO 210 868C 869C. . . . . . . . . . . . . COMPUTE CANDIDATE STEP . . . . . . . . . . 870C 871 140 STEP1 = IV(STEP) 872 DG1 = IV(DG) 873 NWTST1 = IV(NWTSTP) 874 IF (IV(KAGQT) .GE. 0) GO TO 150 875 L = IV(LMAT) 876 CALL DLIVMU(N, V(NWTST1), V(L), G) 877 CALL DLITVM(N, V(NWTST1), V(L), V(NWTST1)) 878 CALL DVVMUP(N, V(STEP1), V(NWTST1), D, 1) 879 V(DST0) = DNRM2(N, V(STEP1),1) 880 CALL DVVMUP(N, V(DG1), V(DG1), D, -1) 881 CALL DLTVMU(N, V(STEP1), V(L), V(DG1)) 882 V(GTHG) = DNRM2(N, V(STEP1),1) 883 IV(KAGQT) = 0 884 150 CALL DDBDOG(V(DG1), G, LV, N, V(NWTST1), V(STEP1), V) 885 IF (IV(IRC) .EQ. 6) GO TO 160 886C 887C *** COMPUTE F(X0 + STEP) *** 888C 889 X01 = IV(X0) 890 STEP1 = IV(STEP) 891 CALL DVAXPY(N, X, ONE, V(STEP1), V(X01)) 892 IV(NFCALL) = IV(NFCALL) + 1 893 IV(1) = 1 894 IV(TOOBIG) = 0 895 GO TO 999 896C 897C. . . . . . . . . . . . . ASSESS CANDIDATE STEP . . . . . . . . . . . 898C 899 160 STEP1 = IV(STEP) 900 LSTGST = IV(STLSTG) 901 X01 = IV(X0) 902 CALL DASSST(D, IV, N, V(STEP1), V(LSTGST), V, X, V(X01)) 903C 904 K = IV(IRC) 905 GO TO (170,200,200,200,170,180,190,190,190,190,190,190,250,220), K 906C 907C *** RECOMPUTE STEP WITH CHANGED RADIUS *** 908C 909 170 V(RADIUS) = V(RADFAC) * V(DSTNRM) 910 GO TO 100 911C 912C *** COMPUTE STEP OF LENGTH V(LMAX0) FOR SINGULAR CONVERGENCE TEST. 913C 914 180 V(RADIUS) = V(LMAX0) 915 GO TO 140 916C 917C *** CONVERGENCE OR FALSE CONVERGENCE *** 918C 919 190 IV(CNVCOD) = K - 4 920 IF (V(F) .GE. V(F0)) GO TO 260 921 IF (IV(XIRC) .EQ. 14) GO TO 260 922 IV(XIRC) = 14 923C 924C. . . . . . . . . . . . PROCESS ACCEPTABLE STEP . . . . . . . . . . . 925C 926 200 IF (IV(IRC) .NE. 3) GO TO 210 927 STEP1 = IV(STEP) 928 TEMP1 = IV(STLSTG) 929C 930C *** SET TEMP1 = HESSIAN * STEP FOR USE IN GRADIENT TESTS *** 931C 932 L = IV(LMAT) 933 CALL DLTVMU(N, V(TEMP1), V(L), V(STEP1)) 934 CALL DLVMUL(N, V(TEMP1), V(L), V(TEMP1)) 935C 936C *** COMPUTE GRADIENT *** 937C 938 210 IV(NGCALL) = IV(NGCALL) + 1 939 IV(1) = 2 940 GO TO 999 941C 942C *** INITIALIZATIONS -- G0 = G - G0, ETC. *** 943C 944 220 G01 = IV(G0) 945 CALL DVAXPY(N, V(G01), NEGONE, V(G01), G) 946 STEP1 = IV(STEP) 947 TEMP1 = IV(STLSTG) 948 IF (IV(IRC) .NE. 3) GO TO 240 949C 950C *** SET V(RADFAC) BY GRADIENT TESTS *** 951C 952C *** SET TEMP1 = DIAG(D)**-1 * (HESSIAN*STEP + (G(X0)-G(X))) *** 953C 954 CALL DVAXPY(N, V(TEMP1), NEGONE, V(G01), V(TEMP1)) 955 CALL DVVMUP(N, V(TEMP1), V(TEMP1), D, -1) 956C 957C *** DO GRADIENT TESTS *** 958C 959 IF (DNRM2(N, V(TEMP1),1) .LE. V(DGNORM) * V(TUNER4)) 960 1 GO TO 230 961 IF (DDOT(N, G,1,V(STEP1),1) 962 1 .GE. V(GTSTEP) * V(TUNER5)) GO TO 240 963 230 V(RADFAC) = V(INCFAC) 964C 965C *** UPDATE H, LOOP *** 966C 967 240 W = IV(NWTSTP) 968 Z = IV(X0) 969 L = IV(LMAT) 970 CALL DWZBFG(V(L), N, V(STEP1), V(W), V(G01), V(Z)) 971C 972C ** USE THE N-VECTORS STARTING AT V(STEP1) AND V(G01) FOR SCRATCH.. 973 CALL DLUPDT(V(TEMP1), V(STEP1), V(L), V(G01), V(L), N, V(W), V(Z)) 974 IV(1) = 2 975 GO TO 70 976C 977C. . . . . . . . . . . . . . MISC. DETAILS . . . . . . . . . . . . . . 978C 979C *** BAD PARAMETERS TO ASSESS *** 980C 981 250 IV(1) = 64 982C 983C *** PRINT SUMMARY OF FINAL ITERATION AND OTHER REQUESTED ITEMS *** 984C 985 260 IV(1) = IV(CNVCOD) 986 IV(CNVCOD) = 0 987 270 CALL DITSUM(D, G, IV, LIV, LV, N, V, X) 988C 989 999 RETURN 990C 991C *** LAST CARD OF DSUMIT FOLLOWS *** 992 END 993 SUBROUTINE DVAXPY(P, W, A, X, Y) 994 save 995C 996C *** SET W = A*X + Y -- W, X, Y = P-VECTORS, A = SCALAR *** 997C 998 INTEGER P 999 DOUBLE PRECISION A, W(P), X(P), Y(P) 1000C 1001 INTEGER I 1002C 1003 DO 10 I = 1, P 1004 10 W(I) = A*X(I) + Y(I) 1005 RETURN 1006 END 1007 SUBROUTINE DVDFLT(ALG, LV, V) 1008 save 1009C 1010C *** SUPPLY ***SOL (VERSION 2.3) DEFAULT VALUES TO V *** 1011C 1012C *** ALG = 1 MEANS REGRESSION CONSTANTS. 1013C *** ALG = 2 MEANS GENERAL UNCONSTRAINED OPTIMIZATION CONSTANTS. 1014C 1015 INTEGER ALG, LV 1016 DOUBLE PRECISION V(LV) 1017C/+ 1018 DOUBLE PRECISION DMAX1 1019C/ 1020 DOUBLE PRECISION D1MACH 1021C 1022 DOUBLE PRECISION MACHEP, MEPCRT, ONE, SQTEPS, THREE 1023C 1024C *** SUBSCRIPTS FOR V *** 1025C 1026 INTEGER AFCTOL, BIAS, COSMIN, DECFAC, DELTA0, DFAC, DINIT, DLTFDC, 1027 1 DLTFDJ, DTINIT, D0INIT, EPSLON, ETA0, FUZZ, HUBERC, 1028 2 INCFAC, LMAX0, LMAXS, PHMNFC, PHMXFC, RDFCMN, RDFCMX, 1029 3 RFCTOL, RLIMIT, RSPTOL, SCTOL, SIGMIN, TUNER1, TUNER2, 1030 4 TUNER3, TUNER4, TUNER5, XCTOL, XFTOL 1031C 1032C/6 1033C DATA ONE/1.D+0/, THREE/3.D+0/ 1034C/7 1035 PARAMETER (ONE=1.D+0, THREE=3.D+0) 1036C/ 1037C 1038C *** V SUBSCRIPT VALUES *** 1039C 1040C/6 1041C DATA AFCTOL/31/, BIAS/43/, COSMIN/47/, DECFAC/22/, DELTA0/44/, 1042C 1 DFAC/41/, DINIT/38/, DLTFDC/42/, DLTFDJ/43/, DTINIT/39/, 1043C 2 D0INIT/40/, EPSLON/19/, ETA0/42/, FUZZ/45/, HUBERC/48/, 1044C 3 INCFAC/23/, LMAX0/35/, LMAXS/36/, PHMNFC/20/, PHMXFC/21/, 1045C 4 RDFCMN/24/, RDFCMX/25/, RFCTOL/32/, RLIMIT/46/, RSPTOL/49/, 1046C 5 SCTOL/37/, SIGMIN/50/, TUNER1/26/, TUNER2/27/, TUNER3/28/, 1047C 6 TUNER4/29/, TUNER5/30/, XCTOL/33/, XFTOL/34/ 1048C/7 1049 PARAMETER (AFCTOL=31, BIAS=43, COSMIN=47, DECFAC=22, DELTA0=44, 1050 1 DFAC=41, DINIT=38, DLTFDC=42, DLTFDJ=43, DTINIT=39, 1051 2 D0INIT=40, EPSLON=19, ETA0=42, FUZZ=45, HUBERC=48, 1052 3 INCFAC=23, LMAX0=35, LMAXS=36, PHMNFC=20, PHMXFC=21, 1053 4 RDFCMN=24, RDFCMX=25, RFCTOL=32, RLIMIT=46, RSPTOL=49, 1054 5 SCTOL=37, SIGMIN=50, TUNER1=26, TUNER2=27, TUNER3=28, 1055 6 TUNER4=29, TUNER5=30, XCTOL=33, XFTOL=34) 1056C/ 1057C 1058C------------------------------- BODY -------------------------------- 1059C 1060 MACHEP = D1MACH(4) 1061 V(AFCTOL) = 1.D-20 1062 IF (MACHEP .GT. 1.D-10) V(AFCTOL) = MACHEP**2 1063 V(DECFAC) = 0.5D+0 1064 SQTEPS = DSQRT(D1MACH(4)) 1065 V(DFAC) = 0.6D+0 1066 V(DELTA0) = SQTEPS 1067 V(DTINIT) = 1.D-6 1068 MEPCRT = MACHEP ** (ONE/THREE) 1069 V(D0INIT) = 1.D+0 1070 V(EPSLON) = 0.1D+0 1071 V(INCFAC) = 2.D+0 1072 V(LMAX0) = 1.D+0 1073 V(LMAXS) = 1.D+0 1074 V(PHMNFC) = -0.1D+0 1075 V(PHMXFC) = 0.1D+0 1076 V(RDFCMN) = 0.1D+0 1077 V(RDFCMX) = 4.D+0 1078 V(RFCTOL) = DMAX1(1.D-10, MEPCRT**2) 1079 V(SCTOL) = V(RFCTOL) 1080 V(TUNER1) = 0.1D+0 1081 V(TUNER2) = 1.D-4 1082 V(TUNER3) = 0.75D+0 1083 V(TUNER4) = 0.5D+0 1084 V(TUNER5) = 0.75D+0 1085 V(XCTOL) = SQTEPS 1086 V(XFTOL) = 1.D+2 * MACHEP 1087C 1088 IF (ALG .GE. 2) GO TO 10 1089C 1090C *** REGRESSION VALUES 1091C 1092 V(COSMIN) = DMAX1(1.D-6, 1.D+2 * MACHEP) 1093 V(DINIT) = 0.D+0 1094 V(DLTFDC) = MEPCRT 1095 V(DLTFDJ) = SQTEPS 1096 V(FUZZ) = 1.5D+0 1097 V(HUBERC) = 0.7D+0 1098 V(RLIMIT) = DSQRT(D1MACH(2))*16. 1099 V(RSPTOL) = 1.D-2 1100 V(SIGMIN) = 1.D-4 1101 GO TO 999 1102C 1103C *** GENERAL OPTIMIZATION VALUES 1104C 1105 10 V(BIAS) = 0.8D+0 1106 V(DINIT) = -1.0D+0 1107 V(ETA0) = 1.0D+3 * MACHEP 1108C 1109 999 RETURN 1110C *** LAST CARD OF DVDFLT FOLLOWS *** 1111 END 1112 SUBROUTINE DVSCPY(P, Y, S) 1113 save 1114C 1115C *** SET P-VECTOR Y TO SCALAR S *** 1116C 1117 INTEGER P 1118 DOUBLE PRECISION S, Y(P) 1119C 1120 INTEGER I 1121C 1122 DO 10 I = 1, P 1123 10 Y(I) = S 1124 RETURN 1125 END 1126 SUBROUTINE DVVMUP(N, X, Y, Z, K) 1127 save 1128C 1129C *** SET X(I) = Y(I) * Z(I)**K, 1 .LE. I .LE. N (FOR K = 1 OR -1) *** 1130C 1131 INTEGER N, K 1132 DOUBLE PRECISION X(N), Y(N), Z(N) 1133 INTEGER I 1134C 1135 IF (K .GE. 0) GO TO 20 1136 DO 10 I = 1, N 1137 10 X(I) = Y(I) / Z(I) 1138 GO TO 999 1139C 1140 20 DO 30 I = 1, N 1141 30 X(I) = Y(I) * Z(I) 1142 999 RETURN 1143C *** LAST CARD OF DVVMUP FOLLOWS *** 1144 END 1145 SUBROUTINE DWZBFG (L, N, S, W, Y, Z) 1146 save 1147C 1148C *** COMPUTE Y AND Z FOR DLUPDT CORRESPONDING TO BFGS UPDATE. 1149C 1150 INTEGER N 1151 DOUBLE PRECISION L, S(N), W(N), Y(N), Z(N) 1152 DIMENSION L(N*(N+1)/2) 1153C 1154C-------------------------- PARAMETER USAGE -------------------------- 1155C 1156C L (I/O) CHOLESKY FACTOR OF HESSIAN, A LOWER TRIANG. MATRIX STORED 1157C COMPACTLY BY ROWS. 1158C N (INPUT) ORDER OF L AND LENGTH OF S, W, Y, Z. 1159C S (INPUT) THE STEP JUST TAKEN. 1160C W (OUTPUT) RIGHT SINGULAR VECTOR OF RANK 1 CORRECTION TO L. 1161C Y (INPUT) CHANGE IN GRADIENTS CORRESPONDING TO S. 1162C Z (OUTPUT) LEFT SINGULAR VECTOR OF RANK 1 CORRECTION TO L. 1163C 1164C------------------------------- NOTES ------------------------------- 1165C 1166C *** ALGORITHM NOTES *** 1167C 1168C WHEN S IS COMPUTED IN CERTAIN WAYS, E.G. BY GQTSTP OR 1169C DBLDOG, IT IS POSSIBLE TO SAVE N**2/2 OPERATIONS SINCE (L**T)*S 1170C OR L*(L**T)*S IS THEN KNOWN. 1171C IF THE BFGS UPDATE TO L*(L**T) WOULD REDUCE ITS DETERMINANT TO 1172C LESS THAN EPS TIMES ITS OLD VALUE, THEN THIS ROUTINE IN EFFECT 1173C REPLACES Y BY THETA*Y + (1 - THETA)*L*(L**T)*S, WHERE THETA 1174C (BETWEEN 0 AND 1) IS CHOSEN TO MAKE THE REDUCTION FACTOR = EPS. 1175C 1176C *** GENERAL *** 1177C 1178C CODED BY DAVID M. GAY (FALL 1979). 1179C THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH SUPPORTED 1180C BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS MCS-7600324 AND 1181C MCS-7906671. 1182C 1183C------------------------ EXTERNAL QUANTITIES ------------------------ 1184C 1185C *** FUNCTIONS AND SUBROUTINES CALLED *** 1186C 1187 EXTERNAL DLIVMU, DLTVMU 1188 DOUBLE PRECISION DDOT 1189C DLIVMU MULTIPLIES L**-1 TIMES A VECTOR. 1190C DLTVMU MULTIPLIES L**T TIMES A VECTOR. 1191C 1192C *** INTRINSIC FUNCTIONS *** 1193C/+ 1194 DOUBLE PRECISION DSQRT 1195C/ 1196C-------------------------- LOCAL VARIABLES -------------------------- 1197C 1198 INTEGER I 1199 DOUBLE PRECISION CS, CY, EPS, EPSRT, ONE, SHS, YS, THETA 1200C 1201C *** DATA INITIALIZATIONS *** 1202C 1203C/6 1204C DATA EPS/0.1D+0/, ONE/1.D+0/ 1205C/7 1206 PARAMETER (EPS=0.1D+0, ONE=1.D+0) 1207C/ 1208C 1209C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ 1210C 1211 CALL DLTVMU(N, W, L, S) 1212 SHS = DDOT(N, W,1,W,1) 1213 YS = DDOT(N, Y,1,S,1) 1214 IF (YS .GE. EPS*SHS) GO TO 10 1215 THETA = (ONE - EPS) * SHS / (SHS - YS) 1216 EPSRT = DSQRT(EPS) 1217 CY = THETA / (SHS * EPSRT) 1218 CS = (ONE + (THETA-ONE)/EPSRT) / SHS 1219 GO TO 20 1220 10 CY = ONE / (DSQRT(YS) * DSQRT(SHS)) 1221 CS = ONE / SHS 1222 20 CALL DLIVMU(N, Z, L, Y) 1223 DO 30 I = 1, N 1224 30 Z(I) = CY * Z(I) - CS * W(I) 1225C 1226 999 RETURN 1227C *** LAST CARD OF DWZBFG FOLLOWS *** 1228 END 1229 1230 SUBROUTINE DASSST(D, IV, P, STEP, STLSTG, V, X, X0) 1231 save 1232C 1233C *** ASSESS CANDIDATE STEP (***SOL VERSION 2.3) *** 1234C 1235 INTEGER P, IV(32) 1236 DOUBLE PRECISION D(P), STEP(P), STLSTG(P), V(37), X(P), X0(P) 1237C 1238C *** PURPOSE *** 1239C 1240C THIS SUBROUTINE IS CALLED BY AN UNCONSTRAINED MINIMIZATION 1241C ROUTINE TO ASSESS THE NEXT CANDIDATE STEP. IT MAY RECOMMEND ONE 1242C OF SEVERAL COURSES OF ACTION, SUCH AS ACCEPTING THE STEP, RECOM- 1243C PUTING IT USING THE SAME OR A NEW QUADRATIC MODEL, OR HALTING DUE 1244C TO CONVERGENCE OR FALSE CONVERGENCE. SEE THE RETURN CODE LISTING 1245C BELOW. 1246C 1247C-------------------------- PARAMETER USAGE -------------------------- 1248C 1249C IV (I/O) INTEGER PARAMETER AND SCRATCH VECTOR -- SEE DESCRIPTION 1250C BELOW OF IV VALUES REFERENCED. 1251C D (IN) SCALE VECTOR USED IN COMPUTING V(RELDX) -- SEE BELOW. 1252C P (IN) NUMBER OF PARAMETERS BEING OPTIMIZED. 1253C STEP (I/O) ON INPUT, STEP IS THE STEP TO BE ASSESSED. IT IS UN- 1254C CHANGED ON OUTPUT UNLESS A PREVIOUS STEP ACHIEVED A 1255C BETTER OBJECTIVE FUNCTION REDUCTION, IN WHICH CASE STLSTG 1256C WILL HAVE BEEN COPIED TO STEP. 1257C STLSTG (I/O) WHEN ASSESS RECOMMENDS RECOMPUTING STEP EVEN THOUGH THE 1258C CURRENT (OR A PREVIOUS) STEP YIELDS AN OBJECTIVE FUNC- 1259C TION DECREASE, IT SAVES IN STLSTG THE STEP THAT GAVE THE 1260C BEST FUNCTION REDUCTION SEEN SO FAR (IN THE CURRENT ITERA- 1261C TION). IF THE RECOMPUTED STEP YIELDS A LARGER FUNCTION 1262C VALUE, THEN STEP IS RESTORED FROM STLSTG AND 1263C X = X0 + STEP IS RECOMPUTED. 1264C V (I/O) REAL PARAMETER AND SCRATCH VECTOR -- SEE DESCRIPTION 1265C BELOW OF V VALUES REFERENCED. 1266C X (I/O) ON INPUT, X = X0 + STEP IS THE POINT AT WHICH THE OBJEC- 1267C TIVE FUNCTION HAS JUST BEEN EVALUATED. IF AN EARLIER 1268C STEP YIELDED A BIGGER FUNCTION DECREASE, THEN X IS 1269C RESTORED TO THE CORRESPONDING EARLIER VALUE. OTHERWISE, 1270C IF THE CURRENT STEP DOES NOT GIVE ANY FUNCTION DECREASE, 1271C THEN X IS RESTORED TO X0. 1272C X0 (IN) INITIAL OBJECTIVE FUNCTION PARAMETER VECTOR (AT THE 1273C START OF THE CURRENT ITERATION). 1274C 1275C *** IV VALUES REFERENCED *** 1276C 1277C IV(IRC) (I/O) ON INPUT FOR THE FIRST STEP TRIED IN A NEW ITERATION, 1278C IV(IRC) SHOULD BE SET TO 3 OR 4 (THE VALUE TO WHICH IT IS 1279C SET WHEN STEP IS DEFINITELY TO BE ACCEPTED). ON INPUT 1280C AFTER STEP HAS BEEN RECOMPUTED, IV(IRC) SHOULD BE 1281C UNCHANGED SINCE THE PREVIOUS RETURN OF ASSESS. 1282C ON OUTPUT, IV(IRC) IS A RETURN CODE HAVING ONE OF THE 1283C FOLLOWING VALUES... 1284C 1 = SWITCH MODELS OR TRY SMALLER STEP. 1285C 2 = SWITCH MODELS OR ACCEPT STEP. 1286C 3 = ACCEPT STEP AND DETERMINE V(RADFAC) BY GRADIENT 1287C TESTS. 1288C 4 = ACCEPT STEP, V(RADFAC) HAS BEEN DETERMINED. 1289C 5 = RECOMPUTE STEP (USING THE SAME MODEL). 1290C 6 = RECOMPUTE STEP WITH RADIUS = V(LMAXS) BUT DO NOT 1291C EVAULATE THE OBJECTIVE FUNCTION. 1292C 7 = X-CONVERGENCE (SEE V(XCTOL)). 1293C 8 = RELATIVE FUNCTION CONVERGENCE (SEE V(RFCTOL)). 1294C 9 = BOTH X- AND RELATIVE FUNCTION CONVERGENCE. 1295C 10 = ABSOLUTE FUNCTION CONVERGENCE (SEE V(AFCTOL)). 1296C 11 = SINGULAR CONVERGENCE (SEE V(LMAXS)). 1297C 12 = FALSE CONVERGENCE (SEE V(XFTOL)). 1298C 13 = IV(IRC) WAS OUT OF RANGE ON INPUT. 1299C RETURN CODE I HAS PRECDENCE OVER I+1 FOR I = 9, 10, 11. 1300C IV(MLSTGD) (I/O) SAVED VALUE OF IV(MODEL). 1301C IV(MODEL) (I/O) ON INPUT, IV(MODEL) SHOULD BE AN INTEGER IDENTIFYING 1302C THE CURRENT QUADRATIC MODEL OF THE OBJECTIVE FUNCTION. 1303C IF A PREVIOUS STEP YIELDED A BETTER FUNCTION REDUCTION, 1304C THEN IV(MODEL) WILL BE SET TO IV(MLSTGD) ON OUTPUT. 1305C IV(NFCALL) (IN) INVOCATION COUNT FOR THE OBJECTIVE FUNCTION. 1306C IV(NFGCAL) (I/O) VALUE OF IV(NFCALL) AT STEP THAT GAVE THE BIGGEST 1307C FUNCTION REDUCTION THIS ITERATION. IV(NFGCAL) REMAINS 1308C UNCHANGED UNTIL A FUNCTION REDUCTION IS OBTAINED. 1309C IV(RADINC) (I/O) THE NUMBER OF RADIUS INCREASES (OR MINUS THE NUMBER 1310C OF DECREASES) SO FAR THIS ITERATION. 1311C IV(RESTOR) (OUT) SET TO 0 UNLESS X AND V(F) HAVE BEEN RESTORED, IN 1312C WHICH CASE ASSESS SETS IV(RESTOR) = 1. 1313C IV(STAGE) (I/O) COUNT OF THE NUMBER OF MODELS TRIED SO FAR IN THE 1314C CURRENT ITERATION. 1315C IV(STGLIM) (IN) MAXIMUM NUMBER OF MODELS TO CONSIDER. 1316C IV(SWITCH) (OUT) SET TO 0 UNLESS A NEW MODEL IS BEING TRIED AND IT 1317C GIVES A SMALLER FUNCTION VALUE THAN THE PREVIOUS MODEL, 1318C IN WHICH CASE ASSESS SETS IV(SWITCH) = 1. 1319C IV(TOOBIG) (IN) IS NONZERO IF STEP WAS TOO BIG (E.G. IF IT CAUSED 1320C OVERFLOW). 1321C IV(XIRC) (I/O) VALUE THAT IV(IRC) WOULD HAVE IN THE ABSENCE OF 1322C CONVERGENCE, FALSE CONVERGENCE, AND OVERSIZED STEPS. 1323C 1324C *** V VALUES REFERENCED *** 1325C 1326C V(AFCTOL) (IN) ABSOLUTE FUNCTION CONVERGENCE TOLERANCE. IF THE 1327C ABSOLUTE VALUE OF THE CURRENT FUNCTION VALUE V(F) IS LESS 1328C THAN V(AFCTOL), THEN ASSESS RETURNS WITH IV(IRC) = 10. 1329C V(DECFAC) (IN) FACTOR BY WHICH TO DECREASE RADIUS WHEN IV(TOOBIG) IS 1330C NONZERO. 1331C V(DSTNRM) (IN) THE 2-NORM OF D*STEP. 1332C V(DSTSAV) (I/O) VALUE OF V(DSTNRM) ON SAVED STEP. 1333C V(DST0) (IN) THE 2-NORM OF D TIMES THE NEWTON STEP (WHEN DEFINED, 1334C I.E., FOR V(NREDUC) .GE. 0). 1335C V(F) (I/O) ON BOTH INPUT AND OUTPUT, V(F) IS THE OBJECTIVE FUNC- 1336C TION VALUE AT X. IF X IS RESTORED TO A PREVIOUS VALUE, 1337C THEN V(F) IS RESTORED TO THE CORRESPONDING VALUE. 1338C V(FDIF) (OUT) THE FUNCTION REDUCTION V(F0) - V(F) (FOR THE OUTPUT 1339C VALUE OF V(F) IF AN EARLIER STEP GAVE A BIGGER FUNCTION 1340C DECREASE, AND FOR THE INPUT VALUE OF V(F) OTHERWISE). 1341C V(FLSTGD) (I/O) SAVED VALUE OF V(F). 1342C V(F0) (IN) OBJECTIVE FUNCTION VALUE AT START OF ITERATION. 1343C V(GTSLST) (I/O) VALUE OF V(GTSTEP) ON SAVED STEP. 1344C V(GTSTEP) (IN) INNER PRODUCT BETWEEN STEP AND GRADIENT. 1345C V(INCFAC) (IN) MINIMUM FACTOR BY WHICH TO INCREASE RADIUS. 1346C V(LMAXS) (IN) MAXIMUM REASONABLE STEP SIZE (AND INITIAL STEP BOUND). 1347C IF THE ACTUAL FUNCTION DECREASE IS NO MORE THAN TWICE 1348C WHAT WAS PREDICTED, IF A RETURN WITH IV(IRC) = 7, 8, 9, 1349C OR 10 DOES NOT OCCUR, IF V(DSTNRM) .GT. V(LMAXS), AND IF 1350C V(PREDUC) .LE. V(SCTOL) * ABS(V(F0)), THEN ASSESS RE- 1351C TURNS WITH IV(IRC) = 11. IF SO DOING APPEARS WORTHWHILE, 1352C THEN ASSESS REPEATS THIS TEST WITH V(PREDUC) COMPUTED FOR 1353C A STEP OF LENGTH V(LMAXS) (BY A RETURN WITH IV(IRC) = 6). 1354C V(NREDUC) (I/O) FUNCTION REDUCTION PREDICTED BY QUADRATIC MODEL FOR 1355C NEWTON STEP. IF ASSESS IS CALLED WITH IV(IRC) = 6, I.E., 1356C IF V(PREDUC) HAS BEEN COMPUTED WITH RADIUS = V(LMAXS) FOR 1357C USE IN THE SINGULAR CONVERVENCE TEST, THEN V(NREDUC) IS 1358C SET TO -V(PREDUC) BEFORE THE LATTER IS RESTORED. 1359C V(PLSTGD) (I/O) VALUE OF V(PREDUC) ON SAVED STEP. 1360C V(PREDUC) (I/O) FUNCTION REDUCTION PREDICTED BY QUADRATIC MODEL FOR 1361C CURRENT STEP. 1362C V(RADFAC) (OUT) FACTOR TO BE USED IN DETERMINING THE NEW RADIUS, 1363C WHICH SHOULD BE V(RADFAC)*DST, WHERE DST IS EITHER THE 1364C OUTPUT VALUE OF V(DSTNRM) OR THE 2-NORM OF 1365C DIAG(NEWD)*STEP FOR THE OUTPUT VALUE OF STEP AND THE 1366C UPDATED VERSION, NEWD, OF THE SCALE VECTOR D. FOR 1367C IV(IRC) = 3, V(RADFAC) = 1.0 IS RETURNED. 1368C V(RDFCMN) (IN) MINIMUM VALUE FOR V(RADFAC) IN TERMS OF THE INPUT 1369C VALUE OF V(DSTNRM) -- SUGGESTED VALUE = 0.1. 1370C V(RDFCMX) (IN) MAXIMUM VALUE FOR V(RADFAC) -- SUGGESTED VALUE = 4.0. 1371C V(RELDX) (OUT) SCALED RELATIVE CHANGE IN X CAUSED BY STEP, COMPUTED 1372C BY FUNCTION DRELST AS 1373C MAX (D(I)*ABS(X(I)-X0(I)), 1 .LE. I .LE. P) / 1374C MAX (D(I)*(ABS(X(I))+ABS(X0(I))), 1 .LE. I .LE. P). 1375C IF AN ACCEPTABLE STEP IS RETURNED, THEN V(RELDX) IS COM- 1376C PUTED USING THE OUTPUT (POSSIBLY RESTORED) VALUES OF X 1377C AND STEP. OTHERWISE IT IS COMPUTED USING THE INPUT 1378C VALUES. 1379C V(RFCTOL) (IN) RELATIVE FUNCTION CONVERGENCE TOLERANCE. IF THE 1380C ACTUAL FUNCTION REDUCTION IS AT MOST TWICE WHAT WAS PRE- 1381C DICTED AND V(NREDUC) .LE. V(RFCTOL)*ABS(V(F0)), THEN 1382C ASSESS RETURNS WITH IV(IRC) = 8 OR 9. 1383C V(STPPAR) (IN) MARQUARDT PARAMETER -- 0 MEANS FULL NEWTON STEP. 1384C V(TUNER1) (IN) TUNING CONSTANT USED TO DECIDE IF THE FUNCTION 1385C REDUCTION WAS MUCH LESS THAN EXPECTED. SUGGESTED 1386C VALUE = 0.1. 1387C V(TUNER2) (IN) TUNING CONSTANT USED TO DECIDE IF THE FUNCTION 1388C REDUCTION WAS LARGE ENOUGH TO ACCEPT STEP. SUGGESTED 1389C VALUE = 10**-4. 1390C V(TUNER3) (IN) TUNING CONSTANT USED TO DECIDE IF THE RADIUS 1391C SHOULD BE INCREASED. SUGGESTED VALUE = 0.75. 1392C V(XCTOL) (IN) X-CONVERGENCE CRITERION. IF STEP IS A NEWTON STEP 1393C (V(STPPAR) = 0) HAVING V(RELDX) .LE. V(XCTOL) AND GIVING 1394C AT MOST TWICE THE PREDICTED FUNCTION DECREASE, THEN 1395C ASSESS RETURNS IV(IRC) = 7 OR 9. 1396C V(XFTOL) (IN) FALSE CONVERGENCE TOLERANCE. IF STEP GAVE NO OR ONLY 1397C A SMALL FUNCTION DECREASE AND V(RELDX) .LE. V(XFTOL), 1398C THEN ASSESS RETURNS WITH IV(IRC) = 12. 1399C 1400C------------------------------- NOTES ------------------------------- 1401C 1402C *** APPLICATION AND USAGE RESTRICTIONS *** 1403C 1404C THIS ROUTINE IS CALLED AS PART OF THE NL2SOL (NONLINEAR 1405C LEAST-SQUARES) PACKAGE. IT MAY BE USED IN ANY UNCONSTRAINED 1406C MINIMIZATION SOLVER THAT USES DOGLEG, GOLDFELD-QUANDT-TROTTER, 1407C OR LEVENBERG-MARQUARDT STEPS. 1408C 1409C *** ALGORITHM NOTES *** 1410C 1411C SEE (1) FOR FURTHER DISCUSSION OF THE ASSESSING AND MODEL 1412C SWITCHING STRATEGIES. WHILE NL2SOL CONSIDERS ONLY TWO MODELS, 1413C ASSESS IS DESIGNED TO HANDLE ANY NUMBER OF MODELS. 1414C 1415C *** USAGE NOTES *** 1416C 1417C ON THE FIRST CALL OF AN ITERATION, ONLY THE I/O VARIABLES 1418C STEP, X, IV(IRC), IV(MODEL), V(F), V(DSTNRM), V(GTSTEP), AND 1419C V(PREDUC) NEED HAVE BEEN INITIALIZED. BETWEEN CALLS, NO I/O 1420C VALUES EXECPT STEP, X, IV(MODEL), V(F) AND THE STOPPING TOLER- 1421C ANCES SHOULD BE CHANGED. 1422C AFTER A RETURN FOR CONVERGENCE OR FALSE CONVERGENCE, ONE CAN 1423C CHANGE THE STOPPING TOLERANCES AND CALL ASSESS AGAIN, IN WHICH 1424C CASE THE STOPPING TESTS WILL BE REPEATED. 1425C 1426C *** REFERENCES *** 1427C 1428C (1) DENNIS, J.E., JR., GAY, D.M., AND WELSCH, R.E. (1981), 1429C AN ADAPTIVE NONLINEAR LEAST-SQUARES ALGORITHM, 1430C ACM TRANS. MATH. SOFTWARE, VOL. 7, NO. 3. 1431C 1432C (2) POWELL, M.J.D. (1970) A FORTRAN SUBROUTINE FOR SOLVING 1433C SYSTEMS OF NONLINEAR ALGEBRAIC EQUATIONS, IN NUMERICAL 1434C METHODS FOR NONLINEAR ALGEBRAIC EQUATIONS, EDITED BY 1435C P. RABINOWITZ, GORDON AND BREACH, LONDON. 1436C 1437C *** HISTORY *** 1438C 1439C JOHN DENNIS DESIGNED MUCH OF THIS ROUTINE, STARTING WITH 1440C IDEAS IN (2). ROY WELSCH SUGGESTED THE MODEL SWITCHING STRATEGY. 1441C DAVID GAY AND STEPHEN PETERS CAST THIS SUBROUTINE INTO A MORE 1442C PORTABLE FORM (WINTER 1977), AND DAVID GAY CAST IT INTO ITS 1443C PRESENT FORM (FALL 1978). 1444C 1445C *** GENERAL *** 1446C 1447C THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH 1448C SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS 1449C MCS-7600324, DCR75-10143, 76-14311DSS, MCS76-11989, AND 1450C MCS-7906671. 1451C 1452C------------------------ EXTERNAL QUANTITIES ------------------------ 1453C 1454C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** 1455C 1456 EXTERNAL DRELST 1457 DOUBLE PRECISION DRELST 1458C 1459C DRELST... COMPUTES V(RELDX) = RELATIVE STEP SIZE. 1460C 1461C *** INTRINSIC FUNCTIONS *** 1462C/+ 1463 DOUBLE PRECISION DABS, DMAX1 1464C/ 1465C *** NO COMMON BLOCKS *** 1466C 1467C-------------------------- LOCAL VARIABLES -------------------------- 1468C 1469 LOGICAL GOODX 1470 INTEGER I, NFC 1471 DOUBLE PRECISION EMAX, EMAXS, GTS, HALF, ONE, RELDX1, RFAC1, TWO, 1472 1 XMAX, ZERO 1473C 1474C *** SUBSCRIPTS FOR IV AND V *** 1475C 1476 INTEGER AFCTOL, DECFAC, DSTNRM, DSTSAV, DST0, F, FDIF, FLSTGD, F0, 1477 1 GTSLST, GTSTEP, INCFAC, IRC, LMAXS, MLSTGD, MODEL, NFCALL, 1478 2 NFGCAL, NREDUC, PLSTGD, PREDUC, RADFAC, RADINC, RDFCMN, 1479 3 RDFCMX, RELDX, RESTOR, RFCTOL, SCTOL, STAGE, STGLIM, 1480 4 STPPAR, SWITCH, TOOBIG, TUNER1, TUNER2, TUNER3, XCTOL, 1481 5 XFTOL, XIRC 1482C 1483C *** DATA INITIALIZATIONS *** 1484C 1485C/6 1486C DATA HALF/0.5D+0/, ONE/1.D+0/, TWO/2.D+0/, ZERO/0.D+0/ 1487C/7 1488 PARAMETER (HALF=0.5D+0, ONE=1.D+0, TWO=2.D+0, ZERO=0.D+0) 1489C/ 1490C 1491C/6 1492C DATA IRC/29/, MLSTGD/32/, MODEL/5/, NFCALL/6/, NFGCAL/7/, 1493C 1 RADINC/8/, RESTOR/9/, STAGE/10/, STGLIM/11/, SWITCH/12/, 1494C 2 TOOBIG/2/, XIRC/13/ 1495C/7 1496 PARAMETER (IRC=29, MLSTGD=32, MODEL=5, NFCALL=6, NFGCAL=7, 1497 1 RADINC=8, RESTOR=9, STAGE=10, STGLIM=11, SWITCH=12, 1498 2 TOOBIG=2, XIRC=13) 1499C/ 1500C/6 1501C DATA AFCTOL/31/, DECFAC/22/, DSTNRM/2/, DST0/3/, DSTSAV/18/, 1502C 1 F/10/, FDIF/11/, FLSTGD/12/, F0/13/, GTSLST/14/, GTSTEP/4/, 1503C 2 INCFAC/23/, LMAXS/36/, NREDUC/6/, PLSTGD/15/, PREDUC/7/, 1504C 3 RADFAC/16/, RDFCMN/24/, RDFCMX/25/, RELDX/17/, RFCTOL/32/, 1505C 4 SCTOL/37/, STPPAR/5/, TUNER1/26/, TUNER2/27/, TUNER3/28/, 1506C 5 XCTOL/33/, XFTOL/34/ 1507C/7 1508 PARAMETER (AFCTOL=31, DECFAC=22, DSTNRM=2, DST0=3, DSTSAV=18, 1509 1 F=10, FDIF=11, FLSTGD=12, F0=13, GTSLST=14, GTSTEP=4, 1510 2 INCFAC=23, LMAXS=36, NREDUC=6, PLSTGD=15, PREDUC=7, 1511 3 RADFAC=16, RDFCMN=24, RDFCMX=25, RELDX=17, RFCTOL=32, 1512 4 SCTOL=37, STPPAR=5, TUNER1=26, TUNER2=27, TUNER3=28, 1513 5 XCTOL=33, XFTOL=34) 1514C/ 1515C 1516C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ 1517C 1518 NFC = IV(NFCALL) 1519 IV(SWITCH) = 0 1520 IV(RESTOR) = 0 1521 RFAC1 = ONE 1522 GOODX = .TRUE. 1523 I = IV(IRC) 1524 IF (I .GE. 1 .AND. I .LE. 12) 1525 1 GO TO (20,30,10,10,40,280,220,220,220,220,220,170), I 1526 IV(IRC) = 13 1527 GO TO 999 1528C 1529C *** INITIALIZE FOR NEW ITERATION *** 1530C 1531 10 IV(STAGE) = 1 1532 IV(RADINC) = 0 1533 V(FLSTGD) = V(F0) 1534 IF (IV(TOOBIG) .EQ. 0) GO TO 90 1535 IV(STAGE) = -1 1536 IV(XIRC) = I 1537 GO TO 60 1538C 1539C *** STEP WAS RECOMPUTED WITH NEW MODEL OR SMALLER RADIUS *** 1540C *** FIRST DECIDE WHICH *** 1541C 1542 20 IF (IV(MODEL) .NE. IV(MLSTGD)) GO TO 30 1543C *** OLD MODEL RETAINED, SMALLER RADIUS TRIED *** 1544C *** DO NOT CONSIDER ANY MORE NEW MODELS THIS ITERATION *** 1545 IV(STAGE) = IV(STGLIM) 1546 IV(RADINC) = -1 1547 GO TO 90 1548C 1549C *** A NEW MODEL IS BEING TRIED. DECIDE WHETHER TO KEEP IT. *** 1550C 1551 30 IV(STAGE) = IV(STAGE) + 1 1552C 1553C *** NOW WE ADD THE POSSIBILTIY THAT STEP WAS RECOMPUTED WITH *** 1554C *** THE SAME MODEL, PERHAPS BECAUSE OF AN OVERSIZED STEP. *** 1555C 1556 40 IF (IV(STAGE) .GT. 0) GO TO 50 1557C 1558C *** STEP WAS RECOMPUTED BECAUSE IT WAS TOO BIG. *** 1559C 1560 IF (IV(TOOBIG) .NE. 0) GO TO 60 1561C 1562C *** RESTORE IV(STAGE) AND PICK UP WHERE WE LEFT OFF. *** 1563C 1564 IV(STAGE) = -IV(STAGE) 1565 I = IV(XIRC) 1566 GO TO (20, 30, 90, 90, 70), I 1567C 1568 50 IF (IV(TOOBIG) .EQ. 0) GO TO 70 1569C 1570C *** HANDLE OVERSIZE STEP *** 1571C 1572 IF (IV(RADINC) .GT. 0) GO TO 80 1573 IV(STAGE) = -IV(STAGE) 1574 IV(XIRC) = IV(IRC) 1575C 1576 60 V(RADFAC) = V(DECFAC) 1577 IV(RADINC) = IV(RADINC) - 1 1578 IV(IRC) = 5 1579 GO TO 999 1580C 1581 70 IF (V(F) .LT. V(FLSTGD)) GO TO 90 1582C 1583C *** THE NEW STEP IS A LOSER. RESTORE OLD MODEL. *** 1584C 1585 IF (IV(MODEL) .EQ. IV(MLSTGD)) GO TO 80 1586 IV(MODEL) = IV(MLSTGD) 1587 IV(SWITCH) = 1 1588C 1589C *** RESTORE STEP, ETC. ONLY IF A PREVIOUS STEP DECREASED V(F). 1590C 1591 80 IF (V(FLSTGD) .GE. V(F0)) GO TO 90 1592 IV(RESTOR) = 1 1593 V(F) = V(FLSTGD) 1594 V(PREDUC) = V(PLSTGD) 1595 V(GTSTEP) = V(GTSLST) 1596 IF (IV(SWITCH) .EQ. 0) RFAC1 = V(DSTNRM) / V(DSTSAV) 1597 V(DSTNRM) = V(DSTSAV) 1598 NFC = IV(NFGCAL) 1599 GOODX = .FALSE. 1600C 1601C 1602C *** COMPUTE RELATIVE CHANGE IN X BY CURRENT STEP *** 1603C 1604 90 RELDX1 = DRELST(P, D, X, X0) 1605C 1606C *** RESTORE X AND STEP IF NECESSARY *** 1607C 1608 IF (GOODX) GO TO 110 1609 DO 100 I = 1, P 1610 STEP(I) = STLSTG(I) 1611 X(I) = X0(I) + STLSTG(I) 1612 100 CONTINUE 1613C 1614 110 V(FDIF) = V(F0) - V(F) 1615 IF (V(FDIF) .GT. V(TUNER2) * V(PREDUC)) GO TO 140 1616C 1617C *** NO (OR ONLY A TRIVIAL) FUNCTION DECREASE 1618C *** -- SO TRY NEW MODEL OR SMALLER RADIUS 1619C 1620 V(RELDX) = RELDX1 1621 IF (V(F) .LT. V(F0)) GO TO 120 1622 IV(MLSTGD) = IV(MODEL) 1623 V(FLSTGD) = V(F) 1624 V(F) = V(F0) 1625 CALL DCOPY(P,X0,1,X,1) 1626 IV(RESTOR) = 1 1627 GO TO 130 1628 120 IV(NFGCAL) = NFC 1629 130 IV(IRC) = 1 1630 IF (IV(STAGE) .LT. IV(STGLIM)) GO TO 160 1631 IV(IRC) = 5 1632 IV(RADINC) = IV(RADINC) - 1 1633 GO TO 160 1634C 1635C *** NONTRIVIAL FUNCTION DECREASE ACHIEVED *** 1636C 1637 140 IV(NFGCAL) = NFC 1638 RFAC1 = ONE 1639 IF (GOODX) V(RELDX) = RELDX1 1640 V(DSTSAV) = V(DSTNRM) 1641 IF (V(FDIF) .GT. V(PREDUC)*V(TUNER1)) GO TO 190 1642C 1643C *** DECREASE WAS MUCH LESS THAN PREDICTED -- EITHER CHANGE MODELS 1644C *** OR ACCEPT STEP WITH DECREASED RADIUS. 1645C 1646 IF (IV(STAGE) .GE. IV(STGLIM)) GO TO 150 1647C *** CONSIDER SWITCHING MODELS *** 1648 IV(IRC) = 2 1649 GO TO 160 1650C 1651C *** ACCEPT STEP WITH DECREASED RADIUS *** 1652C 1653 150 IV(IRC) = 4 1654C 1655C *** SET V(RADFAC) TO FLETCHER*S DECREASE FACTOR *** 1656C 1657 160 IV(XIRC) = IV(IRC) 1658 EMAX = V(GTSTEP) + V(FDIF) 1659 V(RADFAC) = HALF * RFAC1 1660 IF (EMAX .LT. V(GTSTEP)) V(RADFAC) = RFAC1 * DMAX1(V(RDFCMN), 1661 1 HALF * V(GTSTEP)/EMAX) 1662C 1663C *** DO FALSE CONVERGENCE TEST *** 1664C 1665 170 IF (V(RELDX) .LE. V(XFTOL)) GO TO 180 1666 IV(IRC) = IV(XIRC) 1667 IF (V(F) .LT. V(F0)) GO TO 200 1668 GO TO 230 1669C 1670 180 IV(IRC) = 12 1671 GO TO 240 1672C 1673C *** HANDLE GOOD FUNCTION DECREASE *** 1674C 1675 190 IF (V(FDIF) .LT. (-V(TUNER3) * V(GTSTEP))) GO TO 210 1676C 1677C *** INCREASING RADIUS LOOKS WORTHWHILE. SEE IF WE JUST 1678C *** RECOMPUTED STEP WITH A DECREASED RADIUS OR RESTORED STEP 1679C *** AFTER RECOMPUTING IT WITH A LARGER RADIUS. 1680C 1681 IF (IV(RADINC) .LT. 0) GO TO 210 1682 IF (IV(RESTOR) .EQ. 1) GO TO 210 1683C 1684C *** WE DID NOT. TRY A LONGER STEP UNLESS THIS WAS A NEWTON 1685C *** STEP. 1686C 1687 V(RADFAC) = V(RDFCMX) 1688 GTS = V(GTSTEP) 1689 IF (V(FDIF) .LT. (HALF/V(RADFAC) - ONE) * GTS) 1690 1 V(RADFAC) = DMAX1(V(INCFAC), HALF*GTS/(GTS + V(FDIF))) 1691 IV(IRC) = 4 1692 IF (V(STPPAR) .EQ. ZERO) GO TO 230 1693C *** STEP WAS NOT A NEWTON STEP. RECOMPUTE IT WITH 1694C *** A LARGER RADIUS. 1695 IV(IRC) = 5 1696 IV(RADINC) = IV(RADINC) + 1 1697C 1698C *** SAVE VALUES CORRESPONDING TO GOOD STEP *** 1699C 1700 200 V(FLSTGD) = V(F) 1701 IV(MLSTGD) = IV(MODEL) 1702 CALL DCOPY(P, STEP,1,STLSTG,1) 1703 V(DSTSAV) = V(DSTNRM) 1704 IV(NFGCAL) = NFC 1705 V(PLSTGD) = V(PREDUC) 1706 V(GTSLST) = V(GTSTEP) 1707 GO TO 230 1708C 1709C *** ACCEPT STEP WITH RADIUS UNCHANGED *** 1710C 1711 210 V(RADFAC) = ONE 1712 IV(IRC) = 3 1713 GO TO 230 1714C 1715C *** COME HERE FOR A RESTART AFTER CONVERGENCE *** 1716C 1717 220 IV(IRC) = IV(XIRC) 1718 IF (V(DSTSAV) .GE. ZERO) GO TO 240 1719 IV(IRC) = 12 1720 GO TO 240 1721C 1722C *** PERFORM CONVERGENCE TESTS *** 1723C 1724 230 IV(XIRC) = IV(IRC) 1725 240 IF (DABS(V(F)) .LT. V(AFCTOL)) IV(IRC) = 10 1726 IF (HALF * V(FDIF) .GT. V(PREDUC)) GO TO 999 1727 EMAX = V(RFCTOL) * DABS(V(F0)) 1728 EMAXS = V(SCTOL) * DABS(V(F0)) 1729 IF (V(DSTNRM) .GT. V(LMAXS) .AND. V(PREDUC) .LE. EMAXS) 1730 1 IV(IRC) = 11 1731 IF (V(DST0) .LT. ZERO) GO TO 250 1732 I = 0 1733 IF ((V(NREDUC) .GT. ZERO .AND. V(NREDUC) .LE. EMAX) .OR. 1734 1 (V(NREDUC) .EQ. ZERO. AND. V(PREDUC) .EQ. ZERO)) I = 2 1735 IF (V(STPPAR) .EQ. ZERO .AND. V(RELDX) .LE. V(XCTOL) 1736 1 .AND. GOODX) I = I + 1 1737 IF (I .GT. 0) IV(IRC) = I + 6 1738C 1739C *** CONSIDER RECOMPUTING STEP OF LENGTH V(LMAXS) FOR SINGULAR 1740C *** CONVERGENCE TEST. 1741C 1742 250 IF (IV(IRC) .GT. 5 .AND. IV(IRC) .NE. 12) GO TO 999 1743 IF (V(DSTNRM) .GT. V(LMAXS)) GO TO 260 1744 IF (V(PREDUC) .GE. EMAXS) GO TO 999 1745 IF (V(DST0) .LE. ZERO) GO TO 270 1746 IF (HALF * V(DST0) .LE. V(LMAXS)) GO TO 999 1747 GO TO 270 1748 260 IF (HALF * V(DSTNRM) .LE. V(LMAXS)) GO TO 999 1749 XMAX = V(LMAXS) / V(DSTNRM) 1750 IF (XMAX * (TWO - XMAX) * V(PREDUC) .GE. EMAXS) GO TO 999 1751 270 IF (V(NREDUC) .LT. ZERO) GO TO 290 1752C 1753C *** RECOMPUTE V(PREDUC) FOR USE IN SINGULAR CONVERGENCE TEST *** 1754C 1755 V(GTSLST) = V(GTSTEP) 1756 V(DSTSAV) = V(DSTNRM) 1757 IF (IV(IRC) .EQ. 12) V(DSTSAV) = -V(DSTSAV) 1758 V(PLSTGD) = V(PREDUC) 1759 IV(IRC) = 6 1760 CALL DCOPY(P, STEP,1,STLSTG,1) 1761 GO TO 999 1762C 1763C *** PERFORM SINGULAR CONVERGENCE TEST WITH RECOMPUTED V(PREDUC) *** 1764C 1765 280 V(GTSTEP) = V(GTSLST) 1766 V(DSTNRM) = DABS(V(DSTSAV)) 1767 CALL DCOPY(P, STLSTG,1,STEP,1) 1768 IV(IRC) = IV(XIRC) 1769 IF (V(DSTSAV) .LE. ZERO) IV(IRC) = 12 1770 V(NREDUC) = -V(PREDUC) 1771 V(PREDUC) = V(PLSTGD) 1772 290 IF (-V(NREDUC) .LE. V(RFCTOL) * DABS(V(F0))) IV(IRC) = 11 1773C 1774 999 RETURN 1775C 1776C *** LAST CARD OF ASSESS FOLLOWS *** 1777 END 1778 SUBROUTINE DDBDOG(DIG, G, LV, N, NWTSTP, STEP, V) 1779 save 1780C 1781C *** COMPUTE DOUBLE DOGLEG STEP *** 1782C 1783C *** PARAMETER DECLARATIONS *** 1784C 1785 INTEGER LV, N 1786 DOUBLE PRECISION DIG(N), G(N), NWTSTP(N), STEP(N), V(LV) 1787C 1788C *** PURPOSE *** 1789C 1790C THIS SUBROUTINE COMPUTES A CANDIDATE STEP (FOR USE IN AN UNCON- 1791C STRAINED MINIMIZATION CODE) BY THE DOUBLE DOGLEG ALGORITHM OF 1792C DENNIS AND MEI (REF. 1), WHICH IS A VARIATION ON POWELL*S DOGLEG 1793C SCHEME (REF. 2, P. 95). 1794C 1795C-------------------------- PARAMETER USAGE -------------------------- 1796C 1797C DIG (INPUT) DIAG(D)**-2 * G -- SEE ALGORITHM NOTES. 1798C G (INPUT) THE CURRENT GRADIENT VECTOR. 1799C LV (INPUT) LENGTH OF V. 1800C N (INPUT) NUMBER OF COMPONENTS IN DIG, G, NWTSTP, AND STEP. 1801C NWTSTP (INPUT) NEGATIVE NEWTON STEP -- SEE ALGORITHM NOTES. 1802C STEP (OUTPUT) THE COMPUTED STEP. 1803C V (I/O) VALUES ARRAY, THE FOLLOWING COMPONENTS OF WHICH ARE 1804C USED HERE... 1805C V(BIAS) (INPUT) BIAS FOR RELAXED NEWTON STEP, WHICH IS V(BIAS) OF 1806C THE WAY FROM THE FULL NEWTON TO THE FULLY RELAXED NEWTON 1807C STEP. RECOMMENDED VALUE = 0.8 . 1808C V(DGNORM) (INPUT) 2-NORM OF DIAG(D)**-1 * G -- SEE ALGORITHM NOTES. 1809C V(DSTNRM) (OUTPUT) 2-NORM OF DIAG(D) * STEP, WHICH IS V(RADIUS) 1810C UNLESS V(STPPAR) = 0 -- SEE ALGORITHM NOTES. 1811C V(DST0) (INPUT) 2-NORM OF DIAG(D) * NWTSTP -- SEE ALGORITHM NOTES. 1812C V(GRDFAC) (OUTPUT) THE COEFFICIENT OF DIG IN THE STEP RETURNED -- 1813C STEP(I) = V(GRDFAC)*DIG(I) + V(NWTFAC)*NWTSTP(I). 1814C V(GTHG) (INPUT) SQUARE-ROOT OF (DIG**T) * (HESSIAN) * DIG -- SEE 1815C ALGORITHM NOTES. 1816C V(GTSTEP) (OUTPUT) INNER PRODUCT BETWEEN G AND STEP. 1817C V(NREDUC) (OUTPUT) FUNCTION REDUCTION PREDICTED FOR THE FULL NEWTON 1818C STEP. 1819C V(NWTFAC) (OUTPUT) THE COEFFICIENT OF NWTSTP IN THE STEP RETURNED -- 1820C SEE V(GRDFAC) ABOVE. 1821C V(PREDUC) (OUTPUT) FUNCTION REDUCTION PREDICTED FOR THE STEP RETURNED. 1822C V(RADIUS) (INPUT) THE TRUST REGION RADIUS. D TIMES THE STEP RETURNED 1823C HAS 2-NORM V(RADIUS) UNLESS V(STPPAR) = 0. 1824C V(STPPAR) (OUTPUT) CODE TELLING HOW STEP WAS COMPUTED... 0 MEANS A 1825C FULL NEWTON STEP. BETWEEN 0 AND 1 MEANS V(STPPAR) OF THE 1826C WAY FROM THE NEWTON TO THE RELAXED NEWTON STEP. BETWEEN 1827C 1 AND 2 MEANS A TRUE DOUBLE DOGLEG STEP, V(STPPAR) - 1 OF 1828C THE WAY FROM THE RELAXED NEWTON TO THE CAUCHY STEP. 1829C GREATER THAN 2 MEANS 1 / (V(STPPAR) - 1) TIMES THE CAUCHY 1830C STEP. 1831C 1832C------------------------------- NOTES ------------------------------- 1833C 1834C *** ALGORITHM NOTES *** 1835C 1836C LET G AND H BE THE CURRENT GRADIENT AND HESSIAN APPROXIMA- 1837C TION RESPECTIVELY AND LET D BE THE CURRENT SCALE VECTOR. THIS 1838C ROUTINE ASSUMES DIG = DIAG(D)**-2 * G AND NWTSTP = H**-1 * G. 1839C THE STEP COMPUTED IS THE SAME ONE WOULD GET BY REPLACING G AND H 1840C BY DIAG(D)**-1 * G AND DIAG(D)**-1 * H * DIAG(D)**-1, 1841C COMPUTING STEP, AND TRANSLATING STEP BACK TO THE ORIGINAL 1842C VARIABLES, I.E., PREMULTIPLYING IT BY DIAG(D)**-1. 1843C 1844C *** REFERENCES *** 1845C 1846C 1. DENNIS, J.E., AND MEI, H.H.W. (1979), TWO NEW UNCONSTRAINED OPTI- 1847C MIZATION ALGORITHMS WHICH USE FUNCTION AND GRADIENT 1848C VALUES, J. OPTIM. THEORY APPLIC. 28, PP. 453-482. 1849C 2. POWELL, M.J.D. (1970), A HYBRID METHOD FOR NON-LINEAR EQUATIONS, 1850C IN NUMERICAL METHODS FOR NON-LINEAR EQUATIONS, EDITED BY 1851C P. RABINOWITZ, GORDON AND BREACH, LONDON. 1852C 1853C *** GENERAL *** 1854C 1855C CODED BY DAVID M. GAY. 1856C THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH SUPPORTED 1857C BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS MCS-7600324 AND 1858C MCS-7906671. 1859C 1860C------------------------ EXTERNAL QUANTITIES ------------------------ 1861C 1862C *** FUNCTIONS AND SUBROUTINES CALLED *** 1863C 1864 DOUBLE PRECISION DDOT 1865C 1866C 1867C *** INTRINSIC FUNCTIONS *** 1868C/+ 1869 DOUBLE PRECISION DSQRT 1870C/ 1871C-------------------------- LOCAL VARIABLES -------------------------- 1872C 1873 INTEGER I 1874 DOUBLE PRECISION CFACT, CNORM, CTRNWT, GHINVG, FEMNSQ, GNORM, 1875 1 NWTNRM, RELAX, RLAMBD, T, T1, T2 1876 DOUBLE PRECISION HALF, ONE, TWO, ZERO 1877C 1878C *** V SUBSCRIPTS *** 1879C 1880 INTEGER BIAS, DGNORM, DSTNRM, DST0, GRDFAC, GTHG, GTSTEP, 1881 1 NREDUC, NWTFAC, PREDUC, RADIUS, STPPAR 1882C 1883C *** DATA INITIALIZATIONS *** 1884C 1885C/6 1886C DATA HALF/0.5D+0/, ONE/1.D+0/, TWO/2.D+0/, ZERO/0.D+0/ 1887C/7 1888 PARAMETER (HALF=0.5D+0, ONE=1.D+0, TWO=2.D+0, ZERO=0.D+0) 1889C/ 1890C 1891C/6 1892C DATA BIAS/43/, DGNORM/1/, DSTNRM/2/, DST0/3/, GRDFAC/45/, 1893C 1 GTHG/44/, GTSTEP/4/, NREDUC/6/, NWTFAC/46/, PREDUC/7/, 1894C 2 RADIUS/8/, STPPAR/5/ 1895C/7 1896 PARAMETER (BIAS=43, DGNORM=1, DSTNRM=2, DST0=3, GRDFAC=45, 1897 1 GTHG=44, GTSTEP=4, NREDUC=6, NWTFAC=46, PREDUC=7, 1898 2 RADIUS=8, STPPAR=5) 1899C/ 1900C 1901C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ 1902C 1903 NWTNRM = V(DST0) 1904 RLAMBD = ONE 1905 IF (NWTNRM .GT. ZERO) RLAMBD = V(RADIUS) / NWTNRM 1906 GNORM = V(DGNORM) 1907 DO 10 I = 1, N 1908 10 STEP(I) = G(I) / GNORM 1909 GHINVG = DDOT(N, STEP,1,NWTSTP,1) 1910 V(NREDUC) = HALF * GHINVG * GNORM 1911 V(GRDFAC) = ZERO 1912 V(NWTFAC) = ZERO 1913 IF (RLAMBD .LT. ONE) GO TO 30 1914C 1915C *** THE NEWTON STEP IS INSIDE THE TRUST REGION *** 1916C 1917 V(STPPAR) = ZERO 1918 V(DSTNRM) = NWTNRM 1919 V(GTSTEP) = -GHINVG * GNORM 1920 V(PREDUC) = V(NREDUC) 1921 V(NWTFAC) = -ONE 1922 DO 20 I = 1, N 1923 20 STEP(I) = -NWTSTP(I) 1924 GO TO 999 1925C 1926 30 V(DSTNRM) = V(RADIUS) 1927 CFACT = (GNORM / V(GTHG))**2 1928C *** CAUCHY STEP = -CFACT * G. 1929 CNORM = GNORM * CFACT 1930 RELAX = ONE - V(BIAS) * (ONE - CNORM/GHINVG) 1931 IF (RLAMBD .LT. RELAX) GO TO 50 1932C 1933C *** STEP IS BETWEEN RELAXED NEWTON AND FULL NEWTON STEPS *** 1934C 1935 V(STPPAR) = ONE - (RLAMBD - RELAX) / (ONE - RELAX) 1936 T = -RLAMBD 1937 V(GTSTEP) = T * GHINVG * GNORM 1938 V(PREDUC) = RLAMBD * (ONE - HALF*RLAMBD) * GHINVG * GNORM 1939 V(NWTFAC) = T 1940 DO 40 I = 1, N 1941 40 STEP(I) = T * NWTSTP(I) 1942 GO TO 999 1943C 1944 50 IF (CNORM .LT. V(RADIUS)) GO TO 70 1945C 1946C *** THE CAUCHY STEP LIES OUTSIDE THE TRUST REGION -- 1947C *** STEP = SCALED CAUCHY STEP *** 1948C 1949 T = -V(RADIUS) / GNORM 1950 V(GRDFAC) = T 1951 V(STPPAR) = ONE + CNORM / V(RADIUS) 1952 V(GTSTEP) = -V(RADIUS) * GNORM 1953 V(PREDUC) = V(RADIUS)*(GNORM - HALF*V(RADIUS)*(V(GTHG)/GNORM)**2) 1954 DO 60 I = 1, N 1955 60 STEP(I) = T * DIG(I) 1956 GO TO 999 1957C 1958C *** COMPUTE DOGLEG STEP BETWEEN CAUCHY AND RELAXED NEWTON *** 1959C *** FEMUR = RELAXED NEWTON STEP MINUS CAUCHY STEP *** 1960C 1961 70 CTRNWT = CFACT * RELAX * GHINVG / GNORM 1962C *** CTRNWT = INNER PROD. OF CAUCHY AND RELAXED NEWTON STEPS, 1963C *** SCALED BY GNORM**-2. 1964 T1 = CTRNWT - CFACT**2 1965C *** T1 = INNER PROD. OF FEMUR AND CAUCHY STEP, SCALED BY 1966C *** GNORM**-2. 1967 T2 = (V(RADIUS)/GNORM)**2 - CFACT**2 1968 FEMNSQ = (RELAX*NWTNRM/GNORM)**2 - CTRNWT - T1 1969C *** FEMNSQ = SQUARE OF 2-NORM OF FEMUR, SCALED BY GNORM**-2. 1970 T = T2 / (T1 + DSQRT(T1**2 + FEMNSQ*T2)) 1971C *** DOGLEG STEP = CAUCHY STEP + T * FEMUR. 1972 T1 = (T - ONE) * CFACT 1973 V(GRDFAC) = T1 1974 T2 = -T * RELAX 1975 V(NWTFAC) = T2 1976 V(STPPAR) = TWO - T 1977 V(GTSTEP) = GNORM * (T1*GNORM + T2*GHINVG) 1978 V(PREDUC) = -(T1*GNORM) * ((T2 + ONE)*GNORM) 1979 1 - (T2*GNORM) * (ONE + HALF*T2)*GHINVG 1980 2 - HALF * (V(GTHG)*T1)**2 1981 DO 80 I = 1, N 1982 80 STEP(I) = T1*DIG(I) + T2*NWTSTP(I) 1983C 1984 999 RETURN 1985C *** LAST CARD OF DDBDOG FOLLOWS *** 1986 END 1987 SUBROUTINE DITSUM(D, G, IV, LIV, LV, P, V, X) 1988 use cfuncs 1989 save 1990 1991C 1992C *** PRINT ITERATION SUMMARY FOR ***SOL (VERSION 2.3) *** 1993C 1994C *** PARAMETER DECLARATIONS *** 1995C 1996 INTEGER LIV, LV, P 1997 INTEGER IV(LIV) 1998 DOUBLE PRECISION D(P), G(P), V(LV), X(P) 1999C 2000C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2001C 2002C *** LOCAL VARIABLES *** 2003C 2004 INTEGER ALG, I, IV1, M, NF, NG, OL, PU 2005C/6 2006C REAL MODEL1(6), MODEL2(6) 2007C/7 2008 CHARACTER*4 MODEL1(6), MODEL2(6) 2009C/ 2010 DOUBLE PRECISION NRELDF, OLDF, PRELDF, RELDF, ZERO 2011C 2012C *** INTRINSIC FUNCTIONS *** 2013C/+ 2014 INTEGER IABS 2015 DOUBLE PRECISION DABS, DMAX1 2016C/ 2017C *** NO EXTERNAL FUNCTIONS OR SUBROUTINES *** 2018C 2019C *** SUBSCRIPTS FOR IV AND V *** 2020C 2021 INTEGER ALGSAV, DSTNRM, F, FDIF, F0, NEEDHD, NFCALL, NFCOV, NGCOV, 2022 1 NGCALL, NITER, NREDUC, OUTLEV, PREDUC, PRNTIT, PRUNIT, 2023 2 RELDX, SOLPRT, STATPR, STPPAR, SUSED, X0PRT 2024C 2025C *** IV SUBSCRIPT VALUES *** 2026C 2027C/6 2028C DATA ALGSAV/51/, NEEDHD/36/, NFCALL/6/, NFCOV/52/, NGCALL/30/, 2029C 1 NGCOV/53/, NITER/31/, OUTLEV/19/, PRNTIT/39/, PRUNIT/21/, 2030C 2 SOLPRT/22/, STATPR/23/, SUSED/64/, X0PRT/24/ 2031C/7 2032 PARAMETER (ALGSAV=51, NEEDHD=36, NFCALL=6, NFCOV=52, NGCALL=30, 2033 1 NGCOV=53, NITER=31, OUTLEV=19, PRNTIT=39, PRUNIT=21, 2034 2 SOLPRT=22, STATPR=23, SUSED=64, X0PRT=24) 2035C/ 2036C 2037C *** V SUBSCRIPT VALUES *** 2038C 2039C/6 2040C DATA DSTNRM/2/, F/10/, F0/13/, FDIF/11/, NREDUC/6/, PREDUC/7/, 2041C 1 RELDX/17/, STPPAR/5/ 2042C/7 2043 PARAMETER (DSTNRM=2, F=10, F0=13, FDIF=11, NREDUC=6, PREDUC=7, 2044 1 RELDX=17, STPPAR=5) 2045C/ 2046C 2047C/6 2048C DATA ZERO/0.D+0/ 2049C/7 2050 PARAMETER (ZERO=0.D+0) 2051C/ 2052C/6 2053C DATA MODEL1(1)/4H /, MODEL1(2)/4H /, MODEL1(3)/4H /, 2054C 1 MODEL1(4)/4H /, MODEL1(5)/4H G /, MODEL1(6)/4H S /, 2055C 2 MODEL2(1)/4H G /, MODEL2(2)/4H S /, MODEL2(3)/4HG-S /, 2056C 3 MODEL2(4)/4HS-G /, MODEL2(5)/4H-S-G/, MODEL2(6)/4H-G-S/ 2057C/7 2058 DATA MODEL1/' ',' ',' ',' ',' G ',' S '/, 2059 1 MODEL2/' G ',' S ','G-S ','S-G ','-S-G','-G-S'/ 2060C/ 2061C 2062C------------------------------- BODY -------------------------------- 2063C 2064 PU = IV(PRUNIT) 2065 IF (PU .EQ. 0) GO TO 999 2066 IV1 = IV(1) 2067 IF (IV1 .GT. 62) IV1 = IV1 - 51 2068 OL = IV(OUTLEV) 2069 ALG = IV(ALGSAV) 2070 IF (IV1 .LT. 2 .OR. IV1 .GT. 15) GO TO 370 2071 IF (OL .EQ. 0) GO TO 120 2072 IF (IV1 .GE. 12) GO TO 120 2073 IF (IV1 .EQ. 2 .AND. IV(NITER) .EQ. 0) GO TO 390 2074 IF (IV1 .GE. 10 .AND. IV(PRNTIT) .EQ. 0) GO TO 120 2075 IF (IV1 .GT. 2) GO TO 10 2076 IV(PRNTIT) = IV(PRNTIT) + 1 2077 IF (IV(PRNTIT) .LT. IABS(OL)) GO TO 999 2078 10 NF = IV(NFCALL) - IABS(IV(NFCOV)) 2079 IV(PRNTIT) = 0 2080 RELDF = ZERO 2081 PRELDF = ZERO 2082 OLDF = DMAX1(DABS(V(F0)), DABS(V(F))) 2083 IF (OLDF .LE. ZERO) GO TO 20 2084 RELDF = V(FDIF) / OLDF 2085 PRELDF = V(PREDUC) / OLDF 2086 20 IF (OL .GT. 0) GO TO 60 2087C 2088C *** PRINT SHORT SUMMARY LINE *** 2089C 2090 IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 1) call h30() 2091 IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 2) call h40() 2092 IV(NEEDHD) = 0 2093 IF (ALG .EQ. 2) GO TO 50 2094 M = IV(SUSED) 2095 call h100s(IV(NITER), NF, V(F), RELDF, PRELDF, V(RELDX), 2096 1 MODEL1(M), MODEL2(M), V(STPPAR)) 2097 GO TO 120 2098C 2099 50 call h110s(IV(NITER), NF, V(F), RELDF, PRELDF, V(RELDX), 2100 1 V(STPPAR)) 2101 GO TO 120 2102C 2103C *** PRINT LONG SUMMARY LINE *** 2104C 2105 60 IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 1) call h70() 2106 IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 2) call h80() 2107 IV(NEEDHD) = 0 2108 NRELDF = ZERO 2109 IF (OLDF .GT. ZERO) NRELDF = V(NREDUC) / OLDF 2110 IF (ALG .EQ. 2) GO TO 90 2111 M = IV(SUSED) 2112 call h100l(IV(NITER), NF, V(F), RELDF, PRELDF, V(RELDX), 2113 1 MODEL1(M), MODEL2(M), V(STPPAR), V(DSTNRM), NRELDF) 2114 GO TO 120 2115C 2116 90 call h110l(IV(NITER), NF, V(F), RELDF, PRELDF, 2117 1 V(RELDX), V(STPPAR), V(DSTNRM), NRELDF) 2118C 2119 120 IF (IV(STATPR) .LT. 0) GO TO 430 2120 GO TO (999, 999, 130, 150, 170, 190, 210, 230, 250, 270, 290, 310, 2121 1 330, 350, 520), IV1 2122C 2123 130 call cnlprt(' ***** X-CONVERGENCE *****', 26) 2124 GO TO 430 2125C 2126 150 call cnlprt(' ***** RELATIVE FUNCTION CONVERGENCE *****', 42) 2127 GO TO 430 2128C 2129 170 call cnlprt 2130 1(' ***** X- AND RELATIVE FUNCTION CONVERGENCE *****', 49) 2131 GO TO 430 2132C 2133 190 call cnlprt(' ***** ABSOLUTE FUNCTION CONVERGENCE *****', 42) 2134 GO TO 430 2135C 2136 210 call cnlprt(' ***** SINGULAR CONVERGENCE *****', 33) 2137 GO TO 430 2138C 2139 230 call cnlprt(' ***** FALSE CONVERGENCE *****', 30) 2140 GO TO 430 2141C 2142 250 call cnlprt(' ***** FUNCTION EVALUATION LIMIT *****', 38) 2143 GO TO 430 2144C 2145 270 call cnlprt(' ***** ITERATION LIMIT *****', 28) 2146 GO TO 430 2147C 2148 290 call cnlprt(' ***** STOPX *****', 18) 2149 GO TO 430 2150C 2151 310 call cnlprt(' ***** INITIAL F(X) CANNOT BE COMPUTED *****', 44) 2152 GO TO 390 2153C 2154 330 call cnlprt(' ***** BAD PARAMETERS TO ASSESS *****', 37) 2155 GO TO 999 2156C 2157 350 call cnlprt(' ***** GRADIENT COULD NOT BE COMPUTED *****', 43) 2158 IF (IV(NITER) .GT. 0) GO TO 480 2159 GO TO 390 2160C 2161 370 call h380(IV(1)) 2162 GO TO 999 2163C 2164C *** INITIAL CALL ON DITSUM *** 2165C 2166 390 call h400(P, X, D) 2167 IF (IV1 .GE. 12) GO TO 999 2168 IV(NEEDHD) = 0 2169 IV(PRNTIT) = 0 2170 IF (OL .EQ. 0) GO TO 999 2171 IF (OL .LT. 0 .AND. ALG .EQ. 1) call h30() 2172 IF (OL .LT. 0 .AND. ALG .EQ. 2) call h40() 2173 IF (OL .GT. 0 .AND. ALG .EQ. 1) call h70() 2174 IF (OL .GT. 0 .AND. ALG .EQ. 2) call h80() 2175 IF (ALG .EQ. 1) call h410(V(F)) 2176 IF (ALG .EQ. 2) call h420(V(F)) 2177 GO TO 999 2178C 2179C *** PRINT VARIOUS INFORMATION REQUESTED ON SOLUTION *** 2180C 2181 430 IV(NEEDHD) = 1 2182 IF (IV(STATPR) .EQ. 0) GO TO 480 2183 OLDF = DMAX1(DABS(V(F0)), DABS(V(F))) 2184 PRELDF = ZERO 2185 NRELDF = ZERO 2186 IF (OLDF .LE. ZERO) GO TO 440 2187 PRELDF = V(PREDUC) / OLDF 2188 NRELDF = V(NREDUC) / OLDF 2189 440 NF = IV(NFCALL) - IV(NFCOV) 2190 NG = IV(NGCALL) - IV(NGCOV) 2191 call h450(V(F), V(RELDX), NF, NG, PRELDF, NRELDF) 2192C 2193 IF (IV(NFCOV) .GT. 0) call h460(IV(NFCOV)) 2194 IF (IV(NGCOV) .GT. 0) call h470(IV(NGCOV)) 2195C 2196 480 IF (IV(SOLPRT) .EQ. 0) GO TO 999 2197 IV(NEEDHD) = 1 2198 call cnlprt(' I FINAL X(I) D(I) G(I)', 2199 1 48) 2200 call h500(P, X, D, G) 2201 GO TO 999 2202C 2203 520 call cnlprt(' INCONSISTENT DIMENSIONS', 24) 2204 999 RETURN 2205C *** LAST CARD OF DITSUM FOLLOWS *** 2206 END 2207 SUBROUTINE DLITVM(N, X, L, Y) 2208 save 2209C 2210C *** SOLVE (L**T)*X = Y, WHERE L IS AN N X N LOWER TRIANGULAR 2211C *** MATRIX STORED COMPACTLY BY ROWS. X AND Y MAY OCCUPY THE SAME 2212C *** STORAGE. *** 2213C 2214 INTEGER N 2215 DOUBLE PRECISION X(N), L, Y(N) 2216 DIMENSION L(N*(N+1)/2) 2217 INTEGER I, II, IJ, IM1, I0, J, NP1 2218 DOUBLE PRECISION XI, ZERO 2219C/6 2220C DATA ZERO/0.D+0/ 2221C/7 2222 PARAMETER (ZERO=0.D+0) 2223C/ 2224C 2225 DO 10 I = 1, N 2226 10 X(I) = Y(I) 2227 NP1 = N + 1 2228 I0 = N*(N+1)/2 2229 DO 30 II = 1, N 2230 I = NP1 - II 2231 XI = X(I)/L(I0) 2232 X(I) = XI 2233 IF (I .LE. 1) GO TO 999 2234 I0 = I0 - I 2235 IF (XI .EQ. ZERO) GO TO 30 2236 IM1 = I - 1 2237 DO 20 J = 1, IM1 2238 IJ = I0 + J 2239 X(J) = X(J) - XI*L(IJ) 2240 20 CONTINUE 2241 30 CONTINUE 2242 999 RETURN 2243C *** LAST CARD OF DLITVM FOLLOWS *** 2244 END 2245 SUBROUTINE DLIVMU(N, X, L, Y) 2246 save 2247C 2248C *** SOLVE L*X = Y, WHERE L IS AN N X N LOWER TRIANGULAR 2249C *** MATRIX STORED COMPACTLY BY ROWS. X AND Y MAY OCCUPY THE SAME 2250C *** STORAGE. *** 2251C 2252 INTEGER N 2253 DOUBLE PRECISION X(N), L, Y(N) 2254 DIMENSION L(N*(N+1)/2) 2255 DOUBLE PRECISION DDOT 2256 INTEGER I, J, K 2257 DOUBLE PRECISION T, ZERO 2258C/6 2259C DATA ZERO/0.D+0/ 2260C/7 2261 PARAMETER (ZERO=0.D+0) 2262C/ 2263C 2264 DO 10 K = 1, N 2265 IF (Y(K) .NE. ZERO) GO TO 20 2266 X(K) = ZERO 2267 10 CONTINUE 2268 GO TO 999 2269 20 J = K*(K+1)/2 2270 X(K) = Y(K) / L(J) 2271 IF (K .GE. N) GO TO 999 2272 K = K + 1 2273 DO 30 I = K, N 2274 T = DDOT(I-1, L(J+1),1,X,1) 2275 J = J + I 2276 X(I) = (Y(I) - T)/L(J) 2277 30 CONTINUE 2278 999 RETURN 2279C *** LAST CARD OF DLIVMU FOLLOWS *** 2280 END 2281 SUBROUTINE DLTVMU(N, X, L, Y) 2282 save 2283C 2284C *** COMPUTE X = (L**T)*Y, WHERE L IS AN N X N LOWER 2285C *** TRIANGULAR MATRIX STORED COMPACTLY BY ROWS. X AND Y MAY 2286C *** OCCUPY THE SAME STORAGE. *** 2287C 2288 INTEGER N 2289 DOUBLE PRECISION X(N), L, Y(N) 2290 DIMENSION L(N*(N+1)/2) 2291 INTEGER I, IJ, I0, J 2292 DOUBLE PRECISION YI, ZERO 2293C/6 2294C DATA ZERO/0.D+0/ 2295C/7 2296 PARAMETER (ZERO=0.D+0) 2297C/ 2298C 2299 I0 = 0 2300 DO 20 I = 1, N 2301 YI = Y(I) 2302 X(I) = ZERO 2303 DO 10 J = 1, I 2304 IJ = I0 + J 2305 X(J) = X(J) + YI*L(IJ) 2306 10 CONTINUE 2307 I0 = I0 + I 2308 20 CONTINUE 2309 999 RETURN 2310C *** LAST CARD OF DLTVMU FOLLOWS *** 2311 END 2312 SUBROUTINE DLUPDT(BETA, GAMMA, L, LAMBDA, LPLUS, N, W, Z) 2313 save 2314C 2315C *** COMPUTE LPLUS = SECANT UPDATE OF L *** 2316C 2317C *** PARAMETER DECLARATIONS *** 2318C 2319 INTEGER N 2320 DOUBLE PRECISION BETA(N), GAMMA(N), L, LAMBDA(N), LPLUS, 2321 1 W(N), Z(N) 2322 DIMENSION L(N*(N+1)/2), LPLUS(N*(N+1)/2) 2323C 2324C-------------------------- PARAMETER USAGE -------------------------- 2325C 2326C BETA = SCRATCH VECTOR. 2327C GAMMA = SCRATCH VECTOR. 2328C L (INPUT) LOWER TRIANGULAR MATRIX, STORED ROWWISE. 2329C LAMBDA = SCRATCH VECTOR. 2330C LPLUS (OUTPUT) LOWER TRIANGULAR MATRIX, STORED ROWWISE, WHICH MAY 2331C OCCUPY THE SAME STORAGE AS L. 2332C N (INPUT) LENGTH OF VECTOR PARAMETERS AND ORDER OF MATRICES. 2333C W (INPUT, DESTROYED ON OUTPUT) RIGHT SINGULAR VECTOR OF RANK 1 2334C CORRECTION TO L. 2335C Z (INPUT, DESTROYED ON OUTPUT) LEFT SINGULAR VECTOR OF RANK 1 2336C CORRECTION TO L. 2337C 2338C------------------------------- NOTES ------------------------------- 2339C 2340C *** APPLICATION AND USAGE RESTRICTIONS *** 2341C 2342C THIS ROUTINE UPDATES THE CHOLESKY FACTOR L OF A SYMMETRIC 2343C POSITIVE DEFINITE MATRIX TO WHICH A SECANT UPDATE IS BEING 2344C APPLIED -- IT COMPUTES A CHOLESKY FACTOR LPLUS OF 2345C L * (I + Z*W**T) * (I + W*Z**T) * L**T. IT IS ASSUMED THAT W 2346C AND Z HAVE BEEN CHOSEN SO THAT THE UPDATED MATRIX IS STRICTLY 2347C POSITIVE DEFINITE. 2348C 2349C *** ALGORITHM NOTES *** 2350C 2351C THIS CODE USES RECURRENCE 3 OF REF. 1 (WITH D(J) = 1 FOR ALL J) 2352C TO COMPUTE LPLUS OF THE FORM L * (I + Z*W**T) * Q, WHERE Q 2353C IS AN ORTHOGONAL MATRIX THAT MAKES THE RESULT LOWER TRIANGULAR. 2354C LPLUS MAY HAVE SOME NEGATIVE DIAGONAL ELEMENTS. 2355C 2356C *** REFERENCES *** 2357C 2358C 1. GOLDFARB, D. (1976), FACTORIZED VARIABLE METRIC METHODS FOR UNCON- 2359C STRAINED OPTIMIZATION, MATH. COMPUT. 30, PP. 796-811. 2360C 2361C *** GENERAL *** 2362C 2363C CODED BY DAVID M. GAY (FALL 1979). 2364C THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH SUPPORTED 2365C BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS MCS-7600324 AND 2366C MCS-7906671. 2367C 2368C------------------------ EXTERNAL QUANTITIES ------------------------ 2369C 2370C *** INTRINSIC FUNCTIONS *** 2371C/+ 2372 DOUBLE PRECISION DSQRT 2373C/ 2374C-------------------------- LOCAL VARIABLES -------------------------- 2375C 2376 INTEGER I, IJ, J, JJ, JP1, K, NM1, NP1 2377 DOUBLE PRECISION A, B, BJ, ETA, GJ, LJ, LIJ, LJJ, NU, S, THETA, 2378 1 WJ, ZJ 2379 DOUBLE PRECISION ONE, ZERO 2380C 2381C *** DATA INITIALIZATIONS *** 2382C 2383C/6 2384C DATA ONE/1.D+0/, ZERO/0.D+0/ 2385C/7 2386 PARAMETER (ONE=1.D+0, ZERO=0.D+0) 2387C/ 2388C 2389C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ 2390C 2391 NU = ONE 2392 ETA = ZERO 2393 IF (N .LE. 1) GO TO 30 2394 NM1 = N - 1 2395C 2396C *** TEMPORARILY STORE S(J) = SUM OVER K = J+1 TO N OF W(K)**2 IN 2397C *** LAMBDA(J). 2398C 2399 S = ZERO 2400 DO 10 I = 1, NM1 2401 J = N - I 2402 S = S + W(J+1)**2 2403 LAMBDA(J) = S 2404 10 CONTINUE 2405C 2406C *** COMPUTE LAMBDA, GAMMA, AND BETA BY GOLDFARB*S RECURRENCE 3. 2407C 2408 DO 20 J = 1, NM1 2409 WJ = W(J) 2410 A = NU*Z(J) - ETA*WJ 2411 THETA = ONE + A*WJ 2412 S = A*LAMBDA(J) 2413 LJ = DSQRT(THETA**2 + A*S) 2414 IF (THETA .GT. ZERO) LJ = -LJ 2415 LAMBDA(J) = LJ 2416 B = THETA*WJ + S 2417 GAMMA(J) = B * NU / LJ 2418 BETA(J) = (A - B*ETA) / LJ 2419 NU = -NU / LJ 2420 ETA = -(ETA + (A**2)/(THETA - LJ)) / LJ 2421 20 CONTINUE 2422 30 LAMBDA(N) = ONE + (NU*Z(N) - ETA*W(N))*W(N) 2423C 2424C *** UPDATE L, GRADUALLY OVERWRITING W AND Z WITH L*W AND L*Z. 2425C 2426 NP1 = N + 1 2427 JJ = N * (N + 1) / 2 2428 DO 60 K = 1, N 2429 J = NP1 - K 2430 LJ = LAMBDA(J) 2431 LJJ = L(JJ) 2432 LPLUS(JJ) = LJ * LJJ 2433 WJ = W(J) 2434 W(J) = LJJ * WJ 2435 ZJ = Z(J) 2436 Z(J) = LJJ * ZJ 2437 IF (K .EQ. 1) GO TO 50 2438 BJ = BETA(J) 2439 GJ = GAMMA(J) 2440 IJ = JJ + J 2441 JP1 = J + 1 2442 DO 40 I = JP1, N 2443 LIJ = L(IJ) 2444 LPLUS(IJ) = LJ*LIJ + BJ*W(I) + GJ*Z(I) 2445 W(I) = W(I) + LIJ*WJ 2446 Z(I) = Z(I) + LIJ*ZJ 2447 IJ = IJ + I 2448 40 CONTINUE 2449 50 JJ = JJ - J 2450 60 CONTINUE 2451C 2452 999 RETURN 2453C *** LAST CARD OF DLUPDT FOLLOWS *** 2454 END 2455 SUBROUTINE DLVMUL(N, X, L, Y) 2456 save 2457C 2458C *** COMPUTE X = L*Y, WHERE L IS AN N X N LOWER TRIANGULAR 2459C *** MATRIX STORED COMPACTLY BY ROWS. X AND Y MAY OCCUPY THE SAME 2460C *** STORAGE. *** 2461C 2462 INTEGER N 2463 DOUBLE PRECISION X(N), L, Y(N) 2464 DIMENSION L(N*(N+1)/2) 2465 INTEGER I, II, IJ, I0, J, NP1 2466 DOUBLE PRECISION T, ZERO 2467C/6 2468C DATA ZERO/0.D+0/ 2469C/7 2470 PARAMETER (ZERO=0.D+0) 2471C/ 2472C 2473 NP1 = N + 1 2474 I0 = N*(N+1)/2 2475 DO 20 II = 1, N 2476 I = NP1 - II 2477 I0 = I0 - I 2478 T = ZERO 2479 DO 10 J = 1, I 2480 IJ = I0 + J 2481 T = T + L(IJ)*Y(J) 2482 10 CONTINUE 2483 X(I) = T 2484 20 CONTINUE 2485 999 RETURN 2486C *** LAST CARD OF DLVMUL FOLLOWS *** 2487 END 2488 SUBROUTINE DPARCK(ALG, D, IV, LIV, LV, N, V) 2489 save 2490C 2491C *** CHECK ***SOL (VERSION 2.3) PARAMETERS, PRINT CHANGED VALUES *** 2492C 2493C *** ALG = 1 FOR REGRESSION, ALG = 2 FOR GENERAL UNCONSTRAINED OPT. 2494C 2495 INTEGER ALG, LIV, LV, N 2496 INTEGER IV(LIV) 2497 DOUBLE PRECISION D(N), V(LV) 2498C 2499 EXTERNAL DVDFLT 2500 DOUBLE PRECISION D1MACH 2501C DVDFLT -- SUPPLIES DEFAULT PARAMETER VALUES TO V ALONE. 2502C/+ 2503 INTEGER MAX0 2504C/ 2505C 2506C *** LOCAL VARIABLES *** 2507C 2508 INTEGER I, II, IV1, J, K, L, M, MIV1, MIV2, NDFALT, PARSV1, PU 2509 INTEGER IJMP, JLIM(2), MINIV(2), NDFLT(2) 2510C/6 2511C INTEGER VARNM(2), SH(2) 2512C REAL CNGD(3), DFLT(3), VN(2,34), WHICH(3) 2513C/7 2514 CHARACTER*1 VARNM(2), SH(2) 2515 CHARACTER*4 CNGD(3), DFLT(3), VN(2,34), WHICH(3) 2516C/ 2517 DOUBLE PRECISION BIG, MACHEP, TINY, VK, VM(34), VX(34), ZERO 2518C 2519C *** IV AND V SUBSCRIPTS *** 2520C 2521 INTEGER ALGSAV, DINIT, DTYPE, DTYPE0, EPSLON, INITS, IVNEED, 2522 1 LASTIV, LASTV, LMAT, NEXTIV, NEXTV, NVDFLT, OLDN, 2523 2 PARPRT, PARSAV, PERM, PRUNIT, VNEED 2524C 2525C 2526C/6 2527C DATA ALGSAV/51/, DINIT/38/, DTYPE/16/, DTYPE0/54/, EPSLON/19/, 2528C 1 INITS/25/, IVNEED/3/, LASTIV/44/, LASTV/45/, LMAT/42/, 2529C 2 NEXTIV/46/, NEXTV/47/, NVDFLT/50/, OLDN/38/, PARPRT/20/, 2530C 3 PARSAV/49/, PERM/58/, PRUNIT/21/, VNEED/4/ 2531C/7 2532 PARAMETER (ALGSAV=51, DINIT=38, DTYPE=16, DTYPE0=54, EPSLON=19, 2533 1 INITS=25, IVNEED=3, LASTIV=44, LASTV=45, LMAT=42, 2534 2 NEXTIV=46, NEXTV=47, NVDFLT=50, OLDN=38, PARPRT=20, 2535 3 PARSAV=49, PERM=58, PRUNIT=21, VNEED=4) 2536C SAVE BIG, MACHEP, TINY 2537C/ 2538C 2539 DATA BIG/0.D+0/, MACHEP/-1.D+0/, TINY/1.D+0/, ZERO/0.D+0/ 2540C/6 2541C DATA VN(1,1),VN(2,1)/4HEPSL,4HON../ 2542C DATA VN(1,2),VN(2,2)/4HPHMN,4HFC../ 2543C DATA VN(1,3),VN(2,3)/4HPHMX,4HFC../ 2544C DATA VN(1,4),VN(2,4)/4HDECF,4HAC../ 2545C DATA VN(1,5),VN(2,5)/4HINCF,4HAC../ 2546C DATA VN(1,6),VN(2,6)/4HRDFC,4HMN../ 2547C DATA VN(1,7),VN(2,7)/4HRDFC,4HMX../ 2548C DATA VN(1,8),VN(2,8)/4HTUNE,4HR1../ 2549C DATA VN(1,9),VN(2,9)/4HTUNE,4HR2../ 2550C DATA VN(1,10),VN(2,10)/4HTUNE,4HR3../ 2551C DATA VN(1,11),VN(2,11)/4HTUNE,4HR4../ 2552C DATA VN(1,12),VN(2,12)/4HTUNE,4HR5../ 2553C DATA VN(1,13),VN(2,13)/4HAFCT,4HOL../ 2554C DATA VN(1,14),VN(2,14)/4HRFCT,4HOL../ 2555C DATA VN(1,15),VN(2,15)/4HXCTO,4HL.../ 2556C DATA VN(1,16),VN(2,16)/4HXFTO,4HL.../ 2557C DATA VN(1,17),VN(2,17)/4HLMAX,4H0.../ 2558C DATA VN(1,18),VN(2,18)/4HLMAX,4HS.../ 2559C DATA VN(1,19),VN(2,19)/4HSCTO,4HL.../ 2560C DATA VN(1,20),VN(2,20)/4HDINI,4HT.../ 2561C DATA VN(1,21),VN(2,21)/4HDTIN,4HIT../ 2562C DATA VN(1,22),VN(2,22)/4HD0IN,4HIT../ 2563C DATA VN(1,23),VN(2,23)/4HDFAC,4H..../ 2564C DATA VN(1,24),VN(2,24)/4HDLTF,4HDC../ 2565C DATA VN(1,25),VN(2,25)/4HDLTF,4HDJ../ 2566C DATA VN(1,26),VN(2,26)/4HDELT,4HA0../ 2567C DATA VN(1,27),VN(2,27)/4HFUZZ,4H..../ 2568C DATA VN(1,28),VN(2,28)/4HRLIM,4HIT../ 2569C DATA VN(1,29),VN(2,29)/4HCOSM,4HIN../ 2570C DATA VN(1,30),VN(2,30)/4HHUBE,4HRC../ 2571C DATA VN(1,31),VN(2,31)/4HRSPT,4HOL../ 2572C DATA VN(1,32),VN(2,32)/4HSIGM,4HIN../ 2573C DATA VN(1,33),VN(2,33)/4HETA0,4H..../ 2574C DATA VN(1,34),VN(2,34)/4HBIAS,4H..../ 2575C/7 2576 DATA VN(1,1),VN(2,1)/'EPSL','ON..'/ 2577 DATA VN(1,2),VN(2,2)/'PHMN','FC..'/ 2578 DATA VN(1,3),VN(2,3)/'PHMX','FC..'/ 2579 DATA VN(1,4),VN(2,4)/'DECF','AC..'/ 2580 DATA VN(1,5),VN(2,5)/'INCF','AC..'/ 2581 DATA VN(1,6),VN(2,6)/'RDFC','MN..'/ 2582 DATA VN(1,7),VN(2,7)/'RDFC','MX..'/ 2583 DATA VN(1,8),VN(2,8)/'TUNE','R1..'/ 2584 DATA VN(1,9),VN(2,9)/'TUNE','R2..'/ 2585 DATA VN(1,10),VN(2,10)/'TUNE','R3..'/ 2586 DATA VN(1,11),VN(2,11)/'TUNE','R4..'/ 2587 DATA VN(1,12),VN(2,12)/'TUNE','R5..'/ 2588 DATA VN(1,13),VN(2,13)/'AFCT','OL..'/ 2589 DATA VN(1,14),VN(2,14)/'RFCT','OL..'/ 2590 DATA VN(1,15),VN(2,15)/'XCTO','L...'/ 2591 DATA VN(1,16),VN(2,16)/'XFTO','L...'/ 2592 DATA VN(1,17),VN(2,17)/'LMAX','0...'/ 2593 DATA VN(1,18),VN(2,18)/'LMAX','S...'/ 2594 DATA VN(1,19),VN(2,19)/'SCTO','L...'/ 2595 DATA VN(1,20),VN(2,20)/'DINI','T...'/ 2596 DATA VN(1,21),VN(2,21)/'DTIN','IT..'/ 2597 DATA VN(1,22),VN(2,22)/'D0IN','IT..'/ 2598 DATA VN(1,23),VN(2,23)/'DFAC','....'/ 2599 DATA VN(1,24),VN(2,24)/'DLTF','DC..'/ 2600 DATA VN(1,25),VN(2,25)/'DLTF','DJ..'/ 2601 DATA VN(1,26),VN(2,26)/'DELT','A0..'/ 2602 DATA VN(1,27),VN(2,27)/'FUZZ','....'/ 2603 DATA VN(1,28),VN(2,28)/'RLIM','IT..'/ 2604 DATA VN(1,29),VN(2,29)/'COSM','IN..'/ 2605 DATA VN(1,30),VN(2,30)/'HUBE','RC..'/ 2606 DATA VN(1,31),VN(2,31)/'RSPT','OL..'/ 2607 DATA VN(1,32),VN(2,32)/'SIGM','IN..'/ 2608 DATA VN(1,33),VN(2,33)/'ETA0','....'/ 2609 DATA VN(1,34),VN(2,34)/'BIAS','....'/ 2610C/ 2611C 2612 DATA VM(1)/1.0D-3/, VM(2)/-0.99D+0/, VM(3)/1.0D-3/, VM(4)/1.0D-2/, 2613 1 VM(5)/1.2D+0/, VM(6)/1.D-2/, VM(7)/1.2D+0/, VM(8)/0.D+0/, 2614 2 VM(9)/0.D+0/, VM(10)/1.D-3/, VM(11)/-1.D+0/, VM(15)/0.D+0/, 2615 3 VM(16)/0.D+0/, VM(19)/0.D+0/, VM(20)/-10.D+0/, VM(21)/0.D+0/, 2616 4 VM(22)/0.D+0/, VM(23)/0.D+0/, VM(27)/1.01D+0/, 2617 5 VM(28)/1.D+10/, VM(30)/0.D+0/, VM(31)/0.D+0/, VM(32)/0.D+0/, 2618 6 VM(34)/0.D+0/ 2619 DATA VX(1)/0.9D+0/, VX(2)/-1.D-3/, VX(3)/1.D+1/, VX(4)/0.8D+0/, 2620 1 VX(5)/1.D+2/, VX(6)/0.8D+0/, VX(7)/1.D+2/, VX(8)/0.5D+0/, 2621 2 VX(9)/0.5D+0/, VX(10)/1.D+0/, VX(11)/1.D+0/, VX(14)/0.1D+0/, 2622 3 VX(15)/1.D+0/, VX(16)/1.D+0/, VX(19)/1.D+0/, VX(23)/1.D+0/, 2623 4 VX(24)/1.D+0/, VX(25)/1.D+0/, VX(26)/1.D+0/, VX(27)/1.D+10/, 2624 5 VX(29)/1.D+0/, VX(31)/1.D+0/, VX(32)/1.D+0/, VX(33)/1.D+0/, 2625 6 VX(34)/1.D+0/ 2626C 2627C/6 2628C DATA VARNM(1)/1HP/, VARNM(2)/1HN/, SH(1)/1HS/, SH(2)/1HH/ 2629C DATA CNGD(1),CNGD(2),CNGD(3)/4H---C,4HHANG,4HED V/, 2630C 1 DFLT(1),DFLT(2),DFLT(3)/4HNOND,4HEFAU,4HLT V/ 2631C/7 2632 DATA VARNM(1)/'P'/, VARNM(2)/'N'/, SH(1)/'S'/, SH(2)/'H'/ 2633 DATA CNGD(1),CNGD(2),CNGD(3)/'---C','HANG','ED V'/, 2634 1 DFLT(1),DFLT(2),DFLT(3)/'NOND','EFAU','LT V'/ 2635C/ 2636 DATA IJMP/33/, JLIM(1)/0/, JLIM(2)/24/, NDFLT(1)/32/, NDFLT(2)/25/ 2637 DATA MINIV(1)/80/, MINIV(2)/59/ 2638C 2639C............................... BODY ................................ 2640C 2641 IF (ALG .LT. 1 .OR. ALG .GT. 2) GO TO 330 2642 IF (IV(1) .EQ. 0) CALL DDEFLT(ALG, IV, LIV, LV, V) 2643 PU = IV(PRUNIT) 2644 MIV1 = MINIV(ALG) 2645 IF (PERM .LE. LIV) MIV1 = MAX0(MIV1, IV(PERM) - 1) 2646 IF (IVNEED .LE. LIV) MIV2 = MIV1 + MAX0(IV(IVNEED), 0) 2647 IF (LASTIV .LE. LIV) IV(LASTIV) = MIV2 2648 IF (LIV .LT. MIV1) GO TO 290 2649 IV(IVNEED) = 0 2650 IV(LASTV) = MAX0(IV(VNEED), 0) + IV(LMAT) - 1 2651 IF (LIV .LT. MIV2) GO TO 290 2652 IF (LV .LT. IV(LASTV)) GO TO 310 2653 IV(VNEED) = 0 2654 IF (ALG .EQ. IV(ALGSAV)) GO TO 20 2655c IF (PU .NE. 0) WRITE(PU,10) ALG, IV(ALGSAV) 2656c 10 FORMAT(/' THE FIRST PARAMETER TO DDEFLT SHOULD BE',I3, 2657c 1 12H RATHER THAN,I3) 2658 IV(1) = 82 2659 GO TO 999 2660 20 IV1 = IV(1) 2661 IF (IV1 .LT. 12 .OR. IV1 .GT. 14) GO TO 50 2662 IF (N .GE. 1) GO TO 40 2663 IV(1) = 81 2664 IF (PU .EQ. 0) GO TO 999 2665c WRITE(PU,30) VARNM(ALG), N 2666c 30 FORMAT(/8H /// BAD,A1,2H =,I5) 2667 GO TO 999 2668 40 IF (IV1 .NE. 14) IV(NEXTIV) = IV(PERM) 2669 IF (IV1 .NE. 14) IV(NEXTV) = IV(LMAT) 2670 IF (IV1 .EQ. 13) GO TO 999 2671 K = IV(PARSAV) - EPSLON 2672 CALL DVDFLT(ALG, LV-K, V(K+1)) 2673 IV(DTYPE0) = 2 - ALG 2674 IV(OLDN) = N 2675 WHICH(1) = DFLT(1) 2676 WHICH(2) = DFLT(2) 2677 WHICH(3) = DFLT(3) 2678 GO TO 100 2679 50 IF (N .EQ. IV(OLDN)) GO TO 70 2680 IV(1) = 17 2681 IF (PU .EQ. 0) GO TO 999 2682c WRITE(PU,60) VARNM(ALG), IV(OLDN), N 2683c 60 FORMAT(/5H /// ,1A1,14H CHANGED FROM ,I5,4H TO ,I5) 2684 GO TO 999 2685C 2686 70 IF (IV1 .LE. 11 .AND. IV1 .GE. 1) GO TO 90 2687 IV(1) = 80 2688c IF (PU .NE. 0) WRITE(PU,80) IV1 2689c 80 FORMAT(/13H /// IV(1) =,I5,28H SHOULD BE BETWEEN 0 AND 14.) 2690 GO TO 999 2691C 2692 90 WHICH(1) = CNGD(1) 2693 WHICH(2) = CNGD(2) 2694 WHICH(3) = CNGD(3) 2695C 2696 100 IF (IV1 .EQ. 14) IV1 = 12 2697 IF (BIG .GT. TINY) GO TO 110 2698 TINY = D1MACH(1) 2699 MACHEP = D1MACH(4) 2700 BIG = D1MACH(2) 2701 VM(12) = MACHEP 2702 VX(12) = BIG 2703 VM(13) = TINY 2704 VX(13) = BIG 2705 VM(14) = MACHEP 2706 VM(17) = TINY 2707 VX(17) = BIG 2708 VM(18) = TINY 2709 VX(18) = BIG 2710 VX(20) = BIG 2711 VX(21) = BIG 2712 VX(22) = BIG 2713 VM(24) = MACHEP 2714 VM(25) = MACHEP 2715 VM(26) = MACHEP 2716 VX(28) = DSQRT(D1MACH(2))*16. 2717 VM(29) = MACHEP 2718 VX(30) = BIG 2719 VM(33) = MACHEP 2720 110 M = 0 2721 I = 1 2722 J = JLIM(ALG) 2723 K = EPSLON 2724 NDFALT = NDFLT(ALG) 2725 DO 140 L = 1, NDFALT 2726 VK = V(K) 2727 IF (VK .GE. VM(I) .AND. VK .LE. VX(I)) GO TO 130 2728 M = K 2729c IF (PU .NE. 0) WRITE(PU,120) VN(1,I), VN(2,I), K, VK, 2730c 1 VM(I), VX(I) 2731c 120 FORMAT(/6H /// ,2A4,5H.. V(,I2,3H) =,D11.3,7H SHOULD, 2732c 1 11H BE BETWEEN,D11.3,4H AND,D11.3) 2733 130 K = K + 1 2734 I = I + 1 2735 IF (I .EQ. J) I = IJMP 2736 140 CONTINUE 2737C 2738 IF (IV(NVDFLT) .EQ. NDFALT) GO TO 160 2739 IV(1) = 51 2740 IF (PU .EQ. 0) GO TO 999 2741c WRITE(PU,150) IV(NVDFLT), NDFALT 2742c 150 FORMAT(/13H IV(NVDFLT) =,I5,13H RATHER THAN ,I5) 2743 GO TO 999 2744 160 IF ((IV(DTYPE) .GT. 0 .OR. V(DINIT) .GT. ZERO) .AND. IV1 .EQ. 12) 2745 1 GO TO 190 2746 DO 180 I = 1, N 2747 IF (D(I) .GT. ZERO) GO TO 180 2748 M = 18 2749c IF (PU .NE. 0) WRITE(PU,170) I, D(I) 2750c 170 FORMAT(/8H /// D(,I3,3H) =,D11.3,19H SHOULD BE POSITIVE) 2751 180 CONTINUE 2752 190 IF (M .EQ. 0) GO TO 200 2753 IV(1) = M 2754 GO TO 999 2755C 2756 200 IF (PU .EQ. 0 .OR. IV(PARPRT) .EQ. 0) GO TO 999 2757 IF (IV1 .NE. 12 .OR. IV(INITS) .EQ. ALG-1) GO TO 220 2758 M = 1 2759c WRITE(PU,210) SH(ALG), IV(INITS) 2760c 210 FORMAT(/22H NONDEFAULT VALUES..../5H INIT,A1,14H..... IV(25) =, 2761c 1 I3) 2762 220 IF (IV(DTYPE) .EQ. IV(DTYPE0)) GO TO 240 2763c IF (M .EQ. 0) WRITE(PU,250) WHICH 2764 M = 1 2765c WRITE(PU,230) IV(DTYPE) 2766c 230 FORMAT(20H DTYPE..... IV(16) =,I3) 2767 240 I = 1 2768 J = JLIM(ALG) 2769 K = EPSLON 2770 L = IV(PARSAV) 2771 NDFALT = NDFLT(ALG) 2772 DO 280 II = 1, NDFALT 2773 IF (V(K) .EQ. V(L)) GO TO 270 2774c IF (M .EQ. 0) WRITE(PU,250) WHICH 2775c 250 FORMAT(/1H ,3A4,9HALUES..../) 2776 M = 1 2777c WRITE(PU,260) VN(1,I), VN(2,I), K, V(K) 2778c 260 FORMAT(1X,2A4,5H.. V(,I2,3H) =,D15.7) 2779 270 K = K + 1 2780 L = L + 1 2781 I = I + 1 2782 IF (I .EQ. J) I = IJMP 2783 280 CONTINUE 2784C 2785 IV(DTYPE0) = IV(DTYPE) 2786 PARSV1 = IV(PARSAV) 2787 CALL DCOPY(IV(NVDFLT), V(EPSLON),1,V(PARSV1),1) 2788 GO TO 999 2789C 2790 290 IV(1) = 15 2791 IF (PU .EQ. 0) GO TO 999 2792c WRITE(PU,300) LIV, MIV2 2793c 300 FORMAT(/10H /// LIV =,I5,17H MUST BE AT LEAST,I5) 2794 IF (LIV .LT. MIV1) GO TO 999 2795 IF (LV .LT. IV(LASTV)) GO TO 310 2796 GO TO 999 2797C 2798 310 IV(1) = 16 2799 IF (PU .EQ. 0) GO TO 999 2800c WRITE(PU,320) LV, IV(LASTV) 2801c 320 FORMAT(/9H /// LV =,I5,17H MUST BE AT LEAST,I5) 2802 GO TO 999 2803C 2804 330 IV(1) = 67 2805C 2806 999 RETURN 2807C *** LAST CARD OF DPARCK FOLLOWS *** 2808 END 2809 DOUBLE PRECISION FUNCTION DRELST(P, D, X, X0) 2810 save 2811C 2812C *** COMPUTE AND RETURN RELATIVE DIFFERENCE BETWEEN X AND X0 *** 2813C *** NL2SOL VERSION 2.2 *** 2814C 2815 INTEGER P 2816 DOUBLE PRECISION D(P), X(P), X0(P) 2817C/+ 2818 DOUBLE PRECISION DABS 2819C/ 2820 INTEGER I 2821 DOUBLE PRECISION EMAX, T, XMAX, ZERO 2822C/6 2823C DATA ZERO/0.D+0/ 2824C/7 2825 PARAMETER (ZERO=0.D+0) 2826C/ 2827C 2828 EMAX = ZERO 2829 XMAX = ZERO 2830 DO 10 I = 1, P 2831 T = DABS(D(I) * (X(I) - X0(I))) 2832 IF (EMAX .LT. T) EMAX = T 2833 T = D(I) * (DABS(X(I)) + DABS(X0(I))) 2834 IF (XMAX .LT. T) XMAX = T 2835 10 CONTINUE 2836 DRELST = ZERO 2837 IF (XMAX .GT. ZERO) DRELST = EMAX / XMAX 2838 999 RETURN 2839C *** LAST CARD OF DRELST FOLLOWS *** 2840 END 2841 LOGICAL FUNCTION DSTOPX(IDUMMY) 2842 save 2843C *****PARAMETERS... 2844 INTEGER IDUMMY 2845C 2846C .................................................................. 2847C 2848C *****PURPOSE... 2849C THIS FUNCTION MAY SERVE AS THE DSTOPX (ASYNCHRONOUS INTERRUPTION) 2850C FUNCTION FOR THE NL2SOL (NONLINEAR LEAST-SQUARES) PACKAGE AT 2851C THOSE INSTALLATIONS WHICH DO NOT WISH TO IMPLEMENT A 2852C DYNAMIC DSTOPX. 2853C 2854C *****ALGORITHM NOTES... 2855C AT INSTALLATIONS WHERE THE NL2SOL SYSTEM IS USED 2856C INTERACTIVELY, THIS DUMMY DSTOPX SHOULD BE REPLACED BY A 2857C FUNCTION THAT RETURNS .TRUE. IF AND ONLY IF THE INTERRUPT 2858C (BREAK) KEY HAS BEEN PRESSED SINCE THE LAST CALL ON DSTOPX. 2859C 2860C .................................................................. 2861C 2862 DSTOPX = .FALSE. 2863 RETURN 2864 END 2865 SUBROUTINE DSMSNO(N, D, X, CALCF, IV, LIV, LV, V, 2866 1 UIPARM, URPARM, UFPARM) 2867 save 2868C 2869C *** MINIMIZE GENERAL UNCONSTRAINED OBJECTIVE FUNCTION USING 2870C *** FINITE-DIFFERENCE GRADIENTS AND SECANT HESSIAN APPROXIMATIONS. 2871C 2872 INTEGER N, LIV, LV 2873 INTEGER IV(LIV), UIPARM(*) 2874 DOUBLE PRECISION D(N), X(N), V(LV), URPARM(*) 2875C DIMENSION V(77 + N*(N+17)/2), UIPARM(*), URPARM(*) 2876 EXTERNAL CALCF, UFPARM 2877C 2878C *** PURPOSE *** 2879C 2880C THIS ROUTINE INTERACTS WITH SUBROUTINE DSNOIT IN AN ATTEMPT 2881C TO FIND AN N-VECTOR X* THAT MINIMIZES THE (UNCONSTRAINED) 2882C OBJECTIVE FUNCTION COMPUTED BY CALCF. (OFTEN THE X* FOUND IS 2883C A LOCAL MINIMIZER RATHER THAN A GLOBAL ONE.) 2884C 2885C *** PARAMETERS *** 2886C 2887C THE PARAMETERS FOR DSMSNO ARE THE SAME AS THOSE FOR DSUMSL 2888C (WHICH SEE), EXCEPT THAT CALCG IS OMITTED. INSTEAD OF CALLING 2889C CALCG TO OBTAIN THE GRADIENT OF THE OBJECTIVE FUNCTION AT X, 2890C DSMSNO CALLS DSGRD2, WHICH COMPUTES AN APPROXIMATION TO THE 2891C GRADIENT BY FINITE (FORWARD AND CENTRAL) DIFFERENCES USING THE 2892C METHOD OF REF. 1. THE FOLLOWING INPUT COMPONENT IS OF INTEREST 2893C IN THIS REGARD (AND IS NOT DESCRIBED IN DSUMSL). 2894C 2895C V(ETA0)..... V(42) IS AN ESTIMATED BOUND ON THE RELATIVE ERROR IN THE 2896C OBJECTIVE FUNCTION VALUE COMPUTED BY CALCF... 2897C (TRUE VALUE) = (COMPUTED VALUE) * (1 + E), 2898C WHERE ABS(E) .LE. V(ETA0). DEFAULT = MACHEP * 10**3, 2899C WHERE MACHEP IS THE UNIT ROUNDOFF. 2900C 2901C THE OUTPUT VALUES IV(NFCALL) AND IV(NGCALL) HAVE DIFFERENT 2902C MEANINGS FOR DSMSNO THAN FOR DSUMSL... 2903C 2904C IV(NFCALL)... IV(6) IS THE NUMBER OF CALLS SO FAR MADE ON CALCF (I.E., 2905C FUNCTION EVALUATIONS) EXCLUDING THOSE MADE ONLY FOR 2906C COMPUTING GRADIENTS. THE INPUT VALUE IV(MXFCAL) IS A 2907C LIMIT ON IV(NFCALL). 2908C IV(NGCALL)... IV(30) IS THE NUMBER OF FUNCTION EVALUATIONS MADE ONLY 2909C FOR COMPUTING GRADIENTS. THE TOTAL NUMBER OF FUNCTION 2910C EVALUATIONS IS THUS IV(NFCALL) + IV(NGCALL). 2911C 2912C *** REFERENCES *** 2913C 2914C 1. STEWART, G.W. (1967), A MODIFICATION OF DAVIDON*S MINIMIZATION 2915C METHOD TO ACCEPT DIFFERENCE APPROXIMATIONS OF DERIVATIVES, 2916C J. ASSOC. COMPUT. MACH. 14, PP. 72-83. 2917C. 2918C *** GENERAL *** 2919C 2920C CODED BY DAVID M. GAY (WINTER 1980). REVISED SEPT. 1982. 2921C THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH 2922C SUPPORTED IN PART BY THE NATIONAL SCIENCE FOUNDATION UNDER 2923C GRANTS MCS-7600324, DCR75-10143, 76-14311DSS, MCS76-11989, 2924C AND MCS-7906671. 2925C 2926C 2927C---------------------------- DECLARATIONS --------------------------- 2928C 2929 EXTERNAL DSNOIT 2930C 2931C DSNOIT.... OVERSEES COMPUTATION OF FINITE-DIFFERENCE GRADIENT AND 2932C CALLS DSUMIT TO CARRY OUT DSUMSL ALGORITHM. 2933C 2934 INTEGER NF 2935 DOUBLE PRECISION FX 2936C 2937C *** SUBSCRIPTS FOR IV *** 2938C 2939 INTEGER NFCALL, TOOBIG 2940C 2941C/6 2942C DATA NFCALL/6/, TOOBIG/2/ 2943C/7 2944 PARAMETER (NFCALL=6, TOOBIG=2) 2945C/ 2946C 2947C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ 2948C 2949 10 CALL DSNOIT(D, FX, IV, LIV, LV, N, V, X) 2950 IF (IV(1) .GT. 2) GO TO 999 2951C 2952C *** COMPUTE FUNCTION *** 2953C 2954 NF = IV(NFCALL) 2955 CALL CALCF(N, X, NF, FX, UIPARM, URPARM, UFPARM) 2956 IF (NF .LE. 0) IV(TOOBIG) = 1 2957 GO TO 10 2958C 2959C 2960 999 RETURN 2961C *** LAST CARD OF DSMSNO FOLLOWS *** 2962 END 2963 SUBROUTINE DSNOIT(D, FX, IV, LIV, LV, N, V, X) 2964 save 2965C 2966C *** ITERATION DRIVER FOR DSMSNO... 2967C *** MINIMIZE GENERAL UNCONSTRAINED OBJECTIVE FUNCTION USING 2968C *** FINITE-DIFFERENCE GRADIENTS AND SECANT HESSIAN APPROXIMATIONS. 2969C 2970 INTEGER LIV, LV, N 2971 INTEGER IV(LIV) 2972 DOUBLE PRECISION D(N), FX, X(N), V(LV) 2973C DIMENSION V(77 + N*(N+17)/2) 2974C 2975C *** PURPOSE *** 2976C 2977C THIS ROUTINE INTERACTS WITH SUBROUTINE DSUMIT IN AN ATTEMPT 2978C TO FIND AN N-VECTOR X* THAT MINIMIZES THE (UNCONSTRAINED) 2979C OBJECTIVE FUNCTION FX = F(X) COMPUTED BY THE CALLER. (OFTEN 2980C THE X* FOUND IS A LOCAL MINIMIZER RATHER THAN A GLOBAL ONE.) 2981C 2982C *** PARAMETERS *** 2983C 2984C THE PARAMETERS FOR DSNOIT ARE THE SAME AS THOSE FOR DSUMSL 2985C (WHICH SEE), EXCEPT THAT CALCF, CALCG, UIPARM, URPARM, AND UFPARM 2986C ARE OMITTED, AND A PARAMETER FX FOR THE OBJECTIVE FUNCTION 2987C VALUE AT X IS ADDED. INSTEAD OF CALLING CALCG TO OBTAIN THE 2988C GRADIENT OF THE OBJECTIVE FUNCTION AT X, DSNOIT CALLS DSGRD2, 2989C WHICH COMPUTES AN APPROXIMATION TO THE GRADIENT BY FINITE 2990C (FORWARD AND CENTRAL) DIFFERENCES USING THE METHOD OF REF. 1. 2991C THE FOLLOWING INPUT COMPONENT IS OF INTEREST IN THIS REGARD 2992C (AND IS NOT DESCRIBED IN DSUMSL). 2993C 2994C V(ETA0)..... V(42) IS AN ESTIMATED BOUND ON THE RELATIVE ERROR IN THE 2995C OBJECTIVE FUNCTION VALUE COMPUTED BY CALCF... 2996C (TRUE VALUE) = (COMPUTED VALUE) * (1 + E), 2997C WHERE ABS(E) .LE. V(ETA0). DEFAULT = MACHEP * 10**3, 2998C WHERE MACHEP IS THE UNIT ROUNDOFF. 2999C 3000C THE OUTPUT VALUES IV(NFCALL) AND IV(NGCALL) HAVE DIFFERENT 3001C MEANINGS FOR DSMSNO THAN FOR DSUMSL... 3002C 3003C IV(NFCALL)... IV(6) IS THE NUMBER OF CALLS SO FAR MADE ON CALCF (I.E., 3004C FUNCTION EVALUATIONS) EXCLUDING THOSE MADE ONLY FOR 3005C COMPUTING GRADIENTS. THE INPUT VALUE IV(MXFCAL) IS A 3006C LIMIT ON IV(NFCALL). 3007C IV(NGCALL)... IV(30) IS THE NUMBER OF FUNCTION EVALUATIONS MADE ONLY 3008C FOR COMPUTING GRADIENTS. THE TOTAL NUMBER OF FUNCTION 3009C EVALUATIONS IS THUS IV(NFCALL) + IV(NGCALL). 3010C 3011C *** REFERENCES *** 3012C 3013C 1. STEWART, G.W. (1967), A MODIFICATION OF DAVIDON*S MINIMIZATION 3014C METHOD TO ACCEPT DIFFERENCE APPROXIMATIONS OF DERIVATIVES, 3015C J. ASSOC. COMPUT. MACH. 14, PP. 72-83. 3016C. 3017C *** GENERAL *** 3018C 3019C CODED BY DAVID M. GAY (AUGUST 1982). 3020C 3021C---------------------------- DECLARATIONS --------------------------- 3022C 3023 EXTERNAL DDEFLT, DSGRD2, DSUMIT, DVSCPY 3024 DOUBLE PRECISION DDOT 3025C 3026C DDEFLT.... SUPPLIES DEFAULT PARAMETER VALUES. 3027C DSGRD2... COMPUTES FINITE-DIFFERENCE GRADIENT APPROXIMATION. 3028C DSUMIT.... REVERSE-COMMUNICATION ROUTINE THAT DOES DSUMSL ALGORITHM. 3029C DVSCPY... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR. 3030C 3031 INTEGER ALPHA, G1, I, IV1, J, K, W 3032 DOUBLE PRECISION ZERO 3033C 3034C *** SUBSCRIPTS FOR IV *** 3035C 3036 INTEGER ETA0, F, G, LMAT, NEXTV, NFGCAL, NGCALL, 3037 1 NITER, SGIRC, TOOBIG, VNEED 3038C 3039C/6 3040C DATA ETA0/42/, F/10/, G/28/, LMAT/42/, NEXTV/47/, 3041C 1 NFGCAL/7/, NGCALL/30/, NITER/31/, SGIRC/57/, 3042C 2 TOOBIG/2/, VNEED/4/ 3043C/7 3044 PARAMETER (DTYPE=16, ETA0=42, F=10, G=28, LMAT=42, NEXTV=47, 3045 1 NFCALL=6, NFGCAL=7, NGCALL=30, NITER=31, SGIRC=57, 3046 2 TOOBIG=2, VNEED=4) 3047C/ 3048C/6 3049C DATA ZERO/0.D+0/ 3050C/7 3051 PARAMETER (ONE=1.D+0, ZERO=0.D+0) 3052C/ 3053C 3054C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ 3055C 3056 IV1 = IV(1) 3057 IF (IV1 .EQ. 1) GO TO 10 3058 IF (IV1 .EQ. 2) GO TO 50 3059 IF (IV(1) .EQ. 0) CALL DDEFLT(2, IV, LIV, LV, V) 3060 IV(VNEED) = IV(VNEED) + 2*N + 6 3061 IV1 = IV(1) 3062 IF (IV1 .EQ. 14) GO TO 10 3063 IF (IV1 .GT. 2 .AND. IV1 .LT. 12) GO TO 10 3064 G1 = 1 3065 IF (IV1 .EQ. 12) IV(1) = 13 3066 GO TO 20 3067C 3068 10 G1 = IV(G) 3069C 3070 20 CALL DSUMIT(D, FX, V(G1), IV, LIV, LV, N, V, X) 3071c IF (IV(1) - 2) 999, 30, 70 3072 IF (IV(1) .LT. 2) GO TO 999 3073 IF (IV(1) .GT. 2) GO TO 70 3074C 3075C *** COMPUTE GRADIENT *** 3076C 3077 30 IF (IV(NITER) .EQ. 0) CALL DVSCPY(N, V(G1), ZERO) 3078 J = IV(LMAT) 3079 K = G1 - N 3080 DO 40 I = 1, N 3081 V(K) = DDOT(I, V(J),1,V(J),1) 3082 K = K + 1 3083 J = J + I 3084 40 CONTINUE 3085C *** UNDO INCREMENT OF IV(NGCALL) DONE BY DSUMIT *** 3086 IV(NGCALL) = IV(NGCALL) - 1 3087C *** STORE RETURN CODE FROM DSGRD2 IN IV(SGIRC) *** 3088 IV(SGIRC) = 0 3089C *** X MAY HAVE BEEN RESTORED, SO COPY BACK FX... *** 3090 FX = V(F) 3091 GO TO 60 3092C 3093C *** GRADIENT LOOP *** 3094C 3095 50 IF (IV(TOOBIG) .EQ. 0) GO TO 60 3096 IV(NFGCAL) = 0 3097 GO TO 10 3098C 3099 60 G1 = IV(G) 3100 ALPHA = G1 - N 3101 W = ALPHA - 6 3102 CALL DSGRD2(V(ALPHA), D, V(ETA0), FX, V(G1), IV(SGIRC), N, V(W),X) 3103 IF (IV(SGIRC) .EQ. 0) GO TO 10 3104 IV(NGCALL) = IV(NGCALL) + 1 3105 GO TO 999 3106C 3107 70 IF (IV(1) .NE. 14) GO TO 999 3108C 3109C *** STORAGE ALLOCATION *** 3110C 3111 IV(G) = IV(NEXTV) + N + 6 3112 IV(NEXTV) = IV(G) + N 3113 IF (IV1 .NE. 13) GO TO 10 3114C 3115 999 RETURN 3116C *** LAST CARD OF DSNOIT FOLLOWS *** 3117 END 3118 SUBROUTINE DSGRD2 (ALPHA, D, ETA0, FX, G, IRC, N, W, X) 3119 save 3120C 3121C *** COMPUTE FINITE DIFFERENCE GRADIENT BY STWEART*S SCHEME *** 3122C 3123C *** PARAMETERS *** 3124C 3125 INTEGER IRC, N 3126 DOUBLE PRECISION ALPHA(N), D(N), ETA0, FX, G(N), W(6), X(N) 3127C 3128C....................................................................... 3129C 3130C *** PURPOSE *** 3131C 3132C THIS SUBROUTINE USES AN EMBELLISHED FORM OF THE FINITE-DIFFER- 3133C ENCE SCHEME PROPOSED BY STEWART (REF. 1) TO APPROXIMATE THE 3134C GRADIENT OF THE FUNCTION F(X), WHOSE VALUES ARE SUPPLIED BY 3135C REVERSE COMMUNICATION. 3136C 3137C *** PARAMETER DESCRIPTION *** 3138C 3139C ALPHA IN (APPROXIMATE) DIAGONAL ELEMENTS OF THE HESSIAN OF F(X). 3140C D IN SCALE VECTOR SUCH THAT D(I)*X(I), I = 1,...,N, ARE IN 3141C COMPARABLE UNITS. 3142C ETA0 IN ESTIMATED BOUND ON RELATIVE ERROR IN THE FUNCTION VALUE... 3143C (TRUE VALUE) = (COMPUTED VALUE)*(1+E), WHERE 3144C ABS(E) .LE. ETA0. 3145C FX I/O ON INPUT, FX MUST BE THE COMPUTED VALUE OF F(X). ON 3146C OUTPUT WITH IRC = 0, FX HAS BEEN RESTORED TO ITS ORIGINAL 3147C VALUE, THE ONE IT HAD WHEN DSGRD2 WAS LAST CALLED WITH 3148C IRC = 0. 3149C G I/O ON INPUT WITH IRC = 0, G SHOULD CONTAIN AN APPROXIMATION 3150C TO THE GRADIENT OF F NEAR X, E.G., THE GRADIENT AT THE 3151C PREVIOUS ITERATE. WHEN DSGRD2 RETURNS WITH IRC = 0, G IS 3152C THE DESIRED FINITE-DIFFERENCE APPROXIMATION TO THE 3153C GRADIENT AT X. 3154C IRC I/O INPUT/RETURN CODE... BEFORE THE VERY FIRST CALL ON DSGRD2, 3155C THE CALLER MUST SET IRC TO 0. WHENEVER DSGRD2 RETURNS A 3156C NONZERO VALUE FOR IRC, IT HAS PERTURBED SOME COMPONENT OF 3157C X... THE CALLER SHOULD EVALUATE F(X) AND CALL DSGRD2 3158C AGAIN WITH FX = F(X). 3159C N IN THE NUMBER OF VARIABLES (COMPONENTS OF X) ON WHICH F 3160C DEPENDS. 3161C X I/O ON INPUT WITH IRC = 0, X IS THE POINT AT WHICH THE 3162C GRADIENT OF F IS DESIRED. ON OUTPUT WITH IRC NONZERO, X 3163C IS THE POINT AT WHICH F SHOULD BE EVALUATED. ON OUTPUT 3164C WITH IRC = 0, X HAS BEEN RESTORED TO ITS ORIGINAL VALUE 3165C (THE ONE IT HAD WHEN DSGRD2 WAS LAST CALLED WITH IRC = 0) 3166C AND G CONTAINS THE DESIRED GRADIENT APPROXIMATION. 3167C W I/O WORK VECTOR OF LENGTH 6 IN WHICH DSGRD2 SAVES CERTAIN 3168C QUANTITIES WHILE THE CALLER IS EVALUATING F(X) AT A 3169C PERTURBED X. 3170C 3171C *** APPLICATION AND USAGE RESTRICTIONS *** 3172C 3173C THIS ROUTINE IS INTENDED FOR USE WITH QUASI-NEWTON ROUTINES 3174C FOR UNCONSTRAINED MINIMIZATION (IN WHICH CASE ALPHA COMES FROM 3175C THE DIAGONAL OF THE QUASI-NEWTON HESSIAN APPROXIMATION). 3176C 3177C *** ALGORITHM NOTES *** 3178C 3179C THIS CODE DEPARTS FROM THE SCHEME PROPOSED BY STEWART (REF. 1) 3180C IN ITS GUARDING AGAINST OVERLY LARGE OR SMALL STEP SIZES AND ITS 3181C HANDLING OF SPECIAL CASES (SUCH AS ZERO COMPONENTS OF ALPHA OR G). 3182C 3183C *** REFERENCES *** 3184C 3185C 1. STEWART, G.W. (1967), A MODIFICATION OF DAVIDON*S MINIMIZATION 3186C METHOD TO ACCEPT DIFFERENCE APPROXIMATIONS OF DERIVATIVES, 3187C J. ASSOC. COMPUT. MACH. 14, PP. 72-83. 3188C 3189C *** HISTORY *** 3190C 3191C DESIGNED AND CODED BY DAVID M. GAY (SUMMER 1977/SUMMER 1980). 3192C 3193C *** GENERAL *** 3194C 3195C THIS ROUTINE WAS PREPARED IN CONNECTION WITH WORK SUPPORTED BY 3196C THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS MCS76-00324 AND 3197C MCS-7906671. 3198C 3199C....................................................................... 3200C 3201C ***** EXTERNAL FUNCTION ***** 3202C 3203 DOUBLE PRECISION D1MACH 3204C 3205C ***** INTRINSIC FUNCTIONS ***** 3206C/+ 3207 INTEGER IABS 3208 DOUBLE PRECISION DABS, DMAX1, DSQRT 3209C/ 3210C ***** LOCAL VARIABLES ***** 3211C 3212 INTEGER FH, FX0, HSAVE, I, XISAVE 3213 DOUBLE PRECISION AAI, AFX, AFXETA, AGI, ALPHAI, AXI, AXIBAR, 3214 1 DISCON, ETA, GI, H, HMIN 3215 DOUBLE PRECISION C2000, FOUR, HMAX0, HMIN0, H0, MACHEP, ONE, P002, 3216 1 THREE, TWO, ZERO 3217C 3218C/6 3219C DATA C2000/2.0D+3/, FOUR/4.0D+0/, HMAX0/0.02D+0/, HMIN0/5.0D+1/, 3220C 1 ONE/1.0D+0/, P002/0.002D+0/, THREE/3.0D+0/, 3221C 2 TWO/2.0D+0/, ZERO/0.0D+0/ 3222C/7 3223 PARAMETER (C2000=2.0D+3, FOUR=4.0D+0, HMAX0=0.02D+0, HMIN0=5.0D+1, 3224 1 ONE=1.0D+0, P002=0.002D+0, THREE=3.0D+0, 3225 2 TWO=2.0D+0, ZERO=0.0D+0) 3226C/ 3227C/6 3228C DATA FH/3/, FX0/4/, HSAVE/5/, XISAVE/6/ 3229C/7 3230 PARAMETER (FH=3, FX0=4, HSAVE=5, XISAVE=6) 3231C/ 3232C 3233C--------------------------------- BODY ------------------------------ 3234C 3235c IF (IRC) 140, 100, 210 3236 IF (IRC .LT. 0) GO TO 140 3237 IF (IRC .GT. 0) GO TO 210 3238C 3239C *** FRESH START -- GET MACHINE-DEPENDENT CONSTANTS *** 3240C 3241C STORE MACHEP IN W(1) AND H0 IN W(2), WHERE MACHEP IS THE UNIT 3242C ROUNDOFF (THE SMALLEST POSITIVE NUMBER SUCH THAT 3243C 1 + MACHEP .GT. 1 AND 1 - MACHEP .LT. 1), AND H0 IS THE 3244C SQUARE-ROOT OF MACHEP. 3245C 3246 100 W(1) = D1MACH(4) 3247 W(2) = DSQRT(W(1)) 3248C 3249 W(FX0) = FX 3250C 3251C *** INCREMENT I AND START COMPUTING G(I) *** 3252C 3253 110 I = IABS(IRC) + 1 3254 IF (I .GT. N) GO TO 300 3255 IRC = I 3256 AFX = DABS(W(FX0)) 3257 MACHEP = W(1) 3258 H0 = W(2) 3259 HMIN = HMIN0 * MACHEP 3260 W(XISAVE) = X(I) 3261 AXI = DABS(X(I)) 3262 AXIBAR = DMAX1(AXI, ONE/D(I)) 3263 GI = G(I) 3264 AGI = DABS(GI) 3265 ETA = DABS(ETA0) 3266 IF (AFX .GT. ZERO) ETA = DMAX1(ETA, AGI*AXI*MACHEP/AFX) 3267 ALPHAI = ALPHA(I) 3268 IF (ALPHAI .EQ. ZERO) GO TO 170 3269 IF (GI .EQ. ZERO .OR. FX .EQ. ZERO) GO TO 180 3270 AFXETA = AFX*ETA 3271 AAI = DABS(ALPHAI) 3272C 3273C *** COMPUTE H = STEWART*S FORWARD-DIFFERENCE STEP SIZE. 3274C 3275 IF (GI**2 .LE. AFXETA*AAI) GO TO 120 3276 H = TWO*DSQRT(AFXETA/AAI) 3277 H = H*(ONE - AAI*H/(THREE*AAI*H + FOUR*AGI)) 3278 GO TO 130 3279 120 H = TWO*(AFXETA*AGI/(AAI**2))**(ONE/THREE) 3280 H = H*(ONE - TWO*AGI/(THREE*AAI*H + FOUR*AGI)) 3281C 3282C *** ENSURE THAT H IS NOT INSIGNIFICANTLY SMALL *** 3283C 3284 130 H = DMAX1(H, HMIN*AXIBAR) 3285C 3286C *** USE FORWARD DIFFERENCE IF BOUND ON TRUNCATION ERROR IS AT 3287C *** MOST 10**-3. 3288C 3289 IF (AAI*H .LE. P002*AGI) GO TO 160 3290C 3291C *** COMPUTE H = STEWART*S STEP FOR CENTRAL DIFFERENCE. 3292C 3293 DISCON = C2000*AFXETA 3294 H = DISCON/(AGI + DSQRT(GI**2 + AAI*DISCON)) 3295C 3296C *** ENSURE THAT H IS NEITHER TOO SMALL NOR TOO BIG *** 3297C 3298 H = DMAX1(H, HMIN*AXIBAR) 3299 IF (H .GE. HMAX0*AXIBAR) H = AXIBAR * H0**(TWO/THREE) 3300C 3301C *** COMPUTE CENTRAL DIFFERENCE *** 3302C 3303 IRC = -I 3304 GO TO 200 3305C 3306 140 H = -W(HSAVE) 3307 I = IABS(IRC) 3308 IF (H .GT. ZERO) GO TO 150 3309 W(FH) = FX 3310 GO TO 200 3311C 3312 150 G(I) = (W(FH) - FX) / (TWO * H) 3313 X(I) = W(XISAVE) 3314 GO TO 110 3315C 3316C *** COMPUTE FORWARD DIFFERENCES IN VARIOUS CASES *** 3317C 3318 160 IF (H .GE. HMAX0*AXIBAR) H = H0 * AXIBAR 3319 IF (ALPHAI*GI .LT. ZERO) H = -H 3320 GO TO 200 3321 170 H = AXIBAR 3322 GO TO 200 3323 180 H = H0 * AXIBAR 3324C 3325 200 X(I) = W(XISAVE) + H 3326 W(HSAVE) = H 3327 GO TO 999 3328C 3329C *** COMPUTE ACTUAL FORWARD DIFFERENCE *** 3330C 3331 210 G(IRC) = (FX - W(FX0)) / W(HSAVE) 3332 X(IRC) = W(XISAVE) 3333 GO TO 110 3334C 3335C *** RESTORE FX AND INDICATE THAT G HAS BEEN COMPUTED *** 3336C 3337 300 FX = W(FX0) 3338 IRC = 0 3339C 3340 999 RETURN 3341C *** LAST CARD OF DSGRD2 FOLLOWS *** 3342 END 3343