|
Packit |
67cb25 |
/* linalg/sytd.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2001, 2007 Brian Gough
|
|
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 |
/* Factorise a symmetric matrix A into
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* A = Q T Q'
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* where Q is orthogonal and T is symmetric tridiagonal. Only the
|
|
Packit |
67cb25 |
* diagonal and lower triangular part of A is referenced and modified.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* On exit, T is stored in the diagonal and first subdiagonal of
|
|
Packit |
67cb25 |
* A. Since T is symmetric the upper diagonal is not stored.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Q is stored as a packed set of Householder transformations in the
|
|
Packit |
67cb25 |
* lower triangular part of the input matrix below the first subdiagonal.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* The full matrix for Q can be obtained as the product
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Q = Q_1 Q_2 ... Q_(N-2)
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* where
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Q_i = (I - tau_i * v_i * v_i')
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* and where v_i is a Householder vector
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* v_i = [0, ... , 0, 1, A(i+1,i), A(i+2,i), ... , A(N,i)]
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This storage scheme is the same as in LAPACK. See LAPACK's
|
|
Packit |
67cb25 |
* ssytd2.f for details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* See Golub & Van Loan, "Matrix Computations" (3rd ed), Section 8.3
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Note: this description uses 1-based indices. The code below uses
|
|
Packit |
67cb25 |
* 0-based indices
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_vector.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_matrix.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_blas.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <gsl/gsl_linalg.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_symmtd_decomp (gsl_matrix * A, gsl_vector * tau)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (A->size1 != A->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("symmetric tridiagonal decomposition requires square matrix",
|
|
Packit |
67cb25 |
GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (tau->size + 1 != A->size1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of tau must be (matrix size - 1)", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = A->size1;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0 ; i < N - 2; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_view c = gsl_matrix_column (A, i);
|
|
Packit |
67cb25 |
gsl_vector_view v = gsl_vector_subvector (&c.vector, i + 1, N - (i + 1));
|
|
Packit |
67cb25 |
double tau_i = gsl_linalg_householder_transform (&v.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Apply the transformation H^T A H to the remaining columns */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (tau_i != 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_view m = gsl_matrix_submatrix (A, i + 1, i + 1,
|
|
Packit |
67cb25 |
N - (i+1), N - (i+1));
|
|
Packit |
67cb25 |
double ei = gsl_vector_get(&v.vector, 0);
|
|
Packit |
67cb25 |
gsl_vector_view x = gsl_vector_subvector (tau, i, N-(i+1));
|
|
Packit |
67cb25 |
gsl_vector_set (&v.vector, 0, 1.0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* x = tau * A * v */
|
|
Packit |
67cb25 |
gsl_blas_dsymv (CblasLower, tau_i, &m.matrix, &v.vector, 0.0, &x.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* w = x - (1/2) tau * (x' * v) * v */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xv, alpha;
|
|
Packit |
67cb25 |
gsl_blas_ddot(&x.vector, &v.vector, &xv);
|
|
Packit |
67cb25 |
alpha = - (tau_i / 2.0) * xv;
|
|
Packit |
67cb25 |
gsl_blas_daxpy(alpha, &v.vector, &x.vector);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* apply the transformation A = A - v w' - w v' */
|
|
Packit |
67cb25 |
gsl_blas_dsyr2(CblasLower, -1.0, &v.vector, &x.vector, &m.matrix);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set (&v.vector, 0, ei);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set (tau, i, tau_i);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Form the orthogonal matrix Q from the packed QR matrix */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_symmtd_unpack (const gsl_matrix * A,
|
|
Packit |
67cb25 |
const gsl_vector * tau,
|
|
Packit |
67cb25 |
gsl_matrix * Q,
|
|
Packit |
67cb25 |
gsl_vector * diag,
|
|
Packit |
67cb25 |
gsl_vector * sdiag)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (A->size1 != A->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (tau->size + 1 != A->size1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of tau must be (matrix size - 1)", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (Q->size1 != A->size1 || Q->size2 != A->size1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of Q must match size of A", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (diag->size != A->size1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of diagonal must match size of A", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (sdiag->size + 1 != A->size1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of subdiagonal must be (matrix size - 1)", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = A->size1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Initialize Q to the identity */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_set_identity (Q);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = N - 2; i-- > 0;)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_const_view c = gsl_matrix_const_column (A, i);
|
|
Packit |
67cb25 |
gsl_vector_const_view h = gsl_vector_const_subvector (&c.vector, i + 1, N - (i+1));
|
|
Packit |
67cb25 |
double ti = gsl_vector_get (tau, i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_view m = gsl_matrix_submatrix (Q, i + 1, i + 1, N-(i+1), N-(i+1));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_linalg_householder_hm (ti, &h.vector, &m.matrix);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Copy diagonal into diag */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double Aii = gsl_matrix_get (A, i, i);
|
|
Packit |
67cb25 |
gsl_vector_set (diag, i, Aii);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Copy subdiagonal into sd */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N - 1; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double Aji = gsl_matrix_get (A, i+1, i);
|
|
Packit |
67cb25 |
gsl_vector_set (sdiag, i, Aji);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_symmtd_unpack_T (const gsl_matrix * A,
|
|
Packit |
67cb25 |
gsl_vector * diag,
|
|
Packit |
67cb25 |
gsl_vector * sdiag)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (A->size1 != A->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (diag->size != A->size1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of diagonal must match size of A", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (sdiag->size + 1 != A->size1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of subdiagonal must be (matrix size - 1)", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = A->size1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Copy diagonal into diag */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double Aii = gsl_matrix_get (A, i, i);
|
|
Packit |
67cb25 |
gsl_vector_set (diag, i, Aii);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Copy subdiagonal into sdiag */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N - 1; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double Aij = gsl_matrix_get (A, i+1, i);
|
|
Packit |
67cb25 |
gsl_vector_set (sdiag, i, Aij);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|