> On Aug 24, 2016, at 8:07 PM, Harshad Ranadive <harshadranad...@gmail.com> > wrote: > > Thanks a lot Barry for your help. I have resolved the issue of slow linear > solver for excessively large RHS vectors to some degree. > > To summarize: > > 1) In my case, I wanted to solve the same linear system a large number of > times for different RHS vectors. Previously I was using the iterative solver > KSPGMRES with preconditioner PCJACOBI > > 2) I had a huge number of RHS to be solved (~1M/time step) .... this had a > computational cost 9.6 times greater than my explicit code (which did not > need matrix inversion - only reqd. RHS eval) > > 3) Based on Barry's suggestions to use a direct solver I changed the KSP type > to KSPPREONLY and used the preconditioner - PCLU. > > 4) This approach has a small drawback that the linear system needs to be > solved sequentially. Although different systems are now solved by different > processors in parallel.
Note that if you ./configure PETSc with --download-superlu_dist then you can actually solve them in parallel. > > 5) The current computational cost is only 1.16 times the explicit code. > > Thanks and Regards, > Harshad > > > On Fri, Aug 12, 2016 at 1:27 PM, Barry Smith <bsm...@mcs.anl.gov> wrote: > > > On Aug 11, 2016, at 10:14 PM, Harshad Ranadive <harshadranad...@gmail.com> > > wrote: > > > > Hi Barry, > > > > Thanks for this recommendation. > > > > As you mention, the matrix factorization should be on a single processor. > > If the factored matrix A is available on all processors can I then use > > MatMatSolve(A,B,X) in parallel? That is could the RHS block matrix 'B' and > > solution matrix 'X' be distributed in different processors as is done while > > using MatCreateDense(...) ? > > Note sure what you mean. > > You can have different processes handle different right hand sides. So give > the full linear system matrix to each process; each process factors it and > then each process solves a different set of right hand sides. Embarrassingly > parallel except for any communication you need to do to get the matrix and > right hand sides to the right processes. > > If the linear system involves say millions of unknowns this is the way to > go. If the linear system is over say 1 billion unknowns then it might be > worth each linear system in parallel. > > Barry > > > > > Thanks, > > Harshad > > > > > > > > On Fri, Aug 12, 2016 at 2:09 AM, Barry Smith <bsm...@mcs.anl.gov> wrote: > > > > If it is sequential, which it probably should be, then you can you > > MatLUFactorSymbolic(), MatLUFactorNumeric() and MatMatSolve() where you put > > a bunch of your right hand side vectors into a dense array; not all million > > of them but maybe 10 to 100 at a time. > > > > Barry > > > > > On Aug 10, 2016, at 10:18 PM, Harshad Ranadive > > > <harshadranad...@gmail.com> wrote: > > > > > > Hi Barry > > > > > > The matrix A is mostly tridiagonal > > > > > > 1 α 0 ......... 0 > > > > > > α 1 α 0 .......0 > > > > > > > > > 0 α 1 α 0 ....0 > > > > > > > > > .................... > > > 0..............α 1 > > > > > > In some cases (periodic boundaries) there would be an 'α' in > > > right-top-corner and left-bottom corner. > > > > > > I am not using multigrid approach. I just implemented an implicit > > > filtering approach (instead of an explicit existing one) which requires > > > the solution of the above system. > > > > > > Thanks > > > Harshad > > > > > > On Thu, Aug 11, 2016 at 1:07 PM, Barry Smith <bsm...@mcs.anl.gov> wrote: > > > > > > Effectively utilizing multiple right hand sides with the same system > > > can result in roughly 2 or at absolute most 3 times improvement in > > > solve time. A great improvement but when you have a million right hand > > > sides not a giant improvement. > > > > > > The first step is to get the best (most efficient) preconditioner for > > > you problem. Since you have many right hand sides it obviously pays to > > > spend more time building the preconditioner so that each solve is faster. > > > If you provide more information on your linear system we might have > > > suggestions. CFD so is your linear system a Poisson problem? Are you > > > using geometric or algebraic multigrid with PETSc? It not a Poisson > > > problem how can you describe the linear system? > > > > > > Barry > > > > > > > > > > > > > On Aug 10, 2016, at 9:54 PM, Harshad Ranadive > > > > <harshadranad...@gmail.com> wrote: > > > > > > > > Hi All, > > > > > > > > I have currently added the PETSc library with our CFD solver. > > > > > > > > In this I need to use KSPSolve(...) multiple time for the same matrix > > > > A. I have read that PETSc does not support passing multiple RHS vectors > > > > in the form of a matrix and the only solution to this is calling > > > > KSPSolve multiple times as in example given here: > > > > http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex16.c.html > > > > > > > > I have followed this technique, but I find that the performance of the > > > > code is very slow now. I basically have a mesh size of 8-10 Million and > > > > I need to solve the matrix A very large number of times. I have checked > > > > that the statement KSPSolve(..) is taking close to 90% of my > > > > computation time. > > > > > > > > I am setting up the matrix A, KSPCreate, KSPSetup etc just once at the > > > > start. Only the following statements are executed in a repeated loop > > > > > > > > Loop begin: (say million times !!) > > > > > > > > loop over vector length > > > > VecSetValues( ....) > > > > end > > > > > > > > VecAssemblyBegin( ... ) > > > > VecAssemblyEnd (...) > > > > > > > > KSPSolve (...) > > > > > > > > VecGetValues > > > > > > > > Loop end. > > > > > > > > Is there an efficient way of doing this rather than using KSPSolve > > > > multiple times? > > > > > > > > Please note my matrix A never changes during the time steps or across > > > > the mesh ... So essentially if I can get the inverse once would it be > > > > good enough? It has been recommended in the FAQ that matrix inverse > > > > should be avoided but would it be okay to use in my case? > > > > > > > > Also could someone please provide an example of how to use MatLUFactor > > > > and MatCholeskyFactor() to find the matrix inverse... the arguments > > > > below were not clear to me. > > > > IS row > > > > IS col > > > > const MatFactorInfo *info > > > > > > > > Apologies for a long email and thanks to anyone for help. > > > > > > > > Regards > > > > Harshad > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > >
