bug-gsl
[Top][All Lists]
Advanced

[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, 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, eps)
assert float_equal(2/(v_1*v_1), tau, eps)

```

reply via email to

 [Prev in Thread] Current Thread [Next in Thread]