mmult memory leak with petsc

[Richard Schussnig]

Hi everyone,

I am working on incompressible flow problems and stumbled upon an issue when calling PETScWrappers::SparseMatrix::mmult(),
but before I describe the problem in more detail, let me comment on the basic building blocks of the MWE:

(i) parallel::distributed::Triangulation & either PETSc or Trilinos linear algebra packages

(ii) dim-dimensional vector-valued FE space for velocity components & scalar-valued FE space for pressure,

simply constructed via:

FESystem<dim> fe (FE_Q<dim>(vel_degree), dim,
                                 FE_Q<dim>(press_degree), 1);

So, after integrating the weak form -- or just filling the matrices with some entries -- we end up with a block system

A u + B p = f

C u + D p = g.

To construct some preconditioners, we have to perform some matrix-matrix products:

either for the Schur complement

(a) S = D - C inv(diag(A)) B

or some A_gamma

(b) A_gamma = A + gamma * B inv(diag(Mp)) C.

Comletely ignoring now, why that might be necessary or not
(I know that there is the possibility of assembling a grad-div term and using a
Laplacian on the pressure space to account for the reaction term,
but that is not really an option in the stabilized case),

we need those matrix products, and here comes the problem:

using either PETSc or Trilinos I get identical matrix-products when calling mmult(), BUT
when using PETSc, the RAM is slowly but steadily filled (up to 500GB on our local cluster)

I came up with the MWE attached, which does nothing else than initializing the system
and then constructs the matrix product 1000 times in a row.

Am I doing anything wrong or is this supposed to be used differently?
I am using dealii v.9.0.1 installed via candi, so maybe the old version is the reason.

Any suggestions? Any help would be greatly appreciated!

Bonus question:
Is there a way similar to hand the sparsity patterns over to the mmult function?
(For the dealii::SparseMatrix there is, which is why I am asking)

Kind regards,
Richard

[Wolfgang Bangerth]

Richard,

Nice! Can I ask you to play with this some more? I think you can make that
code even more minimal:
* Remove all of the commented out stuff -- for the purposes of reproducing the
problem it shouldn't matter.
* Move the matrix initialization code out of the main loop. You want to show
that it's the mmult that is the problem, but you're having a lot of other code
in there as well that could in principle be the issue. If you move the
initialization of the individual factors out of the loop and only leave
whatever is absolutely necessary for the mmult in the loop, then you've
reduced the number of places where one needs to look.
* I bet you could trim down the list of #includes by a good bit :-)

You seem to be using a pretty old version of deal.II. There are a number of
header files that no longer exist, and some code that doesn't compile for me.
For your reference, attached is a version that compiles on the current master
branch (though with a number of warnings). That said, it seems like the memory
doesn't explode for me -- which raises the question of which version of
deal.II and PETSc you use. For me, this is deal.II dev and PETSc 3.7.5.

> Am I doing anything wrong or is this supposed to be used differently?
> I am using dealii v.9.0.1 installed via candi, so maybe the old version is the
> reason.

Possible -- no need to chase an old bug that has already been fixed if you can
simply upgrade.

> Bonus question:
> Is there a way similar to hand the sparsity patterns over to the mmult function?
> (For the dealii::SparseMatrix there is, which is why I am asking)

DynamicSparsityPattern has compute_mmult_pattern, which should give you a
sparsity pattern you can then use to initialize the resulting PETSc matrix.

Best
W.

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

[Richard Schussnig]

Hi Wolfgang,

thanks a ton for taking a look at my silly example!
- It seems like I am now finally really at a point where an update is inevitable! ; )
... maybe just rewriting my codes to fit a newer version would have saved me loads of time in the long run,
but I am happy that remains a mistery.

So, I will update and comment here afterwards;
also, I will work on the mwe as you suggested in case the upgrade does not get rid of the problem.

Regards,
Richard

[Richard Schussnig]

Hi again,

I took a look at the modified mwe you sent -- immediately noticing that you read though the whole thing,
which is amazing, considering that this is a google group! 

So, thanks again for your time and effort, it is really appreciated!

I have been busy doing the following with our little toy problem:
(i)  Upgraded to deal.ii v9.2.0 with candi, together with PETSc 3.13.1 (and Trilinos 12.18.1).

(ii) Further stripped down the mwe and deleted unnecessary parts and moved all but the mmult() out of the 'main loop'.

So, I still observe the following:

using Trilinos, all is fine.

and using PETSc, using the call

leftMatrix.mmult(result, rightMatrix)

seems fine, whereas calling

leftMatrix.mmult(result, rightMatrix, vector)

still results in a memory leak for me.

Could you verify that again for me please?

If this is doing fine for you I will just multiply the rows of rightMatrix 'by hand'
before using the call without the vector.

Thanks in advance & kind regards,

Richard

[Wolfgang Bangerth]

> I took a look at the modified mwe you sent -- immediately noticing that you
> read though the whole thing,
> which is amazing, considering that this is a google group!
> So, thanks again for your time and effort, it is really appreciated!

:-)

> I have been busy doing the following with our little toy problem:
> (i)  Upgraded to deal.ii v9.2.0 with candi, together with PETSc 3.13.1 (and
> Trilinos 12.18.1).
> (ii) Further stripped down the mwe and deleted unnecessary parts and moved all
> but the mmult() out of the 'main loop'.
>
> So, I still observe the following:
> using Trilinos, all is fine.
>
> and using PETSc, using the call
> leftMatrix.mmult(result, rightMatrix)
> seems fine, whereas calling
> leftMatrix.mmult(result, rightMatrix, vector)
> still results in a memory leak for me.
>
> Could you verify that again for me please?

I'm running this program and watching memory consumption in 'top', looking at
the 'virt' and 'res' columns. It's growing very slowly: after 3.5 minutes, I'm
now at about 2GB and 1GB, respectively. Is this about as much as you get as well?

Looking at the code in question, we end up in this function:
https://github.com/dealii/dealii/blob/master/source/lac/petsc_matrix_base.cc#L494-L561

Interestingly, I believe that your case specifically goes into line 539
(please confirm), which then calls MatDuplicate. But we never delete that
object again -- one could suspect that there needs to be a MatDestroy at the
end of that block. Would you be interested in trying that out?

[Richard Schussnig]

Hi Wolfgang,

I do have approx. the same numbers, but on 4 cores (see screenshot).

Fortunately, I was able to adapt our little program here with success
- which I append here just for the sake of completeness, since I went on github and tried my best to propose that change, see

https://github.com/dealii/dealii/pull/10838

I hope that worked out, but I used git for the first time.

Nonetheless, some concerns of mine remain:

(i) if I allow in my application-code to actually change the row scaling of the right matrix in the product (assuming the vector has all nonzero entries),
I will actually get away with only 1 additional matrix being the product instead of 2 including the product and a temporary one (which was the one not destroyed)
by the cost of scaling and rescaling the right matrix in the product.

(ii) The documentation says, that both Trilinos and PETSc actually rebuilding the sparsity pattern of the product,
so using compute_mmult_pattern in the setup of the product would be of no use, right?

Thanks again & kind regards,
Richard

[Wolfgang Bangerth]

Hi Richard,

> I do have approx. the same numbers, but on 4 cores (see screenshot).
>
> Fortunately, I was able to adapt our little program here with success
> - which I append here just for the sake of completeness, since I went on
> github and tried my best to propose that change, see
> https://github.com/dealii/dealii/pull/10838

> <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fdealii%2Fdealii%2Fpull%2F10838&data=02%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cda1531178e7b4705abac08d843a89035%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637333738540558308&sdata=1VdEvclgxPPonuEzf1nrSpc0TGmge9QTTXlkJ1snjPk%3D&reserved=0>

> I hope that worked out, but I used git for the first time.

It's the first time at some time for everyone -- but we hope it won't be your
last patch!

> Nonetheless, some concerns of mine remain:
> (i) if I allow in my application-code to actually change the row scaling of
> the right matrix in the product (assuming the vector has all nonzero entries),
> I will actually get away with only 1 additional matrix being the product
> instead of 2 including the product and a temporary one (which was the one not
> destroyed)
> by the cost of scaling and rescaling the right matrix in the product.

Right. You might not actually need this function, but it's still worth fixing :-)

> (ii) The documentation says, that both Trilinos and PETSc actually rebuilding
> the sparsity pattern of the product,
> so using compute_mmult_pattern in the setup of the product would be of no use,
> right?

For the moment that's correct. But the PETSc documentation provides a way how
one could re-use a previously computed sparsity pattern. One could add a flag
to the mmult function that switches between either telling PETSc to create a
new sparsity pattern, or re-using the existing one.

If you only call the function once, I have no doubt that there is not much of
a difference between using the deal.II and PETSc ways of creating sparsity
patterns. It's a different story if you plan on calling this function multiple
times using the same sparsity pattern.

[richard....@gmail.com]

Hi Wolfgang,

thank you for spending so much time on that issue; 
I am happy I could at least contribute a tiny bit!

If I happen to implement that functionality in the future, I will simply put in on github ; )

Kind regards,
Richard