how to inverse a sparse matrix in Petsc?
Ben Tay
zonexo at gmail.com
Wed Feb 6 10:24:20 CST 2008
Hi Lisandro,
Thanks for your recommendation. Btw, does the poisson eqn arising from
fractional step gives a matrix which is SPD? Because my grid's are
non-uniform in both x,y directions. Shouldn't that result in a
non-symmetric matrix? But I think it's still PD, positive definite.
Correct me if I'm wrong.
Thanks
Lisandro Dalcin wrote:
> 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>
>>>>>>
>>>>>>
>>>>>
>>>
>>>
>>
>
>
>
More information about the petsc-users
mailing list