Actual source code: mlocalref.c

petsc-3.9.0 2018-04-07
Report Typos and Errors

  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:     PetscMemcpy(idxm,idx,m*sizeof(PetscInt));
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:     PetscMemcpy(idxm,idx,m*sizeof(PetscInt));
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: }