matrix factorization reuse?

[Zhen Tao]

Dear users, 

I am considering a time dependent problem, say, step-26 with mesh adaptive refinement shut down.

Then it turns out that the solve_time_step() is solving exactly the same matrix every time step,

the only thing change is the right hand side.

For problems like this, is it a waste to call a sparse solver (cg in step-26) everytime? pros and cons?

The way I am thinking is to make a LAPACKFullMatrix, call the super expensive inverse() at the beginning,

store it, and then refresh the rhs every time, only do the matrix-vector-multiplication afterwards.

Or set the inverse matrix (a lapack full matrix) as a preconditioner, then call a 1-step pseudo cg every time. (is that possible?)

Therefore my following questions are:

What is the limit of the size of matrix system can I do a lapack LU factorization on a single machine?

Are there any better way to get LU factorization in deal.ii?

Can we incorporate "history experience" to achieve more efficient solvers?

Thanks,

- Zhen Tao

[Jean-Paul Pelteret]

I might be wrong, but it appears to me that SparseDirectUMFPACK factorises the matrix once up front and retains this factorisation for each solve(rhs_and_solution) until you next call factorize(matrix) or initialize(matrix). Might this be what you're looking for?

J-P

[Wolfgang Bangerth]

Zhen,

Jean-Paul already gave the correct answer, namely use a sparse direct solver
and factorize only once.

> Therefore my following questions are:
>
> What is the limit of the size of matrix system can I do a lapack LU
> factorization on a single machine?

For a sparse direct solver, I think you can probably go to around 1M unknowns.
If you use a LAPACK matrix, you will quickly be limited by memory. For
example, a 10,000 x 10,000 matrix already requires almost a GB of memory.

> Are there any better way to get LU factorization in deal.ii?
> Can we incorporate "history experience" to achieve more efficient solvers?

People have tried that and there is a substantial amount of literature on the
topic of Krylov subspace recycling. I always think that it is a great idea and
want to investigate using it, but it has not found its way into mainstream
applications to the best of my knowledge.

Best
W.

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

[Zhen Tao]

Thanks @Jean-Paul & @Wolfgang ! 

SparseDirectUMFPACK is exactly what I need.

[Pascal Kraft]

Hello Wolfgang,

does the same hold true for other direct solvers? I would like to use Mumps (PETSc) to get a LU-factorization and reuse that a couple of times. The PETSc wrappers only provide a constructor and solve. If I reuse solve withe the same matrix argument, will the factorization be done again or does it detect that the inverse is already present? The code for both the PETSc wrappers as well as the Trilinos Wrappers does not look as if that was supported.

Greetings,
Pascal Kraft

[Michael Harmon]

If you want to use Trilinos, I refactored their direct solver so that it can be used the way you want, I will submit the patch, but if you want can send you the changes as well.

[David F]

Hi, I'm interested in this as well. Could you post or link to your solution? Thanks.

[Wolfgang Bangerth]

On 08/30/2016 11:09 AM, David F wrote:
> Hi, I'm interested in this as well. Could you post or link to your solution?

David -- this went into the code base with this patch:
https://github.com/dealii/dealii/pull/2493

Best
Wolfgang

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