1 static void compute_diag (const gsl_matrix * J, gsl_vector * diag);
2 static void update_diag (const gsl_matrix * J, gsl_vector * diag);
3 static double compute_delta (gsl_vector * diag, gsl_vector * x);
4 static double scaled_enorm (const gsl_vector * d, const gsl_vector * f);
5 static double enorm (const gsl_vector * f);
8 enorm (const gsl_vector * f)
10 return gsl_blas_dnrm2 (f);
14 scaled_enorm (const gsl_vector * d, const gsl_vector * f)
17 size_t i, n = f->size;
18 for (i = 0; i < n; i++)
20 double fi = gsl_vector_get (f, i);
21 double di = gsl_vector_get (d, i);
29 compute_delta (gsl_vector * diag, gsl_vector * x)
31 double Dx = scaled_enorm (diag, x);
32 double factor = 100; /* generally recommended value from MINPACK */
34 return (Dx > 0) ? factor * Dx : factor;
38 compute_actual_reduction (double fnorm, double fnorm1)
42 if (0.1 * fnorm1 < fnorm)
44 double u = fnorm1 / fnorm;
56 compute_diag (const gsl_matrix * J, gsl_vector * diag)
58 size_t i, j, n = J->size1, p = J->size2;
60 for (j = 0; j < p; j++)
64 for (i = 0; i < n; i++)
66 double Jij = gsl_matrix_get (J, i, j);
72 gsl_vector_set (diag, j, sqrt (sum));
77 update_diag (const gsl_matrix * J, gsl_vector * diag)
79 size_t i, j, n = diag->size;
81 for (j = 0; j < n; j++)
83 double cnorm, diagj, sum = 0;
84 for (i = 0; i < n; i++)
86 double Jij = gsl_matrix_get (J, i, j);
93 diagj = gsl_vector_get (diag, j);
96 gsl_vector_set (diag, j, cnorm);
101 compute_rptdx (const gsl_matrix * r, const gsl_permutation * p,
102 const gsl_vector * dx, gsl_vector * rptdx)
104 size_t i, j, N = dx->size;
106 for (i = 0; i < N; i++)
110 for (j = i; j < N; j++)
112 size_t pj = gsl_permutation_get (p, j);
114 sum += gsl_matrix_get (r, i, j) * gsl_vector_get (dx, pj);
117 gsl_vector_set (rptdx, i, sum);
123 compute_trial_step (gsl_vector * x, gsl_vector * dx, gsl_vector * x_trial)
125 size_t i, N = x->size;
127 for (i = 0; i < N; i++)
129 double pi = gsl_vector_get (dx, i);
130 double xi = gsl_vector_get (x, i);
131 gsl_vector_set (x_trial, i, xi + pi);