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

   1 /* eigen/hermv.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_complex_math.h>
  26 #include <gsl/gsl_linalg.h>
  27 #include <gsl/gsl_eigen.h>
  28
  29 /* Compute eigenvalues/eigenvectors of complex hermitian matrix using
  30    reduction to real symmetric tridiagonal form, followed by QR
  31    iteration with implicit shifts.
  32
  33    See Golub & Van Loan, "Matrix Computations" (3rd ed), Section 8.3 */
  34
  35 #include "qrstep.c"
  36
  37 gsl_eigen_hermv_workspace *
  38 gsl_eigen_hermv_alloc (const size_t n)
  39 {
  40   gsl_eigen_hermv_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_hermv_workspace *) malloc (sizeof(gsl_eigen_hermv_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       free (w);
  59       GSL_ERROR_NULL ("failed to allocate space for diagonal", GSL_ENOMEM);
  60     }
  61
  62   w->sd = (double *) malloc (n * sizeof (double));
  63
  64   if (w->sd == 0)
  65     {
  66       free (w->d);
  67       free (w);
  68       GSL_ERROR_NULL ("failed to allocate space for subdiagonal", GSL_ENOMEM);
  69     }
  70
  71   w->tau = (double *) malloc (2 * n * sizeof (double));
  72
  73   if (w->tau == 0)
  74     {
  75       free (w->sd);
  76       free (w->d);
  77       free (w);
  78       GSL_ERROR_NULL ("failed to allocate space for tau", GSL_ENOMEM);
  79     }
  80
  81   w->gc = (double *) malloc (n * sizeof (double));
  82
  83   if (w->gc == 0)
  84     {
  85       free (w->tau);
  86       free (w->sd);
  87       free (w->d);
  88       free (w);
  89       GSL_ERROR_NULL ("failed to allocate space for cosines", GSL_ENOMEM);
  90     }
  91
  92   w->gs = (double *) malloc (n * sizeof (double));
  93
  94   if (w->gs == 0)
  95     {
  96       free (w->gc);
  97       free (w->tau);
  98       free (w->sd);
  99       free (w->d);
 100       free (w);
 101       GSL_ERROR_NULL ("failed to allocate space for sines", GSL_ENOMEM);
 102     }
 103
 104   w->size = n;
 105
 106   return w;
 107 }
 108
 109 void
 110 gsl_eigen_hermv_free (gsl_eigen_hermv_workspace * w)
 111 {
 112   free (w->gs);
 113   free (w->gc);
 114   free (w->tau);
 115   free (w->sd);
 116   free (w->d);
 117   free (w);
 118 }
 119
 120 int
 121 gsl_eigen_hermv (gsl_matrix_complex * A, gsl_vector * eval,
 122                        gsl_matrix_complex * evec,
 123                        gsl_eigen_hermv_workspace * w)
 124 {
 125   if (A->size1 != A->size2)
 126     {
 127       GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
 128     }
 129   else if (eval->size != A->size1)
 130     {
 131       GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
 132     }
 133   else if (evec->size1 != A->size1 || evec->size2 != A->size1)
 134     {
 135       GSL_ERROR ("eigenvector matrix must match matrix size", GSL_EBADLEN);
 136     }
 137   else
 138     {
 139       const size_t N = A->size1;
 140       double *const d = w->d;
 141       double *const sd = w->sd;
 142
 143       size_t a, b;
 144
 145       /* handle special case */
 146
 147       if (N == 1)
 148         {
 149           gsl_complex A00 = gsl_matrix_complex_get (A, 0, 0);
 150           gsl_vector_set (eval, 0, GSL_REAL(A00));
 151           gsl_matrix_complex_set (evec, 0, 0, GSL_COMPLEX_ONE);
 152           return GSL_SUCCESS;
 153         }
 154
 155       /* Transform the matrix into a symmetric tridiagonal form */
 156
 157       {
 158         gsl_vector_view d_vec = gsl_vector_view_array (d, N);
 159         gsl_vector_view sd_vec = gsl_vector_view_array (sd, N - 1);
 160         gsl_vector_complex_view tau_vec = gsl_vector_complex_view_array (w->tau, N-1);
 161         gsl_linalg_hermtd_decomp (A, &tau_vec.vector);
 162         gsl_linalg_hermtd_unpack (A, &tau_vec.vector, evec, &d_vec.vector, &sd_vec.vector);
 163       }
 164
 165       /* Make an initial pass through the tridiagonal decomposition
 166          to remove off-diagonal elements which are effectively zero */
 167
 168       chop_small_elements (N, d, sd);
 169
 170       /* Progressively reduce the matrix until it is diagonal */
 171
 172       b = N - 1;
 173
 174       while (b > 0)
 175         {
 176           if (sd[b - 1] == 0.0 || isnan(sd[b - 1]))
 177             {
 178               b--;
 179               continue;
 180             }
 181
 182           /* Find the largest unreduced block (a,b) starting from b
 183              and working backwards */
 184
 185           a = b - 1;
 186
 187           while (a > 0)
 188             {
 189               if (sd[a - 1] == 0.0)
 190                 {
 191                   break;
 192                 }
 193               a--;
 194             }
 195
 196           {
 197             size_t i;
 198             const size_t n_block = b - a + 1;
 199             double *d_block = d + a;
 200             double *sd_block = sd + a;
 201             double * const gc = w->gc;
 202             double * const gs = w->gs;
 203
 204             /* apply QR reduction with implicit deflation to the
 205                unreduced block */
 206
 207             qrstep (n_block, d_block, sd_block, gc, gs);
 208
 209             /* Apply  Givens rotation Gij(c,s) to matrix Q,  Q <- Q G */
 210
 211             for (i = 0; i < n_block - 1; i++)
 212               {
 213                 const double c = gc[i], s = gs[i];
 214                 size_t k;
 215
 216                 for (k = 0; k < N; k++)
 217                   {
 218                     gsl_complex qki = gsl_matrix_complex_get (evec, k, a + i);
 219                     gsl_complex qkj = gsl_matrix_complex_get (evec, k, a + i + 1);
 220                     /* qki <= qki * c - qkj * s */
 221                     /* qkj <= qki * s + qkj * c */
 222                     gsl_complex x1 = gsl_complex_mul_real(qki, c);
 223                     gsl_complex y1 = gsl_complex_mul_real(qkj, -s);
 224
 225                     gsl_complex x2 = gsl_complex_mul_real(qki, s);
 226                     gsl_complex y2 = gsl_complex_mul_real(qkj, c);
 227
 228                     gsl_complex qqki = gsl_complex_add(x1, y1);
 229                     gsl_complex qqkj = gsl_complex_add(x2, y2);
 230
 231                     gsl_matrix_complex_set (evec, k, a + i, qqki);
 232                     gsl_matrix_complex_set (evec, k, a + i + 1, qqkj);
 233                   }
 234               }
 235
 236             /* remove any small off-diagonal elements */
 237
 238             chop_small_elements (n_block, d_block, sd_block);
 239           }
 240         }
 241
 242       {
 243         gsl_vector_view d_vec = gsl_vector_view_array (d, N);
 244         gsl_vector_memcpy (eval, &d_vec.vector);
 245       }
 246
 247       return GSL_SUCCESS;
 248     }
 249 }