* Complex Householder transform functions
2003-12-08 15:38 ` Brian Gough
@ 2003-12-26 4:19 ` Z F
0 siblings, 0 replies; 6+ messages in thread
From: Z F @ 2003-12-26 4:19 UTC (permalink / raw)
To: gsl-discuss
[-- Attachment #1: Type: text/plain, Size: 679 bytes --]
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/
[-- Attachment #2: householder.diff --]
[-- Type: text/plain, Size: 4899 bytes --]
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;
+}
^ permalink raw reply [flat|nested] 6+ messages in thread