|
Packit |
67cb25 |
/* eigen/gensymm.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2007 Patrick Alken
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_eigen.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_linalg.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_blas.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_vector.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_matrix.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* This module computes the eigenvalues of a real generalized
|
|
Packit |
67cb25 |
* symmetric-definite eigensystem A x = \lambda B x, where A and
|
|
Packit |
67cb25 |
* B are symmetric, and B is positive-definite.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_eigen_gensymm_alloc()
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Allocate a workspace for solving the generalized symmetric-definite
|
|
Packit |
67cb25 |
eigenvalue problem. The size of this workspace is O(2n).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: n - size of matrices
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: pointer to workspace
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_gensymm_workspace *
|
|
Packit |
67cb25 |
gsl_eigen_gensymm_alloc(const size_t n)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_eigen_gensymm_workspace *w;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (n == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_NULL ("matrix dimension must be positive integer",
|
|
Packit |
67cb25 |
GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
w = (gsl_eigen_gensymm_workspace *) calloc (1, sizeof (gsl_eigen_gensymm_workspace));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (w == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
w->size = n;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
w->symm_workspace_p = gsl_eigen_symm_alloc(n);
|
|
Packit |
67cb25 |
if (!w->symm_workspace_p)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_eigen_gensymm_free(w);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL("failed to allocate space for symm workspace", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return (w);
|
|
Packit |
67cb25 |
} /* gsl_eigen_gensymm_alloc() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_eigen_gensymm_free()
|
|
Packit |
67cb25 |
Free workspace w
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
gsl_eigen_gensymm_free (gsl_eigen_gensymm_workspace * w)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
RETURN_IF_NULL (w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (w->symm_workspace_p)
|
|
Packit |
67cb25 |
gsl_eigen_symm_free(w->symm_workspace_p);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
free(w);
|
|
Packit |
67cb25 |
} /* gsl_eigen_gensymm_free() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_eigen_gensymm()
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Solve the generalized symmetric-definite eigenvalue problem
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A x = \lambda B x
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for the eigenvalues \lambda.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: A - real symmetric matrix
|
|
Packit |
67cb25 |
B - real symmetric and positive definite matrix
|
|
Packit |
67cb25 |
eval - where to store eigenvalues
|
|
Packit |
67cb25 |
w - workspace
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: success or error
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_eigen_gensymm (gsl_matrix * A, gsl_matrix * B, gsl_vector * eval,
|
|
Packit |
67cb25 |
gsl_eigen_gensymm_workspace * w)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = A->size1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* check matrix and vector sizes */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (N != A->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if ((N != B->size1) || (N != B->size2))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("B matrix dimensions must match A", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (eval->size != N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (w->size != N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size does not match workspace", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute Cholesky factorization of B */
|
|
Packit |
67cb25 |
s = gsl_linalg_cholesky_decomp1(B);
|
|
Packit |
67cb25 |
if (s != GSL_SUCCESS)
|
|
Packit |
67cb25 |
return s; /* B is not positive definite */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* transform to standard symmetric eigenvalue problem */
|
|
Packit |
67cb25 |
gsl_eigen_gensymm_standardize(A, B);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s = gsl_eigen_symm(A, eval, w->symm_workspace_p);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* gsl_eigen_gensymm() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_eigen_gensymm_standardize()
|
|
Packit |
67cb25 |
Reduce the generalized symmetric-definite eigenproblem to
|
|
Packit |
67cb25 |
the standard symmetric eigenproblem by computing
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
C = L^{-1} A L^{-t}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where L L^t is the Cholesky decomposition of B
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: A - (input/output) real symmetric matrix
|
|
Packit |
67cb25 |
B - real symmetric, positive definite matrix in Cholesky form
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: success
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Notes: A is overwritten by L^{-1} A L^{-t}
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_eigen_gensymm_standardize(gsl_matrix *A, const gsl_matrix *B)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = A->size1;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
double a, b, c;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* update lower triangle of A(i:n, i:n) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
a = gsl_matrix_get(A, i, i);
|
|
Packit |
67cb25 |
b = gsl_matrix_get(B, i, i);
|
|
Packit |
67cb25 |
a /= b * b;
|
|
Packit |
67cb25 |
gsl_matrix_set(A, i, i, a);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (i < N - 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_view ai = gsl_matrix_subcolumn(A, i, i + 1, N - i - 1);
|
|
Packit |
67cb25 |
gsl_matrix_view ma =
|
|
Packit |
67cb25 |
gsl_matrix_submatrix(A, i + 1, i + 1, N - i - 1, N - i - 1);
|
|
Packit |
67cb25 |
gsl_vector_const_view bi =
|
|
Packit |
67cb25 |
gsl_matrix_const_subcolumn(B, i, i + 1, N - i - 1);
|
|
Packit |
67cb25 |
gsl_matrix_const_view mb =
|
|
Packit |
67cb25 |
gsl_matrix_const_submatrix(B, i + 1, i + 1, N - i - 1, N - i - 1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_dscal(1.0 / b, &ai.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
c = -0.5 * a;
|
|
Packit |
67cb25 |
gsl_blas_daxpy(c, &bi.vector, &ai.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_dsyr2(CblasLower, -1.0, &ai.vector, &bi.vector, &ma.matrix);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_daxpy(c, &bi.vector, &ai.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_dtrsv(CblasLower,
|
|
Packit |
67cb25 |
CblasNoTrans,
|
|
Packit |
67cb25 |
CblasNonUnit,
|
|
Packit |
67cb25 |
&mb.matrix,
|
|
Packit |
67cb25 |
&ai.vector);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
} /* gsl_eigen_gensymm_standardize() */
|