1 SUBROUTINE SRNSGB(A, ALF, B, C, DA, IN, IV, L, L1, LA, LIV, LV, 2 1 N, NDA, P, V, Y) 3c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. 4c ALL RIGHTS RESERVED. 5c Based on Government Sponsored Research NAS7-03001. 6C>> 2000-12-03 SRNSGB Krogh Removed references to IV(RESTOR). 7C>> 1998-10-29 SRNSGB Krogh Moved external statement up for mangle. 8c>> 1994-10-20 SRNSGB Krogh Changes to use M77CON 9c>> 1990-06-12 SRNSGB CLL @ JPL 10c>> 1990-02-14 CLL @ JPL 11*** from netlib, Fri Feb 9 13:10:09 EST 1990 *** 12C 13C *** ITERATION DRIVER FOR SEPARABLE NONLINEAR LEAST SQUARES, 14C *** WITH SIMPLE BOUNDS ON THE NONLINEAR VARIABLES. 15C 16C *** PARAMETER DECLARATIONS *** 17C 18 INTEGER L, L1, LA, LIV, LV, N, NDA, P 19 INTEGER IN(2,NDA), IV(LIV) 20 REAL A(LA,L1), ALF(P), B(2,P), C(L), DA(LA,NDA), 21 1 V(LV), Y(N) 22C 23C *** PURPOSE *** 24C 25C GIVEN A SET OF N OBSERVATIONS Y(1)....Y(N) OF A DEPENDENT VARIABLE 26C T(1)...T(N), SRNSGB ATTEMPTS TO COMPUTE A LEAST SQUARES FIT 27C TO A FUNCTION ETA (THE MODEL) WHICH IS A LINEAR COMBINATION 28C 29C L 30C ETA(C,ALF,T) = SUM C * PHI(ALF,T) +PHI (ALF,T) 31C J=1 J J L+1 32C 33C OF NONLINEAR FUNCTIONS PHI(J) DEPENDENT ON T AND ALF(1),...,ALF(P) 34C (.E.G. A SUM OF EXPONENTIALS OR GAUSSIANS). THAT IS, IT DETERMINES 35C NONLINEAR PARAMETERS ALF WHICH MINIMIZE 36C 37C 2 N 2 38C NORM(RESIDUAL) = SUM (Y - ETA(C,ALF,T )) , 39C I=1 I I 40C 41C SUBJECT TO THE SIMPLE BOUND CONSTRAINTS B(1,I) .LE. X(I) .LE. B(2,I), 42C I = 1(1)P. 43C 44C THE (L+1)ST TERM IS OPTIONAL. 45C 46C 47C *** PARAMETERS *** 48C 49C A (IN) MATRIX PHI(ALF,T) OF THE MODEL. 50C ALF (I/O) NONLINEAR PARAMETERS. 51C INPUT = INITIAL GUESS, 52C OUTPUT = BEST ESTIMATE FOUND. 53C C (OUT) LINEAR PARAMETERS (ESTIMATED). 54C DA (IN) DERIVATIVES OF COLUMNS OF A WITH RESPECT TO COMPONENTS 55C OF ALF, AS SPECIFIED BY THE IN ARRAY... 56C IN (IN) WHEN SRNSGB IS CALLED WITH IV(1) = 2 OR -2, THEN FOR 57C I = 1(1)NDA, COLUMN I OF DA IS THE PARTIAL 58C DERIVATIVE WITH RESPECT TO ALF(IN(1,I)) OF COLUMN 59C IN(2,I) OF A, UNLESS IV(1,I) IS NOT POSITIVE (IN 60C WHICH CASE COLUMN I OF DA IS IGNORED. IV(1) = -2 61C MEANS THERE ARE MORE COLUMNS OF DA TO COME AND 62C SRNSGB SHOULD RETURN FOR THEM. 63C IV (I/O) INTEGER PARAMETER AND SCRATCH VECTOR. SRNSGB RETURNS 64C WITH IV(1) = 1 WHEN IT WANTS A TO BE EVALUATED AT 65C ALF AND WITH IV(1) = 2 WHEN IT WANTS DA TO BE 66C EVALUATED AT ALF. WHEN CALLED WITH IV(1) = -2 67C (AFTER A RETURN WITH IV(1) = 2), SRNSGB RETURNS 68C WITH IV(1) = -2 TO GET MORE COLUMNS OF DA. 69C L (IN) NUMBER OF LINEAR PARAMETERS TO BE ESTIMATED. 70C L1 (IN) L+1 IF PHI(L+1) IS IN THE MODEL, L IF NOT. 71C LA (IN) LEAD DIMENSION OF A. MUST BE AT LEAST N. 72C LIV (IN) LENGTH OF IV. MUST BE AT LEAST 110 + L + 4*P. 73C LV (IN) LENGTH OF V. MUST BE AT LEAST 74C 105 + 2*N + L*(L+3)/2 + P*(2*P + 21 + N). 75C N (IN) NUMBER OF OBSERVATIONS. 76C NDA (IN) NUMBER OF COLUMNS IN DA AND IN. 77C P (IN) NUMBER OF NONLINEAR PARAMETERS TO BE ESTIMATED. 78C V (I/O) FLOATING-POINT PARAMETER AND SCRATCH VECTOR. 79C Y (IN) RIGHT-HAND SIDE VECTOR. 80C 81C 82C *** EXTERNAL SUBROUTINES *** 83C 84 EXTERNAL SIVSET,SITSUM, SL7ITV, SL7SVX, SL7SVN, SRN2GB, SQ7APL, 85 1 SQ7RFH, SR7MDC, SS7CPR,SV2AXY,SV7CPY,SV7PRM, SV7SCP 86 REAL SL7SVX, SL7SVN, SR7MDC 87C 88c--S replaces "?": ?RNSGB,?ITSUM,?IVSET,?L7ITV,?L7SVN,?L7SVX,?Q7APL 89c--& ?Q7RFH,?R7MDC,?RN2GB,?S7CPR,?V2AXY,?V7CPY,?V7PRM 90c--& ?V7SCP 91C 92C SIVSET.... SUPPLIES DEFAULT PARAMETER VALUES. 93C SITSUM.... PRINTS ITERATION SUMMARY, INITIAL AND FINAL ALF. 94C SL7ITV... APPLIES INVERSE-TRANSPOSE OF COMPACT LOWER TRIANG. MATRIX. 95C SL7SVX... ESTIMATES LARGEST SING. VALUE OF LOWER TRIANG. MATRIX. 96C SL7SVN... ESTIMATES SMALLEST SING. VALUE OF LOWER TRIANG. MATRIX. 97C SRN2GB... UNDERLYING NONLINEAR LEAST-SQUARES SOLVER. 98C SQ7APL... APPLIES HOUSEHOLDER TRANSFORMS STORED BY SQ7RFH. 99C SQ7RFH.... COMPUTES QR FACT. VIA HOUSEHOLDER TRANSFORMS WITH PIVOTING. 100C SR7MDC... RETURNS MACHINE-DEP. CONSTANTS. 101C SS7CPR... PRINTS LINEAR PARAMETERS AT SOLUTION. 102C SV2AXY.... ADDS MULTIPLE OF ONE VECTOR TO ANOTHER. 103C SV7CPY.... COPIES ONE VECTOR TO ANOTHER. 104C SV7PRM.... PERMUTES VECTOR. 105C DV7SCL... SCALES AND COPIES ONE VECTOR TO ANOTHER. 106C 107C *** LOCAL VARIABLES *** 108C 109 INTEGER AR1, CSAVE1, D1, DR1, DR1L, I, I1, 110 1 IPIV1, IER, IV1, J1, JLEN, K, LL1O2, MD, N1, 111 2 NML, NRAN, R1, R1L, RD1 112 REAL SINGTL, T 113 REAL MACHEP, NEGONE, SNGFAC, ZERO 114C 115C *** SUBSCRIPTS FOR IV AND V *** 116C 117 INTEGER AR, CNVCOD, CSAVE, D, FDH, H, IERS, IPIVS, IV1SAV, 118 2 IVNEED, J, LMAT, MODE, NEXTIV, NEXTV, 119 2 NFCALL, NFGCAL, NGCALL, PERM, R, RCOND, 120 3 RDREQ, RDRQSV, REGD, REGD0, RESTOR, TOOBIG, VNEED 121C 122C *** IV SUBSCRIPT VALUES *** 123C 124 PARAMETER (AR=110, CNVCOD=55, CSAVE=105, D=27, FDH=74, H=56, 125 1 IERS=108, IPIVS=109, IV1SAV=104, IVNEED=3, J=70, 126 2 LMAT=42, MODE=35, NEXTIV=46, NEXTV=47, NFCALL=6, 127 3 NFGCAL=7, NGCALL=30, PERM=58, R=61, RCOND=53, RDREQ=57, 128 4 RDRQSV=107, REGD=67, REGD0=82, RESTOR=9, TOOBIG=2, 129 5 VNEED=4) 130 DATA MACHEP/-1.E+0/, NEGONE/-1.E+0/, SNGFAC/1.E+2/, ZERO/0.E+0/ 131C 132C++++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++ 133C 134C 135 IF (IV(1) .EQ. 0) CALL SIVSET(1, IV, LIV, LV, V) 136 N1 = 1 137 NML = N 138 IV1 = IV(1) 139 IF (IV1 .LE. 2) GO TO 20 140C 141C *** CHECK INPUT INTEGERS *** 142C 143 IF (P .LE. 0) GO TO 240 144 IF (L .LT. 0) GO TO 240 145 IF (N .LE. L) GO TO 240 146 IF (LA .LT. N) GO TO 240 147 IF (IV1 .LT. 12) GO TO 20 148 IF (IV1 .EQ. 14) GO TO 20 149 IF (IV1 .EQ. 12) IV(1) = 13 150C 151C *** FRESH START -- COMPUTE STORAGE REQUIREMENTS *** 152C 153 IF (IV(1) .GT. 16) GO TO 240 154 LL1O2 = L*(L+1)/2 155 JLEN = N*P 156 I = L + P 157 IF (IV(1) .NE. 13) GO TO 10 158 IV(IVNEED) = IV(IVNEED) + L 159 IV(VNEED) = IV(VNEED) + P + 2*N + JLEN + LL1O2 + L 160 10 IF (IV(PERM) .LE. AR) IV(PERM) = AR + 1 161 CALL SRN2GB(B, V, V, IV, LIV, LV, N, N, N1, NML, P, V, V, V, ALF) 162 IF (IV(1) .NE. 14) GO TO 999 163C 164C *** STORAGE ALLOCATION *** 165C 166 IV(IPIVS) = IV(NEXTIV) 167 IV(NEXTIV) = IV(NEXTIV) + L 168 IV(D) = IV(NEXTV) 169 IV(REGD0) = IV(D) + P 170 IV(AR) = IV(REGD0) + N 171 IV(CSAVE) = IV(AR) + LL1O2 172 IV(J) = IV(CSAVE) + L 173 IV(R) = IV(J) + JLEN 174 IV(NEXTV) = IV(R) + N 175 IV(IERS) = 0 176 IF (IV1 .EQ. 13) GO TO 999 177C 178C *** SET POINTERS INTO IV AND V *** 179C 180 20 AR1 = IV(AR) 181 D1 = IV(D) 182 DR1 = IV(J) 183 DR1L = DR1 + L 184 R1 = IV(R) 185 R1L = R1 + L 186 RD1 = IV(REGD0) 187 CSAVE1 = IV(CSAVE) 188 NML = N - L 189 IF (IV1 .LE. 2) GO TO 50 190C 191 30 CALL SRN2GB(B, V(D1), V(DR1L), IV, LIV, LV, NML, N, N1, NML, P, 192 1 V(R1L), V(RD1), V, ALF) 193c IF (abs(IV(RESTOR)-2) .EQ. 1 .AND. L .GT. 0) 194c 1 CALL SV7CPY(L, C, V(CSAVE1)) 195 IV1 = IV(1) 196 IF (IV1-2) 40, 150, 230 197C 198C *** NEW FUNCTION VALUE (RESIDUAL) NEEDED *** 199C 200 40 IV(IV1SAV) = IV(1) 201 IV(1) = abs(IV1) 202c IF (IV(RESTOR) .EQ. 2 .AND. L .GT. 0) CALL SV7CPY(L, V(CSAVE1), C 203 GO TO 999 204C 205C *** COMPUTE NEW RESIDUAL OR GRADIENT *** 206C 207 50 IV(1) = IV(IV1SAV) 208 MD = IV(MODE) 209 IF (MD .LE. 0) GO TO 60 210 NML = N 211 DR1L = DR1 212 R1L = R1 213 60 IF (IV(TOOBIG) .NE. 0) GO TO 30 214 IF (abs(IV1) .EQ. 2) GO TO 170 215C 216C *** COMPUTE NEW RESIDUAL *** 217C 218 IF (L1 .LE. L) CALL SV7CPY(N, V(R1), Y) 219 IF (L1 .GT. L) CALL SV2AXY(N, V(R1), NEGONE, A(1,L1), Y) 220 IF (MD .GT. 0) GO TO 120 221 IER = 0 222 IF (L .LE. 0) GO TO 110 223 LL1O2 = L * (L + 1) / 2 224 IPIV1 = IV(IPIVS) 225 CALL SQ7RFH(IER, IV(IPIV1), N, LA, 0, L, A, V(AR1), LL1O2, C) 226C 227C *** DETERMINE NUMERICAL RANK OF A *** 228C 229 IF (MACHEP .LE. ZERO) MACHEP = SR7MDC(3) 230 SINGTL = SNGFAC * real(max(L,N)) * MACHEP 231 K = L 232 IF (IER .NE. 0) K = IER - 1 233 70 IF (K .LE. 0) GO TO 90 234 T = SL7SVX(K, V(AR1), C, C) 235 IF (T .GT. ZERO) T = SL7SVN(K, V(AR1), C, C) / T 236 IF (T .GT. SINGTL) GO TO 80 237 K = K - 1 238 GO TO 70 239C 240C *** RECORD RANK IN IV(IERS)... IV(IERS) = 0 MEANS FULL RANK, 241C *** IV(IERS) .GT. 0 MEANS RANK IV(IERS) - 1. 242C 243 80 IF (K .GE. L) GO TO 100 244 90 IER = K + 1 245 CALL SV7SCP(L-K, C(K+1), ZERO) 246 100 IV(IERS) = IER 247 IF (K .LE. 0) GO TO 110 248C 249C *** APPLY HOUSEHOLDER TRANSFORMATONS TO RESIDUALS... 250C 251 CALL SQ7APL(LA, N, K, A, V(R1), IER) 252C 253C *** COMPUTING C NOW MAY SAVE A FUNCTION EVALUATION AT 254C *** THE LAST ITERATION. 255C 256 CALL SL7ITV(K, C, V(AR1), V(R1)) 257 CALL SV7PRM(L, IV(IPIV1), C) 258C 259 110 IF(IV(1) .LT. 2) GO TO 220 260 GO TO 999 261C 262C 263C *** RESIDUAL COMPUTATION FOR F.D. HESSIAN *** 264C 265 120 IF (L .LE. 0) GO TO 140 266 DO 130 I = 1, L 267 130 CALL SV2AXY(N, V(R1), -C(I), A(1,I), V(R1)) 268 140 IF (IV(1) .GT. 0) GO TO 30 269 IV(1) = 2 270 GO TO 160 271C 272C *** NEW GRADIENT (JACOBIAN) NEEDED *** 273C 274 150 IV(IV1SAV) = IV1 275 IF (IV(NFGCAL) .NE. IV(NFCALL)) IV(1) = 1 276 160 CALL SV7SCP(N*P, V(DR1), ZERO) 277 GO TO 999 278C 279C *** COMPUTE NEW JACOBIAN *** 280C 281 170 continue 282 IF (NDA .LE. 0) GO TO 240 283 DO 180 I = 1, NDA 284 I1 = IN(1,I) - 1 285 IF (I1 .LT. 0) GO TO 180 286 J1 = IN(2,I) 287 K = DR1 + I1*N 288 T = NEGONE 289 IF (J1 .LE. L) T = -C(J1) 290 CALL SV2AXY(N, V(K), T, DA(1,I), V(K)) 291 180 CONTINUE 292 IF (IV1 .EQ. 2) GO TO 190 293 IV(1) = IV1 294 GO TO 999 295 190 IF (L .LE. 0) GO TO 30 296 IF (MD .GT. 0) GO TO 30 297 K = DR1 298 IER = IV(IERS) 299 NRAN = L 300 IF (IER .GT. 0) NRAN = IER - 1 301 IF (NRAN .LE. 0) GO TO 210 302 DO 200 I = 1, P 303 CALL SQ7APL(LA, N, NRAN, A, V(K), IER) 304 K = K + N 305 200 CONTINUE 306 210 CALL SV7CPY(L, V(CSAVE1), C) 307 220 IF (IER .EQ. 0) GO TO 30 308C 309C *** ADJUST SUBSCRIPTS DESCRIBING R AND DR... 310C 311 NRAN = IER - 1 312 DR1L = DR1 + NRAN 313 NML = N - NRAN 314 R1L = R1 + NRAN 315 GO TO 30 316C 317C *** CONVERGENCE OR LIMIT REACHED *** 318C 319 230 IF (IV(REGD) .EQ. 1) IV(REGD) = RD1 320 IF (IV(1) .LE. 11) CALL SS7CPR(C, IV, L, LIV) 321 GO TO 999 322C 323 240 IV(1) = 66 324 CALL SITSUM(V, V, IV, LIV, LV, P, V, ALF) 325C 326 999 RETURN 327C 328C *** LAST CARD OF SRNSGB FOLLOWS *** 329 END 330