1 /* driver for lmdif example. */
2
3 #include <stdio.h>
4 #include <math.h>
5 #include <string.h>
6 #include <assert.h>
7 #include <cminpack.h>
8 #define real __cminpack_real__
9
10 #define TEST_COVAR
11
12 /* the following struct defines the data points */
13 typedef struct {
14 int m;
15 real *y;
16 } fcndata_t;
17
18 int fcn(void *p, int m, int n, const real *x, real *fvec, int iflag);
19
main()20 int main()
21 {
22 int i, j, maxfev, mode, nprint, info, nfev, ldfjac;
23 int ipvt[3];
24 real ftol, xtol, gtol, epsfcn, factor, fnorm;
25 real x[3], fvec[15], diag[3], fjac[15*3], qtf[3],
26 wa1[3], wa2[3], wa3[3], wa4[15];
27 int k;
28 const int m = 15;
29 const int n = 3;
30 /* auxiliary data (e.g. measurements) */
31 real y[15] = {1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1,
32 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39};
33 #ifdef TEST_COVAR
34 real covfac;
35 real fjac1[15*3];
36 #endif
37 fcndata_t data;
38 data.m = m;
39 data.y = y;
40
41 /* the following starting values provide a rough fit. */
42
43 x[0] = 1.;
44 x[1] = 1.;
45 x[2] = 1.;
46
47 ldfjac = 15;
48
49 /* set ftol and xtol to the square root of the machine */
50 /* and gtol to zero. unless high solutions are */
51 /* required, these are the recommended settings. */
52
53 ftol = sqrt(__cminpack_func__(dpmpar)(1));
54 xtol = sqrt(__cminpack_func__(dpmpar)(1));
55 gtol = 0.;
56
57 maxfev = 800;
58 epsfcn = 0.;
59 mode = 1;
60 factor = 1.e2;
61 nprint = 0;
62
63 info = __cminpack_func__(lmdif)(fcn, &data, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn,
64 diag, mode, factor, nprint, &nfev, fjac, ldfjac,
65 ipvt, qtf, wa1, wa2, wa3, wa4);
66
67 fnorm = __cminpack_func__(enorm)(m, fvec);
68
69 printf(" final l2 norm of the residuals%15.7g\n\n", (double)fnorm);
70 printf(" number of function evaluations%10i\n\n", nfev);
71 printf(" exit parameter %10i\n\n", info);
72 printf(" final approximate solution\n");
73 for (j=0; j<n; ++j) {
74 printf("%s%15.7g", j%3==0?"\n ":"", (double)x[j]);
75 }
76 printf("\n");
77 ftol = __cminpack_func__(dpmpar)(1);
78 #ifdef TEST_COVAR
79 /* test the original covar from MINPACK */
80 covfac = fnorm*fnorm/(m-n);
81 memcpy(fjac1, fjac, sizeof(fjac));
82 __cminpack_func__(covar)(n, fjac1, ldfjac, ipvt, ftol, wa1);
83 /*
84 printf(" covariance (using covar)\n");
85 for (i=0; i<n; ++i) {
86 for (j=0; j<n; ++j)
87 printf("%s%15.7g", j%3==1?"\n ":"", (double)fjac1[i*ldfjac+j]*covfac);
88 }
89 printf("\n");
90 */
91 #endif
92 /* test covar1, which also estimates the rank of the Jacobian */
93 k = __cminpack_func__(covar1)(m, n, fnorm*fnorm, fjac, ldfjac, ipvt, ftol, wa1);
94 printf(" covariance\n");
95 for (i=0; i<n; ++i) {
96 for (j=0; j<n; ++j)
97 printf("%s%15.7g", j%3==0?"\n ":"", (double)fjac[i*ldfjac+j]);
98 }
99 printf("\n");
100 (void)k;
101 #ifdef TEST_COVAR
102 if (k == n) {
103 /* comparison only works if covariance matrix has full rank */
104 for (i=0; i<n; ++i) {
105 for (j=0; j<n; ++j) {
106 if (fjac[i*ldfjac+j] != fjac1[i*ldfjac+j]*covfac) {
107 printf("component (%d,%d) of covar and covar1 differ: %g != %g\n", i, j, (double)fjac[i*ldfjac+j], (double)(fjac1[i*ldfjac+j]*covfac));
108 }
109 }
110 }
111 }
112 #endif
113 /* printf(" rank(J) = %d\n", k != 0 ? k : n); */
114 return 0;
115 }
116
fcn(void * p,int m,int n,const real * x,real * fvec,int iflag)117 int fcn(void *p, int m, int n, const real *x, real *fvec, int iflag)
118 {
119
120 /* subroutine fcn for lmdif example. */
121
122 int i;
123 real tmp1, tmp2, tmp3;
124 const real *y = ((fcndata_t*)p)->y;
125 assert(m == 15 && n == 3);
126
127 if (iflag == 0) {
128 /* insert print statements here when nprint is positive. */
129 /* if the nprint parameter to lmdif is positive, the function is
130 called every nprint iterations with iflag=0, so that the
131 function may perform special operations, such as printing
132 residuals. */
133 return 0;
134 }
135
136 /* compute residuals */
137 for (i = 0; i < 15; ++i) {
138 tmp1 = i + 1;
139 tmp2 = 15 - i;
140 tmp3 = (i > 7) ? tmp2 : tmp1;
141 fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3));
142 }
143 return 0;
144 }
145