This is the mail archive of the gsl-discuss@sources.redhat.com mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

gsl_linalg_complex_householder_transform


Hi Brian,
< I think I used the LAPACK method described in Lapack Working Note 72

 thank you for the reference!
 Now my examples for the functions in the section
 Tridiagonal Decomposition of Hermitian Matrices
 work without modifications of the code of gsl.
 The key point is that
 gsl_complex
 gsl_linalg_complex_householder_transform (gsl_vector_complex * v)
 returns a complex number tau such that
 (I - tau^* v v^H) x = -/+ ||x|| e_1
 (see Working Note 72)
 while I thought the relation was (I - tau v v^H) x = -/+ ||x|| e_1
 Maybe this should be added in the comments to the source code of
 this function.

 And you could add to householdercomplex.c the function
gsl_linalg_complex_householder_hv in analogy with householder.c .
======================
int
gsl_linalg_complex_householder_hv(gsl_complex tau, const 
gsl_vector_complex * v, gsl_vector_complex *  w)
{
  size_t i;
  gsl_complex z;

  if (GSL_REAL(tau) == 0.0 && GSL_IMAG(tau) == 0.0)
      return GSL_SUCCESS;

  /* z = v^H w */
  gsl_blas_zdotc(v, w, &z);

  /* w = w - tau * (v^H w) * v   */
  z = gsl_complex_mul(z, tau);
  for (i = 0; i < v->size; i++)
  {
     gsl_complex vi, wi;
         vi = gsl_vector_complex_get(v,i);
         wi = gsl_vector_complex_get(w,i);
         vi = gsl_complex_mul(vi, z);
         wi = gsl_complex_sub(wi, vi);
         gsl_vector_complex_set(w, i, wi);
   }

}

=============================
Bye
 Mario






Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]