gsl_linalg_complex_householder_transform

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

Dear Sirs,
 it seems that  gsl_linalg_complex_householder_transform
 does not follow Golub and Van Loan; is it done
 on purpose or there is a bug?

 Below I present the problem;
 in an attached files I give a modification of the
 code for gsl_linalg_complex_householder_transform
 according to Golub and Van Loan, the function
 gsl_linalg_complex_householder_hv
 and an example. Using the modified code for
 gsl_linalg_complex_householder_transform the example gives the
 expected results (es.1.out); the original code gives results
 not according to Golub and Van Loan (es.2.out).
 Best regards
   Mario

Problem:
 In the function
 double gsl_linalg_householder_transform(gsl_vector * v)
 there is a rescaling
  gsl_blas_dscal (1.0 / (alpha - beta), &x.vector);
 such that the first component of the householder vector is rescaled
 from (alpha - beta) to 1 (and then it is replaced by the norm).

 However this does not happen in
 gsl_complex
 gsl_linalg_complex_householder_transform (gsl_vector_complex * v).
 According to Golub and Van Loan, Problem P5.1.3, (their notation)
 householder matrix:
 P = 1 - beta u u^H
 P x = -e^{i theta} ||x||_2 * e_1
 householder vector:
 u = (x + e^{i theta} ||x||_2 * e_1 * s
 where x(1) = |x(1)| * e^{i theta}
 and where I added the scaling factor s, chosen to be
 s = 1/(x(1) + e^{i theta} ||x||_2) to set first component of the
 householder vector u(1) = 1.
 In the notation of linalg/householdercomplex.c the rescaling factor is
 s = 1/(alpha - e^{i theta} * beta_r);
 beta = (||x||_2 + |x(1)|)/||x||_2  becomes, in gsl notation,
 tau = (beta_r - |alpha|)/beta_r; notice that it is real (but I used
 a gsl_complex for it, not to change notation).

 In the gsl code the scaling factor is
 s = 1/(alpha - beta_r)
 so that u(1) is not rescaled to 1 when alpha = x(1) is not real.

Attachment: house.tar.gz
Description: GNU Zip compressed data