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

   1 /* eigen/symmv.c
   2  *
   3  * Copyright (C) 2001, 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 #include <config.h>
  21 #include <stdlib.h>
  22 #include <gsl/gsl_math.h>
  23 #include <gsl/gsl_vector.h>
  24 #include <gsl/gsl_matrix.h>
  25 #include <gsl/gsl_linalg.h>
  26 #include <gsl/gsl_eigen.h>
  27
  28 /* Compute eigenvalues/eigenvectors of real symmetric matrix using
  29    reduction to tridiagonal form, followed by QR iteration with
  30    implicit shifts.
  31
  32    See Golub & Van Loan, "Matrix Computations" (3rd ed), Section 8.3
  33    */
  34
  35 #include "qrstep.c"
  36
  37 gsl_eigen_symmv_workspace *
  38 gsl_eigen_symmv_alloc (const size_t n)
  39 {
  40   gsl_eigen_symmv_workspace * w ;
  41
  42   if (n == 0)
  43     {
  44       GSL_ERROR_NULL ("matrix dimension must be positive integer", GSL_EINVAL);
  45     }
  46
  47   w= ((gsl_eigen_symmv_workspace *) malloc (sizeof(gsl_eigen_symmv_workspace)));
  48
  49   if (w == 0)
  50     {
  51       GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
  52     }
  53
  54   w->d = (double *) malloc (n * sizeof (double));
  55
  56   if (w->d == 0)
  57     {
  58       GSL_ERROR_NULL ("failed to allocate space for diagonal", GSL_ENOMEM);
  59     }
  60
  61   w->sd = (double *) malloc (n * sizeof (double));
  62
  63   if (w->sd == 0)
  64     {
  65       GSL_ERROR_NULL ("failed to allocate space for subdiagonal", GSL_ENOMEM);
  66     }
  67
  68   w->gc = (double *) malloc (n * sizeof (double));
  69
  70   if (w->gc == 0)
  71     {
  72       GSL_ERROR_NULL ("failed to allocate space for cosines", GSL_ENOMEM);
  73     }
  74
  75   w->gs = (double *) malloc (n * sizeof (double));
  76
  77   if (w->gs == 0)
  78     {
  79       GSL_ERROR_NULL ("failed to allocate space for sines", GSL_ENOMEM);
  80     }
  81
  82   w->size = n;
  83
  84   return w;
  85 }
  86
  87 void
  88 gsl_eigen_symmv_free (gsl_eigen_symmv_workspace * w)
  89 {
  90   free(w->gs);
  91   free(w->gc);
  92   free(w->sd);
  93   free(w->d);
  94   free(w);
  95 }
  96
  97
  98 int
  99 gsl_eigen_symmv (gsl_matrix * A, gsl_vector * eval, gsl_matrix * evec,
 100                        gsl_eigen_symmv_workspace * w)
 101 {
 102   if (A->size1 != A->size2)
 103     {
 104       GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
 105     }
 106   else if (eval->size != A->size1)
 107     {
 108       GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
 109     }
 110   else if (evec->size1 != A->size1 || evec->size2 != A->size1)
 111     {
 112       GSL_ERROR ("eigenvector matrix must match matrix size", GSL_EBADLEN);
 113     }
 114   else
 115     {
 116       double *const d = w->d;
 117       double *const sd = w->sd;
 118       const size_t N = A->size1;
 119       size_t a, b;
 120
 121       /* handle special case */
 122
 123       if (N == 1)
 124         {
 125           double A00 = gsl_matrix_get (A, 0, 0);
 126           gsl_vector_set (eval, 0, A00);
 127           gsl_matrix_set (evec, 0, 0, 1.0);
 128           return GSL_SUCCESS;
 129         }
 130
 131       /* use sd as the temporary workspace for the decomposition when
 132          computing eigenvectors */
 133
 134       {
 135         gsl_vector_view d_vec = gsl_vector_view_array (d, N);
 136         gsl_vector_view sd_vec = gsl_vector_view_array (sd, N - 1);
 137         gsl_vector_view tau = gsl_vector_view_array (sd, N - 1);
 138         gsl_linalg_symmtd_decomp (A, &tau.vector);
 139         gsl_linalg_symmtd_unpack (A, &tau.vector, evec, &d_vec.vector, &sd_vec.vector);
 140       }
 141
 142       /* Make an initial pass through the tridiagonal decomposition
 143          to remove off-diagonal elements which are effectively zero */
 144
 145       chop_small_elements (N, d, sd);
 146
 147       /* Progressively reduce the matrix until it is diagonal */
 148
 149       b = N - 1;
 150
 151       while (b > 0)
 152         {
 153           if (sd[b - 1] == 0.0 || isnan(sd[b - 1]))
 154             {
 155               b--;
 156               continue;
 157             }
 158
 159           /* Find the largest unreduced block (a,b) starting from b
 160              and working backwards */
 161
 162           a = b - 1;
 163
 164           while (a > 0)
 165             {
 166               if (sd[a - 1] == 0.0)
 167                 {
 168                   break;
 169                 }
 170               a--;
 171             }
 172
 173           {
 174             size_t i;
 175             const size_t n_block = b - a + 1;
 176             double *d_block = d + a;
 177             double *sd_block = sd + a;
 178             double * const gc = w->gc;
 179             double * const gs = w->gs;
 180
 181             /* apply QR reduction with implicit deflation to the
 182                unreduced block */
 183
 184             qrstep (n_block, d_block, sd_block, gc, gs);
 185
 186             /* Apply  Givens rotation Gij(c,s) to matrix Q,  Q <- Q G */
 187
 188             for (i = 0; i < n_block - 1; i++)
 189               {
 190                 const double c = gc[i], s = gs[i];
 191                 size_t k;
 192
 193                 for (k = 0; k < N; k++)
 194                   {
 195                     double qki = gsl_matrix_get (evec, k, a + i);
 196                     double qkj = gsl_matrix_get (evec, k, a + i + 1);
 197                     gsl_matrix_set (evec, k, a + i, qki * c - qkj * s);
 198                     gsl_matrix_set (evec, k, a + i + 1, qki * s + qkj * c);
 199                   }
 200               }
 201
 202             /* remove any small off-diagonal elements */
 203
 204             chop_small_elements (N, d, sd);
 205           }
 206         }
 207
 208       {
 209         gsl_vector_view d_vec = gsl_vector_view_array (d, N);
 210         gsl_vector_memcpy (eval, &d_vec.vector);
 211       }
 212
 213       return GSL_SUCCESS;
 214     }
 215 }