Complex Householder transform functions

[Z F <mail4me9999 at yahoo dot com>]

Hello

Just discovered that some functions from householder transform set 
are missing. These functions may be needed for complex SVD.

Anyway, attached is a patch which should be applied to linalg
directory. It will add the

gsl_linalg_complex_householder_mh()
gsl_linalg_complex_householder_hm1()

Also, all complex householder functions should be added to the manual.

There are also some minor typos in sections "Vectors", "Matrices" of
the documentation when gsl_matrix_complex (and the like) are written
without underscore symbol.

Roman and Lazar

__________________________________
Do you Yahoo!?
New Yahoo! Photos - easier uploading and sharing.
http://photos.yahoo.com/

diff -urN linalg_dist/gsl_linalg.h linalg/gsl_linalg.h
--- linalg_dist/gsl_linalg.h	2003-07-25 11:18:21.000000000 -0400
+++ linalg/gsl_linalg.h	2003-12-25 22:51:05.000000000 -0500
@@ -109,10 +109,17 @@
                                        const gsl_vector_complex * v, 
                                        gsl_matrix_complex * A);
 
+int gsl_linalg_complex_householder_mh (gsl_complex tau, 
+                                       const gsl_vector_complex * v, 
+                                       gsl_matrix_complex * A);
+
 int gsl_linalg_complex_householder_hv (gsl_complex tau, 
                                        const gsl_vector_complex * v, 
                                        gsl_vector_complex * w);
 
+int gsl_linalg_complex_householder_hm1 (gsl_complex tau, 
+                                gsl_matrix_complex * A);
+
 /* Singular Value Decomposition
 
  * exceptions: 
diff -urN linalg_dist/householdercomplex.c linalg/householdercomplex.c
--- linalg_dist/householdercomplex.c	2003-07-25 11:18:12.000000000 -0400
+++ linalg/householdercomplex.c	2003-12-25 20:25:07.000000000 -0500
@@ -116,6 +116,58 @@
 }
 
 int
+gsl_linalg_complex_householder_mh (gsl_complex tau, const gsl_vector_complex * v, gsl_matrix_complex * A)
+{
+  /* applies a householder transformation v,tau to matrix m */
+
+  size_t i, j;
+
+  if (GSL_REAL(tau) == 0.0 && GSL_IMAG(tau) == 0.0)
+    {
+      return GSL_SUCCESS;
+    }
+
+  /* w = (v' A)^T */
+
+  for (i = 0; i < A->size1; i++)
+    {
+      gsl_complex tauwi;
+      gsl_complex wi = gsl_matrix_complex_get(A,i,0);  
+
+      for (j = 1; j < A->size2; j++)  /* note, computed for v(0) = 1 above */
+        {
+          gsl_complex Aij = gsl_matrix_complex_get(A,i,j);
+          gsl_complex vj = gsl_vector_complex_get(v,j);
+          gsl_complex Av = gsl_complex_mul (Aij,vj);
+          wi = gsl_complex_add (wi, Av);
+        }
+
+      tauwi = gsl_complex_mul (tau, wi);
+
+      /* A = A - v w^T */
+      
+      {
+        gsl_complex Ai0 = gsl_matrix_complex_get (A, i, 0);
+        gsl_complex Atw = gsl_complex_sub (Ai0, tauwi);
+        /* store A0i - tau  * wi */
+        gsl_matrix_complex_set (A, i, 0, Atw);
+      }
+      
+      for (j = 1; j < A->size2; j++)
+        {
+          gsl_complex vj = gsl_vector_complex_get (v, j);
+          gsl_complex tauvw = gsl_complex_mul(gsl_complex_conjugate(vj), tauwi);
+          gsl_complex Aij = gsl_matrix_complex_get (A, i, j);
+          gsl_complex Atwv = gsl_complex_sub (Aij, tauvw);
+          /* store Aij - tau * vi * wj */
+          gsl_matrix_complex_set (A, i, j, Atwv);
+        }
+    }
+      
+  return GSL_SUCCESS;
+}
+
+int
 gsl_linalg_complex_householder_hm (gsl_complex tau, const gsl_vector_complex * v, gsl_matrix_complex * A)
 {
   /* applies a householder transformation v,tau to matrix m */
@@ -205,3 +257,74 @@
 
   return GSL_SUCCESS;
 }
+
+
+int
+gsl_linalg_complex_householder_hm1 (gsl_complex tau, gsl_matrix_complex * A)
+{
+  /* applies a householder transformation v,tau to a matrix being
+     build up from the identity matrix, using the first column of A as
+     a householder vector */
+
+  size_t i, j;
+  gsl_complex czero = gsl_complex_rect(0.0, 0.0);
+  gsl_complex  cone = gsl_complex_rect(1.0, 0.0);
+
+  if (GSL_REAL(tau) == 0.0 && GSL_IMAG(tau) == 0.0)
+    {
+      gsl_matrix_complex_set (A, 0, 0, cone);
+      
+      for (j = 1; j < A->size2; j++)
+        {
+          gsl_matrix_complex_set (A, 0, j, czero);
+        }
+
+      for (i = 1; i < A->size1; i++)
+        {
+          gsl_matrix_complex_set (A, i, 0, czero);
+        }
+
+      return GSL_SUCCESS;
+    }
+
+  /* w = A' v */
+
+  for (j = 1; j < A->size2; j++)
+    {
+      gsl_complex wj = czero;   /* A0j * v0 */
+
+      for (i = 1; i < A->size1; i++)
+        {
+          gsl_complex vi  = gsl_matrix_complex_get(A, i, 0);
+          gsl_complex Aij = gsl_matrix_complex_get(A,i,j);
+          gsl_complex Av = gsl_complex_mul(Aij,gsl_complex_conjugate(vi));
+          wj =gsl_complex_add(wj,Av);
+        }
+
+      /* A = A - tau v w' */
+
+      gsl_complex tauwj = gsl_complex_mul(tau,wj);
+      gsl_matrix_complex_set (A, 0, j, gsl_complex_negative(tauwj));
+      
+      for (i = 1; i < A->size1; i++)
+        {
+          gsl_complex vi  = gsl_matrix_complex_get (A, i, 0);
+          gsl_complex Aij = gsl_matrix_complex_get (A, i, j);
+          vi = gsl_complex_mul(vi,wj);
+          vi = gsl_complex_mul(vi,tau);
+          vi = gsl_complex_sub(Aij,vi);
+          gsl_matrix_complex_set (A, i, j, vi);
+        }
+    }
+
+  for (i = 1; i < A->size1; i++)
+    {
+      gsl_complex vi = gsl_matrix_complex_get(A, i, 0);
+      vi = gsl_complex_mul(vi,tau);
+      gsl_matrix_complex_set(A, i, 0, gsl_complex_negative(vi));
+    }
+
+  gsl_matrix_complex_set (A, 0, 0, gsl_complex_sub(cone,tau));
+
+  return GSL_SUCCESS;
+}