1 /* driver for hybrj example. */
2
3 #include <stdio.h>
4 #include <math.h>
5 #include <assert.h>
6 #include <minpack.h>
7 #define real __minpack_real__
8
9 void fcn(const int *n, const real *x, real *fvec, real *fjac, const int *ldfjac,
10 int *iflag);
11
main()12 int main()
13 {
14
15 int j, n, ldfjac, maxfev, mode, nprint, info, nfev, njev, lr;
16 real xtol, factor, fnorm;
17 real x[9], fvec[9], fjac[9*9], diag[9], r[45], qtf[9],
18 wa1[9], wa2[9], wa3[9], wa4[9];
19 int one=1;
20
21 n = 9;
22
23 /* the following starting values provide a rough solution. */
24
25 for (j=1; j<=9; j++) {
26 x[j-1] = -1.;
27 }
28
29 ldfjac = 9;
30 lr = 45;
31
32 /* set xtol to the square root of the machine precision. */
33 /* unless high solutions are required, */
34 /* this is the recommended setting. */
35
36 xtol = sqrt(__minpack_func__(dpmpar)(&one));
37
38 maxfev = 1000;
39 mode = 2;
40 for (j=1; j<=9; j++) {
41 diag[j-1] = 1.;
42 }
43 factor = 1.e2;
44 nprint = 0;
45
46 __minpack_func__(hybrj)(&fcn, &n, x, fvec, fjac, &ldfjac, &xtol, &maxfev, diag,
47 &mode, &factor, &nprint, &info, &nfev, &njev, r, &lr, qtf,
48 wa1, wa2, wa3, wa4);
49 fnorm = __minpack_func__(enorm)(&n, fvec);
50
51 printf(" final l2 norm of the residuals%15.7g\n\n", (double)fnorm);
52 printf(" number of function evaluations%10i\n\n", nfev);
53 printf(" number of jacobian evaluations%10i\n\n", njev);
54 printf(" exit parameter %10i\n\n", info);
55 printf(" final approximate solution\n\n");
56 for (j=1; j<=n; j++) {
57 printf("%s%15.7g", j%3==1?"\n ":"", (double)x[j-1]);
58 }
59 printf("\n");
60 return 0;
61 }
62
fcn(const int * n,const real * x,real * fvec,real * fjac,const int * ldfjac,int * iflag)63 void fcn(const int *n, const real *x, real *fvec, real *fjac, const int *ldfjac,
64 int *iflag)
65 {
66
67 /* subroutine fcn for hybrj example. */
68
69 int j, k;
70 real temp, temp1, temp2;
71 assert(*n == 9);
72
73 if (*iflag == 0) {
74 /* insert print statements here when nprint is positive. */
75 /* if the nprint parameter to lmder is positive, the function is
76 called every nprint iterations with iflag=0, so that the
77 function may perform special operations, such as printing
78 residuals. */
79 return;
80 }
81
82 if (*iflag != 2) {
83 /* compute residuals */
84 for (k=1; k <= *n; k++) {
85 temp = (3 - 2*x[k-1])*x[k-1];
86 temp1 = 0;
87 if (k != 1) {
88 temp1 = x[k-1-1];
89 }
90 temp2 = 0;
91 if (k != *n) {
92 temp2 = x[k+1-1];
93 }
94 fvec[k-1] = temp - temp1 - 2*temp2 + 1;
95 }
96 } else {
97 /* compute Jacobian */
98 for (k = 1; k <= *n; k++) {
99 for (j=1; j <= *n; j++) {
100 fjac[k-1 + *ldfjac*(j-1)] = 0;
101 }
102 fjac[k-1 + *ldfjac*(k-1)] = 3 - 4*x[k-1];
103 if (k != 1) {
104 fjac[k-1 + *ldfjac*(k-1-1)] = -1;
105 }
106 if (k != *n) {
107 fjac[k-1 + *ldfjac*(k+1-1)] = -2;
108 }
109 }
110 }
111 return;
112 }
113
114