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@xxxxxxxxx> 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@xxxxxxxxx> 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@xxxxxxxxx>/* 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@xxxxxxxxxxx
<mailto:bsmith@xxxxxxxxxxx>> 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>