Vector eigenvalue problem with SLEPc

[Anton Ermakov]

Dear Deal.II community,

I am working on implementing a vector eigenproblem in deal.ii, starting from the step-36 example. The question I am studying is the normal oscillation of a fluid planet. The eigenproblem is of this type:

 _                                                                      _

 I  I¯  A  0  ¯I    *eigenvalue     +  I¯    0 B ¯ I I  * eigenvector =0

 I  I_   0 0   _I                                  I_   C D   _I I

 ¯                                                                     ¯

where eigenvector consists of a vector displacement and a scalar pressure perturbation: 

eigenvector = transpose([U P])

It seems from the literature (e.g., Wang & Bathe 1997) that it is standard to rewrite this system to eliminate the pressure variable:

eigenvalue * A *  U - B * D^-1 * C * U = 0

or 

eigenvalue * A *  U + M * U = 0

with M = - B * D^-1 * C and P = -D^-1 * C * U 

What I am having trouble with is how to efficiently represent matrix products and inverses that go into matrix M. All the matrices (A, B, C, D) are sparse, but due to the inverse, it seems that M will be a full matrix and it might be unaffordable to store it in the memory. There was a discussion on a similar topic here:

https://groups.google.com/g/dealii/c/5P_fv0zg7jg/m/CPYcYrbsDQAJ

It seems that LinearOperator was used in that discussion to manipulate the matrix products and inverses, which seemed like an approach to follow if you use ARPACK eigensolvers. However, I am using SLEPc and it seems the inputs to SLEPc solvers would have to be objects of class PETScWrappers::MatrixBase. 

So, my question is how to prepare such a matrix M involving matrix products and matrix inverse efficiently and represent it as PETScWrappers::MatrixBase suitable for a SLEPc eigensolver. 

In addition, I am interested in the interior eigenvalues. It seems that shift-invert spectral transformation is available in SLEPc - this is what I plan to use. However, SLEPc also has a polynomial filtering spectral transformation. It does not seem that deal.ii provides a wrapper to that one though. So, I wonder if anyone has had expertise with making an interface in deal.ii to use SLEPc's polynomial filtering transformation.

Thanks a lot,

Anton.

[Wolfgang Bangerth]

On 8/25/21 11:53 AM, Anton Ermakov wrote:
>  I  I¯  A  0  ¯I    *eigenvalue     +  I¯    0 B ¯ I I  * eigenvector =0
>  I  I_   0 0   _I                                  I_   C D   _I I
>  ¯                                                                     ¯
> where eigenvector consists of a vector displacement and a scalar
> pressure perturbation:
>
> eigenvector = transpose([U P])
>
> It seems from the literature (e.g., Wang & Bathe 1997) that it is
> standard to rewrite this system to eliminate the pressure variable:
>
> eigenvalue * A *  U - B * D^-1 * C * U = 0
>
> or
>
> eigenvalue * A *  U + M * U = 0
>
> with M = - B * D^-1 * C and P = -D^-1 * C * U

You can do that, but I believe that it's also possible to work with the
original system where one of the two matrices is only semidefinite.

> What I am having trouble with is how to efficiently represent matrix
> products and inverses that go into matrix M. All the matrices (A, B, C,
> D) are sparse, but due to the inverse, it seems that M will be a full
> matrix and it might be unaffordable to store it in the memory. There was
> a discussion on a similar topic here:
>
> https://groups.google.com/g/dealii/c/5P_fv0zg7jg/m/CPYcYrbsDQAJ
>
> It seems that LinearOperator was used in that discussion to manipulate
> the matrix products and inverses, which seemed like an approach to
> follow if you use ARPACK eigensolvers. However, I am using SLEPc and it
> seems the inputs to SLEPc solvers would have to be objects of class
> PETScWrappers::MatrixBase

> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dealii.org%2Fcurrent%2Fdoxygen%2Fdeal.II%2FclassPETScWrappers_1_1MatrixBase.html&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Ce376b299e56c44690ada08d967f13673%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637655108014800275%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=QQLBoH27O1sZ6GMXDatmucrxAye8edSeN3XNqAybONI%3D&reserved=0>.

>
>
> So, my question is how to prepare such a matrix M involving matrix
> products and matrix inverse efficiently and represent it as
> PETScWrappers::MatrixBase

That might not easily be possible with our PETSc interfaces.

PETSc has the concept of MatShell, which is in essence what deal.II's
LinearOperator is: A class that implements a matrix-vector operation,
without giving access to individual matrix elements. You could try
whether you can implement something derived from MatrixBase that is a
MatShell operation in the same way as the existing classes derived from
MatrixBase model other types of PETSc matrices.

I don't know how difficult that would be, though.

Best
W.


--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Anton Ermakov]

>>> That might not easily be possible with our PETSc interfaces.

PETSc has the concept of MatShell, which is in essence what deal.II's
LinearOperator is: A class that implements a matrix-vector operation,
without giving access to individual matrix elements. You could try
whether you can implement something derived from MatrixBase that is a
MatShell operation in the same way as the existing classes derived from
MatrixBase model other types of PETSc matrices.

I don't know how difficult that would be, though.

Ok. I see. Seems like using ARPACK might be easier then.

>>> You can do that, but I believe that it's also possible to work with the

original system where one of the two matrices is only semidefinite.

I have redone it with the original system, but it seems that PETSc is complaining about missing diagonal entry. 

[0]PETSC ERROR: Object is in wrong state

[0]PETSC ERROR: Matrix is missing diagonal entry 3

I presume this happens due to a zero block in one (or both) of the original matrices. Hmmm, should I then just explicitly see all zero diagonal values to 0.0.

Anton.

[Wolfgang Bangerth]

On 8/25/21 9:37 PM, Anton Ermakov wrote:
>
> >>> You can do that, but I believe that it's also possible to work with the
> original system where one of the two matrices is only semidefinite.
>
> I have redone it with the original system, but it seems that PETSc is
> complaining about missing diagonal entry.
>

> /[0]PETSC ERROR: Object is in wrong state/
>
> /[0]PETSC ERROR: Matrix is missing diagonal entry 3/

>
> I presume this happens due to a zero block in one (or both) of the
> original matrices. Hmmm, should I then just explicitly see all zero
> diagonal values to 0.0.

I don't know whether the error complains about the value missing in the
sparsity pattern, or it being zero.

In any case, you might have to read up in the documentation of SLEPC
about two issues:

* You can solve Ax = lambda M x or M x = mu A x, where mu=1/lambda. The
problem is well posed if one of the two matrices is not singular, but I
bet that SLEPc has a convention where it either wants the matrix on the
left, or the one on the right to be the non-singular one.

* There is probably a preconditioner involved somewhere. If that
preconditioner happens to be something like Jacobi, then it will need to
divide by a diagonal entry. You will need to make sure that the matrix
that the preconditioner is applied to has nonzero diagonals -- or,
alternatively, choose a different preconditioner.