Computing eigenvalues

General Discussion
1

Hello,

could anybody explain me why arpack is not able to compute the eigenvalues in this example?

matrix C = [[ 2, -1, -1, 0],
[-1, 1, 0, 0],
[-1, 0, 2, -1],
[ 0, 0, -1, 1]];

set(C,solver=UMFPACK);

func real[int] FOP(real[int] & u){ real[int] Au = C^-1*u; return Au; }
int nev = 2; //number of computed eigenvalues
real[int] ev(nev);
int mk = EigenValue(4,FOP,mode=3,sym=sym,nev=nev,value=ev);

Thanks,
Ernesto

2

This issue has been raised multiple times: here is a thorough answer.

3

Here is what I get for your specific problem.

matrix C = [[ 2, -1, -1, 0],
[-1, 1, 0, 0],
[-1, 0, 2, -1],
[ 0, 0, -1, 1]];

load "SLEPc"
Mat mC(C);
int nev = 2; // number of computed eigenvalues
real[int] ev(nev);
EPSSolve(mC, sparams = "-eps_nev "+nev+" -eps_largest_magnitude", values = ev);
ev.resize(min(ev.n, nev));
cout << ev << endl;
EPSSolve(mC, sparams = "-eps_nev "+nev+" -eps_smallest_magnitude", values = ev);
ev.resize(min(ev.n, nev));
cout << ev << endl;

load "lapack"
real[int, int] fC(4, 4);
for[i,j,v : C] fC(i, j) = v;
complex[int, int] vectors(4, 4);
complex[int] values(4);
dgeev(fC, values, vectors);
cout << values << endl;

$ mpirun -np 1 FreeFem++-mpi ev.edp -v 0 
2
	3.414213562	  2
2
	-1.498580728e-16	0.5857864376
4
	(3.414213562,0)	(2,0)	(3.474907126e-17,0)	(0.5857864376,0)
4

the basic way with arpack in freefem++
is to do some think like

matrix C = [
[ 2, -1, -1, 0],
[-1, 1, 0, 0],
[-1, 0, 2, -1],
[ 0, 0, -1, 1]];
matrix Id=eye(4);
real s = 1; // do a shift because the matrix C is not invertible
C -= s*Id;
set(C,solver=UMFPACK);

 int nev = 3; //number of computed eigenvalues
real[int] ev(nev);
real[int,int] EV(4,nev);
int mk = EigenValue(C,Id,sigma=s,sym=1,nev=nev,value=ev,rawvector=EV,mode=3);

 cout << ev << endl; 

 cout << EV << endl;
5

Note that this won’t work if you change nev from 3 to 4.

SLEPc: 4
	3.414213562	  2	0.5857864376	-1.498580728e-16
ARPACK: 4
	2.220446049e-16	0.5857864376	  2	  0

Be careful with ARPACK…