Actual source code: mlocalref.c

petsc-3.5.2 2014-09-08
Report Typos and Errors
  2: #include <petsc-private/matimpl.h>          /*I "petscmat.h" I*/

  4: typedef struct {
  5:   Mat Top;
  6:   PetscErrorCode (*SetValues)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode);
  7:   PetscErrorCode (*SetValuesBlocked)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode);
  8: } Mat_LocalRef;

 10: /* These need to be macros because they use sizeof */
 11: #define IndexSpaceGet(buf,nrow,ncol,irowm,icolm) do {                   \
 12:     if (nrow + ncol > (PetscInt)(sizeof(buf)/sizeof(buf[0]))) {         \
 13:       PetscMalloc2(nrow,&irowm,ncol,&icolm); \
 14:     } else {                                                            \
 15:       irowm = &buf[0];                                                  \
 16:       icolm = &buf[nrow];                                               \
 17:     }                                                                   \
 18: } while (0)

 20: #define IndexSpaceRestore(buf,nrow,ncol,irowm,icolm) do {       \
 21:     if (nrow + ncol > (PetscInt)(sizeof(buf)/sizeof(buf[0]))) { \
 22:       PetscFree2(irowm,icolm);             \
 23:     }                                                           \
 24: } while (0)

 26: static void BlockIndicesExpand(PetscInt n,const PetscInt idx[],PetscInt bs,PetscInt idxm[])
 27: {
 28:   PetscInt i,j;
 29:   for (i=0; i<n; i++) {
 30:     for (j=0; j<bs; j++) {
 31:       idxm[i*bs+j] = idx[i]*bs + j;
 32:     }
 33:   }
 34: }

 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,*icolm;

 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: }

 56: static PetscErrorCode MatSetValuesBlockedLocal_LocalRef_Scalar(Mat A,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
 57: {
 58:   Mat_LocalRef   *lr = (Mat_LocalRef*)A->data;
 60:   PetscInt       bs,buf[4096],*irowm,*icolm;

 63:   MatGetBlockSize(A,&bs);
 64:   IndexSpaceGet(buf,nrow*bs,ncol*bs,irowm,icolm);
 65:   BlockIndicesExpand(nrow,irow,bs,irowm);
 66:   BlockIndicesExpand(ncol,icol,bs,icolm);
 67:   ISLocalToGlobalMappingApply(A->rmap->mapping,nrow*bs,irowm,irowm);
 68:   ISLocalToGlobalMappingApply(A->cmap->mapping,ncol*bs,icolm,icolm);
 69:   (*lr->SetValues)(lr->Top,nrow*bs,irowm,ncol*bs,icolm,y,addv);
 70:   IndexSpaceRestore(buf,nrow*bs,ncol*bs,irowm,icolm);
 71:   return(0);
 72: }

 76: static PetscErrorCode MatSetValuesLocal_LocalRef_Scalar(Mat A,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
 77: {
 78:   Mat_LocalRef   *lr = (Mat_LocalRef*)A->data;
 80:   PetscInt       buf[4096],*irowm,*icolm;

 83:   IndexSpaceGet(buf,nrow,ncol,irowm,icolm);
 84:   ISLocalToGlobalMappingApply(A->rmap->mapping,nrow,irow,irowm);
 85:   ISLocalToGlobalMappingApply(A->cmap->mapping,ncol,icol,icolm);
 86:   (*lr->SetValues)(lr->Top,nrow,irowm,ncol,icolm,y,addv);
 87:   IndexSpaceRestore(buf,nrow,ncol,irowm,icolm);
 88:   return(0);
 89: }

 93: /* Compose an IS with an ISLocalToGlobalMapping to map from IS source indices to global indices */
 94: static PetscErrorCode ISL2GCompose(IS is,ISLocalToGlobalMapping ltog,ISLocalToGlobalMapping *cltog)
 95: {
 97:   const PetscInt *idx;
 98:   PetscInt       m,*idxm;
 99:   PetscBool      isblock;

105:   PetscObjectTypeCompare((PetscObject)is,ISBLOCK,&isblock);
106:   if (isblock) {
107:     PetscInt bs,lbs;

109:     ISGetBlockSize(is,&bs);
110:     ISLocalToGlobalMappingGetBlockSize(ltog,&lbs);
111:     if (bs == lbs) {
112:       ISGetLocalSize(is,&m);
113:       m    = m/bs;
114:       ISBlockGetIndices(is,&idx);
115:       PetscMalloc1(m,&idxm);
116:       ISLocalToGlobalMappingApplyBlock(ltog,m,idx,idxm);
117:       ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)is),bs,m,idxm,PETSC_OWN_POINTER,cltog);
118:       ISBlockRestoreIndices(is,&idx);
119:       return(0);
120:     }
121:   }
122:   ISGetLocalSize(is,&m);
123:   ISGetIndices(is,&idx);
124:   PetscMalloc1(m,&idxm);
125:   if (ltog) {
126:     ISLocalToGlobalMappingApply(ltog,m,idx,idxm);
127:   } else {
128:     PetscMemcpy(idxm,idx,m*sizeof(PetscInt));
129:   }
130:   ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)is),1,m,idxm,PETSC_OWN_POINTER,cltog);
131:   ISRestoreIndices(is,&idx);
132:   return(0);
133: }

137: static PetscErrorCode ISL2GComposeBlock(IS is,ISLocalToGlobalMapping ltog,ISLocalToGlobalMapping *cltog)
138: {
140:   const PetscInt *idx;
141:   PetscInt       m,*idxm,bs;

147:   ISBlockGetLocalSize(is,&m);
148:   ISBlockGetIndices(is,&idx);
149:   ISLocalToGlobalMappingGetBlockSize(ltog,&bs);
150:   PetscMalloc1(m,&idxm);
151:   if (ltog) {
152:     ISLocalToGlobalMappingApplyBlock(ltog,m,idx,idxm);
153:   } else {
154:     PetscMemcpy(idxm,idx,m*sizeof(PetscInt));
155:   }
156:   ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)is),bs,m,idxm,PETSC_OWN_POINTER,cltog);
157:   ISBlockRestoreIndices(is,&idx);
158:   return(0);
159: }

163: static PetscErrorCode MatDestroy_LocalRef(Mat B)
164: {

168:   PetscFree(B->data);
169:   return(0);
170: }


175: /*@
176:    MatCreateLocalRef - Gets a logical reference to a local submatrix, for use in assembly

178:    Not Collective

180:    Input Arguments:
181: + A - Full matrix, generally parallel
182: . isrow - Local index set for the rows
183: - iscol - Local index set for the columns

185:    Output Arguments:
186: . newmat - New serial Mat

188:    Level: developer

190:    Notes:
191:    Most will use MatGetLocalSubMatrix() which returns a real matrix corresponding to the local
192:    block if it available, such as with matrix formats that store these blocks separately.

194:    The new matrix forwards MatSetValuesLocal() and MatSetValuesBlockedLocal() to the global system.
195:    In general, it does not define MatMult() or any other functions.  Local submatrices can be nested.

197: .seealso: MatSetValuesLocal(), MatSetValuesBlockedLocal(), MatGetLocalSubMatrix(), MatCreateSubMatrix()
198: @*/
199: PetscErrorCode  MatCreateLocalRef(Mat A,IS isrow,IS iscol,Mat *newmat)
200: {
202:   Mat_LocalRef   *lr;
203:   Mat            B;
204:   PetscInt       m,n;
205:   PetscBool      islr;

212:   *newmat = 0;

214:   MatCreate(PETSC_COMM_SELF,&B);
215:   ISGetLocalSize(isrow,&m);
216:   ISGetLocalSize(iscol,&n);
217:   MatSetSizes(B,m,n,m,n);
218:   PetscObjectChangeTypeName((PetscObject)B,MATLOCALREF);
219:   MatSetUp(B);

221:   B->ops->destroy = MatDestroy_LocalRef;

223:   PetscNewLog(B,&lr);
224:   B->data = (void*)lr;

226:   PetscObjectTypeCompare((PetscObject)A,MATLOCALREF,&islr);
227:   if (islr) {
228:     Mat_LocalRef *alr = (Mat_LocalRef*)A->data;
229:     lr->Top = alr->Top;
230:   } else {
231:     /* This does not increase the reference count because MatLocalRef is not allowed to live longer than its parent */
232:     lr->Top = A;
233:   }
234:   {
235:     ISLocalToGlobalMapping rltog,cltog;
236:     PetscInt               abs,rbs,cbs;

238:     /* We will translate directly to global indices for the top level */
239:     lr->SetValues        = MatSetValues;
240:     lr->SetValuesBlocked = MatSetValuesBlocked;

242:     B->ops->setvalueslocal = MatSetValuesLocal_LocalRef_Scalar;

244:     ISL2GCompose(isrow,A->rmap->mapping,&rltog);
245:     if (isrow == iscol && A->rmap->mapping == A->cmap->mapping) {
246:       PetscObjectReference((PetscObject)rltog);
247:       cltog = rltog;
248:     } else {
249:       ISL2GCompose(iscol,A->cmap->mapping,&cltog);
250:     }
251:     MatSetLocalToGlobalMapping(B,rltog,cltog);
252:     ISLocalToGlobalMappingDestroy(&rltog);
253:     ISLocalToGlobalMappingDestroy(&cltog);

255:     MatGetBlockSize(A,&abs);
256:     ISGetBlockSize(isrow,&rbs);
257:     ISGetBlockSize(iscol,&cbs);
258:     if (rbs == cbs) {           /* submatrix has block structure, so user can insert values with blocked interface */
259:       PetscLayoutSetBlockSize(B->rmap,rbs);
260:       PetscLayoutSetBlockSize(B->cmap,cbs);
261:       if (abs != rbs || abs == 1) {
262:         /* Top-level matrix has different block size, so we have to call its scalar insertion interface */
263:         B->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_LocalRef_Scalar;
264:       } else {
265:         /* Block sizes match so we can forward values to the top level using the block interface */
266:         B->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_LocalRef_Block;

268:         ISL2GComposeBlock(isrow,A->rmap->mapping,&rltog);
269:         if (isrow == iscol && A->rmap->mapping == A->cmap->mapping) {
270:            PetscObjectReference((PetscObject)rltog);
271:           cltog = rltog;
272:         } else {
273:           ISL2GComposeBlock(iscol,A->cmap->mapping,&cltog);
274:         }
275:         MatSetLocalToGlobalMapping(B,rltog,cltog);
276:         ISLocalToGlobalMappingDestroy(&rltog);
277:         ISLocalToGlobalMappingDestroy(&cltog);
278:       }
279:     }
280:   }
281:   *newmat = B;
282:   return(0);
283: }