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;
41: PetscInt buf[4096],*irowm=NULL,*icolm; /* suppress maybe-uninitialized warning */
43: if (!nrow || !ncol) return 0;
44: IndexSpaceGet(buf,nrow,ncol,irowm,icolm);
45: ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,nrow,irow,irowm);
46: ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,ncol,icol,icolm);
47: (*lr->SetValuesBlocked)(lr->Top,nrow,irowm,ncol,icolm,y,addv);
48: IndexSpaceRestore(buf,nrow,ncol,irowm,icolm);
49: return 0;
50: }
52: static PetscErrorCode MatSetValuesBlockedLocal_LocalRef_Scalar(Mat A,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
53: {
54: Mat_LocalRef *lr = (Mat_LocalRef*)A->data;
55: PetscInt rbs,cbs,buf[4096],*irowm,*icolm;
57: MatGetBlockSizes(A,&rbs,&cbs);
58: IndexSpaceGet(buf,nrow*rbs,ncol*cbs,irowm,icolm);
59: BlockIndicesExpand(nrow,irow,rbs,irowm);
60: BlockIndicesExpand(ncol,icol,cbs,icolm);
61: ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,nrow*rbs,irowm,irowm);
62: ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,ncol*cbs,icolm,icolm);
63: (*lr->SetValues)(lr->Top,nrow*rbs,irowm,ncol*cbs,icolm,y,addv);
64: IndexSpaceRestore(buf,nrow*rbs,ncol*cbs,irowm,icolm);
65: return 0;
66: }
68: static PetscErrorCode MatSetValuesLocal_LocalRef_Scalar(Mat A,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
69: {
70: Mat_LocalRef *lr = (Mat_LocalRef*)A->data;
71: PetscInt buf[4096],*irowm,*icolm;
73: IndexSpaceGet(buf,nrow,ncol,irowm,icolm);
74: /* If the row IS defining this submatrix was an ISBLOCK, then the unblocked LGMapApply is the right one to use. If
75: * instead it was (say) an ISSTRIDE with a block size > 1, then we need to use LGMapApplyBlock */
76: if (lr->rowisblock) {
77: ISLocalToGlobalMappingApply(A->rmap->mapping,nrow,irow,irowm);
78: } else {
79: ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,nrow,irow,irowm);
80: }
81: /* As above, but for the column IS. */
82: if (lr->colisblock) {
83: ISLocalToGlobalMappingApply(A->cmap->mapping,ncol,icol,icolm);
84: } else {
85: ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,ncol,icol,icolm);
86: }
87: (*lr->SetValues)(lr->Top,nrow,irowm,ncol,icolm,y,addv);
88: IndexSpaceRestore(buf,nrow,ncol,irowm,icolm);
89: return 0;
90: }
92: /* Compose an IS with an ISLocalToGlobalMapping to map from IS source indices to global indices */
93: static PetscErrorCode ISL2GCompose(IS is,ISLocalToGlobalMapping ltog,ISLocalToGlobalMapping *cltog)
94: {
95: const PetscInt *idx;
96: PetscInt m,*idxm;
97: PetscInt bs;
98: PetscBool isblock;
103: PetscObjectTypeCompare((PetscObject)is,ISBLOCK,&isblock);
104: if (isblock) {
105: PetscInt lbs;
107: ISGetBlockSize(is,&bs);
108: ISLocalToGlobalMappingGetBlockSize(ltog,&lbs);
109: if (bs == lbs) {
110: ISGetLocalSize(is,&m);
111: m = m/bs;
112: ISBlockGetIndices(is,&idx);
113: PetscMalloc1(m,&idxm);
114: ISLocalToGlobalMappingApplyBlock(ltog,m,idx,idxm);
115: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)is),bs,m,idxm,PETSC_OWN_POINTER,cltog);
116: ISBlockRestoreIndices(is,&idx);
117: return 0;
118: }
119: }
120: ISGetLocalSize(is,&m);
121: ISGetIndices(is,&idx);
122: ISGetBlockSize(is,&bs);
123: PetscMalloc1(m,&idxm);
124: if (ltog) {
125: ISLocalToGlobalMappingApply(ltog,m,idx,idxm);
126: } else {
127: PetscArraycpy(idxm,idx,m);
128: }
129: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)is),bs,m,idxm,PETSC_OWN_POINTER,cltog);
130: ISRestoreIndices(is,&idx);
131: return 0;
132: }
134: static PetscErrorCode ISL2GComposeBlock(IS is,ISLocalToGlobalMapping ltog,ISLocalToGlobalMapping *cltog)
135: {
136: const PetscInt *idx;
137: PetscInt m,*idxm,bs;
142: ISBlockGetLocalSize(is,&m);
143: ISBlockGetIndices(is,&idx);
144: ISLocalToGlobalMappingGetBlockSize(ltog,&bs);
145: PetscMalloc1(m,&idxm);
146: if (ltog) {
147: ISLocalToGlobalMappingApplyBlock(ltog,m,idx,idxm);
148: } else {
149: PetscArraycpy(idxm,idx,m);
150: }
151: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)is),bs,m,idxm,PETSC_OWN_POINTER,cltog);
152: ISBlockRestoreIndices(is,&idx);
153: return 0;
154: }
156: static PetscErrorCode MatDestroy_LocalRef(Mat B)
157: {
158: PetscFree(B->data);
159: return 0;
160: }
162: /*@
163: MatCreateLocalRef - Gets a logical reference to a local submatrix, for use in assembly
165: Not Collective
167: Input Parameters:
168: + A - Full matrix, generally parallel
169: . isrow - Local index set for the rows
170: - iscol - Local index set for the columns
172: Output Parameter:
173: . newmat - New serial Mat
175: Level: developer
177: Notes:
178: Most will use MatGetLocalSubMatrix() which returns a real matrix corresponding to the local
179: block if it available, such as with matrix formats that store these blocks separately.
181: The new matrix forwards MatSetValuesLocal() and MatSetValuesBlockedLocal() to the global system.
182: In general, it does not define MatMult() or any other functions. Local submatrices can be nested.
184: .seealso: MatSetValuesLocal(), MatSetValuesBlockedLocal(), MatGetLocalSubMatrix(), MatCreateSubMatrix()
185: @*/
186: PetscErrorCode MatCreateLocalRef(Mat A,IS isrow,IS iscol,Mat *newmat)
187: {
188: Mat_LocalRef *lr;
189: Mat B;
190: PetscInt m,n;
191: PetscBool islr;
198: *newmat = NULL;
200: MatCreate(PETSC_COMM_SELF,&B);
201: ISGetLocalSize(isrow,&m);
202: ISGetLocalSize(iscol,&n);
203: MatSetSizes(B,m,n,m,n);
204: PetscObjectChangeTypeName((PetscObject)B,MATLOCALREF);
205: MatSetUp(B);
207: B->ops->destroy = MatDestroy_LocalRef;
209: PetscNewLog(B,&lr);
210: B->data = (void*)lr;
212: PetscObjectTypeCompare((PetscObject)A,MATLOCALREF,&islr);
213: if (islr) {
214: Mat_LocalRef *alr = (Mat_LocalRef*)A->data;
215: lr->Top = alr->Top;
216: } else {
217: /* This does not increase the reference count because MatLocalRef is not allowed to live longer than its parent */
218: lr->Top = A;
219: }
220: {
221: ISLocalToGlobalMapping rltog,cltog;
222: PetscInt arbs,acbs,rbs,cbs;
224: /* We will translate directly to global indices for the top level */
225: lr->SetValues = MatSetValues;
226: lr->SetValuesBlocked = MatSetValuesBlocked;
228: B->ops->setvalueslocal = MatSetValuesLocal_LocalRef_Scalar;
230: ISL2GCompose(isrow,A->rmap->mapping,&rltog);
231: if (isrow == iscol && A->rmap->mapping == A->cmap->mapping) {
232: PetscObjectReference((PetscObject)rltog);
233: cltog = rltog;
234: } else {
235: ISL2GCompose(iscol,A->cmap->mapping,&cltog);
236: }
237: /* Remember if the ISes we used to pull out the submatrix are of type ISBLOCK. Will be used later in
238: * MatSetValuesLocal_LocalRef_Scalar. */
239: PetscObjectTypeCompare((PetscObject)isrow,ISBLOCK,&lr->rowisblock);
240: PetscObjectTypeCompare((PetscObject)iscol,ISBLOCK,&lr->colisblock);
241: MatSetLocalToGlobalMapping(B,rltog,cltog);
242: ISLocalToGlobalMappingDestroy(&rltog);
243: ISLocalToGlobalMappingDestroy(&cltog);
245: MatGetBlockSizes(A,&arbs,&acbs);
246: ISGetBlockSize(isrow,&rbs);
247: ISGetBlockSize(iscol,&cbs);
248: /* Always support block interface insertion on submatrix */
249: PetscLayoutSetBlockSize(B->rmap,rbs);
250: PetscLayoutSetBlockSize(B->cmap,cbs);
251: if (arbs != rbs || acbs != cbs || (arbs == 1 && acbs == 1)) {
252: /* Top-level matrix has different block size, so we have to call its scalar insertion interface */
253: B->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_LocalRef_Scalar;
254: } else {
255: /* Block sizes match so we can forward values to the top level using the block interface */
256: B->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_LocalRef_Block;
258: ISL2GComposeBlock(isrow,A->rmap->mapping,&rltog);
259: if (isrow == iscol && A->rmap->mapping == A->cmap->mapping) {
260: PetscObjectReference((PetscObject)rltog);
261: cltog = rltog;
262: } else {
263: ISL2GComposeBlock(iscol,A->cmap->mapping,&cltog);
264: }
265: MatSetLocalToGlobalMapping(B,rltog,cltog);
266: ISLocalToGlobalMappingDestroy(&rltog);
267: ISLocalToGlobalMappingDestroy(&cltog);
268: }
269: }
270: *newmat = B;
271: return 0;
272: }