/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:18 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_snlsgu s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include #include #include "p_snlsgu.h" /* program DRSNLSGU *>> 2001-05-24 DRSNLSGU Krogh Minor change for making .f90 version. *>> 1997-06-18 DRSNLSGU Krogh Changes to improve C portability. *>> 1994-11-02 DRSNLSGU Krogh Changes to use M77CON *>> 1994-09-14 DRSNLSGU CLL Set IV(OUTLEV) = 0 for comparing output. *>> 1992-04-13 CLL Rename and reorder common block [D/S]KEY. *>> 1992-02-03 CLL @ JPL *>> 1990-07-02 CLL @ JPL *>> 1990-06-27 CLL @ JPL *>> 1990-04-05 CLL @ JPL *>> 1990-03-29 CLL @ JPL * Demo driver for SNLSGU. A variant of the nonlinear LS code NL2SOL. * SNLSGU solves the "separable" problem. * SNLSGU requires values of the function and the Jacobian matrix. * Note: MDER is the number of ones in the array IND(). * ------------------------------------------------------------------ *--S replaces "?": DR?NLSGU, ?NLSGU, ?CALCA, ?CALCB, ?IVSET, ?KEY * ------------------------------------------------------------------ */ /* PARAMETER translations */ #define COVPRT 14 #define F 10 #define LIV (115 + MA + MB + 2*MDER) #define LV (105 + MDATA*(MB + MDER + 3) + MLEN + (MB*(MB + 3))/2 + MA*(2*MA + 17)) #define MA 2 #define MB 5 #define MDATA 30 #define MDER 4 #define MLEN ((MB + MA)*(MDATA + MB + MA + 1)) #define OUTLEV 19 #define SOLPRT 22 #define STATPR 23 #define X0PRT 24 /* end of PARAMETER translations */ int main( ) { long int iterm, iv[LIV], ivar, na, nb, ndata, _i, _r; static long int ind[MA][MB + 1]; float alf[MA], bet[MB], dof, v[LV]; static float ydata[MDATA]={1.700641e0,1.793512e0,1.838309e0,1.838416e0, 1.792204e0,1.700501e0,1.579804e0,1.426268e0,1.260724e0,1.084901e0, 0.917094e0,0.761920e0,0.627304e0,0.522146e0,0.446645e0,0.404920e0, 0.392033e0,0.409622e0,0.453045e0,0.510765e0,0.584554e0,0.663109e0, 0.747613e0,0.829439e0,0.908496e0,0.983178e0,1.051046e0,1.114072e0, 1.171746e0,1.227823e0}; static int _aini = 1; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Alf = &alf[0] - 1; float *const Bet = &bet[0] - 1; long *const Iv = &iv[0] - 1; float *const V = &v[0] - 1; float *const Ydata = &ydata[0] - 1; /* end of OFFSET VECTORS */ if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */ { static long _itmp0[] = {0,0,1,0,1,0,0,1,0,1,0,0}; for (iterm = 1, _r = 0; iterm <= 6; iterm++) { for (ivar = 1; ivar <= 2; ivar++) { ind[ivar - 1][iterm - 1] = _itmp0[_r++]; } } } _aini = 0; } /* ------------------------------------------------------------------ */ ndata = MDATA; na = MA; nb = MB; Alf[1] = 5.0e0; Alf[2] = 10.0e0; Iv[1] = 0; printf(" Program DRSNLSGU.. Demo driver for SNLSGU.\n A variant of NL2SOL.\n " " SNLSGU handles the Separable problem.\n " " SNLSGU requires values of the function and the Jacobian.\n \n " "Sample problem is a nonlinear curve fit to data.\n " "Model function is B1 + B2 * cosf(A1*t) + B3 * sinf(A1*t) +\n " " B4 * cosf(A2*t) + B5 * sinf(A2*t) + Noise\n " "Data generated using\n (A1, A2, B1, ..., B5) = (6, 9, 1, 0.5, 0.4, 0.2, 0.1)\n " "and Gaussian noise with mean 0 and\n sample standard deviation 0.001\n \n"); sivset( 1, iv, LIV, LV, v ); Iv[X0PRT] = 1; Iv[OUTLEV] = 0; Iv[STATPR] = 1; Iv[SOLPRT] = 1; Iv[COVPRT] = 1; snlsgu( ndata, na, nb, alf, bet, ydata, scalca, scalcb, (long*)ind, nb + 1, iv, LIV, LV, v ); dof = max( ndata - na - nb, 1 ); printf(" \n SIGFAC: sqrtf((2 * V(F))/DOF) =%12.4g\n", sqrtf( 2.0e0*V[F]/dof )); exit(0); } /* end of function */ /* ================================================================== */ /* COMMON translations */ struct t_skey { float c1[MDATA], c2[MDATA], s1[MDATA], s2[MDATA]; long int key; } skey; /* end of COMMON translations */ void /*FUNCTION*/ scalca( long ndata, long na, long nb, float alf[], long *ncount, float *phi) { #define PHI(I_,J_) (*(phi+(I_)*(ndata)+(J_))) long int i; float del, t; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Alf = &alf[0] - 1; float *const C1 = &skey.c1[0] - 1; float *const C2 = &skey.c2[0] - 1; float *const S1 = &skey.s1[0] - 1; float *const S2 = &skey.s2[0] - 1; /* end of OFFSET VECTORS */ /* Test case for separable nonlinear least squares computation. * Computes MDATA x NB matrix PHI as a function of the * nonlinear parameters ALF(). * For J .le. NB the (I,J) term of PHI is the coefficient of the * linear coefficient B(J) in row I of the model. * In this example the model does not have a term that is not * multiplied by a linear coefficient. If such a term is present * then PHI must have an (NB+1)st column to hold this term. * This code Illustrates saving results in common between SCALCA and * SCALCB to avoid recalculation of common subexpressions. * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ t = 0.0e0; del = 1.0e0/29.0e0; skey.key = *ncount; for (i = 1; i <= ndata; i++) { C1[i] = cosf( Alf[1]*t ); S1[i] = sinf( Alf[1]*t ); C2[i] = cosf( Alf[2]*t ); S2[i] = sinf( Alf[2]*t ); PHI(0,i - 1) = 1.0e0; PHI(1,i - 1) = C1[i]; PHI(2,i - 1) = S1[i]; PHI(3,i - 1) = C2[i]; PHI(4,i - 1) = S2[i]; t += del; } return; #undef PHI } /* end of function */ /* ================================================================== */ /* PARAMETER translations */ #define NDER 4 /* end of PARAMETER translations */ void /*FUNCTION*/ scalcb( long ndata, long na, long nb, float alf[], long *ncount, float *der) { #define DER(I_,J_) (*(der+(I_)*(ndata)+(J_))) long int i; float del, t; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Alf = &alf[0] - 1; float *const C1 = &skey.c1[0] - 1; float *const C2 = &skey.c2[0] - 1; float *const S1 = &skey.s1[0] - 1; float *const S2 = &skey.s2[0] - 1; /* end of OFFSET VECTORS */ /* Test case for separable nonlinear least squares computation. * Computes the NDATA x NDER matrix DER. Here NDER is the number of * ones in the indicator array IND(). The columns of DER correspond * to nonzero entries of IND() traversed columnwise. * In this example NDER is 4 and the correspondence is as follows: * Col of DER: 1 2 3 4 * Element of IND(): (2,1) (3,1) (4,2) (5,2) * In this example Row I of DER will be set to contain the values * of the four partial derivatives: * Partial of PHI(I,2) with respect to ALP(1) * Partial of PHI(I,3) with respect to ALP(1) * Partial of PHI(I,4) with respect to ALP(2) * Partial of PHI(I,5) with respect to ALP(2) * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ t = 0.0e0; del = 1.0e0/29.0e0; if (*ncount == skey.key) { for (i = 1; i <= ndata; i++) { DER(0,i - 1) = -S1[i]*t; DER(1,i - 1) = C1[i]*t; DER(2,i - 1) = -S2[i]*t; DER(3,i - 1) = C2[i]*t; t += del; } } else { for (i = 1; i <= ndata; i++) { DER(0,i - 1) = -sinf( Alf[1]*t )*t; DER(1,i - 1) = cosf( Alf[1]*t )*t; DER(2,i - 1) = -sinf( Alf[2]*t )*t; DER(3,i - 1) = cosf( Alf[2]*t )*t; t += del; } } return; #undef DER } /* end of function */