# how to inverse a sparse matrix in Petsc?

Barry Smith bsmith at mcs.anl.gov
Wed Feb 6 07:52:39 CST 2008

```    Whoops, actually the grids in each direction would need to be
uniform.

Barry

On Feb 5, 2008, at 9:39 PM, Ben Tay wrote:

> Sorry Barry, I just would like to confirm that as long as it's a
> constant constant coefficient Poisson eqn with Neumann or Dirchelet
> boundary conditions, I can use FFT. It doesn't matter if the grids
> are uniform or not. Is that correct? Thanks.
>
> Barry Smith wrote:
>>
>> On Feb 5, 2008, at 8:48 PM, Ben Tay wrote:
>>
>>> Thank you Barry for your enlightenment. I'll just continue to use
>>> BoomerAMG for the poisson eqn. I'll also check up on FFTW. Last
>>> time, I recalled that there seemed to be some restrictions for FFT
>>> on solving poisson eqn. It seems that the grids must be constant
>>> in at least 1 dimension.
>>
>>   Yes. Then it decouples into a bunch of tridiagonal solves;
>> Basically if you can do separation of variables you can
>> use FFTs.
>>
>>   Barry
>>
>>> I wonder if that is true? If that's the case, then it's not
>>> possible for me to use it, although it's a constant coefficient
>>> Poisson operator with Neumann or Dirchelet boundary conditions.
>>>
>>> thank you.
>>>
>>> Barry Smith wrote:
>>>>
>>>> On Feb 5, 2008, at 8:04 PM, Ben Tay wrote:
>>>>
>>>>> Hi Lisandro,
>>>>>
>>>>> I'm using the fractional step mtd to solve the NS eqns as well.
>>>>> I've tried the direct mtd and also boomerAMG in solving the
>>>>> poisson eqn. Experience shows that for smaller matrix, direct
>>>>> mtd is slightly faster but if the matrix increases in size,
>>>>> boomerAMG  is faster. Btw, if I'm not wrong, the default solver
>>>>> will be GMRES. I've also tried using the "Struct" interface
>>>>> solely under Hypre. It's even faster for big matrix, although
>>>>> the improvement doesn't seem to be a lot. I need to do more
>>>>> tests to confirm though.
>>>>>
>>>>> I'm now doing 2D simulation with 1400x2000 grids. It's takes
>>>>> quite a while to solve the eqns. I'm wondering if it'll be
>>>>> faster if I get the inverse and then do matrix multiplication.
>>>>> Or just calling KSPSolve is actually doing something similar and
>>>>> there'll not be any speed difference. Hope someone can
>>>>> enlighten...
>>>>>
>>>>> Thanks!
>>>>>
>>>>  Ben,
>>>>
>>>>    Forming the inverse explicitly will be a complete failure.
>>>> Because it is dense it will have (1400x2000)^2 values and
>>>> each multiply will take 2*(1400x2000)^2 floating point
>>>> operations, while boomerAMG should take only O(1400x2000).
>>>>
>>>>    BTW: if this is a constant coefficient Poisson operator with
>>>> Neumann or Dirchelet boundary conditions then
>>>> likely a parallel FFT based algorithm would be fastest. Alas we
>>>> do not yet have this in PETSc. It looks like FFTW finally
>>>> has an updated MPI version so we need to do the PETSc interface
>>>> for that.
>>>>
>>>>
>>>>  Barry
>>>>
>>>>
>>>>> Lisandro Dalcin wrote:
>>>>>> Ben, some time ago I was doing some testing with PETSc for
>>>>>> solving
>>>>>> incompressible NS eqs with fractional step method. I've found
>>>>>> that in
>>>>>> our software and hardware setup, the best way to solve the
>>>>>> pressure
>>>>>> problem was by using HYPRE BoomerAMG. This preconditioner
>>>>>> usually have
>>>>>> some heavy setup, but if your Poison matrix does not change,
>>>>>> then the
>>>>>> sucessive solves at each time step are really fast.
>>>>>>
>>>>>> If you still want to use a direct method, you should use the
>>>>>> combination '-ksp_type preonly -pc_type lu' (by default, this
>>>>>> will
>>>>>> only work on sequential mode, unless you build PETSc with an
>>>>>> external
>>>>>> package like MUMPS). This way, PETSc computes the LU
>>>>>> factorization
>>>>>> only once, and at each time step, the call to KSPSolve end-up
>>>>>> only
>>>>>> doing the triangular solvers.
>>>>>>
>>>>>> The nice thing about PETSc is that, if you next realize the
>>>>>> factorization take a long time (as it usually take in big
>>>>>> problems),
>>>>>> you can switch BoomerAMG by only passing in the command line
>>>>>> '-ksp_type cg -pc_type hypre -pc_hypre_type boomeramg'. And
>>>>>> that's
>>>>>> all, you do not need to change your code. And more, depending
>>>>>> on your
>>>>>> problem you can choose the direct solvers or algebraic
>>>>>> multigrid as
>>>>>> you want, by simply pass the appropriate combination options in
>>>>>> the
>>>>>> command line (or a options file, using the -options_file option).
>>>>>>
>>>>>> Please, if you ever try HYPRE BoomerAMG preconditioners, I
>>>>>> would like
>>>>>>
>>>>>> Regards,
>>>>>>
>>>>>> On 2/5/08, Ben Tay <zonexo at gmail.com> wrote:
>>>>>>
>>>>>>> Hi everyone,
>>>>>>>
>>>>>>> I was reading about the topic abt inversing a sparse matrix. I
>>>>>>> have to
>>>>>>> solve a poisson eqn for my CFD code. Usually, I form a system
>>>>>>> of linear
>>>>>>> eqns and solve Ax=b. The "A" is always the same and only the
>>>>>>> "b" changes
>>>>>>> every timestep. Does it mean that if I'm able to get the
>>>>>>> inverse matrix
>>>>>>> A^(-1), in order to get x at every timestep, I only need to do
>>>>>>> a simple
>>>>>>> matrix multiplication ie x=A^(-1)*b ?
>>>>>>>
>>>>>>> Hi Timothy, if the above is true, can you email me your
>>>>>>> Fortran code
>>>>>>> template? I'm also programming in fortran 90. Thank you very
>>>>>>> much
>>>>>>>
>>>>>>> Regards.
>>>>>>>
>>>>>>> Timothy Stitt wrote:
>>>>>>>
>>>>>>>> Yes Yujie, I was able to put together a parallel code to
>>>>>>>> invert a
>>>>>>>> large sparse matrix with the help of the PETSc developers. If
>>>>>>>> you need
>>>>>>>> any help or maybe a Fortran code template just let me know.
>>>>>>>>
>>>>>>>> Best,
>>>>>>>>
>>>>>>>> Tim.
>>>>>>>>
>>>>>>>>
>>>>>>>>> Hi
>>>>>>>>> There was a discussion between Tim Stitt and petsc
>>>>>>>>> matrix inversion, and it was really helpful. That was in
>>>>>>>>> last Nov.
>>>>>>>>> You can check the emails archive
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> */Yujie <recrusader at gmail.com>/* wrote:
>>>>>>>>>
>>>>>>>>>  what is the difference between sequantial and parallel AIJ
>>>>>>>>> matrix?
>>>>>>>>>  Assuming there is a matrix A, if
>>>>>>>>>  I partitaion this matrix into A1, A2, Ai... An.
>>>>>>>>>  A is a parallel AIJ matrix at the whole view, Ai
>>>>>>>>>  is a sequential AIJ matrix? I want to operate Ai at each
>>>>>>>>> node.
>>>>>>>>>  In addition, whether is it possible to get general inverse
>>>>>>>>> using
>>>>>>>>>  MatMatSolve() if the matrix is not square? Thanks a lot.
>>>>>>>>>
>>>>>>>>>  Regards,
>>>>>>>>>  Yujie
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>  On 2/4/08, *Barry Smith* <bsmith at mcs.anl.gov
>>>>>>>>>  <mailto:bsmith at mcs.anl.gov>> wrote:
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>          For sequential AIJ matrices you can fill the B matrix
>>>>>>>>> with the
>>>>>>>>>      identity and then use
>>>>>>>>>      MatMatSolve().
>>>>>>>>>
>>>>>>>>>          Note since the inverse of a sparse matrix is dense
>>>>>>>>> the B
>>>>>>>>>      matrix is
>>>>>>>>>      a SeqDense matrix.
>>>>>>>>>
>>>>>>>>>          Barry
>>>>>>>>>
>>>>>>>>>      On Feb 4, 2008, at 12:37 AM, Yujie wrote:
>>>>>>>>>
>>>>>>>>>      > Hi,
>>>>>>>>>      > Now, I want to inverse a sparse matrix. I have
>>>>>>>>> browsed the
>>>>>>>>>      manual,
>>>>>>>>>      > however, I can't find some information. could you
>>>>>>>>> give me
>>>>>>>>>      >
>>>>>>>>>      > thanks a lot.
>>>>>>>>>      >
>>>>>>>>>      > Regards,
>>>>>>>>>      > Yujie
>>>>>>>>>      >
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> ------------------------------------------------------------------------
>>>>>>>>> Looking for last minute shopping deals? Find them fast with
>>>>>>>>> Yahoo!
>>>>>>>>> Search.
>>>>>>>>> <http://us.rd.yahoo.com/evt=51734/*http://tools.search.yahoo.com/newsearch/category.php?category=shopping
>>>>>>>>> >
>>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>>
>>>>
>>>
>>
>>
>

```