|
Packit |
67cb25 |
/* linalg/hessenberg.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2006 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 <config.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_linalg.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_matrix.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_vector.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_linalg_hessenberg_decomp()
|
|
Packit |
67cb25 |
Compute the Householder reduction to Hessenberg form of a
|
|
Packit |
67cb25 |
square N-by-N matrix A.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
H = U^t A U
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
See Golub & Van Loan, "Matrix Computations" (3rd ed), algorithm
|
|
Packit |
67cb25 |
7.4.2
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: A - matrix to reduce
|
|
Packit |
67cb25 |
tau - where to store scalar factors in Householder
|
|
Packit |
67cb25 |
matrices; this vector must be of length N,
|
|
Packit |
67cb25 |
where N is the order of A
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: GSL_SUCCESS unless error occurs
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Notes: on output, the upper triangular portion of A (including
|
|
Packit |
67cb25 |
the diagaonal and subdiagonal) contains the Hessenberg matrix.
|
|
Packit |
67cb25 |
The lower triangular portion (below the subdiagonal) contains
|
|
Packit |
67cb25 |
the Householder vectors which can be used to construct
|
|
Packit |
67cb25 |
the similarity transform matrix U.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The matrix U is
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
U = U(1) U(2) ... U(n - 2)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
U(i) = I - tau(i) * v(i) * v(i)^t
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
and the vector v(i) is stored in column i of the matrix A
|
|
Packit |
67cb25 |
underneath the subdiagonal. So the first element of v(i)
|
|
Packit |
67cb25 |
is stored in row i + 2, column i, the second element at
|
|
Packit |
67cb25 |
row i + 3, column i, and so on.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Also note that for the purposes of computing U(i),
|
|
Packit |
67cb25 |
v(1:i) = 0, v(i + 1) = 1, and v(i+2:n) is what is stored in
|
|
Packit |
67cb25 |
column i of A beneath the subdiagonal.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_hessenberg_decomp(gsl_matrix *A, gsl_vector *tau)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = A->size1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (N != A->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("Hessenberg reduction requires square matrix",
|
|
Packit |
67cb25 |
GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (N != tau->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("tau vector must match matrix size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (N < 3)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* nothing to do */
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i; /* looping */
|
|
Packit |
67cb25 |
gsl_vector_view c, /* matrix column */
|
|
Packit |
67cb25 |
hv; /* householder vector */
|
|
Packit |
67cb25 |
gsl_matrix_view m;
|
|
Packit |
67cb25 |
double tau_i; /* beta in algorithm 7.4.2 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N - 2; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* make a copy of A(i + 1:n, i) and store it in the section
|
|
Packit |
67cb25 |
* of 'tau' that we haven't stored coefficients in yet
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
c = gsl_matrix_subcolumn(A, i, i + 1, N - i - 1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
hv = gsl_vector_subvector(tau, i + 1, N - (i + 1));
|
|
Packit |
67cb25 |
gsl_vector_memcpy(&hv.vector, &c.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute householder transformation of A(i+1:n,i) */
|
|
Packit |
67cb25 |
tau_i = gsl_linalg_householder_transform(&hv.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* apply left householder matrix (I - tau_i v v') to A */
|
|
Packit |
67cb25 |
m = gsl_matrix_submatrix(A, i + 1, i, N - (i + 1), N - i);
|
|
Packit |
67cb25 |
gsl_linalg_householder_hm(tau_i, &hv.vector, &m.matrix);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* apply right householder matrix (I - tau_i v v') to A */
|
|
Packit |
67cb25 |
m = gsl_matrix_submatrix(A, 0, i + 1, N, N - (i + 1));
|
|
Packit |
67cb25 |
gsl_linalg_householder_mh(tau_i, &hv.vector, &m.matrix);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* save Householder coefficient */
|
|
Packit |
67cb25 |
gsl_vector_set(tau, i, tau_i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* store Householder vector below the subdiagonal in column
|
|
Packit |
67cb25 |
* i of the matrix. hv(1) does not need to be stored since
|
|
Packit |
67cb25 |
* it is always 1.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
c = gsl_vector_subvector(&c.vector, 1, c.vector.size - 1);
|
|
Packit |
67cb25 |
hv = gsl_vector_subvector(&hv.vector, 1, hv.vector.size - 1);
|
|
Packit |
67cb25 |
gsl_vector_memcpy(&c.vector, &hv.vector);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* gsl_linalg_hessenberg_decomp() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_linalg_hessenberg_unpack()
|
|
Packit |
67cb25 |
Construct the matrix U which transforms a matrix A into
|
|
Packit |
67cb25 |
its upper Hessenberg form:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
H = U^t A U
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
by unpacking the information stored in H from gsl_linalg_hessenberg().
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
U is a product of Householder matrices:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
U = U(1) U(2) ... U(n - 2)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
U(i) = I - tau(i) * v(i) * v(i)^t
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The v(i) are stored in the lower triangular part of H by
|
|
Packit |
67cb25 |
gsl_linalg_hessenberg(). The tau(i) are stored in the vector tau.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: H - Hessenberg matrix computed from
|
|
Packit |
67cb25 |
gsl_linalg_hessenberg()
|
|
Packit |
67cb25 |
tau - tau vector computed from gsl_linalg_hessenberg()
|
|
Packit |
67cb25 |
U - (output) where to store similarity matrix
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: success or error
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_hessenberg_unpack(gsl_matrix * H, gsl_vector * tau,
|
|
Packit |
67cb25 |
gsl_matrix * U)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_set_identity(U);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s = gsl_linalg_hessenberg_unpack_accum(H, tau, U);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
} /* gsl_linalg_hessenberg_unpack() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_linalg_hessenberg_unpack_accum()
|
|
Packit |
67cb25 |
This routine is the same as gsl_linalg_hessenberg_unpack(), except
|
|
Packit |
67cb25 |
instead of storing the similarity matrix in U, it accumulates it,
|
|
Packit |
67cb25 |
so that
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
U -> U * [ U(1) U(2) ... U(n - 2) ]
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
instead of:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
U -> U(1) U(2) ... U(n - 2)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: H - Hessenberg matrix computed from
|
|
Packit |
67cb25 |
gsl_linalg_hessenberg()
|
|
Packit |
67cb25 |
tau - tau vector computed from gsl_linalg_hessenberg()
|
|
Packit |
67cb25 |
V - (input/output) where to accumulate similarity matrix
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: success or error
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Notes: 1) On input, V needs to be initialized. The Householder matrices
|
|
Packit |
67cb25 |
are accumulated into V, so on output,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
V_out = V_in * U(1) * U(2) * ... * U(n - 2)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
so if you just want the product of the Householder matrices,
|
|
Packit |
67cb25 |
initialize V to the identity matrix before calling this
|
|
Packit |
67cb25 |
function.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
2) V does not have to be square, but must have the same
|
|
Packit |
67cb25 |
number of columns as the order of H
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_hessenberg_unpack_accum(gsl_matrix * H, gsl_vector * tau,
|
|
Packit |
67cb25 |
gsl_matrix * V)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = H->size1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (N != H->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("Hessenberg reduction requires square matrix",
|
|
Packit |
67cb25 |
GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (N != tau->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("tau vector must match matrix size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (N != V->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("V matrix has wrong dimension", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t j; /* looping */
|
|
Packit |
67cb25 |
double tau_j; /* householder coefficient */
|
|
Packit |
67cb25 |
gsl_vector_view c, /* matrix column */
|
|
Packit |
67cb25 |
hv; /* householder vector */
|
|
Packit |
67cb25 |
gsl_matrix_view m;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (N < 3)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* nothing to do */
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < (N - 2); ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
c = gsl_matrix_column(H, j);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
tau_j = gsl_vector_get(tau, j);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* get a view to the householder vector in column j, but
|
|
Packit |
67cb25 |
* make sure hv(2) starts at the element below the
|
|
Packit |
67cb25 |
* subdiagonal, since hv(1) was never stored and is always
|
|
Packit |
67cb25 |
* 1
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
hv = gsl_vector_subvector(&c.vector, j + 1, N - (j + 1));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* Only operate on part of the matrix since the first
|
|
Packit |
67cb25 |
* j + 1 entries of the real householder vector are 0
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* V -> V * U(j)
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Note here that V->size1 is not necessarily equal to N
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
m = gsl_matrix_submatrix(V, 0, j + 1, V->size1, N - (j + 1));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* apply right Householder matrix to V */
|
|
Packit |
67cb25 |
gsl_linalg_householder_mh(tau_j, &hv.vector, &m.matrix);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* gsl_linalg_hessenberg_unpack_accum() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_linalg_hessenberg_set_zero()
|
|
Packit |
67cb25 |
Zero out the lower triangular portion of the Hessenberg matrix H.
|
|
Packit |
67cb25 |
This is useful when Householder vectors may be stored in the lower
|
|
Packit |
67cb25 |
part of H, but eigenvalue solvers need some scratch space with zeros.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_hessenberg_set_zero(gsl_matrix * H)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = H->size1;
|
|
Packit |
67cb25 |
size_t i, j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (N < 3)
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < N - 2; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (i = j + 2; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_set(H, i, j, 0.0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
} /* gsl_linalg_hessenberg_set_zero() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_linalg_hessenberg_submatrix()
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This routine does the same thing as gsl_linalg_hessenberg(),
|
|
Packit |
67cb25 |
except that it operates on a submatrix of a larger matrix, but
|
|
Packit |
67cb25 |
updates the larger matrix with the Householder transformations.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
For example, suppose
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
M = [ M_{11} | M_{12} | M_{13} ]
|
|
Packit |
67cb25 |
[ 0 | A | M_{23} ]
|
|
Packit |
67cb25 |
[ 0 | 0 | M_{33} ]
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where M_{11} and M_{33} are already in Hessenberg form, and we
|
|
Packit |
67cb25 |
just want to reduce A to Hessenberg form. Applying the transformations
|
|
Packit |
67cb25 |
to A alone will cause the larger matrix M to lose its similarity
|
|
Packit |
67cb25 |
information. So this routine updates M_{12} and M_{23} as A gets
|
|
Packit |
67cb25 |
reduced.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: M - total matrix
|
|
Packit |
67cb25 |
A - (sub)matrix to reduce
|
|
Packit |
67cb25 |
top - row index of top of A in M
|
|
Packit |
67cb25 |
tau - where to store scalar factors in Householder
|
|
Packit |
67cb25 |
matrices; this vector must be of length N,
|
|
Packit |
67cb25 |
where N is the order of A
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: GSL_SUCCESS unless error occurs
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Notes: on output, the upper triangular portion of A (including
|
|
Packit |
67cb25 |
the diagaonal and subdiagonal) contains the Hessenberg matrix.
|
|
Packit |
67cb25 |
The lower triangular portion (below the subdiagonal) contains
|
|
Packit |
67cb25 |
the Householder vectors which can be used to construct
|
|
Packit |
67cb25 |
the similarity transform matrix U.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The matrix U is
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
U = U(1) U(2) ... U(n - 2)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
U(i) = I - tau(i) * v(i) * v(i)^t
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
and the vector v(i) is stored in column i of the matrix A
|
|
Packit |
67cb25 |
underneath the subdiagonal. So the first element of v(i)
|
|
Packit |
67cb25 |
is stored in row i + 2, column i, the second element at
|
|
Packit |
67cb25 |
row i + 3, column i, and so on.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Also note that for the purposes of computing U(i),
|
|
Packit |
67cb25 |
v(1:i) = 0, v(i + 1) = 1, and v(i+2:n) is what is stored in
|
|
Packit |
67cb25 |
column i of A beneath the subdiagonal.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_hessenberg_submatrix(gsl_matrix *M, gsl_matrix *A, size_t top,
|
|
Packit |
67cb25 |
gsl_vector *tau)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = A->size1;
|
|
Packit |
67cb25 |
const size_t N_M = M->size1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (N != A->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("Hessenberg reduction requires square matrix",
|
|
Packit |
67cb25 |
GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (N != tau->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("tau vector must match matrix size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (N < 3)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* nothing to do */
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i; /* looping */
|
|
Packit |
67cb25 |
gsl_vector_view c, /* matrix column */
|
|
Packit |
67cb25 |
hv; /* householder vector */
|
|
Packit |
67cb25 |
gsl_matrix_view m;
|
|
Packit |
67cb25 |
double tau_i; /* beta in algorithm 7.4.2 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N - 2; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* make a copy of A(i + 1:n, i) and store it in the section
|
|
Packit |
67cb25 |
* of 'tau' that we haven't stored coefficients in yet
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
c = gsl_matrix_subcolumn(A, i, i + 1, N - i - 1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
hv = gsl_vector_subvector(tau, i + 1, N - (i + 1));
|
|
Packit |
67cb25 |
gsl_vector_memcpy(&hv.vector, &c.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute householder transformation of A(i+1:n,i) */
|
|
Packit |
67cb25 |
tau_i = gsl_linalg_householder_transform(&hv.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* apply left householder matrix (I - tau_i v v') to
|
|
Packit |
67cb25 |
* [ A | M_{23} ]
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
m = gsl_matrix_submatrix(M,
|
|
Packit |
67cb25 |
top + i + 1,
|
|
Packit |
67cb25 |
top + i,
|
|
Packit |
67cb25 |
N - (i + 1),
|
|
Packit |
67cb25 |
N_M - top - i);
|
|
Packit |
67cb25 |
gsl_linalg_householder_hm(tau_i, &hv.vector, &m.matrix);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* apply right householder matrix (I - tau_i v v') to
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* [ M_{12} ]
|
|
Packit |
67cb25 |
* [ A ]
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
m = gsl_matrix_submatrix(M,
|
|
Packit |
67cb25 |
0,
|
|
Packit |
67cb25 |
top + i + 1,
|
|
Packit |
67cb25 |
top + N,
|
|
Packit |
67cb25 |
N - (i + 1));
|
|
Packit |
67cb25 |
gsl_linalg_householder_mh(tau_i, &hv.vector, &m.matrix);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* save Householder coefficient */
|
|
Packit |
67cb25 |
gsl_vector_set(tau, i, tau_i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* store Householder vector below the subdiagonal in column
|
|
Packit |
67cb25 |
* i of the matrix. hv(1) does not need to be stored since
|
|
Packit |
67cb25 |
* it is always 1.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
c = gsl_vector_subvector(&c.vector, 1, c.vector.size - 1);
|
|
Packit |
67cb25 |
hv = gsl_vector_subvector(&hv.vector, 1, hv.vector.size - 1);
|
|
Packit |
67cb25 |
gsl_vector_memcpy(&c.vector, &hv.vector);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* gsl_linalg_hessenberg_submatrix() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#ifndef GSL_DISABLE_DEPRECATED
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* To support gsl-1.9 interface: DEPRECATED */
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_hessenberg(gsl_matrix *A, gsl_vector *tau)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return gsl_linalg_hessenberg_decomp(A, tau);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#endif
|