[Bug-gsl] linalg/symmtd

[Mario Pernici]

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)