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 at gmail.com> 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 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> > >>>> > >>> > >>> > >> > > > > > > > > -- Lisandro Dalc?n --------------- Centro Internacional de M?todos Computacionales en Ingenier?a (CIMEC) Instituto de Desarrollo Tecnol?gico para la Industria Qu?mica (INTEC) Consejo Nacional de Investigaciones Cient?ficas y T?cnicas (CONICET) PTLC - G?emes 3450, (3000) Santa Fe, Argentina Tel/Fax: +54-(0)342-451.1594
