Actual source code: mlocalref.c
petsc-3.6.1 2015-08-06
2: #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
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: }
40: static PetscErrorCode MatSetValuesBlockedLocal_LocalRef_Block(Mat A,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
41: {
42: Mat_LocalRef *lr = (Mat_LocalRef*)A->data;
44: PetscInt buf[4096],*irowm,*icolm;
47: if (!nrow || !ncol) return(0);
48: IndexSpaceGet(buf,nrow,ncol,irowm,icolm);
49: ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,nrow,irow,irowm);
50: ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,ncol,icol,icolm);
51: (*lr->SetValuesBlocked)(lr->Top,nrow,irowm,ncol,icolm,y,addv);
52: IndexSpaceRestore(buf,nrow,ncol,irowm,icolm);
53: return(0);
54: }
58: static PetscErrorCode MatSetValuesBlockedLocal_LocalRef_Scalar(Mat A,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
59: {
60: Mat_LocalRef *lr = (Mat_LocalRef*)A->data;
62: PetscInt rbs,cbs,buf[4096],*irowm,*icolm;
65: MatGetBlockSizes(A,&rbs,&cbs);
66: IndexSpaceGet(buf,nrow*rbs,ncol*cbs,irowm,icolm);
67: BlockIndicesExpand(nrow,irow,rbs,irowm);
68: BlockIndicesExpand(ncol,icol,cbs,icolm);
69: ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,nrow*rbs,irowm,irowm);
70: ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,ncol*cbs,icolm,icolm);
71: (*lr->SetValues)(lr->Top,nrow*rbs,irowm,ncol*cbs,icolm,y,addv);
72: IndexSpaceRestore(buf,nrow*rbs,ncol*cbs,irowm,icolm);
73: return(0);
74: }
78: static PetscErrorCode MatSetValuesLocal_LocalRef_Scalar(Mat A,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
79: {
80: Mat_LocalRef *lr = (Mat_LocalRef*)A->data;
82: PetscInt buf[4096],*irowm,*icolm;
85: IndexSpaceGet(buf,nrow,ncol,irowm,icolm);
86: /* If the row IS defining this submatrix was an ISBLOCK, then the unblocked LGMapApply is the right one to use. If
87: * instead it was (say) an ISSTRIDE with a block size > 1, then we need to use LGMapApplyBlock */
88: if (lr->rowisblock) {
89: ISLocalToGlobalMappingApply(A->rmap->mapping,nrow,irow,irowm);
90: } else {
91: ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,nrow,irow,irowm);
92: }
93: /* As above, but for the column IS. */
94: if (lr->colisblock) {
95: ISLocalToGlobalMappingApply(A->cmap->mapping,ncol,icol,icolm);
96: } else {
97: ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,ncol,icol,icolm);
98: }
99: (*lr->SetValues)(lr->Top,nrow,irowm,ncol,icolm,y,addv);
100: IndexSpaceRestore(buf,nrow,ncol,irowm,icolm);
101: return(0);
102: }
106: /* Compose an IS with an ISLocalToGlobalMapping to map from IS source indices to global indices */
107: static PetscErrorCode ISL2GCompose(IS is,ISLocalToGlobalMapping ltog,ISLocalToGlobalMapping *cltog)
108: {
110: const PetscInt *idx;
111: PetscInt m,*idxm;
112: PetscInt bs;
113: PetscBool isblock;
119: PetscObjectTypeCompare((PetscObject)is,ISBLOCK,&isblock);
120: if (isblock) {
121: PetscInt lbs;
123: ISGetBlockSize(is,&bs);
124: ISLocalToGlobalMappingGetBlockSize(ltog,&lbs);
125: if (bs == lbs) {
126: ISGetLocalSize(is,&m);
127: m = m/bs;
128: ISBlockGetIndices(is,&idx);
129: PetscMalloc1(m,&idxm);
130: ISLocalToGlobalMappingApplyBlock(ltog,m,idx,idxm);
131: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)is),bs,m,idxm,PETSC_OWN_POINTER,cltog);
132: ISBlockRestoreIndices(is,&idx);
133: return(0);
134: }
135: }
136: ISGetLocalSize(is,&m);
137: ISGetIndices(is,&idx);
138: ISGetBlockSize(is,&bs);
139: PetscMalloc1(m,&idxm);
140: if (ltog) {
141: ISLocalToGlobalMappingApply(ltog,m,idx,idxm);
142: } else {
143: PetscMemcpy(idxm,idx,m*sizeof(PetscInt));
144: }
145: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)is),bs,m,idxm,PETSC_OWN_POINTER,cltog);
146: ISRestoreIndices(is,&idx);
147: return(0);
148: }
152: static PetscErrorCode ISL2GComposeBlock(IS is,ISLocalToGlobalMapping ltog,ISLocalToGlobalMapping *cltog)
153: {
155: const PetscInt *idx;
156: PetscInt m,*idxm,bs;
162: ISBlockGetLocalSize(is,&m);
163: ISBlockGetIndices(is,&idx);
164: ISLocalToGlobalMappingGetBlockSize(ltog,&bs);
165: PetscMalloc1(m,&idxm);
166: if (ltog) {
167: ISLocalToGlobalMappingApplyBlock(ltog,m,idx,idxm);
168: } else {
169: PetscMemcpy(idxm,idx,m*sizeof(PetscInt));
170: }
171: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)is),bs,m,idxm,PETSC_OWN_POINTER,cltog);
172: ISBlockRestoreIndices(is,&idx);
173: return(0);
174: }
178: static PetscErrorCode MatDestroy_LocalRef(Mat B)
179: {
183: PetscFree(B->data);
184: return(0);
185: }
190: /*@
191: MatCreateLocalRef - Gets a logical reference to a local submatrix, for use in assembly
193: Not Collective
195: Input Arguments:
196: + A - Full matrix, generally parallel
197: . isrow - Local index set for the rows
198: - iscol - Local index set for the columns
200: Output Arguments:
201: . newmat - New serial Mat
203: Level: developer
205: Notes:
206: Most will use MatGetLocalSubMatrix() which returns a real matrix corresponding to the local
207: block if it available, such as with matrix formats that store these blocks separately.
209: The new matrix forwards MatSetValuesLocal() and MatSetValuesBlockedLocal() to the global system.
210: In general, it does not define MatMult() or any other functions. Local submatrices can be nested.
212: .seealso: MatSetValuesLocal(), MatSetValuesBlockedLocal(), MatGetLocalSubMatrix(), MatCreateSubMatrix()
213: @*/
214: PetscErrorCode MatCreateLocalRef(Mat A,IS isrow,IS iscol,Mat *newmat)
215: {
217: Mat_LocalRef *lr;
218: Mat B;
219: PetscInt m,n;
220: PetscBool islr;
227: *newmat = 0;
229: MatCreate(PETSC_COMM_SELF,&B);
230: ISGetLocalSize(isrow,&m);
231: ISGetLocalSize(iscol,&n);
232: MatSetSizes(B,m,n,m,n);
233: PetscObjectChangeTypeName((PetscObject)B,MATLOCALREF);
234: MatSetUp(B);
236: B->ops->destroy = MatDestroy_LocalRef;
238: PetscNewLog(B,&lr);
239: B->data = (void*)lr;
241: PetscObjectTypeCompare((PetscObject)A,MATLOCALREF,&islr);
242: if (islr) {
243: Mat_LocalRef *alr = (Mat_LocalRef*)A->data;
244: lr->Top = alr->Top;
245: } else {
246: /* This does not increase the reference count because MatLocalRef is not allowed to live longer than its parent */
247: lr->Top = A;
248: }
249: {
250: ISLocalToGlobalMapping rltog,cltog;
251: PetscInt arbs,acbs,rbs,cbs;
253: /* We will translate directly to global indices for the top level */
254: lr->SetValues = MatSetValues;
255: lr->SetValuesBlocked = MatSetValuesBlocked;
257: B->ops->setvalueslocal = MatSetValuesLocal_LocalRef_Scalar;
259: ISL2GCompose(isrow,A->rmap->mapping,&rltog);
260: if (isrow == iscol && A->rmap->mapping == A->cmap->mapping) {
261: PetscObjectReference((PetscObject)rltog);
262: cltog = rltog;
263: } else {
264: ISL2GCompose(iscol,A->cmap->mapping,&cltog);
265: }
266: /* Remember if the ISes we used to pull out the submatrix are of type ISBLOCK. Will be used later in
267: * MatSetValuesLocal_LocalRef_Scalar. */
268: PetscObjectTypeCompare((PetscObject)isrow,ISBLOCK,&lr->rowisblock);
269: PetscObjectTypeCompare((PetscObject)iscol,ISBLOCK,&lr->colisblock);
270: MatSetLocalToGlobalMapping(B,rltog,cltog);
271: ISLocalToGlobalMappingDestroy(&rltog);
272: ISLocalToGlobalMappingDestroy(&cltog);
274: MatGetBlockSizes(A,&arbs,&acbs);
275: ISGetBlockSize(isrow,&rbs);
276: ISGetBlockSize(iscol,&cbs);
277: /* Always support block interface insertion on submatrix */
278: PetscLayoutSetBlockSize(B->rmap,rbs);
279: PetscLayoutSetBlockSize(B->cmap,cbs);
280: if (arbs != rbs || acbs != cbs || (arbs == 1 && acbs == 1)) {
281: /* Top-level matrix has different block size, so we have to call its scalar insertion interface */
282: B->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_LocalRef_Scalar;
283: } else {
284: /* Block sizes match so we can forward values to the top level using the block interface */
285: B->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_LocalRef_Block;
287: ISL2GComposeBlock(isrow,A->rmap->mapping,&rltog);
288: if (isrow == iscol && A->rmap->mapping == A->cmap->mapping) {
289: PetscObjectReference((PetscObject)rltog);
290: cltog = rltog;
291: } else {
292: ISL2GComposeBlock(iscol,A->cmap->mapping,&cltog);
293: }
294: MatSetLocalToGlobalMapping(B,rltog,cltog);
295: ISLocalToGlobalMappingDestroy(&rltog);
296: ISLocalToGlobalMappingDestroy(&cltog);
297: }
298: }
299: *newmat = B;
300: return(0);
301: }