# Problem for computing eigenvalues with EPSSolve

**URL:** <https://community.freefem.org/t/problem-for-computing-eigenvalues-with-epssolve/2420>\
**Category:** General Discussion\
**Created:** [April 26, 2023, 7:32am UTC](https://community.freefem.org/t/problem-for-computing-eigenvalues-with-epssolve/2420 "2023-04-26T07:32:46Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Loic](https://avatars.discourse-cdn.com/v4/letter/l/e9c0ed/32.png) [@Loic](https://community.freefem.org/u/Loic)\
**Post date:** [April 26, 2023, 7:32am UTC](https://community.freefem.org/t/problem-for-computing-eigenvalues-with-epssolve/2420/1 "2023-04-26T07:32:46Z")

</div>

Dear all,

I am facing a problem with a generalized eigenvalues problem using `EPSSolve` in FreeFEM of the type `MASS \Phi = \lambda STIFF \Phi`.

`EPSSolve` seems to diverge . Indeed, I have to find 14 eigenvalues. However the two first eigenvalues are big and negative while all my eigenvalues should be positive.

I compare the values with eigenvalues computed in `Python`. With Python, all my eigenvalues are positive. Among my 14 eigenvalues, 12 of them are similar with `Python` and `FreeFEM`.

By the way, the eigenvectors found with the two method are different.

I attach the code, so that you can reproduce the problem.  
[Problem\_EV.zip](https://community.freefem.org/uploads/short-url/xX0Tv5tmlmjPCMY9V0u8B4juLXl.zip) (9.5 KB)

To compare the two method, you have just to do :

`mpiexec -np 1 FreeFem++-mpi FF_COMPUTE_EV.edp`  
`python PY_COMPUTE_EV.py`

Thank you in advance for your help,

Best regards,

Loïc,

---

<div class="post-metadata">

**Author:** ![aszaboa](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/aszaboa/32/2918_2.png) [@aszaboa](https://community.freefem.org/u/aszaboa)\
**Post date:** [April 26, 2023, 7:45am UTC](https://community.freefem.org/t/problem-for-computing-eigenvalues-with-epssolve/2420/2 "2023-04-26T07:45:27Z")

</div>

I think you should look at all the `EPSSolve` options in the SLEPc manual, and that will help you sort out the problem. My guess is that you need to specify `-eps_gen_hermitian` since you are solving a generalized eigenvalue problem.

---

<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:** [April 26, 2023, 8:09am UTC](https://community.freefem.org/t/problem-for-computing-eigenvalues-with-epssolve/2420/3 "2023-04-26T08:09:10Z")

</div>

Also, it’s perfectly normal that the eigenvectors are not the same, as long as they span the same space.

---

<div class="post-metadata">

**Author:** ![Loic](https://avatars.discourse-cdn.com/v4/letter/l/e9c0ed/32.png) [@Loic](https://community.freefem.org/u/Loic)\
**Post date:** [April 26, 2023, 8:17am UTC](https://community.freefem.org/t/problem-for-computing-eigenvalues-with-epssolve/2420/4 "2023-04-26T08:17:36Z")

</div>

Hi,

Thanks for your reply,

However when I use this parameter, I get

```auto
loading matrices
[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Missing or incorrect user input
[0]PETSC ERROR: The inner product is not well defined: indefinite matrix
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.18.3, Dec 28, 2022 
[0]PETSC ERROR: FreeFem++-mpi on a named is244532 by lb265317 Wed Apr 26 10:16:06 2023
[0]PETSC ERROR: Configure options --prefix=/volatile/catB/INSTALL/FREEFEM2/ff-petsc/r --with-debugging=0 COPTFLAGS="-O3 -mtune=native" CXXOPTFLAGS="-O3 -mtune=native" FOPTFLAGS="-O3 -mtune=native" --with-cxx-dialect=11 --with-ssl=0 --with-x=0 --with-fortran-bindings=0 --with-cudac=0 --download-slepc --download-hpddm --download-slepc-commit=d9df183c341775ffc91510ec6827e13db44a25cd --with-cc=mpicc --with-cxx=mpicxx --with-fc=mpif90 --with-scalar-type=real --with-blaslapack-include= --with-blaslapack-lib="-llapack -lblas" --download-metis --download-ptscotch --download-hypre --download-parmetis --download-superlu --download-suitesparse --download-tetgen --download-slepc-configure-arguments=--download-arpack --download-scalapack --download-mumps PETSC_ARCH=fr
[0]PETSC ERROR: #1 BV_SafeSqrt() at /home/catB/lb265317/FreeFem-sources/3rdparty/ff-petsc/petsc-3.18.3/fr/externalpackages/git.slepc/include/slepc/private/bvimpl.h:132
[0]PETSC ERROR: #2 BV_SquareRoot_Default() at /home/catB/lb265317/FreeFem-sources/3rdparty/ff-petsc/petsc-3.18.3/fr/externalpackages/git.slepc/include/slepc/private/bvimpl.h:362
[0]PETSC ERROR: #3 BVOrthogonalizeCGS1() at /home/catB/lb265317/FreeFem-sources/3rdparty/ff-petsc/petsc-3.18.3/fr/externalpackages/git.slepc/src/sys/classes/bv/interface/bvorthog.c:101
[0]PETSC ERROR: #4 BVOrthogonalizeGS() at /home/catB/lb265317/FreeFem-sources/3rdparty/ff-petsc/petsc-3.18.3/fr/externalpackages/git.slepc/src/sys/classes/bv/interface/bvorthog.c:177
[0]PETSC ERROR: #5 BVOrthogonalizeColumn() at /home/catB/lb265317/FreeFem-sources/3rdparty/ff-petsc/petsc-3.18.3/fr/externalpackages/git.slepc/src/sys/classes/bv/interface/bvorthog.c:333
[0]PETSC ERROR: #6 EPSGetStartVector() at /home/catB/lb265317/FreeFem-sources/3rdparty/ff-petsc/petsc-3.18.3/fr/externalpackages/git.slepc/src/eps/interface/epssolve.c:803
[0]PETSC ERROR: #7 EPSSolve_KrylovSchur_Default() at /home/catB/lb265317/FreeFem-sources/3rdparty/ff-petsc/petsc-3.18.3/fr/externalpackages/git.slepc/src/eps/impls/krylov/krylovschur/krylovschur.c:242
[0]PETSC ERROR: #8 EPSSolve() at /home/catB/lb265317/FreeFem-sources/3rdparty/ff-petsc/petsc-3.18.3/fr/externalpackages/git.slepc/src/eps/interface/epssolve.c:146
0

```

Loïc,

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [April 26, 2023, 9:44am UTC](https://community.freefem.org/t/problem-for-computing-eigenvalues-with-epssolve/2420/5 "2023-04-26T09:44:16Z")

</div>

I guess I could dump the matricies that have been read in and first insure that  
the FF edp file read them properly or visualize them or just dump them back  
out and diff the input and output txt files. I know in one case writing out  
a matrix from FF I had to flush the file stream or else last value was dropped.

In phonons anyway you sometimes  
see small values with the wrong sign if there is no  
restoring force- the materials just falls apart.  
In this case, what would large negative eignvalues mean  
in your system? If you look at the eigenvectors for the negative ones  
what does that suggest about your system?

I did note some of these are nearly degenerate and the solver AFAIK has  
no way to pick particular vectors to span the spaces.

---

<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:** [April 26, 2023, 1:49pm UTC](https://community.freefem.org/t/problem-for-computing-eigenvalues-with-epssolve/2420/6 "2023-04-26T13:49:37Z")

</div>

BTW, if I run your FreeFEM script and dump the matrices, then with MATLAB, I can check that the results from SLEPc are correct.

```auto
>> A = PetscBinaryRead('/tmp/A0.bin');
>> B = PetscBinaryRead('/tmp/A1.bin');
>> eig(full(A), full(B))

ans =

    -3.452024300482933e+11
    -5.080069245105576e+11
     1.140810058050668e-02
     1.142933265808896e-02
     5.482857616842862e-03
     5.479991885016102e-03
     2.119836942750551e-03
     2.164645146026448e-03
     1.522260102485258e-03
     3.599930347392806e-04
     5.027322264954978e-04
     1.253948212699970e-03
     8.058577249813228e-04
     7.976357767610110e-04

>> eigs(A, B)

ans =

    -5.080069245105576e+11
    -3.452024300482933e+11
     1.142933265808896e-02
     1.140810058050668e-02
     5.482857616842862e-03
     5.479991885016102e-03

```

---

<div class="post-metadata">

**Author:** ![Loic](https://avatars.discourse-cdn.com/v4/letter/l/e9c0ed/32.png) [@Loic](https://community.freefem.org/u/Loic)\
**Post date:** [April 26, 2023, 2:25pm UTC](https://community.freefem.org/t/problem-for-computing-eigenvalues-with-epssolve/2420/7 "2023-04-26T14:25:48Z")

</div>

Thanks for testing,

It is strange that I find other eigenvalues with Python, especially since both my matrices should de positive definite.

Loïc,

---

<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:** [April 26, 2023, 2:37pm UTC](https://community.freefem.org/t/problem-for-computing-eigenvalues-with-epssolve/2420/9 "2023-04-26T14:37:47Z")

</div>

Your matrix B is rank deficient, so it’s not that simple.

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [April 26, 2023, 2:45pm UTC](https://community.freefem.org/t/problem-for-computing-eigenvalues-with-epssolve/2420/10 "2023-04-26T14:45:02Z")

</div>

Are these the ev’s you got in python? I put your python coefs into  
FF ,

```auto
type or paste code here
`` ff-mpirun -np 1 FF_COMPUTE_EV.edp | grep lambda
   35 : cout << "lambda[" << j << "] = " << ev(j) << " and " << eV(:,j) << endl;
   40 : I am facing a problem with a generalized eigenvalues problem using EPSSolve in FreeFEM of the type MASS \Phi = \lambda STIFF \Phi.
lambda[0] = 177623 and 14	
lambda[1] = 177622 and 14	
lambda[2] = 0.0114293 and 14	
lambda[3] = 0.0114082 and 14	
lambda[4] = 0.00548286 and 14	
lambda[5] = 0.00547999 and 14	
lambda[6] = 0.00216465 and 14	
lambda[7] = 0.00211984 and 14	
lambda[8] = 0.00152226 and 14	
lambda[9] = 0.00125395 and 14	
lambda[10] = 0.000805858 and 14	
lambda[11] = 0.000797636 and 14	
lambda[12] = 0.000502732 and 14	
lambda[13] = 0.000359993 and 14	
marchywka@happy:/home/documents/cpp/proj/freefem/play/junk/Problem_EV`

I went through man solvers same result. Your python coefficients are not the same-
dod you look at the condition number or do SVD in either case? There may be
a sensivitivty problem. The eigenvalues are less than .1 and you have coefficients
as high as 13. 

[https://slepc.upv.es/documentation/current/docs/manualpages/EPS/EPSType.html](https://slepc.upv.es/documentation/current/docs/manualpages/EPS/EPSType.html)
```

---

<div class="post-metadata">

**Author:** ![Loic](https://avatars.discourse-cdn.com/v4/letter/l/e9c0ed/32.png) [@Loic](https://community.freefem.org/u/Loic)\
**Post date:** [April 26, 2023, 3:01pm UTC](https://community.freefem.org/t/problem-for-computing-eigenvalues-with-epssolve/2420/11 "2023-04-26T15:01:58Z")

</div>

Yes, indeed, these are the values I find with Python.

And we can see that 12 of them are similar with those obtained with PETSc.

Loïc,

---

<div class="post-metadata">

**Author:** ![Loic](https://avatars.discourse-cdn.com/v4/letter/l/e9c0ed/32.png) [@Loic](https://community.freefem.org/u/Loic)\
**Post date:** [April 26, 2023, 3:22pm UTC](https://community.freefem.org/t/problem-for-computing-eigenvalues-with-epssolve/2420/12 "2023-04-26T15:22:02Z")

</div>

Ah okay, but theoretically it should not be.

Since, my operators are `MASS = \int_K \Phi_i \Phi_j` and `STIFF= \int_K grad(\Phi_i) grad(\Phi_j)`. The `\Phi_i` define a Finite Element Basis.

Is it possible that this is due to approximations errors ?

Normally, I have the same matrices in Python.

Loïc,

---

<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:** [April 26, 2023, 3:32pm UTC](https://community.freefem.org/t/problem-for-computing-eigenvalues-with-epssolve/2420/13 "2023-04-26T15:32:05Z")

</div>

`STIFF= \int_K grad(\Phi_i) grad(\Phi_j)` this matrix is clearly indefinite. Everything is defined up to a constant.

---

<div class="post-metadata">

**Author:** ![Loic](https://avatars.discourse-cdn.com/v4/letter/l/e9c0ed/32.png) [@Loic](https://community.freefem.org/u/Loic)\
**Post date:** [April 26, 2023, 4:00pm UTC](https://community.freefem.org/t/problem-for-computing-eigenvalues-with-epssolve/2420/14 "2023-04-26T16:00:27Z")

</div>

You mean numerically in my case ? (Since theoretically, it should be)

Loïc,

---

<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:** [April 26, 2023, 4:04pm UTC](https://community.freefem.org/t/problem-for-computing-eigenvalues-with-epssolve/2420/15 "2023-04-26T16:04:26Z")

</div>

I mean that your right-hand side matrix is indefinite, both theoretically and numerically.

---

<div class="post-metadata">

**Author:** ![Loic](https://avatars.discourse-cdn.com/v4/letter/l/e9c0ed/32.png) [@Loic](https://community.freefem.org/u/Loic)\
**Post date:** [April 26, 2023, 4:17pm UTC](https://community.freefem.org/t/problem-for-computing-eigenvalues-with-epssolve/2420/16 "2023-04-26T16:17:44Z")

</div>

To my mind, if I am not mistaken, the stiffness matrix for a Poisson problem for example is symmetric and positive definite.

Loïc,

---

<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:** [April 26, 2023, 4:18pm UTC](https://community.freefem.org/t/problem-for-computing-eigenvalues-with-epssolve/2420/17 "2023-04-26T16:18:32Z")

</div>

With no boundary conditions, of course it is not, again, everything is defined up to a constant.

---

<div class="post-metadata">

**Author:** ![Loic](https://avatars.discourse-cdn.com/v4/letter/l/e9c0ed/32.png) [@Loic](https://community.freefem.org/u/Loic)\
**Post date:** [April 26, 2023, 4:22pm UTC](https://community.freefem.org/t/problem-for-computing-eigenvalues-with-epssolve/2420/18 "2023-04-26T16:22:20Z")

</div>

Ah okay, yes you are right,  
Sorry for this confusion,

Loïc,

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [April 26, 2023, 8:40pm UTC](https://community.freefem.org/t/problem-for-computing-eigenvalues-with-epssolve/2420/19 "2023-04-26T20:40:07Z")

</div>

What is row 9 in the stiffness matrix? The diagonal element is large  
, IIRC 6600, and some of the other terms look like they have similar coefs of opposite  
signs. I guess if there is a delicate cancellation that could cause  
numerical issues. It probably is of no significance but if you zero out the  
terms less than a noise floor, like 1e-6 or so, what do you get?

Since all the solvers I tried gave more or less the same result, it  
does not seem to be a quirk of any one of them with the original  
FF coefficients

---

<div class="post-metadata">

**Author:** ![Loic](https://avatars.discourse-cdn.com/v4/letter/l/e9c0ed/32.png) [@Loic](https://community.freefem.org/u/Loic)\
**Post date:** [April 27, 2023, 5:47am UTC](https://community.freefem.org/t/problem-for-computing-eigenvalues-with-epssolve/2420/20 "2023-04-27T05:47:37Z")

</div>

The Stiffness matrix is just a correlation matrix in norm H1. So my problem doesn’t really make physical sense. It only says that some basis functions have more “energy” than others.

Loïc,

---

<div class="post-metadata">

**Author:** ![Loic](https://avatars.discourse-cdn.com/v4/letter/l/e9c0ed/32.png) [@Loic](https://community.freefem.org/u/Loic)\
**Post date:** [April 27, 2023, 7:02am UTC](https://community.freefem.org/t/problem-for-computing-eigenvalues-with-epssolve/2420/21 "2023-04-27T07:02:19Z")

</div>

HI @prj,

I asked around, and the stiffness matrix is at least positive semi-definite (all the eigenvalues should be greater or equal to zero). The fact that everything is defined up to a constant should not lead to negative eigenvalue, but we can have eigenvalues equal to zero which corresponds to rigid body motion.

The Mass matrix is for sure definite positive (all eigenvalues strictly greater than zero).

Consequently, all eigenvalues should be positive.

By the way when I consider the problem in the other way : STIFF \Phi = \lambda MASS \Phi, I find positive eigenvalues, and the same values than with Python.

So, we think that may be, the problem is the way we compute the eigenvalues. Indeed, when computing generalized eigenvalue problem, there is often a factorization that is is done. May be the fact that some eigenvalues of the stiffness matrix is equal to zero (when I compute STIFF \Phi = \lambda \Phi, I find two eigenvalues equal to zero) and some approximations can cause a problem.

What do you think of this ?

Loïc,

[Next page](https://community.freefem.org/t/problem-for-computing-eigenvalues-with-epssolve/2420.md?page=2)