xref: /petsc/src/mat/impls/localref/mlocalref.c (revision fe998a80077c9ee0917a39496df43fc256e1b478) 
1 
2 #include <petsc-private/matimpl.h>          /*I "petscmat.h" I*/
3 
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;
9 
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       ierr = PetscMalloc2(nrow,&irowm,ncol,&icolm);CHKERRQ(ierr); \
14     } else {                                                            \
15       irowm = &buf[0];                                                  \
16       icolm = &buf[nrow];                                               \
17     }                                                                   \
18 } while (0)
19 
20 #define IndexSpaceRestore(buf,nrow,ncol,irowm,icolm) do {       \
21     if (nrow + ncol > (PetscInt)(sizeof(buf)/sizeof(buf[0]))) { \
22       ierr = PetscFree2(irowm,icolm);CHKERRQ(ierr);             \
23     }                                                           \
24 } while (0)
25 
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 }
35 
36 #undef __FUNCT__
37 #define __FUNCT__ "MatSetValuesBlockedLocal_LocalRef_Block"
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   PetscErrorCode ierr;
42   PetscInt       buf[4096],*irowm,*icolm;
43 
44   PetscFunctionBegin;
45   if (!nrow || !ncol) PetscFunctionReturn(0);
46   IndexSpaceGet(buf,nrow,ncol,irowm,icolm);
47   ierr = ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,nrow,irow,irowm);CHKERRQ(ierr);
48   ierr = ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,ncol,icol,icolm);CHKERRQ(ierr);
49   ierr = (*lr->SetValuesBlocked)(lr->Top,nrow,irowm,ncol,icolm,y,addv);CHKERRQ(ierr);
50   IndexSpaceRestore(buf,nrow,ncol,irowm,icolm);
51   PetscFunctionReturn(0);
52 }
53 
54 #undef __FUNCT__
55 #define __FUNCT__ "MatSetValuesBlockedLocal_LocalRef_Scalar"
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;
59   PetscErrorCode ierr;
60   PetscInt       bs,buf[4096],*irowm,*icolm;
61 
62   PetscFunctionBegin;
63   ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr);
64   IndexSpaceGet(buf,nrow*bs,ncol*bs,irowm,icolm);
65   BlockIndicesExpand(nrow,irow,bs,irowm);
66   BlockIndicesExpand(ncol,icol,bs,icolm);
67   ierr = ISLocalToGlobalMappingApply(A->rmap->mapping,nrow*bs,irowm,irowm);CHKERRQ(ierr);
68   ierr = ISLocalToGlobalMappingApply(A->cmap->mapping,ncol*bs,icolm,icolm);CHKERRQ(ierr);
69   ierr = (*lr->SetValues)(lr->Top,nrow*bs,irowm,ncol*bs,icolm,y,addv);CHKERRQ(ierr);
70   IndexSpaceRestore(buf,nrow*bs,ncol*bs,irowm,icolm);
71   PetscFunctionReturn(0);
72 }
73 
74 #undef __FUNCT__
75 #define __FUNCT__ "MatSetValuesLocal_LocalRef_Scalar"
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;
79   PetscErrorCode ierr;
80   PetscInt       buf[4096],*irowm,*icolm;
81 
82   PetscFunctionBegin;
83   IndexSpaceGet(buf,nrow,ncol,irowm,icolm);
84   ierr = ISLocalToGlobalMappingApply(A->rmap->mapping,nrow,irow,irowm);CHKERRQ(ierr);
85   ierr = ISLocalToGlobalMappingApply(A->cmap->mapping,ncol,icol,icolm);CHKERRQ(ierr);
86   ierr = (*lr->SetValues)(lr->Top,nrow,irowm,ncol,icolm,y,addv);CHKERRQ(ierr);
87   IndexSpaceRestore(buf,nrow,ncol,irowm,icolm);
88   PetscFunctionReturn(0);
89 }
90 
91 #undef __FUNCT__
92 #define __FUNCT__ "ISL2GCompose"
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 {
96   PetscErrorCode ierr;
97   const PetscInt *idx;
98   PetscInt       m,*idxm;
99   PetscBool      isblock;
100 
101   PetscFunctionBegin;
102   PetscValidHeaderSpecific(is,IS_CLASSID,1);
103   PetscValidHeaderSpecific(ltog,IS_LTOGM_CLASSID,2);
104   PetscValidPointer(cltog,3);
105   ierr = PetscObjectTypeCompare((PetscObject)is,ISBLOCK,&isblock);CHKERRQ(ierr);
106   if (isblock) {
107     PetscInt bs,lbs;
108 
109     ierr = ISGetBlockSize(is,&bs);CHKERRQ(ierr);
110     ierr = ISLocalToGlobalMappingGetBlockSize(ltog,&lbs);CHKERRQ(ierr);
111     if (bs == lbs) {
112       ierr = ISGetLocalSize(is,&m);CHKERRQ(ierr);
113       m    = m/bs;
114       ierr = ISBlockGetIndices(is,&idx);CHKERRQ(ierr);
115       ierr = PetscMalloc1(m,&idxm);CHKERRQ(ierr);
116       ierr = ISLocalToGlobalMappingApplyBlock(ltog,m,idx,idxm);CHKERRQ(ierr);
117       ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)is),bs,m,idxm,PETSC_OWN_POINTER,cltog);CHKERRQ(ierr);
118       ierr = ISBlockRestoreIndices(is,&idx);CHKERRQ(ierr);
119       PetscFunctionReturn(0);
120     }
121   }
122   ierr = ISGetLocalSize(is,&m);CHKERRQ(ierr);
123   ierr = ISGetIndices(is,&idx);CHKERRQ(ierr);
124   ierr = PetscMalloc1(m,&idxm);CHKERRQ(ierr);
125   if (ltog) {
126     ierr = ISLocalToGlobalMappingApply(ltog,m,idx,idxm);CHKERRQ(ierr);
127   } else {
128     ierr = PetscMemcpy(idxm,idx,m*sizeof(PetscInt));CHKERRQ(ierr);
129   }
130   ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)is),1,m,idxm,PETSC_OWN_POINTER,cltog);CHKERRQ(ierr);
131   ierr = ISRestoreIndices(is,&idx);CHKERRQ(ierr);
132   PetscFunctionReturn(0);
133 }
134 
135 #undef __FUNCT__
136 #define __FUNCT__ "ISL2GComposeBlock"
137 static PetscErrorCode ISL2GComposeBlock(IS is,ISLocalToGlobalMapping ltog,ISLocalToGlobalMapping *cltog)
138 {
139   PetscErrorCode ierr;
140   const PetscInt *idx;
141   PetscInt       m,*idxm,bs;
142 
143   PetscFunctionBegin;
144   PetscValidHeaderSpecific(is,IS_CLASSID,1);
145   PetscValidHeaderSpecific(ltog,IS_LTOGM_CLASSID,2);
146   PetscValidPointer(cltog,3);
147   ierr = ISBlockGetLocalSize(is,&m);CHKERRQ(ierr);
148   ierr = ISBlockGetIndices(is,&idx);CHKERRQ(ierr);
149   ierr = ISLocalToGlobalMappingGetBlockSize(ltog,&bs);CHKERRQ(ierr);
150   ierr = PetscMalloc1(m,&idxm);CHKERRQ(ierr);
151   if (ltog) {
152     ierr = ISLocalToGlobalMappingApplyBlock(ltog,m,idx,idxm);CHKERRQ(ierr);
153   } else {
154     ierr = PetscMemcpy(idxm,idx,m*sizeof(PetscInt));CHKERRQ(ierr);
155   }
156   ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)is),bs,m,idxm,PETSC_OWN_POINTER,cltog);CHKERRQ(ierr);
157   ierr = ISBlockRestoreIndices(is,&idx);CHKERRQ(ierr);
158   PetscFunctionReturn(0);
159 }
160 
161 #undef __FUNCT__
162 #define __FUNCT__ "MatDestroy_LocalRef"
163 static PetscErrorCode MatDestroy_LocalRef(Mat B)
164 {
165   PetscErrorCode ierr;
166 
167   PetscFunctionBegin;
168   ierr = PetscFree(B->data);CHKERRQ(ierr);
169   PetscFunctionReturn(0);
170 }
171 
172 
173 #undef __FUNCT__
174 #define __FUNCT__ "MatCreateLocalRef"
175 /*@
176    MatCreateLocalRef - Gets a logical reference to a local submatrix, for use in assembly
177 
178    Not Collective
179 
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
184 
185    Output Arguments:
186 . newmat - New serial Mat
187 
188    Level: developer
189 
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.
193 
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.
196 
197 .seealso: MatSetValuesLocal(), MatSetValuesBlockedLocal(), MatGetLocalSubMatrix(), MatCreateSubMatrix()
198 @*/
199 PetscErrorCode  MatCreateLocalRef(Mat A,IS isrow,IS iscol,Mat *newmat)
200 {
201   PetscErrorCode ierr;
202   Mat_LocalRef   *lr;
203   Mat            B;
204   PetscInt       m,n;
205   PetscBool      islr;
206 
207   PetscFunctionBegin;
208   PetscValidHeaderSpecific(A,MAT_CLASSID,1);
209   PetscValidHeaderSpecific(isrow,IS_CLASSID,2);
210   PetscValidHeaderSpecific(iscol,IS_CLASSID,3);
211   PetscValidPointer(newmat,4);
212   *newmat = 0;
213 
214   ierr = MatCreate(PETSC_COMM_SELF,&B);CHKERRQ(ierr);
215   ierr = ISGetLocalSize(isrow,&m);CHKERRQ(ierr);
216   ierr = ISGetLocalSize(iscol,&n);CHKERRQ(ierr);
217   ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr);
218   ierr = PetscObjectChangeTypeName((PetscObject)B,MATLOCALREF);CHKERRQ(ierr);
219   ierr = MatSetUp(B);CHKERRQ(ierr);
220 
221   B->ops->destroy = MatDestroy_LocalRef;
222 
223   ierr    = PetscNewLog(B,&lr);CHKERRQ(ierr);
224   B->data = (void*)lr;
225 
226   ierr = PetscObjectTypeCompare((PetscObject)A,MATLOCALREF,&islr);CHKERRQ(ierr);
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;
237 
238     /* We will translate directly to global indices for the top level */
239     lr->SetValues        = MatSetValues;
240     lr->SetValuesBlocked = MatSetValuesBlocked;
241 
242     B->ops->setvalueslocal = MatSetValuesLocal_LocalRef_Scalar;
243 
244     ierr = ISL2GCompose(isrow,A->rmap->mapping,&rltog);CHKERRQ(ierr);
245     if (isrow == iscol && A->rmap->mapping == A->cmap->mapping) {
246       ierr  = PetscObjectReference((PetscObject)rltog);CHKERRQ(ierr);
247       cltog = rltog;
248     } else {
249       ierr = ISL2GCompose(iscol,A->cmap->mapping,&cltog);CHKERRQ(ierr);
250     }
251     ierr = MatSetLocalToGlobalMapping(B,rltog,cltog);CHKERRQ(ierr);
252     ierr = ISLocalToGlobalMappingDestroy(&rltog);CHKERRQ(ierr);
253     ierr = ISLocalToGlobalMappingDestroy(&cltog);CHKERRQ(ierr);
254 
255     ierr = MatGetBlockSize(A,&abs);CHKERRQ(ierr);
256     ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr);
257     ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr);
258     if (rbs == cbs) {           /* submatrix has block structure, so user can insert values with blocked interface */
259       ierr = PetscLayoutSetBlockSize(B->rmap,rbs);CHKERRQ(ierr);
260       ierr = PetscLayoutSetBlockSize(B->cmap,cbs);CHKERRQ(ierr);
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;
267 
268         ierr = ISL2GComposeBlock(isrow,A->rmap->mapping,&rltog);CHKERRQ(ierr);
269         if (isrow == iscol && A->rmap->mapping == A->cmap->mapping) {
270           ierr  =  PetscObjectReference((PetscObject)rltog);CHKERRQ(ierr);
271           cltog = rltog;
272         } else {
273           ierr = ISL2GComposeBlock(iscol,A->cmap->mapping,&cltog);CHKERRQ(ierr);
274         }
275         ierr = MatSetLocalToGlobalMapping(B,rltog,cltog);CHKERRQ(ierr);
276         ierr = ISLocalToGlobalMappingDestroy(&rltog);CHKERRQ(ierr);
277         ierr = ISLocalToGlobalMappingDestroy(&cltog);CHKERRQ(ierr);
278       }
279     }
280   }
281   *newmat = B;
282   PetscFunctionReturn(0);
283 }
284