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