1 SUBROUTINE SNLSFB(N, P, L, ALF, B, C, Y, SCALCA, INC, IINC, IV, 2 1 LIV, LV, V) 3c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. 4c ALL RIGHTS RESERVED. 5c Based on Government Sponsored Research NAS7-03001. 6c>> 1996-04-27 SNLSFB Krogh Changes to get desired C prototypes. 7c>> 1994-10-20 SNLSFB Krogh Changes to use M77CON 8c>> 1990-07-02 SNLSFB CLL @ JPL 9c>> 1990-06-12 CLL @ JPL 10c>> 1990-02-14 CLL @ JPL 11*** from netlib, Fri Feb 9 13:10:09 EST 1990 *** 12c--S replaces "?": ?NLSFB,?CALCA,?IVSET,?RNSGB,?V2AXY,?V7CPY,?V7SCL 13C 14C *** SOLVE SEPARABLE NONLINEAR LEAST SQUARES USING 15C *** FINITE-DIFFERENCE DERIVATIVES. 16C 17C *** PARAMETER DECLARATIONS *** 18C 19 EXTERNAL SCALCA 20 INTEGER IINC, L, LIV, LV, N, P 21 INTEGER INC(IINC,P), IV(LIV) 22 REAL ALF(P), C(L), B(2,P), V(LV), Y(N) 23C 24C *** PARAMETERS *** 25C 26C N (IN) NUMBER OF OBSERVATIONS. 27C P (IN) NUMBER OF NONLINEAR PARAMETERS TO BE ESTIMATED. 28C L (IN) NUMBER OF LINEAR PARAMETERS TO BE ESTIMATED. 29C ALF (I/O) NONLINEAR PARAMETERS. 30C INPUT = INITIAL GUESS, 31C OUTPUT = BEST ESTIMATE FOUND. 32C B (IN) SIMBLE BOUNDS ON ALF.. B(1,I) .LE. ALF(I) .LE. B(2,I). 33C C (OUT) LINEAR PARAMETERS (ESTIMATED). 34C Y (IN) RIGHT-HAND SIDE VECTOR. 35C SCALCA (IN) SUBROUTINE TO COMPUTE A MATRIX. 36C INC (IN) INCIDENCE MATRIX OF DEPENDENCIES OF COLUMNS OF A ON 37C COMPONENTS OF ALF -- INC(I,J) = 1 MEANS COLUMN I 38C OF A DEPENDS ON ALF(J). 39C IINC (IN) DECLARED LEAD DIMENSION OF INC. MUST BE AT LEAST L+1. 40C IV (I/O) INTEGER PARAMETER AND SCRATCH VECTOR. 41C LIV (IN) LENGTH OF IV. MUST BE AT LEAST 42C 122 + 2*M + 7*P + 2*L + MAX(L+1,6*P), WHERE M IS 43C THE NUMBER OF ONES IN INC. 44C LV (IN) LENGTH OF V. MUST BE AT LEAST 45C 105 + N*(2*L+6+P) + L*(L+3)/2 + P*(2*P + 22). 46c If row L+1 of INC() contains only zeros, meaning 47c the term PHI sub (L+1) is absent from the model, 48c then LV can be 4*N less than just described. 49C V (I/O) FLOATING-POINT PARAMETER AND SCRATCH VECTOR. 50C 51c ------------------------- DECLARATIONS ---------------------------- 52C 53C 54C *** EXTERNAL SUBROUTINES *** 55C 56 EXTERNAL SIVSET, IDSM, SRNSGB,SV2AXY,SV7CPY, SV7SCL 57C 58C SIVSET.... PROVIDES DEFAULT IV AND V VALUES. 59C IDSM...... DETERMINES EFFICIENT ORDER FOR FINITE DIFFERENCES. 60C SRNSGB... CARRIES OUT NL2SOL ALGORITHM. 61C SV2AXY.... ADDS A MULTIPLE OF ONE VECTOR TO ANOTHER. 62C SV7CPY.... COPIES ONE VECTOR TO ANOTHER. 63C SV7SCL... SCALES AND COPIES ONE VECTOR TO ANOTHER. 64C 65C *** LOCAL VARIABLES *** 66C 67 LOGICAL PARTJ 68 INTEGER A0, A1, AJ, ALP1, BWA1, D0, DA0, DA1, DAJ, GPTR1, GRP1, 69 1 GRP2, I, I1, IN0, IN1, IN2, INI, INLEN, IPNTR1, IV1, IWA1, 70 2 IWALEN, J1, JN1, JPNTR1, K, L1, LP1, M, M0, NF, NG, NGRP0, 71 3 NGRP1, NGRP2, RSAVE0, RSAVE1, RSVLEN, X0I, XSAVE0, XSAVE1 72 REAL DELTA, DI, H, XI, XI1 73 REAL NEGONE, ONE, ZERO 74C 75C *** SUBSCRIPTS FOR IV AND V *** 76C 77 INTEGER AMAT, COVREQ, D, DAMAT, DLTFDJ, GPTR, GRP, IN, IVNEED, J, 78 1 L1SAV, MAXGRP, MODE, MSAVE, NEXTIV, NEXTV, NFCALL, NFGCAL, 79 2 PERM, R, RESTOR, TOOBIG, VNEED, XSAVE 80C 81C *** IV SUBSCRIPT VALUES *** 82C 83C/6 84C DATA AMAT/113/, COVREQ/15/, D/27/, DAMAT/114/, DLTFDJ/43/, 85C 1 GPTR/117/, GRP/118/, IN/112/, IVNEED/3/, J/70/, L1SAV/111/, 86C 2 MAXGRP/116/, MODE/35/, MSAVE/115/, NEXTIV/46/, NEXTV/47/, 87C 3 NFCALL/6/, NFGCAL/7/, PERM/58/, R/61/, RESTOR/9/, TOOBIG/2/, 88C 4 VNEED/4/, XSAVE/119/ 89C/7 90 PARAMETER (AMAT=113, COVREQ=15, D=27, DAMAT=114, DLTFDJ=43, 91 1 GPTR=117, GRP=118, IN=112, IVNEED=3, J=70, L1SAV=111, 92 2 MAXGRP=116, MODE=35, MSAVE=115, NEXTIV=46, NEXTV=47, 93 3 NFCALL=6, NFGCAL=7, PERM=58, R=61, RESTOR=9, TOOBIG=2, 94 4 VNEED=4, XSAVE=119) 95C/ 96 DATA NEGONE/-1.E+0/, ONE/1.E+0/, ZERO/0.E+0/ 97C 98C++++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++ 99C 100 LP1 = L + 1 101 IF (IV(1) .EQ. 0) CALL SIVSET(1, IV, LIV, LV, V) 102 IF (P .LE. 0 .OR. L .LT. 0 .OR. IINC .LE. L) GO TO 80 103 IV1 = IV(1) 104 IF (IV1 .EQ. 14) GO TO 120 105 IF (IV1 .GT. 2 .AND. IV1 .LT. 12) GO TO 120 106 IF (IV1 .EQ. 12) IV(1) = 13 107 IF (IV(1) .NE. 13) GO TO 50 108C 109C *** FRESH START *** 110C 111 IF (IV(PERM) .LE. XSAVE) IV(PERM) = XSAVE + 1 112C 113C *** CHECK INC, COUNT ITS NONZEROS 114C 115 L1 = 0 116 M = 0 117 DO 40 I = 1, P 118 IF (B(1,I) .GE. B(2,I)) GO TO 40 119 M0 = M 120 IF (L .EQ. 0) GO TO 20 121 DO 10 K = 1, L 122 IF (INC(K,I) .LT. 0 .OR. INC(K,I) .GT. 1) GO TO 80 123 IF (INC(K,I) .EQ. 1) M = M + 1 124 10 CONTINUE 125 20 IF (INC(LP1,I) .NE. 1) GO TO 30 126 M = M + 1 127 L1 = 1 128 GO TO 40 129 30 IF (M .EQ. M0 .OR. INC(LP1,I) .LT. 0 130 1 .OR. INC(LP1,I) .GT. 1) GO TO 80 131 40 CONTINUE 132C 133C *** NOW L1 = 1 MEANS A HAS COLUMN L+1 *** 134C 135C *** COMPUTE STORAGE REQUIREMENTS *** 136C 137 IWALEN = max(LP1, 6*P) 138 INLEN = 2 * M 139 IV(IVNEED) = IV(IVNEED) + INLEN + 3*P + L + IWALEN + 3 140 RSVLEN = 2 * L1 * N 141 L1 = L + L1 142 IV(VNEED) = IV(VNEED) + 2*N*L1 + RSVLEN + P 143C 144 50 CALL SRNSGB(V, ALF, B, C, V, IV, IV, L, 1, N, LIV, LV, N, M, P, V, 145 1 Y) 146 IF (IV(1) .NE. 14) GO TO 999 147C 148C *** STORAGE ALLOCATION *** 149C 150 IV(IN) = IV(NEXTIV) 151 IV(AMAT) = IV(NEXTV) 152 IV(DAMAT) = IV(AMAT) + N*L1 153 IV(XSAVE) = IV(DAMAT) + N*L1 154 IV(NEXTV) = IV(XSAVE) + P + RSVLEN 155 IV(L1SAV) = L1 156 IV(MSAVE) = M 157C 158C *** DETERMINE HOW MANY GROUPS FOR FINITE DIFFERENCES 159C *** (SET UP TO CALL IDSM) 160C 161 IN1 = IV(IN) 162 JN1 = IN1 + M 163 DO 70 K = 1, P 164 IF (B(1,K) .GE. B(2,K)) GO TO 70 165 DO 60 I = 1, LP1 166 IF (INC(I,K) .EQ. 0) GO TO 60 167 IV(IN1) = I 168 IN1 = IN1 + 1 169 IV(JN1) = K 170 JN1 = JN1 + 1 171 60 CONTINUE 172 70 CONTINUE 173 IN1 = IV(IN) 174 JN1 = IN1 + M 175 IWA1 = IN1 + INLEN 176 NGRP1 = IWA1 + IWALEN 177 BWA1 = NGRP1 + P 178 IPNTR1 = BWA1 + P 179 JPNTR1 = IPNTR1 + L + 2 180 CALL IDSM(LP1, P, M, IV(IN1), IV(JN1), IV(NGRP1), NG, K, I, 181 1 IV(IPNTR1), IV(JPNTR1), IV(IWA1), IWALEN, IV(BWA1)) 182 IF (I .EQ. 1) GO TO 90 183 IV(1) = 69 184 GO TO 50 185 80 IV(1) = 66 186 GO TO 50 187C 188C *** SET UP GRP AND GPTR ARRAYS FOR COMPUTING FINITE DIFFERENCES 189C 190C *** THERE ARE NG GROUPS. GROUP I CONTAINS ALF(GRP(J)) FOR 191C *** GPTR(I) .LE. J .LE. GPTR(I+1)-1. 192C 193 90 IV(MAXGRP) = NG 194 IV(GPTR) = IN1 + 2*L1 195 GPTR1 = IV(GPTR) 196 IV(GRP) = GPTR1 + NG + 1 197 IV(NEXTIV) = IV(GRP) + P 198 GRP1 = IV(GRP) 199 NGRP0 = NGRP1 - 1 200 NGRP2 = NGRP0 + P 201 DO 110 I = 1, NG 202 IV(GPTR1) = GRP1 203 GPTR1 = GPTR1 + 1 204 DO 100 I1 = NGRP1, NGRP2 205 IF (IV(I1) .NE. I) GO TO 100 206 K = I1 - NGRP0 207 IF (B(1,K) .GE. B(2,K)) GO TO 100 208 IV(GRP1) = K 209 GRP1 = GRP1 + 1 210 100 CONTINUE 211 110 CONTINUE 212 IV(GPTR1) = GRP1 213 IF (IV1 .EQ. 13) GO TO 999 214C 215C *** INITIALIZE POINTERS *** 216C 217 120 A1 = IV(AMAT) 218 A0 = A1 - N 219 DA1 = IV(DAMAT) 220 DA0 = DA1 - N 221 IN1 = IV(IN) 222 IN0 = IN1 - 2 223 L1 = IV(L1SAV) 224 IN2 = IN1 + 2*L1 - 1 225 D0 = IV(D) - 1 226 NG = IV(MAXGRP) 227 XSAVE1 = IV(XSAVE) 228 XSAVE0 = XSAVE1 - 1 229 RSAVE1 = XSAVE1 + P 230 RSAVE0 = RSAVE1 + N 231 ALP1 = A1 + L*N 232 DELTA = V(DLTFDJ) 233 IV(COVREQ) = -abs(IV(COVREQ)) 234C 235 130 CALL SRNSGB(V(A1), ALF, B, C, V(DA1), IV(IN1), IV, L, L1, N, LIV, 236 1 LV, N, L1, P, V, Y) 237 IF (IV(1)-2) 140, 150, 999 238C 239C *** NEW FUNCTION VALUE (R VALUE) NEEDED *** 240C 241 140 NF = IV(NFCALL) 242C%% (*scalca)( n, p, l, alf, &nf, &V[a1] ); 243 CALL SCALCA(N, P, L, ALF, NF, V(A1)) 244 IF (NF .LE. 0) IV(TOOBIG) = 1 245 IF (L1 .LE. L) GO TO 130 246 IF (IV(RESTOR) .EQ. 2) CALL SV7CPY(N, V(RSAVE0), V(RSAVE1)) 247 CALL SV7CPY(N, V(RSAVE1), V(ALP1)) 248 GO TO 130 249C 250C *** COMPUTE DR = GRADIENT OF R COMPONENTS *** 251C 252 150 IF (L1 .GT. L .AND. IV(NFGCAL) .EQ. IV(NFCALL)) 253 1 CALL SV7CPY(N, V(RSAVE0), V(RSAVE1)) 254 GPTR1 = IV(GPTR) 255 DO 260 K = 1, NG 256 CALL SV7CPY(P, V(XSAVE1), ALF) 257 GRP1 = IV(GPTR1) 258 GRP2 = IV(GPTR1+1) - 1 259 GPTR1 = GPTR1 + 1 260 DO 180 I1 = GRP1, GRP2 261 I = IV(I1) 262 XI = ALF(I) 263 J1 = D0 + I 264 DI = V(J1) 265 IF (DI .LE. ZERO) DI = ONE 266 H = DELTA * max(abs(XI), ONE/DI) 267 IF (XI .LT. ZERO) GO TO 160 268 XI1 = XI + H 269 IF (XI1 .LE. B(2,I)) GO TO 170 270 XI1 = XI - H 271 IF (XI1 .GE. B(1,I)) GO TO 170 272 GO TO 190 273 160 XI1 = XI - H 274 IF (XI1 .GE. B(1,I)) GO TO 170 275 XI1 = XI + H 276 IF (XI1 .LE. B(2,I)) GO TO 170 277 GO TO 190 278 170 X0I = XSAVE0 + I 279 V(X0I) = XI1 280 180 CONTINUE 281C%% (*scalca)( n, p, l, &V[xsave1], &nf, &V[da1] ); 282 CALL SCALCA(N, P, L, V(XSAVE1), NF, V(DA1)) 283 IF (IV(NFGCAL) .GT. 0) GO TO 200 284 190 IV(TOOBIG) = 1 285 GO TO 130 286 200 JN1 = IN1 287 DO 210 I = IN1, IN2 288 210 IV(I) = 0 289 PARTJ = IV(MODE) .LE. P 290 DO 250 I1 = GRP1, GRP2 291 I = IV(I1) 292 DO 240 J1 = 1, L1 293 IF (INC(J1,I) .EQ. 0) GO TO 240 294 INI = IN0 + 2*J1 295 IV(INI) = I 296 IV(INI+1) = J1 297 X0I = XSAVE0 + I 298 H = ONE / (V(X0I) - ALF(I)) 299 DAJ = DA0 + J1*N 300 IF (PARTJ) GO TO 220 301C *** FULL FINITE DIFFERENCE FOR COV. AND REG. DIAG. *** 302 AJ = A0 + J1*N 303 CALL SV2AXY(N, V(DAJ), NEGONE, V(AJ), V(DAJ)) 304 GO TO 230 305 220 IF (J1 .GT. L) 306 1 CALL SV2AXY(N, V(DAJ), NEGONE, V(RSAVE0), V(DAJ)) 307 230 CALL SV7SCL(N, V(DAJ), H, V(DAJ)) 308 240 CONTINUE 309 250 CONTINUE 310 IF (K .GE. NG) GO TO 270 311 IV(1) = -2 312 CALL SRNSGB(V(A1), ALF, B, C, V(DA1), IV(IN1), IV, L, L1, N, 313 1 LIV, LV, N, L1, P, V, Y) 314 IF (-2 .NE. IV(1)) GO TO 999 315 260 CONTINUE 316 270 IV(1) = 2 317 GO TO 130 318C 319 999 RETURN 320C 321C *** LAST CARD OF SNLSFB FOLLOWS *** 322 END 323