Generalized eigenvalue problem solving with Lagrange multipliers with EPSSOLVE

General Discussion
1

Hello,
I am facing a problem with a generalized eigenvalues problem using EPSSolve in FreeFEM.

I want to solve a variational formulation with constraints, which I handle using Lagrange multipliers. This leads to a matrix problem of the form AX = \lambda BX. (B has some lines of zeros but it shouldn’t be a problem to solve this matricial equation)

I can solve this problem in MATLAB, but I cannot get it to work with the EPSSOLVE function in FreeFEM.

If I specify the parameters for EPSSOLVE as sparams="-eps_type lapack", my code runs (very slowly), but the returned eigenvector is clearly incorrect.
Do you have any idea ?
Here is some code if you want to test (sorry apparently i can’t import a file yet) :

//mpirun -np 1 FreeFem++-mpi -wg main_tmp_filtre.edp -v 0

load "PETSc"
//load "SLEPc"

bool plots =1; // 1 if plots are used, 0 otherwise

///////////////////////////////////////////////////////////////////////
// Declaring and reading external parameters                         //
///////////////////////////////////////////////////////////////////////

// Parameters defining the domain
real L, Lx, Ly;
// Parameters defining the fine mesh
int Ndiscr, Nx, Ny;
// Parameters defining the coarse mesh
real eps;



///////////////////////////////////////////////////////////////////////
// Define internal data structures                                   //
///////////////////////////////////////////////////////////////////////
L= 1.; Lx=L; Ly=L;
Ndiscr =32; Nx=Ndiscr; Ny=Ndiscr;
eps= 1;

real hx=Lx/Nx, hy=Ly/Ny; // Fine mesh Th
mesh Th=square(Nx,Ny,[Lx*x,Ly*y]); //fine mesh
fespace Vh(Th,P1); //fine global FE space

////////// Diffusion coefficient //////////
string ParameterDescription;

func Aeps = 1+0*x*y;

func sigmaeps = 1+0*x*y;

//func FiltrePhi= 1;
func FiltrePhi= (1-x)*(1-x)*x*x*(1-y)*(1-y)*y*y*630; // Fontion unitaire nulle aux  bords pour la norme L2 dans la cellule
////////////////////////////////////////////////////////////////

macro a(dif,u,v) (dif*dx(u)*dx(v) + dif*dy(u)*dy(v))//EOM


// with Filtrage
varf oploc (u, v)
        = int2d(Th)(
            FiltrePhi*sigmaeps(x,y)*u*v)
        + int2d(Th)(eps^2*FiltrePhi*a(Aeps(x,y),u,v))
        ;

varf bloc (u, v) = int2d(Th)(FiltrePhi*u*v); //no boundary condition


varf CLXfiltreloc (unused, v)
    = int2d(Th)(FiltrePhi*1*dx(v)) 
    ;

varf CLYfiltreloc (unused, v)
    = int2d(Th)(FiltrePhi*1*dy(v)) 
    ;

real[int] CLXfiltreloclocS1 = CLXfiltreloc(0, Vh);
real[int] CLYfiltreloclocS1 = CLYfiltreloc(0, Vh);

real[int] Zeros = 0*CLXfiltreloclocS1;
//cout << "CRlocS1z = " << CRlocS1z.n << endl;

matrix AlocS = oploc(Vh, Vh, factorize=1); 
matrix BlocS = bloc(Vh, Vh, solver=CG, eps=1e-20);

matrix AlocStmp(AlocS.n+2,AlocS.m+2);
matrix BlocStmp(BlocS.n+2,BlocS.m+2);
AlocStmp = [ [AlocS , CLXfiltreloclocS1 , CLYfiltreloclocS1] , [CLXfiltreloclocS1' , 0, 0] , [CLYfiltreloclocS1' , 0 ,0] ];
BlocStmp = [ [BlocS , Zeros, Zeros ] , [Zeros' , 0, 0 ] , [Zeros' , 0 , 0] ];

Mat Aloctemp = AlocStmp;
Mat Bloctemp = BlocStmp; 

real[int] Listelambdaloc(1); //to store the nev eigenvalue (nev = 1 ici)
real[int,int] ListeVecteurproprelocS(AlocStmp.n,1) ; //to store the nev eigenvector

string ssparamsfiltre =            // Parameters for the distributed EigenValue solver
" -eps_nev " + 1       + // Number of eigenvalues
" -eps_type krylovschur" +
" -eps_target "+ 0   + // Shift value
//" -st_type sinvert "     +
" -eps_gen_hermitian"      // The problem is symmetric
;

//int kloc = EPSSolve(Aloctemp, Bloctemp, values=Listelambdaloc , array=ListeVecteurproprelocS, sparams= "-eps_type lapack"); //Résolution du pb aux vp
int kloc = EPSSolve(Aloctemp, Bloctemp, values=Listelambdaloc , array=ListeVecteurproprelocS, sparams= ssparamsfiltre); //Résolution du pb aux vp
kloc=min(kloc,1);

cout << "VP" << kloc << " = " << Listelambdaloc(0) << endl;
Vh VecteurproprelocS;

VecteurproprelocS[](0:VecteurproprelocS.n-1) = ListeVecteurproprelocS(0:VecteurproprelocS.n-1,0);
real LagrangeMultiplier1 = ListeVecteurproprelocS(VecteurproprelocS.n,0);
real LagrangeMultiplier2 = ListeVecteurproprelocS(VecteurproprelocS.n+1 ,0);


// Renormalisation
real normeL2VecteurproprelocS=sqrt(int2d(Th)(VecteurproprelocS*VecteurproprelocS));
VecteurproprelocS[] = VecteurproprelocS[]/normeL2VecteurproprelocS;
///////////////////////////////////////////////////////////////////////

cout << "LagrangeMultiplier1 : " << LagrangeMultiplier1 << "LagrangeMultiplier2 : " << LagrangeMultiplier2 << endl;
if (plots) plot(VecteurproprelocS, fill=1, value=1, cmm="Premier vecteur propre", wait=1);
///////////////////////////////////////////////////////////////////////```
1 Like
2

This was the topic of last Friday’s expert tutorial during the FreeFEM days, so there will be a screencast available soon on the YouTube channel that explains this. In the meantime, if you are looking for the smallest eigenpair, you probably want to use:

" -st_type sinvert" +
" -st_pc_factor_mat_solver_type mumps " +
3

It works, Thank you !