how to inverse a sparse matrix in Petsc?

Lisandro Dalcin dalcinl at gmail.com
Wed Feb 6 09:53:48 CST 2008 

Well, after taking into accout Barry's comments, you have have the
following choices.

* You can use a direct method based on LU factorization using
'-ksp_type preonly -pc_type lu' . This way, PETSc will compute the LU
factors the fist time they are needed; after that, every call to
KSPSolve will reuse those factors. This will work only in sequential
with a default PETSc build, but you could also build PETSc with MUMPS,
and it will let you do the parallel factorization. For MUMPS  to
actually work in your matrix, I believe you have to add the following
line:

MatConvert(A, MATAIJMUMPS, MAT_REUSE_MATRIX, &A);

after assembling (ie. MatAssembleBegin/End calls) your Poisson matrix.


* You can use CG with '-ksp_type cg' (I assume your matrix is SPD, as
it is in a standard fractional step method), and a preconditioner. And
then, I believe the best choice for your application will bee
BoomerAMG. It has a rather high setup cost, but solves are fast. Or
your could use ML, it has less setup costs, but the solvers are a bit
slower. So if you make many timesteps, I would say that BoomerAMG will
pay.

Finally, if you use the last option, perhaps you can try Paul Fischer
tricks. I tried to add this to KSP's some time ago, but I stoped for
many reasons (the main one, lack of time). You can take a look at
this:

http://citeseer.ist.psu.edu/492082.html

A similar (equivalent?) approach is this other one (perhaps a bit
easier to implement, depending on your taste)
doi.wiley.com/10.1002/cnm.743


On 2/5/08, Ben Tay <zonexo at gmail.com> 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!
>
> 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
> > to know about your experience.
> >
> > 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.
> >>>
> >>> Waad Subber wrote:
> >>>
> >>>> Hi
> >>>> There was a discussion between Tim Stitt and petsc developers about
> >>>> matrix inversion, and it was really helpful. That was in last Nov.
> >>>> You can check the emails archive
> >>>>
> >>>> http://www-unix.mcs.anl.gov/web-mail-archive/lists/petsc-users/2007/11/threads.html
> >>>>
> >>>>
> >>>> Waad
> >>>>
> >>>> */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
> >>>>         some advice?
> >>>>         >
> >>>>         > 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>
> >>>>
> >>>
> >>>
> >>
> >
> >
> >
>
>


-- 
Lisandro Dalcín
---------------
Centro Internacional de Métodos Computacionales en Ingeniería (CIMEC)
Instituto de Desarrollo Tecnológico para la Industria Química (INTEC)
Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET)
PTLC - Güemes 3450, (3000) Santa Fe, Argentina
Tel/Fax: +54-(0)342-451.1594

More information about the petsc-users mailing list