1 /* driver for lmder1 example. */
2
3
4 #include <stdio.h>
5 #include <math.h>
6 #include <assert.h>
7 #include <minpack.h>
8 #define real __minpack_real__
9
10 void fcn(const int *m, const int *n, const real *x, real *fvec, real *fjac,
11 const int *ldfjac, int *iflag);
12
main()13 int main()
14 {
15 int j, m, n, ldfjac, info, lwa;
16 int ipvt[3];
17 real tol, fnorm;
18 real x[3], fvec[15], fjac[15*3], wa[30];
19 int one=1;
20
21 m = 15;
22 n = 3;
23
24 /* the following starting values provide a rough fit. */
25
26 x[1-1] = 1.;
27 x[2-1] = 1.;
28 x[3-1] = 1.;
29
30 ldfjac = 15;
31 lwa = 30;
32
33 /* set tol to the square root of the machine precision. */
34 /* unless high solutions are required, */
35 /* this is the recommended setting. */
36
37 tol = sqrt(__minpack_func__(dpmpar)(&one));
38
39 __minpack_func__(lmder1)(&fcn, &m, &n, x, fvec, fjac, &ldfjac, &tol,
40 &info, ipvt, wa, &lwa);
41 fnorm = __minpack_func__(enorm)(&m, fvec);
42 printf(" final l2 norm of the residuals%15.7g\n\n", (double)fnorm);
43 printf(" exit parameter %10i\n\n", info);
44 printf(" final approximate solution\n");
45 for (j=1; j<=n; j++) {
46 printf("%s%15.7g", j%3==1?"\n ":"", (double)x[j-1]);
47 }
48 printf("\n");
49
50 return 0;
51 }
52
fcn(const int * m,const int * n,const real * x,real * fvec,real * fjac,const int * ldfjac,int * iflag)53 void fcn(const int *m, const int *n, const real *x, real *fvec, real *fjac,
54 const int *ldfjac, int *iflag)
55 {
56
57 /* subroutine fcn for lmder1 example. */
58
59 int i;
60 real tmp1, tmp2, tmp3, tmp4;
61 real y[15] = {1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1,
62 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39};
63 assert(*m == 15 && *n == 3);
64
65 if (*iflag == 0) {
66 /* insert print statements here when nprint is positive. */
67 /* if the nprint parameter to lmder is positive, the function is
68 called every nprint iterations with iflag=0, so that the
69 function may perform special operations, such as printing
70 residuals. */
71 return;
72 }
73
74 if (*iflag != 2) {
75 /* compute residuals */
76 for (i = 1; i <= 15; i++) {
77 tmp1 = i;
78 tmp2 = 16 - i;
79 tmp3 = tmp1;
80 if (i > 8) {
81 tmp3 = tmp2;
82 }
83 fvec[i-1] = y[i-1] - (x[1-1] + tmp1/(x[2-1]*tmp2 + x[3-1]*tmp3));
84 }
85 } else {
86 /* compute Jacobian */
87 for (i = 1; i <= 15; i++) {
88 tmp1 = i;
89 tmp2 = 16 - i;
90 tmp3 = tmp1;
91 if (i > 8) {
92 tmp3 = tmp2;
93 }
94 tmp4 = (x[2-1]*tmp2 + x[3-1]*tmp3); tmp4 = tmp4*tmp4;
95 fjac[i-1 + *ldfjac*(1-1)] = -1.;
96 fjac[i-1 + *ldfjac*(2-1)] = tmp1*tmp2/tmp4;
97 fjac[i-1 + *ldfjac*(3-1)] = tmp1*tmp3/tmp4;
98 }
99 }
100 }
101