Actual source code: mlocalref.c
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: }
175: /*@
176: MatCreateLocalRef - Gets a logical reference to a local submatrix, for use in assembly
178: Not Collective
180: Input Parameters:
181: + A - Full matrix, generally parallel
182: . isrow - Local index set for the rows
183: - iscol - Local index set for the columns
185: Output Parameter:
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: if (!A->rmap->mapping) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Matrix must have local to global mapping provided before this call");
213: *newmat = NULL;
215: MatCreate(PETSC_COMM_SELF,&B);
216: ISGetLocalSize(isrow,&m);
217: ISGetLocalSize(iscol,&n);
218: MatSetSizes(B,m,n,m,n);
219: PetscObjectChangeTypeName((PetscObject)B,MATLOCALREF);
220: MatSetUp(B);
222: B->ops->destroy = MatDestroy_LocalRef;
224: PetscNewLog(B,&lr);
225: B->data = (void*)lr;
227: PetscObjectTypeCompare((PetscObject)A,MATLOCALREF,&islr);
228: if (islr) {
229: Mat_LocalRef *alr = (Mat_LocalRef*)A->data;
230: lr->Top = alr->Top;
231: } else {
232: /* This does not increase the reference count because MatLocalRef is not allowed to live longer than its parent */
233: lr->Top = A;
234: }
235: {
236: ISLocalToGlobalMapping rltog,cltog;
237: PetscInt arbs,acbs,rbs,cbs;
239: /* We will translate directly to global indices for the top level */
240: lr->SetValues = MatSetValues;
241: lr->SetValuesBlocked = MatSetValuesBlocked;
243: B->ops->setvalueslocal = MatSetValuesLocal_LocalRef_Scalar;
245: ISL2GCompose(isrow,A->rmap->mapping,&rltog);
246: if (isrow == iscol && A->rmap->mapping == A->cmap->mapping) {
247: PetscObjectReference((PetscObject)rltog);
248: cltog = rltog;
249: } else {
250: ISL2GCompose(iscol,A->cmap->mapping,&cltog);
251: }
252: /* Remember if the ISes we used to pull out the submatrix are of type ISBLOCK. Will be used later in
253: * MatSetValuesLocal_LocalRef_Scalar. */
254: PetscObjectTypeCompare((PetscObject)isrow,ISBLOCK,&lr->rowisblock);
255: PetscObjectTypeCompare((PetscObject)iscol,ISBLOCK,&lr->colisblock);
256: MatSetLocalToGlobalMapping(B,rltog,cltog);
257: ISLocalToGlobalMappingDestroy(&rltog);
258: ISLocalToGlobalMappingDestroy(&cltog);
260: MatGetBlockSizes(A,&arbs,&acbs);
261: ISGetBlockSize(isrow,&rbs);
262: ISGetBlockSize(iscol,&cbs);
263: /* Always support block interface insertion on submatrix */
264: PetscLayoutSetBlockSize(B->rmap,rbs);
265: PetscLayoutSetBlockSize(B->cmap,cbs);
266: if (arbs != rbs || acbs != cbs || (arbs == 1 && acbs == 1)) {
267: /* Top-level matrix has different block size, so we have to call its scalar insertion interface */
268: B->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_LocalRef_Scalar;
269: } else {
270: /* Block sizes match so we can forward values to the top level using the block interface */
271: B->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_LocalRef_Block;
273: ISL2GComposeBlock(isrow,A->rmap->mapping,&rltog);
274: if (isrow == iscol && A->rmap->mapping == A->cmap->mapping) {
275: PetscObjectReference((PetscObject)rltog);
276: cltog = rltog;
277: } else {
278: ISL2GComposeBlock(iscol,A->cmap->mapping,&cltog);
279: }
280: MatSetLocalToGlobalMapping(B,rltog,cltog);
281: ISLocalToGlobalMappingDestroy(&rltog);
282: ISLocalToGlobalMappingDestroy(&cltog);
283: }
284: }
285: *newmat = B;
286: return(0);
287: }