1 /*
2 This file defines a "solve the problem redistributely on each subgroup of processor" preconditioner.
3 */
4 #include <petsc/private/pcimpl.h> /*I "petscksp.h" I*/
5 #include <petscksp.h>
6
7 typedef struct _PC_FieldSplitLink *PC_FieldSplitLink;
8 struct _PC_FieldSplitLink {
9 char *splitname;
10 IS is;
11 PC_FieldSplitLink next, previous;
12 };
13
14 typedef struct {
15 KSP ksp;
16 Vec x, b;
17 VecScatter scatter;
18 IS is;
19 PetscInt dcnt, *drows; /* these are the local rows that have only diagonal entry */
20 PetscScalar *diag;
21 Vec work;
22 PetscBool zerodiag;
23
24 PetscInt nsplits;
25 PC_FieldSplitLink splitlinks;
26 } PC_Redistribute;
27
PCFieldSplitSetIS_Redistribute(PC pc,const char splitname[],IS is)28 static PetscErrorCode PCFieldSplitSetIS_Redistribute(PC pc, const char splitname[], IS is)
29 {
30 PC_Redistribute *red = (PC_Redistribute *)pc->data;
31 PC_FieldSplitLink *next = &red->splitlinks;
32
33 PetscFunctionBegin;
34 while (*next) next = &(*next)->next;
35 PetscCall(PetscNew(next));
36 if (splitname) {
37 PetscCall(PetscStrallocpy(splitname, &(*next)->splitname));
38 } else {
39 PetscCall(PetscMalloc1(8, &(*next)->splitname));
40 PetscCall(PetscSNPrintf((*next)->splitname, 7, "%" PetscInt_FMT, red->nsplits++));
41 }
42 PetscCall(PetscObjectReference((PetscObject)is));
43 PetscCall(ISDestroy(&(*next)->is));
44 (*next)->is = is;
45 PetscFunctionReturn(PETSC_SUCCESS);
46 }
47
PCView_Redistribute(PC pc,PetscViewer viewer)48 static PetscErrorCode PCView_Redistribute(PC pc, PetscViewer viewer)
49 {
50 PC_Redistribute *red = (PC_Redistribute *)pc->data;
51 PetscBool isascii, isstring;
52 PetscInt ncnt, N;
53
54 PetscFunctionBegin;
55 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
56 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
57 if (isascii) {
58 PetscCallMPI(MPIU_Allreduce(&red->dcnt, &ncnt, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)pc)));
59 PetscCall(MatGetSize(pc->pmat, &N, NULL));
60 PetscCall(PetscViewerASCIIPrintf(viewer, " Number rows eliminated %" PetscInt_FMT " Percentage rows eliminated %g\n", ncnt, (double)(100 * ((PetscReal)ncnt) / ((PetscReal)N))));
61 PetscCall(PetscViewerASCIIPrintf(viewer, " Redistribute preconditioner: \n"));
62 PetscCall(KSPView(red->ksp, viewer));
63 } else if (isstring) {
64 PetscCall(PetscViewerStringSPrintf(viewer, " Redistribute preconditioner"));
65 PetscCall(KSPView(red->ksp, viewer));
66 }
67 PetscFunctionReturn(PETSC_SUCCESS);
68 }
69
PCSetUp_Redistribute(PC pc)70 static PetscErrorCode PCSetUp_Redistribute(PC pc)
71 {
72 PC_Redistribute *red = (PC_Redistribute *)pc->data;
73 MPI_Comm comm;
74 PetscInt rstart, rend, nrstart, nrend, nz, cnt, *rows, ncnt, dcnt, *drows;
75 PetscLayout map, nmap;
76 PetscMPIInt size, tag, n;
77 PETSC_UNUSED PetscMPIInt imdex;
78 PetscInt *source = NULL;
79 PetscMPIInt *sizes = NULL, nrecvs, nsends;
80 PetscInt j;
81 PetscInt *owner = NULL, *starts = NULL, count, slen;
82 PetscInt *rvalues, *svalues, recvtotal;
83 PetscMPIInt *onodes1, *olengths1;
84 MPI_Request *send_waits = NULL, *recv_waits = NULL;
85 MPI_Status recv_status, *send_status;
86 Vec tvec, diag;
87 Mat tmat;
88 const PetscScalar *d, *values;
89 const PetscInt *cols;
90 PC_FieldSplitLink *next = &red->splitlinks;
91
92 PetscFunctionBegin;
93 if (pc->setupcalled) {
94 PetscCheck(pc->flag == SAME_NONZERO_PATTERN, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PC is not supported for a change in the nonzero structure of the matrix");
95 PetscCall(KSPGetOperators(red->ksp, NULL, &tmat));
96 PetscCall(MatCreateSubMatrix(pc->pmat, red->is, red->is, MAT_REUSE_MATRIX, &tmat));
97 PetscCall(KSPSetOperators(red->ksp, tmat, tmat));
98 } else {
99 PetscInt NN;
100 PC ipc;
101 PetscBool fptr;
102
103 PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
104 PetscCallMPI(MPI_Comm_size(comm, &size));
105 PetscCall(PetscObjectGetNewTag((PetscObject)pc, &tag));
106
107 /* count non-diagonal rows on process */
108 PetscCall(MatGetOwnershipRange(pc->mat, &rstart, &rend));
109 cnt = 0;
110 for (PetscInt i = rstart; i < rend; i++) {
111 PetscCall(MatGetRow(pc->mat, i, &nz, &cols, &values));
112 for (PetscInt j = 0; j < nz; j++) {
113 if (values[j] != 0 && cols[j] != i) {
114 cnt++;
115 break;
116 }
117 }
118 PetscCall(MatRestoreRow(pc->mat, i, &nz, &cols, &values));
119 }
120 PetscCall(PetscMalloc1(cnt, &rows));
121 PetscCall(PetscMalloc1(rend - rstart - cnt, &drows));
122
123 /* list non-diagonal rows on process */
124 cnt = 0;
125 dcnt = 0;
126 for (PetscInt i = rstart; i < rend; i++) {
127 PetscBool diagonly = PETSC_TRUE;
128 PetscCall(MatGetRow(pc->mat, i, &nz, &cols, &values));
129 for (PetscInt j = 0; j < nz; j++) {
130 if (values[j] != 0 && cols[j] != i) {
131 diagonly = PETSC_FALSE;
132 break;
133 }
134 }
135 if (!diagonly) rows[cnt++] = i;
136 else drows[dcnt++] = i - rstart;
137 PetscCall(MatRestoreRow(pc->mat, i, &nz, &cols, &values));
138 }
139
140 /* create PetscLayout for non-diagonal rows on each process */
141 PetscCall(PetscLayoutCreate(comm, &map));
142 PetscCall(PetscLayoutSetLocalSize(map, cnt));
143 PetscCall(PetscLayoutSetBlockSize(map, 1));
144 PetscCall(PetscLayoutSetUp(map));
145 nrstart = map->rstart;
146 nrend = map->rend;
147
148 /* create PetscLayout for load-balanced non-diagonal rows on each process */
149 PetscCall(PetscLayoutCreate(comm, &nmap));
150 PetscCallMPI(MPIU_Allreduce(&cnt, &ncnt, 1, MPIU_INT, MPI_SUM, comm));
151 PetscCall(PetscLayoutSetSize(nmap, ncnt));
152 PetscCall(PetscLayoutSetBlockSize(nmap, 1));
153 PetscCall(PetscLayoutSetUp(nmap));
154
155 PetscCall(MatGetSize(pc->pmat, &NN, NULL));
156 PetscCall(PetscInfo(pc, "Number of diagonal rows eliminated %" PetscInt_FMT ", percentage eliminated %g\n", NN - ncnt, (double)((PetscReal)(NN - ncnt) / (PetscReal)NN)));
157
158 if (size > 1) {
159 /*
160 the following block of code assumes MPI can send messages to self, which is not supported for MPI-uni hence we need to handle
161 the size 1 case as a special case
162
163 this code is taken from VecScatterCreate_PtoS()
164 Determines what rows need to be moved where to
165 load balance the non-diagonal rows
166 */
167 /* count number of contributors to each processor */
168 PetscCall(PetscMalloc2(size, &sizes, cnt, &owner));
169 PetscCall(PetscArrayzero(sizes, size));
170 j = 0;
171 nsends = 0;
172 for (PetscInt i = nrstart; i < nrend; i++) {
173 if (i < nmap->range[j]) j = 0;
174 for (; j < size; j++) {
175 if (i < nmap->range[j + 1]) {
176 if (!sizes[j]++) nsends++;
177 owner[i - nrstart] = j;
178 break;
179 }
180 }
181 }
182 /* inform other processors of number of messages and max length*/
183 PetscCall(PetscGatherNumberOfMessages(comm, NULL, sizes, &nrecvs));
184 PetscCall(PetscGatherMessageLengths(comm, nsends, nrecvs, sizes, &onodes1, &olengths1));
185 PetscCall(PetscSortMPIIntWithArray(nrecvs, onodes1, olengths1));
186 recvtotal = 0;
187 for (PetscMPIInt i = 0; i < nrecvs; i++) recvtotal += olengths1[i];
188
189 /* post receives: rvalues - rows I will own; count - nu */
190 PetscCall(PetscMalloc3(recvtotal, &rvalues, nrecvs, &source, nrecvs, &recv_waits));
191 count = 0;
192 for (PetscMPIInt i = 0; i < nrecvs; i++) {
193 PetscCallMPI(MPIU_Irecv(rvalues + count, olengths1[i], MPIU_INT, onodes1[i], tag, comm, recv_waits + i));
194 count += olengths1[i];
195 }
196
197 /* do sends:
198 1) starts[i] gives the starting index in svalues for stuff going to
199 the ith processor
200 */
201 PetscCall(PetscMalloc3(cnt, &svalues, nsends, &send_waits, size, &starts));
202 starts[0] = 0;
203 for (PetscMPIInt i = 1; i < size; i++) starts[i] = starts[i - 1] + sizes[i - 1];
204 for (PetscInt i = 0; i < cnt; i++) svalues[starts[owner[i]]++] = rows[i];
205 for (PetscInt i = 0; i < cnt; i++) rows[i] = rows[i] - nrstart;
206 red->drows = drows;
207 red->dcnt = dcnt;
208 PetscCall(PetscFree(rows));
209
210 starts[0] = 0;
211 for (PetscMPIInt i = 1; i < size; i++) starts[i] = starts[i - 1] + sizes[i - 1];
212 count = 0;
213 for (PetscMPIInt i = 0; i < size; i++) {
214 if (sizes[i]) PetscCallMPI(MPIU_Isend(svalues + starts[i], sizes[i], MPIU_INT, i, tag, comm, send_waits + count++));
215 }
216
217 /* wait on receives */
218 count = nrecvs;
219 slen = 0;
220 while (count) {
221 PetscCallMPI(MPI_Waitany(nrecvs, recv_waits, &imdex, &recv_status));
222 /* unpack receives into our local space */
223 PetscCallMPI(MPI_Get_count(&recv_status, MPIU_INT, &n));
224 slen += n;
225 count--;
226 }
227 PetscCheck(slen == recvtotal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total message lengths %" PetscInt_FMT " not expected %" PetscInt_FMT, slen, recvtotal);
228 PetscCall(ISCreateGeneral(comm, slen, rvalues, PETSC_COPY_VALUES, &red->is));
229
230 /* free all work space */
231 PetscCall(PetscFree(olengths1));
232 PetscCall(PetscFree(onodes1));
233 PetscCall(PetscFree3(rvalues, source, recv_waits));
234 PetscCall(PetscFree2(sizes, owner));
235 if (nsends) { /* wait on sends */
236 PetscCall(PetscMalloc1(nsends, &send_status));
237 PetscCallMPI(MPI_Waitall(nsends, send_waits, send_status));
238 PetscCall(PetscFree(send_status));
239 }
240 PetscCall(PetscFree3(svalues, send_waits, starts));
241 } else {
242 PetscCall(ISCreateGeneral(comm, cnt, rows, PETSC_OWN_POINTER, &red->is));
243 red->drows = drows;
244 red->dcnt = dcnt;
245 slen = cnt;
246 }
247 PetscCall(PetscLayoutDestroy(&map));
248
249 PetscCall(VecCreateMPI(comm, slen, PETSC_DETERMINE, &red->b));
250 PetscCall(VecDuplicate(red->b, &red->x));
251 PetscCall(MatCreateVecs(pc->pmat, &tvec, NULL));
252 PetscCall(VecScatterCreate(tvec, red->is, red->b, NULL, &red->scatter));
253
254 /* Map the PCFIELDSPLIT fields to redistributed KSP */
255 PetscCall(KSPGetPC(red->ksp, &ipc));
256 PetscCall(PetscObjectHasFunction((PetscObject)ipc, "PCFieldSplitSetIS_C", &fptr));
257 if (fptr && *next) {
258 PetscScalar *atvec;
259 const PetscScalar *ab;
260 PetscInt primes[] = {2, 3, 5, 7, 11, 13, 17, 19};
261 PetscInt cnt = 0;
262
263 PetscCheck(red->nsplits <= (PetscInt)PETSC_STATIC_ARRAY_LENGTH(primes), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No support for this many fields");
264 PetscCall(VecSet(tvec, 1.0));
265 PetscCall(VecGetArray(tvec, &atvec));
266
267 while (*next) {
268 const PetscInt *indices;
269 PetscInt n;
270
271 PetscCall(ISGetIndices((*next)->is, &indices));
272 PetscCall(ISGetLocalSize((*next)->is, &n));
273 for (PetscInt i = 0; i < n; i++) atvec[indices[i] - rstart] *= primes[cnt];
274 PetscCall(ISRestoreIndices((*next)->is, &indices));
275 cnt++;
276 next = &(*next)->next;
277 }
278 PetscCall(VecRestoreArray(tvec, &atvec));
279 PetscCall(VecScatterBegin(red->scatter, tvec, red->b, INSERT_VALUES, SCATTER_FORWARD));
280 PetscCall(VecScatterEnd(red->scatter, tvec, red->b, INSERT_VALUES, SCATTER_FORWARD));
281 cnt = 0;
282 PetscCall(VecGetArrayRead(red->b, &ab));
283 next = &red->splitlinks;
284 while (*next) {
285 PetscInt n = 0;
286 PetscInt *indices;
287 IS ris;
288
289 for (PetscInt i = 0; i < nmap->rend - nmap->rstart; i++) {
290 if (!(((PetscInt)PetscRealPart(ab[i])) % primes[cnt])) n++;
291 }
292 PetscCall(PetscMalloc1(n, &indices));
293 n = 0;
294 for (PetscInt i = 0; i < nmap->rend - nmap->rstart; i++) {
295 if (!(((PetscInt)PetscRealPart(ab[i])) % primes[cnt])) indices[n++] = i + nmap->rstart;
296 }
297 PetscCall(ISCreateGeneral(comm, n, indices, PETSC_OWN_POINTER, &ris));
298 PetscCall(PCFieldSplitSetIS(ipc, (*next)->splitname, ris));
299
300 PetscCall(ISDestroy(&ris));
301 cnt++;
302 next = &(*next)->next;
303 }
304 PetscCall(VecRestoreArrayRead(red->b, &ab));
305 }
306 PetscCall(VecDestroy(&tvec));
307 PetscCall(MatCreateSubMatrix(pc->pmat, red->is, red->is, MAT_INITIAL_MATRIX, &tmat));
308 PetscCall(KSPSetOperators(red->ksp, tmat, tmat));
309 PetscCall(MatDestroy(&tmat));
310 PetscCall(PetscLayoutDestroy(&nmap));
311 }
312
313 /* get diagonal portion of matrix */
314 PetscCall(PetscFree(red->diag));
315 PetscCall(PetscMalloc1(red->dcnt, &red->diag));
316 PetscCall(MatCreateVecs(pc->pmat, &diag, NULL));
317 PetscCall(MatGetDiagonal(pc->pmat, diag));
318 PetscCall(VecGetArrayRead(diag, &d));
319 for (PetscInt i = 0; i < red->dcnt; i++) {
320 if (d[red->drows[i]] != 0) red->diag[i] = 1.0 / d[red->drows[i]];
321 else {
322 red->zerodiag = PETSC_TRUE;
323 red->diag[i] = 0.0;
324 }
325 }
326 PetscCall(VecRestoreArrayRead(diag, &d));
327 PetscCall(VecDestroy(&diag));
328 PetscCall(KSPSetUp(red->ksp));
329 PetscFunctionReturn(PETSC_SUCCESS);
330 }
331
PCApply_Redistribute(PC pc,Vec b,Vec x)332 static PetscErrorCode PCApply_Redistribute(PC pc, Vec b, Vec x)
333 {
334 PC_Redistribute *red = (PC_Redistribute *)pc->data;
335 PetscInt dcnt = red->dcnt, i;
336 const PetscInt *drows = red->drows;
337 PetscScalar *xwork;
338 const PetscScalar *bwork, *diag = red->diag;
339 PetscBool nonzero_guess;
340
341 PetscFunctionBegin;
342 if (!red->work) PetscCall(VecDuplicate(b, &red->work));
343 PetscCall(KSPGetInitialGuessNonzero(red->ksp, &nonzero_guess));
344 if (nonzero_guess) {
345 PetscCall(VecScatterBegin(red->scatter, x, red->x, INSERT_VALUES, SCATTER_FORWARD));
346 PetscCall(VecScatterEnd(red->scatter, x, red->x, INSERT_VALUES, SCATTER_FORWARD));
347 }
348
349 /* compute the rows of solution that have diagonal entries only */
350 PetscCall(VecSet(x, 0.0)); /* x = diag(A)^{-1} b */
351 PetscCall(VecGetArray(x, &xwork));
352 PetscCall(VecGetArrayRead(b, &bwork));
353 if (red->zerodiag) {
354 for (i = 0; i < dcnt; i++) {
355 if (diag[i] == 0.0 && bwork[drows[i]] != 0.0) {
356 PetscCheck(!pc->erroriffailure, PETSC_COMM_SELF, PETSC_ERR_CONV_FAILED, "Linear system is inconsistent, zero matrix row but nonzero right-hand side");
357 PetscCall(PetscInfo(pc, "Linear system is inconsistent, zero matrix row but nonzero right-hand side\n"));
358 pc->failedreasonrank = PC_INCONSISTENT_RHS;
359 }
360 }
361 PetscCall(VecFlag(x, pc->failedreasonrank == PC_INCONSISTENT_RHS));
362 }
363 for (i = 0; i < dcnt; i++) xwork[drows[i]] = diag[i] * bwork[drows[i]];
364 PetscCall(PetscLogFlops(dcnt));
365 PetscCall(VecRestoreArray(red->work, &xwork));
366 PetscCall(VecRestoreArrayRead(b, &bwork));
367 /* update the right-hand side for the reduced system with diagonal rows (and corresponding columns) removed */
368 PetscCall(MatMult(pc->pmat, x, red->work));
369 PetscCall(VecAYPX(red->work, -1.0, b)); /* red->work = b - A x */
370
371 PetscCall(VecScatterBegin(red->scatter, red->work, red->b, INSERT_VALUES, SCATTER_FORWARD));
372 PetscCall(VecScatterEnd(red->scatter, red->work, red->b, INSERT_VALUES, SCATTER_FORWARD));
373 PetscCall(KSPSolve(red->ksp, red->b, red->x));
374 PetscCall(KSPCheckSolve(red->ksp, pc, red->x));
375 PetscCall(VecScatterBegin(red->scatter, red->x, x, INSERT_VALUES, SCATTER_REVERSE));
376 PetscCall(VecScatterEnd(red->scatter, red->x, x, INSERT_VALUES, SCATTER_REVERSE));
377 PetscFunctionReturn(PETSC_SUCCESS);
378 }
379
PCApplyTranspose_Redistribute(PC pc,Vec b,Vec x)380 static PetscErrorCode PCApplyTranspose_Redistribute(PC pc, Vec b, Vec x)
381 {
382 PC_Redistribute *red = (PC_Redistribute *)pc->data;
383 PetscInt dcnt = red->dcnt, i;
384 const PetscInt *drows = red->drows;
385 PetscScalar *xwork;
386 const PetscScalar *bwork, *diag = red->diag;
387 PetscBool set, flg = PETSC_FALSE, nonzero_guess;
388
389 PetscFunctionBegin;
390 PetscCall(MatIsStructurallySymmetricKnown(pc->pmat, &set, &flg));
391 PetscCheck(set || flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PCApplyTranspose() not implemented for structurally unsymmetric Mat");
392 if (!red->work) PetscCall(VecDuplicate(b, &red->work));
393 PetscCall(KSPGetInitialGuessNonzero(red->ksp, &nonzero_guess));
394 if (nonzero_guess) {
395 PetscCall(VecScatterBegin(red->scatter, x, red->x, INSERT_VALUES, SCATTER_FORWARD));
396 PetscCall(VecScatterEnd(red->scatter, x, red->x, INSERT_VALUES, SCATTER_FORWARD));
397 }
398
399 /* compute the rows of solution that have diagonal entries only */
400 PetscCall(VecSet(x, 0.0)); /* x = diag(A)^{-1} b */
401 PetscCall(VecGetArray(x, &xwork));
402 PetscCall(VecGetArrayRead(b, &bwork));
403 if (red->zerodiag) {
404 for (i = 0; i < dcnt; i++) {
405 if (diag[i] == 0.0 && bwork[drows[i]] != 0.0) {
406 PetscCheck(!pc->erroriffailure, PETSC_COMM_SELF, PETSC_ERR_CONV_FAILED, "Linear system is inconsistent, zero matrix row but nonzero right-hand side");
407 PetscCall(PetscInfo(pc, "Linear system is inconsistent, zero matrix row but nonzero right-hand side\n"));
408 pc->failedreasonrank = PC_INCONSISTENT_RHS;
409 }
410 }
411 PetscCall(VecFlag(x, pc->failedreasonrank == PC_INCONSISTENT_RHS));
412 }
413 for (i = 0; i < dcnt; i++) xwork[drows[i]] = diag[i] * bwork[drows[i]];
414 PetscCall(PetscLogFlops(dcnt));
415 PetscCall(VecRestoreArray(red->work, &xwork));
416 PetscCall(VecRestoreArrayRead(b, &bwork));
417 /* update the right-hand side for the reduced system with diagonal rows (and corresponding columns) removed */
418 PetscCall(MatMultTranspose(pc->pmat, x, red->work));
419 PetscCall(VecAYPX(red->work, -1.0, b)); /* red->work = b - A^T x */
420
421 PetscCall(VecScatterBegin(red->scatter, red->work, red->b, INSERT_VALUES, SCATTER_FORWARD));
422 PetscCall(VecScatterEnd(red->scatter, red->work, red->b, INSERT_VALUES, SCATTER_FORWARD));
423 PetscCall(KSPSolveTranspose(red->ksp, red->b, red->x));
424 PetscCall(KSPCheckSolve(red->ksp, pc, red->x));
425 PetscCall(VecScatterBegin(red->scatter, red->x, x, INSERT_VALUES, SCATTER_REVERSE));
426 PetscCall(VecScatterEnd(red->scatter, red->x, x, INSERT_VALUES, SCATTER_REVERSE));
427 PetscFunctionReturn(PETSC_SUCCESS);
428 }
429
PCDestroy_Redistribute(PC pc)430 static PetscErrorCode PCDestroy_Redistribute(PC pc)
431 {
432 PC_Redistribute *red = (PC_Redistribute *)pc->data;
433 PC_FieldSplitLink next = red->splitlinks;
434
435 PetscFunctionBegin;
436 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", NULL));
437
438 while (next) {
439 PC_FieldSplitLink ilink;
440 PetscCall(PetscFree(next->splitname));
441 PetscCall(ISDestroy(&next->is));
442 ilink = next;
443 next = next->next;
444 PetscCall(PetscFree(ilink));
445 }
446 PetscCall(VecScatterDestroy(&red->scatter));
447 PetscCall(ISDestroy(&red->is));
448 PetscCall(VecDestroy(&red->b));
449 PetscCall(VecDestroy(&red->x));
450 PetscCall(KSPDestroy(&red->ksp));
451 PetscCall(VecDestroy(&red->work));
452 PetscCall(PetscFree(red->drows));
453 PetscCall(PetscFree(red->diag));
454 PetscCall(PetscFree(pc->data));
455 PetscFunctionReturn(PETSC_SUCCESS);
456 }
457
PCSetFromOptions_Redistribute(PC pc,PetscOptionItems PetscOptionsObject)458 static PetscErrorCode PCSetFromOptions_Redistribute(PC pc, PetscOptionItems PetscOptionsObject)
459 {
460 PC_Redistribute *red = (PC_Redistribute *)pc->data;
461
462 PetscFunctionBegin;
463 PetscCall(KSPSetFromOptions(red->ksp));
464 PetscFunctionReturn(PETSC_SUCCESS);
465 }
466
467 /*@
468 PCRedistributeGetKSP - Gets the `KSP` created by the `PCREDISTRIBUTE`
469
470 Not Collective
471
472 Input Parameter:
473 . pc - the preconditioner context
474
475 Output Parameter:
476 . innerksp - the inner `KSP`
477
478 Level: advanced
479
480 .seealso: [](ch_ksp), `KSP`, `PCREDISTRIBUTE`
481 @*/
PCRedistributeGetKSP(PC pc,KSP * innerksp)482 PetscErrorCode PCRedistributeGetKSP(PC pc, KSP *innerksp)
483 {
484 PC_Redistribute *red = (PC_Redistribute *)pc->data;
485
486 PetscFunctionBegin;
487 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
488 PetscAssertPointer(innerksp, 2);
489 *innerksp = red->ksp;
490 PetscFunctionReturn(PETSC_SUCCESS);
491 }
492
493 /*MC
494 PCREDISTRIBUTE - Redistributes a matrix for load balancing, removing the rows (and the corresponding columns) that only have a diagonal entry and then
495 applies a `KSP` to that new smaller matrix
496
497 Level: intermediate
498
499 Notes:
500 Options for the redistribute `KSP` and `PC` with the options database prefix `-redistribute_`
501
502 Usually run this with `-ksp_type preonly`
503
504 If you have used `MatZeroRows()` to eliminate (for example, Dirichlet) boundary conditions for a symmetric problem then you can use, for example, `-ksp_type preonly
505 -pc_type redistribute -redistribute_ksp_type cg -redistribute_pc_type bjacobi -redistribute_sub_pc_type icc` to take advantage of the symmetry.
506
507 Supports the function `PCFieldSplitSetIS()`; pass the appropriate reduced field indices to an inner `PCFIELDSPLIT`, set with, for example
508 `-ksp_type preonly -pc_type redistribute -redistribute_pc_type fieldsplit`. Does not support the `PCFIELDSPLIT` options database keys.
509
510 This does NOT call a partitioner to reorder rows to lower communication; the ordering of the rows in the original matrix and redistributed matrix is the same. Rows are moved
511 between MPI processes inside the preconditioner to balance the number of rows on each process.
512
513 The matrix block information is lost with the possible removal of individual rows and columns of the matrix, thus the behavior of the preconditioner on the reduced
514 system may be very different (worse) than running that preconditioner on the full system. This is specifically true for elasticity problems.
515
516 Developer Note:
517 Should add an option to this preconditioner to use a partitioner to redistribute the rows to lower communication.
518
519 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PCRedistributeGetKSP()`, `MatZeroRows()`, `PCFieldSplitSetIS()`, `PCFIELDSPLIT`
520 M*/
521
PCCreate_Redistribute(PC pc)522 PETSC_EXTERN PetscErrorCode PCCreate_Redistribute(PC pc)
523 {
524 PC_Redistribute *red;
525 const char *prefix;
526
527 PetscFunctionBegin;
528 PetscCall(PetscNew(&red));
529 pc->data = (void *)red;
530
531 pc->ops->apply = PCApply_Redistribute;
532 pc->ops->applytranspose = PCApplyTranspose_Redistribute;
533 pc->ops->setup = PCSetUp_Redistribute;
534 pc->ops->destroy = PCDestroy_Redistribute;
535 pc->ops->setfromoptions = PCSetFromOptions_Redistribute;
536 pc->ops->view = PCView_Redistribute;
537
538 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &red->ksp));
539 PetscCall(KSPSetNestLevel(red->ksp, pc->kspnestlevel));
540 PetscCall(KSPSetErrorIfNotConverged(red->ksp, pc->erroriffailure));
541 PetscCall(PetscObjectIncrementTabLevel((PetscObject)red->ksp, (PetscObject)pc, 1));
542 PetscCall(PCGetOptionsPrefix(pc, &prefix));
543 PetscCall(KSPSetOptionsPrefix(red->ksp, prefix));
544 PetscCall(KSPAppendOptionsPrefix(red->ksp, "redistribute_"));
545 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", PCFieldSplitSetIS_Redistribute));
546 PetscFunctionReturn(PETSC_SUCCESS);
547 }
548