Obtaining the eigenvector from EPSSolve

General Discussion
1

I’m able to solve an eigenvalue problem using EPSSolve, but I’m having trouble obtaining the calculated eigenvector and printing it on the screen. In the following code, how could I store the computed eigenvector in my eV object?
Thanks in advance

load "PETSc"

// Mesh
mesh Th = square(9, 9, [x*5.0,y*5.0]);

// FE Space
fespace Vh(Th, P1);
Vh u, v;

// assemble matrix A
matrix[int] Af0(2), Af1(2);
varf a11(u,v) = int2d(Th)( 1.4*( dx(u)*dx(v) + dy(u)*dy(v) ) + 0.02*u*v );
varf a12(u,v) = int2d(Th)( 1*u*v );
varf a21(u,v) = int2d(Th)( 1*u*v );
varf a22(u,v) = int2d(Th)( 0.4*( dx(u)*dx(v) + dy(u)*dy(v) ) + 0.15*u*v );
Af0[0] = a11(Vh,Vh);
Af0[1] = a12(Vh,Vh);
Af1[0] = a21(Vh,Vh);
Af1[1] = a22(Vh,Vh);
matrix A = [ [ Af0[0], Af0[1] ],
	[ Af1[0], Af1[1] ] ];

// assemble matrix B
matrix[int] Bf0(2), Bf1(2);
varf b11(u,v) = int2d(Th)( 0.007*u*v );
varf b12(u,v) = int2d(Th)( 0.200*u*v );
varf b21(u,v) = int2d(Th)( 0*u*v );
varf b22(u,v) = int2d(Th)( 0*u*v );
Bf0[0] = b11(Vh,Vh);
Bf0[1] = b12(Vh,Vh);
Bf1[0] = b21(Vh,Vh);
Bf1[1] = b22(Vh,Vh);
matrix B = [ [ Bf0[0], Bf0[1] ],
	[ Bf1[0], Bf1[1] ] ];

// solve eigenvalue problem
int nev = 1;
real[int] lambda(nev);
Vh[int] eV(nev); // what type ??

Mat dsA(A);
Mat dsB(B);
string ssparams = " -eps_type krylovschur" +
	" -eps_target 0" +
	" -eps_nev" + nev +
	" -st_type sinvert" +
	" -st_pc_type lu" +
	" -eps_smallest_magnitude";

EPSSolve(dsA, dsB, sparams = ssparams, values = lambda, vectors = eV);

cout << " eigenvalue = " << lambda(0) << endl;
cout << " eigenvector = " << eV[0][] << endl;
2

It can’t be of type Vh since your matrix is of dimension 2 * Vh.ndof. So use instead real[int,int] array(2 * Vh.ndof, nev) and then EPSSolve(..., array = array). And then you can pull the respective rows to get each component.

3

I tried
real[int,int] array(2 * Vh.ndof, nev)
but it didn’t work.

Using the following I could store the eigenvector in the variable eV1:

fespace Wh(Th, [P1, P1]);
Wh[int] [eV1, eV2](nev);
...
EPSSolve(..., vectors = eV1); // solve the eigenvalue equation
cout << "eigenvector = " << eV1[0][]; // print the eigenvector

The notation here seems kind of weird because eV1 and eV2 are actually the same arrays and so eV2 might be used in the last two lines of the following code instead of eV1 with the same effect

4

What did not work with real[int,int] array(2 * Vh.ndof, nev)?

5

It says:
Impossible to cast <P3KNMIdE> in <13FEbaseArrayKnIdE> ... at line EPSSolve(..., vectors = array);

6

Of course, that is not expected to work, you need to read more carefully my initial answer.

7

Oh, I mistakenly was trying EPSSolve(..., vectors = array) instead of EPSSolve(..., array = array) that resulted in error.
Statement real[int,int] array(2 * Vh.ndof, nev) works correctly along with EPSSolve(..., array = array), and array is the sought eigenvector that can be printed or used in any other way.