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

   1 /* eigen/genhermv.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 and eigenvectors of a complex
  34  * generalized hermitian-definite eigensystem A x = \lambda B x, where
  35  * A and B are hermitian, and B is positive-definite.
  36  */
  37
  38 static void genhermv_normalize_eigenvectors(gsl_matrix_complex *evec);
  39
  40 /*
  41 gsl_eigen_genhermv_alloc()
  42
  43 Allocate a workspace for solving the generalized hermitian-definite
  44 eigenvalue problem. The size of this workspace is O(5n).
  45
  46 Inputs: n - size of matrices
  47
  48 Return: pointer to workspace
  49 */
  50
  51 gsl_eigen_genhermv_workspace *
  52 gsl_eigen_genhermv_alloc(const size_t n)
  53 {
  54   gsl_eigen_genhermv_workspace *w;
  55
  56   if (n == 0)
  57     {
  58       GSL_ERROR_NULL ("matrix dimension must be positive integer",
  59                       GSL_EINVAL);
  60     }
  61
  62   w = calloc (1, sizeof (gsl_eigen_genhermv_workspace));
  63
  64   if (w == 0)
  65     {
  66       GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
  67     }
  68
  69   w->size = n;
  70
  71   w->hermv_workspace_p = gsl_eigen_hermv_alloc(n);
  72   if (!w->hermv_workspace_p)
  73     {
  74       gsl_eigen_genhermv_free(w);
  75       GSL_ERROR_NULL("failed to allocate space for hermv workspace", GSL_ENOMEM);
  76     }
  77
  78   return (w);
  79 } /* gsl_eigen_genhermv_alloc() */
  80
  81 /*
  82 gsl_eigen_genhermv_free()
  83   Free workspace w
  84 */
  85
  86 void
  87 gsl_eigen_genhermv_free (gsl_eigen_genhermv_workspace * w)
  88 {
  89   if (w->hermv_workspace_p)
  90     gsl_eigen_hermv_free(w->hermv_workspace_p);
  91
  92   free(w);
  93 } /* gsl_eigen_genhermv_free() */
  94
  95 /*
  96 gsl_eigen_genhermv()
  97
  98 Solve the generalized hermitian-definite eigenvalue problem
  99
 100 A x = \lambda B x
 101
 102 for the eigenvalues \lambda and eigenvectors x.
 103
 104 Inputs: A    - complex hermitian matrix
 105         B    - complex hermitian and positive definite matrix
 106         eval - where to store eigenvalues
 107         evec - where to store eigenvectors
 108         w    - workspace
 109
 110 Return: success or error
 111 */
 112
 113 int
 114 gsl_eigen_genhermv (gsl_matrix_complex * A, gsl_matrix_complex * B,
 115                     gsl_vector * eval, gsl_matrix_complex * evec,
 116                     gsl_eigen_genhermv_workspace * w)
 117 {
 118   const size_t N = A->size1;
 119
 120   /* check matrix and vector sizes */
 121
 122   if (N != A->size2)
 123     {
 124       GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
 125     }
 126   else if ((N != B->size1) || (N != B->size2))
 127     {
 128       GSL_ERROR ("B matrix dimensions must match A", GSL_EBADLEN);
 129     }
 130   else if (eval->size != N)
 131     {
 132       GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
 133     }
 134   else if (evec->size1 != evec->size2)
 135     {
 136       GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
 137     }
 138   else if (evec->size1 != N)
 139     {
 140       GSL_ERROR ("eigenvector matrix has wrong size", GSL_EBADLEN);
 141     }
 142   else if (w->size != N)
 143     {
 144       GSL_ERROR ("matrix size does not match workspace", GSL_EBADLEN);
 145     }
 146   else
 147     {
 148       int s;
 149
 150       /* compute Cholesky factorization of B */
 151       s = gsl_linalg_complex_cholesky_decomp(B);
 152       if (s != GSL_SUCCESS)
 153         return s; /* B is not positive definite */
 154
 155       /* transform to standard hermitian eigenvalue problem */
 156       gsl_eigen_genherm_standardize(A, B);
 157
 158       /* compute eigenvalues and eigenvectors */
 159       s = gsl_eigen_hermv(A, eval, evec, w->hermv_workspace_p);
 160       if (s != GSL_SUCCESS)
 161         return s;
 162
 163       /* backtransform eigenvectors: evec -> L^{-H} evec */
 164       gsl_blas_ztrsm(CblasLeft,
 165                      CblasLower,
 166                      CblasConjTrans,
 167                      CblasNonUnit,
 168                      GSL_COMPLEX_ONE,
 169                      B,
 170                      evec);
 171
 172       /* the blas call destroyed the normalization - renormalize */
 173       genhermv_normalize_eigenvectors(evec);
 174
 175       return GSL_SUCCESS;
 176     }
 177 } /* gsl_eigen_genhermv() */
 178
 179 /********************************************
 180  *           INTERNAL ROUTINES              *
 181  ********************************************/
 182
 183 /*
 184 genhermv_normalize_eigenvectors()
 185   Normalize eigenvectors so that their Euclidean norm is 1
 186
 187 Inputs: evec - eigenvectors
 188 */
 189
 190 static void
 191 genhermv_normalize_eigenvectors(gsl_matrix_complex *evec)
 192 {
 193   const size_t N = evec->size1;
 194   size_t i;     /* looping */
 195
 196   for (i = 0; i < N; ++i)
 197     {
 198       gsl_vector_complex_view vi = gsl_matrix_complex_column(evec, i);
 199       double scale = 1.0 / gsl_blas_dznrm2(&vi.vector);
 200
 201       gsl_blas_zdscal(scale, &vi.vector);
 202     }
 203 } /* genhermv_normalize_eigenvectors() */