[petsc-users] MatCreateSBAIJ
Sam Guo
sam.guo at cd-adapco.com
Mon Mar 21 15:36:40 CDT 2022
Barry,
Thanks. Could you elaborate? I try to implement the matrix-vector
multiplication for a symmetric matrix using shell matrix.
Thanks,
Sam
On Mon, Mar 21, 2022 at 12:56 PM Barry Smith <bsmith at petsc.dev> wrote:
>
> The "trick" is that though "more" communication is needed to complete
> the product the communication can still be done in a single VecScatter
> instead of two separate calls to VecScatter. We simply pack both pieces of
> information that needs to be sent into a single vector.
>
> /* 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>);
>
> If you create two symmetric matrices, one with SBAIJ and one with BAIJ and
> compare the time to do the product you will find that the SBAIJ is not
> significantly slower but does save memory.
>
>
> On Mar 21, 2022, at 3:26 PM, Sam Guo <sam.guo at cd-adapco.com> wrote:
>
> Using following example from the MatCreateSBAIJ documentation
>
> 0 1 2 3 4 5 6 7 8 9 10 11
> --------------------------
> row 3 |. . . d d d o o o o o o
> row 4 |. . . d d d o o o o o o
> row 5 |. . . d d d o o o o o o
> --------------------------
>
>
> On a processor that owns rows 3, 4 and 5, rows 0-2 info are still needed. Is is that the processor that owns rows 0-2 will apply B symmetrical and then send the result
>
> to the processor that owns 3-5?
>
>
> On Mon, Mar 21, 2022 at 12:14 PM Mark Adams <mfadams at lbl.gov> wrote:
>
>> 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 <sam.guo at cd-adapco.com> 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 <sam.guo at cd-adapco.com> 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 <mfadams at lbl.gov> 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 <sam.guo at cd-adapco.com> 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
>>>>>>
>>>>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20220321/003556f2/attachment-0001.html>
More information about the petsc-users
mailing list