Eigenvalues of a Hermitian matrix

Eigenvalues of a Hermitian matrix

[Andrew W Steiner]

Hello!
	The following code fails to calculate the eigenvalues of the
second matrix. (It's only 70 lines, and even that would be shortened if
I wouldn't have tried to be so pedantic with my style). Is
my usage correct, or is there possibly an error in my compilation of the
libraries? I am using Red Hat Linux 6.2 with g++ (version egcs-2.91.66).
	Could someone possibly double check this for me? The code is
supposed to print out and calculate the eigenvalues of two simple
matrices but never leaves the second call to gsl_eigen_herm.

Thanks,
Andrew W. Steiner
Nuclear Theory Group
Department of Physics and Astronomy
State University of New York at Stony Brook
http://tonic.physics.sunysb.edu/~asteiner

-----------------------------------------------------------------------
#include <iostream>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_complex.h>

int main(void) {
  double p;
  int i,j,k;
  gsl_complex z, *zp;
  double real, imag;
  gsl_matrix_complex *n; 
  gsl_vector *eval;
  gsl_eigen_herm_workspace *w;

  zp=new gsl_complex;
  n=gsl_matrix_complex_alloc(4,4);
  eval = gsl_vector_alloc (4);
  w=gsl_eigen_herm_alloc(4);

  for(i=0;i<4;i++) {
    for(j=0;j<4;j++) {
      real=0.0;
      imag=0.0;
      if (i==j) real=((double)i);
      GSL_SET_COMPLEX(zp,real,imag);
      gsl_matrix_complex_set(n,i,j,*zp);
      cout << "(" << real << " " << imag << ") ";
    }
    cout << endl;
  }

  gsl_eigen_herm(n, eval, w);

  for(k=0;k<4;k++) cout << gsl_vector_get(eval,k) << " ";
  cout << endl;

  gsl_eigen_herm_free(w);
  gsl_matrix_complex_free(n);
  gsl_vector_free(eval);

  w=gsl_eigen_herm_alloc(4);
  n=gsl_matrix_complex_alloc(4,4);
  eval = gsl_vector_alloc (4);

  p=1.0;

  for(i=0;i<4;i++) {
    for(j=0;j<4;j++) {
      real=0.0;
      imag=0.0;

      if (i==j && j==1) real+=p;
      if ((i==0 && j==2) || i==2 && j==0) real-=p;
      if ((i==1 && j==3) || i==3 && j==1) real+=p;

      GSL_SET_COMPLEX(zp,real,imag);
      gsl_matrix_complex_set(n,i,j,*zp);
      cout << "(" << real << " " << imag << ") ";
    }
    cout << endl;
  }

  gsl_eigen_herm(n, eval, w);

  for(k=0;k<4;k++) cout << gsl_vector_get(eval,k) << " ";
  cout << endl;

  delete zp;

  return 0;
}

Re: Eigenvalues of a Hermitian matrix

[Brian Gough]

Andrew W Steiner writes:
 > 	The following code fails to calculate the eigenvalues of the
 > second matrix. (It's only 70 lines, and even that would be
 > shortened if I wouldn't have tried to be so pedantic with my
 > style). Is my usage correct, or is there possibly an error in my
 > compilation of the libraries? I am using Red Hat Linux 6.2 with g++
 > (version egcs-2.91.66).
 > 	Could someone possibly double check this for me? The code is
 > supposed to print out and calculate the eigenvalues of two simple
 > matrices but never leaves the second call to gsl_eigen_herm.

Thanks for the bug report.  Your program is valid. The eigenvalue
routine was going into an infinite loop.

The following patch to linalg/qrstep.c fixes it.  I will check it into
CVS.

After applying the patch,

$ ./a.out 
....
(0 0) (0 0) (-1 0) (0 0) 
(0 0) (1 0) (0 0) (1 0) 
(-1 0) (0 0) (0 0) (0 0) 
(0 0) (1 0) (0 0) (0 0) 
1 -1 1.61803 -0.618034 


Index: qrstep.c
===================================================================
RCS file: /cvs/gsl/gsl/eigen/qrstep.c,v
retrieving revision 1.2
diff -u -r1.2 qrstep.c
--- qrstep.c    2001/06/12 12:21:46     1.2
+++ qrstep.c    2001/08/02 17:15:23
@@ -60,18 +60,13 @@
 
   double mu;
 
-  if (dt > 0)
+  if (dt >= 0)
     {
       mu = tb - (tab * tab) / (dt + hypot (dt, tab));
     }
-  else if (dt < 0)
-    {
-      mu = tb + (tab * tab) / ((-dt) + hypot (dt, tab));
-    }
   else
     {
-      /* FIXME: what is the best choice of mu for dt == 0 ? */
-      mu = tb;
+      mu = tb + (tab * tab) / ((-dt) + hypot (dt, tab));
     }
 
   return mu;

Re: Eigenvalues of a Hermitian matrix

[Brian Gough]

Brian Gough writes:
 > The following patch to linalg/qrstep.c fixes it.  I will check it into
 > CVS.

There was also a bug in the accumulation of the eigenvectors.  I've
added this to the test suite.

Index: symmv.c
===================================================================
RCS file: /cvs/gsl/gsl/eigen/symmv.c,v
retrieving revision 1.2
retrieving revision 1.3
diff -u -r1.2 -r1.3
--- symmv.c	2001/07/01 21:48:38	1.2
+++ symmv.c	2001/08/02 19:10:53	1.3
@@ -192,10 +192,10 @@
                 
                 for (k = 0; k < N; k++)
                   {
-                    double qki = gsl_matrix_get (evec, k, i);
-                    double qkj = gsl_matrix_get (evec, k, i + 1);
-                    gsl_matrix_set (evec, k, i, qki * c - qkj * s);
-                    gsl_matrix_set (evec, k, i + 1, qki * s + qkj * c);
+                    double qki = gsl_matrix_get (evec, k, a + i);
+                    double qkj = gsl_matrix_get (evec, k, a + i + 1);
+                    gsl_matrix_set (evec, k, a + i, qki * c - qkj * s);
+                    gsl_matrix_set (evec, k, a + i + 1, qki * s + qkj * c);
                   }
               }
             
Index: hermv.c
===================================================================
RCS file: /cvs/gsl/gsl/eigen/hermv.c,v
retrieving revision 1.3
retrieving revision 1.4
diff -u -r1.3 -r1.4
--- hermv.c	2001/07/01 21:48:38	1.3
+++ hermv.c	2001/08/02 19:10:53	1.4
@@ -215,8 +215,8 @@
                 
                 for (k = 0; k < N; k++)
                   {
-                    gsl_complex qki = gsl_matrix_complex_get (evec, k, i);
-                    gsl_complex qkj = gsl_matrix_complex_get (evec, k, i + 1);
+                    gsl_complex qki = gsl_matrix_complex_get (evec, k, a + i);
+                    gsl_complex qkj = gsl_matrix_complex_get (evec, k, a + i + 1);
                     /* qki <= qki * c - qkj * s */
                     /* qkj <= qki * s + qkj * c */
                     gsl_complex x1 = gsl_complex_mul_real(qki, c);
@@ -228,8 +228,8 @@
                     gsl_complex qqki = gsl_complex_add(x1, y1);
                     gsl_complex qqkj = gsl_complex_add(x2, y2);
                     
-                    gsl_matrix_complex_set (evec, k, i, qqki);
-                    gsl_matrix_complex_set (evec, k, i + 1, qqkj);
+                    gsl_matrix_complex_set (evec, k, a + i, qqki);
+                    gsl_matrix_complex_set (evec, k, a + i + 1, qqkj);
                   }
               }
             

Eigenvalues of a Hermitian matrix

[Andrew W Steiner]

Hello!
	I want to calculate the derivative of a function of the
eigenvalues of a particular Hermitian matrix and so I need the eigenvalues
of several matrices whose entries are similar. Is there any way to
provide an inital guess to the gsl_eigen_herm routine if one is known?
	From my understanding of the algorithm, the answer is no. In that
case, could anyone suggest an alternative? 

Thanks,
Andrew W. Steiner
Nuclear Theory Group
Department of Physics and Astronomy
State University of New York at Stony Brook
http://tonic.physics.sunysb.edu/~asteiner

Re: Eigenvalues of a Hermitian matrix

[Brian Gough]

Andrew W Steiner writes:
 > 	I want to calculate the derivative of a function of the
 > eigenvalues of a particular Hermitian matrix and so I need the eigenvalues
 > of several matrices whose entries are similar. Is there any way to
 > provide an inital guess to the gsl_eigen_herm routine if one is known?
 > 	From my understanding of the algorithm, the answer is no. In that
 > case, could anyone suggest an alternative? 

If you premultiply your perturbed matrices by V' (where V is the
unitary matrix of unperturbed eigenvectors) they will be close to
diagonal.  This might improve the computation of the eigenvalues.  The
same applies if you have only a guess for V.

It probably won't make much difference overall, but if it does let me
know.

regards
Brian

Re: Eigenvalues of a Hermitian matrix

[Brian Gough]

Brian Gough writes:
 > 
 > If you premultiply your perturbed matrices by V' (where V is the
          ^^^^^^^^^^^^
Should be "apply a similarity transformation, V' M V", not
premultiply, of course.