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