gsl_linalg_complex_householder_transform

[Mario Pernici]

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)  ]



-------------- next part --------------
/* in linalg/hermtd.c */

int 
gsl_linalg_hermtd_decomp (gsl_matrix_complex * A, gsl_vector_complex * tau)  
{
  if (A->size1 != A->size2)
    {
      GSL_ERROR ("hermitian tridiagonal decomposition requires square matrix",
                 GSL_ENOTSQR);
    }
  else if (tau->size + 1 != A->size1)
    {
      GSL_ERROR ("size of tau must be (matrix size - 1)", GSL_EBADLEN);
    }
  else
    {
      const size_t N = A->size1;
      size_t i;
  
      const gsl_complex zero = gsl_complex_rect (0.0, 0.0);
      const gsl_complex one = gsl_complex_rect (1.0, 0.0);
      const gsl_complex neg_one = gsl_complex_rect (-1.0, 0.0);

      for (i = 0 ; i < N - 2; i++)
        {
          gsl_vector_complex_view c = gsl_matrix_complex_column (A, i);
          gsl_vector_complex_view v = gsl_vector_complex_subvector (&c.vector, i + 1, N - (i + 1));
          gsl_complex c0 = gsl_vector_complex_get(&c.vector, i+1);
          gsl_complex tau_i = gsl_linalg_complex_householder_transform (&v.vector);
	  /* get the phase of the first element */
	  c0 = gsl_complex_div_real(c0, gsl_complex_abs(c0));
          /* Apply the transformation H^T A H to the remaining columns */

          if ( !(GSL_REAL(tau_i) == 0.0 && GSL_IMAG(tau_i) == 0.0)) 
            {
              gsl_matrix_complex_view m = 
                gsl_matrix_complex_submatrix (A, i + 1, i + 1, 
                                              N - (i+1), N - (i+1));
              gsl_complex ei = gsl_vector_complex_get(&v.vector, 0);
              gsl_vector_complex_view x = gsl_vector_complex_subvector (tau, i, N-(i+1));
              gsl_vector_complex_set (&v.vector, 0, one);
              
              /* x = tau * A * v */
              gsl_blas_zhemv (CblasLower, tau_i, &m.matrix, &v.vector, zero, &x.vector);

              /* w = x - (1/2) tau * (x' * v) * v  */
              {
                gsl_complex xv, txv, alpha;
                gsl_blas_zdotc(&x.vector, &v.vector, &xv);
                txv = gsl_complex_mul(tau_i, xv);
                alpha = gsl_complex_mul_real(txv, -0.5);
                gsl_blas_zaxpy(alpha, &v.vector, &x.vector);
              }
              
              /* apply the transformation A = A - v w' - w v' */
              gsl_blas_zher2(CblasLower, neg_one, &v.vector, &x.vector, &m.matrix);

              /* multiply ei by the phase */
	      ei = gsl_complex_mul(ei, c0);
              gsl_vector_complex_set (&v.vector, 0, ei);
            }
          
          gsl_vector_complex_set (tau, i, tau_i);
        }
      
      return GSL_SUCCESS;
    }
}