|
Packit |
67cb25 |
/* linalg/bidiag.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 matrix A into
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* A = U B V^T
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* where U and V are orthogonal and B is upper bidiagonal.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* On exit, B is stored in the diagonal and first superdiagonal of A.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* U is stored as a packed set of Householder transformations in the
|
|
Packit |
67cb25 |
* lower triangular part of the input matrix below the diagonal.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* V is stored as a packed set of Householder transformations in the
|
|
Packit |
67cb25 |
* upper triangular part of the input matrix above the first
|
|
Packit |
67cb25 |
* superdiagonal.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* The full matrix for U can be obtained as the product
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* U = U_1 U_2 .. U_N
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* where
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* U_i = (I - tau_i * u_i * u_i')
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* and where u_i is a Householder vector
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* u_i = [0, .. , 0, 1, A(i+1,i), A(i+3,i), .. , A(M,i)]
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* The full matrix for V can be obtained as the product
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* V = V_1 V_2 .. V_(N-2)
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* where
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* V_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,i+2), A(i,i+3), .. , A(i,N)]
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* See Golub & Van Loan, "Matrix Computations" (3rd ed), Algorithm 5.4.2
|
|
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_bidiag_decomp (gsl_matrix * A, gsl_vector * tau_U, gsl_vector * tau_V)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (A->size1 < A->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("bidiagonal decomposition requires M>=N", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (tau_U->size != A->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of tau_U must be N", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (tau_V->size + 1 != A->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of tau_V must be (N - 1)", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t M = A->size1;
|
|
Packit |
67cb25 |
const size_t N = A->size2;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0 ; i < N; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Apply Householder transformation to current column */
|
|
Packit |
67cb25 |
|
|
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, M - i);
|
|
Packit |
67cb25 |
double tau_i = gsl_linalg_householder_transform (&v.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Apply the transformation to the remaining columns */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (i + 1 < N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_view m =
|
|
Packit |
67cb25 |
gsl_matrix_submatrix (A, i, i + 1, M - i, N - (i + 1));
|
|
Packit |
67cb25 |
gsl_linalg_householder_hm (tau_i, &v.vector, &m.matrix);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set (tau_U, i, tau_i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Apply Householder transformation to current row */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (i + 1 < N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_view r = gsl_matrix_row (A, i);
|
|
Packit |
67cb25 |
gsl_vector_view v = gsl_vector_subvector (&r.vector, i + 1, N - (i + 1));
|
|
Packit |
67cb25 |
double tau_i = gsl_linalg_householder_transform (&v.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Apply the transformation to the remaining rows */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (i + 1 < M)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_view m =
|
|
Packit |
67cb25 |
gsl_matrix_submatrix (A, i+1, i+1, M - (i+1), N - (i+1));
|
|
Packit |
67cb25 |
gsl_linalg_householder_mh (tau_i, &v.vector, &m.matrix);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set (tau_V, i, tau_i);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Form the orthogonal matrices U, V, diagonal d and superdiagonal sd
|
|
Packit |
67cb25 |
from the packed bidiagonal matrix A */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_bidiag_unpack (const gsl_matrix * A,
|
|
Packit |
67cb25 |
const gsl_vector * tau_U,
|
|
Packit |
67cb25 |
gsl_matrix * U,
|
|
Packit |
67cb25 |
const gsl_vector * tau_V,
|
|
Packit |
67cb25 |
gsl_matrix * V,
|
|
Packit |
67cb25 |
gsl_vector * diag,
|
|
Packit |
67cb25 |
gsl_vector * superdiag)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t M = A->size1;
|
|
Packit |
67cb25 |
const size_t N = A->size2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const size_t K = GSL_MIN(M, N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (M < N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix A must have M >= N", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (tau_U->size != K)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (tau_V->size + 1 != K)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of tau must be MIN(M,N) - 1", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (U->size1 != M || U->size2 != N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of U must be M x N", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (V->size1 != N || V->size2 != N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of V must be N x N", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (diag->size != K)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of diagonal must match size of A", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (superdiag->size + 1 != K)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of subdiagonal must be (diagonal size - 1)", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i, j;
|
|
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 superdiagonal into superdiag */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N - 1; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double Aij = gsl_matrix_get (A, i, i+1);
|
|
Packit |
67cb25 |
gsl_vector_set (superdiag, i, Aij);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Initialize V to the identity */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_set_identity (V);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = N - 1; i-- > 0;)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Householder row transformation to accumulate V */
|
|
Packit |
67cb25 |
gsl_vector_const_view r = gsl_matrix_const_row (A, i);
|
|
Packit |
67cb25 |
gsl_vector_const_view h =
|
|
Packit |
67cb25 |
gsl_vector_const_subvector (&r.vector, i + 1, N - (i+1));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double ti = gsl_vector_get (tau_V, i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_view m =
|
|
Packit |
67cb25 |
gsl_matrix_submatrix (V, 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 |
/* Initialize U to the identity */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_set_identity (U);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = N; j-- > 0;)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Householder column transformation to accumulate U */
|
|
Packit |
67cb25 |
gsl_vector_const_view c = gsl_matrix_const_column (A, j);
|
|
Packit |
67cb25 |
gsl_vector_const_view h = gsl_vector_const_subvector (&c.vector, j, M - j);
|
|
Packit |
67cb25 |
double tj = gsl_vector_get (tau_U, j);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_view m =
|
|
Packit |
67cb25 |
gsl_matrix_submatrix (U, j, j, M-j, N-j);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_linalg_householder_hm (tj, &h.vector, &m.matrix);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_bidiag_unpack2 (gsl_matrix * A,
|
|
Packit |
67cb25 |
gsl_vector * tau_U,
|
|
Packit |
67cb25 |
gsl_vector * tau_V,
|
|
Packit |
67cb25 |
gsl_matrix * V)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t M = A->size1;
|
|
Packit |
67cb25 |
const size_t N = A->size2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const size_t K = GSL_MIN(M, N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (M < N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix A must have M >= N", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (tau_U->size != K)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (tau_V->size + 1 != K)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of tau must be MIN(M,N) - 1", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (V->size1 != N || V->size2 != N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of V must be N x N", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i, j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Initialize V to the identity */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_set_identity (V);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = N - 1; i-- > 0;)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Householder row transformation to accumulate V */
|
|
Packit |
67cb25 |
gsl_vector_const_view r = gsl_matrix_const_row (A, i);
|
|
Packit |
67cb25 |
gsl_vector_const_view h =
|
|
Packit |
67cb25 |
gsl_vector_const_subvector (&r.vector, i + 1, N - (i+1));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double ti = gsl_vector_get (tau_V, i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_view m =
|
|
Packit |
67cb25 |
gsl_matrix_submatrix (V, 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 superdiagonal into tau_v */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N - 1; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double Aij = gsl_matrix_get (A, i, i+1);
|
|
Packit |
67cb25 |
gsl_vector_set (tau_V, i, Aij);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Allow U to be unpacked into the same memory as A, copy
|
|
Packit |
67cb25 |
diagonal into tau_U */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = N; j-- > 0;)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Householder column transformation to accumulate U */
|
|
Packit |
67cb25 |
double tj = gsl_vector_get (tau_U, j);
|
|
Packit |
67cb25 |
double Ajj = gsl_matrix_get (A, j, j);
|
|
Packit |
67cb25 |
gsl_matrix_view m = gsl_matrix_submatrix (A, j, j, M-j, N-j);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set (tau_U, j, Ajj);
|
|
Packit |
67cb25 |
gsl_linalg_householder_hm1 (tj, &m.matrix);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_bidiag_unpack_B (const gsl_matrix * A,
|
|
Packit |
67cb25 |
gsl_vector * diag,
|
|
Packit |
67cb25 |
gsl_vector * superdiag)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t M = A->size1;
|
|
Packit |
67cb25 |
const size_t N = A->size2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const size_t K = GSL_MIN(M, N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (diag->size != K)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of diagonal must match size of A", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (superdiag->size + 1 != K)
|
|
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 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Copy diagonal into diag */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < K; 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 superdiagonal into superdiag */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < K - 1; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double Aij = gsl_matrix_get (A, i, i+1);
|
|
Packit |
67cb25 |
gsl_vector_set (superdiag, i, Aij);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|