[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] linalg/symmtd
From: |
Mario Pernici |
Subject: |
[Bug-gsl] linalg/symmtd |
Date: |
Sun, 11 May 2003 22:24:29 +0200 (CEST) |
gsl version: 1.3
Red Hat 7.3
gcc version 2.96
Dear Sirs,
There are errors in the comments in linalg/symmtd.c
See attached the proposed correction, a comment and an example.
Best Regards
Mario Pernici
===================================
Proposed correction:
/* Factorise a symmetric matrix A into
*
* A = Q T Q'
*
*.......................
* The full matrix for Q can be obtained as the product
*
* Q = Q_0 Q_1..Q_(N-3)
* where
*
* Q_i = (I - tau_i * v_i * v_i')
*
* and where v_i is a N-dimensional Householder vector
*
* v_i = [0,..,0, 1, A(i+2,i), A(i+3,i), ... , A(N-1,i)]
*
* This storage scheme is the same as in LAPACK. See LAPACK's
* ssytd2.f for details.
*
* See Golub & Van Loan, "Matrix Computations" (3rd ed), Section 8.3 */
================================
Comment:
If I understand correctly the description in dsytd2.f,
there is a difference with respect to Lapack,
where (rewritten zero-based Q = H(1) H(2) . . . H(n-1) in dsytd2.f)
Q = Q_0 Q_1..Q_(N-2)
so I suppose that
Q_(N-2) = diag(1,1,..,1,-1)
is the householder matrix corresponding to the householder vector
v_(N-2) = [0,...,0,1]
(the last non-trivial householder vector is v_(N-3) =
[0,..,0,1,A(N-1,N-3)])
Both formulations are correct, with
T_LAPACK = Q_(N-2) T_GSL Q_(N-2)'
that is, the tridiagonal matrices differ in sign in the element T(N-2,N-1)
.
Here is an example; I hope you don't mind it is in ruby-gsl!
I added comments so you do not need to run it.
==============================================
require 'assert'
require "GSL"; include GSL; include GSL::Math
require 'gsl_matrix'; require 'gsl_linalg'
eps = 1.0e-15
# Example with N x N matrix a, with N = 3
puts "symmtd_decomp"
a = Matrix.new([1,3,4],
[3,2,8],
[4,8,3])
# a1 contains the diag, subdiag and householder vectors
a1, tau = a.symmtd_decomp
# Householder vector v_0
v_0 = Vector.new([0,1,a1.get(2,0)])
# Q_0
h_0 = v_0.house
# check that tau contains 2/(v_i^T v_i) within precision 1.0e-15
assert float_equal(2/(v_0*v_0), tau[0], eps)
# Q = Q_0 Q_1..Q_(N-3)
q = h_0
# Get the tridiagonal matrix T from a1
a1.set(2,0,0)
t1 = a1.lower(true) + a1.lower(false).t
# check that T = Q^T A Q within precision 1.0e-14
assert Matrix::float_equal(t1, q.t * a * q, 1.0e-14)
# Example with N x N matrix a, with N = 4
a = Matrix.new([1,3,4,5],
[3,2,8,9],
[4,8,3,7],
[5,9,7,6])
# a1 contains the diag, subdiag and householder vectors
a1, tau = a.symmtd_decomp
# Householder vectors
v_0 = Vector.new([0,1,a1.get(2,0), a1.get(3,0)])
v_1 = Vector.new([0,0,1, a1.get(3,1)])
# Q_0 and Q_1
h_0 = v_0.house
h_1 = v_1.house
# Q = Q_0 Q_1..Q_(N-3)
q = h_0 * h_1
# Get the tridiagonal matrix
a1.set(2,0,0)
a1.set(3,0,0); a1.set(3,1,0)
t1 = a1.lower(true) + a1.lower(false).t
# check that T = Q^T A Q within precision 1.0e-14
assert Matrix::float_equal(t1, q.t * a * q, 1.0e-14)
# check that tau contains 2/(v_i^T v_i) within precision 1.0e-15
assert float_equal(2/(v_0*v_0), tau[0], eps)
assert float_equal(2/(v_1*v_1), tau[1], eps)
- [Bug-gsl] linalg/symmtd,
Mario Pernici <=