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
|
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 j, p = J->size2;
for (j = 0; j < p; j++)
{
gsl_vector_const_view v = gsl_matrix_const_column(J, j);
double norm = gsl_blas_dnrm2(&v.vector);
if (norm == 0)
norm = 1.0;
gsl_vector_set (diag, j, norm);
}
}
static void
update_diag (const gsl_matrix * J, gsl_vector * diag)
{
size_t j, p = J->size2;
for (j = 0; j < p; j++)
{
gsl_vector_const_view v = gsl_matrix_const_column(J, j);
double norm = gsl_blas_dnrm2(&v.vector);
double *diagj = gsl_vector_ptr(diag, j);
if (norm == 0)
norm = 1.0;
*diagj = GSL_MAX(*diagj, norm);
}
}
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);
}
}
|
