gsl_linalg_complex_householder_transform

[Mario Pernici <Mario dot Pernici at mi dot infn dot it>]

Dear Sirs,
 it seems that the modification of 
gsl_linalg_complex_householder_transform
 works fine; here I give a similar correction in gsl_linalg_hermtd_decomp.

 As in the previous mail,
 householder matrix:
 P = 1 - beta u u^H
 P x = -e^{i theta} ||x||_2 * e_1

 so that in gsl_linalg_hermtd_decomp
 the phase e^{i theta} must be taken into account when setting
 gsl_vector_complex_set (&v.vector, 0, ei);
 see attached file.

-------------------
 In the manual there are misprints:
 Tridiagonal Decomposition of Hermitian Matrices
A hermitian matrix A can be factorized by similarity transformations
 into the form,
 A = U T U^T
(becomes A = U T U^H)
 where U is an unitary matrix and T is a real symmetric tridiagonal 
matrix.
 ( T is a hermitian tridiagonal matrix).
---------------------

 Below I give an example in ruby-gsl; comments are added, so it should
 be understandable without running it.

require 'assert'
require "GSL"; include GSL; include GSL::Math
require 'gsl_matrix'; require 'gsl_linalg'

# complex numbers 0 and 1
zero = Complex.new2(0,0)
one = Complex.new2(1,0)

# initialize a 4x4 hermitian matrix
a = Matrix_complex.new([2,0,  3,2, 4,1, 5,2],
                       [3,-2, 2,0, 8,1, 4,4],
                       [4,-1, 8,-1,3,0, 3,2],
                       [5,-2, 4,-4,3,-2, 1,0])

# check that it is hermitian
assert a == a.h
# get the output matrix a1 and the vector tau using hermtd_decomp
a1, tau = a.hermtd_decomp

# In the lower part of a1 are stored the elements of the tridiagonal
# matrix and the householder vectors
#
#       [t             ]
#  a1 = [t    t        ]
#       [v_0  t   t    ]
#       [v_0  v_1 t  t ]

# householder vectors
v_0 = Vector_complex.new2(zero, one, a1.get(2,0),a1.get(3,0))
v_1 = Vector_complex.new2(zero,zero,one,a1.get(3,1))

# householder matrices h_i = I - tau_i v_i v_i^H
# where tau_i = 2/(v_i^H v_i)
h_0 = v_0.house
h_1 = v_1.house
# Q
q = h_0 * h_1
# Q is unitary
assert q * q.h =~ Matrix_complex.identity(4)

# Get the tridiagonal matrix t1 out of a1
# first put to zero the elements in the lower part where are stored
# the householder vectors
a1[2,0] = zero; a1[3,0] = zero; a1[3,1] = zero

# get the full hermitian tridiagonal matrix from its lower part
t1 = a1.lower(true) + a1.lower(false).h

# t1 = Q^H * a * Q within precision 1.0e-14
assert t1 =~ [q.h * a * q, 1.0e-14]
p t1
#[ (2, 0)  (-6.3911, -4.26073)  (0, 0)  (0, 0)
#   (-6.3911, 4.26073)  (9.89831, 0)  (-5.12537, -5.45784)  (0, 0)
#   (0, 0)  (-5.12537, 5.45784)  (-0.50685, 0)  (1.10104, -2.20513)
#    (0, 0)  (0, 0)  (1.10104, 2.20513)  (-3.39146, 0)  ]

# check that tau is related to the norms of the householder vectors
n0 = v_0.dznrm2
n1 = v_1.dznrm2
assert float_equal(2/(n0**2),tau[0].real, 1.0e-15)
assert float_equal(2/(n1**2),tau[1].real, 1.0e-15)
p tau
#[(1.4694, 0)  (0.632016, 0)  (-1.54304, 0.69892)  ]

Attachment: hermtd
Description: Text document