Hello, While trying out Stefano’s PCApplyMat_MG() code (*), we stumbled upon weird numerical errors when reusing a Mat for both MatProduct_AB and MatProduct_AtB. This reminded me that there has been a long-standing issue with MatTransposeMatMult(), see https://www.mcs.anl.gov/petsc/petsc-current/src/ksp/pc/impls/hpddm/hpddm.cxx.html#line608 <https://www.mcs.anl.gov/petsc/petsc-current/src/ksp/pc/impls/hpddm/hpddm.cxx.html#line608>, that I never looked into. I’ve now been trying to figure this out, because this has side effects in multiple places (PCMG and PCHPDDM at least), and thus could impact user-code as well? With this commit: https://gitlab.com/petsc/petsc/-/commit/03d8bd538039defc2fcc3e37d523735c4aaceba0 <https://gitlab.com/petsc/petsc/-/commit/03d8bd538039defc2fcc3e37d523735c4aaceba0> + $ mpirun -n 4 src/ksp/ksp/tutorials/ex76 -ksp_converged_reason -pc_type hpddm -pc_hpddm_levels_1_eps_nev 20 -ksp_type preonly -mat_type aij -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO -rhs 2 -pc_hpddm_coarse_correction balanced -C_input_mattransposematmult -D_output_mattransposematmult I’m seeing that C is nonzero, but D is full of zeros. Mat Object: 4 MPI processes type: mpidense 5.7098316584361917e-08 1.0159399260517841e-07 1.5812349976211856e-07 2.0688121715350138e-07 2.4887556933361981e-08 4.8111092300772958e-08 1.4606298643602107e-07 1.7213611729839211e-07 […] Mat Object: 4 MPI processes type: mpidense 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 […]
If one switches to a MatType which has no MatProduct_AtB implementation with B of type MPIDense (reminder: in that case, the product is computed column-by-column), e.g., -mat_type sbaij, one gets the expected result. Mat Object: 4 MPI processes type: mpidense 7.2003197398135299e-01 9.5191869895699011e-01 6.1793966541680234e-02 9.3884397585488877e-01 1.0022337823233585e-02 2.4653068080134588e-01 1.4463931936094099e-01 8.6111517670701687e-01 Is there a bug somewhere with the MatAIJ implementation, or am I doing something which is not allowed by the MatProduct() machinery? Thanks, Pierre (*) https://gitlab.com/petsc/petsc/-/merge_requests/3717 <https://gitlab.com/petsc/petsc/-/merge_requests/3717>
