3 * Copyright (C) 2007 Patrick Alken
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.
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.
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.
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>
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.
39 gsl_eigen_genherm_alloc()
41 Allocate a workspace for solving the generalized hermitian-definite
42 eigenvalue problem. The size of this workspace is O(3n).
44 Inputs: n - size of matrices
46 Return: pointer to workspace
49 gsl_eigen_genherm_workspace *
50 gsl_eigen_genherm_alloc(const size_t n)
52 gsl_eigen_genherm_workspace *w;
56 GSL_ERROR_NULL ("matrix dimension must be positive integer",
60 w = calloc (1, sizeof (gsl_eigen_genherm_workspace));
64 GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
69 w->herm_workspace_p = gsl_eigen_herm_alloc(n);
70 if (!w->herm_workspace_p)
72 gsl_eigen_genherm_free(w);
73 GSL_ERROR_NULL("failed to allocate space for herm workspace", GSL_ENOMEM);
77 } /* gsl_eigen_genherm_alloc() */
80 gsl_eigen_genherm_free()
85 gsl_eigen_genherm_free (gsl_eigen_genherm_workspace * w)
87 if (w->herm_workspace_p)
88 gsl_eigen_herm_free(w->herm_workspace_p);
91 } /* gsl_eigen_genherm_free() */
96 Solve the generalized hermitian-definite eigenvalue problem
100 for the eigenvalues \lambda.
102 Inputs: A - complex hermitian matrix
103 B - complex hermitian and positive definite matrix
104 eval - where to store eigenvalues
107 Return: success or error
111 gsl_eigen_genherm (gsl_matrix_complex * A, gsl_matrix_complex * B,
112 gsl_vector * eval, gsl_eigen_genherm_workspace * w)
114 const size_t N = A->size1;
116 /* check matrix and vector sizes */
120 GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
122 else if ((N != B->size1) || (N != B->size2))
124 GSL_ERROR ("B matrix dimensions must match A", GSL_EBADLEN);
126 else if (eval->size != N)
128 GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
130 else if (w->size != N)
132 GSL_ERROR ("matrix size does not match workspace", GSL_EBADLEN);
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 */
143 /* transform to standard hermitian eigenvalue problem */
144 gsl_eigen_genherm_standardize(A, B);
146 s = gsl_eigen_herm(A, eval, w->herm_workspace_p);
150 } /* gsl_eigen_genherm() */
153 gsl_eigen_genherm_standardize()
154 Reduce the generalized hermitian-definite eigenproblem to
155 the standard hermitian eigenproblem by computing
159 where L L^H is the Cholesky decomposition of B
161 Inputs: A - (input/output) complex hermitian matrix
162 B - complex hermitian, positive definite matrix in Cholesky form
166 Notes: A is overwritten by L^{-1} A L^{-H}
170 gsl_eigen_genherm_standardize(gsl_matrix_complex *A,
171 const gsl_matrix_complex *B)
173 const size_t N = A->size1;
178 GSL_SET_IMAG(&z, 0.0);
180 for (i = 0; i < N; ++i)
182 /* update lower triangle of A(i:n, i:n) */
184 y = gsl_matrix_complex_get(A, i, i);
186 y = gsl_matrix_complex_get(B, i, i);
190 gsl_matrix_complex_set(A, i, i, z);
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);
203 gsl_blas_zdscal(1.0 / b, &ai.vector);
205 GSL_SET_REAL(&z, -0.5 * a);
206 gsl_blas_zaxpy(z, &bi.vector, &ai.vector);
208 gsl_blas_zher2(CblasLower,
214 gsl_blas_zaxpy(z, &bi.vector, &ai.vector);
216 gsl_blas_ztrsv(CblasLower,
225 } /* gsl_eigen_genherm_standardize() */