how to inverse a sparse matrix in Petsc?
Barry Smith
bsmith at mcs.anl.gov
Wed Feb 6 07:52:39 CST 2008
Whoops, actually the grids in each direction would need to be
uniform.
Barry
On Feb 5, 2008, at 9:39 PM, Ben Tay wrote:
> 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