[petsc-dev] PCREDUNDANT

[Pierre Jolivet]

> On 28 Jul 2019, at 3:02 PM, Mark Adams <mfadams at lbl.gov> wrote:
> 
> 
> On Sun, Jul 28, 2019 at 2:54 AM Pierre Jolivet via petsc-dev <petsc-dev at mcs.anl.gov <mailto:petsc-dev at mcs.anl.gov>> wrote:
> Hello,
> I’m facing multiple issues with PCREDUNDANT and MATMPISBAIJ:
> 1) https://www.mcs.anl.gov/petsc/petsc-current/src/mat/impls/sbaij/mpi/mpisbaij.c.html#line3354 <https://www.mcs.anl.gov/petsc/petsc-current/src/mat/impls/sbaij/mpi/mpisbaij.c.html#line3354> shouldn’t this be sum != N? I’m running an example where it says that sum (4) != Nbs (60), with a bs=15.
> 
> Clearly a bug.

OK, fixed here: https://bitbucket.org/petsc/petsc/pull-requests/1915/fixes/diff#Lsrc/mat/impls/sbaij/mpi/mpisbaij.cT3324 <https://bitbucket.org/petsc/petsc/pull-requests/1915/fixes/diff#Lsrc/mat/impls/sbaij/mpi/mpisbaij.cT3324> (once again, this isn’t tested in PETSc for a bs>1, I think, but now my preconditioner is working OK, for a bs>1).

>  
> 2) when I’m using MATMPIBAIJ, I can do stuff like: -prefix_mat_type baij -prefix_pc_type redundant -prefix_redundant_pc_type ilu, and in the KSPView, I have "package used to perform factorization: petsc”, so the underlying MatType is indeed MATSEQBAIJ. 
> However, with MATMPISBAIJ, if I do: -prefix_mat_type sbaij -prefix_pc_type redundant, first, it looks like you are hardwiring a PCLU (MatGetFactor() line 4440 in src/mat/interface/matrix.c
> 
> Using LU as a default for symmetric matrices does seem wrong.

FWIW, I made a similar remark to Jose because the same matrices are being used by SLEPc for a GEVP with a SINVERT ST in a Krylov—Schur EPS (tentative fix here: https://bitbucket.org/slepc/slepc/branch/jose/sinvert-cholesky#diff <https://bitbucket.org/slepc/slepc/branch/jose/sinvert-cholesky#diff>).

> Could not locate a solver package.), then, if I append -prefix_redundant_pc_type cholesky, I end up with an error related to MUMPS: MatGetFactor_sbaij_mumps() line 2625 in src/mat/impls/aij/mpi/mumps/mumps.c Cannot use PETSc SBAIJ matrices with block size > 1 with MUMPS Cholesky, use AIJ matrix instead. Why isn’t this call dispatched to PETSc Cholesky for SeqSBAIJ matrices?
> 
> 
> Generally, we don't like to switch parameters under the covers like this. We would rather you get your inputs right so you know what is going on.

I’m not sure how to interpret your comment, but anyway, I found out that for SBAIJ matrices, the type of the output matrix in MatCreateMPIMatConcatenateSeqMat (used by PCREDUNDANT) was set to MPISBAIJ, whereas for BAIJ, it was BAIJ (and thus SeqBAIJ when used within PCREDUNDANT with “full” redundancy over PETSC_COMM_WORLD).
Here is the fix: https://bitbucket.org/petsc/petsc/pull-requests/1915/fixes/diff#Lsrc/mat/impls/sbaij/mpi/mpisbaij.cT3345 <https://bitbucket.org/petsc/petsc/pull-requests/1915/fixes/diff#Lsrc/mat/impls/sbaij/mpi/mpisbaij.cT3345> (inspired by https://www.mcs.anl.gov/petsc/petsc-current/src/mat/impls/baij/mpi/mpibaij.c.html#line3957 <https://www.mcs.anl.gov/petsc/petsc-current/src/mat/impls/baij/mpi/mpibaij.c.html#line3957>)

Now everything seems OK!
Thanks,
Pierre

>  
> Thanks,
> Pierre
> 
> 1) I don’t think this is tested right now, at least not in src/ksp/ksp/examples/tutorials
> 2) reproducer: src/ksp/ksp/examples/tutorials/ex2.c
> $ mpirun -np 2 ./ex2 -pc_type redundant -mat_type sbaij
> // error because trying to do LU with a symmetric matrix
> $ mpirun -np 2 ./ex2 -pc_type redundant -mat_type sbaij -redundant_pc_type cholesky -ksp_view
> // you’ll see: that MUMPS is being used, but since bs=1, it’s working, but it won’t for the general case
> //                  the MatType is mpisbaij with "1 MPI processes" whereas with baij, it’s seqbaij
> 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20190728/2a44878d/attachment.html>