I'm trying to obtain the eigenvectors and eigenvalues of
a Hermitian matrix (complex nxn matrix) using zheev_ of LAPACK (in C).
Everything works fine up to a certain matrix size, n <= 626, but when
n > 626, a segmentation fault keeps popping up. I'm pretty sure I don't need any more RAM. Here is the part of the code
- Code: Select all
#include <stdio.h>
#include <stdlib.h>
#include <complex.h>
#include <math.h>
#include "ran3.c" //random number generator
int main(void)
{
//parameters for zheev_()
double *HT;
int LWORK,LRWORK;
double *w;
complex double *WORK;
double *RWORK;
long idum = time(0);
///////////////////////////////////////////////////////
//Change this according to the dimension of your matrix
//H: Matrix (nxn), U: Eigenvector (column vectors)
// H[n][n]
complex double H[627][627],U[627][627];
///////////////////////////////////////////////////////
int i,j,k,n,s,ok,size;
char c4,c5;
c4='V';
c5='L';
///////////////////////////////////////////////////////
// Change this according to the dimension of the matrix
n = 627;
///////////////////////////////////////////////////////
// Random nxn Hermitian matrix
for(i=0;i<n;i++) {
for(j=0;j<=i;j++) {
H[i][j]=ran3(&idum)+I*ran3(&idum);
H[j][i]=conj(H[i][j]);
}
}
///////////////////////////////////////////////////////////////
HT = (double *)malloc(sizeof (double)*n*n*2);
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
HT[2*(n*i+j)]=creal(H[j][i]);
HT[2*(n*i+j)+1]=cimag(H[j][i]);
}
w=(double *)malloc(sizeof(double)*n);
LWORK=5*n;
WORK=(complex double *)malloc(sizeof(complex double)*LWORK);
LRWORK=3*n;
RWORK= (double *)malloc(sizeof (double)*LRWORK);
// Solving the linear system
zheev_(&c4,&c5,&n,HT,&n,w,WORK,&LWORK,RWORK,&ok);
if(ok!=0)printf("diagonalization failed\n");
// Eigenvalues in ascending order
for(i=0;i<n;i++)
printf("%lf\t",w[i]);
printf("\n");
// copying eigenvectors
for(i=0;i<size;i++)
for(j=0;j<size;j++)
U[j][i]=HT[2*(n*i+j)]+I*HT[2*(n*i+j)+1];
for(i=0;i<n;i++) {
for(j=0;j<n;j++)
{printf("%lf+I*%lf\t",creal(U[i][j]),cimag(U[i][j]));}
printf("\n");
}
printf("Okay!");
}
Here H is the nxn Hermitian matrix, w is the array containing the eigenvalues, and U is the array containing the eigenvectors. I'm a beginner in C and LAPACK so any suggestions or solutions would be greatly appreciated. Thank you.
