Actual source code: mlocalref.c
petsc-3.3-p7 2013-05-11
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: #if defined(PETSC_USE_DEBUG)
106: {
107: PetscInt i;
108: for (i=0; i<m; i++) {
109: if (idx[i] < 0 || ltog->n <= idx[i]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"is[%D] = %D is not in the local range [0:%D]",i,idx[i],ltog->n);
110: }
111: }
112: #endif
113: PetscMalloc(m*sizeof(PetscInt),&idxm);
114: if (ltog) {
115: ISLocalToGlobalMappingApply(ltog,m,idx,idxm);
116: } else {
117: PetscMemcpy(idxm,idx,m*sizeof(PetscInt));
118: }
119: ISLocalToGlobalMappingCreate(((PetscObject)is)->comm,m,idxm,PETSC_OWN_POINTER,cltog);
120: ISRestoreIndices(is,&idx);
121: return(0);
122: }
126: static PetscErrorCode ISL2GComposeBlock(IS is,ISLocalToGlobalMapping ltog,ISLocalToGlobalMapping *cltog)
127: {
129: const PetscInt *idx;
130: PetscInt m,*idxm;
136: ISBlockGetLocalSize(is,&m);
137: ISBlockGetIndices(is,&idx);
138: #if defined(PETSC_USE_DEBUG)
139: {
140: PetscInt i;
141: for (i=0; i<m; i++) {
142: if (idx[i] < 0 || ltog->n <= idx[i]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"is[%D] = %D is not in the local range [0:%D]",i,idx[i],ltog->n);
143: }
144: }
145: #endif
146: PetscMalloc(m*sizeof(PetscInt),&idxm);
147: if (ltog) {
148: ISLocalToGlobalMappingApply(ltog,m,idx,idxm);
149: } else {
150: PetscMemcpy(idxm,idx,m*sizeof(PetscInt));
151: }
152: ISLocalToGlobalMappingCreate(((PetscObject)is)->comm,m,idxm,PETSC_OWN_POINTER,cltog);
153: ISBlockRestoreIndices(is,&idx);
154: return(0);
155: }
159: static PetscErrorCode MatDestroy_LocalRef(Mat B)
160: {
164: PetscFree(B->data);
165: return(0);
166: }
171: /*@
172: MatCreateLocalRef - Gets a logical reference to a local submatrix, for use in assembly
174: Not Collective
176: Input Arguments:
177: + A - Full matrix, generally parallel
178: . isrow - Local index set for the rows
179: - iscol - Local index set for the columns
181: Output Arguments:
182: . newmat - New serial Mat
184: Level: developer
186: Notes:
187: Most will use MatGetLocalSubMatrix() which returns a real matrix corresponding to the local
188: block if it available, such as with matrix formats that store these blocks separately.
190: The new matrix forwards MatSetValuesLocal() and MatSetValuesBlockedLocal() to the global system.
191: In general, it does not define MatMult() or any other functions. Local submatrices can be nested.
193: .seealso: MatSetValuesLocal(), MatSetValuesBlockedLocal(), MatGetLocalSubMatrix(), MatCreateSubMatrix()
194: @*/
195: PetscErrorCode MatCreateLocalRef(Mat A,IS isrow,IS iscol,Mat *newmat)
196: {
198: Mat_LocalRef *lr;
199: Mat B;
200: PetscInt m,n;
201: PetscBool islr;
208: *newmat = 0;
210: MatCreate(PETSC_COMM_SELF,&B);
211: ISGetLocalSize(isrow,&m);
212: ISGetLocalSize(iscol,&n);
213: MatSetSizes(B,m,n,m,n);
214: PetscObjectChangeTypeName((PetscObject)B,MATLOCALREF);
215: MatSetUp(B);
217: B->ops->destroy = MatDestroy_LocalRef;
219: PetscNewLog(B,Mat_LocalRef,&lr);
220: B->data = (void*)lr;
222: PetscObjectTypeCompare((PetscObject)A,MATLOCALREF,&islr);
223: if (islr) {
224: Mat_LocalRef *alr = (Mat_LocalRef*)A->data;
225: lr->Top = alr->Top;
226: } else {
227: /* This does not increase the reference count because MatLocalRef is not allowed to live longer than its parent */
228: lr->Top = A;
229: }
230: {
231: ISLocalToGlobalMapping rltog,cltog;
232: PetscInt abs,rbs,cbs;
234: /* We will translate directly to global indices for the top level */
235: lr->SetValues = MatSetValues;
236: lr->SetValuesBlocked = MatSetValuesBlocked;
238: B->ops->setvalueslocal = MatSetValuesLocal_LocalRef_Scalar;
239: ISL2GCompose(isrow,A->rmap->mapping,&rltog);
240: if (isrow == iscol && A->rmap->mapping == A->cmap->mapping) {
241: PetscObjectReference((PetscObject)rltog);
242: cltog = rltog;
243: } else {
244: ISL2GCompose(iscol,A->cmap->mapping,&cltog);
245: }
246: MatSetLocalToGlobalMapping(B,rltog,cltog);
247: ISLocalToGlobalMappingDestroy(&rltog);
248: ISLocalToGlobalMappingDestroy(&cltog);
250: MatGetBlockSize(A,&abs);
251: ISGetBlockSize(isrow,&rbs);
252: ISGetBlockSize(iscol,&cbs);
253: if (rbs == cbs) { /* submatrix has block structure, so user can insert values with blocked interface */
254: PetscLayoutSetBlockSize(B->rmap,rbs);
255: PetscLayoutSetBlockSize(B->cmap,cbs);
256: if (abs != rbs || abs == 1) {
257: /* Top-level matrix has different block size, so we have to call its scalar insertion interface */
258: B->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_LocalRef_Scalar;
259: } else {
260: /* Block sizes match so we can forward values to the top level using the block interface */
261: B->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_LocalRef_Block;
262: ISL2GComposeBlock(isrow,A->rmap->bmapping,&rltog);
263: if (isrow == iscol && A->rmap->bmapping == A->cmap->bmapping) {
264: PetscObjectReference((PetscObject)rltog);
265: cltog = rltog;
266: } else {
267: ISL2GComposeBlock(iscol,A->cmap->bmapping,&cltog);
268: }
269: MatSetLocalToGlobalMappingBlock(B,rltog,cltog);
270: ISLocalToGlobalMappingDestroy(&rltog);
271: ISLocalToGlobalMappingDestroy(&cltog);
272: }
273: }
274: }
275: *newmat = B;
276: return(0);
277: }