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> >>>>>>>> >>>>>>>> >>>>>>>> >>>>>>> >>>>>>> >>>>>> >>>>> >>>>> >>>>> >>>> >>> >>> >> > >