load "lapack"

//  laplace with matrix  
verbosity=1;
int nn=5;
mesh Th=square(nn,nn,[2*pi*x,2*pi*y]);
fespace Vh(Th,P2);
Vh<complex> u1,u2;
Vh ur,ui;
int n=u1.n;
complex[int] Bu1(n),Bu2(n);
complex[int] Au1(n),Au2(n);

real sigma=0.0;                     
complex nu=1.0+1i;

varf op(u1,u2)= int2d(Th)( nu*(dx(u1)*dx(u2) + dy(u1)*dy(u2)) + (dx(u1)+dy(u1))*u2 - sigma* u1*u2 )
              //+ on(1,2,3,4,u1=0)
              ;
varf  a(u1,u2)= int2d(Th)( nu*(dx(u1)*dx(u2) + dy(u1)*dy(u2)) + (dx(u1)+dy(u1))*u2 )
              //+ on(1,2,3,4,u1=0)
              ;
                   
varf b(u1,u2) = int2d(Th)(  u1*u2 
              + 1i*u1*dx(u2)
              );

matrix<complex> OP=op(Vh,Vh,solver=sparsesolver); 
matrix<complex> A=a(Vh,Vh,solver=sparsesolver); 
matrix<complex> B=b(Vh,Vh,solver=sparsesolver); 
int nev=10;

complex[int] ev(nev); // to store eigenvalues 
Vh<complex>[int] eV(nev);   // to store eigenvectors  
  
//int k=EigenValue(OP,B,sigma=sigma,value=ev,vector=eV,tol=1e-10,maxit=0,ncv=0); // the default driver is number 4 (M positive definite Hermitian)
int k=EigenValue(OP,B,sigma=sigma,value=ev,vector=eV,tol=1e-10,maxit=0,ncv=0,driver=2); // driver 2 works when M is arbitrary (non Hermitian)
k=min(k,nev); //  some time the number of converged eigen value 
              // can be greater than nev;
for (int kk=0;kk<k;kk++)
  { 
  int i=kk;	
  u1=eV[i];	
  complex v= ev[i];
  Bu1=B*u1[];
  Au1=A*u1[];
  // The Rayleigh quotient lambda = x'Ax/x'Bx 
  // gives the eigenvalue
  // Vincent Laude, 12/01/2017: don't we need the left and right eigenvectors to form the Rayleigh quotient in the general case? 
  complex xAx = u1[]'*Au1 ;
  complex xBx = u1[]'*Bu1 ;
  //  A u = l * B u
  //Bu1 =   v*Bu1;
  //u1[] = Au1 -Bu1; 
  // Vincent Laude, 12/01/2017: why is a zero sum a good indicator of nullity of a vector? Don't we need a norm instead?
  //if(norm(u1[].sum)>1e-5)
  //      cout << "BUG :::   zero ==  " <<u1[].sum << endl;
  cout << " ---- " <<  i<< " " <<  v << " Rayleigh quotient difference: " << abs(v-xAx/xBx) << endl;
  ur=real(eV[i]);
  ui=imag(eV[i]);
  //  plot(ur,cmm="Eigen  Vector (real)  "+i+" valeur =" + v  ,wait=1,value=1);
  //  plot(ui,cmm="Eigen  Vector (imag)  "+i+" valeur =" + v  ,wait=1,value=1);

  // FFCS: add 3d view
  //plot(eV[i],dim=3,fill=1);
  }
 
// Comparison with lapack (full matrices)
// Doesn't seem to work when the on() condition is used (because of tgv)
complex[int, int] Afull(Vh.ndof,Vh.ndof), Bfull(Vh.ndof,Vh.ndof), EVl(Vh.ndof,Vh.ndof), EVr(Vh.ndof,Vh.ndof);
complex[int] eva(Vh.ndof), evb(Vh.ndof);
int[int] I(1), J(1);
complex[int] C(1);
[I,J,C]=A;
for(int i=0; i<I.n; i++) Afull(I(i),J(i))=C(i);
//cout << Afull << endl;
[I,J,C]=B;
for(int i=0; i<I.n; i++) Bfull(I(i),J(i))=C(i);
//zggev(Afull,Bfull,eva,evb,EVl,EVr);
zggev(Afull,Bfull,eva,evb,EVr);
cout << "\n\n Eigenvalues with LAPACK" << endl;
for (int i=0;i<Vh.ndof;i++) cout << i << " " << eva(i)/evb(i) << endl;
cout << endl;


// Comparison with slepc
// Note: only sigma=0.0 case can be tested
/*
// export matrix A in ascii Morse format
{
  ofstream file("matA.txt");
  file.scientific; 
  file << A << endl;
};	
// converts matrix A to PETSC binary format
exec("ffpp2petsc -fin matA.txt -fout matA.petsc");
// export matrix B in ascii Morse format
{
  ofstream file("matB.txt");
  file.scientific; 
  file << B << endl;
};	
// converts matrix B to PETSC binary format
exec("ffpp2petsc -fin matB.txt -fout matB.petsc");

exec("inveigenpb -f1 matA.petsc -f2 matB.petsc -nval "+nev); 
*/

