Operations between sparse matrices with incompatible IndexSets

[giovann...@hotmail.it]

Hi everyone,

   the coupling matrix (non_matching::coupling_matrix) obtained with parallel meshes is (practically) guaranteed to have and IndexSet which is not compatible with the original ones, as the matrix is obtained multiplying base functions from two indipendent different spaces, while all other sensible matrices (e.g. mass matrix) are computed using a functions from a single space. This is a problem as the basic matrix operations ( + and * ) are not possible (at least using TrilinosWarppers::linear_operator).

As a solution, I was thinking about implementing a function which, given a sparse matrix with incompatible indices and the "wished" IndexSet, makes a global communication to generate a "compatible" matrix. Is this reasonable? Are the better approaches? Or maybe Trilinos can already handle this situation and it is only a matter of implementation in deal.II?

Thank you.

Best,

Giovanni

[Bruno Turcksin]

Giovanni,

On Thursday, January 3, 2019 at 9:11:30 AM UTC-5, giovann...@hotmail.it wrote:

As a solution, I was thinking about implementing a function which, given a sparse matrix with incompatible indices and the "wished" IndexSet, makes a global communication to generate a "compatible" matrix. Is this reasonable? Are the better approaches? Or maybe Trilinos can already handle this situation and it is only a matter of implementation in deal.II?

You could try this function (there is a similar one to multiply two matrices). I am not sure if it will work since they don't explicit say how the matrices can or cannot be distributed.

Best,

Bruno

[giovann...@hotmail.it]

Thank you for your quick, helpful, answer Bruno!!

Looking at the documentation of EpetraExt::MatrixMatrix::Multiply which you linked, it's written:

   In a parallel setting, A and B need not have matching distributions, but C needs to have the same row-map as A.

So I looked into it with great hope; currently I am using TrilinosWrappers::SparseMatrix function, which I just discovered

is a wrapper for Epetra_FECrsMatrix ( which is then made into a linear operator),

and inherits from Epetra_CrsMatrix.

I am actually transforming matrices into linear operators, before running the operations. Thus

I tried to run directly the matrix product using "mmult" (which is supposed

to call EpetraExt::MatrixMatrix::Multiply ). Unfortunately this results in an error:

    Parallel partitioning of A and B does not fit.

which is triggered here:

          Assert(inputleft.n() == inputright.m(),
                 ExcDimensionMismatch(inputleft.n(), inputright.m()));
          Assert(inputleft.domain_partitioner().SameAs(
                   inputright.range_partitioner()),
                 ExcMessage("Parallel partitioning of A and B does not fit."));

I now wonder if it would be wise to remove the second check (or substitute for one where inputleft and result are checked)...and then, what about the sum?

[Bruno Turcksin]

Giovanni,

Le jeu. 3 janv. 2019 à 14:10, <giovann...@hotmail.it> a écrit :
> which is triggered here:
>
>           Assert(inputleft.n() == inputright.m(),
>                  ExcDimensionMismatch(inputleft.n(), inputright.m()));
>           Assert(inputleft.domain_partitioner().SameAs(
>                    inputright.range_partitioner()),
>                  ExcMessage("Parallel partitioning of A and B does not fit."));
>

> I now wonder if it would be wise to remove the second check (or substitute for one where inputleft and result are checked)...and then, what about the sum?

I think that the check can be removed. Until recently we were not using Epetra_Ext to do the multiplication so it's possible that the old version did not work in this case.  BTW, you need to be on recent version of deal.II for this function to work correctly all the time. I found some bugs in it and that's why we switched to Epetra_Ext. For the sum, the easiest is to extract the Epetra matrices using trilinos_matrix(), then use the Epetra_Ext function yourself, and finally create a TrilinosWrapper matrix using reinit()

Best,

Bruno

[giovann...@hotmail.it]

The version I am using is the developeing one, maybe a couple of days old, so it should work...Thank you very much for your help! I will now try to follow your suggestions!