htswanalysis/MACS/lib/gsl/gsl-1.11/multifit/lmutil.c

   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);
   6
   7 static double
   8 enorm (const gsl_vector * f)
   9 {
  10   return gsl_blas_dnrm2 (f);
  11 }
  12
  13 static double
  14 scaled_enorm (const gsl_vector * d, const gsl_vector * f)
  15 {
  16   double e2 = 0;
  17   size_t i, n = f->size;
  18   for (i = 0; i < n; i++)
  19     {
  20       double fi = gsl_vector_get (f, i);
  21       double di = gsl_vector_get (d, i);
  22       double u = di * fi;
  23       e2 += u * u;
  24     }
  25   return sqrt (e2);
  26 }
  27
  28 static double
  29 compute_delta (gsl_vector * diag, gsl_vector * x)
  30 {
  31   double Dx = scaled_enorm (diag, x);
  32   double factor = 100;  /* generally recommended value from MINPACK */
  33
  34   return (Dx > 0) ? factor * Dx : factor;
  35 }
  36
  37 static double
  38 compute_actual_reduction (double fnorm, double fnorm1)
  39 {
  40   double actred;
  41
  42   if (0.1 * fnorm1 < fnorm)
  43     {
  44       double u = fnorm1 / fnorm;
  45       actred = 1 - u * u;
  46     }
  47   else
  48     {
  49       actred = -1;
  50     }
  51
  52   return actred;
  53 }
  54
  55 static void
  56 compute_diag (const gsl_matrix * J, gsl_vector * diag)
  57 {
  58   size_t i, j, n = J->size1, p = J->size2;
  59
  60   for (j = 0; j < p; j++)
  61     {
  62       double sum = 0;
  63
  64       for (i = 0; i < n; i++)
  65         {
  66           double Jij = gsl_matrix_get (J, i, j);
  67           sum += Jij * Jij;
  68         }
  69       if (sum == 0)
  70         sum = 1.0;
  71
  72       gsl_vector_set (diag, j, sqrt (sum));
  73     }
  74 }
  75
  76 static void
  77 update_diag (const gsl_matrix * J, gsl_vector * diag)
  78 {
  79   size_t i, j, n = diag->size;
  80
  81   for (j = 0; j < n; j++)
  82     {
  83       double cnorm, diagj, sum = 0;
  84       for (i = 0; i < n; i++)
  85         {
  86           double Jij = gsl_matrix_get (J, i, j);
  87           sum += Jij * Jij;
  88         }
  89       if (sum == 0)
  90         sum = 1.0;
  91
  92       cnorm = sqrt (sum);
  93       diagj = gsl_vector_get (diag, j);
  94
  95       if (cnorm > diagj)
  96         gsl_vector_set (diag, j, cnorm);
  97     }
  98 }
  99
 100 static void
 101 compute_rptdx (const gsl_matrix * r, const gsl_permutation * p,
 102                const gsl_vector * dx, gsl_vector * rptdx)
 103 {
 104   size_t i, j, N = dx->size;
 105
 106   for (i = 0; i < N; i++)
 107     {
 108       double sum = 0;
 109
 110       for (j = i; j < N; j++)
 111         {
 112           size_t pj = gsl_permutation_get (p, j);
 113
 114           sum += gsl_matrix_get (r, i, j) * gsl_vector_get (dx, pj);
 115         }
 116
 117       gsl_vector_set (rptdx, i, sum);
 118     }
 119 }
 120
 121
 122 static void
 123 compute_trial_step (gsl_vector * x, gsl_vector * dx, gsl_vector * x_trial)
 124 {
 125   size_t i, N = x->size;
 126
 127   for (i = 0; i < N; i++)
 128     {
 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);
 132     }
 133 }
 134