Nonsymmetric Eigenproblems

Nonsymmetric Eigenproblems

[Nils Wagner]

Hi all,

I am wondering if there is any effort concerning the handling of
standard and generalized eigenvalue problems with
nonsymmetric matrices.

Afaik, there are only routines for Hermitian (complex)
symmetric (real case) matrices in GSL-1.4.

Any comments or suggestions ?

Nils

Re: Nonsymmetric Eigenproblems

[Kenneth Geisshirt]

Nils Wagner wrote:

> I am wondering if there is any effort concerning the handling of
> standard and generalized eigenvalue problems with
> nonsymmetric matrices.
> 
> Afaik, there are only routines for Hermitian (complex)
> symmetric (real case) matrices in GSL-1.4.
> 
> Any comments or suggestions ?

Back in 1994/5 I ported some routines from Pascal (implementation due to 
my thesis supervisor who ported it from Algol). I haven't used it for 
years. The description in the header file says:

   Eigen is a library for computing eigenvalues and eigenvectors of
   general matrices. There is only one routine exported, namely Eigen.

   The meaning of the arguments to Eigen is:
     1.   The dimension of the general matrix (n).
     2.   A general matrix (A).
     3.   The maximal number of iterations.
     4.   The precision.
     5.   A vector with the eigenvalues.
     6.   A matrix with the eigenvectors.

The source code is around 1000 lines of ISO C. I would be happy to 
contribute it to GSL.

Kneth

-- 
Kenneth Geisshirt, M.Sc., Ph.D.         http://kenneth.geisshirt.dk
Grøndals Parkvej 2A, 3. sal                    kenneth@geisshirt.dk
DK-2720 Vanløse                                     +45 38 87 78 38

Re: Nonsymmetric Eigenproblems

[Brian Gough]

Kenneth Geisshirt writes:
 > The source code is around 1000 lines of ISO C. I would be happy to 
 > contribute it to GSL.

Sounds interesting --- I'd be very interested to see it.

-- 
Brian

Re: Nonsymmetric Eigenproblems

[Kenneth Geisshirt]

Brian Gough wrote:
> Kenneth Geisshirt writes:
>  > The source code is around 1000 lines of ISO C. I would be happy to 
>  > contribute it to GSL.
> 
> Sounds interesting --- I'd be very interested to see it.
> 

I have been reading a bit this weekend (both books and code), and it 
seems that the standard method is to transform a general matrix to a 
hermitian matrix.

Since GSL has routines for hermitian matrices, I think I could implement 
the transformation stuff. I'll be back!

Kneth
-- 
Kenneth Geisshirt, M.Sc., Ph.D.         http://kenneth.geisshirt.dk
Grøndals Parkvej 2A, 3. sal                    kenneth@geisshirt.dk
DK-2720 Vanløse                                     +45 38 87 78 38

Re: Nonsymmetric Eigenproblems

[Brian Gough]

Kenneth Geisshirt writes:
 > I have been reading a bit this weekend (both books and code), and it 
 > seems that the standard method is to transform a general matrix to a 
 > hermitian matrix.
 > 
 > Since GSL has routines for hermitian matrices, I think I could implement 
 > the transformation stuff. I'll be back!

I'm only familiar with the method described in Golub and Van Loan --
be sure to study that. I think it is the best book on the subject.

Complex Householder transform functions

[Z F]

[-- 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;
+}