Hello Yann, MatGetColumnVector() is indeed kind of slow and there is a hard copy under the hood. Better to use MatDenseGetColumnVec()/MatDenseRestoreColumn() which is more efficient (https://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/Mat/MatDenseGetColumn.html <https://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/Mat/MatDenseGetColumn.html>), or do what Barry suggested.
KSPSetMatSolveBlockSize() is a rather advanced option, it’s only used by KSPHPDDM, not KSPPREONLY (I think?), which I’m guessing is what you are using with a full LU. The idea behind this option, borrowed from MUMPS option ICNTL(27), is to solve multiple systems with multiple RHS. So you split AX = B into A[X_0, X_1, …, X_p] = [B_0, B_1, …, B_p] where all blocks are of width at most KSPMatSolveBlockSize. It’s useful when you have a large number of right-hand sides and it’s becoming too memory-demanding to fit everything in memory, e.g., when doing block Krylov methods or out-of-core LU factorization. However, if everything fits, it’s best to just do KSPSetMatSolveBlockSize(ksp, PETSC_DECIDE) or don’t call that routine at all, which will default to a single solve with a single thick block. Direct solvers scale really nicely with the number of right-hand sides, thanks to BLAS3, so the higher the number of RHS, the better the efficiency. Thanks, Pierre > On 17 Sep 2020, at 11:36 PM, Yann Jobic <[email protected]> wrote: > > Hi Stefano, > You're right, i just switched the rows and colomns, and it's working fine. > However, i've done this in order to have better performance when i access the > rows, with MatGetRow. With this order, i need MatGetColumnVector, which is > believe is kind of slow, as stated in the comments of the implementation > (but i couldn't find it again, maybe i was looking at an old implementation?). > I still don't understand why i can not use the transpose of this matrix, as i > give the right parallel decomposition. It should be a bijective operation no ? > I think i'll be using the KSPMatSolve of Pierre, so i don't have to redevelop > this part. > Thanks a lot, > Yann > > > Le 9/17/2020 à 6:29 PM, Stefano Zampini a écrit : >> Yann >> you want to have the number of columns equal to the number of rhs. Not the >> number of rows. >> This is also consistent in terms of logical layouts. Simply create and dump >> the matrix this way and you can read in parallel splitting yhe rows ( ie the >> logical size per process of the vectors) >> Il Gio 17 Set 2020, 18:09 Pierre Jolivet <[email protected] >> <mailto:[email protected]>> ha scritto: >> Hello Yann, >> This is probably not fully answering your question, but the proper >> way to solve a system with N RHS is _not_ to use KSPSolve(), but >> instead KSPMatSolve(), cf. >> >> https://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/KSP/KSPMatSolve.html. >> If you are tracking master (from the GitLab repository), it’s >> available out of the box. If you are using the release tarballs, it >> will be available in 3.14.0 scheduled to be released in a couple of >> days. >> If you want to know more about the current status of block solvers >> in PETSc, please feel free to have a look at this preprint: >> http://jolivet.perso.enseeiht.fr/article.pdf >> If you are using a specific PC which is not “multi-RHS ready”, see >> the list at the top of page 5, please let me know and I’ll tell you >> how easy to support it. >> Thanks, >> Pierre
