1 SUBROUTINE SRN2G(D, DR, IV, LIV, LV, N, ND, N1, N2, P, R, 2 1 RD, V, X) 3c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. 4c ALL RIGHTS RESERVED. 5c Based on Government Sponsored Research NAS7-03001. 6c File: SRN2G.for Ten subrs used by the 7c David Gay & Linda Kaufman nonlinear LS package. 8c Needed for versions that do not allow Bounded variables. 9c SRN2G is called by SNLAFU, SNLAGU, & SRNSG. 10c 11C>> 2015-07-09 SRN2G Krogh Introduced TP to avoid divide by 0, 12C>> 2000-01-07 SRN2G Krogh Moved COV1 = IV(COMAT) up in SN2CVP. 13C>> 1998-10-29 SRN2G Krogh Moved external statement up for mangle. 14c>> 1996-07-09 SRN2G Krogh Changes for conversion to C. 15c>> 1995-01-26 SRN2G Krogh Moved formats up for C conversion. 16c>> 1994-11-02 SRN2G Krogh Changes to use M77CON 17c>> 1993-03-10 SRN2G CLL Moved stmt NN = ... to follow IF (IV1 ... 18c>> 1992-04-27 CLL Comment out unreferenced stmt labels. 19c>> 1992-04-13 CLL Change from Hollerith to '...' syntax in formats. 20c>> 1990-06-29 CLL Changes to formats in SN2CVP. 21c>> 1990-06-12 CLL Revised SRN2G & SG7LIT from DMG 4/19/90 22c>> 1990-03-30 CLL JPL 23c>> 1990-03-14 CLL JPL 24c>> 1990-06-12 CLL 25c>> 1990-04-23 CLL (Recent revision by DMG) 26*** from netlib, Thu Apr 19 11:58:57 EDT 1990 *** 27c--S replaces "?": ?RN2G,?C7VFN,?D7TPR,?D7UPD,?G7LIT,?ITSUM,?IVSET 28c--& ?L7VML,?N2CVP,?N2LRD,?Q7APL,?Q7RAD,?V2NRM,?V7CPY,?V7SCP, ?N2G 29c--& ?A7SST,?F7HES,?G7QTS,?L7MST,?L7SQR,?L7SRT,?L7SVN,?L7SVX,?L7TVM 30c--& ?PARCK,?R7MDC,?RLDST,?S7LUP,?S7LVM,?V2AXY,?L7ITV,?L7IVM,?O7PRD 31c--& ?L7NVR,?L7TSQ,?V7SCL,?N2RDP,?NLAFU,?NLAGU,?RNSG 32C 33C *** REVISED ITERATION DRIVER FOR NL2SOL (VERSION 2.3) *** 34C 35 INTEGER LIV, LV, N, ND, N1, N2, P 36 INTEGER IV(LIV) 37 REAL D(P), DR(ND,P), R(ND), RD(ND), V(LV), X(P) 38C 39C ------------------------- PARAMETER USAGE -------------------------- 40C 41C D........ SCALE VECTOR. 42C DR....... DERIVATIVES OF R AT X. 43C IV....... INTEGER VALUES ARRAY. 44C LIV...... LENGTH OF IV... LIV MUST BE AT LEAST P + 82. 45C LV....... LENGTH OF V... LV MUST BE AT LEAST 105 + P*(2*P+16). 46C N........ TOTAL NUMBER OF RESIDUALS. 47C ND....... MAX. NO. OF RESIDUALS PASSED ON ONE CALL. 48C N1....... LOWEST ROW INDEX FOR RESIDUALS SUPPLIED THIS TIME. 49C N2....... HIGHEST ROW INDEX FOR RESIDUALS SUPPLIED THIS TIME. 50C P........ NUMBER OF PARAMETERS (COMPONENTS OF X) BEING ESTIMATED. 51C R........ RESIDUALS. 52C RD....... RD(I) = SQRT(G(I)**T * H(I)**-1 * G(I)) ON OUTPUT WHEN 53C IV(RDREQ) IS NONZERO. SRN2G SETS IV(REGD) = 1 IF RD 54C IS SUCCESSFULLY COMPUTED, TO 0 IF NO ATTEMPT WAS MADE 55C TO COMPUTE IT, AND TO -1 IF H (THE FINITE-DIFFERENCE HESSIAN) 56C WAS INDEFINITE. IF ND .GE. N, THEN RD IS ALSO USED AS 57C TEMPORARY STORAGE. 58C V........ FLOATING-POINT VALUES ARRAY. 59C X........ PARAMETER VECTOR BEING ESTIMATED (INPUT = INITIAL GUESS, 60C OUTPUT = BEST VALUE FOUND). 61C 62C *** DISCUSSION *** 63C 64C NOTE... NL2SOL AND NL2ITR (MENTIONED BELOW) ARE DESCRIBED IN 65C ACM TRANS. MATH. SOFTWARE, VOL. 7, PP. 369-383 (AN ADAPTIVE 66C NONLINEAR LEAST-SQUARES ALGORITHM, BY J.E. DENNIS, D.M. GAY, 67C AND R.E. WELSCH). 68C 69C THIS ROUTINE CARRIES OUT ITERATIONS FOR SOLVING NONLINEAR 70C LEAST SQUARES PROBLEMS. WHEN ND = N, IT IS SIMILAR TO NL2ITR 71C (WITH J = DR), EXCEPT THAT R(X) AND DR(X) NEED NOT BE INITIALIZED 72C WHEN SRN2G IS CALLED WITH IV(1) = 0 OR 12. SRN2G ALSO ALLOWS 73C R AND DR TO BE SUPPLIED ROW-WISE -- JUST SET ND = 1 AND CALL 74C SRN2G ONCE FOR EACH ROW WHEN PROVIDING RESIDUALS AND JACOBIANS. 75C ANOTHER NEW FEATURE IS THAT CALLING SRN2G WITH IV(1) = 13 76C CAUSES STORAGE ALLOCATION ONLY TO BE PERFORMED -- ON RETURN, SUCH 77C COMPONENTS AS IV(G) (THE FIRST SUBSCRIPT IN G OF THE GRADIENT) 78C AND IV(S) (THE FIRST SUBSCRIPT IN V OF THE S LOWER TRIANGLE OF 79C THE S MATRIX) WILL HAVE BEEN SET (UNLESS LIV OR LV IS TOO SMALL), 80C AND IV(1) WILL HAVE BEEN SET TO 14. CALLING SRN2G WITH IV(1) = 14 81C CAUSES EXECUTION OF THE ALGORITHM TO BEGIN UNDER THE ASSUMPTION 82C THAT STORAGE HAS BEEN ALLOCATED. 83C 84C *** SUPPLYING R AND DR *** 85C 86C SRN2G USES IV AND V IN THE SAME WAY AS NL2SOL, WITH A SMALL 87C NUMBER OF OBVIOUS CHANGES. ONE DIFFERENCE BETWEEN SRN2G AND 88C NL2ITR IS THAT INITIAL FUNCTION AND GRADIENT INFORMATION NEED NOT 89C BE SUPPLIED IN THE VERY FIRST CALL ON SRN2G, THE ONE WITH 90C IV(1) = 0 OR 12. ANOTHER DIFFERENCE IS THAT SRN2G RETURNS WITH 91C IV(1) = -2 WHEN IT WANTS ANOTHER LOOK AT THE OLD JACOBIAN MATRIX 92C AND THE CURRENT RESIDUAL -- THE ONE CORRESPONDING TO X AND 93C IV(NFGCAL). IT THEN RETURNS WITH IV(1) = -3 WHEN IT WANTS TO SEE 94C BOTH THE NEW RESIDUAL AND THE NEW JACOBIAN MATRIX AT ONCE. NOTE 95C THAT IV(NFGCAL) = IV(7) CONTAINS THE VALUE THAT IV(NFCALL) = IV(6) 96C HAD WHEN THE CURRENT RESIDUAL WAS EVALUATED. ALSO NOTE THAT THE 97C VALUE OF X CORRESPONDING TO THE OLD JACOBIAN MATRIX IS STORED IN 98C V, STARTING AT V(IV(X0)) = V(IV(43)). 99C ANOTHER NEW RETURN... SRN2G IV(1) = -1 WHEN IT WANTS BOTH THE 100C RESIDUAL AND THE JACOBIAN TO BE EVALUATED AT X. 101C A NEW RESIDUAL VECTOR MUST BE SUPPLIED WHEN SRN2G RETURNS WITH 102C IV(1) = 1 OR -1. THIS TAKES THE FORM OF VALUES OF R(I,X) PASSED 103C IN R(I-N1+1), I = N1(1)N2. YOU MAY PASS ALL THESE VALUES AT ONCE 104C (I.E., N1 = 1 AND N2 = N) OR IN PIECES BY MAKING SEVERAL CALLS ON 105C SRN2G. EACH TIME SRN2G RETURNS WITH IV(1) = 1, N1 WILL HAVE 106C BEEN SET TO THE INDEX OF THE NEXT RESIDUAL THAT SRN2G EXPECTS TO 107C SEE, AND N2 WILL BE SET TO THE INDEX OF THE HIGHEST RESIDUAL THAT 108C COULD BE GIVEN ON THE NEXT CALL, I.E., N2 = N1 + ND - 1. (THUS 109C WHEN SRN2G FIRST RETURNS WITH IV(1) = 1 FOR A NEW X, IT WILL 110C HAVE SET N1 TO 1 AND N2 TO MIN(ND,N).) THE CALLER MAY PROVIDE 111C FEWER THAN N2-N1+1 RESIDUALS ON THE NEXT CALL BY SETTING N2 TO 112C A SMALLER VALUE. SRN2G ASSUMES IT HAS SEEN ALL THE RESIDUALS 113C FOR THE CURRENT X WHEN IT IS CALLED WITH N2 .GE. N. 114C EXAMPLE... SUPPOSE N = 80 AND THAT R IS TO BE PASSED IN 8 115C BLOCKS OF SIZE 10. THE FOLLOWING CODE WOULD DO THE JOB. 116C 117C N = 80 118C ND = 10 119C ... 120C DO 10 K = 1, 8 121C *** COMPUTE R(I,X) FOR I = 10*K-9 TO 10*K *** 122C *** AND STORE THEM IN R(1),...,R(10) *** 123C CALL SRN2G(..., R, ...) 124C 10 CONTINUE 125C 126C THE SITUATION IS SIMILAR WHEN GRADIENT INFORMATION IS 127C REQUIRED, I.E., WHEN SRN2G RETURNS WITH IV(1) = 2, -1, OR -2. 128C NOTE THAT SRN2G OVERWRITES R, BUT THAT IN THE SPECIAL CASE OF 129C N1 = 1 AND N2 = N ON PREVIOUS CALLS, SRN2G NEVER RETURNS WITH 130C IV(1) = -2. IT SHOULD BE CLEAR THAT THE PARTIAL DERIVATIVE OF 131C R(I,X) WITH RESPECT TO X(L) IS TO BE STORED IN DR(I-N1+1,L), 132C L = 1(1)P, I = N1(1)N2. IT IS ESSENTIAL THAT R(I) AND DR(I,L) 133C ALL CORRESPOND TO THE SAME RESIDUALS WHEN IV(1) = -1 OR -2. 134C 135C *** COVARIANCE MATRIX *** 136C 137C IV(RDREQ) = IV(57) TELLS WHETHER TO COMPUTE A COVARIANCE 138C MATRIX AND/OR REGRESSION DIAGNOSTICS... 0 MEANS NEITHER, 139C 1 MEANS COVARIANCE MATRIX ONLY, 2 MEANS REG. DIAGNOSTICS ONLY, 140C 3 MEANS BOTH. AS WITH NL2SOL, IV(COVREQ) = IV(15) TELLS WHAT 141C HESSIAN APPROXIMATION TO USE IN THIS COMPUTING. 142C 143C *** REGRESSION DIAGNOSTICS *** 144C 145C SEE THE COMMENTS IN SUBROUTINE SN2G. 146C 147C *** GENERAL *** 148C 149C CODED BY DAVID M. GAY. 150C 151C ++++++++++++++++++++++++++++ DECLARATIONS ++++++++++++++++++++++++++ 152C 153C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** 154C 155 EXTERNAL SC7VFN,SIVSET, SD7TPR,SD7UPD,SG7LIT,SITSUM,SL7VML, 156 1 SN2CVP, SN2LRD, SQ7APL,SQ7RAD,SV7CPY, SV7SCP, SV2NRM 157 REAL SD7TPR, SV2NRM 158c ------------------------------------------------------------------ 159C SC7VFN... FINISHES COVARIANCE COMPUTATION. 160C SIVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS. 161C SD7TPR... COMPUTES INNER PRODUCT OF TWO VECTORS. 162C SD7UPD... UPDATES SCALE VECTOR D. 163C SG7LIT.... PERFORMS BASIC MINIMIZATION ALGORITHM. 164C SITSUM.... PRINTS ITERATION SUMMARY, INFO ABOUT INITIAL AND FINAL X. 165C SL7VML.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX. 166C SN2CVP... PRINTS COVARIANCE MATRIX. 167C SN2LRD... COMPUTES REGRESSION DIAGNOSTICS. 168C SQ7APL... APPLIES QR TRANSFORMATIONS STORED BY SQ7RAD. 169C SQ7RAD.... ADDS A NEW BLOCK OF ROWS TO QR DECOMPOSITION. 170C SV7CPY.... COPIES ONE VECTOR TO ANOTHER. 171C SV7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR. 172C 173C *** LOCAL VARIABLES *** 174C 175 INTEGER G1, GI, I, IV1, IVMODE, JTOL1, K, L, LH, NN, QTR1, 176 1 RMAT1, YI, Y1 177 REAL T 178C 179 REAL HALF, ZERO 180C 181C *** SUBSCRIPTS FOR IV AND V *** 182C 183 INTEGER CNVCOD, COVMAT, COVREQ, DINIT, DTYPE, DTINIT, D0INIT, F, 184 1 FDH, G, H, IPIVOT, IVNEED, JCN, JTOL, LMAT, MODE, 185 2 NEXTIV, NEXTV, NF0, NF00, NF1, NFCALL, NFCOV, NFGCAL, 186 3 NGCALL, NGCOV, QTR, RDREQ, REGD, RESTOR, RLIMIT, RMAT, 187 4 TOOBIG, VNEED, Y 188C 189C *** IV SUBSCRIPT VALUES *** 190C 191 PARAMETER (CNVCOD=55, COVMAT=26, COVREQ=15, DTYPE=16, FDH=74, 192 1 G=28, H=56, IPIVOT=76, IVNEED=3, JCN=66, JTOL=59, 193 2 LMAT=42, MODE=35, NEXTIV=46, NEXTV=47, NFCALL=6, 194 3 NFCOV=52, NF0=68, NF00=81, NF1=69, NFGCAL=7, NGCALL=30, 195 4 NGCOV=53, QTR=77, RESTOR=9, RMAT=78, RDREQ=57, REGD=67, 196 5 TOOBIG=2, VNEED=4, Y=48) 197C 198C *** V SUBSCRIPT VALUES *** 199 PARAMETER (DINIT=38, DTINIT=39, D0INIT=40, F=10, RLIMIT=46) 200 PARAMETER (HALF=0.5E+0, ZERO=0.E+0) 201C 202C ++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ 203C 204 LH = P * (P+1) / 2 205 IF (IV(1) .EQ. 0) CALL SIVSET(1, IV, LIV, LV, V) 206 IV1 = IV(1) 207 IF (IV1 .GT. 2) GO TO 10 208 NN = N2 - N1 + 1 209 IV(RESTOR) = 0 210 I = IV1 + 4 211 IF (IV(TOOBIG) .EQ. 0) GO TO (150, 130, 150, 120, 120, 150), I 212 IF (I .NE. 5) IV(1) = 2 213 GO TO 40 214C 215C *** FRESH START OR RESTART -- CHECK INPUT INTEGERS *** 216C 217 10 IF (ND .LE. 0) GO TO 210 218 IF (P .LE. 0) GO TO 210 219 IF (N .LE. 0) GO TO 210 220 IF (IV1 .EQ. 14) GO TO 30 221 IF (IV1 .GT. 16) GO TO 300 222 IF (IV1 .LT. 12) GO TO 40 223 IF (IV1 .EQ. 12) IV(1) = 13 224 IF (IV(1) .NE. 13) GO TO 20 225 IV(IVNEED) = IV(IVNEED) + P 226 IV(VNEED) = IV(VNEED) + P*(P+13)/2 227 20 CALL SG7LIT(D, X, IV, LIV, LV, P, P, V, X, X) 228 IF (IV(1) .NE. 14) GO TO 999 229C 230C *** STORAGE ALLOCATION *** 231C 232 IV(IPIVOT) = IV(NEXTIV) 233 IV(NEXTIV) = IV(IPIVOT) + P 234 IV(Y) = IV(NEXTV) 235 IV(G) = IV(Y) + P 236 IV(JCN) = IV(G) + P 237 IV(RMAT) = IV(JCN) + P 238 IV(QTR) = IV(RMAT) + LH 239 IV(JTOL) = IV(QTR) + P 240 IV(NEXTV) = IV(JTOL) + 2*P 241 IF (IV1 .EQ. 13) GO TO 999 242C 243 30 JTOL1 = IV(JTOL) 244 IF (V(DINIT) .GE. ZERO) CALL SV7SCP(P, D, V(DINIT)) 245 IF (V(DTINIT) .GT. ZERO) CALL SV7SCP(P, V(JTOL1), V(DTINIT)) 246 I = JTOL1 + P 247 IF (V(D0INIT) .GT. ZERO) CALL SV7SCP(P, V(I), V(D0INIT)) 248 IV(NF0) = 0 249 IV(NF1) = 0 250 IF (ND .GE. N) GO TO 40 251C 252C *** SPECIAL CASE HANDLING OF FIRST FUNCTION AND GRADIENT EVALUATION 253C *** -- ASK FOR BOTH RESIDUAL AND JACOBIAN AT ONCE 254C 255 G1 = IV(G) 256 Y1 = IV(Y) 257 CALL SG7LIT(D, V(G1), IV, LIV, LV, P, P, V, X, V(Y1)) 258 IF (IV(1) .NE. 1) GO TO 220 259 V(F) = ZERO 260 CALL SV7SCP(P, V(G1), ZERO) 261 IV(1) = -1 262 QTR1 = IV(QTR) 263 CALL SV7SCP(P, V(QTR1), ZERO) 264 IV(REGD) = 0 265 RMAT1 = IV(RMAT) 266 GO TO 100 267C 268 40 G1 = IV(G) 269 Y1 = IV(Y) 270 CALL SG7LIT(D, V(G1), IV, LIV, LV, P, P, V, X, V(Y1)) 271 IF (IV(1) - 2) 50, 60, 220 272C 273 50 V(F) = ZERO 274 IF (IV(NF1) .EQ. 0) GO TO 260 275 IF (IV(RESTOR) .NE. 2) GO TO 260 276 IV(NF0) = IV(NF1) 277 CALL SV7CPY(N, RD, R) 278 IV(REGD) = 0 279 GO TO 260 280C 281 60 CALL SV7SCP(P, V(G1), ZERO) 282 IF (IV(MODE) .GT. 0) GO TO 230 283 RMAT1 = IV(RMAT) 284 QTR1 = IV(QTR) 285 CALL SV7SCP(P, V(QTR1), ZERO) 286 IV(REGD) = 0 287 IF (ND .LT. N) GO TO 90 288 IF (N1 .NE. 1) GO TO 90 289 IF (IV(MODE) .LT. 0) GO TO 100 290 IF (IV(NF1) .EQ. IV(NFGCAL)) GO TO 70 291 IF (IV(NF0) .NE. IV(NFGCAL)) GO TO 90 292 CALL SV7CPY(N, R, RD) 293 GO TO 80 294 70 CALL SV7CPY(N, RD, R) 295 80 CALL SQ7APL(ND, N, P, DR, RD, 0) 296 CALL SL7VML(P, V(Y1), V(RMAT1), RD) 297 GO TO 110 298C 299 90 IV(1) = -2 300 IF (IV(MODE) .LT. 0) IV(1) = -1 301 100 CALL SV7SCP(P, V(Y1), ZERO) 302 110 CALL SV7SCP(LH, V(RMAT1), ZERO) 303 GO TO 260 304C 305C *** COMPUTE F(X) *** 306C 307 120 T = SV2NRM(NN, R) 308 IF (T .GT. V(RLIMIT)) GO TO 200 309 V(F) = V(F) + HALF * T**2 310 IF (N2 .LT. N) GO TO 270 311 IF (N1 .EQ. 1) IV(NF1) = IV(NFCALL) 312 GO TO 40 313C 314C *** COMPUTE Y *** 315C 316 130 Y1 = IV(Y) 317 YI = Y1 318 DO 140 L = 1, P 319 V(YI) = V(YI) + SD7TPR(NN, DR(1,L), R) 320 YI = YI + 1 321 140 CONTINUE 322 IF (N2 .LT. N) GO TO 270 323 IV(1) = 2 324 IF (N1 .GT. 1) IV(1) = -3 325 GO TO 260 326C 327C *** COMPUTE GRADIENT INFORMATION *** 328C 329 150 IF (IV(MODE) .GT. P) GO TO 240 330 G1 = IV(G) 331 IVMODE = IV(MODE) 332 IF (IVMODE .LT. 0) GO TO 170 333 IF (IVMODE .EQ. 0) GO TO 180 334 IV(1) = 2 335C 336C *** COMPUTE GRADIENT ONLY (FOR USE IN COVARIANCE COMPUTATION) *** 337C 338 GI = G1 339 DO 160 L = 1, P 340 V(GI) = V(GI) + SD7TPR(NN, R, DR(1,L)) 341 GI = GI + 1 342 160 CONTINUE 343 GO TO 190 344C 345C *** COMPUTE INITIAL FUNCTION VALUE WHEN ND .LT. N *** 346C 347 170 IF (N .LE. ND) GO TO 180 348 T = SV2NRM(NN, R) 349 IF (T .GT. V(RLIMIT)) GO TO 200 350 V(F) = V(F) + HALF * T**2 351C 352C *** UPDATE D IF DESIRED *** 353C 354 180 IF (IV(DTYPE) .GT. 0) 355 1 CALL SD7UPD(D, DR, IV, LIV, LV, N, ND, NN, N2, P, V) 356C 357C *** COMPUTE RMAT AND QTR *** 358C 359 QTR1 = IV(QTR) 360 RMAT1 = IV(RMAT) 361 CALL SQ7RAD(NN, ND, P, V(QTR1), .TRUE., V(RMAT1), DR, R) 362 IV(NF1) = 0 363C 364 190 IF (N2 .LT. N) GO TO 270 365 IF (IVMODE .GT. 0) GO TO 40 366 IV(NF00) = IV(NFGCAL) 367C 368C *** COMPUTE G FROM RMAT AND QTR *** 369C 370 CALL SL7VML(P, V(G1), V(RMAT1), V(QTR1)) 371 IV(1) = 2 372 IF (IVMODE .EQ. 0) GO TO 40 373 IF (N .LE. ND) GO TO 40 374C 375C *** FINISH SPECIAL CASE HANDLING OF FIRST FUNCTION AND GRADIENT 376C 377 Y1 = IV(Y) 378 IV(1) = 1 379 CALL SG7LIT(D, V(G1), IV, LIV, LV, P, P, V, X, V(Y1)) 380 IF (IV(1) .NE. 2) GO TO 220 381 GO TO 40 382C 383C *** MISC. DETAILS *** 384C 385C *** X IS OUT OF RANGE (OVERSIZE STEP) *** 386C 387 200 IV(TOOBIG) = 1 388 GO TO 40 389C 390C *** BAD N, ND, OR P *** 391C 392 210 IV(1) = 66 393 GO TO 300 394C 395C *** CONVERGENCE OBTAINED -- SEE WHETHER TO COMPUTE COVARIANCE *** 396C 397 220 IF (IV(COVMAT) .NE. 0) GO TO 290 398 IF (IV(REGD) .NE. 0) GO TO 290 399C 400C *** SEE IF CHOLESKY FACTOR OF HESSIAN IS AVAILABLE *** 401C 402 K = IV(FDH) 403 IF (K .LE. 0) GO TO 280 404 IF (IV(RDREQ) .LE. 0) GO TO 290 405C 406C *** COMPUTE REGRESSION DIAGNOSTICS AND DEFAULT COVARIANCE IF 407C DESIRED *** 408C 409 I = 0 410 IF (mod(IV(RDREQ),4) .GE. 2) I = 1 411 IF (mod(IV(RDREQ),2) .EQ. 1 .AND. abs(IV(COVREQ)) .LE. 1) I = I+2 412 IF (I .EQ. 0) GO TO 250 413 IV(MODE) = P + I 414 IV(NGCALL) = IV(NGCALL) + 1 415 IV(NGCOV) = IV(NGCOV) + 1 416 IV(CNVCOD) = IV(1) 417 IF (I .LT. 2) GO TO 230 418 L = abs(IV(H)) 419 CALL SV7SCP(LH, V(L), ZERO) 420 230 IV(NFCOV) = IV(NFCOV) + 1 421 IV(NFCALL) = IV(NFCALL) + 1 422 IV(NFGCAL) = IV(NFCALL) 423 IV(1) = -1 424 GO TO 260 425C 426 240 L = IV(LMAT) 427 CALL SN2LRD(DR, IV, V(L), LH, LIV, LV, ND, NN, P, R, RD, V) 428 IF (N2 .LT. N) GO TO 270 429 IF (N1 .GT. 1) GO TO 250 430C 431C *** ENSURE WE CAN RESTART -- AND MAKE RETURN STATE OF DR 432C *** INDEPENDENT OF WHETHER REGRESSION DIAGNOSTICS ARE COMPUTED. 433C *** USE STEP VECTOR (ALLOCATED BY SG7LIT) FOR SCRATCH. 434C 435 RMAT1 = IV(RMAT) 436 CALL SV7SCP(LH, V(RMAT1), ZERO) 437 CALL SQ7RAD(NN, ND, P, R, .FALSE., V(RMAT1), DR, R) 438 IV(NF1) = 0 439C 440C *** FINISH COMPUTING COVARIANCE *** 441C 442 250 L = IV(LMAT) 443 CALL SC7VFN(IV, V(L), LH, LIV, LV, N, P, V) 444 GO TO 290 445C 446C *** RETURN FOR MORE FUNCTION OR GRADIENT INFORMATION *** 447C 448 260 N2 = 0 449 270 N1 = N2 + 1 450 N2 = N2 + ND 451 IF (N2 .GT. N) N2 = N 452 GO TO 999 453C 454C *** COME HERE FOR INDEFINITE FINITE-DIFFERENCE HESSIAN *** 455C 456 280 IV(COVMAT) = K 457 IV(REGD) = K 458C 459C *** PRINT SUMMARY OF FINAL ITERATION AND OTHER REQUESTED ITEMS *** 460C 461 290 G1 = IV(G) 462 300 CALL SITSUM(D, V(G1), IV, LIV, LV, P, V, X) 463 IF (IV(1) .LE. 6 .AND. IV(RDREQ) .GT. 0) 464 1 CALL SN2CVP(IV, LIV, LV, P, V) 465C 466 999 RETURN 467C *** LAST LINE OF SRN2G FOLLOWS *** 468 END 469c ================================================================== 470 SUBROUTINE SG7LIT(D, GG, IV, LIV, LV, P, PS, V, X, YY) 471c>> 1990-06-12 CLL @ JPL 472c>> 1990-04-23 CLL (Recent revision by DMG) 473*** from netlib, Mon Apr 23 20:37:24 EDT 1990 *** 474c>> 1990-02-20 CLL @ JPL 475C 476C *** CARRY OUT NL2SOL-LIKE ITERATIONS FOR GENERALIZED LINEAR *** 477C *** REGRESSION PROBLEMS (AND OTHERS OF SIMILAR STRUCTURE) *** 478C 479C *** PARAMETER DECLARATIONS *** 480C 481 INTEGER LIV, LV, P, PS 482 INTEGER IV(LIV) 483 REAL D(P), GG(P), V(LV), X(P), YY(P) 484C 485C ------------------------- PARAMETER USAGE -------------------------- 486C 487C D.... SCALE VECTOR. 488C IV... INTEGER VALUE ARRAY. 489C LIV.. LENGTH OF IV. MUST BE AT LEAST 82. 490C LH... LENGTH OF H = P*(P+1)/2. 491C LV... LENGTH OF V. MUST BE AT LEAST P*(3*P + 19)/2 + 7. 492C GG... GRADIENT AT X (WHEN IV(1) = 2). 493C P.... NUMBER OF PARAMETERS (COMPONENTS IN X). 494C PS... NUMBER OF NONZERO ROWS AND COLUMNS IN S. 495C V.... FLOATING-POINT VALUE ARRAY. 496C X.... PARAMETER VECTOR. 497C YY... PART OF YIELD VECTOR (WHEN IV(1)= 2, SCRATCH OTHERWISE). 498C 499C *** DISCUSSION *** 500C 501C SG7LIT PERFORMS NL2SOL-LIKE ITERATIONS FOR A VARIETY OF 502C REGRESSION PROBLEMS THAT ARE SIMILAR TO NONLINEAR LEAST-SQUARES 503C IN THAT THE HESSIAN IS THE SUM OF TWO TERMS, A READILY-COMPUTED 504C FIRST-ORDER TERM AND A SECOND-ORDER TERM. THE CALLER SUPPLIES 505C THE FIRST-ORDER TERM OF THE HESSIAN IN HC (LOWER TRIANGLE, STORED 506C COMPACTLY BY ROWS IN V, STARTING AT IV(HC)), AND SG7LIT BUILDS AN 507C APPROXIMATION, S, TO THE SECOND-ORDER TERM. THE CALLER ALSO 508C PROVIDES THE FUNCTION VALUE, GRADIENT, AND PART OF THE YIELD 509C VECTOR USED IN UPDATING S. SG7LIT DECIDES DYNAMICALLY WHETHER OR 510C NOT TO USE S WHEN CHOOSING THE NEXT STEP TO TRY... THE HESSIAN 511C APPROXIMATION USED IS EITHER HC ALONE (GAUSS-NEWTON MODEL) OR 512C HC + S (AUGMENTED MODEL). 513C 514C IF PS .LT. P, THEN ROWS AND COLUMNS PS+1...P OF S ARE KEPT 515C CONSTANT. THEY WILL BE ZERO UNLESS THE CALLER SETS IV(INITS) TO 516C 1 OR 2 AND SUPPLIES NONZERO VALUES FOR THEM, OR THE CALLER SETS 517C IV(INITS) TO 3 OR 4 AND THE FINITE-DIFFERENCE INITIAL S THEN 518C COMPUTED HAS NONZERO VALUES IN THESE ROWS. 519C 520C IF IV(INITS) IS 3 OR 4, THEN THE INITIAL S IS COMPUTED BY 521C FINITE DIFFERENCES. 3 MEANS USE FUNCTION DIFFERENCES, 4 MEANS 522C USE GRADIENT DIFFERENCES. FINITE DIFFERENCING IS DONE THE SAME 523C WAY AS IN COMPUTING A COVARIANCE MATRIX (WITH IV(COVREQ) = -1, -2, 524C 1, OR 2). 525C 526C FOR UPDATING S,SG7LIT ASSUMES THAT THE GRADIENT HAS THE FORM 527C OF A SUM OVER I OF RHO(I,X)*GRAD(R(I,X)), WHERE GRAD DENOTES THE 528C GRADIENT WITH RESPECT TO X. THE TRUE SECOND-ORDER TERM THEN IS 529C THE SUM OVER I OF RHO(I,X)*HESSIAN(R(I,X)). IF X = X0 + STEP, 530C THEN WE WISH TO UPDATE S SO THAT S*STEP IS THE SUM OVER I OF 531C RHO(I,X)*(GRAD(R(I,X)) - GRAD(R(I,X0))). THE CALLER MUST SUPPLY 532C PART OF THIS IN YY, NAMELY THE SUM OVER I OF 533C RHO(I,X)*GRAD(R(I,X0)), WHEN CALLING SG7LIT WITH IV(1) = 2 AND 534C IV(MODE) = 0 (WHERE MODE = 38). GG THEN CONTANS THE OTHER PART, 535C SO THAT THE DESIRED YIELD VECTOR IS GG - YY. IF PS .LT. P, THEN 536C THE ABOVE DISCUSSION APPLIES ONLY TO THE FIRST PS COMPONENTS OF 537C GRAD(R(I,X)), STEP, AND YY. 538C 539C PARAMETERS IV, P, V, AND X ARE THE SAME AS THE CORRESPONDING 540C ONES TO NL2SOL (WHICH SEE), EXCEPT THAT V CAN BE SHORTER 541C (SINCE THE PART OF V THAT NL2SOL USES FOR STORING D, J, AND R IS 542C NOT NEEDED). MOREOVER, COMPARED WITH NL2SOL, IV(1) MAY HAVE THE 543C TWO ADDITIONAL OUTPUT VALUES 1 AND 2, WHICH ARE EXPLAINED BELOW, 544C AS IS THE USE OF IV(TOOBIG) AND IV(NFGCAL). THE VALUES IV(D), 545C IV(J), AND IV(R), WHICH ARE OUTPUT VALUES FROM NL2SOL (AND 546C NL2SNO), ARE NOT REFERENCED BY SG7LIT OR THE SUBROUTINES IT CALLS. 547C 548C WHEN SG7LIT IS FIRST CALLED, I.E., WHEN SG7LIT IS CALLED WITH 549C IV(1) = 0 OR 12, V(F), GG, AND HC NEED NOT BE INITIALIZED. TO 550C OBTAIN THESE STARTING VALUES,SG7LIT RETURNS FIRST WITH IV(1) = 1, 551C THEN WITH IV(1) = 2, WITH IV(MODE) = -1 IN BOTH CASES. ON 552C SUBSEQUENT RETURNS WITH IV(1) = 2, IV(MODE) = 0 IMPLIES THAT 553C YY MUST ALSO BE SUPPLIED. (NOTE THAT YY IS USED FOR SCRATCH -- 554c ITS INPUT CONTENTS ARE LOST. BY CONTRAST, HC IS NEVER CHANGED.) 555C ONCE CONVERGENCE HAS BEEN OBTAINED, IV(RDREQ) AND IV(COVREQ) MAY 556C IMPLY THAT A FINITE-DIFFERENCE HESSIAN SHOULD BE COMPUTED FOR USE 557C IN COMPUTING A COVARIANCE MATRIX. IN THIS CASE SG7LIT WILL MAKE A 558C NUMBER OF RETURNS WITH IV(1) = 1 OR 2 AND IV(MODE) POSITIVE. 559C WHEN IV(MODE) IS POSITIVE, YY SHOULD NOT BE CHANGED. 560C 561C IV(1) = 1 MEANS THE CALLER SHOULD SET V(F) (I.E., V(10)) TO F(X), THE 562C FUNCTION VALUE AT X, AND CALL SG7LIT AGAIN, HAVING CHANGED 563C NONE OF THE OTHER PARAMETERS. AN EXCEPTION OCCURS IF F(X) 564C CANNOT BE EVALUATED (E.G. IF OVERFLOW WOULD OCCUR), WHICH 565C MAY HAPPEN BECAUSE OF AN OVERSIZED STEP. IN THIS CASE 566C THE CALLER SHOULD SET IV(TOOBIG) = IV(2) TO 1, WHICH WILL 567C CAUSE SG7LIT TO IGNORE V(F) AND TRY A SMALLER STEP. NOTE 568C THAT THE CURRENT FUNCTION EVALUATION COUNT IS AVAILABLE 569C IN IV(NFCALL) = IV(6). THIS MAY BE USED TO IDENTIFY 570C WHICH COPY OF SAVED INFORMATION SHOULD BE USED IN COM- 571C PUTING GG, HC, AND YY THE NEXT TIME SG7LIT RETURNS WITH 572C IV(1) = 2. SEE MLPIT FOR AN EXAMPLE OF THIS. 573C IV(1) = 2 MEANS THE CALLER SHOULD SET GG TO G(X), THE GRADIENT OF F AT 574C X. THE CALLER SHOULD ALSO SET HC TO THE GAUSS-NEWTON 575C HESSIAN AT X. IF IV(MODE) = 0, THEN THE CALLER SHOULD 576C ALSO COMPUTE THE PART OF THE YIELD VECTOR DESCRIBED ABOVE. 577C THE CALLER SHOULD THEN CALL SG7LIT AGAIN (WITH IV(1) = 2). 578C THE CALLER MAY ALSO CHANGE D AT THIS TIME, BUT SHOULD NOT 579C CHANGE X. NOTE THAT IV(NFGCAL) = IV(7) CONTAINS THE 580C VALUE THAT IV(NFCALL) HAD DURING THE RETURN WITH 581C IV(1) = 1 IN WHICH X HAD THE SAME VALUE AS IT NOW HAS. 582C IV(NFGCAL) IS EITHER IV(NFCALL) OR IV(NFCALL) - 1. MLPIT 583C IS AN EXAMPLE WHERE THIS INFORMATION IS USED. IF GG OR HC 584C CANNOT BE EVALUATED AT X, THEN THE CALLER MAY SET 585C IV(TOOBIG) TO 1, IN WHICH CASE SG7LIT WILL RETURN WITH 586C IV(1) = 15. 587C 588C *** GENERAL *** 589C 590C CODED BY DAVID M. GAY. 591C THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH 592C SUPPORTED IN PART BY D.O.E. GRANT EX-76-A-01-2295 TO MIT/CCREMS. 593C 594C (SEE NL2SOL FOR REFERENCES.) 595C 596c ------------------------------------------------------------------ 597c References to the function STOPX have been commented out of this 598c subroutine. If one wishes to be able to terminate this package 599c gracefully using a keybord "Break" key, one can provide a STOPX 600c function that returns .true. if the Break key has been pressed 601c since the last call to STOPX, and otherwise returns .false., and 602c then uncomment the references to STOPX in this subr. 603c -- CLL 6/12/90 604C ++++++++++++++++++++++++++ DECLARATIONS ++++++++++++++++++++++++++++ 605C 606C *** LOCAL VARIABLES *** 607C 608c integer DUMMY 609 INTEGER DIG1, G01, H1, HC1, I, IPIV1, J, K, L, LMAT1, 610 1 LSTGST, PP1O2, QTR1, RMAT1, RSTRST, STEP1, STPMOD, S1, 611 2 TEMP1, TEMP2, W1, X01 612 REAL E, STTSST, T, T1, TP 613C 614C *** CONSTANTS *** 615C 616 REAL HALF, NEGONE, ONE, ONEP2, ZERO 617C 618C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** 619C 620c external STOPX 621c LOGICAL STOPX 622 EXTERNAL SA7SST, SD7TPR,SF7HES,SG7QTS,SITSUM, SL7MST,SL7SRT, 623 1 SL7SQR, SL7SVX, SL7SVN, SL7TVM,SL7VML,SPARCK, SRLDST, 624 2 SR7MDC, SS7LUP, SS7LVM, SV2AXY,SV7CPY, SV7SCP, 625 3 SV2NRM 626 REAL SD7TPR, SL7SVX, SL7SVN, SRLDST, SR7MDC, SV2NRM 627c ------------------------------------------------------------------ 628C SA7SST.... ASSESSES CANDIDATE STEP. 629C SD7TPR... RETURNS INNER PRODUCT OF TWO VECTORS. 630C SF7HES.... COMPUTE FINITE-DIFFERENCE HESSIAN (FOR COVARIANCE). 631C SG7QTS.... COMPUTES GOLDFELD-QUANDT-TROTTER STEP (AUGMENTED MODEL). 632C SITSUM.... PRINTS ITERATION SUMMARY AND INFO ON INITIAL AND FINAL X. 633C SL7MST... COMPUTES LEVENBERG-MARQUARDT STEP (GAUSS-NEWTON MODEL). 634C SL7SRT.... COMPUTES CHOLESKY FACTOR OF (LOWER TRIANG. OF) SYM. MATRIX. 635C SL7SQR... COMPUTES L * L**T FROM LOWER TRIANGULAR MATRIX L. 636C SL7TVM... COMPUTES L**T * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX. 637C SL7SVX... ESTIMATES LARGEST SING. VALUE OF LOWER TRIANG. MATRIX. 638C SL7SVN... ESTIMATES SMALLEST SING. VALUE OF LOWER TRIANG. MATRIX. 639C SL7VML.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX. 640C SPARCK.... CHECK VALIDITY OF IV AND V INPUT COMPONENTS. 641C SRLDST... COMPUTES V(RELDX) = RELATIVE STEP SIZE. 642C SR7MDC... RETURNS MACHINE-DEPENDENT CONSTANTS. 643C SS7LUP... PERFORMS QUASI-NEWTON UPDATE ON COMPACTLY STORED LOWER TRI- 644C ANGLE OF A SYMMETRIC MATRIX. 645C STOPX.... RETURNS .TRUE. IF THE BREAK KEY HAS BEEN PRESSED. 646c Call to STOPX commented out. -- CLL 6/12/90 647C SV2AXY.... COMPUTES SCALAR TIMES ONE VECTOR PLUS ANOTHER. 648C SV7CPY.... COPIES ONE VECTOR TO ANOTHER. 649C SV7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR. 650C SV2NRM... RETURNS THE 2-NORM OF A VECTOR. 651C 652C *** SUBSCRIPTS FOR IV AND V *** 653C 654 INTEGER CNVCOD, COSMIN, COVMAT, COVREQ, DGNORM, DIG, DSTNRM, F, 655 1 FDH, FDIF, FUZZ, F0, GTSTEP, H, HC, IERR, INCFAC, INITS, 656 2 IPIVOT, IRC, KAGQT, KALM, LMAT, LMAX0, LMAXS, MODE, MODEL, 657 3 MXFCAL, MXITER, NEXTV, NFCALL, NFGCAL, NFCOV, NGCOV, 658 4 NGCALL, NITER, NVSAVE, PHMXFC, PREDUC, QTR, RADFAC, 659 5 RADINC, RADIUS, RAD0, RCOND, RDREQ, REGD, RELDX, RESTOR, 660 6 RMAT, S, SIZE, STEP, STGLIM, STLSTG, STPPAR, SUSED, 661 7 SWITCH, TOOBIG, TUNER4, TUNER5, VNEED, VSAVE, W, WSCALE, 662 8 XIRC, X0 663C 664C *** IV SUBSCRIPT VALUES *** 665C 666 PARAMETER (CNVCOD=55, COVMAT=26, COVREQ=15, DIG=37, FDH=74, H=56, 667 1 HC=71, IERR=75, INITS=25, IPIVOT=76, IRC=29, KAGQT=33, 668 2 KALM=34, LMAT=42, MODE=35, MODEL=5, MXFCAL=17, 669 3 MXITER=18, NEXTV=47, NFCALL=6, NFGCAL=7, NFCOV=52, 670 4 NGCOV=53, NGCALL=30, NITER=31, QTR=77, RADINC=8, 671 5 RDREQ=57, REGD=67, RESTOR=9, RMAT=78, S=62, STEP=40, 672 6 STGLIM=11, STLSTG=41, SUSED=64, SWITCH=12, TOOBIG=2, 673 7 VNEED=4, VSAVE=60, W=65, XIRC=13, X0=43) 674C 675C *** V SUBSCRIPT VALUES *** 676C 677 PARAMETER (COSMIN=47, DGNORM=1, DSTNRM=2, F=10, FDIF=11, FUZZ=45, 678 1 F0=13, GTSTEP=4, INCFAC=23, LMAX0=35, LMAXS=36, 679 2 NVSAVE=9, PHMXFC=21, PREDUC=7, RADFAC=16, RADIUS=8, 680 3 RAD0=9, RCOND=53, RELDX=17, SIZE=55, STPPAR=5, 681 4 TUNER4=29, TUNER5=30, WSCALE=56) 682 PARAMETER (HALF=0.5E+0, NEGONE=-1.E+0, ONE=1.E+0, ONEP2=1.2E+0, 683 1 ZERO=0.E+0) 684C 685C ++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ 686C 687 I = IV(1) 688 IF (I .EQ. 1) GO TO 40 689 IF (I .EQ. 2) GO TO 50 690C 691 IF (I .EQ. 12 .OR. I .EQ. 13) 692 1 IV(VNEED) = IV(VNEED) + P*(3*P + 19)/2 + 7 693 CALL SPARCK(1, D, IV, LIV, LV, P, V) 694 I = IV(1) - 2 695 IF (I .GT. 12) GO TO 999 696 GO TO (290, 290, 290, 290, 290, 290, 170, 120, 170, 10, 10, 20), I 697C 698C *** STORAGE ALLOCATION *** 699C 700 10 PP1O2 = P * (P + 1) / 2 701 IV(S) = IV(LMAT) + PP1O2 702 IV(X0) = IV(S) + PP1O2 703 IV(STEP) = IV(X0) + P 704 IV(STLSTG) = IV(STEP) + P 705 IV(DIG) = IV(STLSTG) + P 706 IV(W) = IV(DIG) + P 707 IV(H) = IV(W) + 4*P + 7 708 IV(NEXTV) = IV(H) + PP1O2 709 IF (IV(1) .NE. 13) GO TO 20 710 IV(1) = 14 711 GO TO 999 712C 713C *** INITIALIZATION *** 714C 715 20 IV(NITER) = 0 716 IV(NFCALL) = 1 717 IV(NGCALL) = 1 718 IV(NFGCAL) = 1 719 IV(MODE) = -1 720 IV(STGLIM) = 2 721 IV(TOOBIG) = 0 722 IV(CNVCOD) = 0 723 IV(COVMAT) = 0 724 IV(NFCOV) = 0 725 IV(NGCOV) = 0 726 IV(RADINC) = 0 727 IV(RESTOR) = 0 728 IV(FDH) = 0 729 V(RAD0) = ZERO 730 V(STPPAR) = ZERO 731 V(RADIUS) = V(LMAX0) / (ONE + V(PHMXFC)) 732C 733C *** SET INITIAL MODEL AND S MATRIX *** 734C 735 IV(MODEL) = 1 736 IF (IV(S) .LT. 0) GO TO 999 737 IF (IV(INITS) .GT. 1) IV(MODEL) = 2 738 S1 = IV(S) 739 IF (IV(INITS) .EQ. 0 .OR. IV(INITS) .GT. 2) 740 1 CALL SV7SCP(P*(P+1)/2, V(S1), ZERO) 741 IV(1) = 1 742 J = IV(IPIVOT) 743 IF (J .LE. 0) GO TO 999 744 DO 30 I = 1, P 745 IV(J) = I 746 J = J + 1 747 30 CONTINUE 748 GO TO 999 749C 750C *** NEW FUNCTION VALUE *** 751C 752 40 IF (IV(MODE) .EQ. 0) GO TO 290 753 IF (IV(MODE) .GT. 0) GO TO 520 754C 755 IV(1) = 2 756 IF (IV(TOOBIG) .EQ. 0) GO TO 999 757 IV(1) = 63 758 GO TO 999 759C 760C *** NEW GRADIENT *** 761C 762 50 IV(KALM) = -1 763 IV(KAGQT) = -1 764 IV(FDH) = 0 765 IF (IV(MODE) .GT. 0) GO TO 520 766C 767C *** MAKE SURE GRADIENT COULD BE COMPUTED *** 768C 769 IF (IV(TOOBIG) .EQ. 0) GO TO 60 770 IV(1) = 65 771 GO TO 999 772 60 IF (IV(HC) .LE. 0 .AND. IV(RMAT) .LE. 0) GO TO 610 773C 774C *** COMPUTE D**-1 * GRADIENT *** 775C 776 DIG1 = IV(DIG) 777 K = DIG1 778 DO 70 I = 1, P 779 V(K) = GG(I) / D(I) 780 K = K + 1 781 70 CONTINUE 782 V(DGNORM) = SV2NRM(P, V(DIG1)) 783C 784 IF (IV(CNVCOD) .NE. 0) GO TO 510 785 IF (IV(MODE) .EQ. 0) GO TO 440 786 IV(MODE) = 0 787 V(F0) = V(F) 788 IF (IV(INITS) .LE. 2) GO TO 100 789C 790C *** ARRANGE FOR FINITE-DIFFERENCE INITIAL S *** 791C 792 IV(XIRC) = IV(COVREQ) 793 IV(COVREQ) = -1 794 IF (IV(INITS) .GT. 3) IV(COVREQ) = 1 795 IV(CNVCOD) = 70 796 GO TO 530 797C 798C *** COME TO NEXT STMT AFTER COMPUTING F.D. HESSIAN FOR INIT. S *** 799C 800 80 IV(CNVCOD) = 0 801 IV(MODE) = 0 802 IV(NFCOV) = 0 803 IV(NGCOV) = 0 804 IV(COVREQ) = IV(XIRC) 805 S1 = IV(S) 806 PP1O2 = PS * (PS + 1) / 2 807 HC1 = IV(HC) 808 IF (HC1 .LE. 0) GO TO 90 809 CALL SV2AXY(PP1O2, V(S1), NEGONE, V(HC1), V(H1)) 810 GO TO 100 811 90 RMAT1 = IV(RMAT) 812 CALL SL7SQR(PS, V(S1), V(RMAT1)) 813 CALL SV2AXY(PP1O2, V(S1), NEGONE, V(S1), V(H1)) 814 100 IV(1) = 2 815C 816C 817C ---------------------------- MAIN LOOP ----------------------------- 818C 819C 820C *** PRINT ITERATION SUMMARY, CHECK ITERATION LIMIT *** 821C 822 110 CALL SITSUM(D, GG, IV, LIV, LV, P, V, X) 823 120 K = IV(NITER) 824 IF (K .LT. IV(MXITER)) GO TO 130 825 IV(1) = 10 826 GO TO 999 827 130 IV(NITER) = K + 1 828C 829C *** UPDATE RADIUS *** 830C 831 IF (K .EQ. 0) GO TO 150 832 STEP1 = IV(STEP) 833 DO 140 I = 1, P 834 V(STEP1) = D(I) * V(STEP1) 835 STEP1 = STEP1 + 1 836 140 CONTINUE 837 STEP1 = IV(STEP) 838 T = V(RADFAC) * SV2NRM(P, V(STEP1)) 839 IF (V(RADFAC) .LT. ONE .OR. T .GT. V(RADIUS)) V(RADIUS) = T 840C 841C *** INITIALIZE FOR START OF NEXT ITERATION *** 842C 843 150 X01 = IV(X0) 844 V(F0) = V(F) 845 IV(IRC) = 4 846 IV(H) = -abs(IV(H)) 847 IV(SUSED) = IV(MODEL) 848C 849C *** COPY X TO X0 *** 850C 851 CALL SV7CPY(P, V(X01), X) 852C 853C *** CHECK STOPX AND FUNCTION EVALUATION LIMIT *** 854C 855 160 continue 856c if (STOPX(DUMMY)) then 857c IV(1) = 11 858c GO TO 190 859c else 860 go to 180 861c endif 862C 863C *** COME HERE WHEN RESTARTING AFTER FUNC. EVAL. LIMIT OR STOPX. 864C 865 170 IF (V(F) .GE. V(F0)) GO TO 180 866 V(RADFAC) = ONE 867 K = IV(NITER) 868 GO TO 130 869C 870 180 IF (IV(NFCALL) .LT. IV(MXFCAL) + IV(NFCOV)) GO TO 200 871 IV(1) = 9 872c 190 continue 873 IF (V(F) .GE. V(F0)) GO TO 999 874C 875C *** IN CASE OF STOPX OR FUNCTION EVALUATION LIMIT WITH 876C *** IMPROVED V(F), EVALUATE THE GRADIENT AT X. 877C 878 IV(CNVCOD) = IV(1) 879 GO TO 430 880C 881C. . . . . . . . . . . . . COMPUTE CANDIDATE STEP . . . . . . . . . . 882C 883 200 STEP1 = IV(STEP) 884 W1 = IV(W) 885 H1 = IV(H) 886 T1 = ONE 887 IF (IV(MODEL) .EQ. 2) GO TO 210 888 T1 = ZERO 889C 890C *** COMPUTE LEVENBERG-MARQUARDT STEP IF POSSIBLE... 891C 892 RMAT1 = IV(RMAT) 893 IF (RMAT1 .LE. 0) GO TO 210 894 QTR1 = IV(QTR) 895 IF (QTR1 .LE. 0) GO TO 210 896 IPIV1 = IV(IPIVOT) 897 CALL SL7MST(D, GG, IV(IERR), IV(IPIV1), IV(KALM), P, V(QTR1), 898 1 V(RMAT1), V(STEP1), V, V(W1)) 899C *** H IS STORED IN THE END OF W AND HAS JUST BEEN OVERWRITTEN, 900C *** SO WE MARK IT INVALID... 901 IV(H) = -abs(H1) 902C *** EVEN IF H WERE STORED ELSEWHERE, IT WOULD BE NECESSARY TO 903C *** MARK INVALID THE INFORMATION SG7QTS MAY HAVE STORED IN V... 904 IV(KAGQT) = -1 905 GO TO 260 906C 907 210 IF (H1 .GT. 0) GO TO 250 908C 909C *** SET H TO D**-1 * (HC + T1*S) * D**-1. *** 910C 911 H1 = -H1 912 IV(H) = H1 913 IV(FDH) = 0 914 J = IV(HC) 915 IF (J .GT. 0) GO TO 220 916 J = H1 917 RMAT1 = IV(RMAT) 918 CALL SL7SQR(P, V(H1), V(RMAT1)) 919 220 S1 = IV(S) 920 DO 240 I = 1, P 921 T = ONE / D(I) 922 DO 230 K = 1, I 923 V(H1) = T * (V(J) + T1*V(S1)) / D(K) 924 J = J + 1 925 H1 = H1 + 1 926 S1 = S1 + 1 927 230 CONTINUE 928 240 CONTINUE 929 H1 = IV(H) 930 IV(KAGQT) = -1 931C 932C *** COMPUTE ACTUAL GOLDFELD-QUANDT-TROTTER STEP *** 933C 934 250 DIG1 = IV(DIG) 935 LMAT1 = IV(LMAT) 936 CALL SG7QTS(D, V(DIG1), V(H1), IV(KAGQT), V(LMAT1), P, V(STEP1), 937 1 V, V(W1)) 938 IF (IV(KALM) .GT. 0) IV(KALM) = 0 939C 940 260 IF (IV(IRC) .NE. 6) GO TO 270 941 IF (IV(RESTOR) .NE. 2) GO TO 290 942 RSTRST = 2 943 GO TO 300 944C 945C *** CHECK WHETHER EVALUATING F(X0 + STEP) LOOKS WORTHWHILE *** 946C 947 270 IV(TOOBIG) = 0 948 IF (V(DSTNRM) .LE. ZERO) GO TO 290 949 IF (IV(IRC) .NE. 5) GO TO 280 950 IF (V(RADFAC) .LE. ONE) GO TO 280 951 IF (V(PREDUC) .GT. ONEP2 * V(FDIF)) GO TO 280 952 IF (IV(RESTOR) .NE. 2) GO TO 290 953 RSTRST = 0 954 GO TO 300 955C 956C *** COMPUTE F(X0 + STEP) *** 957C 958 280 X01 = IV(X0) 959 STEP1 = IV(STEP) 960 CALL SV2AXY(P, X, ONE, V(STEP1), V(X01)) 961 IV(NFCALL) = IV(NFCALL) + 1 962 IV(1) = 1 963 GO TO 999 964C 965C. . . . . . . . . . . . . ASSESS CANDIDATE STEP . . . . . . . . . . . 966C 967 290 RSTRST = 3 968 300 X01 = IV(X0) 969 V(RELDX) = SRLDST(P, D, X, V(X01)) 970 CALL SA7SST(IV, LIV, LV, V) 971 STEP1 = IV(STEP) 972 LSTGST = IV(STLSTG) 973 I = IV(RESTOR) + 1 974 GO TO (340, 310, 320, 330), I 975 310 CALL SV7CPY(P, X, V(X01)) 976 GO TO 340 977 320 CALL SV7CPY(P, V(LSTGST), V(STEP1)) 978 GO TO 340 979 330 CALL SV7CPY(P, V(STEP1), V(LSTGST)) 980 CALL SV2AXY(P, X, ONE, V(STEP1), V(X01)) 981 V(RELDX) = SRLDST(P, D, X, V(X01)) 982 IV(RESTOR) = RSTRST 983C 984C *** IF NECESSARY, SWITCH MODELS *** 985C 986 340 IF (IV(SWITCH) .EQ. 0) GO TO 350 987 IV(H) = -abs(IV(H)) 988 IV(SUSED) = IV(SUSED) + 2 989 L = IV(VSAVE) 990 CALL SV7CPY(NVSAVE, V, V(L)) 991 350 L = IV(IRC) - 4 992 STPMOD = IV(MODEL) 993 IF (L .GT. 0) GO TO (370,380,390,390,390,390,390,390,500,440), L 994C 995C *** DECIDE WHETHER TO CHANGE MODELS *** 996C 997 E = V(PREDUC) - V(FDIF) 998 S1 = IV(S) 999 CALL SS7LVM(PS, YY, V(S1), V(STEP1)) 1000 STTSST = HALF * SD7TPR(PS, V(STEP1), YY) 1001 IF (IV(MODEL) .EQ. 1) STTSST = -STTSST 1002 IF (abs(E + STTSST) * V(FUZZ) .GE. abs(E)) GO TO 360 1003C 1004C *** SWITCH MODELS *** 1005C 1006 IV(MODEL) = 3 - IV(MODEL) 1007 IF (-2 .LT. L) GO TO 400 1008 IV(H) = -abs(IV(H)) 1009 IV(SUSED) = IV(SUSED) + 2 1010 L = IV(VSAVE) 1011 CALL SV7CPY(NVSAVE, V(L), V) 1012 GO TO 160 1013C 1014 360 IF (-3 .LT. L) GO TO 400 1015C 1016C *** RECOMPUTE STEP WITH NEW RADIUS *** 1017C 1018 370 V(RADIUS) = V(RADFAC) * V(DSTNRM) 1019 GO TO 160 1020C 1021C *** COMPUTE STEP OF LENGTH V(LMAXS) FOR SINGULAR CONVERGENCE TEST 1022C 1023 380 V(RADIUS) = V(LMAXS) 1024 GO TO 200 1025C 1026C *** CONVERGENCE OR FALSE CONVERGENCE *** 1027C 1028 390 IV(CNVCOD) = L 1029 IF (V(F) .GE. V(F0)) GO TO 510 1030 IF (IV(XIRC) .EQ. 14) GO TO 510 1031 IV(XIRC) = 14 1032C 1033C. . . . . . . . . . . . PROCESS ACCEPTABLE STEP . . . . . . . . . . . 1034C 1035 400 IV(COVMAT) = 0 1036 IV(REGD) = 0 1037C 1038C *** SEE WHETHER TO SET V(RADFAC) BY GRADIENT TESTS *** 1039C 1040 IF (IV(IRC) .NE. 3) GO TO 430 1041 STEP1 = IV(STEP) 1042 TEMP1 = IV(STLSTG) 1043 TEMP2 = IV(W) 1044C 1045C *** SET TEMP1 = HESSIAN * STEP FOR USE IN GRADIENT TESTS *** 1046C 1047 HC1 = IV(HC) 1048 IF (HC1 .LE. 0) GO TO 410 1049 CALL SS7LVM(P, V(TEMP1), V(HC1), V(STEP1)) 1050 GO TO 420 1051 410 RMAT1 = IV(RMAT) 1052 CALL SL7TVM(P, V(TEMP1), V(RMAT1), V(STEP1)) 1053 CALL SL7VML(P, V(TEMP1), V(RMAT1), V(TEMP1)) 1054C 1055 420 IF (STPMOD .EQ. 1) GO TO 430 1056 S1 = IV(S) 1057 CALL SS7LVM(PS, V(TEMP2), V(S1), V(STEP1)) 1058 CALL SV2AXY(PS, V(TEMP1), ONE, V(TEMP2), V(TEMP1)) 1059C 1060C *** SAVE OLD GRADIENT AND COMPUTE NEW ONE *** 1061C 1062 430 IV(NGCALL) = IV(NGCALL) + 1 1063 G01 = IV(W) 1064 CALL SV7CPY(P, V(G01), GG) 1065 IV(1) = 2 1066 IV(TOOBIG) = 0 1067 GO TO 999 1068C 1069C *** INITIALIZATIONS -- G0 = GG - G0, ETC. *** 1070C 1071 440 G01 = IV(W) 1072 CALL SV2AXY(P, V(G01), NEGONE, V(G01), GG) 1073 STEP1 = IV(STEP) 1074 TEMP1 = IV(STLSTG) 1075 TEMP2 = IV(W) 1076 IF (IV(IRC) .NE. 3) GO TO 470 1077C 1078C *** SET V(RADFAC) BY GRADIENT TESTS *** 1079C 1080C *** SET TEMP1 = D**-1 * (HESSIAN * STEP + (G(X0) - G(X))) *** 1081C 1082 K = TEMP1 1083 L = G01 1084 DO 450 I = 1, P 1085 V(K) = (V(K) - V(L)) / D(I) 1086 K = K + 1 1087 L = L + 1 1088 450 CONTINUE 1089C 1090C *** DO GRADIENT TESTS *** 1091C 1092 IF (SV2NRM(P, V(TEMP1)) .LE. V(DGNORM) * V(TUNER4)) GO TO 460 1093 IF (SD7TPR(P, GG, V(STEP1)) 1094 1 .GE. V(GTSTEP) * V(TUNER5)) GO TO 470 1095 460 V(RADFAC) = V(INCFAC) 1096C 1097C *** COMPUTE YY VECTOR NEEDED FOR UPDATING S *** 1098C 1099 470 CALL SV2AXY(PS, YY, NEGONE, YY, GG) 1100C 1101C *** DETERMINE SIZING FACTOR V(SIZE) *** 1102C 1103C *** SET TEMP1 = S * STEP *** 1104 S1 = IV(S) 1105 CALL SS7LVM(PS, V(TEMP1), V(S1), V(STEP1)) 1106C 1107 T1 = abs(SD7TPR(PS, V(STEP1), V(TEMP1))) 1108 T = abs(SD7TPR(PS, V(STEP1), YY)) 1109 V(SIZE) = ONE 1110 IF (T .LT. T1) V(SIZE) = T / T1 1111C 1112C *** SET G0 TO WCHMTD CHOICE OF FLETCHER AND AL-BAALI *** 1113C 1114 HC1 = IV(HC) 1115 IF (HC1 .LE. 0) GO TO 480 1116 CALL SS7LVM(PS, V(G01), V(HC1), V(STEP1)) 1117 GO TO 490 1118C 1119 480 RMAT1 = IV(RMAT) 1120 CALL SL7TVM(PS, V(G01), V(RMAT1), V(STEP1)) 1121 CALL SL7VML(PS, V(G01), V(RMAT1), V(G01)) 1122C 1123 490 CALL SV2AXY(PS, V(G01), ONE, YY, V(G01)) 1124C 1125C *** UPDATE S *** 1126C 1127 CALL SS7LUP(V(S1), V(COSMIN), PS, V(SIZE), V(STEP1), V(TEMP1), 1128 1 V(TEMP2), V(G01), V(WSCALE), YY) 1129 IV(1) = 2 1130 GO TO 110 1131C 1132C. . . . . . . . . . . . . . MISC. DETAILS . . . . . . . . . . . . . . 1133C 1134C *** BAD PARAMETERS TO ASSESS *** 1135C 1136 500 IV(1) = 64 1137 GO TO 999 1138C 1139C 1140C *** CONVERGENCE OBTAINED -- SEE WHETHER TO COMPUTE COVARIANCE *** 1141C 1142 510 IF (IV(RDREQ) .EQ. 0) GO TO 600 1143 IF (IV(FDH) .NE. 0) GO TO 600 1144 IF (IV(CNVCOD) .GE. 7) GO TO 600 1145 IF (IV(REGD) .GT. 0) GO TO 600 1146 IF (IV(COVMAT) .GT. 0) GO TO 600 1147 IF (abs(IV(COVREQ)) .GE. 3) GO TO 560 1148 IF (IV(RESTOR) .EQ. 0) IV(RESTOR) = 2 1149 GO TO 530 1150C 1151C *** COMPUTE FINITE-DIFFERENCE HESSIAN FOR COMPUTING COVARIANCE *** 1152C 1153 520 IV(RESTOR) = 0 1154 530 CALL SF7HES(D, GG, I, IV, LIV, LV, P, V, X) 1155 GO TO (540, 550, 580), I 1156 540 IV(NFCOV) = IV(NFCOV) + 1 1157 IV(NFCALL) = IV(NFCALL) + 1 1158 IV(1) = 1 1159 GO TO 999 1160C 1161 550 IV(NGCOV) = IV(NGCOV) + 1 1162 IV(NGCALL) = IV(NGCALL) + 1 1163 IV(NFGCAL) = IV(NFCALL) + IV(NGCOV) 1164 IV(1) = 2 1165 GO TO 999 1166C 1167 560 H1 = abs(IV(H)) 1168 IV(H) = -H1 1169 PP1O2 = P * (P + 1) / 2 1170 RMAT1 = IV(RMAT) 1171 IF (RMAT1 .LE. 0) GO TO 570 1172 LMAT1 = IV(LMAT) 1173 CALL SV7CPY(PP1O2, V(LMAT1), V(RMAT1)) 1174 V(RCOND) = ZERO 1175 GO TO 590 1176 570 HC1 = IV(HC) 1177 IV(FDH) = H1 1178 CALL SV7CPY(P*(P+1)/2, V(H1), V(HC1)) 1179C 1180C *** COMPUTE CHOLESKY FACTOR OF FINITE-DIFFERENCE HESSIAN 1181C *** FOR USE IN CALLER*S COVARIANCE CALCULATION... 1182C 1183 580 LMAT1 = IV(LMAT) 1184 H1 = IV(FDH) 1185 IF (H1 .LE. 0) GO TO 600 1186 IF (IV(CNVCOD) .EQ. 70) GO TO 80 1187 CALL SL7SRT(1, P, V(LMAT1), V(H1), I) 1188 IV(FDH) = -1 1189 V(RCOND) = ZERO 1190 IF (I .NE. 0) GO TO 600 1191C 1192 590 IV(FDH) = -1 1193 STEP1 = IV(STEP) 1194 T = SL7SVN(P, V(LMAT1), V(STEP1), V(STEP1)) 1195 IF (T .LE. ZERO) GO TO 600 1196 TP = SL7SVX(P, V(LMAT1), V(STEP1), V(STEP1)) 1197 IF (TP .NE. ZERO) then 1198 T = T / TP 1199 IF (T .GT. SR7MDC(4)) IV(FDH) = H1 1200 V(RCOND) = T 1201 END IF 1202C 1203 600 IV(MODE) = 0 1204 IV(1) = IV(CNVCOD) 1205 IV(CNVCOD) = 0 1206 GO TO 999 1207C 1208C *** SPECIAL RETURN FOR MISSING HESSIAN INFORMATION -- BOTH 1209C *** IV(HC) .LE. 0 AND IV(RMAT) .LE. 0 1210C 1211 610 IV(1) = 1400 1212C 1213 999 RETURN 1214C 1215C *** LAST LINE OF SG7LIT FOLLOWS *** 1216 END 1217c ================================================================== 1218 SUBROUTINE SN2LRD(DR, IV, L, LH, LIV, LV, ND, NN, P, R, RD, V) 1219C 1220C *** COMPUTE REGRESSION DIAGNOSTIC AND DEFAULT COVARIANCE MATRIX FOR 1221C SRN2G *** 1222C 1223C *** PARAMETERS *** 1224C 1225 INTEGER LH, LIV, LV, ND, NN, P 1226 INTEGER IV(LIV) 1227 REAL DR(ND,P), L(LH), R(NN), RD(NN), V(LV) 1228C 1229C *** CODED BY DAVID M. GAY (WINTER 1982, FALL 1983) *** 1230C 1231C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** 1232C 1233 EXTERNAL SD7TPR, SL7ITV, SL7IVM,SO7PRD, SV7SCP 1234 REAL SD7TPR 1235C 1236C *** LOCAL VARIABLES *** 1237C 1238 INTEGER COV, I, J, M, STEP1 1239 REAL A, S, T 1240C 1241C *** CONSTANTS *** 1242C 1243 REAL NEGONE, ONE, ONEV(1), ZERO 1244C 1245C *** IV SUBSCRIPTS *** 1246C 1247 INTEGER D, H, MODE, RDREQ, STEP 1248 PARAMETER (D=27, H=56, MODE=35, RDREQ=57, STEP=40) 1249 PARAMETER (NEGONE=-1.E+0, ONE=1.E+0, ZERO=0.E+0) 1250 DATA ONEV(1)/1.E+0/ 1251C 1252C +++++++++++++++++++++++++++++++ BODY +++++++++++++++++++++++++++++++ 1253C 1254 STEP1 = IV(STEP) 1255 I = IV(RDREQ) 1256 IF (I .LE. 0) GO TO 999 1257 IF (MOD(I,4) .LT. 2) GO TO 30 1258 CALL SV7SCP(NN, RD, NEGONE) 1259 DO 20 I = 1, NN 1260 A = R(I)**2 1261 M = STEP1 1262 DO 10 J = 1, P 1263 V(M) = DR(I,J) 1264 M = M + 1 1265 10 CONTINUE 1266 CALL SL7IVM(P, V(STEP1), L, V(STEP1)) 1267 S = SD7TPR(P, V(STEP1), V(STEP1)) 1268 T = ONE - S 1269 IF (T .LE. ZERO) GO TO 20 1270 A = A * S / T 1271 RD(I) = sqrt(A) 1272 20 CONTINUE 1273C 1274 30 IF (IV(MODE) - P .LT. 2) GO TO 999 1275C 1276C *** COMPUTE DEFAULT COVARIANCE MATRIX *** 1277C 1278 COV = abs(IV(H)) 1279 DO 50 I = 1, NN 1280 M = STEP1 1281 DO 40 J = 1, P 1282 V(M) = DR(I,J) 1283 M = M + 1 1284 40 CONTINUE 1285 CALL SL7IVM(P, V(STEP1), L, V(STEP1)) 1286 CALL SL7ITV(P, V(STEP1), L, V(STEP1)) 1287 CALL SO7PRD(1, LH, P, V(COV), ONEV, V(STEP1), V(STEP1)) 1288 50 CONTINUE 1289C 1290 999 RETURN 1291C *** LAST CARD OF SN2LRD FOLLOWS *** 1292 END 1293 SUBROUTINE SC7VFN(IV, L, LH, LIV, LV, N, P, V) 1294C 1295C *** FINISH COVARIANCE COMPUTATION FOR SRN2G, SRNSG *** 1296C 1297 INTEGER LH, LIV, LV, N, P 1298 INTEGER IV(LIV) 1299 REAL L(LH), V(LV) 1300C 1301 EXTERNAL SL7NVR, SL7TSQ, SV7SCL 1302C 1303C *** LOCAL VARIABLES *** 1304C 1305 INTEGER COV, I 1306 REAL HALF 1307C 1308C *** SUBSCRIPTS FOR IV AND V *** 1309C 1310 INTEGER CNVCOD, COVMAT, F, FDH, H, MODE, RDREQ, REGD 1311C 1312 PARAMETER (CNVCOD=55, COVMAT=26, F=10, FDH=74, H=56, MODE=35, 1313 1 RDREQ=57, REGD=67) 1314 DATA HALF/0.5E+0/ 1315C 1316C *** BODY *** 1317C 1318 IV(1) = IV(CNVCOD) 1319 I = IV(MODE) - P 1320 IV(MODE) = 0 1321 IV(CNVCOD) = 0 1322 IF (IV(FDH) .LE. 0) GO TO 999 1323 IF ((I-2)**2 .EQ. 1) IV(REGD) = 1 1324 IF (MOD(IV(RDREQ),2) .NE. 1) GO TO 999 1325C 1326C *** FINISH COMPUTING COVARIANCE MATRIX = INVERSE OF F.D. HESSIAN. 1327C 1328 COV = abs(IV(H)) 1329 IV(FDH) = 0 1330C 1331 IF (IV(COVMAT) .NE. 0) GO TO 999 1332 IF (I .GE. 2) GO TO 10 1333 CALL SL7NVR(P, V(COV), L) 1334 CALL SL7TSQ(P, V(COV), V(COV)) 1335C 1336 10 CALL SV7SCL(LH, V(COV), V(F)/(HALF * real(max(1,N-P))), V(COV)) 1337 IV(COVMAT) = COV 1338C 1339 999 RETURN 1340C *** LAST LINE OF SC7VFN FOLLOWS *** 1341 END 1342 SUBROUTINE SF7HES(D, GG, IRT, IV, LIV, LV, P, V, X) 1343C 1344C *** COMPUTE FINITE-DIFFERENCE HESSIAN, STORE IT IN V STARTING 1345C *** AT V(IV(FDH)) = V(-IV(H)). 1346C 1347C *** IF IV(COVREQ) .GE. 0 THEN SF7HES USES GRADIENT DIFFERENCES, 1348C *** OTHERWISE FUNCTION DIFFERENCES. STORAGE IN V IS AS IN SG7LIT. 1349C 1350C IRT VALUES... 1351C 1 = COMPUTE FUNCTION VALUE, I.E., V(F). 1352C 2 = COMPUTE G. 1353C 3 = DONE. 1354C 1355C 1356C *** PARAMETER DECLARATIONS *** 1357C 1358 INTEGER IRT, LIV, LV, P 1359 INTEGER IV(LIV) 1360 REAL D(P), GG(P), V(LV), X(P) 1361C 1362C *** LOCAL VARIABLES *** 1363C 1364 INTEGER GSAVE1, HES, HMI, HPI, HPM, I, K, KIND, L, M, MM1, MM1O2, 1365 1 PP1O2, STPI, STPM, STP0 1366 REAL DEL, HALF, NEGPT5, ONE, TWO, ZERO 1367C 1368C *** EXTERNAL SUBROUTINES *** 1369C 1370 EXTERNAL SV7CPY 1371C 1372C SV7CPY.... COPY ONE VECTOR TO ANOTHER. 1373C 1374C *** SUBSCRIPTS FOR IV AND V *** 1375C 1376 INTEGER COVREQ, DELTA, DELTA0, DLTFDC, F, FDH, FX, H, KAGQT, MODE, 1377 1 NFGCAL, SAVEI, SWITCH, TOOBIG, W, XMSAVE 1378 PARAMETER (HALF=0.5E+0, NEGPT5=-0.5E+0, ONE=1.E+0, TWO=2.E+0, 1379 1 ZERO=0.E+0) 1380 PARAMETER (COVREQ=15, DELTA=52, DELTA0=44, DLTFDC=42, F=10, 1381 1 FDH=74, FX=53, H=56, KAGQT=33, MODE=35, NFGCAL=7, 1382 2 SAVEI=63, SWITCH=12, TOOBIG=2, W=65, XMSAVE=51) 1383C 1384C ++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ 1385C 1386 IRT = 4 1387 KIND = IV(COVREQ) 1388 M = IV(MODE) 1389 IF (M .GT. 0) GO TO 10 1390 IV(H) = -abs(IV(H)) 1391 IV(FDH) = 0 1392 IV(KAGQT) = -1 1393 V(FX) = V(F) 1394 10 IF (M .GT. P) GO TO 999 1395 IF (KIND .LT. 0) GO TO 110 1396C 1397C *** COMPUTE FINITE-DIFFERENCE HESSIAN USING BOTH FUNCTION AND 1398C *** GRADIENT VALUES. 1399C 1400 GSAVE1 = IV(W) + P 1401 IF (M .GT. 0) GO TO 20 1402C *** FIRST CALL ON SF7HES. SET GSAVE = G, TAKE FIRST STEP *** 1403 CALL SV7CPY(P, V(GSAVE1), GG) 1404 IV(SWITCH) = IV(NFGCAL) 1405 GO TO 90 1406C 1407 20 DEL = V(DELTA) 1408 X(M) = V(XMSAVE) 1409 IF (IV(TOOBIG) .EQ. 0) GO TO 40 1410C 1411C *** HANDLE OVERSIZE V(DELTA) *** 1412C 1413 IF (DEL*X(M) .GT. ZERO) GO TO 30 1414C *** WE ALREADY TRIED SHRINKING V(DELTA), SO QUIT *** 1415 IV(FDH) = -2 1416 GO TO 220 1417C 1418C *** TRY SHRINKING V(DELTA) *** 1419 30 DEL = NEGPT5 * DEL 1420 GO TO 100 1421C 1422 40 HES = -IV(H) 1423C 1424C *** SET GG = (GG - GSAVE)/DEL *** 1425C 1426 DO 50 I = 1, P 1427 GG(I) = (GG(I) - V(GSAVE1)) / DEL 1428 GSAVE1 = GSAVE1 + 1 1429 50 CONTINUE 1430C 1431C *** ADD GG AS NEW COL. TO FINITE-DIFF. HESSIAN MATRIX *** 1432C 1433 K = HES + M*(M-1)/2 1434 L = K + M - 2 1435 IF (M .EQ. 1) GO TO 70 1436C 1437C *** SET H(I,M) = 0.5 * (H(I,M) + GG(I)) FOR I = 1 TO M-1 *** 1438C 1439 MM1 = M - 1 1440 DO 60 I = 1, MM1 1441 V(K) = HALF * (V(K) + GG(I)) 1442 K = K + 1 1443 60 CONTINUE 1444C 1445C *** ADD H(I,M) = GG(I) FOR I = M TO P *** 1446C 1447 70 L = L + 1 1448 DO 80 I = M, P 1449 V(L) = GG(I) 1450 L = L + I 1451 80 CONTINUE 1452C 1453 90 M = M + 1 1454 IV(MODE) = M 1455 IF (M .GT. P) GO TO 210 1456C 1457C *** CHOOSE NEXT FINITE-DIFFERENCE STEP, RETURN TO GET GG THERE *** 1458C 1459 DEL = V(DELTA0) * max(ONE/D(M), abs(X(M))) 1460 IF (X(M) .LT. ZERO) DEL = -DEL 1461 V(XMSAVE) = X(M) 1462 100 X(M) = X(M) + DEL 1463 V(DELTA) = DEL 1464 IRT = 2 1465 GO TO 999 1466C 1467C *** COMPUTE FINITE-DIFFERENCE HESSIAN USING FUNCTION VALUES ONLY. 1468C 1469 110 STP0 = IV(W) + P - 1 1470 MM1 = M - 1 1471 MM1O2 = M*MM1/2 1472 IF (M .GT. 0) GO TO 120 1473C *** FIRST CALL ON SF7HES. *** 1474 IV(SAVEI) = 0 1475 GO TO 200 1476C 1477 120 I = IV(SAVEI) 1478 HES = -IV(H) 1479 IF (I .GT. 0) GO TO 180 1480 IF (IV(TOOBIG) .EQ. 0) GO TO 140 1481C 1482C *** HANDLE OVERSIZE STEP *** 1483C 1484 STPM = STP0 + M 1485 DEL = V(STPM) 1486 IF (DEL*X(XMSAVE) .GT. ZERO) GO TO 130 1487C *** WE ALREADY TRIED SHRINKING THE STEP, SO QUIT *** 1488 IV(FDH) = -2 1489 GO TO 220 1490C 1491C *** TRY SHRINKING THE STEP *** 1492 130 DEL = NEGPT5 * DEL 1493 X(M) = X(XMSAVE) + DEL 1494 V(STPM) = DEL 1495 IRT = 1 1496 GO TO 999 1497C 1498C *** SAVE F(X + STP(M)*E(M)) IN H(P,M) *** 1499C 1500 140 PP1O2 = P * (P-1) / 2 1501 HPM = HES + PP1O2 + MM1 1502 V(HPM) = V(F) 1503C 1504C *** START COMPUTING ROW M OF THE FINITE-DIFFERENCE HESSIAN H. *** 1505C 1506 HMI = HES + MM1O2 1507 IF (MM1 .EQ. 0) GO TO 160 1508 HPI = HES + PP1O2 1509 DO 150 I = 1, MM1 1510 V(HMI) = V(FX) - (V(F) + V(HPI)) 1511 HMI = HMI + 1 1512 HPI = HPI + 1 1513 150 CONTINUE 1514 160 V(HMI) = V(F) - TWO*V(FX) 1515C 1516C *** COMPUTE FUNCTION VALUES NEEDED TO COMPLETE ROW M OF H. *** 1517C 1518 I = 1 1519C 1520 170 IV(SAVEI) = I 1521 STPI = STP0 + I 1522 V(DELTA) = X(I) 1523 X(I) = X(I) + V(STPI) 1524 IF (I .EQ. M) X(I) = V(XMSAVE) - V(STPI) 1525 IRT = 1 1526 GO TO 999 1527C 1528 180 X(I) = V(DELTA) 1529 IF (IV(TOOBIG) .EQ. 0) GO TO 190 1530C *** PUNT IN THE EVENT OF AN OVERSIZE STEP *** 1531 IV(FDH) = -2 1532 GO TO 220 1533C 1534C *** FINISH COMPUTING H(M,I) *** 1535C 1536 190 STPI = STP0 + I 1537 HMI = HES + MM1O2 + I - 1 1538 STPM = STP0 + M 1539 V(HMI) = (V(HMI) + V(F)) / (V(STPI)*V(STPM)) 1540 I = I + 1 1541 IF (I .LE. M) GO TO 170 1542 IV(SAVEI) = 0 1543 X(M) = V(XMSAVE) 1544C 1545 200 M = M + 1 1546 IV(MODE) = M 1547 IF (M .GT. P) GO TO 210 1548C 1549C *** PREPARE TO COMPUTE ROW M OF THE FINITE-DIFFERENCE HESSIAN H. 1550C *** COMPUTE M-TH STEP SIZE STP(M), THEN RETURN TO OBTAIN 1551C *** F(X + STP(M)*E(M)), WHERE E(M) = M-TH STD. UNIT VECTOR. 1552C 1553 DEL = V(DLTFDC) * max(ONE/D(M), abs(X(M))) 1554 IF (X(M) .LT. ZERO) DEL = -DEL 1555 V(XMSAVE) = X(M) 1556 X(M) = X(M) + DEL 1557 STPM = STP0 + M 1558 V(STPM) = DEL 1559 IRT = 1 1560 GO TO 999 1561C 1562C *** RESTORE V(F), ETC. *** 1563C 1564 210 IV(FDH) = HES 1565 220 V(F) = V(FX) 1566 IRT = 3 1567 IF (KIND .LT. 0) GO TO 999 1568 IV(NFGCAL) = IV(SWITCH) 1569 GSAVE1 = IV(W) + P 1570 CALL SV7CPY(P, GG, V(GSAVE1)) 1571 GO TO 999 1572C 1573 999 RETURN 1574C *** LAST CARD OF SF7HES FOLLOWS *** 1575 END 1576 SUBROUTINE SN2CVP(IV, LIV, LV, P, V) 1577C 1578C *** PRINT COVARIANCE MATRIX FOR SRN2G *** 1579C 6/27/90 CLL changed 'SCALE' to 'VARFAC' in output labels. 1580c ------------------------------------------------------------------ 1581c%% long int j, k; 1582 INTEGER J 1583 INTEGER LIV, LV, P 1584 INTEGER IV(LIV) 1585 REAL V(LV) 1586C 1587C *** LOCAL VARIABLES *** 1588C 1589 INTEGER COV1, I, II, I1, PU 1590 REAL T 1591C 1592C *** IV SUBSCRIPTS *** 1593C 1594 INTEGER COVMAT, COVPRT, COVREQ, NEEDHD, NFCOV, NGCOV, PRUNIT, 1595 1 RCOND, REGD, STATPR 1596C 1597 PARAMETER (COVMAT=26, COVPRT=14, COVREQ=15, NEEDHD=36, NFCOV=52, 1598 1 NGCOV=53, PRUNIT=21, RCOND=53, REGD=67, STATPR=23) 1599C *** BODY *** 1600C 1601 10 FORMAT(/1X,I4, 1602 * ' EXTRA FUNC. EVALS FOR COVARIANCE AND DIAGNOSTICS.') 1603 20 FORMAT(1X,I4, 1604 * ' EXTRA GRAD. EVALS FOR COVARIANCE AND DIAGNOSTICS.') 1605 40 FORMAT(/' RECIPROCAL CONDITION OF F.D. HESSIAN = AT MOST',g10.2) 1606 60 FORMAT(/' RECIPROCAL CONDITION OF (J**T)*J = AT LEAST',g10.2) 1607 90 FORMAT(/' ++++++ INDEFINITE COVARIANCE MATRIX ++++++') 1608 100 FORMAT(/' ++++++ OVERSIZE STEPS IN COMPUTING COVARIANCE +++++') 1609 120 FORMAT(/' ++++++ COVARIANCE MATRIX NOT COMPUTED ++++++') 1610c 1611c++(~.C.) Default UNITNO='(PU,' 1612c++(.C.) Default UNITNO='(*,' 1613c++ Replace "(PU," = UNITNO 1614c 1615 IF (IV(1) .GT. 8) GO TO 999 1616 PU = IV(PRUNIT) 1617 IF (PU .EQ. 0) GO TO 999 1618 COV1 = IV(COVMAT) 1619 IF (-2 .EQ. COV1) WRITE(PU,100) 1620 IF (IV(STATPR) .EQ. 0) GO TO 30 1621 IF (IV(NFCOV) .GT. 0) WRITE(PU,10) IV(NFCOV) 1622 IF (IV(NGCOV) .GT. 0) WRITE(PU,20) IV(NGCOV) 1623C 1624 30 IF (IV(COVPRT) .LE. 0) GO TO 999 1625 IF (IV(REGD) .LE. 0 .AND. COV1 .LE. 0) GO TO 70 1626 IV(NEEDHD) = 1 1627 T = V(RCOND)**2 1628 IF (abs(IV(COVREQ)) .GT. 2) GO TO 50 1629C 1630 WRITE(PU,40) T 1631 GO TO 70 1632C 1633 50 WRITE(PU,60) T 1634C 1635 70 IF (MOD(IV(COVPRT),2) .EQ. 0) GO TO 999 1636 IV(NEEDHD) = 1 1637 IF (COV1) 80,110,130 1638 80 IF (-1 .EQ. COV1) WRITE(PU,90) 1639 GO TO 999 1640C 1641 110 WRITE(PU,120) 1642 GO TO 999 1643C 1644 130 I = abs(IV(COVREQ)) 1645 IF (I .LE. 1) WRITE(PU,'(/ 1646 * '' COVARIANCE = VARFAC * H**-1 * (J**T * J) * H**-1''/ 1647 * '' WHERE H = F.D. HESSIAN''/1x)') 1648 IF (I .EQ. 2) WRITE(PU,'(1x/'' COVARIANCE = VARFAC * H**-1,'', 1649 * '' WHERE H = FINITE-DIFFERENCE HESSIAN''/1x)') 1650 IF (I.GT.2) WRITE(PU,'(/'' COVARIANCE = VARFAC * J**T * J''/1x)') 1651 II = COV1 - 1 1652 DO 170 I = 1, P 1653 I1 = II + 1 1654 II = II + I 1655C%% printf( " ROW%3ld ", i ); 1656C%% for (j = i1; j <= ii; j+=5){ 1657C%% for (k = j; k <= (j <= ii - 5 ? j+4 : ii); k++) 1658C%% printf( "%12.3g", V[k] ); 1659C%% printf( "\n"); 1660C%% if (j <= ii - 5) printf( " ");} 1661 WRITE(PU,'('' ROW'',I3,2X,5g12.3/(9X,5g12.3))')I,(V(J),J=I1,II) 1662 170 CONTINUE 1663C 1664 999 RETURN 1665C *** LAST CARD OF SN2CVP FOLLOWS *** 1666 END 1667 SUBROUTINE SN2RDP(IV, LIV, N, RD) 1668C 1669C *** PRINT REGRESSION DIAGNOSTICS FOR MLPSL AND NL2S1 *** 1670C 1671c ------------------------------------------------------------------ 1672c++ Code for .C. is inactive 1673c%% long int j,k; 1674c++ End 1675 INTEGER LIV, N 1676 INTEGER IV(LIV) 1677 REAL RD(N) 1678 INTEGER PU 1679C 1680C *** IV SUBSCRIPTS *** 1681C 1682 INTEGER COVPRT, NEEDHD, PRUNIT, RDREQ, REGD 1683C 1684C DATA COVPRT/14/, NEEDHD/36/, PRUNIT/21/, RDREQ/57/, REGD/67/ 1685 PARAMETER (COVPRT=14, NEEDHD=36, PRUNIT=21, RDREQ=57, REGD=67) 1686C 1687C ++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ 1688C 1689 PU = IV(PRUNIT) 1690 IF (PU .EQ. 0) GO TO 999 1691 IF (IV(COVPRT) .LT. 2) GO TO 999 1692 IF (IV(REGD) .LE. 0) GO TO 999 1693 IV(NEEDHD) = 1 1694 WRITE(PU,'('' REGRESSION DIAGNOSTIC = SQRT(G(I)**T * H(I)**-1 *'', 1695 * ''G(I))...''/)') 1696C%% for(j=0; j < n; j+=6){ 1697C%% for (k = j; k < (j < n - 6 ? j+6 : n); k++) 1698C%% printf( "%12.3g", rd[k] ); 1699C%% printf( "\n" );} 1700 WRITE(PU,'(6g12.3)') RD 1701C 1702 999 RETURN 1703C *** LAST CARD OF SN2RDP FOLLOWS *** 1704 END 1705 SUBROUTINE SO7PRD(L, LS, PP, SS, WW, YY, Z) 1706C 1707C *** FOR I = 1..L, SET SS = SS + WW(I)*YY(.,I)*(Z(.,I)**T), I.E., 1708C *** ADD WW(I) TIMES THE OUTER PRODUCT OF Y(.,I) AND Z(.,I). 1709C 1710c ------------------------------------------------------------------ 1711 INTEGER L, LS, PP 1712 REAL SS(LS), WW(L), YY(PP,L), Z(PP,L) 1713C DIMENSION SS(PP*(PP+1)/2) 1714C 1715 INTEGER I, J, K, M 1716 REAL WK, YI, ZERO 1717 parameter(ZERO = 0.0E0) 1718C 1719 DO 30 K = 1, L 1720 WK = WW(K) 1721 IF (WK .EQ. ZERO) GO TO 30 1722 M = 1 1723 DO 20 I = 1, PP 1724 YI = WK * YY(I,K) 1725 DO 10 J = 1, I 1726 SS(M) = SS(M) + YI*Z(J,K) 1727 M = M + 1 1728 10 CONTINUE 1729 20 CONTINUE 1730 30 CONTINUE 1731C 1732 RETURN 1733C *** LAST CARD OF SO7PRD FOLLOWS *** 1734 END 1735 SUBROUTINE SL7NVR(N, LIN, L) 1736C 1737C *** COMPUTE LIN = L**-1, BOTH N X N LOWER TRIANG. STORED *** 1738C *** COMPACTLY BY ROWS. LIN AND L MAY SHARE THE SAME STORAGE. *** 1739C 1740C *** PARAMETERS *** 1741C 1742c ------------------------------------------------------------------ 1743 INTEGER N 1744 REAL L(*), LIN(*) 1745C DIMENSION L(N*(N+1)/2), LIN(N*(N+1)/2) 1746C 1747C *** LOCAL VARIABLES *** 1748C 1749 INTEGER I, II, IM1, JJ, J0, J1, K, K0, NP1 1750 REAL ONE, T, ZERO 1751 PARAMETER (ONE=1.E+0, ZERO=0.E+0) 1752C 1753C *** BODY *** 1754C 1755 NP1 = N + 1 1756 J0 = N*(NP1)/2 1757 DO 30 II = 1, N 1758 I = NP1 - II 1759 LIN(J0) = ONE/L(J0) 1760 IF (I .LE. 1) GO TO 999 1761 J1 = J0 1762 IM1 = I - 1 1763 DO 20 JJ = 1, IM1 1764 T = ZERO 1765 J0 = J1 1766 K0 = J1 - JJ 1767 DO 10 K = 1, JJ 1768 T = T - L(K0)*LIN(J0) 1769 J0 = J0 - 1 1770 K0 = K0 + K - I 1771 10 CONTINUE 1772 LIN(J0) = T/L(K0) 1773 20 CONTINUE 1774 J0 = J0 - 1 1775 30 CONTINUE 1776 999 RETURN 1777C *** LAST CARD OF SL7NVR FOLLOWS *** 1778 END 1779 SUBROUTINE SL7TSQ(N, A, L) 1780C 1781C *** SET A TO LOWER TRIANGLE OF (L**T) * L *** 1782C 1783C *** L = N X N LOWER TRIANG. MATRIX STORED ROWWISE. *** 1784C *** A IS ALSO STORED ROWWISE AND MAY SHARE STORAGE WITH L. *** 1785C 1786c ------------------------------------------------------------------ 1787 INTEGER N 1788 REAL A(*), L(*) 1789C DIMENSION A(N*(N+1)/2), L(N*(N+1)/2) 1790C 1791 INTEGER I, II, IIM1, I1, J, K, M 1792 REAL LII, LJ 1793C 1794 II = 0 1795 DO 50 I = 1, N 1796 I1 = II + 1 1797 II = II + I 1798 M = 1 1799 IF (I .EQ. 1) GO TO 30 1800 IIM1 = II - 1 1801 DO 20 J = I1, IIM1 1802 LJ = L(J) 1803 DO 10 K = I1, J 1804 A(M) = A(M) + LJ*L(K) 1805 M = M + 1 1806 10 CONTINUE 1807 20 CONTINUE 1808 30 LII = L(II) 1809 DO 40 J = I1, II 1810 40 A(J) = LII * L(J) 1811 50 CONTINUE 1812C 1813 RETURN 1814C *** LAST CARD OF SL7TSQ FOLLOWS *** 1815 END 1816