htswanalysis/MACS/lib/gsl/gsl-1.11/eigen/qrstep.c

   1 /* eigen/qrstep.c
   2  *
   3  * Copyright (C) 2007 Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* remove off-diagonal elements which are neglegible compared with the
  21    neighboring diagonal elements */
  22
  23 static void
  24 chop_small_elements (const size_t N, const double d[], double sd[])
  25 {
  26   double d_i = d[0];
  27
  28   size_t i;
  29
  30   for (i = 0; i < N - 1; i++)
  31     {
  32       double sd_i = sd[i];
  33       double d_ip1 = d[i + 1];
  34
  35       if (fabs (sd_i) < GSL_DBL_EPSILON * (fabs (d_i) + fabs (d_ip1)))
  36         {
  37           sd[i] = 0.0;
  38         }
  39       d_i = d_ip1;
  40     }
  41 }
  42
  43 /* Generate a Givens rotation (cos,sin) which takes v=(x,y) to (|v|,0)
  44
  45    From Golub and Van Loan, "Matrix Computations", Section 5.1.8 */
  46
  47 inline static void
  48 create_givens (const double a, const double b, double *c, double *s)
  49 {
  50   if (b == 0)
  51     {
  52       *c = 1;
  53       *s = 0;
  54     }
  55   else if (fabs (b) > fabs (a))
  56     {
  57       double t = -a / b;
  58       double s1 = 1.0 / sqrt (1 + t * t);
  59       *s = s1;
  60       *c = s1 * t;
  61     }
  62   else
  63     {
  64       double t = -b / a;
  65       double c1 = 1.0 / sqrt (1 + t * t);
  66       *c = c1;
  67       *s = c1 * t;
  68     }
  69 }
  70
  71 inline static double
  72 trailing_eigenvalue (const size_t n, const double d[], const double sd[])
  73 {
  74   double ta = d[n - 2];
  75   double tb = d[n - 1];
  76   double tab = sd[n - 2];
  77
  78   double dt = (ta - tb) / 2.0;
  79
  80   double mu;
  81
  82   if (dt > 0)
  83     {
  84       mu = tb - tab * (tab / (dt + hypot (dt, tab)));
  85     }
  86   else if (dt == 0)
  87     {
  88       mu = tb - fabs(tab);
  89     }
  90   else
  91     {
  92       mu = tb + tab * (tab / ((-dt) + hypot (dt, tab)));
  93     }
  94
  95   return mu;
  96 }
  97
  98 static void
  99 qrstep (const size_t n, double d[], double sd[], double gc[], double gs[])
 100 {
 101   double x, z;
 102   double ak, bk, zk, ap, bp, aq, bq;
 103   size_t k;
 104
 105   double mu = trailing_eigenvalue (n, d, sd);
 106
 107   x = d[0] - mu;
 108   z = sd[0];
 109
 110   ak = 0;
 111   bk = 0;
 112   zk = 0;
 113
 114   ap = d[0];
 115   bp = sd[0];
 116
 117   aq = d[1];
 118
 119   if (n == 2)
 120     {
 121       double c, s;
 122       create_givens (x, z, &c, &s);
 123
 124       if (gc != NULL)
 125         gc[0] = c;
 126       if (gs != NULL)
 127         gs[0] = s;
 128
 129       {
 130         double ap1 = c * (c * ap - s * bp) + s * (s * aq - c * bp);
 131         double bp1 = c * (s * ap + c * bp) - s * (s * bp + c * aq);
 132
 133         double aq1 = s * (s * ap + c * bp) + c * (s * bp + c * aq);
 134
 135         ak = ap1;
 136         bk = bp1;
 137
 138         ap = aq1;
 139       }
 140
 141       d[0] = ak;
 142       sd[0] = bk;
 143       d[1] = ap;
 144
 145       return;
 146     }
 147
 148   bq = sd[1];
 149
 150   for (k = 0; k < n - 1; k++)
 151     {
 152       double c, s;
 153       create_givens (x, z, &c, &s);
 154
 155       /* store Givens rotation */
 156       if (gc != NULL)
 157         gc[k] = c;
 158       if (gs != NULL)
 159         gs[k] = s;
 160
 161       /* compute G' T G */
 162
 163       {
 164         double bk1 = c * bk - s * zk;
 165
 166         double ap1 = c * (c * ap - s * bp) + s * (s * aq - c * bp);
 167         double bp1 = c * (s * ap + c * bp) - s * (s * bp + c * aq);
 168         double zp1 = -s * bq;
 169
 170         double aq1 = s * (s * ap + c * bp) + c * (s * bp + c * aq);
 171         double bq1 = c * bq;
 172
 173         ak = ap1;
 174         bk = bp1;
 175         zk = zp1;
 176
 177         ap = aq1;
 178         bp = bq1;
 179
 180         if (k < n - 2)
 181           aq = d[k + 2];
 182         if (k < n - 3)
 183           bq = sd[k + 2];
 184
 185         d[k] = ak;
 186
 187         if (k > 0)
 188           sd[k - 1] = bk1;
 189
 190         if (k < n - 2)
 191           sd[k + 1] = bp;
 192
 193         x = bk;
 194         z = zk;
 195       }
 196     }
 197
 198   /* k = n - 1 */
 199   d[k] = ap;
 200   sd[k - 1] = bk;
 201 }