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

   1 /* eigen/genherm.c
   2  *
   3  * Copyright (C) 2007 Patrick Alken
   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 <stdlib.h>
  21
  22 #include <config.h>
  23 #include <gsl/gsl_eigen.h>
  24 #include <gsl/gsl_linalg.h>
  25 #include <gsl/gsl_math.h>
  26 #include <gsl/gsl_blas.h>
  27 #include <gsl/gsl_vector.h>
  28 #include <gsl/gsl_matrix.h>
  29 #include <gsl/gsl_complex.h>
  30 #include <gsl/gsl_complex_math.h>
  31
  32 /*
  33  * This module computes the eigenvalues of a complex generalized
  34  * hermitian-definite eigensystem A x = \lambda B x, where A and
  35  * B are hermitian, and B is positive-definite.
  36  */
  37
  38 /*
  39 gsl_eigen_genherm_alloc()
  40
  41 Allocate a workspace for solving the generalized hermitian-definite
  42 eigenvalue problem. The size of this workspace is O(3n).
  43
  44 Inputs: n - size of matrices
  45
  46 Return: pointer to workspace
  47 */
  48
  49 gsl_eigen_genherm_workspace *
  50 gsl_eigen_genherm_alloc(const size_t n)
  51 {
  52   gsl_eigen_genherm_workspace *w;
  53
  54   if (n == 0)
  55     {
  56       GSL_ERROR_NULL ("matrix dimension must be positive integer",
  57                       GSL_EINVAL);
  58     }
  59
  60   w = calloc (1, sizeof (gsl_eigen_genherm_workspace));
  61
  62   if (w == 0)
  63     {
  64       GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
  65     }
  66
  67   w->size = n;
  68
  69   w->herm_workspace_p = gsl_eigen_herm_alloc(n);
  70   if (!w->herm_workspace_p)
  71     {
  72       gsl_eigen_genherm_free(w);
  73       GSL_ERROR_NULL("failed to allocate space for herm workspace", GSL_ENOMEM);
  74     }
  75
  76   return (w);
  77 } /* gsl_eigen_genherm_alloc() */
  78
  79 /*
  80 gsl_eigen_genherm_free()
  81   Free workspace w
  82 */
  83
  84 void
  85 gsl_eigen_genherm_free (gsl_eigen_genherm_workspace * w)
  86 {
  87   if (w->herm_workspace_p)
  88     gsl_eigen_herm_free(w->herm_workspace_p);
  89
  90   free(w);
  91 } /* gsl_eigen_genherm_free() */
  92
  93 /*
  94 gsl_eigen_genherm()
  95
  96 Solve the generalized hermitian-definite eigenvalue problem
  97
  98 A x = \lambda B x
  99
 100 for the eigenvalues \lambda.
 101
 102 Inputs: A    - complex hermitian matrix
 103         B    - complex hermitian and positive definite matrix
 104         eval - where to store eigenvalues
 105         w    - workspace
 106
 107 Return: success or error
 108 */
 109
 110 int
 111 gsl_eigen_genherm (gsl_matrix_complex * A, gsl_matrix_complex * B,
 112                    gsl_vector * eval, gsl_eigen_genherm_workspace * w)
 113 {
 114   const size_t N = A->size1;
 115
 116   /* check matrix and vector sizes */
 117
 118   if (N != A->size2)
 119     {
 120       GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
 121     }
 122   else if ((N != B->size1) || (N != B->size2))
 123     {
 124       GSL_ERROR ("B matrix dimensions must match A", GSL_EBADLEN);
 125     }
 126   else if (eval->size != N)
 127     {
 128       GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
 129     }
 130   else if (w->size != N)
 131     {
 132       GSL_ERROR ("matrix size does not match workspace", GSL_EBADLEN);
 133     }
 134   else
 135     {
 136       int s;
 137
 138       /* compute Cholesky factorization of B */
 139       s = gsl_linalg_complex_cholesky_decomp(B);
 140       if (s != GSL_SUCCESS)
 141         return s; /* B is not positive definite */
 142
 143       /* transform to standard hermitian eigenvalue problem */
 144       gsl_eigen_genherm_standardize(A, B);
 145
 146       s = gsl_eigen_herm(A, eval, w->herm_workspace_p);
 147
 148       return s;
 149     }
 150 } /* gsl_eigen_genherm() */
 151
 152 /*
 153 gsl_eigen_genherm_standardize()
 154   Reduce the generalized hermitian-definite eigenproblem to
 155 the standard hermitian eigenproblem by computing
 156
 157 C = L^{-1} A L^{-H}
 158
 159 where L L^H is the Cholesky decomposition of B
 160
 161 Inputs: A - (input/output) complex hermitian matrix
 162         B - complex hermitian, positive definite matrix in Cholesky form
 163
 164 Return: success
 165
 166 Notes: A is overwritten by L^{-1} A L^{-H}
 167 */
 168
 169 int
 170 gsl_eigen_genherm_standardize(gsl_matrix_complex *A,
 171                               const gsl_matrix_complex *B)
 172 {
 173   const size_t N = A->size1;
 174   size_t i;
 175   double a, b;
 176   gsl_complex y, z;
 177
 178   GSL_SET_IMAG(&z, 0.0);
 179
 180   for (i = 0; i < N; ++i)
 181     {
 182       /* update lower triangle of A(i:n, i:n) */
 183
 184       y = gsl_matrix_complex_get(A, i, i);
 185       a = GSL_REAL(y);
 186       y = gsl_matrix_complex_get(B, i, i);
 187       b = GSL_REAL(y);
 188       a /= b * b;
 189       GSL_SET_REAL(&z, a);
 190       gsl_matrix_complex_set(A, i, i, z);
 191
 192       if (i < N - 1)
 193         {
 194           gsl_vector_complex_view ai =
 195             gsl_matrix_complex_subcolumn(A, i, i + 1, N - i - 1);
 196           gsl_matrix_complex_view ma =
 197             gsl_matrix_complex_submatrix(A, i + 1, i + 1, N - i - 1, N - i - 1);
 198           gsl_vector_complex_const_view bi =
 199             gsl_matrix_complex_const_subcolumn(B, i, i + 1, N - i - 1);
 200           gsl_matrix_complex_const_view mb =
 201             gsl_matrix_complex_const_submatrix(B, i + 1, i + 1, N - i - 1, N - i - 1);
 202
 203           gsl_blas_zdscal(1.0 / b, &ai.vector);
 204
 205           GSL_SET_REAL(&z, -0.5 * a);
 206           gsl_blas_zaxpy(z, &bi.vector, &ai.vector);
 207
 208           gsl_blas_zher2(CblasLower,
 209                          GSL_COMPLEX_NEGONE,
 210                          &ai.vector,
 211                          &bi.vector,
 212                          &ma.matrix);
 213
 214           gsl_blas_zaxpy(z, &bi.vector, &ai.vector);
 215
 216           gsl_blas_ztrsv(CblasLower,
 217                          CblasNoTrans,
 218                          CblasNonUnit,
 219                          &mb.matrix,
 220                          &ai.vector);
 221         }
 222     }
 223
 224   return GSL_SUCCESS;
 225 } /* gsl_eigen_genherm_standardize() */