/*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_dnlsgb s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include #include #include "p_dnlsgb.h" /* program DRDNLSGB *>> 2001-05-24 DRDNLSGB Krogh Minor change for making .f90 version. *>> 1997-06-18 DRDNLSGB Krogh Changes to improve C portability. *>> 1994-11-02 DRDNLSGB Krogh Changes to use M77CON *>> 1994-09-14 DRDNLSGB 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 DNLSGB. A variant of the nonlinear LS code NL2SOL. * DNLSGB solves the "separable" problem. * DNLSGB handles bounds on the nonlinear variables. * DNLSGB requires values of the function and the Jacobian matrix. * Note: MDER is the number of ones in the array IND(). * ------------------------------------------------------------------ *--D replaces "?": DR?NLSGB, ?NLSGB, ?CALCA, ?CALCB, ?IVSET, ?KEY * ------------------------------------------------------------------ */ /* PARAMETER translations */ #define F 10 #define LIV (115 + 4*MA + MB + 2*MDER) #define LV (105 + MDATA*(MB + MDER + MA + 3) + (MB*(MB + 3))/2 + MA*(2*MA + 21)) #define MA 2 #define MB 5 #define MDATA 30 #define MDER 4 #define OUTLEV 19 #define SOLPRT 22 #define STATPR 23 #define X0PRT 24 /* end of PARAMETER translations */ int main( ) { long int iterm, iv[LIV], ivar, j, na, nb, ndata, _i, _r; static long int ind[MA][MB + 1]; double alf[MA], bet[MB], dof, v[LV]; static double bnd[MA][2]; 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}; static int _aini = 1; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Alf = &alf[0] - 1; double *const Bet = &bet[0] - 1; long *const Iv = &iv[0] - 1; double *const V = &v[0] - 1; double *const Ydata = &ydata[0] - 1; /* end of OFFSET VECTORS */ if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */ for (j = 1; j <= 2; j++) { bnd[j - 1][0] = 5.0e0; } for (j = 1; j <= 2; j++) { bnd[j - 1][1] = 10.0e0; } { 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 DRDNLSGB.. Demo driver for DNLSGB.\n A variant of NL2SOL.\n " " DNLSGB handles the Separable problem.\n " " DNLSGB requires values of the function and the Jacobian.\n " " DNLSGB handles bounds on the nonlinear variables.\n \n " "Sample problem is a nonlinear curve fit to data.\n " "Model function is B1 + B2 * cos(A1*t) + B3 * sin(A1*t) +\n " " B4 * cos(A2*t) + B5 * sin(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"); divset( 1, iv, LIV, LV, v ); Iv[X0PRT] = 1; Iv[OUTLEV] = 0; Iv[STATPR] = 1; Iv[SOLPRT] = 1; dnlsgb( ndata, na, nb, alf, bnd, bet, ydata, dcalca, dcalcb, (long*)ind, nb + 1, iv, LIV, LV, v ); dof = max( ndata - na - nb, 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*/ dcalca( long ndata, long na, long nb, double alf[], long *ncount, double *phi) { #define PHI(I_,J_) (*(phi+(I_)*(ndata)+(J_))) long int i; double del, t; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Alf = &alf[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 */ /* 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 DCALCA and * DCALCB 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( Alf[1]*t ); S1[i] = sin( Alf[1]*t ); C2[i] = cos( Alf[2]*t ); S2[i] = sin( 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*/ dcalcb( long ndata, long na, long nb, double alf[], long *ncount, double *der) { #define DER(I_,J_) (*(der+(I_)*(ndata)+(J_))) long int i; double del, t; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Alf = &alf[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 */ /* 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 == dkey.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) = -sin( Alf[1]*t )*t; DER(1,i - 1) = cos( Alf[1]*t )*t; DER(2,i - 1) = -sin( Alf[2]*t )*t; DER(3,i - 1) = cos( Alf[2]*t )*t; t += del; } } return; #undef DER } /* end of function */