SUBROUTINE DNLSFB(N, P, L, ALF, B, C, Y, DCALCA, INC, IINC, IV, 1 LIV, LV, V) c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. c ALL RIGHTS RESERVED. c Based on Government Sponsored Research NAS7-03001. c>> 1996-04-27 DNLSFB Krogh Changes to get desired C prototypes. c>> 1994-10-20 DNLSFB Krogh Changes to use M77CON c>> 1990-07-02 DNLSFB CLL @ JPL c>> 1990-06-12 CLL @ JPL c>> 1990-02-14 CLL @ JPL *** from netlib, Fri Feb 9 13:10:09 EST 1990 *** c--D replaces "?": ?NLSFB,?CALCA,?IVSET,?RNSGB,?V2AXY,?V7CPY,?V7SCL C C *** SOLVE SEPARABLE NONLINEAR LEAST SQUARES USING C *** FINITE-DIFFERENCE DERIVATIVES. C C *** PARAMETER DECLARATIONS *** C EXTERNAL DCALCA INTEGER IINC, L, LIV, LV, N, P INTEGER INC(IINC,P), IV(LIV) DOUBLE PRECISION ALF(P), C(L), B(2,P), V(LV), Y(N) C C *** PARAMETERS *** C C N (IN) NUMBER OF OBSERVATIONS. C P (IN) NUMBER OF NONLINEAR PARAMETERS TO BE ESTIMATED. C L (IN) NUMBER OF LINEAR PARAMETERS TO BE ESTIMATED. C ALF (I/O) NONLINEAR PARAMETERS. C INPUT = INITIAL GUESS, C OUTPUT = BEST ESTIMATE FOUND. C B (IN) SIMBLE BOUNDS ON ALF.. B(1,I) .LE. ALF(I) .LE. B(2,I). C C (OUT) LINEAR PARAMETERS (ESTIMATED). C Y (IN) RIGHT-HAND SIDE VECTOR. C DCALCA (IN) SUBROUTINE TO COMPUTE A MATRIX. C INC (IN) INCIDENCE MATRIX OF DEPENDENCIES OF COLUMNS OF A ON C COMPONENTS OF ALF -- INC(I,J) = 1 MEANS COLUMN I C OF A DEPENDS ON ALF(J). C IINC (IN) DECLARED LEAD DIMENSION OF INC. MUST BE AT LEAST L+1. C IV (I/O) INTEGER PARAMETER AND SCRATCH VECTOR. C LIV (IN) LENGTH OF IV. MUST BE AT LEAST C 122 + 2*M + 7*P + 2*L + MAX(L+1,6*P), WHERE M IS C THE NUMBER OF ONES IN INC. C LV (IN) LENGTH OF V. MUST BE AT LEAST C 105 + N*(2*L+6+P) + L*(L+3)/2 + P*(2*P + 22). c If row L+1 of INC() contains only zeros, meaning c the term PHI sub (L+1) is absent from the model, c then LV can be 4*N less than just described. C V (I/O) FLOATING-POINT PARAMETER AND SCRATCH VECTOR. C c ------------------------- DECLARATIONS ---------------------------- C C C *** EXTERNAL SUBROUTINES *** C EXTERNAL DIVSET, IDSM, DRNSGB,DV2AXY,DV7CPY, DV7SCL C C DIVSET.... PROVIDES DEFAULT IV AND V VALUES. C IDSM...... DETERMINES EFFICIENT ORDER FOR FINITE DIFFERENCES. C DRNSGB... CARRIES OUT NL2SOL ALGORITHM. C DV2AXY.... ADDS A MULTIPLE OF ONE VECTOR TO ANOTHER. C DV7CPY.... COPIES ONE VECTOR TO ANOTHER. C DV7SCL... SCALES AND COPIES ONE VECTOR TO ANOTHER. C C *** LOCAL VARIABLES *** C LOGICAL PARTJ INTEGER A0, A1, AJ, ALP1, BWA1, D0, DA0, DA1, DAJ, GPTR1, GRP1, 1 GRP2, I, I1, IN0, IN1, IN2, INI, INLEN, IPNTR1, IV1, IWA1, 2 IWALEN, J1, JN1, JPNTR1, K, L1, LP1, M, M0, NF, NG, NGRP0, 3 NGRP1, NGRP2, RSAVE0, RSAVE1, RSVLEN, X0I, XSAVE0, XSAVE1 DOUBLE PRECISION DELTA, DI, H, XI, XI1 DOUBLE PRECISION NEGONE, ONE, ZERO C C *** SUBSCRIPTS FOR IV AND V *** C INTEGER AMAT, COVREQ, D, DAMAT, DLTFDJ, GPTR, GRP, IN, IVNEED, J, 1 L1SAV, MAXGRP, MODE, MSAVE, NEXTIV, NEXTV, NFCALL, NFGCAL, 2 PERM, R, RESTOR, TOOBIG, VNEED, XSAVE C C *** IV SUBSCRIPT VALUES *** C C/6 C DATA AMAT/113/, COVREQ/15/, D/27/, DAMAT/114/, DLTFDJ/43/, C 1 GPTR/117/, GRP/118/, IN/112/, IVNEED/3/, J/70/, L1SAV/111/, C 2 MAXGRP/116/, MODE/35/, MSAVE/115/, NEXTIV/46/, NEXTV/47/, C 3 NFCALL/6/, NFGCAL/7/, PERM/58/, R/61/, RESTOR/9/, TOOBIG/2/, C 4 VNEED/4/, XSAVE/119/ C/7 PARAMETER (AMAT=113, COVREQ=15, D=27, DAMAT=114, DLTFDJ=43, 1 GPTR=117, GRP=118, IN=112, IVNEED=3, J=70, L1SAV=111, 2 MAXGRP=116, MODE=35, MSAVE=115, NEXTIV=46, NEXTV=47, 3 NFCALL=6, NFGCAL=7, PERM=58, R=61, RESTOR=9, TOOBIG=2, 4 VNEED=4, XSAVE=119) C/ DATA NEGONE/-1.D+0/, ONE/1.D+0/, ZERO/0.D+0/ C C++++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++ C LP1 = L + 1 IF (IV(1) .EQ. 0) CALL DIVSET(1, IV, LIV, LV, V) IF (P .LE. 0 .OR. L .LT. 0 .OR. IINC .LE. L) GO TO 80 IV1 = IV(1) IF (IV1 .EQ. 14) GO TO 120 IF (IV1 .GT. 2 .AND. IV1 .LT. 12) GO TO 120 IF (IV1 .EQ. 12) IV(1) = 13 IF (IV(1) .NE. 13) GO TO 50 C C *** FRESH START *** C IF (IV(PERM) .LE. XSAVE) IV(PERM) = XSAVE + 1 C C *** CHECK INC, COUNT ITS NONZEROS C L1 = 0 M = 0 DO 40 I = 1, P IF (B(1,I) .GE. B(2,I)) GO TO 40 M0 = M IF (L .EQ. 0) GO TO 20 DO 10 K = 1, L IF (INC(K,I) .LT. 0 .OR. INC(K,I) .GT. 1) GO TO 80 IF (INC(K,I) .EQ. 1) M = M + 1 10 CONTINUE 20 IF (INC(LP1,I) .NE. 1) GO TO 30 M = M + 1 L1 = 1 GO TO 40 30 IF (M .EQ. M0 .OR. INC(LP1,I) .LT. 0 1 .OR. INC(LP1,I) .GT. 1) GO TO 80 40 CONTINUE C C *** NOW L1 = 1 MEANS A HAS COLUMN L+1 *** C C *** COMPUTE STORAGE REQUIREMENTS *** C IWALEN = max(LP1, 6*P) INLEN = 2 * M IV(IVNEED) = IV(IVNEED) + INLEN + 3*P + L + IWALEN + 3 RSVLEN = 2 * L1 * N L1 = L + L1 IV(VNEED) = IV(VNEED) + 2*N*L1 + RSVLEN + P C 50 CALL DRNSGB(V, ALF, B, C, V, IV, IV, L, 1, N, LIV, LV, N, M, P, V, 1 Y) IF (IV(1) .NE. 14) GO TO 999 C C *** STORAGE ALLOCATION *** C IV(IN) = IV(NEXTIV) IV(AMAT) = IV(NEXTV) IV(DAMAT) = IV(AMAT) + N*L1 IV(XSAVE) = IV(DAMAT) + N*L1 IV(NEXTV) = IV(XSAVE) + P + RSVLEN IV(L1SAV) = L1 IV(MSAVE) = M C C *** DETERMINE HOW MANY GROUPS FOR FINITE DIFFERENCES C *** (SET UP TO CALL IDSM) C IN1 = IV(IN) JN1 = IN1 + M DO 70 K = 1, P IF (B(1,K) .GE. B(2,K)) GO TO 70 DO 60 I = 1, LP1 IF (INC(I,K) .EQ. 0) GO TO 60 IV(IN1) = I IN1 = IN1 + 1 IV(JN1) = K JN1 = JN1 + 1 60 CONTINUE 70 CONTINUE IN1 = IV(IN) JN1 = IN1 + M IWA1 = IN1 + INLEN NGRP1 = IWA1 + IWALEN BWA1 = NGRP1 + P IPNTR1 = BWA1 + P JPNTR1 = IPNTR1 + L + 2 CALL IDSM(LP1, P, M, IV(IN1), IV(JN1), IV(NGRP1), NG, K, I, 1 IV(IPNTR1), IV(JPNTR1), IV(IWA1), IWALEN, IV(BWA1)) IF (I .EQ. 1) GO TO 90 IV(1) = 69 GO TO 50 80 IV(1) = 66 GO TO 50 C C *** SET UP GRP AND GPTR ARRAYS FOR COMPUTING FINITE DIFFERENCES C C *** THERE ARE NG GROUPS. GROUP I CONTAINS ALF(GRP(J)) FOR C *** GPTR(I) .LE. J .LE. GPTR(I+1)-1. C 90 IV(MAXGRP) = NG IV(GPTR) = IN1 + 2*L1 GPTR1 = IV(GPTR) IV(GRP) = GPTR1 + NG + 1 IV(NEXTIV) = IV(GRP) + P GRP1 = IV(GRP) NGRP0 = NGRP1 - 1 NGRP2 = NGRP0 + P DO 110 I = 1, NG IV(GPTR1) = GRP1 GPTR1 = GPTR1 + 1 DO 100 I1 = NGRP1, NGRP2 IF (IV(I1) .NE. I) GO TO 100 K = I1 - NGRP0 IF (B(1,K) .GE. B(2,K)) GO TO 100 IV(GRP1) = K GRP1 = GRP1 + 1 100 CONTINUE 110 CONTINUE IV(GPTR1) = GRP1 IF (IV1 .EQ. 13) GO TO 999 C C *** INITIALIZE POINTERS *** C 120 A1 = IV(AMAT) A0 = A1 - N DA1 = IV(DAMAT) DA0 = DA1 - N IN1 = IV(IN) IN0 = IN1 - 2 L1 = IV(L1SAV) IN2 = IN1 + 2*L1 - 1 D0 = IV(D) - 1 NG = IV(MAXGRP) XSAVE1 = IV(XSAVE) XSAVE0 = XSAVE1 - 1 RSAVE1 = XSAVE1 + P RSAVE0 = RSAVE1 + N ALP1 = A1 + L*N DELTA = V(DLTFDJ) IV(COVREQ) = -abs(IV(COVREQ)) C 130 CALL DRNSGB(V(A1), ALF, B, C, V(DA1), IV(IN1), IV, L, L1, N, LIV, 1 LV, N, L1, P, V, Y) IF (IV(1)-2) 140, 150, 999 C C *** NEW FUNCTION VALUE (R VALUE) NEEDED *** C 140 NF = IV(NFCALL) C%% (*dcalca)( n, p, l, alf, &nf, &V[a1] ); CALL DCALCA(N, P, L, ALF, NF, V(A1)) IF (NF .LE. 0) IV(TOOBIG) = 1 IF (L1 .LE. L) GO TO 130 IF (IV(RESTOR) .EQ. 2) CALL DV7CPY(N, V(RSAVE0), V(RSAVE1)) CALL DV7CPY(N, V(RSAVE1), V(ALP1)) GO TO 130 C C *** COMPUTE DR = GRADIENT OF R COMPONENTS *** C 150 IF (L1 .GT. L .AND. IV(NFGCAL) .EQ. IV(NFCALL)) 1 CALL DV7CPY(N, V(RSAVE0), V(RSAVE1)) GPTR1 = IV(GPTR) DO 260 K = 1, NG CALL DV7CPY(P, V(XSAVE1), ALF) GRP1 = IV(GPTR1) GRP2 = IV(GPTR1+1) - 1 GPTR1 = GPTR1 + 1 DO 180 I1 = GRP1, GRP2 I = IV(I1) XI = ALF(I) J1 = D0 + I DI = V(J1) IF (DI .LE. ZERO) DI = ONE H = DELTA * max(abs(XI), ONE/DI) IF (XI .LT. ZERO) GO TO 160 XI1 = XI + H IF (XI1 .LE. B(2,I)) GO TO 170 XI1 = XI - H IF (XI1 .GE. B(1,I)) GO TO 170 GO TO 190 160 XI1 = XI - H IF (XI1 .GE. B(1,I)) GO TO 170 XI1 = XI + H IF (XI1 .LE. B(2,I)) GO TO 170 GO TO 190 170 X0I = XSAVE0 + I V(X0I) = XI1 180 CONTINUE C%% (*dcalca)( n, p, l, &V[xsave1], &nf, &V[da1] ); CALL DCALCA(N, P, L, V(XSAVE1), NF, V(DA1)) IF (IV(NFGCAL) .GT. 0) GO TO 200 190 IV(TOOBIG) = 1 GO TO 130 200 JN1 = IN1 DO 210 I = IN1, IN2 210 IV(I) = 0 PARTJ = IV(MODE) .LE. P DO 250 I1 = GRP1, GRP2 I = IV(I1) DO 240 J1 = 1, L1 IF (INC(J1,I) .EQ. 0) GO TO 240 INI = IN0 + 2*J1 IV(INI) = I IV(INI+1) = J1 X0I = XSAVE0 + I H = ONE / (V(X0I) - ALF(I)) DAJ = DA0 + J1*N IF (PARTJ) GO TO 220 C *** FULL FINITE DIFFERENCE FOR COV. AND REG. DIAG. *** AJ = A0 + J1*N CALL DV2AXY(N, V(DAJ), NEGONE, V(AJ), V(DAJ)) GO TO 230 220 IF (J1 .GT. L) 1 CALL DV2AXY(N, V(DAJ), NEGONE, V(RSAVE0), V(DAJ)) 230 CALL DV7SCL(N, V(DAJ), H, V(DAJ)) 240 CONTINUE 250 CONTINUE IF (K .GE. NG) GO TO 270 IV(1) = -2 CALL DRNSGB(V(A1), ALF, B, C, V(DA1), IV(IN1), IV, L, L1, N, 1 LIV, LV, N, L1, P, V, Y) IF (-2 .NE. IV(1)) GO TO 999 260 CONTINUE 270 IV(1) = 2 GO TO 130 C 999 RETURN C C *** LAST CARD OF DNLSFB FOLLOWS *** END