PETSc stores parallel matrices as two serial matrices. One for the diagonal (d or A) block and one for the rest (o or B). I would guess that for symmetric matrices it has a symmetric matrix for the diagonal and a full AIJ matrix for the (upper) off-diagonal. So the multtranspose is applying B symmetrically. This lower off-diagonal and the diagonal block can be done without communication. Then the off processor values are collected, and the upper off-diagonal is applied.
On Mon, Mar 21, 2022 at 2:35 PM Sam Guo <[email protected]> wrote: > I am most interested in how the lower triangular part is redistributed. It > seems that SBAJI saves memory but requires more communication than BAIJ. > > On Mon, Mar 21, 2022 at 11:27 AM Sam Guo <[email protected]> wrote: > >> Mark, thanks for the quick response. I am more interested in parallel >> implementation of MatMult for SBAIJ. I found following >> >> 1094: *PetscErrorCode >> <https://petsc.org/main/docs/manualpages/Sys/PetscErrorCode.html#PetscErrorCode> >> MatMult_MPISBAIJ(Mat >> <https://petsc.org/main/docs/manualpages/Mat/Mat.html#Mat> A,Vec >> <https://petsc.org/main/docs/manualpages/Vec/Vec.html#Vec> xx,Vec >> <https://petsc.org/main/docs/manualpages/Vec/Vec.html#Vec> yy)*1095: {1096: >> Mat_MPISBAIJ *a = (Mat_MPISBAIJ*)A->data;1097: PetscErrorCode >> <https://petsc.org/main/docs/manualpages/Sys/PetscErrorCode.html#PetscErrorCode> >> ierr;1098: PetscInt >> <https://petsc.org/main/docs/manualpages/Sys/PetscInt.html#PetscInt> >> mbs=a->mbs,bs=A->rmap->bs;1099: PetscScalar >> <https://petsc.org/main/docs/manualpages/Sys/PetscScalar.html#PetscScalar> >> *from;1100: const PetscScalar >> <https://petsc.org/main/docs/manualpages/Sys/PetscScalar.html#PetscScalar> >> *x; >> 1103: /* diagonal part */1104: >> (*a->A->ops->mult)(a->A,xx,a->slvec1a);1105: VecSet >> <https://petsc.org/main/docs/manualpages/Vec/VecSet.html#VecSet>(a->slvec1b,0.0); >> 1107: /* subdiagonal part */1108: >> (*a->B->ops->multtranspose)(a->B,xx,a->slvec0b); >> 1110: /* copy x into the vec slvec0 */1111: VecGetArray >> <https://petsc.org/main/docs/manualpages/Vec/VecGetArray.html#VecGetArray>(a->slvec0,&from);1112: >> VecGetArrayRead >> <https://petsc.org/main/docs/manualpages/Vec/VecGetArrayRead.html#VecGetArrayRead>(xx,&x); >> 1114: PetscArraycpy >> <https://petsc.org/main/docs/manualpages/Sys/PetscArraycpy.html#PetscArraycpy>(from,x,bs*mbs);1115: >> VecRestoreArray >> <https://petsc.org/main/docs/manualpages/Vec/VecRestoreArray.html#VecRestoreArray>(a->slvec0,&from);1116: >> VecRestoreArrayRead >> <https://petsc.org/main/docs/manualpages/Vec/VecRestoreArrayRead.html#VecRestoreArrayRead>(xx,&x); >> 1118: VecScatterBegin >> <https://petsc.org/main/docs/manualpages/PetscSF/VecScatterBegin.html#VecScatterBegin>(a->sMvctx,a->slvec0,a->slvec1,ADD_VALUES >> >> <https://petsc.org/main/docs/manualpages/Sys/ADD_VALUES.html#ADD_VALUES>,SCATTER_FORWARD >> >> <https://petsc.org/main/docs/manualpages/Vec/SCATTER_FORWARD.html#SCATTER_FORWARD>);1119: >> VecScatterEnd >> <https://petsc.org/main/docs/manualpages/PetscSF/VecScatterEnd.html#VecScatterEnd>(a->sMvctx,a->slvec0,a->slvec1,ADD_VALUES >> >> <https://petsc.org/main/docs/manualpages/Sys/ADD_VALUES.html#ADD_VALUES>,SCATTER_FORWARD >> >> <https://petsc.org/main/docs/manualpages/Vec/SCATTER_FORWARD.html#SCATTER_FORWARD>);1120: >> /* supperdiagonal part */1121: >> (*a->B->ops->multadd)(a->B,a->slvec1b,a->slvec1a,yy);1122: return(0);1123: >> } >> >> I try to understand the algorithm. >> >> >> Thanks, >> >> Sam >> >> >> On Mon, Mar 21, 2022 at 11:14 AM Mark Adams <[email protected]> wrote: >> >>> This code looks fine to me and the code is >>> in src/mat/impls/sbaij/seq/sbaij2.c >>> >>> On Mon, Mar 21, 2022 at 2:02 PM Sam Guo <[email protected]> wrote: >>> >>>> Dear PETSc dev team, >>>> The documentation about MatCreateSBAIJ has following >>>> "It is recommended that one use the MatCreate >>>> <https://petsc.org/main/docs/manualpages/Mat/MatCreate.html#MatCreate> >>>> (), MatSetType >>>> <https://petsc.org/main/docs/manualpages/Mat/MatSetType.html#MatSetType>() >>>> and/or MatSetFromOptions >>>> <https://petsc.org/main/docs/manualpages/Mat/MatSetFromOptions.html#MatSetFromOptions>(), >>>> MatXXXXSetPreallocation() paradigm instead of this routine directly. >>>> [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation >>>> <https://petsc.org/main/docs/manualpages/Mat/MatSeqAIJSetPreallocation.html#MatSeqAIJSetPreallocation> >>>> ]" >>>> I currently call MatCreateSBAIJ directly as follows: >>>> MatCreateSBAIJ (with d_nnz and o_nnz) >>>> MatSetValues (to add row by row) >>>> MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); >>>> MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); >>>> MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE); >>>> >>>> Two questions: >>>> (1) I am wondering whether what I am doing is the most efficient. >>>> >>>> (2) I try to find out how the matrix vector multiplication is >>>> implemented in PETSc for SBAIJ storage. >>>> >>>> Thanks, >>>> Sam >>>> >>>
