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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
|
/* Code to test the conjgrad minimiser */
/*
* Copyright 1999, 2018 Graeme W. Gill
* All rights reserved.
*
* This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
* see the License.txt file for licencing details.
*/
#include <stdio.h>
#include "numlib.h"
/* Final approximate solution: */
double expect[9] = {
-0.5706545E+00,
-0.6816283E+00,
-0.7017325E+00,
-0.7042129E+00,
-0.7013690E+00,
-0.6918656E+00,
-0.6657920E+00,
-0.5960342E+00,
-0.4164121E+00 };
double fcn( /* Return function value */
void *fdata, /* Opaque data pointer */
double tp[] /* Multivriate input value */
);
static double dfcn(
void *fdata,
double dp[],
double tp[]
);
#define N 9
static void progress(void *pdata, int perc) {
printf("%c% 3d%%",cr_char,perc);
if (perc == 100)
printf("\n");
fflush(stdout);
}
int main(void)
{
double cp[N]; /* Function input values */
double s[N]; /* Search area */
double err;
int j;
int nprint = 0; /* Itteration debugging print = off */
int rc;
error_program = "tpowell"; /* Set global error reporting string */
check_if_not_interactive();
/* The following starting values provide a rough solution. */
for (j = 0; j < N; j++) {
cp[j] = -1.f;
s[j] = 0.9;
}
nprint = 0;
/* Set tol to the square root of the machine precision. */
/* Unless high precision solutions are required, */
/* this is the recommended setting. */
rc = conjgrad(
&err,
N, /* Dimentionality */
cp, /* Initial starting point */
s, /* Size of initial search area */
0.00000001, /* Tollerance of error change to stop on */
1000, /* Maximum iterations allowed */
fcn, /* Error function to evaluate */
dfcn, /* Partial derivative function */
NULL, /* Opaque data needed by function */
progress, /* Progress callback */
NULL /* Context for callback */
);
fprintf(stdout,"Status = %d, final approximate solution err = %f:\n",rc,err);
for (j = 0; j < N; j++) {
fprintf(stdout,"cp[%d] = %e, expect %e\n",j,cp[j],expect[j]);
}
return 0;
} /* main() */
/* Function being minimized */
double fcn( /* Return function value */
void *fdata, /* Opaque data pointer */
double tp[] /* Multivriate input value */
) {
double err, tt;
double temp, temp1, temp2;
int k;
/* Function Body */
err = 0.0;
for (k = 0; k < N; ++k) {
temp = (3.0 - 2.0 * tp[k]) * tp[k];
temp1 = 0.0;
if (k != 0) {
temp1 = tp[k-1];
}
temp2 = 0.0;
if (k != ((N)-1))
temp2 = tp[k+1];
tt = temp - temp1 - 2.0 * temp2 + 1.0;
err += tt * tt;
}
err = sqrt(err);
//printf("Returning %16.14f\n",err);
return err;
}
/* Compute aprox. forward difference */
#define JEPS 1.0e-8 /* Aprox. sqrt of machine precision */
static double dfcn(
void *fdata,
double dp[],
double tp[]
) {
int i;
double d, h, temp;
d = fcn(fdata, tp);
for (i = 0; i < N; i++) {
temp = tp[i];
h = JEPS * fabs(temp);
if (h == 0.0)
h = JEPS;
tp[i] = temp + h; /* Add delta */
h = tp[i] - temp; /* Actual delta with fp precision limits */
dp[i] = (fcn(fdata, tp) - d)/h;
tp[i] = temp; /* Restore value */
}
return DFUNC_NRV;
}
|
