Eigenvalues of a Hermitian matrix

[Andrew W Steiner <asteiner at tonic dot physics dot sunysb dot edu>]

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