Eigenvalue problem with a mesh-independent input space

General Discussion
1

Hello everyone,

I would like to solve an eigenvalue problem in the form of sigma^2 B x = A x using PETSc The matrices A and B are symmetric and positive definite. The catch is that I would like the vector ‘x’ be independent of the mesh size, lets say a size (3,1) vector. I can specify the action of the matrices B and A using the ‘‘‘func’’’ feature of FreeFEM; but to do that I need a PETSc Mat of the same size. The issue is that I want the (3,1) input vector to be the same on all processes, but naturally if I define a B matrix with size (3,1) each process will have its separate unknowns. I think defining a (3,1) vector on rank 0 would resolve the issue but I cannot seem to have zero size matrices on the other processes (I cannot seem to create a PETSc MAT from a FreeFEM Mat that way). Any suggestions are welcome how to resolve this issue!

2

Your notation are not consistent, what is a size (3,1) vector and how can you have both this and a B matrix with size (3,1)? Anyway, here is how to define a Mat of size m \times n, exclusively stored on rank p:

load "PETSc"

int p = 0;
int m = 4;
int n = 3;
Mat A(mpirank == p ? m : 0, mpirank == p ? n : 0);
3

Thanks for the prompt reply, I really appreciate it. The method you suggested works well to define the matrices, but there is a following step in which I am stuck.

I want to use EPSSolve so that the eigenproblem is defined for a few unknowns, like 3, but the application of the operator requires solution of large linear systems of equations. Thus, I think the shell capability FreeFEM/PETSc should be used. However, I cannot manage to make it work the way I want.

The MWE below illustrates the behavior I do not understand. I would imagine that the function applyM is being called to apply multiplication with the matrix M; however, that does not seem to be the case. What am I missing?

load "PETSc"

include "getARGV.idp"

/* matrix creation */
int n = 3;
matrix MinFF(n,n);
for(int i=0; i<n; i++){
    for(int j=0; j<n; j++){
        MinFF(i,j) = 10*i + j;
    }
}

Mat Min(n,n);
Min = MinFF;

/* shell function for successive operations */
int lab;
func real[int] applyM(real[int]& inPETSc) {

    cout << "Rank " << mpirank << "/" << mpisize << ", M appliead" << endl;
    real[int] outPETSc(n);
    MatMult(Min,inPETSc,outPETSc);
    return outPETSc;

}
Mat Mmat(Min, applyM);
ObjectView(Min);
ObjectView(Mmat);

/* testing the B matrix application - seems OK */
real[int] testOut(n);
real[int] testIn(n);
testIn = [1,2,3];

real evtol = getARGV("-evtol", 1e-6);
int nev = getARGV("-nev", 2);
int ncv = getARGV("-ncv", 6);
int nIter = getARGV("-nIter", 10);
string params = " -eps_non_hermitian -eps_tol " + evtol + " -eps_nev " + nev + " -eps_ncv " + ncv + " -eps_max_it " + nIter + 
      " -eps_largest_magnitude -eps_type krylovschur -eps_monitor_conv -eps_monitor_all -eps_converged_reason";

real[int] eval(nev);      // array to store eigenvalues
real[int] evalErrors(nev);      // array to store eigenvalues
real[int,int] eigenvec(n,nev); // array to store singular values

int k = EPSSolve(Mmat, array = eigenvec, values  = eval, sparams = params, errorestimate = evalErrors);
4

I need to look into this some more, but as a temporary workaround, you can replace Mat Mmat(Min, applyM); by Mat Mmat(Min, Min, applyM); and your functor will be called.

5

Fixed in Fix wrongful call to MatDuplicate() · FreeFem/FreeFem-sources@1d0be0f · GitHub, thanks for the report.

1 Like
6

Thanks for the quick fix. I compiled FreeFEM with the fix, and the MWE works indeed!