how to inverse a sparse matrix in Petsc?

Barry Smith bsmith at mcs.anl.gov
Tue Feb 5 21:21:46 CST 2008 

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
>>>> 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