1 /*
2 Creates a matrix class for using the Neumann-Neumann type preconditioners.
3 This stores the matrices in globally unassembled form. Each processor
4 assembles only its local Neumann problem and the parallel matrix vector
5 product is handled "implicitly".
6
7 Currently this allows for only one subdomain per processor.
8 */
9
10 #include <petsc/private/matisimpl.h> /*I "petscmat.h" I*/
11 #include <../src/mat/impls/aij/mpi/mpiaij.h>
12 #include <petsc/private/sfimpl.h>
13 #include <petsc/private/vecimpl.h>
14 #include <petsc/private/hashseti.h>
15
16 #define MATIS_MAX_ENTRIES_INSERTION 2048
17
18 /* copied from src/mat/impls/localref/mlocalref.c */
19 #define IndexSpaceGet(buf, nrow, ncol, irowm, icolm) \
20 do { \
21 if (nrow + ncol > (PetscInt)PETSC_STATIC_ARRAY_LENGTH(buf)) { \
22 PetscCall(PetscMalloc2(nrow, &irowm, ncol, &icolm)); \
23 } else { \
24 irowm = &buf[0]; \
25 icolm = &buf[nrow]; \
26 } \
27 } while (0)
28
29 #define IndexSpaceRestore(buf, nrow, ncol, irowm, icolm) \
30 do { \
31 if (nrow + ncol > (PetscInt)PETSC_STATIC_ARRAY_LENGTH(buf)) PetscCall(PetscFree2(irowm, icolm)); \
32 } while (0)
33
BlockIndicesExpand(PetscInt n,const PetscInt idx[],PetscInt bs,PetscInt idxm[])34 static void BlockIndicesExpand(PetscInt n, const PetscInt idx[], PetscInt bs, PetscInt idxm[])
35 {
36 for (PetscInt i = 0; i < n; i++) {
37 for (PetscInt j = 0; j < bs; j++) idxm[i * bs + j] = idx[i] * bs + j;
38 }
39 }
40
41 static PetscErrorCode MatSetValuesLocal_IS(Mat, PetscInt, const PetscInt *, PetscInt, const PetscInt *, const PetscScalar *, InsertMode);
42 static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat, PetscInt, const PetscInt *, PetscInt, const PetscInt *, const PetscScalar *, InsertMode);
43 static PetscErrorCode MatISSetUpScatters_Private(Mat);
44
MatISContainerDestroyPtAP_Private(PetscCtxRt ptr)45 static PetscErrorCode MatISContainerDestroyPtAP_Private(PetscCtxRt ptr)
46 {
47 MatISPtAP ptap = *(MatISPtAP *)ptr;
48
49 PetscFunctionBegin;
50 PetscCall(MatDestroySubMatrices(ptap->ris1 ? 2 : 1, &ptap->lP));
51 PetscCall(ISDestroy(&ptap->cis0));
52 PetscCall(ISDestroy(&ptap->cis1));
53 PetscCall(ISDestroy(&ptap->ris0));
54 PetscCall(ISDestroy(&ptap->ris1));
55 PetscCall(PetscFree(ptap));
56 PetscFunctionReturn(PETSC_SUCCESS);
57 }
58
MatPtAPNumeric_IS_XAIJ(Mat A,Mat P,Mat C)59 static PetscErrorCode MatPtAPNumeric_IS_XAIJ(Mat A, Mat P, Mat C)
60 {
61 MatISPtAP ptap;
62 Mat_IS *matis = (Mat_IS *)A->data;
63 Mat lA, lC;
64 MatReuse reuse;
65 IS ris[2], cis[2];
66 PetscContainer c;
67 PetscInt n;
68
69 PetscFunctionBegin;
70 PetscCall(PetscObjectQuery((PetscObject)C, "_MatIS_PtAP", (PetscObject *)&c));
71 PetscCheck(c, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Missing PtAP information");
72 PetscCall(PetscContainerGetPointer(c, &ptap));
73 ris[0] = ptap->ris0;
74 ris[1] = ptap->ris1;
75 cis[0] = ptap->cis0;
76 cis[1] = ptap->cis1;
77 n = ptap->ris1 ? 2 : 1;
78 reuse = ptap->lP ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX;
79 PetscCall(MatCreateSubMatrices(P, n, ris, cis, reuse, &ptap->lP));
80
81 PetscCall(MatISGetLocalMat(A, &lA));
82 PetscCall(MatISGetLocalMat(C, &lC));
83 if (ptap->ris1) { /* unsymmetric A mapping */
84 Mat lPt;
85
86 PetscCall(MatTranspose(ptap->lP[1], MAT_INITIAL_MATRIX, &lPt));
87 PetscCall(MatMatMatMult(lPt, lA, ptap->lP[0], reuse, ptap->fill, &lC));
88 if (matis->storel2l) PetscCall(PetscObjectCompose((PetscObject)A, "_MatIS_PtAP_l2l", (PetscObject)lPt));
89 PetscCall(MatDestroy(&lPt));
90 } else {
91 PetscCall(MatPtAP(lA, ptap->lP[0], reuse, ptap->fill, &lC));
92 if (matis->storel2l) PetscCall(PetscObjectCompose((PetscObject)C, "_MatIS_PtAP_l2l", (PetscObject)ptap->lP[0]));
93 }
94 if (reuse == MAT_INITIAL_MATRIX) {
95 PetscCall(MatISSetLocalMat(C, lC));
96 PetscCall(MatDestroy(&lC));
97 }
98 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
99 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
100 PetscFunctionReturn(PETSC_SUCCESS);
101 }
102
MatGetNonzeroColumnsLocal_Private(Mat PT,IS * cis)103 static PetscErrorCode MatGetNonzeroColumnsLocal_Private(Mat PT, IS *cis)
104 {
105 Mat Po, Pd;
106 IS zd, zo;
107 const PetscInt *garray;
108 PetscInt *aux, i, bs;
109 PetscInt dc, stc, oc, ctd, cto;
110 PetscBool ismpiaij, ismpibaij, isseqaij, isseqbaij;
111 MPI_Comm comm;
112
113 PetscFunctionBegin;
114 PetscValidHeaderSpecific(PT, MAT_CLASSID, 1);
115 PetscAssertPointer(cis, 2);
116 PetscCall(PetscObjectGetComm((PetscObject)PT, &comm));
117 bs = 1;
118 PetscCall(PetscObjectBaseTypeCompare((PetscObject)PT, MATMPIAIJ, &ismpiaij));
119 PetscCall(PetscObjectBaseTypeCompare((PetscObject)PT, MATMPIBAIJ, &ismpibaij));
120 PetscCall(PetscObjectBaseTypeCompare((PetscObject)PT, MATSEQAIJ, &isseqaij));
121 PetscCall(PetscObjectTypeCompare((PetscObject)PT, MATSEQBAIJ, &isseqbaij));
122 if (isseqaij || isseqbaij) {
123 Pd = PT;
124 Po = NULL;
125 garray = NULL;
126 } else if (ismpiaij) {
127 PetscCall(MatMPIAIJGetSeqAIJ(PT, &Pd, &Po, &garray));
128 } else if (ismpibaij) {
129 PetscCall(MatMPIBAIJGetSeqBAIJ(PT, &Pd, &Po, &garray));
130 PetscCall(MatGetBlockSize(PT, &bs));
131 } else SETERRQ(comm, PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)PT)->type_name);
132
133 /* identify any null columns in Pd or Po */
134 /* We use a tolerance comparison since it may happen that, with geometric multigrid,
135 some of the columns are not really zero, but very close to */
136 zo = zd = NULL;
137 if (Po) PetscCall(MatFindNonzeroRowsOrCols_Basic(Po, PETSC_TRUE, PETSC_SMALL, &zo));
138 PetscCall(MatFindNonzeroRowsOrCols_Basic(Pd, PETSC_TRUE, PETSC_SMALL, &zd));
139
140 PetscCall(MatGetLocalSize(PT, NULL, &dc));
141 PetscCall(MatGetOwnershipRangeColumn(PT, &stc, NULL));
142 if (Po) PetscCall(MatGetLocalSize(Po, NULL, &oc));
143 else oc = 0;
144 PetscCall(PetscMalloc1((dc + oc) / bs, &aux));
145 if (zd) {
146 const PetscInt *idxs;
147 PetscInt nz;
148
149 /* this will throw an error if bs is not valid */
150 PetscCall(ISSetBlockSize(zd, bs));
151 PetscCall(ISGetLocalSize(zd, &nz));
152 PetscCall(ISGetIndices(zd, &idxs));
153 ctd = nz / bs;
154 for (i = 0; i < ctd; i++) aux[i] = (idxs[bs * i] + stc) / bs;
155 PetscCall(ISRestoreIndices(zd, &idxs));
156 } else {
157 ctd = dc / bs;
158 for (i = 0; i < ctd; i++) aux[i] = i + stc / bs;
159 }
160 if (zo) {
161 const PetscInt *idxs;
162 PetscInt nz;
163
164 /* this will throw an error if bs is not valid */
165 PetscCall(ISSetBlockSize(zo, bs));
166 PetscCall(ISGetLocalSize(zo, &nz));
167 PetscCall(ISGetIndices(zo, &idxs));
168 cto = nz / bs;
169 for (i = 0; i < cto; i++) aux[i + ctd] = garray[idxs[bs * i] / bs];
170 PetscCall(ISRestoreIndices(zo, &idxs));
171 } else {
172 cto = oc / bs;
173 for (i = 0; i < cto; i++) aux[i + ctd] = garray[i];
174 }
175 PetscCall(ISCreateBlock(comm, bs, ctd + cto, aux, PETSC_OWN_POINTER, cis));
176 PetscCall(ISDestroy(&zd));
177 PetscCall(ISDestroy(&zo));
178 PetscFunctionReturn(PETSC_SUCCESS);
179 }
180
MatPtAPSymbolic_IS_XAIJ(Mat A,Mat P,PetscReal fill,Mat C)181 static PetscErrorCode MatPtAPSymbolic_IS_XAIJ(Mat A, Mat P, PetscReal fill, Mat C)
182 {
183 Mat PT, lA;
184 MatISPtAP ptap;
185 ISLocalToGlobalMapping Crl2g, Ccl2g, rl2g, cl2g;
186 PetscContainer c;
187 MatType lmtype;
188 const PetscInt *garray;
189 PetscInt ibs, N, dc;
190 MPI_Comm comm;
191
192 PetscFunctionBegin;
193 PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
194 PetscCall(MatSetType(C, MATIS));
195 PetscCall(MatISGetLocalMat(A, &lA));
196 PetscCall(MatGetType(lA, &lmtype));
197 PetscCall(MatISSetLocalMatType(C, lmtype));
198 PetscCall(MatGetSize(P, NULL, &N));
199 PetscCall(MatGetLocalSize(P, NULL, &dc));
200 PetscCall(MatSetSizes(C, dc, dc, N, N));
201 /* Not sure about this
202 PetscCall(MatGetBlockSizes(P,NULL,&ibs));
203 PetscCall(MatSetBlockSize(*C,ibs));
204 */
205
206 PetscCall(PetscNew(&ptap));
207 PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &c));
208 PetscCall(PetscContainerSetPointer(c, ptap));
209 PetscCall(PetscContainerSetCtxDestroy(c, MatISContainerDestroyPtAP_Private));
210 PetscCall(PetscObjectCompose((PetscObject)C, "_MatIS_PtAP", (PetscObject)c));
211 PetscCall(PetscContainerDestroy(&c));
212 ptap->fill = fill;
213
214 PetscCall(MatISGetLocalToGlobalMapping(A, &rl2g, &cl2g));
215
216 PetscCall(ISLocalToGlobalMappingGetBlockSize(cl2g, &ibs));
217 PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &N));
218 PetscCall(ISLocalToGlobalMappingGetBlockIndices(cl2g, &garray));
219 PetscCall(ISCreateBlock(comm, ibs, N / ibs, garray, PETSC_COPY_VALUES, &ptap->ris0));
220 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(cl2g, &garray));
221
222 PetscCall(MatCreateSubMatrix(P, ptap->ris0, NULL, MAT_INITIAL_MATRIX, &PT));
223 PetscCall(MatGetNonzeroColumnsLocal_Private(PT, &ptap->cis0));
224 PetscCall(ISLocalToGlobalMappingCreateIS(ptap->cis0, &Ccl2g));
225 PetscCall(MatDestroy(&PT));
226
227 Crl2g = NULL;
228 if (rl2g != cl2g) { /* unsymmetric A mapping */
229 PetscBool same, lsame = PETSC_FALSE;
230 PetscInt N1, ibs1;
231
232 PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &N1));
233 PetscCall(ISLocalToGlobalMappingGetBlockSize(rl2g, &ibs1));
234 PetscCall(ISLocalToGlobalMappingGetBlockIndices(rl2g, &garray));
235 PetscCall(ISCreateBlock(comm, ibs, N1 / ibs, garray, PETSC_COPY_VALUES, &ptap->ris1));
236 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(rl2g, &garray));
237 if (ibs1 == ibs && N1 == N) { /* check if the l2gmaps are the same */
238 const PetscInt *i1, *i2;
239
240 PetscCall(ISBlockGetIndices(ptap->ris0, &i1));
241 PetscCall(ISBlockGetIndices(ptap->ris1, &i2));
242 PetscCall(PetscArraycmp(i1, i2, N, &lsame));
243 }
244 PetscCallMPI(MPIU_Allreduce(&lsame, &same, 1, MPI_C_BOOL, MPI_LAND, comm));
245 if (same) {
246 PetscCall(ISDestroy(&ptap->ris1));
247 } else {
248 PetscCall(MatCreateSubMatrix(P, ptap->ris1, NULL, MAT_INITIAL_MATRIX, &PT));
249 PetscCall(MatGetNonzeroColumnsLocal_Private(PT, &ptap->cis1));
250 PetscCall(ISLocalToGlobalMappingCreateIS(ptap->cis1, &Crl2g));
251 PetscCall(MatDestroy(&PT));
252 }
253 }
254 /* Not sure about this
255 if (!Crl2g) {
256 PetscCall(MatGetBlockSize(C,&ibs));
257 PetscCall(ISLocalToGlobalMappingSetBlockSize(Ccl2g,ibs));
258 }
259 */
260 PetscCall(MatSetLocalToGlobalMapping(C, Crl2g ? Crl2g : Ccl2g, Ccl2g));
261 PetscCall(ISLocalToGlobalMappingDestroy(&Crl2g));
262 PetscCall(ISLocalToGlobalMappingDestroy(&Ccl2g));
263
264 C->ops->ptapnumeric = MatPtAPNumeric_IS_XAIJ;
265 PetscFunctionReturn(PETSC_SUCCESS);
266 }
267
MatProductSymbolic_PtAP_IS_XAIJ(Mat C)268 static PetscErrorCode MatProductSymbolic_PtAP_IS_XAIJ(Mat C)
269 {
270 Mat_Product *product = C->product;
271 Mat A = product->A, P = product->B;
272 PetscReal fill = product->fill;
273
274 PetscFunctionBegin;
275 PetscCall(MatPtAPSymbolic_IS_XAIJ(A, P, fill, C));
276 C->ops->productnumeric = MatProductNumeric_PtAP;
277 PetscFunctionReturn(PETSC_SUCCESS);
278 }
279
MatProductSetFromOptions_IS_XAIJ_PtAP(Mat C)280 static PetscErrorCode MatProductSetFromOptions_IS_XAIJ_PtAP(Mat C)
281 {
282 PetscFunctionBegin;
283 C->ops->productsymbolic = MatProductSymbolic_PtAP_IS_XAIJ;
284 PetscFunctionReturn(PETSC_SUCCESS);
285 }
286
MatProductSetFromOptions_IS_XAIJ(Mat C)287 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat C)
288 {
289 Mat_Product *product = C->product;
290
291 PetscFunctionBegin;
292 if (product->type == MATPRODUCT_PtAP) PetscCall(MatProductSetFromOptions_IS_XAIJ_PtAP(C));
293 PetscFunctionReturn(PETSC_SUCCESS);
294 }
295
MatISContainerDestroyFields_Private(PetscCtxRt ptr)296 static PetscErrorCode MatISContainerDestroyFields_Private(PetscCtxRt ptr)
297 {
298 MatISLocalFields lf = *(MatISLocalFields *)ptr;
299 PetscInt i;
300
301 PetscFunctionBegin;
302 for (i = 0; i < lf->nr; i++) PetscCall(ISDestroy(&lf->rf[i]));
303 for (i = 0; i < lf->nc; i++) PetscCall(ISDestroy(&lf->cf[i]));
304 PetscCall(PetscFree2(lf->rf, lf->cf));
305 PetscCall(PetscFree(lf));
306 PetscFunctionReturn(PETSC_SUCCESS);
307 }
308
MatConvert_SeqXAIJ_IS(Mat A,MatType type,MatReuse reuse,Mat * newmat)309 static PetscErrorCode MatConvert_SeqXAIJ_IS(Mat A, MatType type, MatReuse reuse, Mat *newmat)
310 {
311 Mat B, lB;
312
313 PetscFunctionBegin;
314 if (reuse != MAT_REUSE_MATRIX) {
315 ISLocalToGlobalMapping rl2g, cl2g;
316 PetscInt bs;
317 IS is;
318
319 PetscCall(MatGetBlockSize(A, &bs));
320 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->rmap->n / bs, 0, 1, &is));
321 if (bs > 1) {
322 IS is2;
323 PetscInt i, *aux;
324
325 PetscCall(ISGetLocalSize(is, &i));
326 PetscCall(ISGetIndices(is, (const PetscInt **)&aux));
327 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), bs, i, aux, PETSC_COPY_VALUES, &is2));
328 PetscCall(ISRestoreIndices(is, (const PetscInt **)&aux));
329 PetscCall(ISDestroy(&is));
330 is = is2;
331 }
332 PetscCall(ISSetIdentity(is));
333 PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
334 PetscCall(ISDestroy(&is));
335 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->cmap->n / bs, 0, 1, &is));
336 if (bs > 1) {
337 IS is2;
338 PetscInt i, *aux;
339
340 PetscCall(ISGetLocalSize(is, &i));
341 PetscCall(ISGetIndices(is, (const PetscInt **)&aux));
342 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), bs, i, aux, PETSC_COPY_VALUES, &is2));
343 PetscCall(ISRestoreIndices(is, (const PetscInt **)&aux));
344 PetscCall(ISDestroy(&is));
345 is = is2;
346 }
347 PetscCall(ISSetIdentity(is));
348 PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
349 PetscCall(ISDestroy(&is));
350 PetscCall(MatCreateIS(PetscObjectComm((PetscObject)A), bs, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N, rl2g, cl2g, &B));
351 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
352 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
353 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &lB));
354 if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
355 } else {
356 B = *newmat;
357 PetscCall(PetscObjectReference((PetscObject)A));
358 lB = A;
359 }
360 PetscCall(MatISSetLocalMat(B, lB));
361 PetscCall(MatDestroy(&lB));
362 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
363 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
364 if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B));
365 PetscFunctionReturn(PETSC_SUCCESS);
366 }
367
MatISScaleDisassembling_Private(Mat A)368 static PetscErrorCode MatISScaleDisassembling_Private(Mat A)
369 {
370 Mat_IS *matis = (Mat_IS *)A->data;
371 PetscScalar *aa;
372 const PetscInt *ii, *jj;
373 PetscInt i, n, m;
374 PetscInt *ecount, **eneighs;
375 PetscBool flg;
376
377 PetscFunctionBegin;
378 PetscCall(MatGetRowIJ(matis->A, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg));
379 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
380 PetscCall(ISLocalToGlobalMappingGetNodeInfo(matis->rmapping, &n, &ecount, &eneighs));
381 PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected %" PetscInt_FMT " != %" PetscInt_FMT, m, n);
382 PetscCall(MatSeqAIJGetArray(matis->A, &aa));
383 for (i = 0; i < n; i++) {
384 if (ecount[i] > 1) {
385 PetscInt j;
386
387 for (j = ii[i]; j < ii[i + 1]; j++) {
388 PetscInt i2 = jj[j], p, p2;
389 PetscReal scal = 0.0;
390
391 for (p = 0; p < ecount[i]; p++) {
392 for (p2 = 0; p2 < ecount[i2]; p2++) {
393 if (eneighs[i][p] == eneighs[i2][p2]) {
394 scal += 1.0;
395 break;
396 }
397 }
398 }
399 if (scal) aa[j] /= scal;
400 }
401 }
402 }
403 PetscCall(ISLocalToGlobalMappingRestoreNodeInfo(matis->rmapping, &n, &ecount, &eneighs));
404 PetscCall(MatSeqAIJRestoreArray(matis->A, &aa));
405 PetscCall(MatRestoreRowIJ(matis->A, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg));
406 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot restore IJ structure");
407 PetscFunctionReturn(PETSC_SUCCESS);
408 }
409
410 typedef enum {
411 MAT_IS_DISASSEMBLE_L2G_NATURAL,
412 MAT_IS_DISASSEMBLE_L2G_MAT,
413 MAT_IS_DISASSEMBLE_L2G_ND
414 } MatISDisassemblel2gType;
415
MatMPIXAIJComputeLocalToGlobalMapping_Private(Mat A,ISLocalToGlobalMapping * l2g)416 static PetscErrorCode MatMPIXAIJComputeLocalToGlobalMapping_Private(Mat A, ISLocalToGlobalMapping *l2g)
417 {
418 Mat Ad, Ao;
419 IS is, ndmap, ndsub;
420 MPI_Comm comm;
421 const PetscInt *garray, *ndmapi;
422 PetscInt bs, i, cnt, nl, *ncount, *ndmapc;
423 PetscBool ismpiaij, ismpibaij;
424 const char *const MatISDisassemblel2gTypes[] = {"NATURAL", "MAT", "ND", "MatISDisassemblel2gType", "MAT_IS_DISASSEMBLE_L2G_", NULL};
425 MatISDisassemblel2gType mode = MAT_IS_DISASSEMBLE_L2G_NATURAL;
426 MatPartitioning part;
427 PetscSF sf;
428 PetscObject dm;
429
430 PetscFunctionBegin;
431 PetscOptionsBegin(PetscObjectComm((PetscObject)A), ((PetscObject)A)->prefix, "MatIS l2g disassembling options", "Mat");
432 PetscCall(PetscOptionsEnum("-mat_is_disassemble_l2g_type", "Type of local-to-global mapping to be used for disassembling", "MatISDisassemblel2gType", MatISDisassemblel2gTypes, (PetscEnum)mode, (PetscEnum *)&mode, NULL));
433 PetscOptionsEnd();
434 if (mode == MAT_IS_DISASSEMBLE_L2G_MAT) {
435 PetscCall(MatGetLocalToGlobalMapping(A, l2g, NULL));
436 PetscFunctionReturn(PETSC_SUCCESS);
437 }
438 PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
439 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
440 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIBAIJ, &ismpibaij));
441 PetscCall(MatGetBlockSize(A, &bs));
442 switch (mode) {
443 case MAT_IS_DISASSEMBLE_L2G_ND:
444 PetscCall(MatPartitioningCreate(comm, &part));
445 PetscCall(MatPartitioningSetAdjacency(part, A));
446 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)part, ((PetscObject)A)->prefix));
447 PetscCall(MatPartitioningSetFromOptions(part));
448 PetscCall(MatPartitioningApplyND(part, &ndmap));
449 PetscCall(MatPartitioningDestroy(&part));
450 PetscCall(ISBuildTwoSided(ndmap, NULL, &ndsub));
451 PetscCall(MatMPIAIJSetUseScalableIncreaseOverlap(A, PETSC_TRUE));
452 PetscCall(MatIncreaseOverlap(A, 1, &ndsub, 1));
453 PetscCall(ISLocalToGlobalMappingCreateIS(ndsub, l2g));
454
455 /* it may happen that a separator node is not properly shared */
456 PetscCall(ISLocalToGlobalMappingGetNodeInfo(*l2g, &nl, &ncount, NULL));
457 PetscCall(PetscSFCreate(comm, &sf));
458 PetscCall(ISLocalToGlobalMappingGetIndices(*l2g, &garray));
459 PetscCall(PetscSFSetGraphLayout(sf, A->rmap, nl, NULL, PETSC_OWN_POINTER, garray));
460 PetscCall(ISLocalToGlobalMappingRestoreIndices(*l2g, &garray));
461 PetscCall(PetscCalloc1(A->rmap->n, &ndmapc));
462 PetscCall(PetscSFReduceBegin(sf, MPIU_INT, ncount, ndmapc, MPI_REPLACE));
463 PetscCall(PetscSFReduceEnd(sf, MPIU_INT, ncount, ndmapc, MPI_REPLACE));
464 PetscCall(ISLocalToGlobalMappingRestoreNodeInfo(*l2g, NULL, &ncount, NULL));
465 PetscCall(ISGetIndices(ndmap, &ndmapi));
466 for (i = 0, cnt = 0; i < A->rmap->n; i++)
467 if (ndmapi[i] < 0 && ndmapc[i] < 2) cnt++;
468
469 PetscCallMPI(MPIU_Allreduce(&cnt, &i, 1, MPIU_INT, MPI_MAX, comm));
470 if (i) { /* we detected isolated separator nodes */
471 Mat A2, A3;
472 IS *workis, is2;
473 PetscScalar *vals;
474 PetscInt gcnt = i, *dnz, *onz, j, *lndmapi;
475 ISLocalToGlobalMapping ll2g;
476 PetscBool flg;
477 const PetscInt *ii, *jj;
478
479 /* communicate global id of separators */
480 MatPreallocateBegin(comm, A->rmap->n, A->cmap->n, dnz, onz);
481 for (i = 0, cnt = 0; i < A->rmap->n; i++) dnz[i] = ndmapi[i] < 0 ? i + A->rmap->rstart : -1;
482
483 PetscCall(PetscMalloc1(nl, &lndmapi));
484 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, dnz, lndmapi, MPI_REPLACE));
485
486 /* compute adjacency of isolated separators node */
487 PetscCall(PetscMalloc1(gcnt, &workis));
488 for (i = 0, cnt = 0; i < A->rmap->n; i++) {
489 if (ndmapi[i] < 0 && ndmapc[i] < 2) PetscCall(ISCreateStride(comm, 1, i + A->rmap->rstart, 1, &workis[cnt++]));
490 }
491 for (i = cnt; i < gcnt; i++) PetscCall(ISCreateStride(comm, 0, 0, 1, &workis[i]));
492 for (i = 0; i < gcnt; i++) {
493 PetscCall(PetscObjectSetName((PetscObject)workis[i], "ISOLATED"));
494 PetscCall(ISViewFromOptions(workis[i], NULL, "-view_isolated_separators"));
495 }
496
497 /* no communications since all the ISes correspond to locally owned rows */
498 PetscCall(MatIncreaseOverlap(A, gcnt, workis, 1));
499
500 /* end communicate global id of separators */
501 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, dnz, lndmapi, MPI_REPLACE));
502
503 /* communicate new layers : create a matrix and transpose it */
504 PetscCall(PetscArrayzero(dnz, A->rmap->n));
505 PetscCall(PetscArrayzero(onz, A->rmap->n));
506 for (i = 0, j = 0; i < A->rmap->n; i++) {
507 if (ndmapi[i] < 0 && ndmapc[i] < 2) {
508 const PetscInt *idxs;
509 PetscInt s;
510
511 PetscCall(ISGetLocalSize(workis[j], &s));
512 PetscCall(ISGetIndices(workis[j], &idxs));
513 PetscCall(MatPreallocateSet(i + A->rmap->rstart, s, idxs, dnz, onz));
514 j++;
515 }
516 }
517 PetscCheck(j == cnt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected local count %" PetscInt_FMT " != %" PetscInt_FMT, j, cnt);
518
519 for (i = 0; i < gcnt; i++) {
520 PetscCall(PetscObjectSetName((PetscObject)workis[i], "EXTENDED"));
521 PetscCall(ISViewFromOptions(workis[i], NULL, "-view_isolated_separators"));
522 }
523
524 for (i = 0, j = 0; i < A->rmap->n; i++) j = PetscMax(j, dnz[i] + onz[i]);
525 PetscCall(PetscMalloc1(j, &vals));
526 for (i = 0; i < j; i++) vals[i] = 1.0;
527
528 PetscCall(MatCreate(comm, &A2));
529 PetscCall(MatSetType(A2, MATMPIAIJ));
530 PetscCall(MatSetSizes(A2, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
531 PetscCall(MatMPIAIJSetPreallocation(A2, 0, dnz, 0, onz));
532 PetscCall(MatSetOption(A2, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
533 for (i = 0, j = 0; i < A2->rmap->n; i++) {
534 PetscInt row = i + A2->rmap->rstart, s = dnz[i] + onz[i];
535 const PetscInt *idxs;
536
537 if (s) {
538 PetscCall(ISGetIndices(workis[j], &idxs));
539 PetscCall(MatSetValues(A2, 1, &row, s, idxs, vals, INSERT_VALUES));
540 PetscCall(ISRestoreIndices(workis[j], &idxs));
541 j++;
542 }
543 }
544 PetscCheck(j == cnt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected local count %" PetscInt_FMT " != %" PetscInt_FMT, j, cnt);
545 PetscCall(PetscFree(vals));
546 PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
547 PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
548 PetscCall(MatTranspose(A2, MAT_INPLACE_MATRIX, &A2));
549
550 /* extract submatrix corresponding to the coupling "owned separators" x "isolated separators" */
551 for (i = 0, j = 0; i < nl; i++)
552 if (lndmapi[i] >= 0) lndmapi[j++] = lndmapi[i];
553 PetscCall(ISCreateGeneral(comm, j, lndmapi, PETSC_USE_POINTER, &is));
554 PetscCall(MatMPIAIJGetLocalMatCondensed(A2, MAT_INITIAL_MATRIX, &is, NULL, &A3));
555 PetscCall(ISDestroy(&is));
556 PetscCall(MatDestroy(&A2));
557
558 /* extend local to global map to include connected isolated separators */
559 PetscCall(PetscObjectQuery((PetscObject)A3, "_petsc_GetLocalMatCondensed_iscol", (PetscObject *)&is));
560 PetscCheck(is, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing column map");
561 PetscCall(ISLocalToGlobalMappingCreateIS(is, &ll2g));
562 PetscCall(MatGetRowIJ(A3, 0, PETSC_FALSE, PETSC_FALSE, &i, &ii, &jj, &flg));
563 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
564 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ii[i], jj, PETSC_COPY_VALUES, &is));
565 PetscCall(MatRestoreRowIJ(A3, 0, PETSC_FALSE, PETSC_FALSE, &i, &ii, &jj, &flg));
566 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
567 PetscCall(ISLocalToGlobalMappingApplyIS(ll2g, is, &is2));
568 PetscCall(ISDestroy(&is));
569 PetscCall(ISLocalToGlobalMappingDestroy(&ll2g));
570
571 /* add new nodes to the local-to-global map */
572 PetscCall(ISLocalToGlobalMappingDestroy(l2g));
573 PetscCall(ISExpand(ndsub, is2, &is));
574 PetscCall(ISDestroy(&is2));
575 PetscCall(ISLocalToGlobalMappingCreateIS(is, l2g));
576 PetscCall(ISDestroy(&is));
577
578 PetscCall(MatDestroy(&A3));
579 PetscCall(PetscFree(lndmapi));
580 MatPreallocateEnd(dnz, onz);
581 for (i = 0; i < gcnt; i++) PetscCall(ISDestroy(&workis[i]));
582 PetscCall(PetscFree(workis));
583 }
584 PetscCall(ISRestoreIndices(ndmap, &ndmapi));
585 PetscCall(PetscSFDestroy(&sf));
586 PetscCall(PetscFree(ndmapc));
587 PetscCall(ISDestroy(&ndmap));
588 PetscCall(ISDestroy(&ndsub));
589 PetscCall(ISLocalToGlobalMappingSetBlockSize(*l2g, bs));
590 PetscCall(ISLocalToGlobalMappingViewFromOptions(*l2g, NULL, "-mat_is_nd_l2g_view"));
591 break;
592 case MAT_IS_DISASSEMBLE_L2G_NATURAL:
593 PetscCall(PetscObjectQuery((PetscObject)A, "__PETSc_dm", &dm));
594 if (dm) { /* if a matrix comes from a DM, most likely we can use the l2gmap if any */
595 PetscCall(MatGetLocalToGlobalMapping(A, l2g, NULL));
596 PetscCall(PetscObjectReference((PetscObject)*l2g));
597 if (*l2g) PetscFunctionReturn(PETSC_SUCCESS);
598 }
599 if (ismpiaij) {
600 PetscCall(MatMPIAIJGetSeqAIJ(A, &Ad, &Ao, &garray));
601 } else if (ismpibaij) {
602 PetscCall(MatMPIBAIJGetSeqBAIJ(A, &Ad, &Ao, &garray));
603 } else SETERRQ(comm, PETSC_ERR_SUP, "Type %s", ((PetscObject)A)->type_name);
604 if (A->rmap->n) {
605 PetscInt dc, oc, stc, *aux;
606
607 PetscCall(MatGetLocalSize(Ad, NULL, &dc));
608 PetscCall(MatGetLocalSize(Ao, NULL, &oc));
609 PetscCheck(!oc || garray, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "garray not present");
610 PetscCall(MatGetOwnershipRangeColumn(A, &stc, NULL));
611 PetscCall(PetscMalloc1((dc + oc) / bs, &aux));
612 for (i = 0; i < dc / bs; i++) aux[i] = i + stc / bs;
613 for (i = 0; i < oc / bs; i++) aux[i + dc / bs] = (ismpiaij ? garray[i * bs] / bs : garray[i]);
614 PetscCall(ISCreateBlock(comm, bs, (dc + oc) / bs, aux, PETSC_OWN_POINTER, &is));
615 } else {
616 PetscCall(ISCreateBlock(comm, 1, 0, NULL, PETSC_OWN_POINTER, &is));
617 }
618 PetscCall(ISLocalToGlobalMappingCreateIS(is, l2g));
619 PetscCall(ISDestroy(&is));
620 break;
621 default:
622 SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unsupported l2g disassembling type %d", mode);
623 }
624 PetscFunctionReturn(PETSC_SUCCESS);
625 }
626
MatConvert_XAIJ_IS(Mat A,MatType type,MatReuse reuse,Mat * newmat)627 PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat A, MatType type, MatReuse reuse, Mat *newmat)
628 {
629 Mat lA, Ad, Ao, B = NULL;
630 ISLocalToGlobalMapping rl2g, cl2g;
631 IS is;
632 MPI_Comm comm;
633 void *ptrs[2];
634 const char *names[2] = {"_convert_csr_aux", "_convert_csr_data"};
635 const PetscInt *garray;
636 PetscScalar *dd, *od, *aa, *data;
637 const PetscInt *di, *dj, *oi, *oj;
638 const PetscInt *odi, *odj, *ooi, *ooj;
639 PetscInt *aux, *ii, *jj;
640 PetscInt rbs, cbs, lc, dr, dc, oc, str, stc, nnz, i, jd, jo, cum;
641 PetscBool flg, ismpiaij, ismpibaij, was_inplace = PETSC_FALSE, cong;
642 PetscMPIInt size;
643
644 PetscFunctionBegin;
645 PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
646 PetscCallMPI(MPI_Comm_size(comm, &size));
647 if (size == 1) {
648 PetscCall(MatConvert_SeqXAIJ_IS(A, type, reuse, newmat));
649 PetscFunctionReturn(PETSC_SUCCESS);
650 }
651 PetscCall(MatGetBlockSizes(A, &rbs, &cbs));
652 PetscCall(MatHasCongruentLayouts(A, &cong));
653 if (reuse != MAT_REUSE_MATRIX && cong && rbs == cbs) {
654 PetscCall(MatMPIXAIJComputeLocalToGlobalMapping_Private(A, &rl2g));
655 PetscCall(MatCreate(comm, &B));
656 PetscCall(MatSetType(B, MATIS));
657 PetscCall(MatSetSizes(B, A->rmap->n, A->rmap->n, A->rmap->N, A->rmap->N));
658 PetscCall(MatSetLocalToGlobalMapping(B, rl2g, rl2g));
659 PetscCall(MatSetBlockSizes(B, rbs, rbs));
660 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
661 if (reuse == MAT_INPLACE_MATRIX) was_inplace = PETSC_TRUE;
662 reuse = MAT_REUSE_MATRIX;
663 }
664 if (reuse == MAT_REUSE_MATRIX) {
665 Mat *newlA, lA;
666 IS rows, cols;
667 const PetscInt *ridx, *cidx;
668 PetscInt nr, nc;
669
670 if (!B) B = *newmat;
671 PetscCall(MatISGetLocalToGlobalMapping(B, &rl2g, &cl2g));
672 PetscCall(ISLocalToGlobalMappingGetBlockIndices(rl2g, &ridx));
673 PetscCall(ISLocalToGlobalMappingGetBlockIndices(cl2g, &cidx));
674 PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &nr));
675 PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &nc));
676 PetscCall(ISLocalToGlobalMappingGetBlockSize(rl2g, &rbs));
677 PetscCall(ISLocalToGlobalMappingGetBlockSize(cl2g, &cbs));
678 PetscCall(ISCreateBlock(comm, rbs, nr / rbs, ridx, PETSC_USE_POINTER, &rows));
679 if (rl2g != cl2g) {
680 PetscCall(ISCreateBlock(comm, cbs, nc / cbs, cidx, PETSC_USE_POINTER, &cols));
681 } else {
682 PetscCall(PetscObjectReference((PetscObject)rows));
683 cols = rows;
684 }
685 PetscCall(MatISGetLocalMat(B, &lA));
686 PetscCall(MatCreateSubMatrices(A, 1, &rows, &cols, MAT_INITIAL_MATRIX, &newlA));
687 PetscCall(MatConvert(newlA[0], MATSEQAIJ, MAT_INPLACE_MATRIX, &newlA[0]));
688 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(rl2g, &ridx));
689 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(cl2g, &cidx));
690 PetscCall(ISDestroy(&rows));
691 PetscCall(ISDestroy(&cols));
692 if (!lA->preallocated) { /* first time */
693 PetscCall(MatDuplicate(newlA[0], MAT_COPY_VALUES, &lA));
694 PetscCall(MatISSetLocalMat(B, lA));
695 PetscCall(PetscObjectDereference((PetscObject)lA));
696 }
697 PetscCall(MatCopy(newlA[0], lA, SAME_NONZERO_PATTERN));
698 PetscCall(MatDestroySubMatrices(1, &newlA));
699 PetscCall(MatISScaleDisassembling_Private(B));
700 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
701 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
702 if (was_inplace) PetscCall(MatHeaderReplace(A, &B));
703 else *newmat = B;
704 PetscFunctionReturn(PETSC_SUCCESS);
705 }
706 /* general case, just compress out the column space */
707 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
708 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIBAIJ, &ismpibaij));
709 if (ismpiaij) {
710 cbs = 1; /* We cannot guarantee the off-process matrix will respect the column block size */
711 PetscCall(MatMPIAIJGetSeqAIJ(A, &Ad, &Ao, &garray));
712 } else if (ismpibaij) {
713 PetscCall(MatMPIBAIJGetSeqBAIJ(A, &Ad, &Ao, &garray));
714 PetscCall(MatConvert(Ad, MATSEQAIJ, MAT_INITIAL_MATRIX, &Ad));
715 PetscCall(MatConvert(Ao, MATSEQAIJ, MAT_INITIAL_MATRIX, &Ao));
716 } else SETERRQ(comm, PETSC_ERR_SUP, "Type %s", ((PetscObject)A)->type_name);
717 PetscCall(MatSeqAIJGetArray(Ad, &dd));
718 PetscCall(MatSeqAIJGetArray(Ao, &od));
719
720 /* access relevant information from MPIAIJ */
721 PetscCall(MatGetOwnershipRange(A, &str, NULL));
722 PetscCall(MatGetOwnershipRangeColumn(A, &stc, NULL));
723 PetscCall(MatGetLocalSize(Ad, &dr, &dc));
724 PetscCall(MatGetLocalSize(Ao, NULL, &oc));
725 PetscCheck(!oc || garray, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "garray not present");
726
727 PetscCall(MatGetRowIJ(Ad, 0, PETSC_FALSE, PETSC_FALSE, &i, &di, &dj, &flg));
728 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
729 PetscCall(MatGetRowIJ(Ao, 0, PETSC_FALSE, PETSC_FALSE, &i, &oi, &oj, &flg));
730 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
731 nnz = di[dr] + oi[dr];
732 /* store original pointers to be restored later */
733 odi = di;
734 odj = dj;
735 ooi = oi;
736 ooj = oj;
737
738 /* generate l2g maps for rows and cols */
739 PetscCall(ISCreateStride(comm, dr / rbs, str / rbs, 1, &is));
740 if (rbs > 1) {
741 IS is2;
742
743 PetscCall(ISGetLocalSize(is, &i));
744 PetscCall(ISGetIndices(is, (const PetscInt **)&aux));
745 PetscCall(ISCreateBlock(comm, rbs, i, aux, PETSC_COPY_VALUES, &is2));
746 PetscCall(ISRestoreIndices(is, (const PetscInt **)&aux));
747 PetscCall(ISDestroy(&is));
748 is = is2;
749 }
750 PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
751 PetscCall(ISDestroy(&is));
752 if (dr) {
753 PetscCall(PetscMalloc1((dc + oc) / cbs, &aux));
754 for (i = 0; i < dc / cbs; i++) aux[i] = i + stc / cbs;
755 for (i = 0; i < oc / cbs; i++) aux[i + dc / cbs] = garray[i];
756 PetscCall(ISCreateBlock(comm, cbs, (dc + oc) / cbs, aux, PETSC_OWN_POINTER, &is));
757 lc = dc + oc;
758 } else {
759 PetscCall(ISCreateBlock(comm, cbs, 0, NULL, PETSC_OWN_POINTER, &is));
760 lc = 0;
761 }
762 PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
763 PetscCall(ISDestroy(&is));
764
765 /* create MATIS object */
766 PetscCall(MatCreate(comm, &B));
767 PetscCall(MatSetSizes(B, dr, dc, PETSC_DECIDE, PETSC_DECIDE));
768 PetscCall(MatSetType(B, MATIS));
769 PetscCall(MatSetBlockSizes(B, rbs, cbs));
770 PetscCall(MatSetLocalToGlobalMapping(B, rl2g, cl2g));
771 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
772 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
773
774 /* merge local matrices */
775 PetscCall(PetscMalloc1(nnz + dr + 1, &aux));
776 PetscCall(PetscMalloc1(nnz, &data));
777 ii = aux;
778 jj = aux + dr + 1;
779 aa = data;
780 *ii = *(di++) + *(oi++);
781 for (jd = 0, jo = 0, cum = 0; *ii < nnz; cum++) {
782 for (; jd < *di; jd++) {
783 *jj++ = *dj++;
784 *aa++ = *dd++;
785 }
786 for (; jo < *oi; jo++) {
787 *jj++ = *oj++ + dc;
788 *aa++ = *od++;
789 }
790 *(++ii) = *(di++) + *(oi++);
791 }
792 for (; cum < dr; cum++) *(++ii) = nnz;
793
794 PetscCall(MatRestoreRowIJ(Ad, 0, PETSC_FALSE, PETSC_FALSE, &i, &odi, &odj, &flg));
795 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot restore IJ structure");
796 PetscCall(MatRestoreRowIJ(Ao, 0, PETSC_FALSE, PETSC_FALSE, &i, &ooi, &ooj, &flg));
797 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot restore IJ structure");
798 PetscCall(MatSeqAIJRestoreArray(Ad, &dd));
799 PetscCall(MatSeqAIJRestoreArray(Ao, &od));
800
801 ii = aux;
802 jj = aux + dr + 1;
803 aa = data;
804 PetscCall(MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, dr, lc, ii, jj, aa, &lA));
805
806 /* create containers to destroy the data */
807 ptrs[0] = aux;
808 ptrs[1] = data;
809 for (i = 0; i < 2; i++) PetscCall(PetscObjectContainerCompose((PetscObject)lA, names[i], ptrs[i], PetscCtxDestroyDefault));
810 if (ismpibaij) { /* destroy converted local matrices */
811 PetscCall(MatDestroy(&Ad));
812 PetscCall(MatDestroy(&Ao));
813 }
814
815 /* finalize matrix */
816 PetscCall(MatISSetLocalMat(B, lA));
817 PetscCall(MatDestroy(&lA));
818 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
819 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
820 if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B));
821 else *newmat = B;
822 PetscFunctionReturn(PETSC_SUCCESS);
823 }
824
MatConvert_Nest_IS(Mat A,MatType type,MatReuse reuse,Mat * newmat)825 PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat A, MatType type, MatReuse reuse, Mat *newmat)
826 {
827 Mat **nest, *snest, **rnest, lA, B;
828 IS *iscol, *isrow, *islrow, *islcol;
829 ISLocalToGlobalMapping rl2g, cl2g;
830 MPI_Comm comm;
831 PetscInt *lr, *lc, *l2gidxs;
832 PetscInt i, j, nr, nc, rbs, cbs;
833 PetscBool convert, lreuse, *istrans;
834 PetscBool3 allow_repeated = PETSC_BOOL3_UNKNOWN;
835
836 PetscFunctionBegin;
837 PetscCall(MatNestGetSubMats(A, &nr, &nc, &nest));
838 lreuse = PETSC_FALSE;
839 rnest = NULL;
840 if (reuse == MAT_REUSE_MATRIX) {
841 PetscBool ismatis, isnest;
842
843 PetscCall(PetscObjectTypeCompare((PetscObject)*newmat, MATIS, &ismatis));
844 PetscCheck(ismatis, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_USER, "Cannot reuse matrix of type %s", ((PetscObject)*newmat)->type_name);
845 PetscCall(MatISGetLocalMat(*newmat, &lA));
846 PetscCall(PetscObjectTypeCompare((PetscObject)lA, MATNEST, &isnest));
847 if (isnest) {
848 PetscCall(MatNestGetSubMats(lA, &i, &j, &rnest));
849 lreuse = (PetscBool)(i == nr && j == nc);
850 if (!lreuse) rnest = NULL;
851 }
852 }
853 PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
854 PetscCall(PetscCalloc2(nr, &lr, nc, &lc));
855 PetscCall(PetscCalloc6(nr, &isrow, nc, &iscol, nr, &islrow, nc, &islcol, nr * nc, &snest, nr * nc, &istrans));
856 PetscCall(MatNestGetISs(A, isrow, iscol));
857 for (i = 0; i < nr; i++) {
858 for (j = 0; j < nc; j++) {
859 PetscBool ismatis, sallow;
860 PetscInt l1, l2, lb1, lb2, ij = i * nc + j;
861
862 /* Null matrix pointers are allowed in MATNEST */
863 if (!nest[i][j]) continue;
864
865 /* Nested matrices should be of type MATIS */
866 PetscCall(PetscObjectTypeCompare((PetscObject)nest[i][j], MATTRANSPOSEVIRTUAL, &istrans[ij]));
867 if (istrans[ij]) {
868 Mat T, lT;
869 PetscCall(MatTransposeGetMat(nest[i][j], &T));
870 PetscCall(PetscObjectTypeCompare((PetscObject)T, MATIS, &ismatis));
871 PetscCheck(ismatis, comm, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") (transposed) is not of type MATIS", i, j);
872 PetscCall(MatISGetAllowRepeated(T, &sallow));
873 PetscCall(MatISGetLocalMat(T, &lT));
874 PetscCall(MatCreateTranspose(lT, &snest[ij]));
875 } else {
876 PetscCall(PetscObjectTypeCompare((PetscObject)nest[i][j], MATIS, &ismatis));
877 PetscCheck(ismatis, comm, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") is not of type MATIS", i, j);
878 PetscCall(MatISGetAllowRepeated(nest[i][j], &sallow));
879 PetscCall(MatISGetLocalMat(nest[i][j], &snest[ij]));
880 }
881 if (allow_repeated == PETSC_BOOL3_UNKNOWN) allow_repeated = PetscBoolToBool3(sallow);
882 PetscCheck(sallow == PetscBool3ToBool(allow_repeated), comm, PETSC_ERR_SUP, "Cannot mix repeated and non repeated maps");
883
884 /* Check compatibility of local sizes */
885 PetscCall(MatGetSize(snest[ij], &l1, &l2));
886 PetscCall(MatGetBlockSizes(snest[ij], &lb1, &lb2));
887 if (!l1 || !l2) continue;
888 PetscCheck(!lr[i] || l1 == lr[i], PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid local size %" PetscInt_FMT " != %" PetscInt_FMT, i, j, lr[i], l1);
889 PetscCheck(!lc[j] || l2 == lc[j], PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid local size %" PetscInt_FMT " != %" PetscInt_FMT, i, j, lc[j], l2);
890 lr[i] = l1;
891 lc[j] = l2;
892
893 /* check compatibility for local matrix reusage */
894 if (rnest && !rnest[i][j] != !snest[ij]) lreuse = PETSC_FALSE;
895 }
896 }
897
898 if (PetscDefined(USE_DEBUG)) {
899 /* Check compatibility of l2g maps for rows */
900 for (i = 0; i < nr; i++) {
901 rl2g = NULL;
902 for (j = 0; j < nc; j++) {
903 PetscInt n1, n2;
904
905 if (!nest[i][j]) continue;
906 if (istrans[i * nc + j]) {
907 Mat T;
908
909 PetscCall(MatTransposeGetMat(nest[i][j], &T));
910 PetscCall(MatISGetLocalToGlobalMapping(T, NULL, &cl2g));
911 } else {
912 PetscCall(MatISGetLocalToGlobalMapping(nest[i][j], &cl2g, NULL));
913 }
914 PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &n1));
915 if (!n1) continue;
916 if (!rl2g) {
917 rl2g = cl2g;
918 } else {
919 const PetscInt *idxs1, *idxs2;
920 PetscBool same;
921
922 PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &n2));
923 PetscCheck(n1 == n2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid row l2gmap size %" PetscInt_FMT " != %" PetscInt_FMT, i, j, n1, n2);
924 PetscCall(ISLocalToGlobalMappingGetIndices(cl2g, &idxs1));
925 PetscCall(ISLocalToGlobalMappingGetIndices(rl2g, &idxs2));
926 PetscCall(PetscArraycmp(idxs1, idxs2, n1, &same));
927 PetscCall(ISLocalToGlobalMappingRestoreIndices(cl2g, &idxs1));
928 PetscCall(ISLocalToGlobalMappingRestoreIndices(rl2g, &idxs2));
929 PetscCheck(same, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid row l2gmap", i, j);
930 }
931 }
932 }
933 /* Check compatibility of l2g maps for columns */
934 for (i = 0; i < nc; i++) {
935 rl2g = NULL;
936 for (j = 0; j < nr; j++) {
937 PetscInt n1, n2;
938
939 if (!nest[j][i]) continue;
940 if (istrans[j * nc + i]) {
941 Mat T;
942
943 PetscCall(MatTransposeGetMat(nest[j][i], &T));
944 PetscCall(MatISGetLocalToGlobalMapping(T, &cl2g, NULL));
945 } else {
946 PetscCall(MatISGetLocalToGlobalMapping(nest[j][i], NULL, &cl2g));
947 }
948 PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &n1));
949 if (!n1) continue;
950 if (!rl2g) {
951 rl2g = cl2g;
952 } else {
953 const PetscInt *idxs1, *idxs2;
954 PetscBool same;
955
956 PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &n2));
957 PetscCheck(n1 == n2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid column l2gmap size %" PetscInt_FMT " != %" PetscInt_FMT, j, i, n1, n2);
958 PetscCall(ISLocalToGlobalMappingGetIndices(cl2g, &idxs1));
959 PetscCall(ISLocalToGlobalMappingGetIndices(rl2g, &idxs2));
960 PetscCall(PetscArraycmp(idxs1, idxs2, n1, &same));
961 PetscCall(ISLocalToGlobalMappingRestoreIndices(cl2g, &idxs1));
962 PetscCall(ISLocalToGlobalMappingRestoreIndices(rl2g, &idxs2));
963 PetscCheck(same, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid column l2gmap", j, i);
964 }
965 }
966 }
967 }
968
969 B = NULL;
970 if (reuse != MAT_REUSE_MATRIX) {
971 PetscInt stl;
972
973 /* Create l2g map for the rows of the new matrix and index sets for the local MATNEST */
974 for (i = 0, stl = 0; i < nr; i++) stl += lr[i];
975 PetscCall(PetscMalloc1(stl, &l2gidxs));
976 for (i = 0, stl = 0; i < nr; i++) {
977 Mat usedmat;
978 Mat_IS *matis;
979 const PetscInt *idxs;
980
981 /* local IS for local NEST */
982 PetscCall(ISCreateStride(PETSC_COMM_SELF, lr[i], stl, 1, &islrow[i]));
983
984 /* l2gmap */
985 j = 0;
986 usedmat = nest[i][j];
987 while (!usedmat && j < nc - 1) usedmat = nest[i][++j];
988 PetscCheck(usedmat, comm, PETSC_ERR_SUP, "Cannot find valid row mat");
989
990 if (istrans[i * nc + j]) {
991 Mat T;
992 PetscCall(MatTransposeGetMat(usedmat, &T));
993 usedmat = T;
994 }
995 matis = (Mat_IS *)usedmat->data;
996 PetscCall(ISGetIndices(isrow[i], &idxs));
997 if (istrans[i * nc + j]) {
998 PetscCall(PetscSFBcastBegin(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
999 PetscCall(PetscSFBcastEnd(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
1000 } else {
1001 PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
1002 PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
1003 }
1004 PetscCall(ISRestoreIndices(isrow[i], &idxs));
1005 stl += lr[i];
1006 }
1007 PetscCall(ISLocalToGlobalMappingCreate(comm, 1, stl, l2gidxs, PETSC_OWN_POINTER, &rl2g));
1008
1009 /* Create l2g map for columns of the new matrix and index sets for the local MATNEST */
1010 for (i = 0, stl = 0; i < nc; i++) stl += lc[i];
1011 PetscCall(PetscMalloc1(stl, &l2gidxs));
1012 for (i = 0, stl = 0; i < nc; i++) {
1013 Mat usedmat;
1014 Mat_IS *matis;
1015 const PetscInt *idxs;
1016
1017 /* local IS for local NEST */
1018 PetscCall(ISCreateStride(PETSC_COMM_SELF, lc[i], stl, 1, &islcol[i]));
1019
1020 /* l2gmap */
1021 j = 0;
1022 usedmat = nest[j][i];
1023 while (!usedmat && j < nr - 1) usedmat = nest[++j][i];
1024 PetscCheck(usedmat, comm, PETSC_ERR_SUP, "Cannot find valid column mat");
1025 if (istrans[j * nc + i]) {
1026 Mat T;
1027 PetscCall(MatTransposeGetMat(usedmat, &T));
1028 usedmat = T;
1029 }
1030 matis = (Mat_IS *)usedmat->data;
1031 PetscCall(ISGetIndices(iscol[i], &idxs));
1032 if (istrans[j * nc + i]) {
1033 PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
1034 PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
1035 } else {
1036 PetscCall(PetscSFBcastBegin(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
1037 PetscCall(PetscSFBcastEnd(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
1038 }
1039 PetscCall(ISRestoreIndices(iscol[i], &idxs));
1040 stl += lc[i];
1041 }
1042 PetscCall(ISLocalToGlobalMappingCreate(comm, 1, stl, l2gidxs, PETSC_OWN_POINTER, &cl2g));
1043
1044 /* Create MATIS */
1045 PetscCall(MatCreate(comm, &B));
1046 PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
1047 PetscCall(MatGetBlockSizes(A, &rbs, &cbs));
1048 PetscCall(MatSetBlockSizes(B, rbs, cbs));
1049 PetscCall(MatSetType(B, MATIS));
1050 PetscCall(MatISSetLocalMatType(B, MATNEST));
1051 PetscCall(MatISSetAllowRepeated(B, PetscBool3ToBool(allow_repeated)));
1052 { /* hack : avoid setup of scatters */
1053 Mat_IS *matis = (Mat_IS *)B->data;
1054 matis->islocalref = B;
1055 }
1056 PetscCall(MatSetLocalToGlobalMapping(B, rl2g, cl2g));
1057 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
1058 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
1059 PetscCall(MatCreateNest(PETSC_COMM_SELF, nr, islrow, nc, islcol, snest, &lA));
1060 PetscCall(MatNestSetVecType(lA, VECNEST));
1061 for (i = 0; i < nr * nc; i++) {
1062 if (istrans[i]) PetscCall(MatDestroy(&snest[i]));
1063 }
1064 PetscCall(MatISSetLocalMat(B, lA));
1065 PetscCall(MatDestroy(&lA));
1066 { /* hack : setup of scatters done here */
1067 Mat_IS *matis = (Mat_IS *)B->data;
1068
1069 matis->islocalref = NULL;
1070 PetscCall(MatISSetUpScatters_Private(B));
1071 }
1072 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
1073 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
1074 if (reuse == MAT_INPLACE_MATRIX) {
1075 PetscCall(MatHeaderReplace(A, &B));
1076 } else {
1077 *newmat = B;
1078 }
1079 } else {
1080 if (lreuse) {
1081 PetscCall(MatISGetLocalMat(*newmat, &lA));
1082 for (i = 0; i < nr; i++) {
1083 for (j = 0; j < nc; j++) {
1084 if (snest[i * nc + j]) {
1085 PetscCall(MatNestSetSubMat(lA, i, j, snest[i * nc + j]));
1086 if (istrans[i * nc + j]) PetscCall(MatDestroy(&snest[i * nc + j]));
1087 }
1088 }
1089 }
1090 } else {
1091 PetscInt stl;
1092 for (i = 0, stl = 0; i < nr; i++) {
1093 PetscCall(ISCreateStride(PETSC_COMM_SELF, lr[i], stl, 1, &islrow[i]));
1094 stl += lr[i];
1095 }
1096 for (i = 0, stl = 0; i < nc; i++) {
1097 PetscCall(ISCreateStride(PETSC_COMM_SELF, lc[i], stl, 1, &islcol[i]));
1098 stl += lc[i];
1099 }
1100 PetscCall(MatCreateNest(PETSC_COMM_SELF, nr, islrow, nc, islcol, snest, &lA));
1101 for (i = 0; i < nr * nc; i++) {
1102 if (istrans[i]) PetscCall(MatDestroy(&snest[i]));
1103 }
1104 PetscCall(MatISSetLocalMat(*newmat, lA));
1105 PetscCall(MatDestroy(&lA));
1106 }
1107 PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY));
1108 PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY));
1109 }
1110
1111 /* Create local matrix in MATNEST format */
1112 convert = PETSC_FALSE;
1113 PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)A)->prefix, "-mat_is_convert_local_nest", &convert, NULL));
1114 if (convert) {
1115 Mat M;
1116 MatISLocalFields lf;
1117
1118 PetscCall(MatISGetLocalMat(*newmat, &lA));
1119 PetscCall(MatConvert(lA, MATAIJ, MAT_INITIAL_MATRIX, &M));
1120 PetscCall(MatISSetLocalMat(*newmat, M));
1121 PetscCall(MatDestroy(&M));
1122
1123 /* attach local fields to the matrix */
1124 PetscCall(PetscNew(&lf));
1125 PetscCall(PetscMalloc2(nr, &lf->rf, nc, &lf->cf));
1126 for (i = 0; i < nr; i++) {
1127 PetscInt n, st;
1128
1129 PetscCall(ISGetLocalSize(islrow[i], &n));
1130 PetscCall(ISStrideGetInfo(islrow[i], &st, NULL));
1131 PetscCall(ISCreateStride(comm, n, st, 1, &lf->rf[i]));
1132 }
1133 for (i = 0; i < nc; i++) {
1134 PetscInt n, st;
1135
1136 PetscCall(ISGetLocalSize(islcol[i], &n));
1137 PetscCall(ISStrideGetInfo(islcol[i], &st, NULL));
1138 PetscCall(ISCreateStride(comm, n, st, 1, &lf->cf[i]));
1139 }
1140 lf->nr = nr;
1141 lf->nc = nc;
1142 PetscCall(PetscObjectContainerCompose((PetscObject)*newmat, "_convert_nest_lfields", lf, MatISContainerDestroyFields_Private));
1143 }
1144
1145 /* Free workspace */
1146 for (i = 0; i < nr; i++) PetscCall(ISDestroy(&islrow[i]));
1147 for (i = 0; i < nc; i++) PetscCall(ISDestroy(&islcol[i]));
1148 PetscCall(PetscFree6(isrow, iscol, islrow, islcol, snest, istrans));
1149 PetscCall(PetscFree2(lr, lc));
1150 PetscFunctionReturn(PETSC_SUCCESS);
1151 }
1152
MatDiagonalScale_IS(Mat A,Vec l,Vec r)1153 static PetscErrorCode MatDiagonalScale_IS(Mat A, Vec l, Vec r)
1154 {
1155 Mat_IS *matis = (Mat_IS *)A->data;
1156 Vec ll, rr;
1157 const PetscScalar *Y, *X;
1158 PetscScalar *x, *y;
1159
1160 PetscFunctionBegin;
1161 if (l) {
1162 ll = matis->y;
1163 PetscCall(VecGetArrayRead(l, &Y));
1164 PetscCall(VecGetArray(ll, &y));
1165 PetscCall(PetscSFBcastBegin(matis->sf, MPIU_SCALAR, Y, y, MPI_REPLACE));
1166 } else {
1167 ll = NULL;
1168 }
1169 if (r) {
1170 rr = matis->x;
1171 PetscCall(VecGetArrayRead(r, &X));
1172 PetscCall(VecGetArray(rr, &x));
1173 PetscCall(PetscSFBcastBegin(matis->csf, MPIU_SCALAR, X, x, MPI_REPLACE));
1174 } else {
1175 rr = NULL;
1176 }
1177 if (ll) {
1178 PetscCall(PetscSFBcastEnd(matis->sf, MPIU_SCALAR, Y, y, MPI_REPLACE));
1179 PetscCall(VecRestoreArrayRead(l, &Y));
1180 PetscCall(VecRestoreArray(ll, &y));
1181 }
1182 if (rr) {
1183 PetscCall(PetscSFBcastEnd(matis->csf, MPIU_SCALAR, X, x, MPI_REPLACE));
1184 PetscCall(VecRestoreArrayRead(r, &X));
1185 PetscCall(VecRestoreArray(rr, &x));
1186 }
1187 PetscCall(MatDiagonalScale(matis->A, ll, rr));
1188 PetscFunctionReturn(PETSC_SUCCESS);
1189 }
1190
MatGetInfo_IS(Mat A,MatInfoType flag,MatInfo * ginfo)1191 static PetscErrorCode MatGetInfo_IS(Mat A, MatInfoType flag, MatInfo *ginfo)
1192 {
1193 Mat_IS *matis = (Mat_IS *)A->data;
1194 MatInfo info;
1195 PetscLogDouble isend[6], irecv[6];
1196 PetscInt bs;
1197
1198 PetscFunctionBegin;
1199 PetscCall(MatGetBlockSize(A, &bs));
1200 if (matis->A->ops->getinfo) {
1201 PetscCall(MatGetInfo(matis->A, MAT_LOCAL, &info));
1202 isend[0] = info.nz_used;
1203 isend[1] = info.nz_allocated;
1204 isend[2] = info.nz_unneeded;
1205 isend[3] = info.memory;
1206 isend[4] = info.mallocs;
1207 } else {
1208 isend[0] = 0.;
1209 isend[1] = 0.;
1210 isend[2] = 0.;
1211 isend[3] = 0.;
1212 isend[4] = 0.;
1213 }
1214 isend[5] = matis->A->num_ass;
1215 if (flag == MAT_LOCAL) {
1216 ginfo->nz_used = isend[0];
1217 ginfo->nz_allocated = isend[1];
1218 ginfo->nz_unneeded = isend[2];
1219 ginfo->memory = isend[3];
1220 ginfo->mallocs = isend[4];
1221 ginfo->assemblies = isend[5];
1222 } else if (flag == MAT_GLOBAL_MAX) {
1223 PetscCallMPI(MPIU_Allreduce(isend, irecv, 6, MPIU_PETSCLOGDOUBLE, MPI_MAX, PetscObjectComm((PetscObject)A)));
1224
1225 ginfo->nz_used = irecv[0];
1226 ginfo->nz_allocated = irecv[1];
1227 ginfo->nz_unneeded = irecv[2];
1228 ginfo->memory = irecv[3];
1229 ginfo->mallocs = irecv[4];
1230 ginfo->assemblies = irecv[5];
1231 } else if (flag == MAT_GLOBAL_SUM) {
1232 PetscCallMPI(MPIU_Allreduce(isend, irecv, 5, MPIU_PETSCLOGDOUBLE, MPI_SUM, PetscObjectComm((PetscObject)A)));
1233
1234 ginfo->nz_used = irecv[0];
1235 ginfo->nz_allocated = irecv[1];
1236 ginfo->nz_unneeded = irecv[2];
1237 ginfo->memory = irecv[3];
1238 ginfo->mallocs = irecv[4];
1239 ginfo->assemblies = A->num_ass;
1240 }
1241 ginfo->block_size = bs;
1242 ginfo->fill_ratio_given = 0;
1243 ginfo->fill_ratio_needed = 0;
1244 ginfo->factor_mallocs = 0;
1245 PetscFunctionReturn(PETSC_SUCCESS);
1246 }
1247
MatTranspose_IS(Mat A,MatReuse reuse,Mat * B)1248 static PetscErrorCode MatTranspose_IS(Mat A, MatReuse reuse, Mat *B)
1249 {
1250 Mat C, lC, lA;
1251
1252 PetscFunctionBegin;
1253 if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B));
1254 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
1255 ISLocalToGlobalMapping rl2g, cl2g;
1256 PetscBool allow_repeated;
1257
1258 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
1259 PetscCall(MatSetSizes(C, A->cmap->n, A->rmap->n, A->cmap->N, A->rmap->N));
1260 PetscCall(MatSetBlockSizes(C, A->cmap->bs, A->rmap->bs));
1261 PetscCall(MatSetType(C, MATIS));
1262 PetscCall(MatISGetAllowRepeated(A, &allow_repeated));
1263 PetscCall(MatISSetAllowRepeated(C, allow_repeated));
1264 PetscCall(MatGetLocalToGlobalMapping(A, &rl2g, &cl2g));
1265 PetscCall(MatSetLocalToGlobalMapping(C, cl2g, rl2g));
1266 } else C = *B;
1267
1268 /* perform local transposition */
1269 PetscCall(MatISGetLocalMat(A, &lA));
1270 PetscCall(MatTranspose(lA, MAT_INITIAL_MATRIX, &lC));
1271 PetscCall(MatSetLocalToGlobalMapping(lC, lA->cmap->mapping, lA->rmap->mapping));
1272 PetscCall(MatISSetLocalMat(C, lC));
1273 PetscCall(MatDestroy(&lC));
1274
1275 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
1276 *B = C;
1277 } else {
1278 PetscCall(MatHeaderMerge(A, &C));
1279 }
1280 PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY));
1281 PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY));
1282 PetscFunctionReturn(PETSC_SUCCESS);
1283 }
1284
MatDiagonalSet_IS(Mat A,Vec D,InsertMode insmode)1285 static PetscErrorCode MatDiagonalSet_IS(Mat A, Vec D, InsertMode insmode)
1286 {
1287 Mat_IS *is = (Mat_IS *)A->data;
1288
1289 PetscFunctionBegin;
1290 PetscCheck(!is->allow_repeated || insmode == ADD_VALUES, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "INSERT_VALUES with repeated entries not supported");
1291 if (D) { /* MatShift_IS pass D = NULL */
1292 PetscCall(VecScatterBegin(is->rctx, D, is->y, INSERT_VALUES, SCATTER_FORWARD));
1293 PetscCall(VecScatterEnd(is->rctx, D, is->y, INSERT_VALUES, SCATTER_FORWARD));
1294 }
1295 PetscCall(VecPointwiseDivide(is->y, is->y, is->counter));
1296 PetscCall(MatDiagonalSet(is->A, is->y, insmode));
1297 PetscFunctionReturn(PETSC_SUCCESS);
1298 }
1299
MatShift_IS(Mat A,PetscScalar a)1300 static PetscErrorCode MatShift_IS(Mat A, PetscScalar a)
1301 {
1302 Mat_IS *is = (Mat_IS *)A->data;
1303
1304 PetscFunctionBegin;
1305 PetscCall(VecSet(is->y, a));
1306 PetscCall(MatDiagonalSet_IS(A, NULL, ADD_VALUES));
1307 PetscFunctionReturn(PETSC_SUCCESS);
1308 }
1309
MatSetValuesLocal_SubMat_IS(Mat A,PetscInt m,const PetscInt * rows,PetscInt n,const PetscInt * cols,const PetscScalar * values,InsertMode addv)1310 static PetscErrorCode MatSetValuesLocal_SubMat_IS(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
1311 {
1312 PetscInt buf[2 * MATIS_MAX_ENTRIES_INSERTION], *rows_l = NULL, *cols_l = NULL;
1313
1314 PetscFunctionBegin;
1315 IndexSpaceGet(buf, m, n, rows_l, cols_l);
1316 PetscCall(ISLocalToGlobalMappingApply(A->rmap->mapping, m, rows, rows_l));
1317 PetscCall(ISLocalToGlobalMappingApply(A->cmap->mapping, n, cols, cols_l));
1318 PetscCall(MatSetValuesLocal_IS(A, m, rows_l, n, cols_l, values, addv));
1319 IndexSpaceRestore(buf, m, n, rows_l, cols_l);
1320 PetscFunctionReturn(PETSC_SUCCESS);
1321 }
1322
MatSetValuesBlockedLocal_SubMat_IS(Mat A,PetscInt m,const PetscInt * rows,PetscInt n,const PetscInt * cols,const PetscScalar * values,InsertMode addv)1323 static PetscErrorCode MatSetValuesBlockedLocal_SubMat_IS(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
1324 {
1325 PetscInt buf[2 * MATIS_MAX_ENTRIES_INSERTION], *rows_l = NULL, *cols_l = NULL, rbs, cbs;
1326
1327 PetscFunctionBegin;
1328 /* We cannot guarantee the local matrix will have the same block size of the original matrix */
1329 PetscCall(ISLocalToGlobalMappingGetBlockSize(A->rmap->mapping, &rbs));
1330 PetscCall(ISLocalToGlobalMappingGetBlockSize(A->cmap->mapping, &cbs));
1331 IndexSpaceGet(buf, m * rbs, n * cbs, rows_l, cols_l);
1332 BlockIndicesExpand(m, rows, rbs, rows_l);
1333 BlockIndicesExpand(n, cols, cbs, cols_l);
1334 PetscCall(ISLocalToGlobalMappingApply(A->rmap->mapping, m * rbs, rows_l, rows_l));
1335 PetscCall(ISLocalToGlobalMappingApply(A->cmap->mapping, n * cbs, cols_l, cols_l));
1336 PetscCall(MatSetValuesLocal_IS(A, m * rbs, rows_l, n * cbs, cols_l, values, addv));
1337 IndexSpaceRestore(buf, m * rbs, n * cbs, rows_l, cols_l);
1338 PetscFunctionReturn(PETSC_SUCCESS);
1339 }
1340
MatZeroRowsLocal_SubMat_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)1341 static PetscErrorCode MatZeroRowsLocal_SubMat_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
1342 {
1343 PetscInt *rows_l;
1344 Mat_IS *is = (Mat_IS *)A->data;
1345
1346 PetscFunctionBegin;
1347 PetscCall(PetscMalloc1(n, &rows_l));
1348 PetscCall(ISLocalToGlobalMappingApply(A->rmap->mapping, n, rows, rows_l));
1349 PetscCall(MatZeroRowsLocal(is->islocalref, n, rows_l, diag, x, b));
1350 PetscCall(PetscFree(rows_l));
1351 PetscFunctionReturn(PETSC_SUCCESS);
1352 }
1353
MatZeroRowsColumnsLocal_SubMat_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)1354 static PetscErrorCode MatZeroRowsColumnsLocal_SubMat_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
1355 {
1356 PetscInt *rows_l;
1357 Mat_IS *is = (Mat_IS *)A->data;
1358
1359 PetscFunctionBegin;
1360 PetscCall(PetscMalloc1(n, &rows_l));
1361 PetscCall(ISLocalToGlobalMappingApply(A->rmap->mapping, n, rows, rows_l));
1362 PetscCall(MatZeroRowsColumnsLocal(is->islocalref, n, rows_l, diag, x, b));
1363 PetscCall(PetscFree(rows_l));
1364 PetscFunctionReturn(PETSC_SUCCESS);
1365 }
1366
MatCreateSubMatrix_IS(Mat mat,IS irow,IS icol,MatReuse scall,Mat * newmat)1367 static PetscErrorCode MatCreateSubMatrix_IS(Mat mat, IS irow, IS icol, MatReuse scall, Mat *newmat)
1368 {
1369 Mat locmat, newlocmat;
1370 Mat_IS *newmatis;
1371 const PetscInt *idxs;
1372 PetscInt i, m, n;
1373
1374 PetscFunctionBegin;
1375 if (scall == MAT_REUSE_MATRIX) {
1376 PetscBool ismatis;
1377
1378 PetscCall(PetscObjectTypeCompare((PetscObject)*newmat, MATIS, &ismatis));
1379 PetscCheck(ismatis, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_ARG_WRONG, "Cannot reuse matrix! Not of MATIS type");
1380 newmatis = (Mat_IS *)(*newmat)->data;
1381 PetscCheck(newmatis->getsub_ris, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_ARG_WRONG, "Cannot reuse matrix! Misses local row IS");
1382 PetscCheck(newmatis->getsub_cis, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_ARG_WRONG, "Cannot reuse matrix! Misses local col IS");
1383 }
1384 /* irow and icol may not have duplicate entries */
1385 if (PetscDefined(USE_DEBUG)) {
1386 Vec rtest, ltest;
1387 const PetscScalar *array;
1388
1389 PetscCall(MatCreateVecs(mat, <est, &rtest));
1390 PetscCall(ISGetLocalSize(irow, &n));
1391 PetscCall(ISGetIndices(irow, &idxs));
1392 for (i = 0; i < n; i++) PetscCall(VecSetValue(rtest, idxs[i], 1.0, ADD_VALUES));
1393 PetscCall(VecAssemblyBegin(rtest));
1394 PetscCall(VecAssemblyEnd(rtest));
1395 PetscCall(VecGetLocalSize(rtest, &n));
1396 PetscCall(VecGetOwnershipRange(rtest, &m, NULL));
1397 PetscCall(VecGetArrayRead(rtest, &array));
1398 for (i = 0; i < n; i++) PetscCheck(array[i] == 0. || array[i] == 1., PETSC_COMM_SELF, PETSC_ERR_SUP, "Index %" PetscInt_FMT " counted %" PetscInt_FMT " times! Irow may not have duplicate entries", i + m, (PetscInt)PetscRealPart(array[i]));
1399 PetscCall(VecRestoreArrayRead(rtest, &array));
1400 PetscCall(ISRestoreIndices(irow, &idxs));
1401 PetscCall(ISGetLocalSize(icol, &n));
1402 PetscCall(ISGetIndices(icol, &idxs));
1403 for (i = 0; i < n; i++) PetscCall(VecSetValue(ltest, idxs[i], 1.0, ADD_VALUES));
1404 PetscCall(VecAssemblyBegin(ltest));
1405 PetscCall(VecAssemblyEnd(ltest));
1406 PetscCall(VecGetLocalSize(ltest, &n));
1407 PetscCall(VecGetOwnershipRange(ltest, &m, NULL));
1408 PetscCall(VecGetArrayRead(ltest, &array));
1409 for (i = 0; i < n; i++) PetscCheck(array[i] == 0. || array[i] == 1., PETSC_COMM_SELF, PETSC_ERR_SUP, "Index %" PetscInt_FMT " counted %" PetscInt_FMT " times! Icol may not have duplicate entries", i + m, (PetscInt)PetscRealPart(array[i]));
1410 PetscCall(VecRestoreArrayRead(ltest, &array));
1411 PetscCall(ISRestoreIndices(icol, &idxs));
1412 PetscCall(VecDestroy(&rtest));
1413 PetscCall(VecDestroy(<est));
1414 }
1415 if (scall == MAT_INITIAL_MATRIX) {
1416 Mat_IS *matis = (Mat_IS *)mat->data;
1417 ISLocalToGlobalMapping rl2g;
1418 IS is;
1419 PetscInt *lidxs, *lgidxs, *newgidxs;
1420 PetscInt ll, newloc, irbs, icbs, arbs, acbs, rbs, cbs;
1421 PetscBool cong;
1422 MPI_Comm comm;
1423
1424 PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
1425 PetscCall(MatGetBlockSizes(mat, &arbs, &acbs));
1426 PetscCall(ISGetBlockSize(irow, &irbs));
1427 PetscCall(ISGetBlockSize(icol, &icbs));
1428 rbs = arbs == irbs ? irbs : 1;
1429 cbs = acbs == icbs ? icbs : 1;
1430 PetscCall(ISGetLocalSize(irow, &m));
1431 PetscCall(ISGetLocalSize(icol, &n));
1432 PetscCall(MatCreate(comm, newmat));
1433 PetscCall(MatSetType(*newmat, MATIS));
1434 PetscCall(MatISSetAllowRepeated(*newmat, matis->allow_repeated));
1435 PetscCall(MatSetSizes(*newmat, m, n, PETSC_DECIDE, PETSC_DECIDE));
1436 PetscCall(MatSetBlockSizes(*newmat, rbs, cbs));
1437 /* communicate irow to their owners in the layout */
1438 PetscCall(ISGetIndices(irow, &idxs));
1439 PetscCall(PetscLayoutMapLocal(mat->rmap, m, idxs, &ll, &lidxs, &lgidxs));
1440 PetscCall(ISRestoreIndices(irow, &idxs));
1441 PetscCall(PetscArrayzero(matis->sf_rootdata, matis->sf->nroots));
1442 for (i = 0; i < ll; i++) matis->sf_rootdata[lidxs[i]] = lgidxs[i] + 1;
1443 PetscCall(PetscFree(lidxs));
1444 PetscCall(PetscFree(lgidxs));
1445 PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
1446 PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
1447 for (i = 0, newloc = 0; i < matis->sf->nleaves; i++)
1448 if (matis->sf_leafdata[i]) newloc++;
1449 PetscCall(PetscMalloc1(newloc, &newgidxs));
1450 PetscCall(PetscMalloc1(newloc, &lidxs));
1451 for (i = 0, newloc = 0; i < matis->sf->nleaves; i++)
1452 if (matis->sf_leafdata[i]) {
1453 lidxs[newloc] = i;
1454 newgidxs[newloc++] = matis->sf_leafdata[i] - 1;
1455 }
1456 PetscCall(ISCreateGeneral(comm, newloc, newgidxs, PETSC_OWN_POINTER, &is));
1457 PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
1458 PetscCall(ISLocalToGlobalMappingSetBlockSize(rl2g, rbs));
1459 PetscCall(ISDestroy(&is));
1460 /* local is to extract local submatrix */
1461 newmatis = (Mat_IS *)(*newmat)->data;
1462 PetscCall(ISCreateGeneral(comm, newloc, lidxs, PETSC_OWN_POINTER, &newmatis->getsub_ris));
1463 PetscCall(MatHasCongruentLayouts(mat, &cong));
1464 if (cong && irow == icol && matis->csf == matis->sf) {
1465 PetscCall(MatSetLocalToGlobalMapping(*newmat, rl2g, rl2g));
1466 PetscCall(PetscObjectReference((PetscObject)newmatis->getsub_ris));
1467 newmatis->getsub_cis = newmatis->getsub_ris;
1468 } else {
1469 ISLocalToGlobalMapping cl2g;
1470
1471 /* communicate icol to their owners in the layout */
1472 PetscCall(ISGetIndices(icol, &idxs));
1473 PetscCall(PetscLayoutMapLocal(mat->cmap, n, idxs, &ll, &lidxs, &lgidxs));
1474 PetscCall(ISRestoreIndices(icol, &idxs));
1475 PetscCall(PetscArrayzero(matis->csf_rootdata, matis->csf->nroots));
1476 for (i = 0; i < ll; i++) matis->csf_rootdata[lidxs[i]] = lgidxs[i] + 1;
1477 PetscCall(PetscFree(lidxs));
1478 PetscCall(PetscFree(lgidxs));
1479 PetscCall(PetscSFBcastBegin(matis->csf, MPIU_INT, matis->csf_rootdata, matis->csf_leafdata, MPI_REPLACE));
1480 PetscCall(PetscSFBcastEnd(matis->csf, MPIU_INT, matis->csf_rootdata, matis->csf_leafdata, MPI_REPLACE));
1481 for (i = 0, newloc = 0; i < matis->csf->nleaves; i++)
1482 if (matis->csf_leafdata[i]) newloc++;
1483 PetscCall(PetscMalloc1(newloc, &newgidxs));
1484 PetscCall(PetscMalloc1(newloc, &lidxs));
1485 for (i = 0, newloc = 0; i < matis->csf->nleaves; i++)
1486 if (matis->csf_leafdata[i]) {
1487 lidxs[newloc] = i;
1488 newgidxs[newloc++] = matis->csf_leafdata[i] - 1;
1489 }
1490 PetscCall(ISCreateGeneral(comm, newloc, newgidxs, PETSC_OWN_POINTER, &is));
1491 PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
1492 PetscCall(ISLocalToGlobalMappingSetBlockSize(cl2g, cbs));
1493 PetscCall(ISDestroy(&is));
1494 /* local is to extract local submatrix */
1495 PetscCall(ISCreateGeneral(comm, newloc, lidxs, PETSC_OWN_POINTER, &newmatis->getsub_cis));
1496 PetscCall(MatSetLocalToGlobalMapping(*newmat, rl2g, cl2g));
1497 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
1498 }
1499 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
1500 } else {
1501 PetscCall(MatISGetLocalMat(*newmat, &newlocmat));
1502 }
1503 PetscCall(MatISGetLocalMat(mat, &locmat));
1504 newmatis = (Mat_IS *)(*newmat)->data;
1505 PetscCall(MatCreateSubMatrix(locmat, newmatis->getsub_ris, newmatis->getsub_cis, scall, &newlocmat));
1506 if (scall == MAT_INITIAL_MATRIX) {
1507 PetscCall(MatISSetLocalMat(*newmat, newlocmat));
1508 PetscCall(MatDestroy(&newlocmat));
1509 }
1510 PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY));
1511 PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY));
1512 PetscFunctionReturn(PETSC_SUCCESS);
1513 }
1514
MatCopy_IS(Mat A,Mat B,MatStructure str)1515 static PetscErrorCode MatCopy_IS(Mat A, Mat B, MatStructure str)
1516 {
1517 Mat_IS *a = (Mat_IS *)A->data, *b;
1518 PetscBool ismatis;
1519
1520 PetscFunctionBegin;
1521 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATIS, &ismatis));
1522 PetscCheck(ismatis, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Need to be implemented");
1523 b = (Mat_IS *)B->data;
1524 PetscCall(MatCopy(a->A, b->A, str));
1525 PetscCall(PetscObjectStateIncrease((PetscObject)B));
1526 PetscFunctionReturn(PETSC_SUCCESS);
1527 }
1528
MatISSetUpSF_IS(Mat B)1529 static PetscErrorCode MatISSetUpSF_IS(Mat B)
1530 {
1531 Mat_IS *matis = (Mat_IS *)B->data;
1532 const PetscInt *gidxs;
1533 PetscInt nleaves;
1534
1535 PetscFunctionBegin;
1536 if (matis->sf) PetscFunctionReturn(PETSC_SUCCESS);
1537 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)B), &matis->sf));
1538 PetscCall(ISLocalToGlobalMappingGetIndices(matis->rmapping, &gidxs));
1539 PetscCall(ISLocalToGlobalMappingGetSize(matis->rmapping, &nleaves));
1540 PetscCall(PetscSFSetGraphLayout(matis->sf, B->rmap, nleaves, NULL, PETSC_OWN_POINTER, gidxs));
1541 PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->rmapping, &gidxs));
1542 PetscCall(PetscMalloc2(matis->sf->nroots, &matis->sf_rootdata, matis->sf->nleaves, &matis->sf_leafdata));
1543 if (matis->rmapping != matis->cmapping) { /* setup SF for columns */
1544 PetscCall(ISLocalToGlobalMappingGetSize(matis->cmapping, &nleaves));
1545 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)B), &matis->csf));
1546 PetscCall(ISLocalToGlobalMappingGetIndices(matis->cmapping, &gidxs));
1547 PetscCall(PetscSFSetGraphLayout(matis->csf, B->cmap, nleaves, NULL, PETSC_OWN_POINTER, gidxs));
1548 PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->cmapping, &gidxs));
1549 PetscCall(PetscMalloc2(matis->csf->nroots, &matis->csf_rootdata, matis->csf->nleaves, &matis->csf_leafdata));
1550 } else {
1551 matis->csf = matis->sf;
1552 matis->csf_leafdata = matis->sf_leafdata;
1553 matis->csf_rootdata = matis->sf_rootdata;
1554 }
1555 PetscFunctionReturn(PETSC_SUCCESS);
1556 }
1557
1558 /*@
1559 MatISGetAllowRepeated - Get the flag to allow repeated entries in the local to global map
1560
1561 Not Collective
1562
1563 Input Parameter:
1564 . A - the matrix
1565
1566 Output Parameter:
1567 . flg - the boolean flag
1568
1569 Level: intermediate
1570
1571 .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateIS()`, `MatSetLocalToGlobalMapping()`, `MatISSetAllowRepeated()`
1572 @*/
MatISGetAllowRepeated(Mat A,PetscBool * flg)1573 PetscErrorCode MatISGetAllowRepeated(Mat A, PetscBool *flg)
1574 {
1575 PetscBool ismatis;
1576
1577 PetscFunctionBegin;
1578 PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
1579 PetscAssertPointer(flg, 2);
1580 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATIS, &ismatis));
1581 PetscCheck(ismatis, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)A)->type_name);
1582 *flg = ((Mat_IS *)A->data)->allow_repeated;
1583 PetscFunctionReturn(PETSC_SUCCESS);
1584 }
1585
1586 /*@
1587 MatISSetAllowRepeated - Set the flag to allow repeated entries in the local to global map
1588
1589 Logically Collective
1590
1591 Input Parameters:
1592 + A - the matrix
1593 - flg - the boolean flag
1594
1595 Level: intermediate
1596
1597 Notes:
1598 The default value is `PETSC_FALSE`.
1599 When called AFTER calling `MatSetLocalToGlobalMapping()` it will recreate the local matrices
1600 if `flg` is different from the previously set value.
1601 Specifically, when `flg` is true it will just recreate the local matrices, while if
1602 `flg` is false will assemble the local matrices summing up repeated entries.
1603
1604 .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateIS()`, `MatSetLocalToGlobalMapping()`, `MatISGetAllowRepeated()`
1605 @*/
MatISSetAllowRepeated(Mat A,PetscBool flg)1606 PetscErrorCode MatISSetAllowRepeated(Mat A, PetscBool flg)
1607 {
1608 PetscFunctionBegin;
1609 PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
1610 PetscValidType(A, 1);
1611 PetscValidLogicalCollectiveBool(A, flg, 2);
1612 PetscTryMethod(A, "MatISSetAllowRepeated_C", (Mat, PetscBool), (A, flg));
1613 PetscFunctionReturn(PETSC_SUCCESS);
1614 }
1615
MatISSetAllowRepeated_IS(Mat A,PetscBool flg)1616 static PetscErrorCode MatISSetAllowRepeated_IS(Mat A, PetscBool flg)
1617 {
1618 Mat_IS *matis = (Mat_IS *)A->data;
1619 Mat lA = NULL;
1620 ISLocalToGlobalMapping lrmap, lcmap;
1621
1622 PetscFunctionBegin;
1623 if (flg == matis->allow_repeated) PetscFunctionReturn(PETSC_SUCCESS);
1624 if (!matis->A) { /* matrix has not been preallocated yet */
1625 matis->allow_repeated = flg;
1626 PetscFunctionReturn(PETSC_SUCCESS);
1627 }
1628 PetscCheck(!matis->islocalref, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not implemented for local references");
1629 if (matis->allow_repeated) { /* we will assemble the old local matrix if needed */
1630 lA = matis->A;
1631 PetscCall(PetscObjectReference((PetscObject)lA));
1632 }
1633 /* In case flg is True, we only recreate the local matrix */
1634 matis->allow_repeated = flg;
1635 PetscCall(MatSetLocalToGlobalMapping(A, A->rmap->mapping, A->cmap->mapping));
1636 if (lA) { /* assemble previous local matrix if needed */
1637 Mat nA = matis->A;
1638
1639 PetscCall(MatGetLocalToGlobalMapping(nA, &lrmap, &lcmap));
1640 if (!lrmap && !lcmap) {
1641 PetscCall(MatISSetLocalMat(A, lA));
1642 } else {
1643 Mat P = NULL, R = NULL;
1644 MatProductType ptype;
1645
1646 if (lrmap == lcmap) {
1647 ptype = MATPRODUCT_PtAP;
1648 PetscCall(MatCreateFromISLocalToGlobalMapping(lcmap, nA, PETSC_TRUE, PETSC_FALSE, NULL, &P));
1649 PetscCall(MatProductCreate(lA, P, NULL, &nA));
1650 } else {
1651 if (lcmap) PetscCall(MatCreateFromISLocalToGlobalMapping(lcmap, nA, PETSC_TRUE, PETSC_FALSE, NULL, &P));
1652 if (lrmap) PetscCall(MatCreateFromISLocalToGlobalMapping(lrmap, nA, PETSC_FALSE, PETSC_TRUE, NULL, &R));
1653 if (R && P) {
1654 ptype = MATPRODUCT_ABC;
1655 PetscCall(MatProductCreate(R, lA, P, &nA));
1656 } else if (R) {
1657 ptype = MATPRODUCT_AB;
1658 PetscCall(MatProductCreate(R, lA, NULL, &nA));
1659 } else {
1660 ptype = MATPRODUCT_AB;
1661 PetscCall(MatProductCreate(lA, P, NULL, &nA));
1662 }
1663 }
1664 PetscCall(MatProductSetType(nA, ptype));
1665 PetscCall(MatProductSetFromOptions(nA));
1666 PetscCall(MatProductSymbolic(nA));
1667 PetscCall(MatProductNumeric(nA));
1668 PetscCall(MatProductClear(nA));
1669 PetscCall(MatConvert(nA, matis->lmattype, MAT_INPLACE_MATRIX, &nA));
1670 PetscCall(MatISSetLocalMat(A, nA));
1671 PetscCall(MatDestroy(&nA));
1672 PetscCall(MatDestroy(&P));
1673 PetscCall(MatDestroy(&R));
1674 }
1675 }
1676 PetscCall(MatDestroy(&lA));
1677 PetscFunctionReturn(PETSC_SUCCESS);
1678 }
1679
1680 /*@
1681 MatISStoreL2L - Store local-to-local operators during the Galerkin process of computing `MatPtAP()`
1682
1683 Logically Collective
1684
1685 Input Parameters:
1686 + A - the matrix
1687 - store - the boolean flag
1688
1689 Level: advanced
1690
1691 .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateIS()`, `MatISSetPreallocation()`, `MatPtAP()`
1692 @*/
MatISStoreL2L(Mat A,PetscBool store)1693 PetscErrorCode MatISStoreL2L(Mat A, PetscBool store)
1694 {
1695 PetscFunctionBegin;
1696 PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
1697 PetscValidType(A, 1);
1698 PetscValidLogicalCollectiveBool(A, store, 2);
1699 PetscTryMethod(A, "MatISStoreL2L_C", (Mat, PetscBool), (A, store));
1700 PetscFunctionReturn(PETSC_SUCCESS);
1701 }
1702
MatISStoreL2L_IS(Mat A,PetscBool store)1703 static PetscErrorCode MatISStoreL2L_IS(Mat A, PetscBool store)
1704 {
1705 Mat_IS *matis = (Mat_IS *)A->data;
1706
1707 PetscFunctionBegin;
1708 matis->storel2l = store;
1709 if (!store) PetscCall(PetscObjectCompose((PetscObject)A, "_MatIS_PtAP_l2l", NULL));
1710 PetscFunctionReturn(PETSC_SUCCESS);
1711 }
1712
1713 /*@
1714 MatISFixLocalEmpty - Compress out zero local rows from the local matrices
1715
1716 Logically Collective
1717
1718 Input Parameters:
1719 + A - the matrix
1720 - fix - the boolean flag
1721
1722 Level: advanced
1723
1724 Note:
1725 When `fix` is `PETSC_TRUE`, new local matrices and l2g maps are generated during the final assembly process.
1726
1727 .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatCreate()`, `MatCreateIS()`, `MatISSetPreallocation()`, `MatAssemblyEnd()`, `MAT_FINAL_ASSEMBLY`
1728 @*/
MatISFixLocalEmpty(Mat A,PetscBool fix)1729 PetscErrorCode MatISFixLocalEmpty(Mat A, PetscBool fix)
1730 {
1731 PetscFunctionBegin;
1732 PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
1733 PetscValidType(A, 1);
1734 PetscValidLogicalCollectiveBool(A, fix, 2);
1735 PetscTryMethod(A, "MatISFixLocalEmpty_C", (Mat, PetscBool), (A, fix));
1736 PetscFunctionReturn(PETSC_SUCCESS);
1737 }
1738
MatISFixLocalEmpty_IS(Mat A,PetscBool fix)1739 static PetscErrorCode MatISFixLocalEmpty_IS(Mat A, PetscBool fix)
1740 {
1741 Mat_IS *matis = (Mat_IS *)A->data;
1742
1743 PetscFunctionBegin;
1744 matis->locempty = fix;
1745 PetscFunctionReturn(PETSC_SUCCESS);
1746 }
1747
1748 /*@
1749 MatISSetPreallocation - Preallocates memory for a `MATIS` parallel matrix.
1750
1751 Collective
1752
1753 Input Parameters:
1754 + B - the matrix
1755 . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix
1756 (same value is used for all local rows)
1757 . d_nnz - array containing the number of nonzeros in the various rows of the
1758 DIAGONAL portion of the local submatrix (possibly different for each row)
1759 or `NULL`, if `d_nz` is used to specify the nonzero structure.
1760 The size of this array is equal to the number of local rows, i.e `m`.
1761 For matrices that will be factored, you must leave room for (and set)
1762 the diagonal entry even if it is zero.
1763 . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local
1764 submatrix (same value is used for all local rows).
1765 - o_nnz - array containing the number of nonzeros in the various rows of the
1766 OFF-DIAGONAL portion of the local submatrix (possibly different for
1767 each row) or `NULL`, if `o_nz` is used to specify the nonzero
1768 structure. The size of this array is equal to the number
1769 of local rows, i.e `m`.
1770
1771 If the *_nnz parameter is given then the *_nz parameter is ignored
1772
1773 Level: intermediate
1774
1775 Note:
1776 This function has the same interface as the `MATMPIAIJ` preallocation routine in order to simplify the transition
1777 from the asssembled format to the unassembled one. It overestimates the preallocation of `MATIS` local
1778 matrices; for exact preallocation, the user should set the preallocation directly on local matrix objects.
1779
1780 .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateIS()`, `MatMPIAIJSetPreallocation()`, `MatISGetLocalMat()`, `MATIS`
1781 @*/
MatISSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])1782 PetscErrorCode MatISSetPreallocation(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[])
1783 {
1784 PetscFunctionBegin;
1785 PetscValidHeaderSpecific(B, MAT_CLASSID, 1);
1786 PetscValidType(B, 1);
1787 PetscTryMethod(B, "MatISSetPreallocation_C", (Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (B, d_nz, d_nnz, o_nz, o_nnz));
1788 PetscFunctionReturn(PETSC_SUCCESS);
1789 }
1790
MatISSetPreallocation_IS(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])1791 static PetscErrorCode MatISSetPreallocation_IS(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[])
1792 {
1793 Mat_IS *matis = (Mat_IS *)B->data;
1794 PetscInt bs, i, nlocalcols;
1795
1796 PetscFunctionBegin;
1797 PetscCall(MatSetUp(B));
1798 if (!d_nnz)
1799 for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] = d_nz;
1800 else
1801 for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] = d_nnz[i];
1802
1803 if (!o_nnz)
1804 for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] += o_nz;
1805 else
1806 for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] += o_nnz[i];
1807
1808 PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
1809 PetscCall(MatGetSize(matis->A, NULL, &nlocalcols));
1810 PetscCall(MatGetBlockSize(matis->A, &bs));
1811 PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
1812
1813 for (i = 0; i < matis->sf->nleaves; i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i], nlocalcols);
1814 PetscCall(MatSeqAIJSetPreallocation(matis->A, 0, matis->sf_leafdata));
1815 #if defined(PETSC_HAVE_HYPRE)
1816 PetscCall(MatHYPRESetPreallocation(matis->A, 0, matis->sf_leafdata, 0, NULL));
1817 #endif
1818
1819 for (i = 0; i < matis->sf->nleaves / bs; i++) {
1820 PetscInt b;
1821
1822 matis->sf_leafdata[i] = matis->sf_leafdata[i * bs] / bs;
1823 for (b = 1; b < bs; b++) matis->sf_leafdata[i] = PetscMax(matis->sf_leafdata[i], matis->sf_leafdata[i * bs + b] / bs);
1824 }
1825 PetscCall(MatSeqBAIJSetPreallocation(matis->A, bs, 0, matis->sf_leafdata));
1826
1827 nlocalcols /= bs;
1828 for (i = 0; i < matis->sf->nleaves / bs; i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i], nlocalcols - i);
1829 PetscCall(MatSeqSBAIJSetPreallocation(matis->A, bs, 0, matis->sf_leafdata));
1830
1831 /* for other matrix types */
1832 PetscCall(MatSetUp(matis->A));
1833 PetscFunctionReturn(PETSC_SUCCESS);
1834 }
1835
MatConvert_IS_XAIJ(Mat mat,MatType mtype,MatReuse reuse,Mat * M)1836 PETSC_INTERN PetscErrorCode MatConvert_IS_XAIJ(Mat mat, MatType mtype, MatReuse reuse, Mat *M)
1837 {
1838 Mat_IS *matis = (Mat_IS *)mat->data;
1839 Mat local_mat = NULL, MT;
1840 PetscInt rbs, cbs, rows, cols, lrows, lcols;
1841 PetscInt local_rows, local_cols;
1842 PetscBool isseqdense, isseqsbaij, isseqaij, isseqbaij;
1843 PetscMPIInt size;
1844 const PetscScalar *array;
1845
1846 PetscFunctionBegin;
1847 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)mat), &size));
1848 if (size == 1 && mat->rmap->N == matis->A->rmap->N && mat->cmap->N == matis->A->cmap->N && !matis->allow_repeated) {
1849 Mat B;
1850 IS irows = NULL, icols = NULL;
1851 PetscInt rbs, cbs;
1852
1853 PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->rmapping, &rbs));
1854 PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->cmapping, &cbs));
1855 if (reuse != MAT_REUSE_MATRIX) { /* check if l2g maps are one-to-one */
1856 IS rows, cols;
1857 const PetscInt *ridxs, *cidxs;
1858 PetscInt i, nw;
1859 PetscBT work;
1860
1861 PetscCall(ISLocalToGlobalMappingGetBlockIndices(matis->rmapping, &ridxs));
1862 PetscCall(ISLocalToGlobalMappingGetSize(matis->rmapping, &nw));
1863 nw = nw / rbs;
1864 PetscCall(PetscBTCreate(nw, &work));
1865 for (i = 0; i < nw; i++) PetscCall(PetscBTSet(work, ridxs[i]));
1866 for (i = 0; i < nw; i++)
1867 if (!PetscBTLookup(work, i)) break;
1868 if (i == nw) {
1869 PetscCall(ISCreateBlock(PETSC_COMM_SELF, rbs, nw, ridxs, PETSC_USE_POINTER, &rows));
1870 PetscCall(ISSetPermutation(rows));
1871 PetscCall(ISInvertPermutation(rows, PETSC_DECIDE, &irows));
1872 PetscCall(ISDestroy(&rows));
1873 }
1874 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(matis->rmapping, &ridxs));
1875 PetscCall(PetscBTDestroy(&work));
1876 if (irows && matis->rmapping != matis->cmapping) {
1877 PetscCall(ISLocalToGlobalMappingGetBlockIndices(matis->cmapping, &cidxs));
1878 PetscCall(ISLocalToGlobalMappingGetSize(matis->cmapping, &nw));
1879 nw = nw / cbs;
1880 PetscCall(PetscBTCreate(nw, &work));
1881 for (i = 0; i < nw; i++) PetscCall(PetscBTSet(work, cidxs[i]));
1882 for (i = 0; i < nw; i++)
1883 if (!PetscBTLookup(work, i)) break;
1884 if (i == nw) {
1885 PetscCall(ISCreateBlock(PETSC_COMM_SELF, cbs, nw, cidxs, PETSC_USE_POINTER, &cols));
1886 PetscCall(ISSetPermutation(cols));
1887 PetscCall(ISInvertPermutation(cols, PETSC_DECIDE, &icols));
1888 PetscCall(ISDestroy(&cols));
1889 }
1890 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(matis->cmapping, &cidxs));
1891 PetscCall(PetscBTDestroy(&work));
1892 } else if (irows) {
1893 PetscCall(PetscObjectReference((PetscObject)irows));
1894 icols = irows;
1895 }
1896 } else {
1897 PetscCall(PetscObjectQuery((PetscObject)*M, "_MatIS_IS_XAIJ_irows", (PetscObject *)&irows));
1898 PetscCall(PetscObjectQuery((PetscObject)*M, "_MatIS_IS_XAIJ_icols", (PetscObject *)&icols));
1899 if (irows) PetscCall(PetscObjectReference((PetscObject)irows));
1900 if (icols) PetscCall(PetscObjectReference((PetscObject)icols));
1901 }
1902 if (!irows || !icols) {
1903 PetscCall(ISDestroy(&icols));
1904 PetscCall(ISDestroy(&irows));
1905 goto general_assembly;
1906 }
1907 PetscCall(MatConvert(matis->A, mtype, MAT_INITIAL_MATRIX, &B));
1908 if (reuse != MAT_INPLACE_MATRIX) {
1909 PetscCall(MatCreateSubMatrix(B, irows, icols, reuse, M));
1910 PetscCall(PetscObjectCompose((PetscObject)*M, "_MatIS_IS_XAIJ_irows", (PetscObject)irows));
1911 PetscCall(PetscObjectCompose((PetscObject)*M, "_MatIS_IS_XAIJ_icols", (PetscObject)icols));
1912 } else {
1913 Mat C;
1914
1915 PetscCall(MatCreateSubMatrix(B, irows, icols, MAT_INITIAL_MATRIX, &C));
1916 PetscCall(MatHeaderReplace(mat, &C));
1917 }
1918 PetscCall(MatDestroy(&B));
1919 PetscCall(ISDestroy(&icols));
1920 PetscCall(ISDestroy(&irows));
1921 PetscFunctionReturn(PETSC_SUCCESS);
1922 }
1923 general_assembly:
1924 PetscCall(MatGetSize(mat, &rows, &cols));
1925 PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->rmapping, &rbs));
1926 PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->cmapping, &cbs));
1927 PetscCall(MatGetLocalSize(mat, &lrows, &lcols));
1928 PetscCall(MatGetSize(matis->A, &local_rows, &local_cols));
1929 PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQDENSE, &isseqdense));
1930 PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQAIJ, &isseqaij));
1931 PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQBAIJ, &isseqbaij));
1932 PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQSBAIJ, &isseqsbaij));
1933 PetscCheck(isseqdense || isseqaij || isseqbaij || isseqsbaij, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)matis->A)->type_name);
1934 if (PetscDefined(USE_DEBUG)) {
1935 PetscBool lb[4], bb[4];
1936
1937 lb[0] = isseqdense;
1938 lb[1] = isseqaij;
1939 lb[2] = isseqbaij;
1940 lb[3] = isseqsbaij;
1941 PetscCallMPI(MPIU_Allreduce(lb, bb, 4, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)mat)));
1942 PetscCheck(bb[0] || bb[1] || bb[2] || bb[3], PETSC_COMM_SELF, PETSC_ERR_SUP, "Local matrices must have the same type");
1943 }
1944
1945 if (reuse != MAT_REUSE_MATRIX) {
1946 PetscCount ncoo;
1947 PetscInt *coo_i, *coo_j;
1948
1949 PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &MT));
1950 PetscCall(MatSetSizes(MT, lrows, lcols, rows, cols));
1951 PetscCall(MatSetType(MT, mtype));
1952 PetscCall(MatSetBlockSizes(MT, rbs, cbs));
1953 if (!isseqaij && !isseqdense) {
1954 PetscCall(MatConvert(matis->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &local_mat));
1955 } else {
1956 PetscCall(PetscObjectReference((PetscObject)matis->A));
1957 local_mat = matis->A;
1958 }
1959 PetscCall(MatSetLocalToGlobalMapping(MT, matis->rmapping, matis->cmapping));
1960 if (isseqdense) {
1961 PetscInt nr, nc;
1962
1963 PetscCall(MatGetSize(local_mat, &nr, &nc));
1964 ncoo = nr * nc;
1965 PetscCall(PetscMalloc2(ncoo, &coo_i, ncoo, &coo_j));
1966 for (PetscInt j = 0; j < nc; j++) {
1967 for (PetscInt i = 0; i < nr; i++) {
1968 coo_i[j * nr + i] = i;
1969 coo_j[j * nr + i] = j;
1970 }
1971 }
1972 } else {
1973 const PetscInt *ii, *jj;
1974 PetscInt nr;
1975 PetscBool done;
1976
1977 PetscCall(MatGetRowIJ(local_mat, 0, PETSC_FALSE, PETSC_FALSE, &nr, &ii, &jj, &done));
1978 PetscCheck(done, PetscObjectComm((PetscObject)local_mat), PETSC_ERR_PLIB, "Error in MatGetRowIJ");
1979 ncoo = ii[nr];
1980 PetscCall(PetscMalloc2(ncoo, &coo_i, ncoo, &coo_j));
1981 PetscCall(PetscArraycpy(coo_j, jj, ncoo));
1982 for (PetscInt i = 0; i < nr; i++) {
1983 for (PetscInt j = ii[i]; j < ii[i + 1]; j++) coo_i[j] = i;
1984 }
1985 PetscCall(MatRestoreRowIJ(local_mat, 0, PETSC_FALSE, PETSC_FALSE, &nr, &ii, &jj, &done));
1986 PetscCheck(done, PetscObjectComm((PetscObject)local_mat), PETSC_ERR_PLIB, "Error in MatRestoreRowIJ");
1987 }
1988 PetscCall(MatSetPreallocationCOOLocal(MT, ncoo, coo_i, coo_j));
1989 PetscCall(PetscFree2(coo_i, coo_j));
1990 } else {
1991 PetscInt mrbs, mcbs, mrows, mcols, mlrows, mlcols;
1992
1993 /* some checks */
1994 MT = *M;
1995 PetscCall(MatGetBlockSizes(MT, &mrbs, &mcbs));
1996 PetscCall(MatGetSize(MT, &mrows, &mcols));
1997 PetscCall(MatGetLocalSize(MT, &mlrows, &mlcols));
1998 PetscCheck(mrows == rows, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of rows (%" PetscInt_FMT " != %" PetscInt_FMT ")", rows, mrows);
1999 PetscCheck(mcols == cols, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of cols (%" PetscInt_FMT " != %" PetscInt_FMT ")", cols, mcols);
2000 PetscCheck(mlrows == lrows, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of local rows (%" PetscInt_FMT " != %" PetscInt_FMT ")", lrows, mlrows);
2001 PetscCheck(mlcols == lcols, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of local cols (%" PetscInt_FMT " != %" PetscInt_FMT ")", lcols, mlcols);
2002 PetscCheck(mrbs == rbs, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong row block size (%" PetscInt_FMT " != %" PetscInt_FMT ")", rbs, mrbs);
2003 PetscCheck(mcbs == cbs, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong col block size (%" PetscInt_FMT " != %" PetscInt_FMT ")", cbs, mcbs);
2004 PetscCall(MatZeroEntries(MT));
2005 if (!isseqaij && !isseqdense) {
2006 PetscCall(MatConvert(matis->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &local_mat));
2007 } else {
2008 PetscCall(PetscObjectReference((PetscObject)matis->A));
2009 local_mat = matis->A;
2010 }
2011 }
2012
2013 /* Set values */
2014 if (isseqdense) {
2015 PetscCall(MatDenseGetArrayRead(local_mat, &array));
2016 PetscCall(MatSetValuesCOO(MT, array, INSERT_VALUES));
2017 PetscCall(MatDenseRestoreArrayRead(local_mat, &array));
2018 } else {
2019 PetscCall(MatSeqAIJGetArrayRead(local_mat, &array));
2020 PetscCall(MatSetValuesCOO(MT, array, INSERT_VALUES));
2021 PetscCall(MatSeqAIJRestoreArrayRead(local_mat, &array));
2022 }
2023 PetscCall(MatDestroy(&local_mat));
2024 PetscCall(MatAssemblyBegin(MT, MAT_FINAL_ASSEMBLY));
2025 PetscCall(MatAssemblyEnd(MT, MAT_FINAL_ASSEMBLY));
2026 if (reuse == MAT_INPLACE_MATRIX) {
2027 PetscCall(MatHeaderReplace(mat, &MT));
2028 } else if (reuse == MAT_INITIAL_MATRIX) {
2029 *M = MT;
2030 }
2031 PetscFunctionReturn(PETSC_SUCCESS);
2032 }
2033
MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat * newmat)2034 static PetscErrorCode MatDuplicate_IS(Mat mat, MatDuplicateOption op, Mat *newmat)
2035 {
2036 Mat_IS *matis = (Mat_IS *)mat->data;
2037 PetscInt rbs, cbs, m, n, M, N;
2038 Mat B, localmat;
2039
2040 PetscFunctionBegin;
2041 PetscCall(ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping, &rbs));
2042 PetscCall(ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping, &cbs));
2043 PetscCall(MatGetSize(mat, &M, &N));
2044 PetscCall(MatGetLocalSize(mat, &m, &n));
2045 PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &B));
2046 PetscCall(MatSetSizes(B, m, n, M, N));
2047 PetscCall(MatSetBlockSize(B, rbs == cbs ? rbs : 1));
2048 PetscCall(MatSetType(B, MATIS));
2049 PetscCall(MatISSetLocalMatType(B, matis->lmattype));
2050 PetscCall(MatISSetAllowRepeated(B, matis->allow_repeated));
2051 PetscCall(MatSetLocalToGlobalMapping(B, mat->rmap->mapping, mat->cmap->mapping));
2052 PetscCall(MatDuplicate(matis->A, op, &localmat));
2053 PetscCall(MatSetLocalToGlobalMapping(localmat, matis->A->rmap->mapping, matis->A->cmap->mapping));
2054 PetscCall(MatISSetLocalMat(B, localmat));
2055 PetscCall(MatDestroy(&localmat));
2056 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2057 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
2058 *newmat = B;
2059 PetscFunctionReturn(PETSC_SUCCESS);
2060 }
2061
MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool * flg)2062 static PetscErrorCode MatIsHermitian_IS(Mat A, PetscReal tol, PetscBool *flg)
2063 {
2064 Mat_IS *matis = (Mat_IS *)A->data;
2065 PetscBool local_sym;
2066
2067 PetscFunctionBegin;
2068 PetscCall(MatIsHermitian(matis->A, tol, &local_sym));
2069 PetscCallMPI(MPIU_Allreduce(&local_sym, flg, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
2070 PetscFunctionReturn(PETSC_SUCCESS);
2071 }
2072
MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool * flg)2073 static PetscErrorCode MatIsSymmetric_IS(Mat A, PetscReal tol, PetscBool *flg)
2074 {
2075 Mat_IS *matis = (Mat_IS *)A->data;
2076 PetscBool local_sym;
2077
2078 PetscFunctionBegin;
2079 if (matis->rmapping != matis->cmapping) {
2080 *flg = PETSC_FALSE;
2081 PetscFunctionReturn(PETSC_SUCCESS);
2082 }
2083 PetscCall(MatIsSymmetric(matis->A, tol, &local_sym));
2084 PetscCallMPI(MPIU_Allreduce(&local_sym, flg, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
2085 PetscFunctionReturn(PETSC_SUCCESS);
2086 }
2087
MatIsStructurallySymmetric_IS(Mat A,PetscBool * flg)2088 static PetscErrorCode MatIsStructurallySymmetric_IS(Mat A, PetscBool *flg)
2089 {
2090 Mat_IS *matis = (Mat_IS *)A->data;
2091 PetscBool local_sym;
2092
2093 PetscFunctionBegin;
2094 if (matis->rmapping != matis->cmapping) {
2095 *flg = PETSC_FALSE;
2096 PetscFunctionReturn(PETSC_SUCCESS);
2097 }
2098 PetscCall(MatIsStructurallySymmetric(matis->A, &local_sym));
2099 PetscCallMPI(MPIU_Allreduce(&local_sym, flg, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
2100 PetscFunctionReturn(PETSC_SUCCESS);
2101 }
2102
MatDestroy_IS(Mat A)2103 static PetscErrorCode MatDestroy_IS(Mat A)
2104 {
2105 Mat_IS *b = (Mat_IS *)A->data;
2106
2107 PetscFunctionBegin;
2108 PetscCall(PetscFree(b->bdiag));
2109 PetscCall(PetscFree(b->lmattype));
2110 PetscCall(MatDestroy(&b->A));
2111 PetscCall(VecScatterDestroy(&b->cctx));
2112 PetscCall(VecScatterDestroy(&b->rctx));
2113 PetscCall(VecDestroy(&b->x));
2114 PetscCall(VecDestroy(&b->y));
2115 PetscCall(VecDestroy(&b->counter));
2116 PetscCall(ISDestroy(&b->getsub_ris));
2117 PetscCall(ISDestroy(&b->getsub_cis));
2118 if (b->sf != b->csf) {
2119 PetscCall(PetscSFDestroy(&b->csf));
2120 PetscCall(PetscFree2(b->csf_rootdata, b->csf_leafdata));
2121 } else b->csf = NULL;
2122 PetscCall(PetscSFDestroy(&b->sf));
2123 PetscCall(PetscFree2(b->sf_rootdata, b->sf_leafdata));
2124 PetscCall(ISLocalToGlobalMappingDestroy(&b->rmapping));
2125 PetscCall(ISLocalToGlobalMappingDestroy(&b->cmapping));
2126 PetscCall(MatDestroy(&b->dA));
2127 PetscCall(MatDestroy(&b->assembledA));
2128 PetscCall(PetscFree(A->data));
2129 PetscCall(PetscObjectChangeTypeName((PetscObject)A, NULL));
2130 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMatType_C", NULL));
2131 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalMat_C", NULL));
2132 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMat_C", NULL));
2133 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISRestoreLocalMat_C", NULL));
2134 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetPreallocation_C", NULL));
2135 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISStoreL2L_C", NULL));
2136 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISFixLocalEmpty_C", NULL));
2137 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalToGlobalMapping_C", NULL));
2138 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpiaij_C", NULL));
2139 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpibaij_C", NULL));
2140 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpisbaij_C", NULL));
2141 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqaij_C", NULL));
2142 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqbaij_C", NULL));
2143 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqsbaij_C", NULL));
2144 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_aij_C", NULL));
2145 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOOLocal_C", NULL));
2146 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", NULL));
2147 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", NULL));
2148 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetAllowRepeated_C", NULL));
2149 PetscFunctionReturn(PETSC_SUCCESS);
2150 }
2151
MatMult_IS(Mat A,Vec x,Vec y)2152 static PetscErrorCode MatMult_IS(Mat A, Vec x, Vec y)
2153 {
2154 Mat_IS *is = (Mat_IS *)A->data;
2155 PetscScalar zero = 0.0;
2156
2157 PetscFunctionBegin;
2158 /* scatter the global vector x into the local work vector */
2159 PetscCall(VecScatterBegin(is->cctx, x, is->x, INSERT_VALUES, SCATTER_FORWARD));
2160 PetscCall(VecScatterEnd(is->cctx, x, is->x, INSERT_VALUES, SCATTER_FORWARD));
2161
2162 /* multiply the local matrix */
2163 PetscCall(MatMult(is->A, is->x, is->y));
2164
2165 /* scatter product back into global memory */
2166 PetscCall(VecSet(y, zero));
2167 PetscCall(VecScatterBegin(is->rctx, is->y, y, ADD_VALUES, SCATTER_REVERSE));
2168 PetscCall(VecScatterEnd(is->rctx, is->y, y, ADD_VALUES, SCATTER_REVERSE));
2169 PetscFunctionReturn(PETSC_SUCCESS);
2170 }
2171
MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)2172 static PetscErrorCode MatMultAdd_IS(Mat A, Vec v1, Vec v2, Vec v3)
2173 {
2174 Vec temp_vec;
2175
2176 PetscFunctionBegin; /* v3 = v2 + A * v1.*/
2177 if (v3 != v2) {
2178 PetscCall(MatMult(A, v1, v3));
2179 PetscCall(VecAXPY(v3, 1.0, v2));
2180 } else {
2181 PetscCall(VecDuplicate(v2, &temp_vec));
2182 PetscCall(MatMult(A, v1, temp_vec));
2183 PetscCall(VecAXPY(temp_vec, 1.0, v2));
2184 PetscCall(VecCopy(temp_vec, v3));
2185 PetscCall(VecDestroy(&temp_vec));
2186 }
2187 PetscFunctionReturn(PETSC_SUCCESS);
2188 }
2189
MatMultTranspose_IS(Mat A,Vec y,Vec x)2190 static PetscErrorCode MatMultTranspose_IS(Mat A, Vec y, Vec x)
2191 {
2192 Mat_IS *is = (Mat_IS *)A->data;
2193
2194 PetscFunctionBegin;
2195 /* scatter the global vector x into the local work vector */
2196 PetscCall(VecScatterBegin(is->rctx, y, is->y, INSERT_VALUES, SCATTER_FORWARD));
2197 PetscCall(VecScatterEnd(is->rctx, y, is->y, INSERT_VALUES, SCATTER_FORWARD));
2198
2199 /* multiply the local matrix */
2200 PetscCall(MatMultTranspose(is->A, is->y, is->x));
2201
2202 /* scatter product back into global vector */
2203 PetscCall(VecSet(x, 0));
2204 PetscCall(VecScatterBegin(is->cctx, is->x, x, ADD_VALUES, SCATTER_REVERSE));
2205 PetscCall(VecScatterEnd(is->cctx, is->x, x, ADD_VALUES, SCATTER_REVERSE));
2206 PetscFunctionReturn(PETSC_SUCCESS);
2207 }
2208
MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)2209 static PetscErrorCode MatMultTransposeAdd_IS(Mat A, Vec v1, Vec v2, Vec v3)
2210 {
2211 Vec temp_vec;
2212
2213 PetscFunctionBegin; /* v3 = v2 + A' * v1.*/
2214 if (v3 != v2) {
2215 PetscCall(MatMultTranspose(A, v1, v3));
2216 PetscCall(VecAXPY(v3, 1.0, v2));
2217 } else {
2218 PetscCall(VecDuplicate(v2, &temp_vec));
2219 PetscCall(MatMultTranspose(A, v1, temp_vec));
2220 PetscCall(VecAXPY(temp_vec, 1.0, v2));
2221 PetscCall(VecCopy(temp_vec, v3));
2222 PetscCall(VecDestroy(&temp_vec));
2223 }
2224 PetscFunctionReturn(PETSC_SUCCESS);
2225 }
2226
ISLocalToGlobalMappingView_Multi(ISLocalToGlobalMapping mapping,PetscInt lsize,PetscInt gsize,const PetscInt vblocks[],PetscViewer viewer)2227 static PetscErrorCode ISLocalToGlobalMappingView_Multi(ISLocalToGlobalMapping mapping, PetscInt lsize, PetscInt gsize, const PetscInt vblocks[], PetscViewer viewer)
2228 {
2229 PetscInt tr[3], n;
2230 const PetscInt *indices;
2231
2232 PetscFunctionBegin;
2233 tr[0] = IS_LTOGM_FILE_CLASSID;
2234 tr[1] = 1;
2235 tr[2] = gsize;
2236 PetscCall(PetscViewerBinaryWrite(viewer, tr, 3, PETSC_INT));
2237 PetscCall(PetscViewerBinaryWriteAll(viewer, vblocks, lsize, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
2238 PetscCall(ISLocalToGlobalMappingGetSize(mapping, &n));
2239 PetscCall(ISLocalToGlobalMappingGetIndices(mapping, &indices));
2240 PetscCall(PetscViewerBinaryWriteAll(viewer, indices, n, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
2241 PetscCall(ISLocalToGlobalMappingRestoreIndices(mapping, &indices));
2242 PetscFunctionReturn(PETSC_SUCCESS);
2243 }
2244
MatView_IS(Mat A,PetscViewer viewer)2245 static PetscErrorCode MatView_IS(Mat A, PetscViewer viewer)
2246 {
2247 Mat_IS *a = (Mat_IS *)A->data;
2248 PetscViewer sviewer;
2249 PetscBool isascii, isbinary, viewl2g = PETSC_FALSE, native;
2250 PetscViewerFormat format;
2251 ISLocalToGlobalMapping rmap, cmap;
2252
2253 PetscFunctionBegin;
2254 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2255 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
2256 PetscCall(PetscViewerGetFormat(viewer, &format));
2257 native = (PetscBool)(format == PETSC_VIEWER_NATIVE);
2258 if (native) {
2259 rmap = A->rmap->mapping;
2260 cmap = A->cmap->mapping;
2261 } else {
2262 rmap = a->rmapping;
2263 cmap = a->cmapping;
2264 }
2265 if (isascii) {
2266 if (format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS);
2267 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL || format == PETSC_VIEWER_ASCII_MATLAB) viewl2g = PETSC_TRUE;
2268 } else if (isbinary) {
2269 PetscInt tr[6], nr, nc, lsize = 0;
2270 char lmattype[64] = {'\0'};
2271 PetscMPIInt size;
2272 PetscBool skipHeader, vbs = PETSC_FALSE;
2273 IS is;
2274 const PetscInt *vblocks = NULL;
2275
2276 PetscCall(PetscViewerSetUp(viewer));
2277 PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)A)->prefix, "-mat_is_view_variableblocksizes", &vbs, NULL));
2278 if (vbs) {
2279 PetscCall(MatGetVariableBlockSizes(a->A, &lsize, &vblocks));
2280 PetscCall(PetscMPIIntCast(lsize, &size));
2281 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &size, 1, MPI_INT, MPI_SUM, PetscObjectComm((PetscObject)viewer)));
2282 } else {
2283 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)viewer), &size));
2284 }
2285 tr[0] = MAT_FILE_CLASSID;
2286 tr[1] = A->rmap->N;
2287 tr[2] = A->cmap->N;
2288 tr[3] = -size; /* AIJ stores nnz here */
2289 tr[4] = (PetscInt)(rmap == cmap);
2290 tr[5] = a->allow_repeated;
2291 PetscCall(PetscSNPrintf(lmattype, sizeof(lmattype), "%s", a->lmattype));
2292
2293 PetscCall(PetscViewerBinaryWrite(viewer, tr, PETSC_STATIC_ARRAY_LENGTH(tr), PETSC_INT));
2294 PetscCall(PetscViewerBinaryWrite(viewer, lmattype, sizeof(lmattype), PETSC_CHAR));
2295
2296 /* first dump l2g info (we need the header for proper loading on different number of processes) */
2297 PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
2298 PetscCall(PetscViewerBinarySetSkipHeader(viewer, PETSC_FALSE));
2299 if (vbs) {
2300 PetscCall(ISLocalToGlobalMappingView_Multi(rmap, lsize, size, vblocks, viewer));
2301 if (cmap != rmap) PetscCall(ISLocalToGlobalMappingView_Multi(cmap, lsize, size, vblocks, viewer));
2302 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)viewer), lsize, vblocks, PETSC_USE_POINTER, &is));
2303 PetscCall(ISView(is, viewer));
2304 PetscCall(ISView(is, viewer));
2305 PetscCall(ISDestroy(&is));
2306 } else {
2307 PetscCall(ISLocalToGlobalMappingView(rmap, viewer));
2308 if (cmap != rmap) PetscCall(ISLocalToGlobalMappingView(cmap, viewer));
2309
2310 /* then the sizes of the local matrices */
2311 PetscCall(MatGetSize(a->A, &nr, &nc));
2312 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)viewer), 1, &nr, PETSC_USE_POINTER, &is));
2313 PetscCall(ISView(is, viewer));
2314 PetscCall(ISDestroy(&is));
2315 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)viewer), 1, &nc, PETSC_USE_POINTER, &is));
2316 PetscCall(ISView(is, viewer));
2317 PetscCall(ISDestroy(&is));
2318 }
2319 PetscCall(PetscViewerBinarySetSkipHeader(viewer, skipHeader));
2320 }
2321 if (format == PETSC_VIEWER_ASCII_MATLAB) {
2322 char name[64];
2323 PetscMPIInt size, rank;
2324
2325 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)viewer), &size));
2326 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
2327 if (size > 1) PetscCall(PetscSNPrintf(name, sizeof(name), "lmat_%d", rank));
2328 else PetscCall(PetscSNPrintf(name, sizeof(name), "lmat"));
2329 PetscCall(PetscObjectSetName((PetscObject)a->A, name));
2330 }
2331
2332 /* Dump the local matrices */
2333 if (isbinary) { /* ViewerGetSubViewer does not work in parallel */
2334 PetscBool isaij;
2335 PetscInt nr, nc;
2336 Mat lA, B;
2337 Mat_MPIAIJ *b;
2338
2339 /* We create a temporary MPIAIJ matrix that stores the unassembled operator */
2340 PetscCall(PetscObjectBaseTypeCompare((PetscObject)a->A, MATAIJ, &isaij));
2341 if (!isaij) PetscCall(MatConvert(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &lA));
2342 else {
2343 PetscCall(PetscObjectReference((PetscObject)a->A));
2344 lA = a->A;
2345 }
2346 PetscCall(MatCreate(PetscObjectComm((PetscObject)viewer), &B));
2347 PetscCall(MatSetType(B, MATMPIAIJ));
2348 PetscCall(MatGetSize(lA, &nr, &nc));
2349 PetscCall(MatSetSizes(B, nr, nc, PETSC_DECIDE, PETSC_DECIDE));
2350 PetscCall(MatMPIAIJSetPreallocation(B, 0, NULL, 0, NULL));
2351
2352 b = (Mat_MPIAIJ *)B->data;
2353 PetscCall(MatDestroy(&b->A));
2354 b->A = lA;
2355
2356 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2357 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
2358 PetscCall(MatView(B, viewer));
2359 PetscCall(MatDestroy(&B));
2360 } else {
2361 PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
2362 PetscCall(MatView(a->A, sviewer));
2363 PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
2364 }
2365
2366 /* with ASCII, we dump the l2gmaps at the end */
2367 if (viewl2g) {
2368 if (format == PETSC_VIEWER_ASCII_MATLAB) {
2369 PetscCall(PetscObjectSetName((PetscObject)rmap, "row"));
2370 PetscCall(ISLocalToGlobalMappingView(rmap, viewer));
2371 PetscCall(PetscObjectSetName((PetscObject)cmap, "col"));
2372 PetscCall(ISLocalToGlobalMappingView(cmap, viewer));
2373 } else {
2374 PetscCall(ISLocalToGlobalMappingView(rmap, viewer));
2375 PetscCall(ISLocalToGlobalMappingView(cmap, viewer));
2376 }
2377 }
2378 PetscFunctionReturn(PETSC_SUCCESS);
2379 }
2380
ISLocalToGlobalMappingHasRepeatedLocal_Private(ISLocalToGlobalMapping map,PetscBool * has)2381 static PetscErrorCode ISLocalToGlobalMappingHasRepeatedLocal_Private(ISLocalToGlobalMapping map, PetscBool *has)
2382 {
2383 const PetscInt *idxs;
2384 PetscHSetI ht;
2385 PetscInt n, bs;
2386
2387 PetscFunctionBegin;
2388 PetscCall(ISLocalToGlobalMappingGetSize(map, &n));
2389 PetscCall(ISLocalToGlobalMappingGetBlockSize(map, &bs));
2390 PetscCall(ISLocalToGlobalMappingGetBlockIndices(map, &idxs));
2391 PetscCall(PetscHSetICreate(&ht));
2392 *has = PETSC_FALSE;
2393 for (PetscInt i = 0; i < n / bs; i++) {
2394 PetscBool missing = PETSC_TRUE;
2395 if (idxs[i] < 0) continue;
2396 PetscCall(PetscHSetIQueryAdd(ht, idxs[i], &missing));
2397 if (!missing) {
2398 *has = PETSC_TRUE;
2399 break;
2400 }
2401 }
2402 PetscCall(PetscHSetIDestroy(&ht));
2403 PetscFunctionReturn(PETSC_SUCCESS);
2404 }
2405
MatLoad_IS(Mat A,PetscViewer viewer)2406 static PetscErrorCode MatLoad_IS(Mat A, PetscViewer viewer)
2407 {
2408 ISLocalToGlobalMapping rmap, cmap;
2409 MPI_Comm comm = PetscObjectComm((PetscObject)A);
2410 PetscBool isbinary, samel, allow, isbaij;
2411 PetscInt tr[6], M, N, nr, nc, Asize, isn;
2412 const PetscInt *idx;
2413 PetscMPIInt size;
2414 char lmattype[64];
2415 Mat dA, lA;
2416 IS is;
2417
2418 PetscFunctionBegin;
2419 PetscCheckSameComm(A, 1, viewer, 2);
2420 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
2421 PetscCheck(isbinary, PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Invalid viewer of type %s", ((PetscObject)viewer)->type_name);
2422
2423 PetscCall(PetscViewerBinaryRead(viewer, tr, PETSC_STATIC_ARRAY_LENGTH(tr), NULL, PETSC_INT));
2424 PetscCheck(tr[0] == MAT_FILE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a matrix next in file");
2425 PetscCheck(tr[1] >= 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file");
2426 PetscCheck(tr[2] >= 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file");
2427 PetscCheck(tr[3] < 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file");
2428 PetscCheck(tr[4] == 0 || tr[4] == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file");
2429 PetscCheck(tr[5] == 0 || tr[5] == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file");
2430 M = tr[1];
2431 N = tr[2];
2432 Asize = -tr[3];
2433 samel = (PetscBool)tr[4];
2434 allow = (PetscBool)tr[5];
2435 PetscCall(PetscViewerBinaryRead(viewer, lmattype, sizeof(lmattype), NULL, PETSC_CHAR));
2436
2437 /* if we are loading from a larger set of processes, allow repeated entries */
2438 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)viewer), &size));
2439 if (Asize > size) allow = PETSC_TRUE;
2440
2441 /* set global sizes if not set already */
2442 if (A->rmap->N < 0) A->rmap->N = M;
2443 if (A->cmap->N < 0) A->cmap->N = N;
2444 PetscCall(PetscLayoutSetUp(A->rmap));
2445 PetscCall(PetscLayoutSetUp(A->cmap));
2446 PetscCheck(M == A->rmap->N, comm, PETSC_ERR_ARG_SIZ, "Matrix rows should be %" PetscInt_FMT ", found %" PetscInt_FMT, M, A->rmap->N);
2447 PetscCheck(N == A->cmap->N, comm, PETSC_ERR_ARG_SIZ, "Matrix columns should be %" PetscInt_FMT ", found %" PetscInt_FMT, N, A->cmap->N);
2448
2449 /* load l2g maps */
2450 PetscCall(ISLocalToGlobalMappingCreate(comm, 0, 0, NULL, PETSC_USE_POINTER, &rmap));
2451 PetscCall(ISLocalToGlobalMappingLoad(rmap, viewer));
2452 if (!samel) {
2453 PetscCall(ISLocalToGlobalMappingCreate(comm, 0, 0, NULL, PETSC_USE_POINTER, &cmap));
2454 PetscCall(ISLocalToGlobalMappingLoad(cmap, viewer));
2455 } else {
2456 PetscCall(PetscObjectReference((PetscObject)rmap));
2457 cmap = rmap;
2458 }
2459
2460 /* load sizes of local matrices */
2461 PetscCall(ISCreate(comm, &is));
2462 PetscCall(ISSetType(is, ISGENERAL));
2463 PetscCall(ISLoad(is, viewer));
2464 PetscCall(ISGetLocalSize(is, &isn));
2465 PetscCall(ISGetIndices(is, &idx));
2466 nr = 0;
2467 for (PetscInt i = 0; i < isn; i++) nr += idx[i];
2468 PetscCall(ISRestoreIndices(is, &idx));
2469 PetscCall(ISDestroy(&is));
2470 PetscCall(ISCreate(comm, &is));
2471 PetscCall(ISSetType(is, ISGENERAL));
2472 PetscCall(ISLoad(is, viewer));
2473 PetscCall(ISGetLocalSize(is, &isn));
2474 PetscCall(ISGetIndices(is, &idx));
2475 nc = 0;
2476 for (PetscInt i = 0; i < isn; i++) nc += idx[i];
2477 PetscCall(ISRestoreIndices(is, &idx));
2478 PetscCall(ISDestroy(&is));
2479
2480 /* now load the unassembled operator */
2481 PetscCall(MatCreate(comm, &dA));
2482 PetscCall(MatSetType(dA, MATMPIAIJ));
2483 PetscCall(MatSetSizes(dA, nr, nc, PETSC_DECIDE, PETSC_DECIDE));
2484 PetscCall(MatLoad(dA, viewer));
2485 PetscCall(MatMPIAIJGetSeqAIJ(dA, &lA, NULL, NULL));
2486 PetscCall(PetscObjectReference((PetscObject)lA));
2487 PetscCall(MatDestroy(&dA));
2488
2489 /* and convert to the desired format */
2490 PetscCall(PetscStrcmpAny(lmattype, &isbaij, MATSBAIJ, MATSEQSBAIJ, ""));
2491 if (isbaij) PetscCall(MatSetOption(lA, MAT_SYMMETRIC, PETSC_TRUE));
2492 PetscCall(MatConvert(lA, lmattype, MAT_INPLACE_MATRIX, &lA));
2493
2494 /* check if we actually have repeated entries */
2495 if (allow) {
2496 PetscBool rhas, chas, hasrepeated;
2497
2498 PetscCall(ISLocalToGlobalMappingHasRepeatedLocal_Private(rmap, &rhas));
2499 if (rmap != cmap) PetscCall(ISLocalToGlobalMappingHasRepeatedLocal_Private(cmap, &chas));
2500 else chas = rhas;
2501 hasrepeated = (PetscBool)(rhas || chas);
2502 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &hasrepeated, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A)));
2503 if (!hasrepeated) allow = PETSC_FALSE;
2504 }
2505
2506 /* assemble the MATIS object */
2507 PetscCall(MatISSetAllowRepeated(A, allow));
2508 PetscCall(MatSetLocalToGlobalMapping(A, rmap, cmap));
2509 PetscCall(MatISSetLocalMat(A, lA));
2510 PetscCall(MatDestroy(&lA));
2511 PetscCall(ISLocalToGlobalMappingDestroy(&rmap));
2512 PetscCall(ISLocalToGlobalMappingDestroy(&cmap));
2513 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
2514 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
2515 PetscFunctionReturn(PETSC_SUCCESS);
2516 }
2517
MatInvertBlockDiagonal_IS(Mat mat,const PetscScalar ** values)2518 static PetscErrorCode MatInvertBlockDiagonal_IS(Mat mat, const PetscScalar **values)
2519 {
2520 Mat_IS *is = (Mat_IS *)mat->data;
2521 MPI_Datatype nodeType;
2522 const PetscScalar *lv;
2523 PetscInt bs;
2524 PetscMPIInt mbs;
2525
2526 PetscFunctionBegin;
2527 PetscCall(MatGetBlockSize(mat, &bs));
2528 PetscCall(MatSetBlockSize(is->A, bs));
2529 PetscCall(MatInvertBlockDiagonal(is->A, &lv));
2530 if (!is->bdiag) PetscCall(PetscMalloc1(bs * mat->rmap->n, &is->bdiag));
2531 PetscCall(PetscMPIIntCast(bs, &mbs));
2532 PetscCallMPI(MPI_Type_contiguous(mbs, MPIU_SCALAR, &nodeType));
2533 PetscCallMPI(MPI_Type_commit(&nodeType));
2534 PetscCall(PetscSFReduceBegin(is->sf, nodeType, lv, is->bdiag, MPI_REPLACE));
2535 PetscCall(PetscSFReduceEnd(is->sf, nodeType, lv, is->bdiag, MPI_REPLACE));
2536 PetscCallMPI(MPI_Type_free(&nodeType));
2537 if (values) *values = is->bdiag;
2538 PetscFunctionReturn(PETSC_SUCCESS);
2539 }
2540
MatISSetUpScatters_Private(Mat A)2541 static PetscErrorCode MatISSetUpScatters_Private(Mat A)
2542 {
2543 Vec cglobal, rglobal;
2544 IS from;
2545 Mat_IS *is = (Mat_IS *)A->data;
2546 PetscScalar sum;
2547 const PetscInt *garray;
2548 PetscInt nr, rbs, nc, cbs;
2549 VecType rtype;
2550
2551 PetscFunctionBegin;
2552 PetscCall(ISLocalToGlobalMappingGetSize(is->rmapping, &nr));
2553 PetscCall(ISLocalToGlobalMappingGetBlockSize(is->rmapping, &rbs));
2554 PetscCall(ISLocalToGlobalMappingGetSize(is->cmapping, &nc));
2555 PetscCall(ISLocalToGlobalMappingGetBlockSize(is->cmapping, &cbs));
2556 PetscCall(VecDestroy(&is->x));
2557 PetscCall(VecDestroy(&is->y));
2558 PetscCall(VecDestroy(&is->counter));
2559 PetscCall(VecScatterDestroy(&is->rctx));
2560 PetscCall(VecScatterDestroy(&is->cctx));
2561 PetscCall(MatCreateVecs(is->A, &is->x, &is->y));
2562 PetscCall(VecBindToCPU(is->y, PETSC_TRUE));
2563 PetscCall(VecGetRootType_Private(is->y, &rtype));
2564 PetscCall(PetscFree(A->defaultvectype));
2565 PetscCall(PetscStrallocpy(rtype, &A->defaultvectype));
2566 PetscCall(MatCreateVecs(A, &cglobal, &rglobal));
2567 PetscCall(VecBindToCPU(rglobal, PETSC_TRUE));
2568 PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->rmapping, &garray));
2569 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), rbs, nr / rbs, garray, PETSC_USE_POINTER, &from));
2570 PetscCall(VecScatterCreate(rglobal, from, is->y, NULL, &is->rctx));
2571 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping, &garray));
2572 PetscCall(ISDestroy(&from));
2573 if (is->rmapping != is->cmapping) {
2574 PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->cmapping, &garray));
2575 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), cbs, nc / cbs, garray, PETSC_USE_POINTER, &from));
2576 PetscCall(VecScatterCreate(cglobal, from, is->x, NULL, &is->cctx));
2577 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping, &garray));
2578 PetscCall(ISDestroy(&from));
2579 } else {
2580 PetscCall(PetscObjectReference((PetscObject)is->rctx));
2581 is->cctx = is->rctx;
2582 }
2583 PetscCall(VecDestroy(&cglobal));
2584
2585 /* interface counter vector (local) */
2586 PetscCall(VecDuplicate(is->y, &is->counter));
2587 PetscCall(VecBindToCPU(is->counter, PETSC_TRUE));
2588 PetscCall(VecSet(is->y, 1.));
2589 PetscCall(VecScatterBegin(is->rctx, is->y, rglobal, ADD_VALUES, SCATTER_REVERSE));
2590 PetscCall(VecScatterEnd(is->rctx, is->y, rglobal, ADD_VALUES, SCATTER_REVERSE));
2591 PetscCall(VecScatterBegin(is->rctx, rglobal, is->counter, INSERT_VALUES, SCATTER_FORWARD));
2592 PetscCall(VecScatterEnd(is->rctx, rglobal, is->counter, INSERT_VALUES, SCATTER_FORWARD));
2593 PetscCall(VecBindToCPU(is->y, PETSC_FALSE));
2594 PetscCall(VecBindToCPU(is->counter, PETSC_FALSE));
2595
2596 /* special functions for block-diagonal matrices */
2597 PetscCall(VecSum(rglobal, &sum));
2598 A->ops->invertblockdiagonal = NULL;
2599 if ((PetscInt)(PetscRealPart(sum)) == A->rmap->N && A->rmap->N == A->cmap->N && is->rmapping == is->cmapping) A->ops->invertblockdiagonal = MatInvertBlockDiagonal_IS;
2600 PetscCall(VecDestroy(&rglobal));
2601
2602 /* setup SF for general purpose shared indices based communications */
2603 PetscCall(MatISSetUpSF_IS(A));
2604 PetscFunctionReturn(PETSC_SUCCESS);
2605 }
2606
MatISFilterL2GMap(Mat A,ISLocalToGlobalMapping map,ISLocalToGlobalMapping * nmap,ISLocalToGlobalMapping * lmap)2607 static PetscErrorCode MatISFilterL2GMap(Mat A, ISLocalToGlobalMapping map, ISLocalToGlobalMapping *nmap, ISLocalToGlobalMapping *lmap)
2608 {
2609 Mat_IS *matis = (Mat_IS *)A->data;
2610 IS is;
2611 ISLocalToGlobalMappingType l2gtype;
2612 const PetscInt *idxs;
2613 PetscHSetI ht;
2614 PetscInt *nidxs;
2615 PetscInt i, n, bs, c;
2616 PetscBool flg[] = {PETSC_FALSE, PETSC_FALSE};
2617
2618 PetscFunctionBegin;
2619 PetscCall(ISLocalToGlobalMappingGetSize(map, &n));
2620 PetscCall(ISLocalToGlobalMappingGetBlockSize(map, &bs));
2621 PetscCall(ISLocalToGlobalMappingGetBlockIndices(map, &idxs));
2622 PetscCall(PetscHSetICreate(&ht));
2623 PetscCall(PetscMalloc1(n / bs, &nidxs));
2624 for (i = 0, c = 0; i < n / bs; i++) {
2625 PetscBool missing = PETSC_TRUE;
2626 if (idxs[i] < 0) {
2627 flg[0] = PETSC_TRUE;
2628 continue;
2629 }
2630 if (!matis->allow_repeated) PetscCall(PetscHSetIQueryAdd(ht, idxs[i], &missing));
2631 if (!missing) flg[1] = PETSC_TRUE;
2632 else nidxs[c++] = idxs[i];
2633 }
2634 PetscCall(PetscHSetIDestroy(&ht));
2635 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, flg, 2, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A)));
2636 if (!flg[0] && !flg[1]) { /* Entries are all non negative and unique */
2637 *nmap = NULL;
2638 *lmap = NULL;
2639 PetscCall(PetscFree(nidxs));
2640 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(map, &idxs));
2641 PetscFunctionReturn(PETSC_SUCCESS);
2642 }
2643
2644 /* New l2g map without negative indices (and repeated indices if not allowed) */
2645 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), bs, c, nidxs, PETSC_USE_POINTER, &is));
2646 PetscCall(ISLocalToGlobalMappingCreateIS(is, nmap));
2647 PetscCall(ISDestroy(&is));
2648 PetscCall(ISLocalToGlobalMappingGetType(map, &l2gtype));
2649 PetscCall(ISLocalToGlobalMappingSetType(*nmap, l2gtype));
2650
2651 /* New local l2g map for repeated indices if not allowed */
2652 PetscCall(ISGlobalToLocalMappingApplyBlock(*nmap, IS_GTOLM_MASK, n / bs, idxs, NULL, nidxs));
2653 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, n / bs, nidxs, PETSC_USE_POINTER, &is));
2654 PetscCall(ISLocalToGlobalMappingCreateIS(is, lmap));
2655 PetscCall(ISDestroy(&is));
2656 PetscCall(PetscFree(nidxs));
2657 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(map, &idxs));
2658 PetscFunctionReturn(PETSC_SUCCESS);
2659 }
2660
MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)2661 static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A, ISLocalToGlobalMapping rmapping, ISLocalToGlobalMapping cmapping)
2662 {
2663 Mat_IS *is = (Mat_IS *)A->data;
2664 ISLocalToGlobalMapping localrmapping = NULL, localcmapping = NULL;
2665 PetscInt nr, rbs, nc, cbs;
2666 PetscBool cong, freem[] = {PETSC_FALSE, PETSC_FALSE};
2667
2668 PetscFunctionBegin;
2669 if (rmapping) PetscCheckSameComm(A, 1, rmapping, 2);
2670 if (cmapping) PetscCheckSameComm(A, 1, cmapping, 3);
2671
2672 PetscCall(ISLocalToGlobalMappingDestroy(&is->rmapping));
2673 PetscCall(ISLocalToGlobalMappingDestroy(&is->cmapping));
2674 PetscCall(PetscLayoutSetUp(A->rmap));
2675 PetscCall(PetscLayoutSetUp(A->cmap));
2676 PetscCall(MatHasCongruentLayouts(A, &cong));
2677
2678 /* If NULL, local space matches global space */
2679 if (!rmapping) {
2680 IS is;
2681
2682 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->rmap->N, 0, 1, &is));
2683 PetscCall(ISLocalToGlobalMappingCreateIS(is, &rmapping));
2684 PetscCall(ISLocalToGlobalMappingSetBlockSize(rmapping, A->rmap->bs));
2685 PetscCall(ISDestroy(&is));
2686 freem[0] = PETSC_TRUE;
2687 if (!cmapping && cong && A->rmap->bs == A->cmap->bs) cmapping = rmapping;
2688 } else if (!is->islocalref) { /* check if the l2g map has negative or repeated entries */
2689 PetscCall(MatISFilterL2GMap(A, rmapping, &is->rmapping, &localrmapping));
2690 if (rmapping == cmapping) {
2691 PetscCall(PetscObjectReference((PetscObject)is->rmapping));
2692 is->cmapping = is->rmapping;
2693 PetscCall(PetscObjectReference((PetscObject)localrmapping));
2694 localcmapping = localrmapping;
2695 }
2696 }
2697 if (!cmapping) {
2698 IS is;
2699
2700 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->cmap->N, 0, 1, &is));
2701 PetscCall(ISLocalToGlobalMappingCreateIS(is, &cmapping));
2702 PetscCall(ISLocalToGlobalMappingSetBlockSize(cmapping, A->cmap->bs));
2703 PetscCall(ISDestroy(&is));
2704 freem[1] = PETSC_TRUE;
2705 } else if (cmapping != rmapping && !is->islocalref) { /* check if the l2g map has negative or repeated entries */
2706 PetscCall(MatISFilterL2GMap(A, cmapping, &is->cmapping, &localcmapping));
2707 }
2708 if (!is->rmapping) {
2709 PetscCall(PetscObjectReference((PetscObject)rmapping));
2710 is->rmapping = rmapping;
2711 }
2712 if (!is->cmapping) {
2713 PetscCall(PetscObjectReference((PetscObject)cmapping));
2714 is->cmapping = cmapping;
2715 }
2716
2717 /* Clean up */
2718 is->lnnzstate = 0;
2719 PetscCall(MatDestroy(&is->dA));
2720 PetscCall(MatDestroy(&is->assembledA));
2721 PetscCall(MatDestroy(&is->A));
2722 if (is->csf != is->sf) {
2723 PetscCall(PetscSFDestroy(&is->csf));
2724 PetscCall(PetscFree2(is->csf_rootdata, is->csf_leafdata));
2725 } else is->csf = NULL;
2726 PetscCall(PetscSFDestroy(&is->sf));
2727 PetscCall(PetscFree2(is->sf_rootdata, is->sf_leafdata));
2728 PetscCall(PetscFree(is->bdiag));
2729
2730 /* check if the two mappings are actually the same for square matrices since MATIS has some optimization for this case
2731 (DOLFIN passes 2 different objects) */
2732 PetscCall(ISLocalToGlobalMappingGetSize(is->rmapping, &nr));
2733 PetscCall(ISLocalToGlobalMappingGetBlockSize(is->rmapping, &rbs));
2734 PetscCall(ISLocalToGlobalMappingGetSize(is->cmapping, &nc));
2735 PetscCall(ISLocalToGlobalMappingGetBlockSize(is->cmapping, &cbs));
2736 if (is->rmapping != is->cmapping && cong) {
2737 PetscBool same = PETSC_FALSE;
2738 if (nr == nc && cbs == rbs) {
2739 const PetscInt *idxs1, *idxs2;
2740
2741 PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->rmapping, &idxs1));
2742 PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->cmapping, &idxs2));
2743 PetscCall(PetscArraycmp(idxs1, idxs2, nr / rbs, &same));
2744 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping, &idxs1));
2745 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping, &idxs2));
2746 }
2747 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &same, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
2748 if (same) {
2749 PetscCall(ISLocalToGlobalMappingDestroy(&is->cmapping));
2750 PetscCall(PetscObjectReference((PetscObject)is->rmapping));
2751 is->cmapping = is->rmapping;
2752 }
2753 }
2754 PetscCall(PetscLayoutSetBlockSize(A->rmap, rbs));
2755 PetscCall(PetscLayoutSetBlockSize(A->cmap, cbs));
2756 /* Pass the user defined maps to the layout */
2757 PetscCall(PetscLayoutSetISLocalToGlobalMapping(A->rmap, rmapping));
2758 PetscCall(PetscLayoutSetISLocalToGlobalMapping(A->cmap, cmapping));
2759 if (freem[0]) PetscCall(ISLocalToGlobalMappingDestroy(&rmapping));
2760 if (freem[1]) PetscCall(ISLocalToGlobalMappingDestroy(&cmapping));
2761
2762 if (!is->islocalref) {
2763 /* Create the local matrix A */
2764 PetscCall(MatCreate(PETSC_COMM_SELF, &is->A));
2765 PetscCall(MatSetType(is->A, is->lmattype));
2766 PetscCall(MatSetSizes(is->A, nr, nc, nr, nc));
2767 PetscCall(MatSetBlockSizes(is->A, rbs, cbs));
2768 PetscCall(MatSetOptionsPrefix(is->A, "is_"));
2769 PetscCall(MatAppendOptionsPrefix(is->A, ((PetscObject)A)->prefix));
2770 PetscCall(PetscLayoutSetUp(is->A->rmap));
2771 PetscCall(PetscLayoutSetUp(is->A->cmap));
2772 PetscCall(MatSetLocalToGlobalMapping(is->A, localrmapping, localcmapping));
2773 PetscCall(ISLocalToGlobalMappingDestroy(&localrmapping));
2774 PetscCall(ISLocalToGlobalMappingDestroy(&localcmapping));
2775 /* setup scatters and local vectors for MatMult */
2776 PetscCall(MatISSetUpScatters_Private(A));
2777 }
2778 A->preallocated = PETSC_TRUE;
2779 PetscFunctionReturn(PETSC_SUCCESS);
2780 }
2781
MatSetUp_IS(Mat A)2782 static PetscErrorCode MatSetUp_IS(Mat A)
2783 {
2784 Mat_IS *is = (Mat_IS *)A->data;
2785 ISLocalToGlobalMapping rmap, cmap;
2786
2787 PetscFunctionBegin;
2788 if (!is->sf) {
2789 PetscCall(MatGetLocalToGlobalMapping(A, &rmap, &cmap));
2790 PetscCall(MatSetLocalToGlobalMapping(A, rmap, cmap));
2791 }
2792 PetscFunctionReturn(PETSC_SUCCESS);
2793 }
2794
MatSetValues_IS(Mat mat,PetscInt m,const PetscInt * rows,PetscInt n,const PetscInt * cols,const PetscScalar * values,InsertMode addv)2795 static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2796 {
2797 Mat_IS *is = (Mat_IS *)mat->data;
2798 PetscInt buf[2 * MATIS_MAX_ENTRIES_INSERTION], *rows_l = NULL, *cols_l = NULL;
2799
2800 PetscFunctionBegin;
2801 IndexSpaceGet(buf, m, n, rows_l, cols_l);
2802 PetscCall(ISGlobalToLocalMappingApply(is->rmapping, IS_GTOLM_MASK, m, rows, &m, rows_l));
2803 if (m != n || rows != cols || is->cmapping != is->rmapping) {
2804 PetscCall(ISGlobalToLocalMappingApply(is->cmapping, IS_GTOLM_MASK, n, cols, &n, cols_l));
2805 PetscCall(MatSetValues(is->A, m, rows_l, n, cols_l, values, addv));
2806 } else {
2807 PetscCall(MatSetValues(is->A, m, rows_l, m, rows_l, values, addv));
2808 }
2809 IndexSpaceRestore(buf, m, n, rows_l, cols_l);
2810 PetscFunctionReturn(PETSC_SUCCESS);
2811 }
2812
MatSetValuesBlocked_IS(Mat mat,PetscInt m,const PetscInt * rows,PetscInt n,const PetscInt * cols,const PetscScalar * values,InsertMode addv)2813 static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2814 {
2815 Mat_IS *is = (Mat_IS *)mat->data;
2816 PetscInt buf[2 * MATIS_MAX_ENTRIES_INSERTION], *rows_l = NULL, *cols_l = NULL;
2817
2818 PetscFunctionBegin;
2819 IndexSpaceGet(buf, m, n, rows_l, cols_l);
2820 PetscCall(ISGlobalToLocalMappingApplyBlock(is->rmapping, IS_GTOLM_MASK, m, rows, &m, rows_l));
2821 if (m != n || rows != cols || is->cmapping != is->rmapping) {
2822 PetscCall(ISGlobalToLocalMappingApplyBlock(is->cmapping, IS_GTOLM_MASK, n, cols, &n, cols_l));
2823 PetscCall(MatSetValuesBlocked(is->A, m, rows_l, n, cols_l, values, addv));
2824 } else {
2825 PetscCall(MatSetValuesBlocked(is->A, m, rows_l, m, rows_l, values, addv));
2826 }
2827 IndexSpaceRestore(buf, m, n, rows_l, cols_l);
2828 PetscFunctionReturn(PETSC_SUCCESS);
2829 }
2830
MatSetValuesLocal_IS(Mat A,PetscInt m,const PetscInt * rows,PetscInt n,const PetscInt * cols,const PetscScalar * values,InsertMode addv)2831 static PetscErrorCode MatSetValuesLocal_IS(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2832 {
2833 Mat_IS *is = (Mat_IS *)A->data;
2834
2835 PetscFunctionBegin;
2836 if (is->A->rmap->mapping || is->A->cmap->mapping) {
2837 PetscCall(MatSetValuesLocal(is->A, m, rows, n, cols, values, addv));
2838 } else {
2839 PetscCall(MatSetValues(is->A, m, rows, n, cols, values, addv));
2840 }
2841 PetscFunctionReturn(PETSC_SUCCESS);
2842 }
2843
MatSetValuesBlockedLocal_IS(Mat A,PetscInt m,const PetscInt * rows,PetscInt n,const PetscInt * cols,const PetscScalar * values,InsertMode addv)2844 static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2845 {
2846 Mat_IS *is = (Mat_IS *)A->data;
2847
2848 PetscFunctionBegin;
2849 if (is->A->rmap->mapping || is->A->cmap->mapping) {
2850 PetscCall(MatSetValuesBlockedLocal(is->A, m, rows, n, cols, values, addv));
2851 } else {
2852 PetscCall(MatSetValuesBlocked(is->A, m, rows, n, cols, values, addv));
2853 }
2854 PetscFunctionReturn(PETSC_SUCCESS);
2855 }
2856
MatISZeroRowsColumnsLocal_Private(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,PetscBool columns)2857 static PetscErrorCode MatISZeroRowsColumnsLocal_Private(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, PetscBool columns)
2858 {
2859 Mat_IS *is = (Mat_IS *)A->data;
2860
2861 PetscFunctionBegin;
2862 if (!n) PetscFunctionReturn(PETSC_SUCCESS);
2863 is->pure_neumann = PETSC_FALSE;
2864
2865 if (columns) {
2866 PetscCall(MatZeroRowsColumns(is->A, n, rows, diag, NULL, NULL));
2867 } else {
2868 PetscCall(MatZeroRows(is->A, n, rows, diag, NULL, NULL));
2869 }
2870 if (diag != 0.) {
2871 const PetscScalar *array;
2872
2873 PetscCall(VecGetArrayRead(is->counter, &array));
2874 for (PetscInt i = 0; i < n; i++) PetscCall(MatSetValue(is->A, rows[i], rows[i], diag / (array[rows[i]]), INSERT_VALUES));
2875 PetscCall(VecRestoreArrayRead(is->counter, &array));
2876 PetscCall(MatAssemblyBegin(is->A, MAT_FINAL_ASSEMBLY));
2877 PetscCall(MatAssemblyEnd(is->A, MAT_FINAL_ASSEMBLY));
2878 }
2879 PetscFunctionReturn(PETSC_SUCCESS);
2880 }
2881
MatZeroRowsColumns_Private_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b,PetscBool columns)2882 static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b, PetscBool columns)
2883 {
2884 Mat_IS *matis = (Mat_IS *)A->data;
2885 PetscInt nr, nl, len;
2886 PetscInt *lrows = NULL;
2887
2888 PetscFunctionBegin;
2889 if (PetscUnlikelyDebug(columns || diag != 0. || (x && b))) {
2890 PetscBool cong;
2891
2892 PetscCall(PetscLayoutCompare(A->rmap, A->cmap, &cong));
2893 cong = (PetscBool)(cong && matis->sf == matis->csf);
2894 PetscCheck(cong || !columns, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Columns can be zeroed if and only if A->rmap and A->cmap are congruent and the l2g maps are the same for MATIS");
2895 PetscCheck(cong || diag == 0., PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Nonzero diagonal value supported if and only if A->rmap and A->cmap are congruent and the l2g maps are the same for MATIS");
2896 PetscCheck(cong || !x || !b, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "A->rmap and A->cmap need to be congruent, and the l2g maps be the same");
2897 }
2898 PetscCall(MatGetSize(matis->A, &nl, NULL));
2899 /* get locally owned rows */
2900 PetscCall(PetscLayoutMapLocal(A->rmap, n, rows, &len, &lrows, NULL));
2901 /* fix right-hand side if needed */
2902 if (x && b) {
2903 const PetscScalar *xx;
2904 PetscScalar *bb;
2905
2906 PetscCall(VecGetArrayRead(x, &xx));
2907 PetscCall(VecGetArray(b, &bb));
2908 for (PetscInt i = 0; i < len; ++i) bb[lrows[i]] = diag * xx[lrows[i]];
2909 PetscCall(VecRestoreArrayRead(x, &xx));
2910 PetscCall(VecRestoreArray(b, &bb));
2911 }
2912 /* get rows associated to the local matrices */
2913 PetscCall(PetscArrayzero(matis->sf_leafdata, nl));
2914 PetscCall(PetscArrayzero(matis->sf_rootdata, A->rmap->n));
2915 for (PetscInt i = 0; i < len; i++) matis->sf_rootdata[lrows[i]] = 1;
2916 PetscCall(PetscFree(lrows));
2917 PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
2918 PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
2919 PetscCall(PetscMalloc1(nl, &lrows));
2920 nr = 0;
2921 for (PetscInt i = 0; i < nl; i++)
2922 if (matis->sf_leafdata[i]) lrows[nr++] = i;
2923 PetscCall(MatISZeroRowsColumnsLocal_Private(A, nr, lrows, diag, columns));
2924 PetscCall(PetscFree(lrows));
2925 PetscFunctionReturn(PETSC_SUCCESS);
2926 }
2927
MatZeroRows_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)2928 static PetscErrorCode MatZeroRows_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
2929 {
2930 PetscFunctionBegin;
2931 PetscCall(MatZeroRowsColumns_Private_IS(A, n, rows, diag, x, b, PETSC_FALSE));
2932 PetscFunctionReturn(PETSC_SUCCESS);
2933 }
2934
MatZeroRowsColumns_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)2935 static PetscErrorCode MatZeroRowsColumns_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
2936 {
2937 PetscFunctionBegin;
2938 PetscCall(MatZeroRowsColumns_Private_IS(A, n, rows, diag, x, b, PETSC_TRUE));
2939 PetscFunctionReturn(PETSC_SUCCESS);
2940 }
2941
MatAssemblyBegin_IS(Mat A,MatAssemblyType type)2942 static PetscErrorCode MatAssemblyBegin_IS(Mat A, MatAssemblyType type)
2943 {
2944 Mat_IS *is = (Mat_IS *)A->data;
2945
2946 PetscFunctionBegin;
2947 PetscCall(MatAssemblyBegin(is->A, type));
2948 PetscFunctionReturn(PETSC_SUCCESS);
2949 }
2950
MatAssemblyEnd_IS(Mat A,MatAssemblyType type)2951 static PetscErrorCode MatAssemblyEnd_IS(Mat A, MatAssemblyType type)
2952 {
2953 Mat_IS *is = (Mat_IS *)A->data;
2954 PetscBool lnnz;
2955
2956 PetscFunctionBegin;
2957 PetscCall(MatAssemblyEnd(is->A, type));
2958 /* fix for local empty rows/cols */
2959 if (is->locempty && type == MAT_FINAL_ASSEMBLY) {
2960 Mat newlA;
2961 ISLocalToGlobalMapping rl2g, cl2g;
2962 IS nzr, nzc;
2963 PetscInt nr, nc, nnzr, nnzc;
2964 PetscBool lnewl2g, newl2g;
2965
2966 PetscCall(MatGetSize(is->A, &nr, &nc));
2967 PetscCall(MatFindNonzeroRowsOrCols_Basic(is->A, PETSC_FALSE, PETSC_SMALL, &nzr));
2968 if (!nzr) PetscCall(ISCreateStride(PetscObjectComm((PetscObject)is->A), nr, 0, 1, &nzr));
2969 PetscCall(MatFindNonzeroRowsOrCols_Basic(is->A, PETSC_TRUE, PETSC_SMALL, &nzc));
2970 if (!nzc) PetscCall(ISCreateStride(PetscObjectComm((PetscObject)is->A), nc, 0, 1, &nzc));
2971 PetscCall(ISGetSize(nzr, &nnzr));
2972 PetscCall(ISGetSize(nzc, &nnzc));
2973 if (nnzr != nr || nnzc != nc) { /* need new global l2g map */
2974 lnewl2g = PETSC_TRUE;
2975 PetscCallMPI(MPIU_Allreduce(&lnewl2g, &newl2g, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A)));
2976
2977 /* extract valid submatrix */
2978 PetscCall(MatCreateSubMatrix(is->A, nzr, nzc, MAT_INITIAL_MATRIX, &newlA));
2979 } else { /* local matrix fully populated */
2980 lnewl2g = PETSC_FALSE;
2981 PetscCallMPI(MPIU_Allreduce(&lnewl2g, &newl2g, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A)));
2982 PetscCall(PetscObjectReference((PetscObject)is->A));
2983 newlA = is->A;
2984 }
2985
2986 /* attach new global l2g map if needed */
2987 if (newl2g) {
2988 IS zr, zc;
2989 const PetscInt *ridxs, *cidxs, *zridxs, *zcidxs;
2990 PetscInt *nidxs, i;
2991
2992 PetscCall(ISComplement(nzr, 0, nr, &zr));
2993 PetscCall(ISComplement(nzc, 0, nc, &zc));
2994 PetscCall(PetscMalloc1(PetscMax(nr, nc), &nidxs));
2995 PetscCall(ISLocalToGlobalMappingGetIndices(is->rmapping, &ridxs));
2996 PetscCall(ISLocalToGlobalMappingGetIndices(is->cmapping, &cidxs));
2997 PetscCall(ISGetIndices(zr, &zridxs));
2998 PetscCall(ISGetIndices(zc, &zcidxs));
2999 PetscCall(ISGetLocalSize(zr, &nnzr));
3000 PetscCall(ISGetLocalSize(zc, &nnzc));
3001
3002 PetscCall(PetscArraycpy(nidxs, ridxs, nr));
3003 for (i = 0; i < nnzr; i++) nidxs[zridxs[i]] = -1;
3004 PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A), 1, nr, nidxs, PETSC_COPY_VALUES, &rl2g));
3005 PetscCall(PetscArraycpy(nidxs, cidxs, nc));
3006 for (i = 0; i < nnzc; i++) nidxs[zcidxs[i]] = -1;
3007 PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A), 1, nc, nidxs, PETSC_COPY_VALUES, &cl2g));
3008
3009 PetscCall(ISRestoreIndices(zr, &zridxs));
3010 PetscCall(ISRestoreIndices(zc, &zcidxs));
3011 PetscCall(ISLocalToGlobalMappingRestoreIndices(is->rmapping, &ridxs));
3012 PetscCall(ISLocalToGlobalMappingRestoreIndices(is->cmapping, &cidxs));
3013 PetscCall(ISDestroy(&nzr));
3014 PetscCall(ISDestroy(&nzc));
3015 PetscCall(ISDestroy(&zr));
3016 PetscCall(ISDestroy(&zc));
3017 PetscCall(PetscFree(nidxs));
3018 PetscCall(MatSetLocalToGlobalMapping(A, rl2g, cl2g));
3019 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
3020 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
3021 }
3022 PetscCall(MatISSetLocalMat(A, newlA));
3023 PetscCall(MatDestroy(&newlA));
3024 PetscCall(ISDestroy(&nzr));
3025 PetscCall(ISDestroy(&nzc));
3026 is->locempty = PETSC_FALSE;
3027 }
3028 lnnz = (PetscBool)(is->A->nonzerostate == is->lnnzstate);
3029 is->lnnzstate = is->A->nonzerostate;
3030 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &lnnz, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
3031 if (!lnnz) A->nonzerostate++;
3032 PetscFunctionReturn(PETSC_SUCCESS);
3033 }
3034
MatISGetLocalMat_IS(Mat mat,Mat * local)3035 static PetscErrorCode MatISGetLocalMat_IS(Mat mat, Mat *local)
3036 {
3037 Mat_IS *is = (Mat_IS *)mat->data;
3038
3039 PetscFunctionBegin;
3040 *local = is->A;
3041 PetscFunctionReturn(PETSC_SUCCESS);
3042 }
3043
MatISRestoreLocalMat_IS(Mat mat,Mat * local)3044 static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat, Mat *local)
3045 {
3046 PetscFunctionBegin;
3047 *local = NULL;
3048 PetscFunctionReturn(PETSC_SUCCESS);
3049 }
3050
3051 /*@
3052 MatISGetLocalMat - Gets the local matrix stored inside a `MATIS` matrix.
3053
3054 Not Collective.
3055
3056 Input Parameter:
3057 . mat - the matrix
3058
3059 Output Parameter:
3060 . local - the local matrix
3061
3062 Level: intermediate
3063
3064 Notes:
3065 This can be called if you have precomputed the nonzero structure of the
3066 matrix and want to provide it to the inner matrix object to improve the performance
3067 of the `MatSetValues()` operation.
3068
3069 Call `MatISRestoreLocalMat()` when finished with the local matrix.
3070
3071 .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatISRestoreLocalMat()`
3072 @*/
MatISGetLocalMat(Mat mat,Mat * local)3073 PetscErrorCode MatISGetLocalMat(Mat mat, Mat *local)
3074 {
3075 PetscFunctionBegin;
3076 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1);
3077 PetscAssertPointer(local, 2);
3078 PetscUseMethod(mat, "MatISGetLocalMat_C", (Mat, Mat *), (mat, local));
3079 PetscFunctionReturn(PETSC_SUCCESS);
3080 }
3081
3082 /*@
3083 MatISRestoreLocalMat - Restores the local matrix obtained with `MatISGetLocalMat()`
3084
3085 Not Collective.
3086
3087 Input Parameters:
3088 + mat - the matrix
3089 - local - the local matrix
3090
3091 Level: intermediate
3092
3093 .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatISGetLocalMat()`
3094 @*/
MatISRestoreLocalMat(Mat mat,Mat * local)3095 PetscErrorCode MatISRestoreLocalMat(Mat mat, Mat *local)
3096 {
3097 PetscFunctionBegin;
3098 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1);
3099 PetscAssertPointer(local, 2);
3100 PetscUseMethod(mat, "MatISRestoreLocalMat_C", (Mat, Mat *), (mat, local));
3101 PetscFunctionReturn(PETSC_SUCCESS);
3102 }
3103
MatISSetLocalMatType_IS(Mat mat,MatType mtype)3104 static PetscErrorCode MatISSetLocalMatType_IS(Mat mat, MatType mtype)
3105 {
3106 Mat_IS *is = (Mat_IS *)mat->data;
3107
3108 PetscFunctionBegin;
3109 if (is->A) PetscCall(MatSetType(is->A, mtype));
3110 PetscCall(PetscFree(is->lmattype));
3111 PetscCall(PetscStrallocpy(mtype, &is->lmattype));
3112 PetscFunctionReturn(PETSC_SUCCESS);
3113 }
3114
3115 /*@
3116 MatISSetLocalMatType - Specifies the type of local matrix inside the `MATIS`
3117
3118 Logically Collective.
3119
3120 Input Parameters:
3121 + mat - the matrix
3122 - mtype - the local matrix type
3123
3124 Level: intermediate
3125
3126 .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatSetType()`, `MatType`
3127 @*/
MatISSetLocalMatType(Mat mat,MatType mtype)3128 PetscErrorCode MatISSetLocalMatType(Mat mat, MatType mtype)
3129 {
3130 PetscFunctionBegin;
3131 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1);
3132 PetscUseMethod(mat, "MatISSetLocalMatType_C", (Mat, MatType), (mat, mtype));
3133 PetscFunctionReturn(PETSC_SUCCESS);
3134 }
3135
MatISSetLocalMat_IS(Mat mat,Mat local)3136 static PetscErrorCode MatISSetLocalMat_IS(Mat mat, Mat local)
3137 {
3138 Mat_IS *is = (Mat_IS *)mat->data;
3139 PetscInt nrows, ncols, orows, ocols;
3140 MatType mtype, otype;
3141 PetscBool sametype = PETSC_TRUE;
3142
3143 PetscFunctionBegin;
3144 if (is->A && !is->islocalref) {
3145 PetscCall(MatGetSize(is->A, &orows, &ocols));
3146 PetscCall(MatGetSize(local, &nrows, &ncols));
3147 PetscCheck(orows == nrows && ocols == ncols, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local MATIS matrix should be of size %" PetscInt_FMT "x%" PetscInt_FMT " (passed a %" PetscInt_FMT "x%" PetscInt_FMT " matrix)", orows, ocols, nrows, ncols);
3148 PetscCall(MatGetType(local, &mtype));
3149 PetscCall(MatGetType(is->A, &otype));
3150 PetscCall(PetscStrcmp(mtype, otype, &sametype));
3151 }
3152 PetscCall(PetscObjectReference((PetscObject)local));
3153 PetscCall(MatDestroy(&is->A));
3154 is->A = local;
3155 PetscCall(MatGetType(is->A, &mtype));
3156 PetscCall(MatISSetLocalMatType(mat, mtype));
3157 if (!sametype && !is->islocalref) PetscCall(MatISSetUpScatters_Private(mat));
3158 is->lnnzstate = 0;
3159 PetscFunctionReturn(PETSC_SUCCESS);
3160 }
3161
3162 /*@
3163 MatISSetLocalMat - Replace the local matrix stored inside a `MATIS` object.
3164
3165 Not Collective
3166
3167 Input Parameters:
3168 + mat - the matrix
3169 - local - the local matrix
3170
3171 Level: intermediate
3172
3173 .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatISSetLocalMatType`, `MatISGetLocalMat()`
3174 @*/
MatISSetLocalMat(Mat mat,Mat local)3175 PetscErrorCode MatISSetLocalMat(Mat mat, Mat local)
3176 {
3177 PetscFunctionBegin;
3178 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1);
3179 PetscValidHeaderSpecific(local, MAT_CLASSID, 2);
3180 PetscUseMethod(mat, "MatISSetLocalMat_C", (Mat, Mat), (mat, local));
3181 PetscFunctionReturn(PETSC_SUCCESS);
3182 }
3183
MatZeroEntries_IS(Mat A)3184 static PetscErrorCode MatZeroEntries_IS(Mat A)
3185 {
3186 Mat_IS *a = (Mat_IS *)A->data;
3187
3188 PetscFunctionBegin;
3189 PetscCall(MatZeroEntries(a->A));
3190 PetscFunctionReturn(PETSC_SUCCESS);
3191 }
3192
MatScale_IS(Mat A,PetscScalar a)3193 static PetscErrorCode MatScale_IS(Mat A, PetscScalar a)
3194 {
3195 Mat_IS *is = (Mat_IS *)A->data;
3196
3197 PetscFunctionBegin;
3198 PetscCall(MatScale(is->A, a));
3199 PetscFunctionReturn(PETSC_SUCCESS);
3200 }
3201
MatGetDiagonal_IS(Mat A,Vec v)3202 static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
3203 {
3204 Mat_IS *is = (Mat_IS *)A->data;
3205
3206 PetscFunctionBegin;
3207 /* get diagonal of the local matrix */
3208 PetscCall(MatGetDiagonal(is->A, is->y));
3209
3210 /* scatter diagonal back into global vector */
3211 PetscCall(VecSet(v, 0));
3212 PetscCall(VecScatterBegin(is->rctx, is->y, v, ADD_VALUES, SCATTER_REVERSE));
3213 PetscCall(VecScatterEnd(is->rctx, is->y, v, ADD_VALUES, SCATTER_REVERSE));
3214 PetscFunctionReturn(PETSC_SUCCESS);
3215 }
3216
MatSetOption_IS(Mat A,MatOption op,PetscBool flg)3217 static PetscErrorCode MatSetOption_IS(Mat A, MatOption op, PetscBool flg)
3218 {
3219 Mat_IS *a = (Mat_IS *)A->data;
3220
3221 PetscFunctionBegin;
3222 PetscCall(MatSetOption(a->A, op, flg));
3223 PetscFunctionReturn(PETSC_SUCCESS);
3224 }
3225
MatAXPY_IS(Mat Y,PetscScalar a,Mat X,MatStructure str)3226 static PetscErrorCode MatAXPY_IS(Mat Y, PetscScalar a, Mat X, MatStructure str)
3227 {
3228 Mat_IS *y = (Mat_IS *)Y->data;
3229 Mat_IS *x;
3230
3231 PetscFunctionBegin;
3232 if (PetscDefined(USE_DEBUG)) {
3233 PetscBool ismatis;
3234 PetscCall(PetscObjectTypeCompare((PetscObject)X, MATIS, &ismatis));
3235 PetscCheck(ismatis, PetscObjectComm((PetscObject)Y), PETSC_ERR_SUP, "Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS");
3236 }
3237 x = (Mat_IS *)X->data;
3238 PetscCall(MatAXPY(y->A, a, x->A, str));
3239 PetscFunctionReturn(PETSC_SUCCESS);
3240 }
3241
MatGetLocalSubMatrix_IS(Mat A,IS row,IS col,Mat * submat)3242 static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A, IS row, IS col, Mat *submat)
3243 {
3244 Mat lA;
3245 Mat_IS *matis = (Mat_IS *)A->data;
3246 ISLocalToGlobalMapping rl2g, cl2g;
3247 IS is;
3248 const PetscInt *rg, *rl;
3249 PetscInt nrg, rbs, cbs;
3250 PetscInt N, M, nrl, i, *idxs;
3251
3252 PetscFunctionBegin;
3253 PetscCall(ISGetBlockSize(row, &rbs));
3254 PetscCall(ISGetBlockSize(col, &cbs));
3255 PetscCall(ISLocalToGlobalMappingGetIndices(A->rmap->mapping, &rg));
3256 PetscCall(ISGetLocalSize(row, &nrl));
3257 PetscCall(ISGetIndices(row, &rl));
3258 PetscCall(ISLocalToGlobalMappingGetSize(A->rmap->mapping, &nrg));
3259 if (PetscDefined(USE_DEBUG)) {
3260 for (i = 0; i < nrl; i++) PetscCheck(rl[i] < nrg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Local row index %" PetscInt_FMT " -> %" PetscInt_FMT " greater than maximum possible %" PetscInt_FMT, i, rl[i], nrg);
3261 }
3262 if (nrg % rbs) nrg = rbs * (nrg / rbs + 1);
3263 PetscCall(PetscMalloc1(nrg, &idxs));
3264 /* map from [0,nrl) to row */
3265 for (i = 0; i < nrl; i++) idxs[i] = rl[i];
3266 for (i = nrl; i < nrg; i++) idxs[i] = -1;
3267 PetscCall(ISRestoreIndices(row, &rl));
3268 PetscCall(ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping, &rg));
3269 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), nrg, idxs, PETSC_OWN_POINTER, &is));
3270 PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
3271 PetscCall(ISLocalToGlobalMappingSetBlockSize(rl2g, rbs));
3272 PetscCall(ISDestroy(&is));
3273 /* compute new l2g map for columns */
3274 if (col != row || matis->rmapping != matis->cmapping || matis->A->rmap->mapping != matis->A->cmap->mapping) {
3275 const PetscInt *cg, *cl;
3276 PetscInt ncg;
3277 PetscInt ncl;
3278
3279 PetscCall(ISLocalToGlobalMappingGetIndices(A->cmap->mapping, &cg));
3280 PetscCall(ISGetLocalSize(col, &ncl));
3281 PetscCall(ISGetIndices(col, &cl));
3282 PetscCall(ISLocalToGlobalMappingGetSize(A->cmap->mapping, &ncg));
3283 if (PetscDefined(USE_DEBUG)) {
3284 for (i = 0; i < ncl; i++) PetscCheck(cl[i] < ncg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Local column index %" PetscInt_FMT " -> %" PetscInt_FMT " greater than maximum possible %" PetscInt_FMT, i, cl[i], ncg);
3285 }
3286 if (ncg % cbs) ncg = cbs * (ncg / cbs + 1);
3287 PetscCall(PetscMalloc1(ncg, &idxs));
3288 /* map from [0,ncl) to col */
3289 for (i = 0; i < ncl; i++) idxs[i] = cl[i];
3290 for (i = ncl; i < ncg; i++) idxs[i] = -1;
3291 PetscCall(ISRestoreIndices(col, &cl));
3292 PetscCall(ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping, &cg));
3293 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), ncg, idxs, PETSC_OWN_POINTER, &is));
3294 PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
3295 PetscCall(ISLocalToGlobalMappingSetBlockSize(cl2g, cbs));
3296 PetscCall(ISDestroy(&is));
3297 } else {
3298 PetscCall(PetscObjectReference((PetscObject)rl2g));
3299 cl2g = rl2g;
3300 }
3301
3302 /* create the MATIS submatrix */
3303 PetscCall(MatGetSize(A, &M, &N));
3304 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), submat));
3305 PetscCall(MatSetSizes(*submat, PETSC_DECIDE, PETSC_DECIDE, M, N));
3306 PetscCall(MatSetType(*submat, MATIS));
3307 matis = (Mat_IS *)((*submat)->data);
3308 matis->islocalref = A;
3309 PetscCall(MatSetLocalToGlobalMapping(*submat, rl2g, cl2g));
3310 PetscCall(MatISGetLocalMat(A, &lA));
3311 PetscCall(MatISSetLocalMat(*submat, lA));
3312 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
3313 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
3314
3315 /* remove unsupported ops */
3316 PetscCall(PetscMemzero((*submat)->ops, sizeof(struct _MatOps)));
3317 (*submat)->ops->destroy = MatDestroy_IS;
3318 (*submat)->ops->setvalueslocal = MatSetValuesLocal_SubMat_IS;
3319 (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS;
3320 (*submat)->ops->zerorowslocal = MatZeroRowsLocal_SubMat_IS;
3321 (*submat)->ops->zerorowscolumnslocal = MatZeroRowsColumnsLocal_SubMat_IS;
3322 (*submat)->ops->getlocalsubmatrix = MatGetLocalSubMatrix_IS;
3323 PetscFunctionReturn(PETSC_SUCCESS);
3324 }
3325
MatSetFromOptions_IS(Mat A,PetscOptionItems PetscOptionsObject)3326 static PetscErrorCode MatSetFromOptions_IS(Mat A, PetscOptionItems PetscOptionsObject)
3327 {
3328 Mat_IS *a = (Mat_IS *)A->data;
3329 char type[256];
3330 PetscBool flg;
3331
3332 PetscFunctionBegin;
3333 PetscOptionsHeadBegin(PetscOptionsObject, "MATIS options");
3334 PetscCall(PetscOptionsDeprecated("-matis_keepassembled", "-mat_is_keepassembled", "3.21", NULL));
3335 PetscCall(PetscOptionsDeprecated("-matis_fixempty", "-mat_is_fixempty", "3.21", NULL));
3336 PetscCall(PetscOptionsDeprecated("-matis_storel2l", "-mat_is_storel2l", "3.21", NULL));
3337 PetscCall(PetscOptionsDeprecated("-matis_localmat_type", "-mat_is_localmat_type", "3.21", NULL));
3338 PetscCall(PetscOptionsBool("-mat_is_keepassembled", "Store an assembled version if needed", NULL, a->keepassembled, &a->keepassembled, NULL));
3339 PetscCall(PetscOptionsBool("-mat_is_fixempty", "Fix local matrices in case of empty local rows/columns", "MatISFixLocalEmpty", a->locempty, &a->locempty, NULL));
3340 PetscCall(PetscOptionsBool("-mat_is_storel2l", "Store local-to-local matrices generated from PtAP operations", "MatISStoreL2L", a->storel2l, &a->storel2l, NULL));
3341 PetscCall(PetscOptionsBool("-mat_is_allow_repeated", "Allow local repeated entries", "MatISSetAllowRepeated", a->allow_repeated, &a->allow_repeated, NULL));
3342 PetscCall(PetscOptionsFList("-mat_is_localmat_type", "Matrix type", "MatISSetLocalMatType", MatList, a->lmattype, type, 256, &flg));
3343 if (flg) PetscCall(MatISSetLocalMatType(A, type));
3344 if (a->A) PetscCall(MatSetFromOptions(a->A));
3345 PetscOptionsHeadEnd();
3346 PetscFunctionReturn(PETSC_SUCCESS);
3347 }
3348
3349 /*@
3350 MatCreateIS - Creates a "process" unassembled matrix.
3351
3352 Collective.
3353
3354 Input Parameters:
3355 + comm - MPI communicator that will share the matrix
3356 . bs - block size of the matrix
3357 . m - local size of left vector used in matrix vector products
3358 . n - local size of right vector used in matrix vector products
3359 . M - global size of left vector used in matrix vector products
3360 . N - global size of right vector used in matrix vector products
3361 . rmap - local to global map for rows
3362 - cmap - local to global map for cols
3363
3364 Output Parameter:
3365 . A - the resulting matrix
3366
3367 Level: intermediate
3368
3369 Notes:
3370 `m` and `n` are NOT related to the size of the map; they represent the size of the local parts of the distributed vectors
3371 used in `MatMult()` operations. The local sizes of `rmap` and `cmap` define the size of the local matrices.
3372
3373 If `rmap` (`cmap`) is `NULL`, then the local row (column) spaces matches the global space.
3374
3375 .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatSetLocalToGlobalMapping()`
3376 @*/
MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat * A)3377 PetscErrorCode MatCreateIS(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt M, PetscInt N, ISLocalToGlobalMapping rmap, ISLocalToGlobalMapping cmap, Mat *A)
3378 {
3379 PetscFunctionBegin;
3380 PetscCall(MatCreate(comm, A));
3381 PetscCall(MatSetSizes(*A, m, n, M, N));
3382 if (bs > 0) PetscCall(MatSetBlockSize(*A, bs));
3383 PetscCall(MatSetType(*A, MATIS));
3384 PetscCall(MatSetLocalToGlobalMapping(*A, rmap, cmap));
3385 PetscFunctionReturn(PETSC_SUCCESS);
3386 }
3387
MatHasOperation_IS(Mat A,MatOperation op,PetscBool * has)3388 static PetscErrorCode MatHasOperation_IS(Mat A, MatOperation op, PetscBool *has)
3389 {
3390 Mat_IS *a = (Mat_IS *)A->data;
3391 MatOperation tobefiltered[] = {MATOP_MULT_ADD, MATOP_MULT_TRANSPOSE_ADD, MATOP_GET_DIAGONAL_BLOCK, MATOP_INCREASE_OVERLAP};
3392
3393 PetscFunctionBegin;
3394 *has = PETSC_FALSE;
3395 if (!((void **)A->ops)[op] || !a->A) PetscFunctionReturn(PETSC_SUCCESS);
3396 *has = PETSC_TRUE;
3397 for (PetscInt i = 0; i < (PetscInt)PETSC_STATIC_ARRAY_LENGTH(tobefiltered); i++)
3398 if (op == tobefiltered[i]) PetscFunctionReturn(PETSC_SUCCESS);
3399 PetscCall(MatHasOperation(a->A, op, has));
3400 PetscFunctionReturn(PETSC_SUCCESS);
3401 }
3402
MatSetValuesCOO_IS(Mat A,const PetscScalar v[],InsertMode imode)3403 static PetscErrorCode MatSetValuesCOO_IS(Mat A, const PetscScalar v[], InsertMode imode)
3404 {
3405 Mat_IS *a = (Mat_IS *)A->data;
3406
3407 PetscFunctionBegin;
3408 PetscCall(MatSetValuesCOO(a->A, v, imode));
3409 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
3410 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
3411 PetscFunctionReturn(PETSC_SUCCESS);
3412 }
3413
MatSetPreallocationCOOLocal_IS(Mat A,PetscCount ncoo,PetscInt coo_i[],PetscInt coo_j[])3414 static PetscErrorCode MatSetPreallocationCOOLocal_IS(Mat A, PetscCount ncoo, PetscInt coo_i[], PetscInt coo_j[])
3415 {
3416 Mat_IS *a = (Mat_IS *)A->data;
3417
3418 PetscFunctionBegin;
3419 PetscCheck(a->A, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to provide l2g map first via MatSetLocalToGlobalMapping");
3420 if (a->A->rmap->mapping || a->A->cmap->mapping) {
3421 PetscCall(MatSetPreallocationCOOLocal(a->A, ncoo, coo_i, coo_j));
3422 } else {
3423 PetscCall(MatSetPreallocationCOO(a->A, ncoo, coo_i, coo_j));
3424 }
3425 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", MatSetValuesCOO_IS));
3426 A->preallocated = PETSC_TRUE;
3427 PetscFunctionReturn(PETSC_SUCCESS);
3428 }
3429
MatSetPreallocationCOO_IS(Mat A,PetscCount ncoo,PetscInt coo_i[],PetscInt coo_j[])3430 static PetscErrorCode MatSetPreallocationCOO_IS(Mat A, PetscCount ncoo, PetscInt coo_i[], PetscInt coo_j[])
3431 {
3432 Mat_IS *a = (Mat_IS *)A->data;
3433 PetscInt ncoo_i;
3434
3435 PetscFunctionBegin;
3436 PetscCheck(a->A, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to provide l2g map first via MatSetLocalToGlobalMapping");
3437 PetscCall(PetscIntCast(ncoo, &ncoo_i));
3438 PetscCall(ISGlobalToLocalMappingApply(a->rmapping, IS_GTOLM_MASK, ncoo_i, coo_i, NULL, coo_i));
3439 PetscCall(ISGlobalToLocalMappingApply(a->cmapping, IS_GTOLM_MASK, ncoo_i, coo_j, NULL, coo_j));
3440 PetscCall(MatSetPreallocationCOO(a->A, ncoo, coo_i, coo_j));
3441 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", MatSetValuesCOO_IS));
3442 A->preallocated = PETSC_TRUE;
3443 PetscFunctionReturn(PETSC_SUCCESS);
3444 }
3445
MatISGetAssembled_Private(Mat A,Mat * tA)3446 static PetscErrorCode MatISGetAssembled_Private(Mat A, Mat *tA)
3447 {
3448 Mat_IS *a = (Mat_IS *)A->data;
3449 PetscObjectState Astate, aAstate = PETSC_INT_MIN;
3450 PetscObjectState Annzstate, aAnnzstate = PETSC_INT_MIN;
3451
3452 PetscFunctionBegin;
3453 PetscCall(PetscObjectStateGet((PetscObject)A, &Astate));
3454 Annzstate = A->nonzerostate;
3455 if (a->assembledA) {
3456 PetscCall(PetscObjectStateGet((PetscObject)a->assembledA, &aAstate));
3457 aAnnzstate = a->assembledA->nonzerostate;
3458 }
3459 if (aAnnzstate != Annzstate) PetscCall(MatDestroy(&a->assembledA));
3460 if (Astate != aAstate || !a->assembledA) {
3461 MatType aAtype;
3462 PetscMPIInt size;
3463 PetscInt rbs, cbs, bs;
3464
3465 /* the assembled form is used as temporary storage for parallel operations
3466 like createsubmatrices and the like, do not waste device memory */
3467 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
3468 PetscCall(ISLocalToGlobalMappingGetBlockSize(a->cmapping, &cbs));
3469 PetscCall(ISLocalToGlobalMappingGetBlockSize(a->rmapping, &rbs));
3470 bs = rbs == cbs ? rbs : 1;
3471 if (a->assembledA) PetscCall(MatGetType(a->assembledA, &aAtype));
3472 else if (size > 1) aAtype = bs > 1 ? MATMPIBAIJ : MATMPIAIJ;
3473 else aAtype = bs > 1 ? MATSEQBAIJ : MATSEQAIJ;
3474
3475 PetscCall(MatConvert(A, aAtype, a->assembledA ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, &a->assembledA));
3476 PetscCall(PetscObjectStateSet((PetscObject)a->assembledA, Astate));
3477 a->assembledA->nonzerostate = Annzstate;
3478 }
3479 PetscCall(PetscObjectReference((PetscObject)a->assembledA));
3480 *tA = a->assembledA;
3481 if (!a->keepassembled) PetscCall(MatDestroy(&a->assembledA));
3482 PetscFunctionReturn(PETSC_SUCCESS);
3483 }
3484
MatISRestoreAssembled_Private(Mat A,Mat * tA)3485 static PetscErrorCode MatISRestoreAssembled_Private(Mat A, Mat *tA)
3486 {
3487 PetscFunctionBegin;
3488 PetscCall(MatDestroy(tA));
3489 *tA = NULL;
3490 PetscFunctionReturn(PETSC_SUCCESS);
3491 }
3492
MatGetDiagonalBlock_IS(Mat A,Mat * dA)3493 static PetscErrorCode MatGetDiagonalBlock_IS(Mat A, Mat *dA)
3494 {
3495 Mat_IS *a = (Mat_IS *)A->data;
3496 PetscObjectState Astate, dAstate = PETSC_INT_MIN;
3497
3498 PetscFunctionBegin;
3499 PetscCall(PetscObjectStateGet((PetscObject)A, &Astate));
3500 if (a->dA) PetscCall(PetscObjectStateGet((PetscObject)a->dA, &dAstate));
3501 if (Astate != dAstate) {
3502 Mat tA;
3503 MatType ltype;
3504
3505 PetscCall(MatDestroy(&a->dA));
3506 PetscCall(MatISGetAssembled_Private(A, &tA));
3507 PetscCall(MatGetDiagonalBlock(tA, &a->dA));
3508 PetscCall(MatPropagateSymmetryOptions(tA, a->dA));
3509 PetscCall(MatGetType(a->A, <ype));
3510 PetscCall(MatConvert(a->dA, ltype, MAT_INPLACE_MATRIX, &a->dA));
3511 PetscCall(PetscObjectReference((PetscObject)a->dA));
3512 PetscCall(MatISRestoreAssembled_Private(A, &tA));
3513 PetscCall(PetscObjectStateSet((PetscObject)a->dA, Astate));
3514 }
3515 *dA = a->dA;
3516 PetscFunctionReturn(PETSC_SUCCESS);
3517 }
3518
MatCreateSubMatrices_IS(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse reuse,Mat * submat[])3519 static PetscErrorCode MatCreateSubMatrices_IS(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse reuse, Mat *submat[])
3520 {
3521 Mat tA;
3522
3523 PetscFunctionBegin;
3524 PetscCall(MatISGetAssembled_Private(A, &tA));
3525 PetscCall(MatCreateSubMatrices(tA, n, irow, icol, reuse, submat));
3526 /* MatCreateSubMatrices_MPIAIJ is a mess at the moment */
3527 #if 0
3528 {
3529 Mat_IS *a = (Mat_IS*)A->data;
3530 MatType ltype;
3531 VecType vtype;
3532 char *flg;
3533
3534 PetscCall(MatGetType(a->A,<ype));
3535 PetscCall(MatGetVecType(a->A,&vtype));
3536 PetscCall(PetscStrstr(vtype,"cuda",&flg));
3537 if (!flg) PetscCall(PetscStrstr(vtype,"hip",&flg));
3538 if (!flg) PetscCall(PetscStrstr(vtype,"kokkos",&flg));
3539 if (flg) {
3540 for (PetscInt i = 0; i < n; i++) {
3541 Mat sA = (*submat)[i];
3542
3543 PetscCall(MatConvert(sA,ltype,MAT_INPLACE_MATRIX,&sA));
3544 (*submat)[i] = sA;
3545 }
3546 }
3547 }
3548 #endif
3549 PetscCall(MatISRestoreAssembled_Private(A, &tA));
3550 PetscFunctionReturn(PETSC_SUCCESS);
3551 }
3552
MatIncreaseOverlap_IS(Mat A,PetscInt n,IS is[],PetscInt ov)3553 static PetscErrorCode MatIncreaseOverlap_IS(Mat A, PetscInt n, IS is[], PetscInt ov)
3554 {
3555 Mat tA;
3556
3557 PetscFunctionBegin;
3558 PetscCall(MatISGetAssembled_Private(A, &tA));
3559 PetscCall(MatIncreaseOverlap(tA, n, is, ov));
3560 PetscCall(MatISRestoreAssembled_Private(A, &tA));
3561 PetscFunctionReturn(PETSC_SUCCESS);
3562 }
3563
3564 /*@
3565 MatISGetLocalToGlobalMapping - Gets the local-to-global numbering of the `MATIS` object
3566
3567 Not Collective
3568
3569 Input Parameter:
3570 . A - the matrix
3571
3572 Output Parameters:
3573 + rmapping - row mapping
3574 - cmapping - column mapping
3575
3576 Level: advanced
3577
3578 Note:
3579 The returned map can be different from the one used to construct the `MATIS` object, since it will not contain negative or repeated indices.
3580
3581 .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatSetLocalToGlobalMapping()`
3582 @*/
MatISGetLocalToGlobalMapping(Mat A,ISLocalToGlobalMapping * rmapping,ISLocalToGlobalMapping * cmapping)3583 PetscErrorCode MatISGetLocalToGlobalMapping(Mat A, ISLocalToGlobalMapping *rmapping, ISLocalToGlobalMapping *cmapping)
3584 {
3585 PetscFunctionBegin;
3586 PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
3587 PetscValidType(A, 1);
3588 if (rmapping) PetscAssertPointer(rmapping, 2);
3589 if (cmapping) PetscAssertPointer(cmapping, 3);
3590 PetscUseMethod(A, "MatISGetLocalToGlobalMapping_C", (Mat, ISLocalToGlobalMapping *, ISLocalToGlobalMapping *), (A, rmapping, cmapping));
3591 PetscFunctionReturn(PETSC_SUCCESS);
3592 }
3593
MatISGetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping * r,ISLocalToGlobalMapping * c)3594 static PetscErrorCode MatISGetLocalToGlobalMapping_IS(Mat A, ISLocalToGlobalMapping *r, ISLocalToGlobalMapping *c)
3595 {
3596 Mat_IS *a = (Mat_IS *)A->data;
3597
3598 PetscFunctionBegin;
3599 if (r) *r = a->rmapping;
3600 if (c) *c = a->cmapping;
3601 PetscFunctionReturn(PETSC_SUCCESS);
3602 }
3603
MatSetBlockSizes_IS(Mat A,PetscInt rbs,PetscInt cbs)3604 static PetscErrorCode MatSetBlockSizes_IS(Mat A, PetscInt rbs, PetscInt cbs)
3605 {
3606 Mat_IS *a = (Mat_IS *)A->data;
3607
3608 PetscFunctionBegin;
3609 if (a->A) PetscCall(MatSetBlockSizes(a->A, rbs, cbs));
3610 PetscFunctionReturn(PETSC_SUCCESS);
3611 }
3612
3613 /*MC
3614 MATIS - MATIS = "is" - A matrix type to be used for non-overlapping domain decomposition methods (e.g. `PCBDDC` or `KSPFETIDP`).
3615 This stores the matrices in globally unassembled form and the parallel matrix vector product is handled "implicitly".
3616
3617 Options Database Keys:
3618 + -mat_type is - Set the matrix type to `MATIS`.
3619 . -mat_is_allow_repeated - Allow repeated entries in the local part of the local to global maps.
3620 . -mat_is_fixempty - Fix local matrices in case of empty local rows/columns.
3621 - -mat_is_storel2l - Store the local-to-local operators generated by the Galerkin process of `MatPtAP()`.
3622
3623 Level: intermediate
3624
3625 Notes:
3626 Options prefix for the inner matrix are given by `-is_mat_xxx`
3627
3628 You must call `MatSetLocalToGlobalMapping()` before using this matrix type.
3629
3630 You can do matrix preallocation on the local matrix after you obtain it with
3631 `MatISGetLocalMat()`; otherwise, you could use `MatISSetPreallocation()` or `MatXAIJSetPreallocation()`
3632
3633 .seealso: [](ch_matrices), `Mat`, `MATIS`, `Mat`, `MatISGetLocalMat()`, `MatSetLocalToGlobalMapping()`, `MatISSetPreallocation()`, `MatCreateIS()`, `PCBDDC`, `KSPFETIDP`
3634 M*/
MatCreate_IS(Mat A)3635 PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
3636 {
3637 Mat_IS *a;
3638
3639 PetscFunctionBegin;
3640 PetscCall(PetscNew(&a));
3641 PetscCall(PetscStrallocpy(MATAIJ, &a->lmattype));
3642 A->data = (void *)a;
3643
3644 /* matrix ops */
3645 PetscCall(PetscMemzero(A->ops, sizeof(struct _MatOps)));
3646 A->ops->mult = MatMult_IS;
3647 A->ops->multadd = MatMultAdd_IS;
3648 A->ops->multtranspose = MatMultTranspose_IS;
3649 A->ops->multtransposeadd = MatMultTransposeAdd_IS;
3650 A->ops->destroy = MatDestroy_IS;
3651 A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
3652 A->ops->setvalues = MatSetValues_IS;
3653 A->ops->setvaluesblocked = MatSetValuesBlocked_IS;
3654 A->ops->setvalueslocal = MatSetValuesLocal_IS;
3655 A->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_IS;
3656 A->ops->zerorows = MatZeroRows_IS;
3657 A->ops->zerorowscolumns = MatZeroRowsColumns_IS;
3658 A->ops->assemblybegin = MatAssemblyBegin_IS;
3659 A->ops->assemblyend = MatAssemblyEnd_IS;
3660 A->ops->view = MatView_IS;
3661 A->ops->load = MatLoad_IS;
3662 A->ops->zeroentries = MatZeroEntries_IS;
3663 A->ops->scale = MatScale_IS;
3664 A->ops->getdiagonal = MatGetDiagonal_IS;
3665 A->ops->setoption = MatSetOption_IS;
3666 A->ops->ishermitian = MatIsHermitian_IS;
3667 A->ops->issymmetric = MatIsSymmetric_IS;
3668 A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS;
3669 A->ops->duplicate = MatDuplicate_IS;
3670 A->ops->copy = MatCopy_IS;
3671 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_IS;
3672 A->ops->createsubmatrix = MatCreateSubMatrix_IS;
3673 A->ops->axpy = MatAXPY_IS;
3674 A->ops->diagonalset = MatDiagonalSet_IS;
3675 A->ops->shift = MatShift_IS;
3676 A->ops->transpose = MatTranspose_IS;
3677 A->ops->getinfo = MatGetInfo_IS;
3678 A->ops->diagonalscale = MatDiagonalScale_IS;
3679 A->ops->setfromoptions = MatSetFromOptions_IS;
3680 A->ops->setup = MatSetUp_IS;
3681 A->ops->hasoperation = MatHasOperation_IS;
3682 A->ops->getdiagonalblock = MatGetDiagonalBlock_IS;
3683 A->ops->createsubmatrices = MatCreateSubMatrices_IS;
3684 A->ops->increaseoverlap = MatIncreaseOverlap_IS;
3685 A->ops->setblocksizes = MatSetBlockSizes_IS;
3686
3687 /* special MATIS functions */
3688 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMatType_C", MatISSetLocalMatType_IS));
3689 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalMat_C", MatISGetLocalMat_IS));
3690 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISRestoreLocalMat_C", MatISRestoreLocalMat_IS));
3691 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMat_C", MatISSetLocalMat_IS));
3692 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetPreallocation_C", MatISSetPreallocation_IS));
3693 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetAllowRepeated_C", MatISSetAllowRepeated_IS));
3694 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISStoreL2L_C", MatISStoreL2L_IS));
3695 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISFixLocalEmpty_C", MatISFixLocalEmpty_IS));
3696 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalToGlobalMapping_C", MatISGetLocalToGlobalMapping_IS));
3697 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpiaij_C", MatConvert_IS_XAIJ));
3698 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpibaij_C", MatConvert_IS_XAIJ));
3699 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpisbaij_C", MatConvert_IS_XAIJ));
3700 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqaij_C", MatConvert_IS_XAIJ));
3701 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqbaij_C", MatConvert_IS_XAIJ));
3702 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqsbaij_C", MatConvert_IS_XAIJ));
3703 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_aij_C", MatConvert_IS_XAIJ));
3704 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOOLocal_C", MatSetPreallocationCOOLocal_IS));
3705 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", MatSetPreallocationCOO_IS));
3706 PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATIS));
3707 PetscFunctionReturn(PETSC_SUCCESS);
3708 }
3709