how to inverse a sparse matrix in Petsc?
Ben Tay
zonexo at gmail.com
Tue Feb 5 20:48:38 CST 2008
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. 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