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

   1 /* eigen/gensymmv.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
  30 /*
  31  * This module computes the eigenvalues and eigenvectors of a real
  32  * generalized symmetric-definite eigensystem A x = \lambda B x, where
  33  * A and B are symmetric, and B is positive-definite.
  34  */
  35
  36 static void gensymmv_normalize_eigenvectors(gsl_matrix *evec);
  37
  38 /*
  39 gsl_eigen_gensymmv_alloc()
  40
  41 Allocate a workspace for solving the generalized symmetric-definite
  42 eigenvalue problem. The size of this workspace is O(4n).
  43
  44 Inputs: n - size of matrices
  45
  46 Return: pointer to workspace
  47 */
  48
  49 gsl_eigen_gensymmv_workspace *
  50 gsl_eigen_gensymmv_alloc(const size_t n)
  51 {
  52   gsl_eigen_gensymmv_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_gensymmv_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->symmv_workspace_p = gsl_eigen_symmv_alloc(n);
  70   if (!w->symmv_workspace_p)
  71     {
  72       gsl_eigen_gensymmv_free(w);
  73       GSL_ERROR_NULL("failed to allocate space for symmv workspace", GSL_ENOMEM);
  74     }
  75
  76   return (w);
  77 } /* gsl_eigen_gensymmv_alloc() */
  78
  79 /*
  80 gsl_eigen_gensymmv_free()
  81   Free workspace w
  82 */
  83
  84 void
  85 gsl_eigen_gensymmv_free (gsl_eigen_gensymmv_workspace * w)
  86 {
  87   if (w->symmv_workspace_p)
  88     gsl_eigen_symmv_free(w->symmv_workspace_p);
  89
  90   free(w);
  91 } /* gsl_eigen_gensymmv_free() */
  92
  93 /*
  94 gsl_eigen_gensymmv()
  95
  96 Solve the generalized symmetric-definite eigenvalue problem
  97
  98 A x = \lambda B x
  99
 100 for the eigenvalues \lambda and eigenvectors x.
 101
 102 Inputs: A    - real symmetric matrix
 103         B    - real symmetric and positive definite matrix
 104         eval - where to store eigenvalues
 105         evec - where to store eigenvectors
 106         w    - workspace
 107
 108 Return: success or error
 109 */
 110
 111 int
 112 gsl_eigen_gensymmv (gsl_matrix * A, gsl_matrix * B, gsl_vector * eval,
 113                     gsl_matrix * evec, gsl_eigen_gensymmv_workspace * w)
 114 {
 115   const size_t N = A->size1;
 116
 117   /* check matrix and vector sizes */
 118
 119   if (N != A->size2)
 120     {
 121       GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
 122     }
 123   else if ((N != B->size1) || (N != B->size2))
 124     {
 125       GSL_ERROR ("B matrix dimensions must match A", GSL_EBADLEN);
 126     }
 127   else if (eval->size != N)
 128     {
 129       GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
 130     }
 131   else if (evec->size1 != evec->size2)
 132     {
 133       GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
 134     }
 135   else if (evec->size1 != N)
 136     {
 137       GSL_ERROR ("eigenvector matrix has wrong size", GSL_EBADLEN);
 138     }
 139   else if (w->size != N)
 140     {
 141       GSL_ERROR ("matrix size does not match workspace", GSL_EBADLEN);
 142     }
 143   else
 144     {
 145       int s;
 146
 147       /* compute Cholesky factorization of B */
 148       s = gsl_linalg_cholesky_decomp(B);
 149       if (s != GSL_SUCCESS)
 150         return s; /* B is not positive definite */
 151
 152       /* transform to standard symmetric eigenvalue problem */
 153       gsl_eigen_gensymm_standardize(A, B);
 154
 155       /* compute eigenvalues and eigenvectors */
 156       s = gsl_eigen_symmv(A, eval, evec, w->symmv_workspace_p);
 157       if (s != GSL_SUCCESS)
 158         return s;
 159
 160       /* backtransform eigenvectors: evec -> L^{-T} evec */
 161       gsl_blas_dtrsm(CblasLeft,
 162                      CblasLower,
 163                      CblasTrans,
 164                      CblasNonUnit,
 165                      1.0,
 166                      B,
 167                      evec);
 168
 169       /* the blas call destroyed the normalization - renormalize */
 170       gensymmv_normalize_eigenvectors(evec);
 171
 172       return GSL_SUCCESS;
 173     }
 174 } /* gsl_eigen_gensymmv() */
 175
 176 /********************************************
 177  *           INTERNAL ROUTINES              *
 178  ********************************************/
 179
 180 /*
 181 gensymmv_normalize_eigenvectors()
 182   Normalize eigenvectors so that their Euclidean norm is 1
 183
 184 Inputs: evec - eigenvectors
 185 */
 186
 187 static void
 188 gensymmv_normalize_eigenvectors(gsl_matrix *evec)
 189 {
 190   const size_t N = evec->size1;
 191   size_t i;     /* looping */
 192
 193   for (i = 0; i < N; ++i)
 194     {
 195       gsl_vector_view vi = gsl_matrix_column(evec, i);
 196       double scale = 1.0 / gsl_blas_dnrm2(&vi.vector);
 197
 198       gsl_blas_dscal(scale, &vi.vector);
 199     }
 200 } /* gensymmv_normalize_eigenvectors() */