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

   1 /* eigen/gensymm.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 of a real generalized
  32  * symmetric-definite eigensystem A x = \lambda B x, where A and
  33  * B are symmetric, and B is positive-definite.
  34  */
  35
  36 /*
  37 gsl_eigen_gensymm_alloc()
  38
  39 Allocate a workspace for solving the generalized symmetric-definite
  40 eigenvalue problem. The size of this workspace is O(2n).
  41
  42 Inputs: n - size of matrices
  43
  44 Return: pointer to workspace
  45 */
  46
  47 gsl_eigen_gensymm_workspace *
  48 gsl_eigen_gensymm_alloc(const size_t n)
  49 {
  50   gsl_eigen_gensymm_workspace *w;
  51
  52   if (n == 0)
  53     {
  54       GSL_ERROR_NULL ("matrix dimension must be positive integer",
  55                       GSL_EINVAL);
  56     }
  57
  58   w = calloc (1, sizeof (gsl_eigen_gensymm_workspace));
  59
  60   if (w == 0)
  61     {
  62       GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
  63     }
  64
  65   w->size = n;
  66
  67   w->symm_workspace_p = gsl_eigen_symm_alloc(n);
  68   if (!w->symm_workspace_p)
  69     {
  70       gsl_eigen_gensymm_free(w);
  71       GSL_ERROR_NULL("failed to allocate space for symm workspace", GSL_ENOMEM);
  72     }
  73
  74   return (w);
  75 } /* gsl_eigen_gensymm_alloc() */
  76
  77 /*
  78 gsl_eigen_gensymm_free()
  79   Free workspace w
  80 */
  81
  82 void
  83 gsl_eigen_gensymm_free (gsl_eigen_gensymm_workspace * w)
  84 {
  85   if (w->symm_workspace_p)
  86     gsl_eigen_symm_free(w->symm_workspace_p);
  87
  88   free(w);
  89 } /* gsl_eigen_gensymm_free() */
  90
  91 /*
  92 gsl_eigen_gensymm()
  93
  94 Solve the generalized symmetric-definite eigenvalue problem
  95
  96 A x = \lambda B x
  97
  98 for the eigenvalues \lambda.
  99
 100 Inputs: A    - real symmetric matrix
 101         B    - real symmetric and positive definite matrix
 102         eval - where to store eigenvalues
 103         w    - workspace
 104
 105 Return: success or error
 106 */
 107
 108 int
 109 gsl_eigen_gensymm (gsl_matrix * A, gsl_matrix * B, gsl_vector * eval,
 110                    gsl_eigen_gensymm_workspace * w)
 111 {
 112   const size_t N = A->size1;
 113
 114   /* check matrix and vector sizes */
 115
 116   if (N != A->size2)
 117     {
 118       GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
 119     }
 120   else if ((N != B->size1) || (N != B->size2))
 121     {
 122       GSL_ERROR ("B matrix dimensions must match A", GSL_EBADLEN);
 123     }
 124   else if (eval->size != N)
 125     {
 126       GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
 127     }
 128   else if (w->size != N)
 129     {
 130       GSL_ERROR ("matrix size does not match workspace", GSL_EBADLEN);
 131     }
 132   else
 133     {
 134       int s;
 135
 136       /* compute Cholesky factorization of B */
 137       s = gsl_linalg_cholesky_decomp(B);
 138       if (s != GSL_SUCCESS)
 139         return s; /* B is not positive definite */
 140
 141       /* transform to standard symmetric eigenvalue problem */
 142       gsl_eigen_gensymm_standardize(A, B);
 143
 144       s = gsl_eigen_symm(A, eval, w->symm_workspace_p);
 145
 146       return s;
 147     }
 148 } /* gsl_eigen_gensymm() */
 149
 150 /*
 151 gsl_eigen_gensymm_standardize()
 152   Reduce the generalized symmetric-definite eigenproblem to
 153 the standard symmetric eigenproblem by computing
 154
 155 C = L^{-1} A L^{-t}
 156
 157 where L L^t is the Cholesky decomposition of B
 158
 159 Inputs: A - (input/output) real symmetric matrix
 160         B - real symmetric, positive definite matrix in Cholesky form
 161
 162 Return: success
 163
 164 Notes: A is overwritten by L^{-1} A L^{-t}
 165 */
 166
 167 int
 168 gsl_eigen_gensymm_standardize(gsl_matrix *A, const gsl_matrix *B)
 169 {
 170   const size_t N = A->size1;
 171   size_t i;
 172   double a, b, c;
 173
 174   for (i = 0; i < N; ++i)
 175     {
 176       /* update lower triangle of A(i:n, i:n) */
 177
 178       a = gsl_matrix_get(A, i, i);
 179       b = gsl_matrix_get(B, i, i);
 180       a /= b * b;
 181       gsl_matrix_set(A, i, i, a);
 182
 183       if (i < N - 1)
 184         {
 185           gsl_vector_view ai = gsl_matrix_subcolumn(A, i, i + 1, N - i - 1);
 186           gsl_matrix_view ma =
 187             gsl_matrix_submatrix(A, i + 1, i + 1, N - i - 1, N - i - 1);
 188           gsl_vector_const_view bi =
 189             gsl_matrix_const_subcolumn(B, i, i + 1, N - i - 1);
 190           gsl_matrix_const_view mb =
 191             gsl_matrix_const_submatrix(B, i + 1, i + 1, N - i - 1, N - i - 1);
 192
 193           gsl_blas_dscal(1.0 / b, &ai.vector);
 194
 195           c = -0.5 * a;
 196           gsl_blas_daxpy(c, &bi.vector, &ai.vector);
 197
 198           gsl_blas_dsyr2(CblasLower, -1.0, &ai.vector, &bi.vector, &ma.matrix);
 199
 200           gsl_blas_daxpy(c, &bi.vector, &ai.vector);
 201
 202           gsl_blas_dtrsv(CblasLower,
 203                          CblasNoTrans,
 204                          CblasNonUnit,
 205                          &mb.matrix,
 206                          &ai.vector);
 207         }
 208     }
 209
 210   return GSL_SUCCESS;
 211 } /* gsl_eigen_gensymm_standardize() */