# Generalized eigenvalue problem solving with Lagrange multipliers with EPSSOLVE

**URL:** <https://community.freefem.org/t/generalized-eigenvalue-problem-solving-with-lagrange-multipliers-with-epssolve/3650>\
**Category:** General Discussion\
**Created:** [December 18, 2024, 1:44pm UTC](https://community.freefem.org/t/generalized-eigenvalue-problem-solving-with-lagrange-multipliers-with-epssolve/3650 "2024-12-18T13:44:30Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![alberic](https://avatars.discourse-cdn.com/v4/letter/a/a4c791/32.png) [@alberic](https://community.freefem.org/u/alberic)\
**Post date:** [December 18, 2024, 1:44pm UTC](https://community.freefem.org/t/generalized-eigenvalue-problem-solving-with-lagrange-multipliers-with-epssolve/3650/1 "2024-12-18T13:44:30Z")

</div>

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) :

````C++
//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);
///////////////////////////////////////////////////////////////////////```
````

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [December 18, 2024, 4:57pm UTC](https://community.freefem.org/t/generalized-eigenvalue-problem-solving-with-lagrange-multipliers-with-epssolve/3650/2 "2024-12-18T16:57:34Z")

</div>

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:

```auto
" -st_type sinvert" +
" -st_pc_factor_mat_solver_type mumps " +

```

---

<div class="post-metadata">

**Author:** ![alberic](https://avatars.discourse-cdn.com/v4/letter/a/a4c791/32.png) [@alberic](https://community.freefem.org/u/alberic)\
**Post date:** [December 19, 2024, 11:43am UTC](https://community.freefem.org/t/generalized-eigenvalue-problem-solving-with-lagrange-multipliers-with-epssolve/3650/3 "2024-12-19T11:43:50Z")

</div>

It works, Thank you !