/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:17 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_snlagu s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include #include #include "p_snlagu.h" /* program DRSNLAGU *>> 1997-06-18 DRSNLAGU Krogh Changes to improve C portability. *>> 1994-11-02 DRSNLAGU Krogh Changes to use M77CON *>> 1994-09-14 DRSNLAGU 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-06-14 CLL @ JPL *>> 1990-03-28 CLL @ JPL * Demo driver for SNLAGU. A variant of the nonlinear LS code NL2SOL. * SNLAGU requires values of the function and the Jacobian matrix. * ------------------------------------------------------------------ *--S replaces "?": DR?NLAGU, ?NLAGU, ?CALCR, ?CALCJ, ?IVSET, ?KEY * ------------------------------------------------------------------ */ /* PARAMETER translations */ #define COVPRT 14 #define F 10 #define LIV (82 + MC) #define LV (105 + MC*(MDATA + 2*MC + 17) + 2*MDATA) #define MC 7 #define MDATA 30 #define OUTLEV 19 #define SOLPRT 22 #define STATPR 23 #define X0PRT 24 /* end of PARAMETER translations */ int main( ) { long int iv[LIV], nc, ndata; float coef[MC], dof, v[LV]; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Coef = &coef[0] - 1; long *const Iv = &iv[0] - 1; float *const V = &v[0] - 1; /* end of OFFSET VECTORS */ /* ------------------------------------------------------------------ */ ndata = MDATA; nc = MC; Coef[1] = 5.0e0; Coef[2] = 10.0e0; Coef[3] = 0.5e0; Coef[4] = 0.5e0; Coef[5] = 0.5e0; Coef[6] = 0.5e0; Coef[7] = 0.5e0; Iv[1] = 0; printf(" Program DRSNLAGU.. Demo driver for SNLAGU.\n A variant of NL2SOL.\n " " SNLAGU requires values of the function and the Jacobian.\n \n " "Sample problem is a nonlinear curve fit to data.\n " "Model function is C3 + C4 * cosf(C1*t) + C5 * sinf(C1*t) +\n " " C6 * cosf(C2*t) + C7 * sinf(C2*t) + Noise\n " "Data generated using\n (C1, ..., C7) = (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; snlagu( ndata, nc, coef, scalcr, scalcj, iv, LIV, LV, v ); dof = max( ndata - nc, 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*/ scalcr( long ndata, long nc, float c[], long *ncount, float rvec[]) { long int i; float del, t; 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}; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const C = &c[0] - 1; float *const C1 = &skey.c1[0] - 1; float *const C2 = &skey.c2[0] - 1; float *const Rvec = &rvec[0] - 1; float *const S1 = &skey.s1[0] - 1; float *const S2 = &skey.s2[0] - 1; float *const Ydata = &ydata[0] - 1; /* end of OFFSET VECTORS */ /* Function evaluation to test nonlinear least squares computation. * Illustrates saving results in common between SCALCR and SCALCJ * 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( C[1]*t ); S1[i] = sinf( C[1]*t ); C2[i] = cosf( C[2]*t ); S2[i] = sinf( C[2]*t ); Rvec[i] = C[3] + C[4]*C1[i] + C[5]*S1[i] + C[6]*C2[i] + C[7]* S2[i] - Ydata[i]; t += del; } return; } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ scalcj( long ndata, long nc, float c[], long *ncount, float *ajac) { #define AJAC(I_,J_) (*(ajac+(I_)*(ndata)+(J_))) long int i; float del, t; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const C = &c[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 */ /* Jacobian evaluation to test nonlinear least squares computation. * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ t = 0.0e0; del = 1.0e0/29.0e0; if (*ncount == skey.key) { for (i = 1; i <= ndata; i++) { AJAC(0,i - 1) = -C[4]*S1[i]*t + C[5]*C1[i]*t; AJAC(1,i - 1) = -C[6]*S2[i]*t + C[7]*C2[i]*t; AJAC(2,i - 1) = 1.0e0; AJAC(3,i - 1) = C1[i]; AJAC(4,i - 1) = S1[i]; AJAC(5,i - 1) = C2[i]; AJAC(6,i - 1) = S2[i]; t += del; } } else { for (i = 1; i <= ndata; i++) { AJAC(2,i - 1) = 1.0e0; AJAC(3,i - 1) = cosf( C[1]*t ); AJAC(4,i - 1) = sinf( C[1]*t ); AJAC(5,i - 1) = cosf( C[2]*t ); AJAC(6,i - 1) = sinf( C[2]*t ); AJAC(0,i - 1) = -C[4]*AJAC(4,i - 1)*t + C[5]*AJAC(3,i - 1)* t; AJAC(1,i - 1) = -C[6]*AJAC(6,i - 1)*t + C[7]*AJAC(5,i - 1)* t; t += del; } } return; #undef AJAC } /* end of function */