/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:47 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "snlsgu.h" #include /* PARAMETER translations */ #define AMAT 113 #define DAMAT 114 #define IN 112 #define IVNEED 3 #define L1SAV 111 #define MSAVE 115 #define NEXTIV 46 #define NEXTV 47 #define NFCALL 6 #define NFGCAL 7 #define PERM 58 #define TOOBIG 2 #define VNEED 4 /* end of PARAMETER translations */ void /*FUNCTION*/ snlsgu( long n, long p, long l, float alf[], float c[], float y[], void (*scalca)(long,long,long,float[],long*,float[]), void (*scalcb)(long,long,long,float[],long*,float[]), long *inc, long iinc, long iv[], long liv, long lv, float v[]) { #define INC(I_,J_) (*(inc+(I_)*(iinc)+(J_))) long int a1, da1, i, in1, iv1, k, l1, lp1, m, m0, nf; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Alf = &alf[0] - 1; float *const C = &c[0] - 1; long *const Iv = &iv[0] - 1; float *const V = &v[0] - 1; float *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 SNLSGU Krogh Changes to get desired C prototypes. *>> 1994-10-20 SNLSGU Krogh Changes to use M77CON *>> 1990-07-02 SNLSGU CLL @ JPL *>> 1990-06-12 CLL @ JPL *>> 1990-02-20 CLL @ JPL *--S replaces "?": ?NLSGU,?CALCA,?CALCB,?IVSET,?RNSG * * *** SOLVE SEPARABLE NONLINEAR LEAST SQUARES USING *** * *** ANALYTICALLY COMPUTED DERIVATIVES. *** * * *** PARAMETER DECLARATIONS *** * */ /* *** PURPOSE *** * * GIVEN A SET OF N OBSERVATIONS Y(1)....Y(N) OF A DEPENDENT VARIABLE * T(1)...T(N), SNLSGU ATTEMPTS TO COMPUTE A LEAST SQUARES FIT * TO A FUNCTION ETA (THE MODEL) WHICH IS A LINEAR COMBINATION * * L * ETA(C,ALF,T) = SUM C * PHI(ALF,T) +PHI (ALF,T) * J=1 J J L+1 * * OF NONLINEAR FUNCTIONS PHI(J) DEPENDENT ON T AND ALF(1),...,ALF(P) * (.E.G. A SUM OF EXPONENTIALS OR GAUSSIANS). THAT IS, IT DETERMINES * NONLINEAR PARAMETERS ALF WHICH MINIMIZE * * 2 N 2 * NORM(RESIDUAL) = SUM (Y - ETA(C,ALF,T )). * I=1 I I * * THE (L+1)ST TERM IS OPTIONAL. * * ------------------------- PARAMETER USAGE ------------------------- * * INPUT PARAMETERS * * N INTEGER NUMBER OF OBSERVATIONS (MUST BE .GE. MAX(L,P)). * * P INTEGER NUMBER OF NONLINEAR PARAMETERS (MUST BE .GE. 1). * * L INTEGER NUMBER OF LINEAR PARAMETERS (MUST BE .GE. 0). * * ALF D.P. ARRAY P VECTOR = INITIAL ESTIMATE OF THE NONLINEAR * PARAMETERS. * * SCALCA SUBROUTINE USER PROVIDED FUNCTION TO CALCULATE THE MODEL * (I.E., TO CALCULATE PHI) -- SEE THE NOTE BELOW * ON THE CALLING SEQUENCE FOR SCALCA. * SCALCA MUST BE DECLARED EXTERNAL IN THE CALLING * PROGRAM. * * SCALCB SUBROUTINE USER PROVIDED FUNCTION TO CALCULATE THE DERIVA- * TIVE OF THE MODEL (I.E., OF PHI) WITH RESPECT TO * ALF -- SEE THE NOTE BELOW ON THE CALLING * SEQUENCE FOR SCALCB. SCALCB MUST BE DECLARED * EXTERNAL IN THE CALLING PROGRAM. * * Y D.P. ARRAY VECTOR OF OBSERVATIONS. * * INC INTEGER ARRAY A 2 DIM. ARRAY OF DIMENSION AT LEAST (L+1,P) * INDICATING THE POSITION OF THE NONLINEAR PARA- * METERS IN THE MODEL. SET INC(J,K) = 1 IF ALF(K) * APPEARS IN PHI(J). OTHERWISE SET INC(J,K) = 0. * IF PHI((L+1)) IS NOT IN THE MODEL, SET THE L+1ST * ROW OF INC TO ALL ZEROS. EVERY COLUMN OF INC * MUST CONTAIN AT LEAST ONE 1. * * IINC INTEGER DECLARED ROW DIMENSION OF INC, WHICH MUST BE AT * LEAST L+1. * * IV INTEGER ARRAY OF LENGTH AT LEAST LIV THAT CONTAINS * VARIOUS PARAMETERS FOR THE SUBROUTINE, SUCH AS * THE ITERATION AND FUNCTION EVALUATION LIMITS AND * SWITCHES THAT CONTROL PRINTING. THE INPUT COM- * PONENTS OF IV ARE DESCRIBED IN DETAIL IN THE * PORT OPTIMIZATION DOCUMENTATION. * IF IV(1)=0 ON INPUT, THEN DEFAULT PARAMETERS * ARE SUPPLIED TO IV AND V. THE CALLER MAY SUPPLY * NONDEFAULT PARAMETERS TO IV AND V BY EXECUTING A * CALL SIVSET(1,IV,LIV,LV,V) AND THEN ASSIGNING * NONDEFAULT VALUES TO THE APPROPRIATE COMPONENTS * OF IV AND V BEFORE CALLING SNLSGU. * * LIV INTEGER LENGTH OF IV. MUST BE AT LEAST 115+P+L + 2*M, * WHERE M IS THE NUMBER OF ONES IN INC. * * LV INTEGER LENGTH OF V. MUST BE AT LEAST * 105 + N*(L+M+3) + JLEN + L*(L+3)/2 + P*(2*P+17), * WHERE M IS AS FOR LIV (SEE ABOVE) AND * JLEN = (L+P)*(N+L+P+1), UNLESS NEITHER A * COVARIANCE MATRIX NOR REGRESSION DIAGNOSTICS ARE * REQUESTED, IN WHICH CASE JLEN = N*P. * 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 N less than just described. * * V D.P. ARRAY WORK AND PARAMETER ARRAY OF LENGTH AT LEAST LV * THAT CONTAINS SUCH INPUT COMPONENTS AS THE * CONVERGENCE TOLERANCES. THE INPUT COMPONENTS OF * V MAY BE SUPPLIED AS FOR IV (SEE ABOVE). NOTE * THAT V(35) CONTAINS THE INITIAL STEP BOUND, * WHICH, IF TOO LARGE, MAY LEAD TO OVERFLOW. * * OUTPUT PARAMETERS * * ALF D.P. ARRAY FINAL NONLINEAR PARAMETERS. * * C D.P. ARRAY L VECTOR OF LINEAR PARAMETERS -- NOTE THAT NO * INITIAL GUESS FOR C IS REQUIRED. * * IV IV(1) CONTAINS A RETURN CODE DESCRIBED IN THE */ /* PORT OPTIMIZATION DOCUMENTATION. IF IV(1) LIES * BETWEEN 3 AND 7, THEN THE ALGORITHM HAS * CONVERGED (BUT IV(1) = 7 INDICATES POSSIBLE * TROUBLE WITH THE MODEL). IV(1) = 9 OR 10 MEANS * FUNCTION EVALUATION OR ITERATION LIMIT REACHED. * IV(1) = 66 MEANS BAD PARAMETERS (INCLUDING A * COLUMN OF ZEROS IN INC). NOTE THAT THE * ALGORITHM CAN BE RESTARTED AFTER ANY RETURN WITH * IV(1) .LT. 12 -- SEE THE PORT DOCUMENTATION. * * V VARIOUS ITEMS OF INTEREST, INCLUDING THE NORM OF * THE GRADIENT(1) AND THE FUNCTION VALUE(10). SEE * THE PORT DOCUMENTATION FOR A COMPLETE LIST. * * * * PARAMETERS FOR SCALCA(N,P,L,ALF,NF,PHI) * * N,L,P,ALF ARE INPUT PARAMETERS AS DESCRIBED ABOVE * * PHI D.P. ARRAY N*(L+1) or N*L array. SCALCA must store into * PHI(i,j) the value of of the nonlinear function * PHI sub j evaluated with the current values of * ALF() and at the ith point of the data set. * Note that if the model does not have the * term PHI sub L+1 then one must not store to the * (L+1)st column of PHI() as this would overwrite * other work space. * * NF INTEGER CURRENT INVOCATION COUNT FOR SCALCA. IF PHI CANNOT * BE EVALUATED AT ALF (E.G. BECAUSE AN ARGUMENT TO * AN INTRINSIC FUNCTION IS OUT OF RANGE), THEN SCALCA * SHOULD SIMPLY SET NF TO 0 AND RETURN. THIS * TELLS THE ALGORITHM TO TRY A SMALLER STEP. * * N.B. THE DEPENDENT VARIABLE T IS NOT EXPLICITLY PASSED. IF REQUIRED, * IT MAY BE PASSED IN NAMED COMMON. * * * PARAMETERS FOR SCALCB(N,P,L,ALF,NF,DER) * * N,P,L,ALF,NF ARE AS FOR SCALCA * * DER D.P. ARRAY N*M ARRAY, WHERE M IS THE NUMBER OF ONES IN INC. * SCALCB MUST SET DER TO THE DERIVATIVES OF THE MODEL * WITH RESPECT TO ALF. IF THE MODEL HAS K TERMS THAT * DEPEND ON ALF(I), THEN DER WILL HAVE K CONSECUTIVE * COLUMNS OF DERIVATIVES WITH RESPECT TO ALF(I). THE * COLUMNS OF DER CORRESPOND TO THE ONES IN INC WHEN * ONE TRAVELS THROUGH INC BY COLUMNS. FOR EXAMPLE, * IF INC HAS THE FORM... * 1 1 0 * 0 1 0 * 1 0 0 * 0 0 1 * THEN THE FIRST TWO COLUMNS OF DER ARE FOR THE * DERIVATIVES OF COLUMNS 1 AND 3 OF PHI WITH RESPECT * TO ALF(1), COLUMNS 3 AND 4 OF DER ARE FOR THE * DERIVATIVES OF COLUMNS 1 AND 2 OF PHI WITH RESPECT * TO ALF(2), AND COLUMN 5 OF DER IS FOR THE DERIVA- * TIVE OF COLUMN 4 OF PHI WITH RESPECT TO ALF(3). * MORE SPECIFICALLY, DER(I,2) IS FOR THE DERIVATIVE * OF PHI(I,3) WITH RESPECT TO ALF(1) AND DER(I,5) IS * FOR THE DERIVATIVE OF PHI(I,4) WITH RESPECT TO * ALF(3) (FOR I = 1,2,...,N). * THE VALUE OF ALF PASSED TO SCALCB IS THE SAME AS * THAT PASSED TO SCALCA THE LAST TIME IT WAS CALLED. * (IF DER CANNOT BE EVALUATED, THEN SCALCB SHOULD SET * NF TO 0. THIS WILL CAUSE AN ERROR RETURN.) * * N.B. DER IS FOR DERIVATIVES WITH RESPECT TO ALF, NOT T. * * ----------------------------- NOTES ------------------------------- * * THIS PROGRAM WAS WRITTEN BY LINDA KAUFMAN AT BELL LABS, MURRAY * HILL, N.J. IN 1977 AND EXTENSIVELY REVISED BY HER AND DAVID GAY IN * 1980, 1981, 1983, 1984. THE WORK OF DAVID GAY WAS SUPPORTED IN PART * BY NATIONAL SCIENCE FOUNDATION GRANT MCS-7906671. * * ------------------------- DECLARATIONS ---------------------------- * * * *** EXTERNAL SUBROUTINES *** * */ /* SIVSET.... PROVIDES DEFAULT IV AND V VALUES. * SRNSG... CARRIES OUT NL2SOL ALGORITHM. * * *** LOCAL VARIABLES *** * */ /* *** SUBSCRIPTS FOR IV AND V *** * */ /* *** IV SUBSCRIPT VALUES *** * * /6 * DATA AMAT/113/, D/27/, DAMAT/114/, IN/112/, IVNEED/3/, J/70/, * 1 L1SAV/111/, MSAVE/115/, NEXTIV/46/, NEXTV/47/, NFCALL/6/, * 2 NFGCAL/7/, PERM/58/, R/61/, TOOBIG/2/, VNEED/4/ * /7 */ /*/ * *++++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++ * */ if (Iv[1] == 0) sivset( 1, iv, liv, lv, v ); if ((p <= 0 || l < 0) || iinc <= l) goto L_50; iv1 = Iv[1]; if (iv1 == 14) goto L_90; if (iv1 > 2 && iv1 < 12) goto L_90; if (iv1 == 12) Iv[1] = 13; if (Iv[1] != 13) goto L_60; if (Iv[PERM] <= MSAVE) Iv[PERM] = MSAVE + 1; lp1 = l + 1; l1 = 0; m = 0; for (i = 1; i <= p; i++) { 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_50; 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; L_30: if ((m == m0 || INC(i - 1,lp1 - 1) < 0) || INC(i - 1,lp1 - 1) > 1) goto L_50; } Iv[IVNEED] += 2*m; l1 += l; Iv[VNEED] += n*(l1 + m); goto L_60; L_50: Iv[1] = 66; L_60: srnsg( v, alf, 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[NEXTIV] = Iv[IN] + 2*m; Iv[AMAT] = Iv[NEXTV]; Iv[DAMAT] = Iv[AMAT] + n*l1; Iv[NEXTV] = Iv[DAMAT] + n*m; Iv[L1SAV] = l1; Iv[MSAVE] = m; /* *** SET UP IN ARRAY *** * */ in1 = Iv[IN]; for (i = 1; i <= p; i++) { for (k = 1; k <= lp1; k++) { if (INC(i - 1,k - 1) == 0) goto L_70; Iv[in1] = i; Iv[in1 + 1] = k; in1 += 2; L_70: ; } } if (iv1 == 13) goto L_999; L_90: a1 = Iv[AMAT]; da1 = Iv[DAMAT]; in1 = Iv[IN]; l1 = Iv[L1SAV]; m = Iv[MSAVE]; L_100: srnsg( &V[a1], alf, c, &V[da1], (long(*)[2])&Iv[in1], iv, l, l1, n, liv, lv, n, m, p, v, y ); switch (IARITHIF(Iv[1] - 2)) { case -1: goto L_110; case 0: goto L_120; case 1: goto L_999; } /* *** NEW FUNCTION VALUE (R VALUE) NEEDED *** * */ L_110: nf = Iv[NFCALL]; (*scalca)( n, p, l, alf, &nf, &V[a1] ); if (nf <= 0) Iv[TOOBIG] = 1; /* CALL SCALCA(N, P, L, ALF, NF, V(A1)) */ goto L_100; /* *** COMPUTE DR = GRADIENT OF R COMPONENTS *** * */ L_120: ; (*scalcb)( n, p, l, alf, &Iv[NFGCAL], &V[da1] ); if (Iv[NFGCAL] == 0) Iv[TOOBIG] = 1; /* CALL SCALCB(N, P, L, ALF, IV(NFGCAL), V(DA1)) */ goto L_100; L_999: return; /* *** LAST CARD OF SNLSGU FOLLOWS *** */ #undef INC } /* end of function */