/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:11 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_dnlagb s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include #include #include "p_dnlagb.h" /* program DRDNLAGB *>> 2001-05-24 DRDNLAGB Krogh Minor change for making .f90 version. *>> 1997-06-18 DRDNLAGB Krogh Changes to improve C portability. *>> 1994-11-02 DRDNLAGB Krogh Changes to use M77CON *>> 1994-09-14 DRDNLAGB 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 DNLAGB. A variant of the nonlinear LS code NL2SOL. * DNLAGB handles bounded variables. Requires values of the function * and the Jacobian matrix. * ------------------------------------------------------------------ *--D replaces "?": DR?NLAGB, ?NLAGB, ?CALCR, ?CALCJ, ?IVSET, ?KEY * ------------------------------------------------------------------ */ /* PARAMETER translations */ #define F 10 #define LIV (82 + 4*MC) #define LV (105 + MC*(MDATA + 2*MC + 21) + 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], j, nc, ndata, _i, _r; double coef[MC], dof, v[LV]; static double bnd[MC][2]; static int _aini = 1; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Coef = &coef[0] - 1; long *const Iv = &iv[0] - 1; double *const V = &v[0] - 1; /* end of OFFSET VECTORS */ if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */ { static double _itmp0[] = {5.0e0,5.0e0,0.0e0,0.0e0,0.0e0, 0.0e0,0.0e0}; for (j = 1, _r = 0; j <= 7; j++) { bnd[j - 1][0] = _itmp0[_r++]; } } { static double _itmp1[] = {10.0e0,10.0e0,1.0e0,1.0e0,1.0e0, 1.0e0,1.0e0}; for (j = 1, _r = 0; j <= 7; j++) { bnd[j - 1][1] = _itmp1[_r++]; } } _aini = 0; } /* ------------------------------------------------------------------ */ 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; printf(" Program DRDNLAGB.. Demo driver for DNLAGB.\n A variant of NL2SOL.\n " " DNLAGB requires values of the function and the Jacobian.\n " " DNLAGB handles bounds on variables.\n \n " "Sample problem is a nonlinear curve fit to data.\n " "Model function is C3 + C4 * cos(C1*t) + C5 * sin(C1*t) +\n " " C6 * cos(C2*t) + C7 * sin(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 " "Setting lower bounds = ( 5, 5,0,0,0,0,0).\n " "Setting upper bounds = (10,10,1,1,1,1,1).\n \n"); divset( 1, iv, LIV, LV, v ); Iv[X0PRT] = 1; Iv[OUTLEV] = 0; Iv[STATPR] = 1; Iv[SOLPRT] = 1; dnlagb( ndata, nc, coef, bnd, dcalcr, dcalcj, iv, LIV, LV, v ); dof = max( ndata - nc, 1 ); printf(" \n SIGFAC: sqrt((2 * V(F))/DOF) =%12.4g\n", sqrt( 2.0e0*V[F]/dof )); exit(0); } /* end of function */ /* ================================================================== */ /* COMMON translations */ struct t_dkey { double c1[MDATA], c2[MDATA], s1[MDATA], s2[MDATA]; long int key; } dkey; /* end of COMMON translations */ void /*FUNCTION*/ dcalcr( long ndata, long nc, double c[], long *ncount, double rvec[]) { long int i; double del, t; static double 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 */ double *const C = &c[0] - 1; double *const C1 = &dkey.c1[0] - 1; double *const C2 = &dkey.c2[0] - 1; double *const Rvec = &rvec[0] - 1; double *const S1 = &dkey.s1[0] - 1; double *const S2 = &dkey.s2[0] - 1; double *const Ydata = &ydata[0] - 1; /* end of OFFSET VECTORS */ /* Function evaluation to test nonlinear least squares computation. * Illustrates saving results in common between DCALCR and DCALCJ * to avoid recalculation of common subexpressions. * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ t = 0.0e0; del = 1.0e0/29.0e0; dkey.key = *ncount; for (i = 1; i <= ndata; i++) { C1[i] = cos( C[1]*t ); S1[i] = sin( C[1]*t ); C2[i] = cos( C[2]*t ); S2[i] = sin( 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*/ dcalcj( long ndata, long nc, double c[], long *ncount, double *ajac) { #define AJAC(I_,J_) (*(ajac+(I_)*(ndata)+(J_))) long int i; double del, t; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const C = &c[0] - 1; double *const C1 = &dkey.c1[0] - 1; double *const C2 = &dkey.c2[0] - 1; double *const S1 = &dkey.s1[0] - 1; double *const S2 = &dkey.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 == dkey.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) = cos( C[1]*t ); AJAC(4,i - 1) = sin( C[1]*t ); AJAC(5,i - 1) = cos( C[2]*t ); AJAC(6,i - 1) = sin( 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 */