Actual source code: mlocalref.c
petsc-3.13.6 2020-09-29
2: #include <petsc/private/matimpl.h>
4: typedef struct {
5: Mat Top;
6: PetscBool rowisblock;
7: PetscBool colisblock;
8: PetscErrorCode (*SetValues)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode);
9: PetscErrorCode (*SetValuesBlocked)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode);
10: } Mat_LocalRef;
12: /* These need to be macros because they use sizeof */
13: #define IndexSpaceGet(buf,nrow,ncol,irowm,icolm) do { \
14: if (nrow + ncol > (PetscInt)(sizeof(buf)/sizeof(buf[0]))) { \
15: PetscMalloc2(nrow,&irowm,ncol,&icolm); \
16: } else { \
17: irowm = &buf[0]; \
18: icolm = &buf[nrow]; \
19: } \
20: } while (0)
22: #define IndexSpaceRestore(buf,nrow,ncol,irowm,icolm) do { \
23: if (nrow + ncol > (PetscInt)(sizeof(buf)/sizeof(buf[0]))) { \
24: PetscFree2(irowm,icolm); \
25: } \
26: } while (0)
28: static void BlockIndicesExpand(PetscInt n,const PetscInt idx[],PetscInt bs,PetscInt idxm[])
29: {
30: PetscInt i,j;
31: for (i=0; i<n; i++) {
32: for (j=0; j<bs; j++) {
33: idxm[i*bs+j] = idx[i]*bs + j;
34: }
35: }
36: }
38: static PetscErrorCode MatSetValuesBlockedLocal_LocalRef_Block(Mat A,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
39: {
40: Mat_LocalRef *lr = (Mat_LocalRef*)A->data;
42: PetscInt buf[4096],*irowm=NULL,*icolm; /* suppress maybe-uninitialized warning */
45: if (!nrow || !ncol) return(0);
46: IndexSpaceGet(buf,nrow,ncol,irowm,icolm);
47: ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,nrow,irow,irowm);
48: ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,ncol,icol,icolm);
49: (*lr->SetValuesBlocked)(lr->Top,nrow,irowm,ncol,icolm,y,addv);
50: IndexSpaceRestore(buf,nrow,ncol,irowm,icolm);
51: return(0);
52: }
54: static PetscErrorCode MatSetValuesBlockedLocal_LocalRef_Scalar(Mat A,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
55: {
56: Mat_LocalRef *lr = (Mat_LocalRef*)A->data;
58: PetscInt rbs,cbs,buf[4096],*irowm,*icolm;
61: MatGetBlockSizes(A,&rbs,&cbs);
62: IndexSpaceGet(buf,nrow*rbs,ncol*cbs,irowm,icolm);
63: BlockIndicesExpand(nrow,irow,rbs,irowm);
64: BlockIndicesExpand(ncol,icol,cbs,icolm);
65: ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,nrow*rbs,irowm,irowm);
66: ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,ncol*cbs,icolm,icolm);
67: (*lr->SetValues)(lr->Top,nrow*rbs,irowm,ncol*cbs,icolm,y,addv);
68: IndexSpaceRestore(buf,nrow*rbs,ncol*cbs,irowm,icolm);
69: return(0);
70: }
72: static PetscErrorCode MatSetValuesLocal_LocalRef_Scalar(Mat A,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
73: {
74: Mat_LocalRef *lr = (Mat_LocalRef*)A->data;
76: PetscInt buf[4096],*irowm,*icolm;
79: IndexSpaceGet(buf,nrow,ncol,irowm,icolm);
80: /* If the row IS defining this submatrix was an ISBLOCK, then the unblocked LGMapApply is the right one to use. If
81: * instead it was (say) an ISSTRIDE with a block size > 1, then we need to use LGMapApplyBlock */
82: if (lr->rowisblock) {
83: ISLocalToGlobalMappingApply(A->rmap->mapping,nrow,irow,irowm);
84: } else {
85: ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,nrow,irow,irowm);
86: }
87: /* As above, but for the column IS. */
88: if (lr->colisblock) {
89: ISLocalToGlobalMappingApply(A->cmap->mapping,ncol,icol,icolm);
90: } else {
91: ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,ncol,icol,icolm);
92: }
93: (*lr->SetValues)(lr->Top,nrow,irowm,ncol,icolm,y,addv);
94: IndexSpaceRestore(buf,nrow,ncol,irowm,icolm);
95: return(0);
96: }
98: /* Compose an IS with an ISLocalToGlobalMapping to map from IS source indices to global indices */
99: static PetscErrorCode ISL2GCompose(IS is,ISLocalToGlobalMapping ltog,ISLocalToGlobalMapping *cltog)
100: {
102: const PetscInt *idx;
103: PetscInt m,*idxm;
104: PetscInt bs;
105: PetscBool isblock;
111: PetscObjectTypeCompare((PetscObject)is,ISBLOCK,&isblock);
112: if (isblock) {
113: PetscInt lbs;
115: ISGetBlockSize(is,&bs);
116: ISLocalToGlobalMappingGetBlockSize(ltog,&lbs);
117: if (bs == lbs) {
118: ISGetLocalSize(is,&m);
119: m = m/bs;
120: ISBlockGetIndices(is,&idx);
121: PetscMalloc1(m,&idxm);
122: ISLocalToGlobalMappingApplyBlock(ltog,m,idx,idxm);
123: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)is),bs,m,idxm,PETSC_OWN_POINTER,cltog);
124: ISBlockRestoreIndices(is,&idx);
125: return(0);
126: }
127: }
128: ISGetLocalSize(is,&m);
129: ISGetIndices(is,&idx);
130: ISGetBlockSize(is,&bs);
131: PetscMalloc1(m,&idxm);
132: if (ltog) {
133: ISLocalToGlobalMappingApply(ltog,m,idx,idxm);
134: } else {
135: PetscArraycpy(idxm,idx,m);
136: }
137: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)is),bs,m,idxm,PETSC_OWN_POINTER,cltog);
138: ISRestoreIndices(is,&idx);
139: return(0);
140: }
142: static PetscErrorCode ISL2GComposeBlock(IS is,ISLocalToGlobalMapping ltog,ISLocalToGlobalMapping *cltog)
143: {
145: const PetscInt *idx;
146: PetscInt m,*idxm,bs;
152: ISBlockGetLocalSize(is,&m);
153: ISBlockGetIndices(is,&idx);
154: ISLocalToGlobalMappingGetBlockSize(ltog,&bs);
155: PetscMalloc1(m,&idxm);
156: if (ltog) {
157: ISLocalToGlobalMappingApplyBlock(ltog,m,idx,idxm);
158: } else {
159: PetscArraycpy(idxm,idx,m);
160: }
161: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)is),bs,m,idxm,PETSC_OWN_POINTER,cltog);
162: ISBlockRestoreIndices(is,&idx);
163: return(0);
164: }
166: static PetscErrorCode MatDestroy_LocalRef(Mat B)
167: {
171: PetscFree(B->data);
172: return(0);
173: }
176: /*@
177: MatCreateLocalRef - Gets a logical reference to a local submatrix, for use in assembly
179: Not Collective
181: Input Arguments:
182: + A - Full matrix, generally parallel
183: . isrow - Local index set for the rows
184: - iscol - Local index set for the columns
186: Output Arguments:
187: . newmat - New serial Mat
189: Level: developer
191: Notes:
192: Most will use MatGetLocalSubMatrix() which returns a real matrix corresponding to the local
193: block if it available, such as with matrix formats that store these blocks separately.
195: The new matrix forwards MatSetValuesLocal() and MatSetValuesBlockedLocal() to the global system.
196: In general, it does not define MatMult() or any other functions. Local submatrices can be nested.
198: .seealso: MatSetValuesLocal(), MatSetValuesBlockedLocal(), MatGetLocalSubMatrix(), MatCreateSubMatrix()
199: @*/
200: PetscErrorCode MatCreateLocalRef(Mat A,IS isrow,IS iscol,Mat *newmat)
201: {
203: Mat_LocalRef *lr;
204: Mat B;
205: PetscInt m,n;
206: PetscBool islr;
213: if (!A->rmap->mapping) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Matrix must have local to global mapping provided before this call");
214: *newmat = 0;
216: MatCreate(PETSC_COMM_SELF,&B);
217: ISGetLocalSize(isrow,&m);
218: ISGetLocalSize(iscol,&n);
219: MatSetSizes(B,m,n,m,n);
220: PetscObjectChangeTypeName((PetscObject)B,MATLOCALREF);
221: MatSetUp(B);
223: B->ops->destroy = MatDestroy_LocalRef;
225: PetscNewLog(B,&lr);
226: B->data = (void*)lr;
228: PetscObjectTypeCompare((PetscObject)A,MATLOCALREF,&islr);
229: if (islr) {
230: Mat_LocalRef *alr = (Mat_LocalRef*)A->data;
231: lr->Top = alr->Top;
232: } else {
233: /* This does not increase the reference count because MatLocalRef is not allowed to live longer than its parent */
234: lr->Top = A;
235: }
236: {
237: ISLocalToGlobalMapping rltog,cltog;
238: PetscInt arbs,acbs,rbs,cbs;
240: /* We will translate directly to global indices for the top level */
241: lr->SetValues = MatSetValues;
242: lr->SetValuesBlocked = MatSetValuesBlocked;
244: B->ops->setvalueslocal = MatSetValuesLocal_LocalRef_Scalar;
246: ISL2GCompose(isrow,A->rmap->mapping,&rltog);
247: if (isrow == iscol && A->rmap->mapping == A->cmap->mapping) {
248: PetscObjectReference((PetscObject)rltog);
249: cltog = rltog;
250: } else {
251: ISL2GCompose(iscol,A->cmap->mapping,&cltog);
252: }
253: /* Remember if the ISes we used to pull out the submatrix are of type ISBLOCK. Will be used later in
254: * MatSetValuesLocal_LocalRef_Scalar. */
255: PetscObjectTypeCompare((PetscObject)isrow,ISBLOCK,&lr->rowisblock);
256: PetscObjectTypeCompare((PetscObject)iscol,ISBLOCK,&lr->colisblock);
257: MatSetLocalToGlobalMapping(B,rltog,cltog);
258: ISLocalToGlobalMappingDestroy(&rltog);
259: ISLocalToGlobalMappingDestroy(&cltog);
261: MatGetBlockSizes(A,&arbs,&acbs);
262: ISGetBlockSize(isrow,&rbs);
263: ISGetBlockSize(iscol,&cbs);
264: /* Always support block interface insertion on submatrix */
265: PetscLayoutSetBlockSize(B->rmap,rbs);
266: PetscLayoutSetBlockSize(B->cmap,cbs);
267: if (arbs != rbs || acbs != cbs || (arbs == 1 && acbs == 1)) {
268: /* Top-level matrix has different block size, so we have to call its scalar insertion interface */
269: B->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_LocalRef_Scalar;
270: } else {
271: /* Block sizes match so we can forward values to the top level using the block interface */
272: B->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_LocalRef_Block;
274: ISL2GComposeBlock(isrow,A->rmap->mapping,&rltog);
275: if (isrow == iscol && A->rmap->mapping == A->cmap->mapping) {
276: PetscObjectReference((PetscObject)rltog);
277: cltog = rltog;
278: } else {
279: ISL2GComposeBlock(iscol,A->cmap->mapping,&cltog);
280: }
281: MatSetLocalToGlobalMapping(B,rltog,cltog);
282: ISLocalToGlobalMappingDestroy(&rltog);
283: ISLocalToGlobalMappingDestroy(&cltog);
284: }
285: }
286: *newmat = B;
287: return(0);
288: }