Actual source code: mlocalref.c
petsc-3.5.4 2015-05-23
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: }