[petsc-users] Cholesky: matrix ordering with aij? Bug with debugging=1 for ordering on sbaij?

Thu Nov 27 13:22:10 CST 2014

> On Nov 27, 2014, at 12:52 PM, Eric Chamberland <Eric.Chamberland at giref.ulaval.ca> wrote:
> 
> Hi,
> 
> Using petsc 3.5.2.
> 
> ------------------------
> Question #1:
> ------------------------
> 
> I configured a cholesky like this:
> 
> KSP Object:(o_slin) 1 MPI processes
>  type: preonly
>  maximum iterations=4000, initial guess is zero
>  tolerances:  relative=1e-12, absolute=1e-20, divergence=1e+30
>  left preconditioning
>  using NONE norm type for convergence test
> PC Object:(o_slin) 1 MPI processes
>  type: cholesky
>    Cholesky: out-of-place factorization
>    tolerance for zero pivot 2.22045e-14
>    matrix ordering: nd
>    factor fill ratio given 5, needed 2.59172
>      Factored matrix follows:
>        Mat Object:         1 MPI processes
>          type: seqsbaij
>          rows=693, cols=693
>          package used to perform factorization: petsc
>          total: nonzeros=10328, allocated nonzeros=10328
>          total number of mallocs used during MatSetValues calls =0
>              block size is 1
>  linear system matrix = precond matrix:
>  Mat Object:  (o_slin)   1 MPI processes
>    type: seqaij
>    rows=693, cols=693
>    total: nonzeros=7277, allocated nonzeros=7345
>    total number of mallocs used during MatSetValues calls =0
>      not using I-node routines
> 
> with a "mat_type aij" for the assembled matrix and "pc_factor_mat_ordering_type nd".  The PC creates an "sbaij".
> 
> When/where is the "matrix ordering: nd" used since it is documented to not works on sbaij matrix?

   The factorization is done on the aij matrix and the result put into a sbaij matrix. Since with aij it is easy to reorder rows/columns it can do the factorization in any order it wants. sbaij only stores the upper half of the matrix so factorizing in any ordering besides "natural" is nasty and we don't do it.

> 
> Compare the number of non-zeros for the factored matrix obtained with a "natural ordering":
> 
> KSP Object:(o_slin) 1 MPI processes
>  type: preonly
>  maximum iterations=4000, initial guess is zero
>  tolerances:  relative=1e-12, absolute=1e-20, divergence=1e+30
>  left preconditioning
>  using NONE norm type for convergence test
> PC Object:(o_slin) 1 MPI processes
>  type: cholesky
>    Cholesky: out-of-place factorization
>    tolerance for zero pivot 2.22045e-14
>    matrix ordering: natural
>    factor fill ratio given 5, needed 44.1099
>      Factored matrix follows:
>        Mat Object:         1 MPI processes
>          type: seqsbaij
>          rows=693, cols=693
>          package used to perform factorization: petsc
>          total: nonzeros=175778, allocated nonzeros=175778
>          total number of mallocs used during MatSetValues calls =0
>              block size is 1
>  linear system matrix = precond matrix:
>  Mat Object:  (o_slin)   1 MPI processes
>    type: seqaij
>    rows=693, cols=693
>    total: nonzeros=7277, allocated nonzeros=7345
>    total number of mallocs used during MatSetValues calls =0
>      not using I-node routines
> 
> there is really something happening with "nd"...
> 
> ------------------------
> Question #2:
> ------------------------
> Using petsc 3.5.2 compiled with "--with-debugging=no", if a try to use "mat_type sbaij", I get this friendly error:
> 
> [0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
> [0]PETSC ERROR: No support for this operation for this object type
> [0]PETSC ERROR: Matrix reordering is not supported for sbaij matrix. Use aij format
> [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
> [0]PETSC ERROR: Petsc Release Version 3.5.2, Sep, 08, 2014
> [0]PETSC ERROR: /home/mefpp_ericc/depots_prepush/GIREF/bin/probDiffusion.dev on a arch-linux2-c-opt named melkor by eric Thu Nov 27 13:35:01 2014
> [0]PETSC ERROR: Configure options --with-mpi-compilers=1 --with-make-np=12 --prefix=/opt/petsc-3.5.2_opt_openmpi-1.6.5_mkl_mt --with-mpi-dir=/opt/openmpi-1.6.5 --with-debugging=no --download-ml=yes --with-mkl_pardiso=1 --with-mkl_pardiso-dir=/opt/intel/composerxe/mkl --download-mumps=yes --download-superlu=yes --download-superlu_dist=yes --download-parmetis=yes --download-metis=yes --download-hypre=yes --with-shared-libraries=1 --with-scalapack=1 --with-scalapack-include=/opt/intel/composerxe/mkl/include --with-scalapack-lib="-L/opt/intel/composerxe/mkl/lib/intel64 -lmkl_scalapack_lp64 -lmkl_blacs_intelmpi_lp64" --with-blas-lapack-dir=/opt/intel/composerxe/mkl/lib/intel64
> [0]PETSC ERROR: #1 MatCholeskyFactorSymbolic_SeqSBAIJ() line 245 in /home/mefpp_ericc/petsc-3.5.2/src/mat/impls/sbaij/seq/sbaijfact.c
> [0]PETSC ERROR: #2 MatCholeskyFactorSymbolic() line 3012 in /home/mefpp_ericc/petsc-3.5.2/src/mat/interface/matrix.c
> [0]PETSC ERROR: #3 PCSetUp_Cholesky() line 124 in /home/mefpp_ericc/petsc-3.5.2/src/ksp/pc/impls/factor/cholesky/cholesky.c
> [0]PETSC ERROR: #4 PCSetUp() line 902 in /home/mefpp_ericc/petsc-3.5.2/src/ksp/pc/interface/precon.c
> [0]PETSC ERROR: #5 KSPSetUp() line 305 in /home/mefpp_ericc/petsc-3.5.2/src/ksp/ksp/interface/itfunc.c
> 
> But if i use petsc compiled with "--with-debugging=yes", I got this one:
> 
> [0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
> [0]PETSC ERROR: Invalid argument
> [0]PETSC ERROR: Index set is not a permutation
> [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
> [0]PETSC ERROR: Petsc Release Version 3.5.2, Sep, 08, 2014
> [0]PETSC ERROR: probDiffusion.dev on a arch-linux2-c-debug named melkor by eric Thu Nov 27 13:45:04 2014
> [0]PETSC ERROR: Configure options --with-mpi-compilers=1 --with-make-np=12 --prefix=/opt/petsc-3.5.2_debug_openmpi-1.6.5_mkl_mt --with-mpi-dir=/opt/openmpi-1.6.5 --with-debugging=yes --download-ml=yes --with-mkl_pardiso=1 --with-mkl_pardiso-dir=/opt/intel/composerxe/mkl --download-mumps=yes --download-superlu=yes --download-superlu_dist=yes --download-parmetis=yes --download-metis=yes --download-hypre=yes --with-shared-libraries=1 --with-scalapack=1 --with-scalapack-include=/opt/intel/composerxe/mkl/include --with-scalapack-lib="-L/opt/intel/composerxe/mkl/lib/intel64 -lmkl_scalapack_lp64 -lmkl_blacs_intelmpi_lp64" --with-blas-lapack-dir=/opt/intel/composerxe/mkl/lib/intel64
> [0]PETSC ERROR: #1 ISSetPermutation() line 183 in /home/mefpp_ericc/petsc-3.5.2/src/vec/is/is/interface/index.c
> [0]PETSC ERROR: #2 MatGetOrdering() line 261 in /home/mefpp_ericc/petsc-3.5.2/src/mat/order/sorder.c
> [0]PETSC ERROR: #3 PCSetUp_Cholesky() line 107 in /home/mefpp_ericc/petsc-3.5.2/src/ksp/pc/impls/factor/cholesky/cholesky.c
> [0]PETSC ERROR: #4 PCSetUp() line 902 in /home/mefpp_ericc/petsc-3.5.2/src/ksp/pc/interface/precon.c
> [0]PETSC ERROR: #5 KSPSetUp() line 305 in /home/mefpp_ericc/petsc-3.5.2/src/ksp/ksp/interface/itfunc.c
> 
> Is this the expected behavior?

   Well kind of. When debugging is turned off the the error checking that an ordering is actually a permutation is turned off hence the code gets further (until it cannot do the factorization). Hence two different error messages.

   Now I suspect that it is generating a bad ordering with the sbaij matrix so probably it would be best if we errored out in the MatGetOrdering() always. But even in the current situation nothing is really going wrong.

  Barry

> 
> Thanks,
> 
> Eric