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
|
static void compute_diag (const gsl_matrix * J, gsl_vector * diag);
static void update_diag (const gsl_matrix * J, gsl_vector * diag);
static double compute_delta (gsl_vector * diag, gsl_vector * x);
static double scaled_enorm (const gsl_vector * d, const gsl_vector * f);
static double enorm (const gsl_vector * f);
static double
enorm (const gsl_vector * f)
{
return gsl_blas_dnrm2 (f);
}
static double
scaled_enorm (const gsl_vector * d, const gsl_vector * f)
{
double e2 = 0;
size_t i, n = f->size;
for (i = 0; i < n; i++)
{
double fi = gsl_vector_get (f, i);
double di = gsl_vector_get (d, i);
double u = di * fi;
e2 += u * u;
}
return sqrt (e2);
}
static double
compute_delta (gsl_vector * diag, gsl_vector * x)
{
double Dx = scaled_enorm (diag, x);
double factor = 100; /* generally recommended value from MINPACK */
return (Dx > 0) ? factor * Dx : factor;
}
static double
compute_actual_reduction (double fnorm, double fnorm1)
{
double actred;
if (0.1 * fnorm1 < fnorm)
{
double u = fnorm1 / fnorm;
actred = 1 - u * u;
}
else
{
actred = -1;
}
return actred;
}
static void
compute_diag (const gsl_matrix * J, gsl_vector * diag)
{
size_t i, j, n = J->size1, p = J->size2;
for (j = 0; j < p; j++)
{
double sum = 0;
for (i = 0; i < n; i++)
{
double Jij = gsl_matrix_get (J, i, j);
sum += Jij * Jij;
}
if (sum == 0)
sum = 1.0;
gsl_vector_set (diag, j, sqrt (sum));
}
}
static void
update_diag (const gsl_matrix * J, gsl_vector * diag)
{
size_t i, j, n = diag->size;
for (j = 0; j < n; j++)
{
double cnorm, diagj, sum = 0;
for (i = 0; i < n; i++)
{
double Jij = gsl_matrix_get (J, i, j);
sum += Jij * Jij;
}
if (sum == 0)
sum = 1.0;
cnorm = sqrt (sum);
diagj = gsl_vector_get (diag, j);
if (cnorm > diagj)
gsl_vector_set (diag, j, cnorm);
}
}
static void
compute_rptdx (const gsl_matrix * r, const gsl_permutation * p,
const gsl_vector * dx, gsl_vector * rptdx)
{
size_t i, j, N = dx->size;
for (i = 0; i < N; i++)
{
double sum = 0;
for (j = i; j < N; j++)
{
size_t pj = gsl_permutation_get (p, j);
sum += gsl_matrix_get (r, i, j) * gsl_vector_get (dx, pj);
}
gsl_vector_set (rptdx, i, sum);
}
}
static void
compute_trial_step (gsl_vector * x, gsl_vector * dx, gsl_vector * x_trial)
{
size_t i, N = x->size;
for (i = 0; i < N; i++)
{
double pi = gsl_vector_get (dx, i);
double xi = gsl_vector_get (x, i);
gsl_vector_set (x_trial, i, xi + pi);
}
}
|
