1 #include <../src/ksp/pc/impls/gamg/gamg.h> /*I "petscpc.h" I*/
2 #include <petscsf.h>
3
4 static PetscFunctionList PCGAMGClassicalProlongatorList = NULL;
5 static PetscBool PCGAMGClassicalPackageInitialized = PETSC_FALSE;
6
7 typedef struct {
8 PetscReal interp_threshold; /* interpolation threshold */
9 char prolongtype[256];
10 PetscInt nsmooths; /* number of jacobi smoothings on the prolongator */
11 } PC_GAMG_Classical;
12
13 /*@
14 PCGAMGClassicalSetType - Sets the type of classical interpolation to use with `PCGAMG`
15
16 Collective
17
18 Input Parameters:
19 + pc - the preconditioner context
20 - type - the interpolation to use, see `PCGAMGClassicalType()`
21
22 Options Database Key:
23 . -pc_gamg_classical_type <direct,standard> - set type of classical AMG prolongation
24
25 Level: intermediate
26
27 .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGClassicalType`, `PCGAMGClassicalGetType()`
28 @*/
PCGAMGClassicalSetType(PC pc,PCGAMGClassicalType type)29 PetscErrorCode PCGAMGClassicalSetType(PC pc, PCGAMGClassicalType type)
30 {
31 PetscFunctionBegin;
32 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
33 PetscTryMethod(pc, "PCGAMGClassicalSetType_C", (PC, PCGAMGClassicalType), (pc, type));
34 PetscFunctionReturn(PETSC_SUCCESS);
35 }
36
37 /*@
38 PCGAMGClassicalGetType - Gets the type of classical interpolation to use with `PCGAMG`
39
40 Collective
41
42 Input Parameter:
43 . pc - the preconditioner context
44
45 Output Parameter:
46 . type - the type used, see `PCGAMGClassicalType()`
47
48 Level: intermediate
49
50 .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGClassicalType`, `PCGAMGClassicalSetType()`
51 @*/
PCGAMGClassicalGetType(PC pc,PCGAMGClassicalType * type)52 PetscErrorCode PCGAMGClassicalGetType(PC pc, PCGAMGClassicalType *type)
53 {
54 PetscFunctionBegin;
55 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
56 PetscUseMethod(pc, "PCGAMGClassicalGetType_C", (PC, PCGAMGClassicalType *), (pc, type));
57 PetscFunctionReturn(PETSC_SUCCESS);
58 }
59
PCGAMGClassicalSetType_GAMG(PC pc,PCGAMGClassicalType type)60 static PetscErrorCode PCGAMGClassicalSetType_GAMG(PC pc, PCGAMGClassicalType type)
61 {
62 PC_MG *mg = (PC_MG *)pc->data;
63 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
64 PC_GAMG_Classical *cls = (PC_GAMG_Classical *)pc_gamg->subctx;
65
66 PetscFunctionBegin;
67 PetscCall(PetscStrncpy(cls->prolongtype, type, sizeof(cls->prolongtype)));
68 PetscFunctionReturn(PETSC_SUCCESS);
69 }
70
PCGAMGClassicalGetType_GAMG(PC pc,PCGAMGClassicalType * type)71 static PetscErrorCode PCGAMGClassicalGetType_GAMG(PC pc, PCGAMGClassicalType *type)
72 {
73 PC_MG *mg = (PC_MG *)pc->data;
74 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
75 PC_GAMG_Classical *cls = (PC_GAMG_Classical *)pc_gamg->subctx;
76
77 PetscFunctionBegin;
78 *type = cls->prolongtype;
79 PetscFunctionReturn(PETSC_SUCCESS);
80 }
81
PCGAMGCreateGraph_Classical(PC pc,Mat A,Mat * G)82 static PetscErrorCode PCGAMGCreateGraph_Classical(PC pc, Mat A, Mat *G)
83 {
84 PetscInt s, f, n, idx, lidx, gidx;
85 PetscInt r, c, ncols;
86 const PetscInt *rcol;
87 const PetscScalar *rval;
88 PetscInt *gcol;
89 PetscScalar *gval;
90 PetscReal rmax;
91 PetscInt cmax = 0;
92 PC_MG *mg = (PC_MG *)pc->data;
93 PC_GAMG *gamg = (PC_GAMG *)mg->innerctx;
94 PetscInt *gsparse, *lsparse;
95 PetscScalar *Amax;
96 MatType mtype;
97
98 PetscFunctionBegin;
99 PetscCall(MatGetOwnershipRange(A, &s, &f));
100 n = f - s;
101 PetscCall(PetscMalloc3(n, &lsparse, n, &gsparse, n, &Amax));
102
103 for (r = 0; r < n; r++) {
104 lsparse[r] = 0;
105 gsparse[r] = 0;
106 }
107
108 for (r = s; r < f; r++) {
109 /* determine the maximum off-diagonal in each row */
110 rmax = 0.;
111 PetscCall(MatGetRow(A, r, &ncols, &rcol, &rval));
112 for (c = 0; c < ncols; c++) {
113 if (PetscRealPart(-rval[c]) > rmax && rcol[c] != r) rmax = PetscRealPart(-rval[c]);
114 }
115 Amax[r - s] = rmax;
116 if (ncols > cmax) cmax = ncols;
117 lidx = 0;
118 gidx = 0;
119 /* create the local and global sparsity patterns */
120 for (c = 0; c < ncols; c++) {
121 if (PetscRealPart(-rval[c]) > gamg->threshold[0] * PetscRealPart(Amax[r - s]) || rcol[c] == r) {
122 if (rcol[c] < f && rcol[c] >= s) {
123 lidx++;
124 } else {
125 gidx++;
126 }
127 }
128 }
129 PetscCall(MatRestoreRow(A, r, &ncols, &rcol, &rval));
130 lsparse[r - s] = lidx;
131 gsparse[r - s] = gidx;
132 }
133 PetscCall(PetscMalloc2(cmax, &gval, cmax, &gcol));
134
135 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), G));
136 PetscCall(MatGetType(A, &mtype));
137 PetscCall(MatSetType(*G, mtype));
138 PetscCall(MatSetSizes(*G, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
139 PetscCall(MatMPIAIJSetPreallocation(*G, 0, lsparse, 0, gsparse));
140 PetscCall(MatSeqAIJSetPreallocation(*G, 0, lsparse));
141 for (r = s; r < f; r++) {
142 PetscCall(MatGetRow(A, r, &ncols, &rcol, &rval));
143 idx = 0;
144 for (c = 0; c < ncols; c++) {
145 /* classical strength of connection */
146 if (PetscRealPart(-rval[c]) > gamg->threshold[0] * PetscRealPart(Amax[r - s]) || rcol[c] == r) {
147 gcol[idx] = rcol[c];
148 gval[idx] = rval[c];
149 idx++;
150 }
151 }
152 PetscCall(MatSetValues(*G, 1, &r, idx, gcol, gval, INSERT_VALUES));
153 PetscCall(MatRestoreRow(A, r, &ncols, &rcol, &rval));
154 }
155 PetscCall(MatAssemblyBegin(*G, MAT_FINAL_ASSEMBLY));
156 PetscCall(MatAssemblyEnd(*G, MAT_FINAL_ASSEMBLY));
157
158 PetscCall(PetscFree2(gval, gcol));
159 PetscCall(PetscFree3(lsparse, gsparse, Amax));
160 PetscFunctionReturn(PETSC_SUCCESS);
161 }
162
PCGAMGCoarsen_Classical(PC pc,Mat * G,PetscCoarsenData ** agg_lists)163 static PetscErrorCode PCGAMGCoarsen_Classical(PC pc, Mat *G, PetscCoarsenData **agg_lists)
164 {
165 MatCoarsen crs;
166 MPI_Comm fcomm = ((PetscObject)pc)->comm;
167 const char *prefix;
168
169 PetscFunctionBegin;
170 PetscCheck(G, fcomm, PETSC_ERR_ARG_WRONGSTATE, "Must set Graph in PC in PCGAMG before coarsening");
171
172 PetscCall(MatCoarsenCreate(fcomm, &crs));
173 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix));
174 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)crs, prefix));
175 PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)crs, "pc_gamg_"));
176 PetscCall(MatCoarsenSetFromOptions(crs));
177 PetscCall(MatCoarsenSetAdjacency(crs, *G));
178 PetscCall(MatCoarsenSetStrictAggs(crs, PETSC_TRUE));
179 PetscCall(MatCoarsenApply(crs));
180 PetscCall(MatCoarsenGetData(crs, agg_lists));
181 PetscCall(MatCoarsenDestroy(&crs));
182 PetscFunctionReturn(PETSC_SUCCESS);
183 }
184
PCGAMGProlongator_Classical_Direct(PC pc,Mat A,PetscCoarsenData * agg_lists,Mat * P)185 static PetscErrorCode PCGAMGProlongator_Classical_Direct(PC pc, Mat A, PetscCoarsenData *agg_lists, Mat *P)
186 {
187 PC_MG *mg = (PC_MG *)pc->data;
188 PC_GAMG *gamg = (PC_GAMG *)mg->innerctx;
189 PetscBool iscoarse, isMPIAIJ, isSEQAIJ;
190 PetscInt fn, cn, fs, fe, cs, ce, i, j, ncols, col, row_f, row_c, cmax = 0, idx, noff;
191 PetscInt *lcid, *gcid, *lsparse, *gsparse, *colmap, *pcols;
192 const PetscInt *rcol;
193 PetscReal *Amax_pos, *Amax_neg;
194 PetscScalar g_pos, g_neg, a_pos, a_neg, diag, invdiag, alpha, beta, pij;
195 PetscScalar *pvals;
196 const PetscScalar *rval;
197 Mat lA, gA = NULL;
198 MatType mtype;
199 Vec C, lvec;
200 PetscLayout clayout;
201 PetscSF sf;
202 Mat_MPIAIJ *mpiaij;
203
204 PetscFunctionBegin;
205 PetscCall(MatGetOwnershipRange(A, &fs, &fe));
206 fn = fe - fs;
207 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPIAIJ, &isMPIAIJ));
208 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQAIJ, &isSEQAIJ));
209 PetscCheck(isMPIAIJ || isSEQAIJ, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Classical AMG requires MPIAIJ matrix");
210 if (isMPIAIJ) {
211 mpiaij = (Mat_MPIAIJ *)A->data;
212 lA = mpiaij->A;
213 gA = mpiaij->B;
214 lvec = mpiaij->lvec;
215 PetscCall(VecGetSize(lvec, &noff));
216 colmap = mpiaij->garray;
217 PetscCall(MatGetLayouts(A, NULL, &clayout));
218 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)A), &sf));
219 PetscCall(PetscSFSetGraphLayout(sf, clayout, noff, NULL, PETSC_COPY_VALUES, colmap));
220 PetscCall(PetscMalloc1(noff, &gcid));
221 } else {
222 lA = A;
223 }
224 PetscCall(PetscMalloc5(fn, &lsparse, fn, &gsparse, fn, &lcid, fn, &Amax_pos, fn, &Amax_neg));
225
226 /* count the number of coarse unknowns */
227 cn = 0;
228 for (i = 0; i < fn; i++) {
229 /* filter out singletons */
230 PetscCall(PetscCDIsEmptyAt(agg_lists, i, &iscoarse));
231 lcid[i] = -1;
232 if (!iscoarse) cn++;
233 }
234
235 /* create the coarse vector */
236 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)A), cn, PETSC_DECIDE, &C));
237 PetscCall(VecGetOwnershipRange(C, &cs, &ce));
238
239 cn = 0;
240 for (i = 0; i < fn; i++) {
241 PetscCall(PetscCDIsEmptyAt(agg_lists, i, &iscoarse));
242 if (!iscoarse) {
243 lcid[i] = cs + cn;
244 cn++;
245 } else {
246 lcid[i] = -1;
247 }
248 }
249
250 if (gA) {
251 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, lcid, gcid, MPI_REPLACE));
252 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, lcid, gcid, MPI_REPLACE));
253 }
254
255 /* determine the largest off-diagonal entries in each row */
256 for (i = fs; i < fe; i++) {
257 Amax_pos[i - fs] = 0.;
258 Amax_neg[i - fs] = 0.;
259 PetscCall(MatGetRow(A, i, &ncols, &rcol, &rval));
260 for (j = 0; j < ncols; j++) {
261 if ((PetscRealPart(-rval[j]) > Amax_neg[i - fs]) && i != rcol[j]) Amax_neg[i - fs] = PetscAbsScalar(rval[j]);
262 if ((PetscRealPart(rval[j]) > Amax_pos[i - fs]) && i != rcol[j]) Amax_pos[i - fs] = PetscAbsScalar(rval[j]);
263 }
264 if (ncols > cmax) cmax = ncols;
265 PetscCall(MatRestoreRow(A, i, &ncols, &rcol, &rval));
266 }
267 PetscCall(PetscMalloc2(cmax, &pcols, cmax, &pvals));
268 PetscCall(VecDestroy(&C));
269
270 /* count the on and off processor sparsity patterns for the prolongator */
271 for (i = 0; i < fn; i++) {
272 /* on */
273 lsparse[i] = 0;
274 gsparse[i] = 0;
275 if (lcid[i] >= 0) {
276 lsparse[i] = 1;
277 gsparse[i] = 0;
278 } else {
279 PetscCall(MatGetRow(lA, i, &ncols, &rcol, &rval));
280 for (j = 0; j < ncols; j++) {
281 col = rcol[j];
282 if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0] * Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0] * Amax_neg[i])) lsparse[i] += 1;
283 }
284 PetscCall(MatRestoreRow(lA, i, &ncols, &rcol, &rval));
285 /* off */
286 if (gA) {
287 PetscCall(MatGetRow(gA, i, &ncols, &rcol, &rval));
288 for (j = 0; j < ncols; j++) {
289 col = rcol[j];
290 if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0] * Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0] * Amax_neg[i])) gsparse[i] += 1;
291 }
292 PetscCall(MatRestoreRow(gA, i, &ncols, &rcol, &rval));
293 }
294 }
295 }
296
297 /* preallocate and create the prolongator */
298 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), P));
299 PetscCall(MatGetType(A, &mtype));
300 PetscCall(MatSetType(*P, mtype));
301 PetscCall(MatSetSizes(*P, fn, cn, PETSC_DETERMINE, PETSC_DETERMINE));
302 PetscCall(MatMPIAIJSetPreallocation(*P, 0, lsparse, 0, gsparse));
303 PetscCall(MatSeqAIJSetPreallocation(*P, 0, lsparse));
304
305 /* loop over local fine nodes -- get the diagonal, the sum of positive and negative strong and weak weights, and set up the row */
306 for (i = 0; i < fn; i++) {
307 /* determine on or off */
308 row_f = i + fs;
309 row_c = lcid[i];
310 if (row_c >= 0) {
311 pij = 1.;
312 PetscCall(MatSetValues(*P, 1, &row_f, 1, &row_c, &pij, INSERT_VALUES));
313 } else {
314 g_pos = 0.;
315 g_neg = 0.;
316 a_pos = 0.;
317 a_neg = 0.;
318 diag = 0.;
319
320 /* local connections */
321 PetscCall(MatGetRow(lA, i, &ncols, &rcol, &rval));
322 for (j = 0; j < ncols; j++) {
323 col = rcol[j];
324 if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0] * Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0] * Amax_neg[i])) {
325 if (PetscRealPart(rval[j]) > 0.) {
326 g_pos += rval[j];
327 } else {
328 g_neg += rval[j];
329 }
330 }
331 if (col != i) {
332 if (PetscRealPart(rval[j]) > 0.) {
333 a_pos += rval[j];
334 } else {
335 a_neg += rval[j];
336 }
337 } else {
338 diag = rval[j];
339 }
340 }
341 PetscCall(MatRestoreRow(lA, i, &ncols, &rcol, &rval));
342
343 /* ghosted connections */
344 if (gA) {
345 PetscCall(MatGetRow(gA, i, &ncols, &rcol, &rval));
346 for (j = 0; j < ncols; j++) {
347 col = rcol[j];
348 if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0] * Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0] * Amax_neg[i])) {
349 if (PetscRealPart(rval[j]) > 0.) {
350 g_pos += rval[j];
351 } else {
352 g_neg += rval[j];
353 }
354 }
355 if (PetscRealPart(rval[j]) > 0.) {
356 a_pos += rval[j];
357 } else {
358 a_neg += rval[j];
359 }
360 }
361 PetscCall(MatRestoreRow(gA, i, &ncols, &rcol, &rval));
362 }
363
364 if (g_neg == 0.) {
365 alpha = 0.;
366 } else {
367 alpha = -a_neg / g_neg;
368 }
369
370 if (g_pos == 0.) {
371 diag += a_pos;
372 beta = 0.;
373 } else {
374 beta = -a_pos / g_pos;
375 }
376 if (diag == 0.) {
377 invdiag = 0.;
378 } else invdiag = 1. / diag;
379 /* on */
380 PetscCall(MatGetRow(lA, i, &ncols, &rcol, &rval));
381 idx = 0;
382 for (j = 0; j < ncols; j++) {
383 col = rcol[j];
384 if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0] * Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0] * Amax_neg[i])) {
385 row_f = i + fs;
386 row_c = lcid[col];
387 /* set the values for on-processor ones */
388 if (PetscRealPart(rval[j]) < 0.) {
389 pij = rval[j] * alpha * invdiag;
390 } else {
391 pij = rval[j] * beta * invdiag;
392 }
393 if (PetscAbsScalar(pij) != 0.) {
394 pvals[idx] = pij;
395 pcols[idx] = row_c;
396 idx++;
397 }
398 }
399 }
400 PetscCall(MatRestoreRow(lA, i, &ncols, &rcol, &rval));
401 /* off */
402 if (gA) {
403 PetscCall(MatGetRow(gA, i, &ncols, &rcol, &rval));
404 for (j = 0; j < ncols; j++) {
405 col = rcol[j];
406 if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0] * Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0] * Amax_neg[i])) {
407 row_f = i + fs;
408 row_c = gcid[col];
409 /* set the values for on-processor ones */
410 if (PetscRealPart(rval[j]) < 0.) {
411 pij = rval[j] * alpha * invdiag;
412 } else {
413 pij = rval[j] * beta * invdiag;
414 }
415 if (PetscAbsScalar(pij) != 0.) {
416 pvals[idx] = pij;
417 pcols[idx] = row_c;
418 idx++;
419 }
420 }
421 }
422 PetscCall(MatRestoreRow(gA, i, &ncols, &rcol, &rval));
423 }
424 PetscCall(MatSetValues(*P, 1, &row_f, idx, pcols, pvals, INSERT_VALUES));
425 }
426 }
427
428 PetscCall(MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY));
429 PetscCall(MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY));
430
431 PetscCall(PetscFree5(lsparse, gsparse, lcid, Amax_pos, Amax_neg));
432
433 PetscCall(PetscFree2(pcols, pvals));
434 if (gA) {
435 PetscCall(PetscSFDestroy(&sf));
436 PetscCall(PetscFree(gcid));
437 }
438 PetscFunctionReturn(PETSC_SUCCESS);
439 }
440
PCGAMGTruncateProlongator_Private(PC pc,Mat * P)441 static PetscErrorCode PCGAMGTruncateProlongator_Private(PC pc, Mat *P)
442 {
443 PetscInt j, i, ps, pf, pn, pcs, pcf, pcn, idx, cmax;
444 const PetscScalar *pval;
445 const PetscInt *pcol;
446 PetscScalar *pnval;
447 PetscInt *pncol;
448 PetscInt ncols;
449 Mat Pnew;
450 PetscInt *lsparse, *gsparse;
451 PetscReal pmax_pos, pmax_neg, ptot_pos, ptot_neg, pthresh_pos, pthresh_neg;
452 PC_MG *mg = (PC_MG *)pc->data;
453 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
454 PC_GAMG_Classical *cls = (PC_GAMG_Classical *)pc_gamg->subctx;
455 MatType mtype;
456
457 PetscFunctionBegin;
458 /* trim and rescale with reallocation */
459 PetscCall(MatGetOwnershipRange(*P, &ps, &pf));
460 PetscCall(MatGetOwnershipRangeColumn(*P, &pcs, &pcf));
461 pn = pf - ps;
462 pcn = pcf - pcs;
463 PetscCall(PetscMalloc2(pn, &lsparse, pn, &gsparse));
464 /* allocate */
465 cmax = 0;
466 for (i = ps; i < pf; i++) {
467 lsparse[i - ps] = 0;
468 gsparse[i - ps] = 0;
469 PetscCall(MatGetRow(*P, i, &ncols, &pcol, &pval));
470 if (ncols > cmax) cmax = ncols;
471 pmax_pos = 0.;
472 pmax_neg = 0.;
473 for (j = 0; j < ncols; j++) {
474 if (PetscRealPart(pval[j]) > pmax_pos) {
475 pmax_pos = PetscRealPart(pval[j]);
476 } else if (PetscRealPart(pval[j]) < pmax_neg) {
477 pmax_neg = PetscRealPart(pval[j]);
478 }
479 }
480 for (j = 0; j < ncols; j++) {
481 if (PetscRealPart(pval[j]) >= pmax_pos * cls->interp_threshold || PetscRealPart(pval[j]) <= pmax_neg * cls->interp_threshold) {
482 if (pcol[j] >= pcs && pcol[j] < pcf) {
483 lsparse[i - ps]++;
484 } else {
485 gsparse[i - ps]++;
486 }
487 }
488 }
489 PetscCall(MatRestoreRow(*P, i, &ncols, &pcol, &pval));
490 }
491
492 PetscCall(PetscMalloc2(cmax, &pnval, cmax, &pncol));
493
494 PetscCall(MatGetType(*P, &mtype));
495 PetscCall(MatCreate(PetscObjectComm((PetscObject)*P), &Pnew));
496 PetscCall(MatSetType(Pnew, mtype));
497 PetscCall(MatSetSizes(Pnew, pn, pcn, PETSC_DETERMINE, PETSC_DETERMINE));
498 PetscCall(MatSeqAIJSetPreallocation(Pnew, 0, lsparse));
499 PetscCall(MatMPIAIJSetPreallocation(Pnew, 0, lsparse, 0, gsparse));
500
501 for (i = ps; i < pf; i++) {
502 PetscCall(MatGetRow(*P, i, &ncols, &pcol, &pval));
503 pmax_pos = 0.;
504 pmax_neg = 0.;
505 for (j = 0; j < ncols; j++) {
506 if (PetscRealPart(pval[j]) > pmax_pos) {
507 pmax_pos = PetscRealPart(pval[j]);
508 } else if (PetscRealPart(pval[j]) < pmax_neg) {
509 pmax_neg = PetscRealPart(pval[j]);
510 }
511 }
512 pthresh_pos = 0.;
513 pthresh_neg = 0.;
514 ptot_pos = 0.;
515 ptot_neg = 0.;
516 for (j = 0; j < ncols; j++) {
517 if (PetscRealPart(pval[j]) >= cls->interp_threshold * pmax_pos) {
518 pthresh_pos += PetscRealPart(pval[j]);
519 } else if (PetscRealPart(pval[j]) <= cls->interp_threshold * pmax_neg) {
520 pthresh_neg += PetscRealPart(pval[j]);
521 }
522 if (PetscRealPart(pval[j]) > 0.) {
523 ptot_pos += PetscRealPart(pval[j]);
524 } else {
525 ptot_neg += PetscRealPart(pval[j]);
526 }
527 }
528 if (PetscAbsReal(pthresh_pos) > 0.) ptot_pos /= pthresh_pos;
529 if (PetscAbsReal(pthresh_neg) > 0.) ptot_neg /= pthresh_neg;
530 idx = 0;
531 for (j = 0; j < ncols; j++) {
532 if (PetscRealPart(pval[j]) >= pmax_pos * cls->interp_threshold) {
533 pnval[idx] = ptot_pos * pval[j];
534 pncol[idx] = pcol[j];
535 idx++;
536 } else if (PetscRealPart(pval[j]) <= pmax_neg * cls->interp_threshold) {
537 pnval[idx] = ptot_neg * pval[j];
538 pncol[idx] = pcol[j];
539 idx++;
540 }
541 }
542 PetscCall(MatRestoreRow(*P, i, &ncols, &pcol, &pval));
543 PetscCall(MatSetValues(Pnew, 1, &i, idx, pncol, pnval, INSERT_VALUES));
544 }
545
546 PetscCall(MatAssemblyBegin(Pnew, MAT_FINAL_ASSEMBLY));
547 PetscCall(MatAssemblyEnd(Pnew, MAT_FINAL_ASSEMBLY));
548 PetscCall(MatDestroy(P));
549
550 *P = Pnew;
551 PetscCall(PetscFree2(lsparse, gsparse));
552 PetscCall(PetscFree2(pnval, pncol));
553 PetscFunctionReturn(PETSC_SUCCESS);
554 }
555
PCGAMGProlongator_Classical_Standard(PC pc,Mat A,PetscCoarsenData * agg_lists,Mat * P)556 static PetscErrorCode PCGAMGProlongator_Classical_Standard(PC pc, Mat A, PetscCoarsenData *agg_lists, Mat *P)
557 {
558 Mat lA, *lAs;
559 MatType mtype;
560 Vec cv;
561 PetscInt *gcid, *lcid, *lsparse, *gsparse, *picol;
562 PetscInt fs, fe, cs, ce, nl, i, j, k, li, lni, ci, ncols, maxcols, fn, cn, cid;
563 PetscMPIInt size;
564 const PetscInt *lidx, *icol, *gidx;
565 PetscBool iscoarse;
566 PetscScalar vi, pentry, pjentry;
567 PetscScalar *pcontrib, *pvcol;
568 const PetscScalar *vcol;
569 PetscReal diag, jdiag, jwttotal;
570 PetscInt pncols;
571 PetscSF sf;
572 PetscLayout clayout;
573 IS lis;
574
575 PetscFunctionBegin;
576 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
577 PetscCall(MatGetOwnershipRange(A, &fs, &fe));
578 fn = fe - fs;
579 PetscCall(ISCreateStride(PETSC_COMM_SELF, fe - fs, fs, 1, &lis));
580 if (size > 1) {
581 PetscCall(MatGetLayouts(A, NULL, &clayout));
582 /* increase the overlap by two to get neighbors of neighbors */
583 PetscCall(MatIncreaseOverlap(A, 1, &lis, 2));
584 PetscCall(ISSort(lis));
585 /* get the local part of A */
586 PetscCall(MatCreateSubMatrices(A, 1, &lis, &lis, MAT_INITIAL_MATRIX, &lAs));
587 lA = lAs[0];
588 /* build an SF out of it */
589 PetscCall(ISGetLocalSize(lis, &nl));
590 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)A), &sf));
591 PetscCall(ISGetIndices(lis, &lidx));
592 PetscCall(PetscSFSetGraphLayout(sf, clayout, nl, NULL, PETSC_COPY_VALUES, lidx));
593 PetscCall(ISRestoreIndices(lis, &lidx));
594 } else {
595 lA = A;
596 nl = fn;
597 }
598 /* create a communication structure for the overlapped portion and transmit coarse indices */
599 PetscCall(PetscMalloc3(fn, &lsparse, fn, &gsparse, nl, &pcontrib));
600 /* create coarse vector */
601 cn = 0;
602 for (i = 0; i < fn; i++) {
603 PetscCall(PetscCDIsEmptyAt(agg_lists, i, &iscoarse));
604 if (!iscoarse) cn++;
605 }
606 PetscCall(PetscMalloc1(fn, &gcid));
607 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)A), cn, PETSC_DECIDE, &cv));
608 PetscCall(VecGetOwnershipRange(cv, &cs, &ce));
609 cn = 0;
610 for (i = 0; i < fn; i++) {
611 PetscCall(PetscCDIsEmptyAt(agg_lists, i, &iscoarse));
612 if (!iscoarse) {
613 gcid[i] = cs + cn;
614 cn++;
615 } else {
616 gcid[i] = -1;
617 }
618 }
619 if (size > 1) {
620 PetscCall(PetscMalloc1(nl, &lcid));
621 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, gcid, lcid, MPI_REPLACE));
622 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, gcid, lcid, MPI_REPLACE));
623 } else {
624 lcid = gcid;
625 }
626 /* count to preallocate the prolongator */
627 PetscCall(ISGetIndices(lis, &gidx));
628 maxcols = 0;
629 /* count the number of unique contributing coarse cells for each fine */
630 for (i = 0; i < nl; i++) {
631 pcontrib[i] = 0.;
632 PetscCall(MatGetRow(lA, i, &ncols, &icol, NULL));
633 if (gidx[i] >= fs && gidx[i] < fe) {
634 li = gidx[i] - fs;
635 lsparse[li] = 0;
636 gsparse[li] = 0;
637 cid = lcid[i];
638 if (cid >= 0) {
639 lsparse[li] = 1;
640 } else {
641 for (j = 0; j < ncols; j++) {
642 if (lcid[icol[j]] >= 0) {
643 pcontrib[icol[j]] = 1.;
644 } else {
645 ci = icol[j];
646 PetscCall(MatRestoreRow(lA, i, &ncols, &icol, NULL));
647 PetscCall(MatGetRow(lA, ci, &ncols, &icol, NULL));
648 for (k = 0; k < ncols; k++) {
649 if (lcid[icol[k]] >= 0) pcontrib[icol[k]] = 1.;
650 }
651 PetscCall(MatRestoreRow(lA, ci, &ncols, &icol, NULL));
652 PetscCall(MatGetRow(lA, i, &ncols, &icol, NULL));
653 }
654 }
655 for (j = 0; j < ncols; j++) {
656 if (lcid[icol[j]] >= 0 && pcontrib[icol[j]] != 0.) {
657 lni = lcid[icol[j]];
658 if (lni >= cs && lni < ce) {
659 lsparse[li]++;
660 } else {
661 gsparse[li]++;
662 }
663 pcontrib[icol[j]] = 0.;
664 } else {
665 ci = icol[j];
666 PetscCall(MatRestoreRow(lA, i, &ncols, &icol, NULL));
667 PetscCall(MatGetRow(lA, ci, &ncols, &icol, NULL));
668 for (k = 0; k < ncols; k++) {
669 if (lcid[icol[k]] >= 0 && pcontrib[icol[k]] != 0.) {
670 lni = lcid[icol[k]];
671 if (lni >= cs && lni < ce) {
672 lsparse[li]++;
673 } else {
674 gsparse[li]++;
675 }
676 pcontrib[icol[k]] = 0.;
677 }
678 }
679 PetscCall(MatRestoreRow(lA, ci, &ncols, &icol, NULL));
680 PetscCall(MatGetRow(lA, i, &ncols, &icol, NULL));
681 }
682 }
683 }
684 if (lsparse[li] + gsparse[li] > maxcols) maxcols = lsparse[li] + gsparse[li];
685 }
686 PetscCall(MatRestoreRow(lA, i, &ncols, &icol, &vcol));
687 }
688 PetscCall(PetscMalloc2(maxcols, &picol, maxcols, &pvcol));
689 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), P));
690 PetscCall(MatGetType(A, &mtype));
691 PetscCall(MatSetType(*P, mtype));
692 PetscCall(MatSetSizes(*P, fn, cn, PETSC_DETERMINE, PETSC_DETERMINE));
693 PetscCall(MatMPIAIJSetPreallocation(*P, 0, lsparse, 0, gsparse));
694 PetscCall(MatSeqAIJSetPreallocation(*P, 0, lsparse));
695 for (i = 0; i < nl; i++) {
696 diag = 0.;
697 if (gidx[i] >= fs && gidx[i] < fe) {
698 pncols = 0;
699 cid = lcid[i];
700 if (cid >= 0) {
701 pncols = 1;
702 picol[0] = cid;
703 pvcol[0] = 1.;
704 } else {
705 PetscCall(MatGetRow(lA, i, &ncols, &icol, &vcol));
706 for (j = 0; j < ncols; j++) {
707 pentry = vcol[j];
708 if (lcid[icol[j]] >= 0) {
709 /* coarse neighbor */
710 pcontrib[icol[j]] += pentry;
711 } else if (icol[j] != i) {
712 /* the neighbor is a strongly connected fine node */
713 ci = icol[j];
714 vi = vcol[j];
715 PetscCall(MatRestoreRow(lA, i, &ncols, &icol, &vcol));
716 PetscCall(MatGetRow(lA, ci, &ncols, &icol, &vcol));
717 jwttotal = 0.;
718 jdiag = 0.;
719 for (k = 0; k < ncols; k++) {
720 if (ci == icol[k]) jdiag = PetscRealPart(vcol[k]);
721 }
722 for (k = 0; k < ncols; k++) {
723 if (lcid[icol[k]] >= 0 && jdiag * PetscRealPart(vcol[k]) < 0.) {
724 pjentry = vcol[k];
725 jwttotal += PetscRealPart(pjentry);
726 }
727 }
728 if (jwttotal != 0.) {
729 jwttotal = PetscRealPart(vi) / jwttotal;
730 for (k = 0; k < ncols; k++) {
731 if (lcid[icol[k]] >= 0 && jdiag * PetscRealPart(vcol[k]) < 0.) {
732 pjentry = vcol[k] * jwttotal;
733 pcontrib[icol[k]] += pjentry;
734 }
735 }
736 } else {
737 diag += PetscRealPart(vi);
738 }
739 PetscCall(MatRestoreRow(lA, ci, &ncols, &icol, &vcol));
740 PetscCall(MatGetRow(lA, i, &ncols, &icol, &vcol));
741 } else {
742 diag += PetscRealPart(vcol[j]);
743 }
744 }
745 if (diag != 0.) {
746 diag = 1. / diag;
747 for (j = 0; j < ncols; j++) {
748 if (lcid[icol[j]] >= 0 && pcontrib[icol[j]] != 0.) {
749 /* the neighbor is a coarse node */
750 if (PetscAbsScalar(pcontrib[icol[j]]) > 0.0) {
751 lni = lcid[icol[j]];
752 pvcol[pncols] = -pcontrib[icol[j]] * diag;
753 picol[pncols] = lni;
754 pncols++;
755 }
756 pcontrib[icol[j]] = 0.;
757 } else {
758 /* the neighbor is a strongly connected fine node */
759 ci = icol[j];
760 PetscCall(MatRestoreRow(lA, i, &ncols, &icol, &vcol));
761 PetscCall(MatGetRow(lA, ci, &ncols, &icol, &vcol));
762 for (k = 0; k < ncols; k++) {
763 if (lcid[icol[k]] >= 0 && pcontrib[icol[k]] != 0.) {
764 if (PetscAbsScalar(pcontrib[icol[k]]) > 0.0) {
765 lni = lcid[icol[k]];
766 pvcol[pncols] = -pcontrib[icol[k]] * diag;
767 picol[pncols] = lni;
768 pncols++;
769 }
770 pcontrib[icol[k]] = 0.;
771 }
772 }
773 PetscCall(MatRestoreRow(lA, ci, &ncols, &icol, &vcol));
774 PetscCall(MatGetRow(lA, i, &ncols, &icol, &vcol));
775 }
776 pcontrib[icol[j]] = 0.;
777 }
778 PetscCall(MatRestoreRow(lA, i, &ncols, &icol, &vcol));
779 }
780 }
781 ci = gidx[i];
782 if (pncols > 0) PetscCall(MatSetValues(*P, 1, &ci, pncols, picol, pvcol, INSERT_VALUES));
783 }
784 }
785 PetscCall(ISRestoreIndices(lis, &gidx));
786 PetscCall(PetscFree2(picol, pvcol));
787 PetscCall(PetscFree3(lsparse, gsparse, pcontrib));
788 PetscCall(ISDestroy(&lis));
789 PetscCall(PetscFree(gcid));
790 if (size > 1) {
791 PetscCall(PetscFree(lcid));
792 PetscCall(MatDestroyMatrices(1, &lAs));
793 PetscCall(PetscSFDestroy(&sf));
794 }
795 PetscCall(VecDestroy(&cv));
796 PetscCall(MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY));
797 PetscCall(MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY));
798 PetscFunctionReturn(PETSC_SUCCESS);
799 }
800
PCGAMGOptProlongator_Classical_Jacobi(PC pc,Mat A,Mat * P)801 static PetscErrorCode PCGAMGOptProlongator_Classical_Jacobi(PC pc, Mat A, Mat *P)
802 {
803 PetscInt f, s, n, cf, cs, i, idx;
804 PetscInt *coarserows;
805 PetscInt ncols;
806 const PetscInt *pcols;
807 const PetscScalar *pvals;
808 Mat Pnew;
809 Vec diag;
810 PC_MG *mg = (PC_MG *)pc->data;
811 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
812 PC_GAMG_Classical *cls = (PC_GAMG_Classical *)pc_gamg->subctx;
813
814 PetscFunctionBegin;
815 if (cls->nsmooths == 0) {
816 PetscCall(PCGAMGTruncateProlongator_Private(pc, P));
817 PetscFunctionReturn(PETSC_SUCCESS);
818 }
819 PetscCall(MatGetOwnershipRange(*P, &s, &f));
820 n = f - s;
821 PetscCall(MatGetOwnershipRangeColumn(*P, &cs, &cf));
822 PetscCall(PetscMalloc1(n, &coarserows));
823 /* identify the rows corresponding to coarse unknowns */
824 idx = 0;
825 for (i = s; i < f; i++) {
826 PetscCall(MatGetRow(*P, i, &ncols, &pcols, &pvals));
827 /* assume, for now, that it's a coarse unknown if it has a single unit entry */
828 if (ncols == 1) {
829 if (pvals[0] == 1.) {
830 coarserows[idx] = i;
831 idx++;
832 }
833 }
834 PetscCall(MatRestoreRow(*P, i, &ncols, &pcols, &pvals));
835 }
836 PetscCall(MatCreateVecs(A, &diag, NULL));
837 PetscCall(MatGetDiagonal(A, diag));
838 PetscCall(VecReciprocal(diag));
839 for (i = 0; i < cls->nsmooths; i++) {
840 PetscCall(MatMatMult(A, *P, MAT_INITIAL_MATRIX, PETSC_CURRENT, &Pnew));
841 PetscCall(MatZeroRows(Pnew, idx, coarserows, 0., NULL, NULL));
842 PetscCall(MatDiagonalScale(Pnew, diag, NULL));
843 PetscCall(MatAYPX(Pnew, -1.0, *P, DIFFERENT_NONZERO_PATTERN));
844 PetscCall(MatDestroy(P));
845 *P = Pnew;
846 Pnew = NULL;
847 }
848 PetscCall(VecDestroy(&diag));
849 PetscCall(PetscFree(coarserows));
850 PetscCall(PCGAMGTruncateProlongator_Private(pc, P));
851 PetscFunctionReturn(PETSC_SUCCESS);
852 }
853
PCGAMGProlongator_Classical(PC pc,Mat A,PetscCoarsenData * agg_lists,Mat * P)854 static PetscErrorCode PCGAMGProlongator_Classical(PC pc, Mat A, PetscCoarsenData *agg_lists, Mat *P)
855 {
856 PetscErrorCode (*f)(PC, Mat, PetscCoarsenData *, Mat *);
857 PC_MG *mg = (PC_MG *)pc->data;
858 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
859 PC_GAMG_Classical *cls = (PC_GAMG_Classical *)pc_gamg->subctx;
860
861 PetscFunctionBegin;
862 PetscCall(PetscFunctionListFind(PCGAMGClassicalProlongatorList, cls->prolongtype, &f));
863 PetscCheck(f, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Cannot find PCGAMG Classical prolongator type");
864 PetscCall((*f)(pc, A, agg_lists, P));
865 PetscFunctionReturn(PETSC_SUCCESS);
866 }
867
PCGAMGDestroy_Classical(PC pc)868 static PetscErrorCode PCGAMGDestroy_Classical(PC pc)
869 {
870 PC_MG *mg = (PC_MG *)pc->data;
871 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
872
873 PetscFunctionBegin;
874 PetscCall(PetscFree(pc_gamg->subctx));
875 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGClassicalSetType_C", NULL));
876 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGClassicalGetType_C", NULL));
877 PetscFunctionReturn(PETSC_SUCCESS);
878 }
879
PCGAMGSetFromOptions_Classical(PC pc,PetscOptionItems PetscOptionsObject)880 static PetscErrorCode PCGAMGSetFromOptions_Classical(PC pc, PetscOptionItems PetscOptionsObject)
881 {
882 PC_MG *mg = (PC_MG *)pc->data;
883 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
884 PC_GAMG_Classical *cls = (PC_GAMG_Classical *)pc_gamg->subctx;
885 char tname[256];
886 PetscBool flg;
887
888 PetscFunctionBegin;
889 PetscOptionsHeadBegin(PetscOptionsObject, "GAMG-Classical options");
890 PetscCall(PetscOptionsFList("-pc_gamg_classical_type", "Type of Classical AMG prolongation", "PCGAMGClassicalSetType", PCGAMGClassicalProlongatorList, cls->prolongtype, tname, sizeof(tname), &flg));
891 if (flg) PetscCall(PCGAMGClassicalSetType(pc, tname));
892 PetscCall(PetscOptionsReal("-pc_gamg_classical_interp_threshold", "Threshold for classical interpolator entries", "", cls->interp_threshold, &cls->interp_threshold, NULL));
893 PetscCall(PetscOptionsInt("-pc_gamg_classical_nsmooths", "Threshold for classical interpolator entries", "", cls->nsmooths, &cls->nsmooths, NULL));
894 PetscOptionsHeadEnd();
895 PetscFunctionReturn(PETSC_SUCCESS);
896 }
897
PCGAMGSetData_Classical(PC pc,Mat A)898 static PetscErrorCode PCGAMGSetData_Classical(PC pc, Mat A)
899 {
900 PC_MG *mg = (PC_MG *)pc->data;
901 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
902
903 PetscFunctionBegin;
904 /* no data for classical AMG */
905 pc_gamg->data = NULL;
906 pc_gamg->data_cell_cols = 0;
907 pc_gamg->data_cell_rows = 0;
908 pc_gamg->data_sz = 0;
909 PetscFunctionReturn(PETSC_SUCCESS);
910 }
911
PCGAMGClassicalFinalizePackage(void)912 static PetscErrorCode PCGAMGClassicalFinalizePackage(void)
913 {
914 PetscFunctionBegin;
915 PCGAMGClassicalPackageInitialized = PETSC_FALSE;
916 PetscCall(PetscFunctionListDestroy(&PCGAMGClassicalProlongatorList));
917 PetscFunctionReturn(PETSC_SUCCESS);
918 }
919
PCGAMGClassicalInitializePackage(void)920 static PetscErrorCode PCGAMGClassicalInitializePackage(void)
921 {
922 PetscFunctionBegin;
923 if (PCGAMGClassicalPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
924 PetscCall(PetscFunctionListAdd(&PCGAMGClassicalProlongatorList, PCGAMGCLASSICALDIRECT, PCGAMGProlongator_Classical_Direct));
925 PetscCall(PetscFunctionListAdd(&PCGAMGClassicalProlongatorList, PCGAMGCLASSICALSTANDARD, PCGAMGProlongator_Classical_Standard));
926 PetscCall(PetscRegisterFinalize(PCGAMGClassicalFinalizePackage));
927 PetscFunctionReturn(PETSC_SUCCESS);
928 }
929
PCCreateGAMG_Classical(PC pc)930 PetscErrorCode PCCreateGAMG_Classical(PC pc)
931 {
932 PC_MG *mg = (PC_MG *)pc->data;
933 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
934 PC_GAMG_Classical *pc_gamg_classical;
935
936 PetscFunctionBegin;
937 PetscCall(PCGAMGClassicalInitializePackage());
938 if (pc_gamg->subctx) {
939 /* call base class */
940 PetscCall(PCDestroy_GAMG(pc));
941 }
942
943 /* create sub context for SA */
944 PetscCall(PetscNew(&pc_gamg_classical));
945 pc_gamg->subctx = pc_gamg_classical;
946 pc->ops->setfromoptions = PCGAMGSetFromOptions_Classical;
947 /* reset does not do anything; setup not virtual */
948
949 /* set internal function pointers */
950 pc_gamg->ops->destroy = PCGAMGDestroy_Classical;
951 pc_gamg->ops->creategraph = PCGAMGCreateGraph_Classical;
952 pc_gamg->ops->coarsen = PCGAMGCoarsen_Classical;
953 pc_gamg->ops->prolongator = PCGAMGProlongator_Classical;
954 pc_gamg->ops->optprolongator = PCGAMGOptProlongator_Classical_Jacobi;
955 pc_gamg->ops->setfromoptions = PCGAMGSetFromOptions_Classical;
956
957 pc_gamg->ops->createdefaultdata = PCGAMGSetData_Classical;
958 pc_gamg_classical->interp_threshold = 0.2;
959 pc_gamg_classical->nsmooths = 0;
960 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGClassicalSetType_C", PCGAMGClassicalSetType_GAMG));
961 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGClassicalGetType_C", PCGAMGClassicalGetType_GAMG));
962 PetscCall(PCGAMGClassicalSetType(pc, PCGAMGCLASSICALSTANDARD));
963 PetscFunctionReturn(PETSC_SUCCESS);
964 }
965