/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:42 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dnlsfb.h" #include /* PARAMETER translations */ #define AMAT 113 #define COVREQ 15 #define D 27 #define DAMAT 114 #define DLTFDJ 43 #define GPTR 117 #define GRP 118 #define IN 112 #define IVNEED 3 #define L1SAV 111 #define MAXGRP 116 #define MODE 35 #define MSAVE 115 #define NEXTIV 46 #define NEXTV 47 #define NFCALL 6 #define NFGCAL 7 #define PERM 58 #define RESTOR 9 #define TOOBIG 2 #define VNEED 4 #define XSAVE 119 /* end of PARAMETER translations */ void /*FUNCTION*/ dnlsfb( long n, long p, long l, double alf[], double b[][2], double c[], double y[], void (*dcalca)(long,long,long,double[],long*,double[]), long *inc, long iinc, long iv[], long liv, long lv, double v[]) { #define INC(I_,J_) (*(inc+(I_)*(iinc)+(J_))) LOGICAL32 partj; long int a0, a1, aj, alp1, bwa1, d0, da0, da1, daj, gptr1, grp1, grp2, i, i1, in0, in1, in2, ini, inlen, ipntr1, iv1, iwa1, iwalen, j1, jn1, jpntr1, k, l1, lp1, m, m0, nf, ng, ngrp0, ngrp1, ngrp2, rsave0, rsave1, rsvlen, x0i, xsave0, xsave1; double delta, di, h, xi, xi1; static double negone = -1.e0; static double one = 1.e0; static double zero = 0.e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Alf = &alf[0] - 1; double *const C = &c[0] - 1; long *const Iv = &iv[0] - 1; double *const V = &v[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 1996-04-27 DNLSFB Krogh Changes to get desired C prototypes. *>> 1994-10-20 DNLSFB Krogh Changes to use M77CON *>> 1990-07-02 DNLSFB CLL @ JPL *>> 1990-06-12 CLL @ JPL *>> 1990-02-14 CLL @ JPL *** from netlib, Fri Feb 9 13:10:09 EST 1990 *** *--D replaces "?": ?NLSFB,?CALCA,?IVSET,?RNSGB,?V2AXY,?V7CPY,?V7SCL * * *** SOLVE SEPARABLE NONLINEAR LEAST SQUARES USING * *** FINITE-DIFFERENCE DERIVATIVES. * * *** PARAMETER DECLARATIONS *** * */ /* *** PARAMETERS *** * * N (IN) NUMBER OF OBSERVATIONS. * P (IN) NUMBER OF NONLINEAR PARAMETERS TO BE ESTIMATED. * L (IN) NUMBER OF LINEAR PARAMETERS TO BE ESTIMATED. * ALF (I/O) NONLINEAR PARAMETERS. * INPUT = INITIAL GUESS, * OUTPUT = BEST ESTIMATE FOUND. * B (IN) SIMBLE BOUNDS ON ALF.. B(1,I) .LE. ALF(I) .LE. B(2,I). * C (OUT) LINEAR PARAMETERS (ESTIMATED). * Y (IN) RIGHT-HAND SIDE VECTOR. * DCALCA (IN) SUBROUTINE TO COMPUTE A MATRIX. * INC (IN) INCIDENCE MATRIX OF DEPENDENCIES OF COLUMNS OF A ON * COMPONENTS OF ALF -- INC(I,J) = 1 MEANS COLUMN I * OF A DEPENDS ON ALF(J). * IINC (IN) DECLARED LEAD DIMENSION OF INC. MUST BE AT LEAST L+1. * IV (I/O) INTEGER PARAMETER AND SCRATCH VECTOR. * LIV (IN) LENGTH OF IV. MUST BE AT LEAST * 122 + 2*M + 7*P + 2*L + MAX(L+1,6*P), WHERE M IS * THE NUMBER OF ONES IN INC. * LV (IN) LENGTH OF V. MUST BE AT LEAST * 105 + N*(2*L+6+P) + L*(L+3)/2 + P*(2*P + 22). * If row L+1 of INC() contains only zeros, meaning * the term PHI sub (L+1) is absent from the model, * then LV can be 4*N less than just described. * V (I/O) FLOATING-POINT PARAMETER AND SCRATCH VECTOR. * * ------------------------- DECLARATIONS ---------------------------- * * * *** EXTERNAL SUBROUTINES *** * */ /* DIVSET.... PROVIDES DEFAULT IV AND V VALUES. * IDSM...... DETERMINES EFFICIENT ORDER FOR FINITE DIFFERENCES. * DRNSGB... CARRIES OUT NL2SOL ALGORITHM. * DV2AXY.... ADDS A MULTIPLE OF ONE VECTOR TO ANOTHER. * DV7CPY.... COPIES ONE VECTOR TO ANOTHER. * DV7SCL... SCALES AND COPIES ONE VECTOR TO ANOTHER. * * *** LOCAL VARIABLES *** * */ /* *** SUBSCRIPTS FOR IV AND V *** * */ /* *** IV SUBSCRIPT VALUES *** * * /6 * DATA 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/ * /7 */ /*/ */ /*++++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++ * */ lp1 = l + 1; if (Iv[1] == 0) divset( 1, iv, liv, lv, v ); if ((p <= 0 || l < 0) || iinc <= l) goto L_80; iv1 = Iv[1]; if (iv1 == 14) goto L_120; if (iv1 > 2 && iv1 < 12) goto L_120; if (iv1 == 12) Iv[1] = 13; if (Iv[1] != 13) goto L_50; /* *** FRESH START *** * */ if (Iv[PERM] <= XSAVE) Iv[PERM] = XSAVE + 1; /* *** CHECK INC, COUNT ITS NONZEROS * */ l1 = 0; m = 0; for (i = 1; i <= p; i++) { if (b[i - 1][0] >= b[i - 1][1]) goto L_40; m0 = m; if (l == 0) goto L_20; for (k = 1; k <= l; k++) { if (INC(i - 1,k - 1) < 0 || INC(i - 1,k - 1) > 1) goto L_80; if (INC(i - 1,k - 1) == 1) m += 1; } L_20: if (INC(i - 1,lp1 - 1) != 1) goto L_30; m += 1; l1 = 1; goto L_40; L_30: if ((m == m0 || INC(i - 1,lp1 - 1) < 0) || INC(i - 1,lp1 - 1) > 1) goto L_80; L_40: ; } /* *** NOW L1 = 1 MEANS A HAS COLUMN L+1 *** * * *** COMPUTE STORAGE REQUIREMENTS *** * */ iwalen = max( lp1, 6*p ); inlen = 2*m; Iv[IVNEED] += inlen + 3*p + l + iwalen + 3; rsvlen = 2*l1*n; l1 += l; Iv[VNEED] += 2*n*l1 + rsvlen + p; L_50: drnsgb( v, alf, b, c, v, (long(*)[2])iv, iv, l, 1, n, liv, lv, n, m, p, v, y ); if (Iv[1] != 14) goto L_999; /* *** STORAGE ALLOCATION *** * */ 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; /* *** DETERMINE HOW MANY GROUPS FOR FINITE DIFFERENCES * *** (SET UP TO CALL IDSM) * */ in1 = Iv[IN]; jn1 = in1 + m; for (k = 1; k <= p; k++) { if (b[k - 1][0] >= b[k - 1][1]) goto L_70; for (i = 1; i <= lp1; i++) { if (INC(k - 1,i - 1) == 0) goto L_60; Iv[in1] = i; in1 += 1; Iv[jn1] = k; jn1 += 1; L_60: ; } L_70: ; } in1 = Iv[IN]; jn1 = in1 + m; iwa1 = in1 + inlen; ngrp1 = iwa1 + iwalen; bwa1 = ngrp1 + p; ipntr1 = bwa1 + p; jpntr1 = ipntr1 + l + 2; idsm( lp1, p, m, &Iv[in1], &Iv[jn1], &Iv[ngrp1], &ng, &k, &i, &Iv[ipntr1], &Iv[jpntr1], &Iv[iwa1], iwalen, &Iv[bwa1] ); if (i == 1) goto L_90; Iv[1] = 69; goto L_50; L_80: Iv[1] = 66; goto L_50; /* *** SET UP GRP AND GPTR ARRAYS FOR COMPUTING FINITE DIFFERENCES * * *** THERE ARE NG GROUPS. GROUP I CONTAINS ALF(GRP(J)) FOR * *** GPTR(I) .LE. J .LE. GPTR(I+1)-1. * */ L_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; for (i = 1; i <= ng; i++) { Iv[gptr1] = grp1; gptr1 += 1; for (i1 = ngrp1; i1 <= ngrp2; i1++) { if (Iv[i1] != i) goto L_100; k = i1 - ngrp0; if (b[k - 1][0] >= b[k - 1][1]) goto L_100; Iv[grp1] = k; grp1 += 1; L_100: ; } } Iv[gptr1] = grp1; if (iv1 == 13) goto L_999; /* *** INITIALIZE POINTERS *** * */ L_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] = -labs( Iv[COVREQ] ); L_130: drnsgb( &V[a1], alf, b, c, &V[da1], (long(*)[2])&Iv[in1], iv, l, l1, n, liv, lv, n, l1, p, v, y ); switch (IARITHIF(Iv[1] - 2)) { case -1: goto L_140; case 0: goto L_150; case 1: goto L_999; } /* *** NEW FUNCTION VALUE (R VALUE) NEEDED *** * */ L_140: nf = Iv[NFCALL]; (*dcalca)( n, p, l, alf, &nf, &V[a1] ); if (nf <= 0) Iv[TOOBIG] = 1; /* CALL DCALCA(N, P, L, ALF, NF, V(A1)) */ if (l1 <= l) goto L_130; if (Iv[RESTOR] == 2) dv7cpy( n, &V[rsave0], &V[rsave1] ); dv7cpy( n, &V[rsave1], &V[alp1] ); goto L_130; /* *** COMPUTE DR = GRADIENT OF R COMPONENTS *** * */ L_150: if (l1 > l && Iv[NFGCAL] == Iv[NFCALL]) dv7cpy( n, &V[rsave0], &V[rsave1] ); gptr1 = Iv[GPTR]; for (k = 1; k <= ng; k++) { dv7cpy( p, &V[xsave1], alf ); grp1 = Iv[gptr1]; grp2 = Iv[gptr1 + 1] - 1; gptr1 += 1; for (i1 = grp1; i1 <= grp2; i1++) { i = Iv[i1]; xi = Alf[i]; j1 = d0 + i; di = V[j1]; if (di <= zero) di = one; h = delta*fmax( fabs( xi ), one/di ); if (xi < zero) goto L_160; xi1 = xi + h; if (xi1 <= b[i - 1][1]) goto L_170; xi1 = xi - h; if (xi1 >= b[i - 1][0]) goto L_170; goto L_190; L_160: xi1 = xi - h; if (xi1 >= b[i - 1][0]) goto L_170; xi1 = xi + h; if (xi1 <= b[i - 1][1]) goto L_170; goto L_190; L_170: x0i = xsave0 + i; V[x0i] = xi1; } (*dcalca)( n, p, l, &V[xsave1], &nf, &V[da1] ); if (Iv[NFGCAL] > 0) goto L_200; /* CALL DCALCA(N, P, L, V(XSAVE1), NF, V(DA1)) */ L_190: Iv[TOOBIG] = 1; goto L_130; L_200: jn1 = in1; for (i = in1; i <= in2; i++) { Iv[i] = 0; } partj = Iv[MODE] <= p; for (i1 = grp1; i1 <= grp2; i1++) { i = Iv[i1]; for (j1 = 1; j1 <= l1; j1++) { if (INC(i - 1,j1 - 1) == 0) goto L_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) goto L_220; /* *** FULL FINITE DIFFERENCE FOR COV. AND REG. DIAG. *** */ aj = a0 + j1*n; dv2axy( n, &V[daj], negone, &V[aj], &V[daj] ); goto L_230; L_220: if (j1 > l) dv2axy( n, &V[daj], negone, &V[rsave0], &V[daj] ); L_230: dv7scl( n, &V[daj], h, &V[daj] ); L_240: ; } } if (k >= ng) goto L_270; Iv[1] = -2; drnsgb( &V[a1], alf, b, c, &V[da1], (long(*)[2])&Iv[in1], iv, l, l1, n, liv, lv, n, l1, p, v, y ); if (-2 != Iv[1]) goto L_999; } L_270: Iv[1] = 2; goto L_130; L_999: return; /* *** LAST CARD OF DNLSFB FOLLOWS *** */ #undef INC } /* end of function */