Actual source code: mlocalref.c
petsc-3.4.5 2014-06-29
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,PetscInt,&irowm,ncol,PetscInt,&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: ISLocalToGlobalMappingApply(A->rmap->bmapping,nrow,irow,irowm);
48: ISLocalToGlobalMappingApply(A->cmap->bmapping,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 = A->rmap->bs,buf[4096],*irowm,*icolm;
63: IndexSpaceGet(buf,nrow*bs,ncol*bs,irowm,icolm);
64: BlockIndicesExpand(nrow,irow,bs,irowm);
65: BlockIndicesExpand(ncol,icol,bs,icolm);
66: ISLocalToGlobalMappingApply(A->rmap->mapping,nrow*bs,irowm,irowm);
67: ISLocalToGlobalMappingApply(A->cmap->mapping,ncol*bs,icolm,icolm);
68: (*lr->SetValues)(lr->Top,nrow*bs,irowm,ncol*bs,icolm,y,addv);
69: IndexSpaceRestore(buf,nrow*bs,ncol*bs,irowm,icolm);
70: return(0);
71: }
75: static PetscErrorCode MatSetValuesLocal_LocalRef_Scalar(Mat A,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
76: {
77: Mat_LocalRef *lr = (Mat_LocalRef*)A->data;
79: PetscInt buf[4096],*irowm,*icolm;
82: IndexSpaceGet(buf,nrow,ncol,irowm,icolm);
83: ISLocalToGlobalMappingApply(A->rmap->mapping,nrow,irow,irowm);
84: ISLocalToGlobalMappingApply(A->cmap->mapping,ncol,icol,icolm);
85: (*lr->SetValues)(lr->Top,nrow,irowm,ncol,icolm,y,addv);
86: IndexSpaceRestore(buf,nrow,ncol,irowm,icolm);
87: return(0);
88: }
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: {
96: const PetscInt *idx;
97: PetscInt m,*idxm;
103: ISGetLocalSize(is,&m);
104: ISGetIndices(is,&idx);
105: PetscMalloc(m*sizeof(PetscInt),&idxm);
106: if (ltog) {
107: ISLocalToGlobalMappingApply(ltog,m,idx,idxm);
108: } else {
109: PetscMemcpy(idxm,idx,m*sizeof(PetscInt));
110: }
111: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)is),m,idxm,PETSC_OWN_POINTER,cltog);
112: ISRestoreIndices(is,&idx);
113: return(0);
114: }
118: static PetscErrorCode ISL2GComposeBlock(IS is,ISLocalToGlobalMapping ltog,ISLocalToGlobalMapping *cltog)
119: {
121: const PetscInt *idx;
122: PetscInt m,*idxm;
128: ISBlockGetLocalSize(is,&m);
129: ISBlockGetIndices(is,&idx);
130: PetscMalloc(m*sizeof(PetscInt),&idxm);
131: if (ltog) {
132: ISLocalToGlobalMappingApply(ltog,m,idx,idxm);
133: } else {
134: PetscMemcpy(idxm,idx,m*sizeof(PetscInt));
135: }
136: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)is),m,idxm,PETSC_OWN_POINTER,cltog);
137: ISBlockRestoreIndices(is,&idx);
138: return(0);
139: }
143: static PetscErrorCode MatDestroy_LocalRef(Mat B)
144: {
148: PetscFree(B->data);
149: return(0);
150: }
155: /*@
156: MatCreateLocalRef - Gets a logical reference to a local submatrix, for use in assembly
158: Not Collective
160: Input Arguments:
161: + A - Full matrix, generally parallel
162: . isrow - Local index set for the rows
163: - iscol - Local index set for the columns
165: Output Arguments:
166: . newmat - New serial Mat
168: Level: developer
170: Notes:
171: Most will use MatGetLocalSubMatrix() which returns a real matrix corresponding to the local
172: block if it available, such as with matrix formats that store these blocks separately.
174: The new matrix forwards MatSetValuesLocal() and MatSetValuesBlockedLocal() to the global system.
175: In general, it does not define MatMult() or any other functions. Local submatrices can be nested.
177: .seealso: MatSetValuesLocal(), MatSetValuesBlockedLocal(), MatGetLocalSubMatrix(), MatCreateSubMatrix()
178: @*/
179: PetscErrorCode MatCreateLocalRef(Mat A,IS isrow,IS iscol,Mat *newmat)
180: {
182: Mat_LocalRef *lr;
183: Mat B;
184: PetscInt m,n;
185: PetscBool islr;
192: *newmat = 0;
194: MatCreate(PETSC_COMM_SELF,&B);
195: ISGetLocalSize(isrow,&m);
196: ISGetLocalSize(iscol,&n);
197: MatSetSizes(B,m,n,m,n);
198: PetscObjectChangeTypeName((PetscObject)B,MATLOCALREF);
199: MatSetUp(B);
201: B->ops->destroy = MatDestroy_LocalRef;
203: PetscNewLog(B,Mat_LocalRef,&lr);
204: B->data = (void*)lr;
206: PetscObjectTypeCompare((PetscObject)A,MATLOCALREF,&islr);
207: if (islr) {
208: Mat_LocalRef *alr = (Mat_LocalRef*)A->data;
209: lr->Top = alr->Top;
210: } else {
211: /* This does not increase the reference count because MatLocalRef is not allowed to live longer than its parent */
212: lr->Top = A;
213: }
214: {
215: ISLocalToGlobalMapping rltog,cltog;
216: PetscInt abs,rbs,cbs;
218: /* We will translate directly to global indices for the top level */
219: lr->SetValues = MatSetValues;
220: lr->SetValuesBlocked = MatSetValuesBlocked;
222: B->ops->setvalueslocal = MatSetValuesLocal_LocalRef_Scalar;
224: ISL2GCompose(isrow,A->rmap->mapping,&rltog);
225: if (isrow == iscol && A->rmap->mapping == A->cmap->mapping) {
226: PetscObjectReference((PetscObject)rltog);
227: cltog = rltog;
228: } else {
229: ISL2GCompose(iscol,A->cmap->mapping,&cltog);
230: }
231: MatSetLocalToGlobalMapping(B,rltog,cltog);
232: ISLocalToGlobalMappingDestroy(&rltog);
233: ISLocalToGlobalMappingDestroy(&cltog);
235: MatGetBlockSize(A,&abs);
236: ISGetBlockSize(isrow,&rbs);
237: ISGetBlockSize(iscol,&cbs);
238: if (rbs == cbs) { /* submatrix has block structure, so user can insert values with blocked interface */
239: PetscLayoutSetBlockSize(B->rmap,rbs);
240: PetscLayoutSetBlockSize(B->cmap,cbs);
241: if (abs != rbs || abs == 1) {
242: /* Top-level matrix has different block size, so we have to call its scalar insertion interface */
243: B->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_LocalRef_Scalar;
244: } else {
245: /* Block sizes match so we can forward values to the top level using the block interface */
246: B->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_LocalRef_Block;
248: ISL2GComposeBlock(isrow,A->rmap->bmapping,&rltog);
249: if (isrow == iscol && A->rmap->bmapping == A->cmap->bmapping) {
250: PetscObjectReference((PetscObject)rltog);
251: cltog = rltog;
252: } else {
253: ISL2GComposeBlock(iscol,A->cmap->bmapping,&cltog);
254: }
255: MatSetLocalToGlobalMappingBlock(B,rltog,cltog);
256: ISLocalToGlobalMappingDestroy(&rltog);
257: ISLocalToGlobalMappingDestroy(&cltog);
258: }
259: }
260: }
261: *newmat = B;
262: return(0);
263: }