Actual source code: agg.c
1: /*
2: GAMG geometric-algebric multigrid PC - Mark Adams 2011
3: */
5: #include <../src/ksp/pc/impls/gamg/gamg.h>
6: #include <petscblaslapack.h>
7: #include <petscdm.h>
8: #include <petsc/private/kspimpl.h>
10: typedef struct {
11: PetscInt nsmooths; // number of smoothing steps to construct prolongation
12: PetscInt aggressive_coarsening_levels; // number of aggressive coarsening levels (square or MISk)
13: PetscInt aggressive_mis_k; // the k in MIS-k
14: PetscBool use_aggressive_square_graph;
15: PetscBool use_minimum_degree_ordering;
16: PetscBool use_low_mem_filter;
17: PetscBool graph_symmetrize;
18: MatCoarsen crs;
19: } PC_GAMG_AGG;
21: /*@
22: PCGAMGSetNSmooths - Set number of smoothing steps (1 is typical) used to construct the prolongation operator
24: Logically Collective
26: Input Parameters:
27: + pc - the preconditioner context
28: - n - the number of smooths, default is 1
30: Options Database Key:
31: . -pc_gamg_agg_nsmooths nsmooth - number of smoothing steps to use
33: Level: intermediate
35: Note:
36: This is a different concept from the number smoothing steps used during the linear solution process which
37: can be set with `-mg_levels_ksp_max_it`
39: Developer Note:
40: This should be named `PCGAMGAGGSetNSmooths()`.
42: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCMG`, `PCGAMG`
43: @*/
44: PetscErrorCode PCGAMGSetNSmooths(PC pc, PetscInt n)
45: {
46: PetscFunctionBegin;
49: PetscTryMethod(pc, "PCGAMGSetNSmooths_C", (PC, PetscInt), (pc, n));
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: static PetscErrorCode PCGAMGSetNSmooths_AGG(PC pc, PetscInt n)
54: {
55: PC_MG *mg = (PC_MG *)pc->data;
56: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
57: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
59: PetscFunctionBegin;
60: pc_gamg_agg->nsmooths = n;
61: PetscFunctionReturn(PETSC_SUCCESS);
62: }
64: /*@
65: PCGAMGSetAggressiveLevels - Use aggressive coarsening on first n levels
67: Logically Collective
69: Input Parameters:
70: + pc - the preconditioner context
71: - n - 0, 1 or more, the default is 1
73: Options Database Key:
74: . -pc_gamg_aggressive_coarsening n - the number of coarsenings to do aggressively
76: Level: intermediate
78: Note:
79: By default, aggressive coarsening squares the matrix (computes $A^T A$) before coarsening.
80: Calling `PCGAMGSetAggressiveSquareGraph()` with a value of `PETSC_FALSE` changes the aggressive coarsening strategy to use MIS-k, see `PCGAMGMISkSetAggressive()`.
82: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGMISkSetAggressive()`,
83: `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`, `PCGAMGSetLowMemoryFilter()`
84: @*/
85: PetscErrorCode PCGAMGSetAggressiveLevels(PC pc, PetscInt n)
86: {
87: PetscFunctionBegin;
90: PetscTryMethod(pc, "PCGAMGSetAggressiveLevels_C", (PC, PetscInt), (pc, n));
91: PetscFunctionReturn(PETSC_SUCCESS);
92: }
94: /*@
95: PCGAMGMISkSetAggressive - Number (k) distance in MIS coarsening (> 2 is aggressive)
97: Logically Collective
99: Input Parameters:
100: + pc - the preconditioner context
101: - n - 1 or more (default = 2)
103: Options Database Key:
104: . -pc_gamg_aggressive_mis_k n - the distance to use
106: Level: intermediate
108: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`,
109: `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`, `PCGAMGSetLowMemoryFilter()`
110: @*/
111: PetscErrorCode PCGAMGMISkSetAggressive(PC pc, PetscInt n)
112: {
113: PetscFunctionBegin;
116: PetscTryMethod(pc, "PCGAMGMISkSetAggressive_C", (PC, PetscInt), (pc, n));
117: PetscFunctionReturn(PETSC_SUCCESS);
118: }
120: /*@
121: PCGAMGSetAggressiveSquareGraph - Use graph square ($A^T A$) for aggressive coarsening. Coarsening is slower than the alternative (MIS-2), which is faster and uses less memory
123: Logically Collective
125: Input Parameters:
126: + pc - the preconditioner context
127: - b - default true
129: Options Database Key:
130: . -pc_gamg_aggressive_square_graph (true|false) - whether to use the graph square to aggressively coarsen
132: Level: intermediate
134: Notes:
135: If `b` is `PETSC_FALSE` then MIS-k is used for aggressive coarsening, see `PCGAMGMISkSetAggressive()`
137: Squaring the matrix to perform the aggressive coarsening is slower and requires more memory than using MIS-k, but may result in a better preconditioner
138: that converges faster.
140: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`, `PCGAMGMISkSetAggressive()`, `PCGAMGMISkSetMinDegreeOrdering()`, `PCGAMGSetLowMemoryFilter()`
141: @*/
142: PetscErrorCode PCGAMGSetAggressiveSquareGraph(PC pc, PetscBool b)
143: {
144: PetscFunctionBegin;
147: PetscTryMethod(pc, "PCGAMGSetAggressiveSquareGraph_C", (PC, PetscBool), (pc, b));
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: /*@
152: PCGAMGMISkSetMinDegreeOrdering - Use minimum degree ordering in greedy MIS algorithm
154: Logically Collective
156: Input Parameters:
157: + pc - the preconditioner context
158: - b - default false
160: Options Database Key:
161: . -pc_gamg_mis_k_minimum_degree_ordering (true|false) - use the minimum degree ordering
163: Level: intermediate
165: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`,
166: `PCGAMGSetAggressiveLevels()`, `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGSetLowMemoryFilter()`
167: @*/
168: PetscErrorCode PCGAMGMISkSetMinDegreeOrdering(PC pc, PetscBool b)
169: {
170: PetscFunctionBegin;
173: PetscTryMethod(pc, "PCGAMGMISkSetMinDegreeOrdering_C", (PC, PetscBool), (pc, b));
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: /*@
178: PCGAMGSetLowMemoryFilter - Use low memory graph/matrix filter
180: Logically Collective
182: Input Parameters:
183: + pc - the preconditioner context
184: - b - default false
186: Options Database Key:
187: . -pc_gamg_low_memory_threshold_filter (true|false) - use the low memory filter
189: Level: intermediate
191: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`,
192: `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`
193: @*/
194: PetscErrorCode PCGAMGSetLowMemoryFilter(PC pc, PetscBool b)
195: {
196: PetscFunctionBegin;
199: PetscTryMethod(pc, "PCGAMGSetLowMemoryFilter_C", (PC, PetscBool), (pc, b));
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: /*@
204: PCGAMGSetGraphSymmetrize - Symmetrize graph used for coarsening. Defaults to true, but if matrix has symmetric attribute, then not needed since the graph is already known to be symmetric
206: Logically Collective
208: Input Parameters:
209: + pc - the preconditioner context
210: - b - default true
212: Options Database Key:
213: . -pc_gamg_graph_symmetrize (true|false) - symmetrize the graph
215: Level: intermediate
217: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`, `MatCreateGraph()`,
218: `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`
219: @*/
220: PetscErrorCode PCGAMGSetGraphSymmetrize(PC pc, PetscBool b)
221: {
222: PetscFunctionBegin;
225: PetscTryMethod(pc, "PCGAMGSetGraphSymmetrize_C", (PC, PetscBool), (pc, b));
226: PetscFunctionReturn(PETSC_SUCCESS);
227: }
229: /*@
230: PCGAMGSetProlongatorFilter - Set threshold for filtering small entries from the prolongator (a kernel-preserving correction is applied afterward)
232: Logically Collective
234: Input Parameters:
235: + pc - the preconditioner context
236: - thr - threshold value; entries with absolute value below this are dropped (0 disables filtering)
238: Options Database Key:
239: . -pc_gamg_prolongator_filter thr - threshold for filtering small entries from prolongator (0=disabled, 0.0025=typical)
241: Level: intermediate
243: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGGetProlongatorFilter()`, `PCGAMGSetLowMemoryFilter()`
244: @*/
245: PetscErrorCode PCGAMGSetProlongatorFilter(PC pc, PetscReal thr)
246: {
247: PetscFunctionBegin;
250: PetscTryMethod(pc, "PCGAMGSetProlongatorFilter_C", (PC, PetscReal), (pc, thr));
251: PetscFunctionReturn(PETSC_SUCCESS);
252: }
254: /*@
255: PCGAMGGetProlongatorFilter - Get threshold for filtering small entries from the prolongator
257: Not Collective
259: Input Parameter:
260: . pc - the preconditioner context
262: Output Parameter:
263: . thr - threshold value; entries with absolute value below this are dropped (0 disables filtering)
265: Level: intermediate
267: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetProlongatorFilter()`, `PCGAMGSetLowMemoryFilter()`
268: @*/
269: PetscErrorCode PCGAMGGetProlongatorFilter(PC pc, PetscReal *thr)
270: {
271: PetscFunctionBegin;
273: PetscAssertPointer(thr, 2);
274: PetscUseMethod(pc, "PCGAMGGetProlongatorFilter_C", (PC, PetscReal *), (pc, thr));
275: PetscFunctionReturn(PETSC_SUCCESS);
276: }
278: static PetscErrorCode PCGAMGSetAggressiveLevels_AGG(PC pc, PetscInt n)
279: {
280: PC_MG *mg = (PC_MG *)pc->data;
281: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
282: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
284: PetscFunctionBegin;
285: pc_gamg_agg->aggressive_coarsening_levels = n;
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }
289: static PetscErrorCode PCGAMGMISkSetAggressive_AGG(PC pc, PetscInt n)
290: {
291: PC_MG *mg = (PC_MG *)pc->data;
292: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
293: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
295: PetscFunctionBegin;
296: pc_gamg_agg->aggressive_mis_k = n;
297: PetscFunctionReturn(PETSC_SUCCESS);
298: }
300: static PetscErrorCode PCGAMGSetAggressiveSquareGraph_AGG(PC pc, PetscBool b)
301: {
302: PC_MG *mg = (PC_MG *)pc->data;
303: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
304: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
306: PetscFunctionBegin;
307: pc_gamg_agg->use_aggressive_square_graph = b;
308: PetscFunctionReturn(PETSC_SUCCESS);
309: }
311: static PetscErrorCode PCGAMGSetLowMemoryFilter_AGG(PC pc, PetscBool b)
312: {
313: PC_MG *mg = (PC_MG *)pc->data;
314: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
315: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
317: PetscFunctionBegin;
318: pc_gamg_agg->use_low_mem_filter = b;
319: PetscFunctionReturn(PETSC_SUCCESS);
320: }
322: static PetscErrorCode PCGAMGSetGraphSymmetrize_AGG(PC pc, PetscBool b)
323: {
324: PC_MG *mg = (PC_MG *)pc->data;
325: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
326: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
328: PetscFunctionBegin;
329: pc_gamg_agg->graph_symmetrize = b;
330: PetscFunctionReturn(PETSC_SUCCESS);
331: }
333: static PetscErrorCode PCGAMGMISkSetMinDegreeOrdering_AGG(PC pc, PetscBool b)
334: {
335: PC_MG *mg = (PC_MG *)pc->data;
336: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
337: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
339: PetscFunctionBegin;
340: pc_gamg_agg->use_minimum_degree_ordering = b;
341: PetscFunctionReturn(PETSC_SUCCESS);
342: }
344: static PetscErrorCode PCGAMGSetProlongatorFilter_AGG(PC pc, PetscReal thr)
345: {
346: PC_MG *mg = (PC_MG *)pc->data;
347: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
349: PetscFunctionBegin;
350: PetscCheck(thr >= 0.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Filter threshold %g must be non-negative", (double)thr);
351: pc_gamg->prolongator_filter = thr;
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: static PetscErrorCode PCGAMGGetProlongatorFilter_AGG(PC pc, PetscReal *thr)
356: {
357: PC_MG *mg = (PC_MG *)pc->data;
358: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
360: PetscFunctionBegin;
361: *thr = pc_gamg->prolongator_filter;
362: PetscFunctionReturn(PETSC_SUCCESS);
363: }
365: static PetscErrorCode PCSetFromOptions_GAMG_AGG(PC pc, PetscOptionItems PetscOptionsObject)
366: {
367: PC_MG *mg = (PC_MG *)pc->data;
368: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
369: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
370: PetscBool n_aggressive_flg, old_sq_provided = PETSC_FALSE, new_sq_provided = PETSC_FALSE, new_sqr_graph = pc_gamg_agg->use_aggressive_square_graph;
371: PetscInt nsq_graph_old = 0;
372: PetscReal thr = pc_gamg->prolongator_filter;
373: PetscBool flg;
375: PetscFunctionBegin;
376: PetscOptionsHeadBegin(PetscOptionsObject, "GAMG-AGG options");
377: PetscCall(PetscOptionsInt("-pc_gamg_agg_nsmooths", "number of smoothing steps to construct prolongation, usually 1", "PCGAMGSetNSmooths", pc_gamg_agg->nsmooths, &pc_gamg_agg->nsmooths, NULL));
378: // aggressive coarsening logic with deprecated -pc_gamg_square_graph
379: PetscCall(PetscOptionsInt("-pc_gamg_aggressive_coarsening", "Number of aggressive coarsening (MIS-2) levels from finest", "PCGAMGSetAggressiveLevels", pc_gamg_agg->aggressive_coarsening_levels, &pc_gamg_agg->aggressive_coarsening_levels, &n_aggressive_flg));
380: if (!n_aggressive_flg)
381: PetscCall(PetscOptionsInt("-pc_gamg_square_graph", "Number of aggressive coarsening (MIS-2) levels from finest (deprecated alias for -pc_gamg_aggressive_coarsening)", "PCGAMGSetAggressiveLevels", nsq_graph_old, &nsq_graph_old, &old_sq_provided));
382: PetscCall(PetscOptionsBool("-pc_gamg_aggressive_square_graph", "Use square graph $(A^T A)$ for aggressive coarsening, if false, MIS-k (k=2) is used, see PCGAMGMISkSetAggressive()", "PCGAMGSetAggressiveSquareGraph", new_sqr_graph, &pc_gamg_agg->use_aggressive_square_graph, &new_sq_provided));
383: if (!new_sq_provided && old_sq_provided) {
384: pc_gamg_agg->aggressive_coarsening_levels = nsq_graph_old; // could be zero
385: pc_gamg_agg->use_aggressive_square_graph = PETSC_TRUE;
386: }
387: if (new_sq_provided && old_sq_provided)
388: PetscCall(PetscInfo(pc, "Warning: both -pc_gamg_square_graph and -pc_gamg_aggressive_coarsening are used. -pc_gamg_square_graph is deprecated, Number of aggressive levels is %" PetscInt_FMT "\n", pc_gamg_agg->aggressive_coarsening_levels));
389: PetscCall(PetscOptionsBool("-pc_gamg_mis_k_minimum_degree_ordering", "Use minimum degree ordering for greedy MIS", "PCGAMGMISkSetMinDegreeOrdering", pc_gamg_agg->use_minimum_degree_ordering, &pc_gamg_agg->use_minimum_degree_ordering, NULL));
390: PetscCall(PetscOptionsBool("-pc_gamg_low_memory_threshold_filter", "Use the (built-in) low memory graph/matrix filter", "PCGAMGSetLowMemoryFilter", pc_gamg_agg->use_low_mem_filter, &pc_gamg_agg->use_low_mem_filter, NULL));
391: PetscCall(PetscOptionsInt("-pc_gamg_aggressive_mis_k", "Number of levels of multigrid to use.", "PCGAMGMISkSetAggressive", pc_gamg_agg->aggressive_mis_k, &pc_gamg_agg->aggressive_mis_k, NULL));
392: PetscCall(PetscOptionsBool("-pc_gamg_graph_symmetrize", "Symmetrize graph for coarsening", "PCGAMGSetGraphSymmetrize", pc_gamg_agg->graph_symmetrize, &pc_gamg_agg->graph_symmetrize, NULL));
393: PetscCall(PetscOptionsReal("-pc_gamg_prolongator_filter", "Threshold for filtering small entries from prolongator (0=disabled)", "PCGAMGSetProlongatorFilter", thr, &thr, &flg));
394: if (flg) PetscCall(PCGAMGSetProlongatorFilter(pc, thr));
396: PetscOptionsHeadEnd();
397: PetscFunctionReturn(PETSC_SUCCESS);
398: }
400: static PetscErrorCode PCDestroy_GAMG_AGG(PC pc)
401: {
402: PC_MG *mg = (PC_MG *)pc->data;
403: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
404: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
406: PetscFunctionBegin;
407: PetscCall(MatCoarsenDestroy(&pc_gamg_agg->crs));
408: PetscCall(PetscFree(pc_gamg->subctx));
409: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", NULL));
410: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", NULL));
411: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetAggressive_C", NULL));
412: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetMinDegreeOrdering_C", NULL));
413: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetLowMemoryFilter_C", NULL));
414: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveSquareGraph_C", NULL));
415: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetGraphSymmetrize_C", NULL));
416: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetProlongatorFilter_C", NULL));
417: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGGetProlongatorFilter_C", NULL));
418: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
419: PetscFunctionReturn(PETSC_SUCCESS);
420: }
422: /*
423: PCSetCoordinates_AGG
425: Collective
427: Input Parameter:
428: . pc - the preconditioner context
429: . ndm - dimension of data (used for dof/vertex for Stokes)
430: . a_nloc - number of vertices local
431: . coords - [a_nloc][ndm] - interleaved coordinate data: {x_0, y_0, z_0, x_1, y_1, ...}
432: */
434: static PetscErrorCode PCSetCoordinates_AGG(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords)
435: {
436: PC_MG *mg = (PC_MG *)pc->data;
437: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
438: PetscInt arrsz, kk, ii, jj, nloc, ndatarows, ndf;
439: Mat mat = pc->pmat;
441: PetscFunctionBegin;
444: nloc = a_nloc;
446: /* SA: null space vectors */
447: PetscCall(MatGetBlockSize(mat, &ndf)); /* this does not work for Stokes */
448: if (coords && ndf == 1) pc_gamg->data_cell_cols = 1; /* scalar w/ coords and SA (not needed) */
449: else if (coords) {
450: PetscCheck(ndm <= ndf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "degrees of motion %" PetscInt_FMT " > block size %" PetscInt_FMT, ndm, ndf);
451: pc_gamg->data_cell_cols = (ndm == 2 ? 3 : 6); /* displacement elasticity */
452: if (ndm != ndf) PetscCheck(pc_gamg->data_cell_cols == ndf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Don't know how to create null space for ndm=%" PetscInt_FMT ", ndf=%" PetscInt_FMT ". Use MatSetNearNullSpace().", ndm, ndf);
453: } else pc_gamg->data_cell_cols = ndf; /* no data, force SA with constant null space vectors */
454: pc_gamg->data_cell_rows = ndatarows = ndf;
455: PetscCheck(pc_gamg->data_cell_cols > 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "pc_gamg->data_cell_cols %" PetscInt_FMT " <= 0", pc_gamg->data_cell_cols);
456: arrsz = nloc * pc_gamg->data_cell_rows * pc_gamg->data_cell_cols;
458: if (!pc_gamg->data || (pc_gamg->data_sz != arrsz)) {
459: PetscCall(PetscFree(pc_gamg->data));
460: PetscCall(PetscMalloc1(arrsz + 1, &pc_gamg->data));
461: }
462: /* copy data in - column-oriented */
463: for (kk = 0; kk < nloc; kk++) {
464: const PetscInt M = nloc * pc_gamg->data_cell_rows; /* stride into data */
465: PetscReal *data = &pc_gamg->data[kk * ndatarows]; /* start of cell */
467: if (pc_gamg->data_cell_cols == 1) *data = 1.0;
468: else {
469: /* translational modes */
470: for (ii = 0; ii < ndatarows; ii++) {
471: for (jj = 0; jj < ndatarows; jj++) {
472: if (ii == jj) data[ii * M + jj] = 1.0;
473: else data[ii * M + jj] = 0.0;
474: }
475: }
477: /* rotational modes */
478: if (coords) {
479: if (ndm == 2) {
480: data += 2 * M;
481: data[0] = -coords[2 * kk + 1];
482: data[1] = coords[2 * kk];
483: } else {
484: data += 3 * M;
485: data[0] = 0.0;
486: data[M + 0] = coords[3 * kk + 2];
487: data[2 * M + 0] = -coords[3 * kk + 1];
488: data[1] = -coords[3 * kk + 2];
489: data[M + 1] = 0.0;
490: data[2 * M + 1] = coords[3 * kk];
491: data[2] = coords[3 * kk + 1];
492: data[M + 2] = -coords[3 * kk];
493: data[2 * M + 2] = 0.0;
494: }
495: }
496: }
497: }
498: pc_gamg->data_sz = arrsz;
499: PetscFunctionReturn(PETSC_SUCCESS);
500: }
502: /*
503: PCSetData_AGG - called if data is not set with PCSetCoordinates.
504: Looks in Mat for near null space.
505: Does not work for Stokes
507: Input Parameter:
508: . pc -
509: . a_A - matrix to get (near) null space out of.
510: */
511: static PetscErrorCode PCSetData_AGG(PC pc, Mat a_A)
512: {
513: PC_MG *mg = (PC_MG *)pc->data;
514: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
515: MatNullSpace mnull;
517: PetscFunctionBegin;
518: PetscCall(MatGetNearNullSpace(a_A, &mnull));
519: if (!mnull) {
520: DM dm;
522: PetscCall(PCGetDM(pc, &dm));
523: if (!dm) PetscCall(MatGetDM(a_A, &dm));
524: if (dm) {
525: PetscObject deformation;
526: PetscInt Nf;
528: PetscCall(DMGetNumFields(dm, &Nf));
529: if (Nf) {
530: PetscCall(DMGetField(dm, 0, NULL, &deformation));
531: if (deformation) {
532: PetscCall(PetscObjectQuery(deformation, "nearnullspace", (PetscObject *)&mnull));
533: if (!mnull) PetscCall(PetscObjectQuery(deformation, "nullspace", (PetscObject *)&mnull));
534: }
535: }
536: }
537: }
539: if (!mnull) {
540: PetscInt bs, NN, MM;
542: PetscCall(MatGetBlockSize(a_A, &bs));
543: PetscCall(MatGetLocalSize(a_A, &MM, &NN));
544: PetscCheck(MM % bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MM %" PetscInt_FMT " must be divisible by bs %" PetscInt_FMT, MM, bs);
545: PetscCall(PCSetCoordinates_AGG(pc, bs, MM / bs, NULL));
546: } else {
547: PetscReal *nullvec;
548: PetscBool has_const;
549: PetscInt i, j, mlocal, nvec, bs;
550: const Vec *vecs;
551: const PetscScalar *v;
553: PetscCall(MatGetLocalSize(a_A, &mlocal, NULL));
554: PetscCall(MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs));
555: for (i = 0; i < nvec; i++) {
556: PetscCall(VecGetLocalSize(vecs[i], &j));
557: PetscCheck(j == mlocal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Attached null space vector size %" PetscInt_FMT " != matrix size %" PetscInt_FMT, j, mlocal);
558: }
559: pc_gamg->data_sz = (nvec + !!has_const) * mlocal;
560: PetscCall(PetscMalloc1((nvec + !!has_const) * mlocal, &nullvec));
561: if (has_const)
562: for (i = 0; i < mlocal; i++) nullvec[i] = 1.0;
563: for (i = 0; i < nvec; i++) {
564: PetscCall(VecGetArrayRead(vecs[i], &v));
565: for (j = 0; j < mlocal; j++) nullvec[(i + !!has_const) * mlocal + j] = PetscRealPart(v[j]);
566: PetscCall(VecRestoreArrayRead(vecs[i], &v));
567: }
568: pc_gamg->data = nullvec;
569: pc_gamg->data_cell_cols = (nvec + !!has_const);
570: PetscCall(MatGetBlockSize(a_A, &bs));
571: pc_gamg->data_cell_rows = bs;
572: }
573: PetscFunctionReturn(PETSC_SUCCESS);
574: }
576: /*
577: formProl0 - collect null space data for each aggregate, do QR, put R in coarse grid data and Q in P_0
579: Input Parameter:
580: . agg_llists - list of arrays with aggregates -- list from selected vertices of aggregate unselected vertices
581: . bs - row block size
582: . nSAvec - column bs of new P
583: . my0crs - global index of start of locals
584: . data_stride - bs*(nloc nodes + ghost nodes) [data_stride][nSAvec]
585: . data_in[data_stride*nSAvec] - local data on fine grid
586: . flid_fgid[data_stride/bs] - make local to global IDs, includes ghosts in 'locals_llist'
588: Output Parameter:
589: . a_data_out - in with fine grid data (w/ghosts), out with coarse grid data
590: . a_Prol - prolongation operator
591: */
592: static PetscErrorCode formProl0(PetscCoarsenData *agg_llists, PetscInt bs, PetscInt nSAvec, PetscInt my0crs, PetscInt data_stride, PetscReal data_in[], const PetscInt flid_fgid[], PetscReal **a_data_out, Mat a_Prol)
593: {
594: PetscInt Istart, my0, Iend, nloc, clid, flid = 0, aggID, kk, jj, ii, mm, nSelected, minsz, nghosts, out_data_stride;
595: MPI_Comm comm;
596: PetscReal *out_data;
597: PetscCDIntNd *pos;
598: PetscHMapI fgid_flid;
600: PetscFunctionBegin;
601: PetscCall(PetscObjectGetComm((PetscObject)a_Prol, &comm));
602: PetscCall(MatGetOwnershipRange(a_Prol, &Istart, &Iend));
603: nloc = (Iend - Istart) / bs;
604: my0 = Istart / bs;
605: PetscCheck((Iend - Istart) % bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Iend %" PetscInt_FMT " - Istart %" PetscInt_FMT " must be divisible by bs %" PetscInt_FMT, Iend, Istart, bs);
606: Iend /= bs;
607: nghosts = data_stride / bs - nloc;
609: PetscCall(PetscHMapICreateWithSize(2 * nghosts + 1, &fgid_flid));
611: for (kk = 0; kk < nghosts; kk++) PetscCall(PetscHMapISet(fgid_flid, flid_fgid[nloc + kk], nloc + kk));
613: /* count selected -- same as number of cols of P */
614: for (nSelected = mm = 0; mm < nloc; mm++) {
615: PetscBool ise;
617: PetscCall(PetscCDIsEmptyAt(agg_llists, mm, &ise));
618: if (!ise) nSelected++;
619: }
620: PetscCall(MatGetOwnershipRangeColumn(a_Prol, &ii, &jj));
621: PetscCheck((ii / nSAvec) == my0crs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ii %" PetscInt_FMT " /nSAvec %" PetscInt_FMT " != my0crs %" PetscInt_FMT, ii, nSAvec, my0crs);
622: PetscCheck(nSelected == (jj - ii) / nSAvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "nSelected %" PetscInt_FMT " != (jj %" PetscInt_FMT " - ii %" PetscInt_FMT ")/nSAvec %" PetscInt_FMT, nSelected, jj, ii, nSAvec);
624: /* aloc space for coarse point data (output) */
625: out_data_stride = nSelected * nSAvec;
627: PetscCall(PetscMalloc1(out_data_stride * nSAvec, &out_data));
628: for (ii = 0; ii < out_data_stride * nSAvec; ii++) out_data[ii] = PETSC_MAX_REAL;
629: *a_data_out = out_data; /* output - stride nSelected*nSAvec */
631: /* find points and set prolongation */
632: minsz = 100;
633: for (mm = clid = 0; mm < nloc; mm++) {
634: PetscCall(PetscCDCountAt(agg_llists, mm, &jj));
635: if (jj > 0) {
636: const PetscInt lid = mm, cgid = my0crs + clid;
637: PetscInt cids[100]; /* max bs */
638: PetscBLASInt asz, M, N, info;
639: PetscBLASInt Mdata, LDA, LWORK;
640: PetscScalar *qqc, *qqr, *TAU, *WORK;
641: PetscInt *fids;
642: PetscReal *data;
644: PetscCall(PetscBLASIntCast(jj, &asz));
645: PetscCall(PetscBLASIntCast(asz * bs, &M));
646: PetscCall(PetscBLASIntCast(nSAvec, &N));
647: PetscCall(PetscBLASIntCast(M + ((N - M > 0) ? N - M : 0), &Mdata));
648: PetscCall(PetscBLASIntCast(Mdata, &LDA));
649: PetscCall(PetscBLASIntCast(N * bs, &LWORK));
650: /* count agg */
651: if (asz < minsz) minsz = asz;
653: /* get block */
654: PetscCall(PetscMalloc5(Mdata * N, &qqc, M * N, &qqr, N, &TAU, LWORK, &WORK, M, &fids));
656: aggID = 0;
657: PetscCall(PetscCDGetHeadPos(agg_llists, lid, &pos));
658: while (pos) {
659: PetscInt gid1;
661: PetscCall(PetscCDIntNdGetID(pos, &gid1));
662: PetscCall(PetscCDGetNextPos(agg_llists, lid, &pos));
664: if (gid1 >= my0 && gid1 < Iend) flid = gid1 - my0;
665: else {
666: PetscCall(PetscHMapIGet(fgid_flid, gid1, &flid));
667: PetscCheck(flid >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cannot find gid1 in table");
668: }
669: /* copy in B_i matrix - column-oriented */
670: data = &data_in[flid * bs];
671: for (ii = 0; ii < bs; ii++) {
672: for (jj = 0; jj < N; jj++) {
673: PetscReal d = data[jj * data_stride + ii];
675: qqc[jj * Mdata + aggID * bs + ii] = d;
676: }
677: }
678: /* set fine IDs */
679: for (kk = 0; kk < bs; kk++) fids[aggID * bs + kk] = flid_fgid[flid] * bs + kk;
680: aggID++;
681: }
683: /* pad with zeros */
684: for (ii = asz * bs; ii < Mdata; ii++) {
685: for (jj = 0; jj < N; jj++, kk++) qqc[jj * Mdata + ii] = .0;
686: }
688: /* QR */
689: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
690: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&Mdata, &N, qqc, &LDA, TAU, WORK, &LWORK, &info));
691: PetscCall(PetscFPTrapPop());
692: PetscCheck(info == 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in xGEQRF LAPACK routine %" PetscBLASInt_FMT, info);
693: /* get R - column-oriented - output B_{i+1} */
694: {
695: PetscReal *data = &out_data[clid * nSAvec];
697: for (jj = 0; jj < nSAvec; jj++) {
698: for (ii = 0; ii < nSAvec; ii++) {
699: PetscCheck(data[jj * out_data_stride + ii] == PETSC_MAX_REAL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "data[jj*out_data_stride + ii] != %e", (double)PETSC_MAX_REAL);
700: if (ii <= jj) data[jj * out_data_stride + ii] = PetscRealPart(qqc[jj * Mdata + ii]);
701: else data[jj * out_data_stride + ii] = 0.;
702: }
703: }
704: }
706: /* get Q - row-oriented */
707: PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&Mdata, &N, &N, qqc, &LDA, TAU, WORK, &LWORK, &info));
708: PetscCheck(info == 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in ORGQR LAPACK routine argument %" PetscBLASInt_FMT, -info);
710: for (ii = 0; ii < M; ii++) {
711: for (jj = 0; jj < N; jj++) qqr[N * ii + jj] = qqc[jj * Mdata + ii];
712: }
714: /* add diagonal block of P0 */
715: for (kk = 0; kk < N; kk++) cids[kk] = N * cgid + kk; /* global col IDs in P0 */
716: PetscCall(MatSetValues(a_Prol, M, fids, N, cids, qqr, INSERT_VALUES));
717: PetscCall(PetscFree5(qqc, qqr, TAU, WORK, fids));
718: clid++;
719: } /* coarse agg */
720: } /* for all fine nodes */
721: PetscCall(MatAssemblyBegin(a_Prol, MAT_FINAL_ASSEMBLY));
722: PetscCall(MatAssemblyEnd(a_Prol, MAT_FINAL_ASSEMBLY));
723: PetscCall(PetscHMapIDestroy(&fgid_flid));
724: PetscFunctionReturn(PETSC_SUCCESS);
725: }
727: static PetscErrorCode PCView_GAMG_AGG(PC pc, PetscViewer viewer)
728: {
729: PC_MG *mg = (PC_MG *)pc->data;
730: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
731: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
733: PetscFunctionBegin;
734: PetscCall(PetscViewerASCIIPrintf(viewer, " AGG specific options\n"));
735: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of levels of aggressive coarsening %" PetscInt_FMT "\n", pc_gamg_agg->aggressive_coarsening_levels));
736: if (pc_gamg_agg->aggressive_coarsening_levels > 0) {
737: PetscCall(PetscViewerASCIIPrintf(viewer, " %s aggressive coarsening\n", !pc_gamg_agg->use_aggressive_square_graph ? "MIS-k" : "Square graph"));
738: if (!pc_gamg_agg->use_aggressive_square_graph) PetscCall(PetscViewerASCIIPrintf(viewer, " MIS-%" PetscInt_FMT " coarsening on aggressive levels\n", pc_gamg_agg->aggressive_mis_k));
739: }
740: PetscCall(PetscViewerASCIIPushTab(viewer));
741: PetscCall(PetscViewerASCIIPushTab(viewer));
742: PetscCall(PetscViewerASCIIPushTab(viewer));
743: PetscCall(PetscViewerASCIIPushTab(viewer));
744: if (pc_gamg_agg->crs) PetscCall(MatCoarsenView(pc_gamg_agg->crs, viewer));
745: else PetscCall(PetscViewerASCIIPrintf(viewer, "Coarsening algorithm not yet selected\n"));
746: PetscCall(PetscViewerASCIIPopTab(viewer));
747: PetscCall(PetscViewerASCIIPopTab(viewer));
748: PetscCall(PetscViewerASCIIPopTab(viewer));
749: PetscCall(PetscViewerASCIIPopTab(viewer));
750: PetscCall(PetscViewerASCIIPrintf(viewer, " Number smoothing steps to construct prolongation %" PetscInt_FMT "\n", pc_gamg_agg->nsmooths));
751: PetscFunctionReturn(PETSC_SUCCESS);
752: }
754: static PetscErrorCode PCGAMGCreateGraph_AGG(PC pc, Mat Amat, Mat *a_Gmat)
755: {
756: PC_MG *mg = (PC_MG *)pc->data;
757: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
758: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
759: const PetscReal vfilter = pc_gamg->threshold[pc_gamg->current_level];
760: PetscBool ishem, ismis;
761: const char *prefix;
762: MatInfo info0, info1;
763: PetscInt bs;
765: PetscFunctionBegin;
766: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
767: /* Note: depending on the algorithm that will be used for computing the coarse grid points this should pass PETSC_TRUE or PETSC_FALSE as the first argument */
768: /* MATCOARSENHEM requires numerical weights for edges so ensure they are computed */
769: PetscCall(MatCoarsenDestroy(&pc_gamg_agg->crs));
770: PetscCall(MatCoarsenCreate(PetscObjectComm((PetscObject)pc), &pc_gamg_agg->crs));
771: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix));
772: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc_gamg_agg->crs, prefix));
773: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)pc_gamg_agg->crs, "pc_gamg_"));
774: PetscCall(MatCoarsenSetFromOptions(pc_gamg_agg->crs));
775: PetscCall(MatGetBlockSize(Amat, &bs));
776: // check for valid indices wrt bs
777: for (int ii = 0; ii < pc_gamg_agg->crs->strength_index_size; ii++) {
778: PetscCheck(pc_gamg_agg->crs->strength_index[ii] >= 0 && pc_gamg_agg->crs->strength_index[ii] < bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Indices (%" PetscInt_FMT ") must be non-negative and < block size (%" PetscInt_FMT "), NB, can not use -mat_coarsen_strength_index with -mat_coarsen_strength_index",
779: pc_gamg_agg->crs->strength_index[ii], bs);
780: }
781: PetscCall(PetscObjectTypeCompare((PetscObject)pc_gamg_agg->crs, MATCOARSENHEM, &ishem));
782: if (ishem) {
783: if (pc_gamg_agg->aggressive_coarsening_levels) PetscCall(PetscInfo(pc, "HEM and aggressive coarsening ignored: HEM using %" PetscInt_FMT " iterations\n", pc_gamg_agg->crs->max_it));
784: pc_gamg_agg->aggressive_coarsening_levels = 0; // aggressive and HEM does not make sense
785: PetscCall(MatCoarsenSetMaximumIterations(pc_gamg_agg->crs, pc_gamg_agg->crs->max_it)); // for code coverage
786: PetscCall(MatCoarsenSetThreshold(pc_gamg_agg->crs, vfilter)); // for code coverage
787: } else {
788: PetscCall(PetscObjectTypeCompare((PetscObject)pc_gamg_agg->crs, MATCOARSENMIS, &ismis));
789: if (ismis && pc_gamg_agg->aggressive_coarsening_levels && !pc_gamg_agg->use_aggressive_square_graph) {
790: PetscCall(PetscInfo(pc, "MIS and aggressive coarsening and no square graph: force square graph\n"));
791: pc_gamg_agg->use_aggressive_square_graph = PETSC_TRUE;
792: }
793: }
794: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
795: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_GRAPH], 0, 0, 0, 0));
796: PetscCall(MatGetInfo(Amat, MAT_LOCAL, &info0)); /* global reduction */
798: if (ishem || pc_gamg_agg->use_low_mem_filter) {
799: PetscCall(MatCreateGraph(Amat, pc_gamg_agg->graph_symmetrize, (vfilter >= 0 || ishem) ? PETSC_TRUE : PETSC_FALSE, vfilter, pc_gamg_agg->crs->strength_index_size, pc_gamg_agg->crs->strength_index, a_Gmat));
800: } else {
801: // make scalar graph, symmetrize if not known to be symmetric, scale, but do not filter (expensive)
802: PetscCall(MatCreateGraph(Amat, pc_gamg_agg->graph_symmetrize, PETSC_TRUE, -1, pc_gamg_agg->crs->strength_index_size, pc_gamg_agg->crs->strength_index, a_Gmat));
803: if (vfilter >= 0) {
804: PetscInt Istart, Iend, ncols, nnz0, nnz1, NN, MM, nloc;
805: Mat tGmat, Gmat = *a_Gmat;
806: MPI_Comm comm;
807: const PetscScalar *vals;
808: const PetscInt *idx;
809: PetscInt *d_nnz, *o_nnz, kk, *garray = NULL, *AJ, maxcols = 0;
810: MatScalar *AA; // this is checked in graph
811: PetscBool isseqaij;
812: Mat a, b, c;
813: MatType jtype;
815: PetscCall(PetscObjectGetComm((PetscObject)Gmat, &comm));
816: PetscCall(PetscObjectBaseTypeCompare((PetscObject)Gmat, MATSEQAIJ, &isseqaij));
817: PetscCall(MatGetType(Gmat, &jtype));
818: PetscCall(MatCreate(comm, &tGmat));
819: PetscCall(MatSetType(tGmat, jtype));
821: /* TODO GPU: this can be called when filter = 0 -> Probably provide MatAIJThresholdCompress that compresses the entries below a threshold?
822: Also, if the matrix is symmetric, can we skip this
823: operation? It can be very expensive on large matrices. */
825: // global sizes
826: PetscCall(MatGetSize(Gmat, &MM, &NN));
827: PetscCall(MatGetOwnershipRange(Gmat, &Istart, &Iend));
828: nloc = Iend - Istart;
829: PetscCall(PetscMalloc2(nloc, &d_nnz, nloc, &o_nnz));
830: if (isseqaij) {
831: a = Gmat;
832: b = NULL;
833: } else {
834: Mat_MPIAIJ *d = (Mat_MPIAIJ *)Gmat->data;
836: a = d->A;
837: b = d->B;
838: garray = d->garray;
839: }
840: /* Determine upper bound on non-zeros needed in new filtered matrix */
841: for (PetscInt row = 0; row < nloc; row++) {
842: PetscCall(MatGetRow(a, row, &ncols, NULL, NULL));
843: d_nnz[row] = ncols;
844: if (ncols > maxcols) maxcols = ncols;
845: PetscCall(MatRestoreRow(a, row, &ncols, NULL, NULL));
846: }
847: if (b) {
848: for (PetscInt row = 0; row < nloc; row++) {
849: PetscCall(MatGetRow(b, row, &ncols, NULL, NULL));
850: o_nnz[row] = ncols;
851: if (ncols > maxcols) maxcols = ncols;
852: PetscCall(MatRestoreRow(b, row, &ncols, NULL, NULL));
853: }
854: }
855: PetscCall(MatSetSizes(tGmat, nloc, nloc, MM, MM));
856: PetscCall(MatSetBlockSizes(tGmat, 1, 1));
857: PetscCall(MatSeqAIJSetPreallocation(tGmat, 0, d_nnz));
858: PetscCall(MatMPIAIJSetPreallocation(tGmat, 0, d_nnz, 0, o_nnz));
859: PetscCall(MatSetOption(tGmat, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
860: PetscCall(PetscFree2(d_nnz, o_nnz));
861: PetscCall(PetscMalloc2(maxcols, &AA, maxcols, &AJ));
862: nnz0 = nnz1 = 0;
863: for (c = a, kk = 0; c && kk < 2; c = b, kk++) {
864: for (PetscInt row = 0, grow = Istart, ncol_row, jj; row < nloc; row++, grow++) {
865: PetscCall(MatGetRow(c, row, &ncols, &idx, &vals));
866: for (ncol_row = jj = 0; jj < ncols; jj++, nnz0++) {
867: PetscScalar sv = PetscAbs(PetscRealPart(vals[jj]));
868: if (PetscRealPart(sv) > vfilter) {
869: PetscInt cid = idx[jj] + Istart; //diag
871: nnz1++;
872: if (c != a) cid = garray[idx[jj]];
873: AA[ncol_row] = vals[jj];
874: AJ[ncol_row] = cid;
875: ncol_row++;
876: }
877: }
878: PetscCall(MatRestoreRow(c, row, &ncols, &idx, &vals));
879: PetscCall(MatSetValues(tGmat, 1, &grow, ncol_row, AJ, AA, INSERT_VALUES));
880: }
881: }
882: PetscCall(PetscFree2(AA, AJ));
883: PetscCall(MatAssemblyBegin(tGmat, MAT_FINAL_ASSEMBLY));
884: PetscCall(MatAssemblyEnd(tGmat, MAT_FINAL_ASSEMBLY));
885: PetscCall(MatPropagateSymmetryOptions(Gmat, tGmat)); /* Normal Mat options are not relevant ? */
886: PetscCall(PetscInfo(pc, "\t %g%% nnz after filtering, with threshold %g, %g nnz ave. (N=%" PetscInt_FMT ", max row size %" PetscInt_FMT "\n", (!nnz0) ? 1. : 100. * (double)nnz1 / (double)nnz0, (double)vfilter, (!nloc) ? 1. : (double)nnz0 / (double)nloc, MM, maxcols));
887: PetscCall(MatViewFromOptions(tGmat, NULL, "-mat_filter_graph_view"));
888: PetscCall(MatDestroy(&Gmat));
889: *a_Gmat = tGmat;
890: }
891: }
893: PetscCall(MatGetInfo(*a_Gmat, MAT_LOCAL, &info1)); /* global reduction */
894: if (info0.nz_used > 0) PetscCall(PetscInfo(pc, "Filtering left %g %% edges in graph (%e %e)\n", 100.0 * info1.nz_used * (double)(bs * bs) / info0.nz_used, info0.nz_used, info1.nz_used));
895: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_GRAPH], 0, 0, 0, 0));
896: PetscFunctionReturn(PETSC_SUCCESS);
897: }
899: typedef PetscInt NState;
900: static const NState NOT_DONE = -2;
901: static const NState DELETED = -1;
902: static const NState REMOVED = -3;
903: #define IS_SELECTED(s) (s != DELETED && s != NOT_DONE && s != REMOVED)
905: /*
906: fixAggregatesWithSquare - greedy grab of with G1 (unsquared graph) -- AIJ specific -- change to fixAggregatesWithSquare -- TODD
907: - AGG-MG specific: clears singletons out of 'selected_2'
909: Input Parameter:
910: . Gmat_2 - global matrix of squared graph (data not defined)
911: . Gmat_1 - base graph to grab with base graph
912: Input/Output Parameter:
913: . aggs_2 - linked list of aggs with gids)
914: */
915: static PetscErrorCode fixAggregatesWithSquare(PC pc, Mat Gmat_2, Mat Gmat_1, PetscCoarsenData *aggs_2)
916: {
917: PetscBool isMPI;
918: Mat_SeqAIJ *matA_1, *matB_1 = NULL;
919: MPI_Comm comm;
920: PetscInt lid, *ii, *idx, ix, Iend, my0, kk, n, j;
921: Mat_MPIAIJ *mpimat_2 = NULL, *mpimat_1 = NULL;
922: const PetscInt nloc = Gmat_2->rmap->n;
923: PetscScalar *cpcol_1_state, *cpcol_2_state, *cpcol_2_par_orig, *lid_parent_gid;
924: PetscInt *lid_cprowID_1 = NULL;
925: NState *lid_state;
926: Vec ghost_par_orig2;
927: PetscMPIInt rank;
929: PetscFunctionBegin;
930: PetscCall(PetscObjectGetComm((PetscObject)Gmat_2, &comm));
931: PetscCallMPI(MPI_Comm_rank(comm, &rank));
932: PetscCall(MatGetOwnershipRange(Gmat_1, &my0, &Iend));
934: /* get submatrices */
935: PetscCall(PetscStrbeginswith(((PetscObject)Gmat_1)->type_name, MATMPIAIJ, &isMPI));
936: PetscCall(PetscInfo(pc, "isMPI = %s\n", isMPI ? "yes" : "no"));
937: PetscCall(PetscMalloc3(nloc, &lid_state, nloc, &lid_parent_gid, nloc, &lid_cprowID_1));
938: for (lid = 0; lid < nloc; lid++) lid_cprowID_1[lid] = -1;
939: if (isMPI) {
940: /* grab matrix objects */
941: mpimat_2 = (Mat_MPIAIJ *)Gmat_2->data;
942: mpimat_1 = (Mat_MPIAIJ *)Gmat_1->data;
943: matA_1 = (Mat_SeqAIJ *)mpimat_1->A->data;
944: matB_1 = (Mat_SeqAIJ *)mpimat_1->B->data;
946: /* force compressed row storage for B matrix in AuxMat */
947: PetscCall(MatCheckCompressedRow(mpimat_1->B, matB_1->nonzerorowcnt, &matB_1->compressedrow, matB_1->i, Gmat_1->rmap->n, -1.0));
948: for (ix = 0; ix < matB_1->compressedrow.nrows; ix++) {
949: PetscInt lid = matB_1->compressedrow.rindex[ix];
951: PetscCheck(lid <= nloc && lid >= -1, PETSC_COMM_SELF, PETSC_ERR_USER, "lid %" PetscInt_FMT " out of range. nloc = %" PetscInt_FMT, lid, nloc);
952: if (lid != -1) lid_cprowID_1[lid] = ix;
953: }
954: } else {
955: PetscBool isAIJ;
957: PetscCall(PetscStrbeginswith(((PetscObject)Gmat_1)->type_name, MATSEQAIJ, &isAIJ));
958: PetscCheck(isAIJ, PETSC_COMM_SELF, PETSC_ERR_USER, "Require AIJ matrix.");
959: matA_1 = (Mat_SeqAIJ *)Gmat_1->data;
960: }
961: if (nloc > 0) PetscCheck(!matB_1 || matB_1->compressedrow.use, PETSC_COMM_SELF, PETSC_ERR_PLIB, "matB_1 && !matB_1->compressedrow.use: PETSc bug???");
962: /* get state of locals and selected gid for deleted */
963: for (lid = 0; lid < nloc; lid++) {
964: lid_parent_gid[lid] = -1.0;
965: lid_state[lid] = DELETED;
966: }
968: /* set lid_state */
969: for (lid = 0; lid < nloc; lid++) {
970: PetscCDIntNd *pos;
972: PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
973: if (pos) {
974: PetscInt gid1;
976: PetscCall(PetscCDIntNdGetID(pos, &gid1));
977: PetscCheck(gid1 == lid + my0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "gid1 %" PetscInt_FMT " != lid %" PetscInt_FMT " + my0 %" PetscInt_FMT, gid1, lid, my0);
978: lid_state[lid] = gid1;
979: }
980: }
982: /* map local to selected local, DELETED means a ghost owns it */
983: for (lid = 0; lid < nloc; lid++) {
984: NState state = lid_state[lid];
986: if (IS_SELECTED(state)) {
987: PetscCDIntNd *pos;
989: PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
990: while (pos) {
991: PetscInt gid1;
993: PetscCall(PetscCDIntNdGetID(pos, &gid1));
994: PetscCall(PetscCDGetNextPos(aggs_2, lid, &pos));
995: if (gid1 >= my0 && gid1 < Iend) lid_parent_gid[gid1 - my0] = (PetscScalar)(lid + my0);
996: }
997: }
998: }
999: /* get 'cpcol_1/2_state' & cpcol_2_par_orig - uses mpimat_1/2->lvec for temp space */
1000: if (isMPI) {
1001: Vec tempVec;
1003: /* get 'cpcol_1_state' */
1004: PetscCall(MatCreateVecs(Gmat_1, &tempVec, NULL));
1005: for (kk = 0, j = my0; kk < nloc; kk++, j++) {
1006: PetscScalar v = (PetscScalar)lid_state[kk];
1008: PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
1009: }
1010: PetscCall(VecAssemblyBegin(tempVec));
1011: PetscCall(VecAssemblyEnd(tempVec));
1012: PetscCall(VecScatterBegin(mpimat_1->Mvctx, tempVec, mpimat_1->lvec, INSERT_VALUES, SCATTER_FORWARD));
1013: PetscCall(VecScatterEnd(mpimat_1->Mvctx, tempVec, mpimat_1->lvec, INSERT_VALUES, SCATTER_FORWARD));
1014: PetscCall(VecGetArray(mpimat_1->lvec, &cpcol_1_state));
1015: /* get 'cpcol_2_state' */
1016: PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, mpimat_2->lvec, INSERT_VALUES, SCATTER_FORWARD));
1017: PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, mpimat_2->lvec, INSERT_VALUES, SCATTER_FORWARD));
1018: PetscCall(VecGetArray(mpimat_2->lvec, &cpcol_2_state));
1019: /* get 'cpcol_2_par_orig' */
1020: for (kk = 0, j = my0; kk < nloc; kk++, j++) {
1021: PetscScalar v = lid_parent_gid[kk];
1023: PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
1024: }
1025: PetscCall(VecAssemblyBegin(tempVec));
1026: PetscCall(VecAssemblyEnd(tempVec));
1027: PetscCall(VecDuplicate(mpimat_2->lvec, &ghost_par_orig2));
1028: PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghost_par_orig2, INSERT_VALUES, SCATTER_FORWARD));
1029: PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghost_par_orig2, INSERT_VALUES, SCATTER_FORWARD));
1030: PetscCall(VecGetArray(ghost_par_orig2, &cpcol_2_par_orig));
1032: PetscCall(VecDestroy(&tempVec));
1033: } /* ismpi */
1034: for (lid = 0; lid < nloc; lid++) {
1035: NState state = lid_state[lid];
1037: if (IS_SELECTED(state)) {
1038: /* steal locals */
1039: ii = matA_1->i;
1040: n = ii[lid + 1] - ii[lid];
1041: idx = matA_1->j + ii[lid];
1042: for (j = 0; j < n; j++) {
1043: PetscInt lidj = idx[j], sgid;
1044: NState statej = lid_state[lidj];
1046: if (statej == DELETED && (sgid = (PetscInt)PetscRealPart(lid_parent_gid[lidj])) != lid + my0) { /* steal local */
1047: lid_parent_gid[lidj] = (PetscScalar)(lid + my0); /* send this if sgid is not local */
1048: if (sgid >= my0 && sgid < Iend) { /* I'm stealing this local from a local sgid */
1049: PetscInt hav = 0, slid = sgid - my0, gidj = lidj + my0;
1050: PetscCDIntNd *pos, *last = NULL;
1052: /* looking for local from local so id_llist_2 works */
1053: PetscCall(PetscCDGetHeadPos(aggs_2, slid, &pos));
1054: while (pos) {
1055: PetscInt gid;
1057: PetscCall(PetscCDIntNdGetID(pos, &gid));
1058: if (gid == gidj) {
1059: PetscCheck(last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "last cannot be null");
1060: PetscCall(PetscCDRemoveNextNode(aggs_2, slid, last));
1061: PetscCall(PetscCDAppendNode(aggs_2, lid, pos));
1062: hav = 1;
1063: break;
1064: } else last = pos;
1065: PetscCall(PetscCDGetNextPos(aggs_2, slid, &pos));
1066: }
1067: if (hav != 1) {
1068: PetscCheck(hav, PETSC_COMM_SELF, PETSC_ERR_PLIB, "failed to find adj in 'selected' lists - structurally unsymmetric matrix");
1069: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "found node %" PetscInt_FMT " times???", hav);
1070: }
1071: } else { /* I'm stealing this local, owned by a ghost */
1072: PetscCheck(sgid == -1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mat has an un-symmetric graph. Use '-%spc_gamg_sym_graph true' to symmetrize the graph or '-%spc_gamg_threshold -1' if the matrix is structurally symmetric.",
1073: ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "");
1074: PetscCall(PetscCDAppendID(aggs_2, lid, lidj + my0));
1075: }
1076: }
1077: } /* local neighbors */
1078: } else if (state == DELETED /* && lid_cprowID_1 */) {
1079: PetscInt sgidold = (PetscInt)PetscRealPart(lid_parent_gid[lid]);
1081: /* see if I have a selected ghost neighbor that will steal me */
1082: if ((ix = lid_cprowID_1[lid]) != -1) {
1083: ii = matB_1->compressedrow.i;
1084: n = ii[ix + 1] - ii[ix];
1085: idx = matB_1->j + ii[ix];
1086: for (j = 0; j < n; j++) {
1087: PetscInt cpid = idx[j];
1088: NState statej = (NState)PetscRealPart(cpcol_1_state[cpid]);
1090: if (IS_SELECTED(statej) && sgidold != statej) { /* ghost will steal this, remove from my list */
1091: lid_parent_gid[lid] = (PetscScalar)statej; /* send who selected */
1092: if (sgidold >= my0 && sgidold < Iend) { /* this was mine */
1093: PetscInt hav = 0, oldslidj = sgidold - my0;
1094: PetscCDIntNd *pos, *last = NULL;
1096: /* remove from 'oldslidj' list */
1097: PetscCall(PetscCDGetHeadPos(aggs_2, oldslidj, &pos));
1098: while (pos) {
1099: PetscInt gid;
1101: PetscCall(PetscCDIntNdGetID(pos, &gid));
1102: if (lid + my0 == gid) {
1103: /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove lid from oldslidj list *\/ */
1104: PetscCheck(last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "last cannot be null");
1105: PetscCall(PetscCDRemoveNextNode(aggs_2, oldslidj, last));
1106: /* ghost (PetscScalar)statej will add this later */
1107: hav = 1;
1108: break;
1109: } else last = pos;
1110: PetscCall(PetscCDGetNextPos(aggs_2, oldslidj, &pos));
1111: }
1112: if (hav != 1) {
1113: PetscCheck(hav, PETSC_COMM_SELF, PETSC_ERR_PLIB, "failed to find (hav=%" PetscInt_FMT ") adj in 'selected' lists - structurally unsymmetric matrix", hav);
1114: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "found node %" PetscInt_FMT " times???", hav);
1115: }
1116: } else {
1117: /* TODO: ghosts remove this later */
1118: }
1119: }
1120: }
1121: }
1122: } /* selected/deleted */
1123: } /* node loop */
1125: if (isMPI) {
1126: PetscScalar *cpcol_2_parent, *cpcol_2_gid;
1127: Vec tempVec, ghostgids2, ghostparents2;
1128: PetscInt cpid, nghost_2;
1129: PetscHMapI gid_cpid;
1131: PetscCall(VecGetSize(mpimat_2->lvec, &nghost_2));
1132: PetscCall(MatCreateVecs(Gmat_2, &tempVec, NULL));
1134: /* get 'cpcol_2_parent' */
1135: for (kk = 0, j = my0; kk < nloc; kk++, j++) PetscCall(VecSetValues(tempVec, 1, &j, &lid_parent_gid[kk], INSERT_VALUES));
1136: PetscCall(VecAssemblyBegin(tempVec));
1137: PetscCall(VecAssemblyEnd(tempVec));
1138: PetscCall(VecDuplicate(mpimat_2->lvec, &ghostparents2));
1139: PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghostparents2, INSERT_VALUES, SCATTER_FORWARD));
1140: PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghostparents2, INSERT_VALUES, SCATTER_FORWARD));
1141: PetscCall(VecGetArray(ghostparents2, &cpcol_2_parent));
1143: /* get 'cpcol_2_gid' */
1144: for (kk = 0, j = my0; kk < nloc; kk++, j++) {
1145: PetscScalar v = (PetscScalar)j;
1147: PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
1148: }
1149: PetscCall(VecAssemblyBegin(tempVec));
1150: PetscCall(VecAssemblyEnd(tempVec));
1151: PetscCall(VecDuplicate(mpimat_2->lvec, &ghostgids2));
1152: PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghostgids2, INSERT_VALUES, SCATTER_FORWARD));
1153: PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghostgids2, INSERT_VALUES, SCATTER_FORWARD));
1154: PetscCall(VecGetArray(ghostgids2, &cpcol_2_gid));
1155: PetscCall(VecDestroy(&tempVec));
1157: /* look for deleted ghosts and add to table */
1158: PetscCall(PetscHMapICreateWithSize(2 * nghost_2 + 1, &gid_cpid));
1159: for (cpid = 0; cpid < nghost_2; cpid++) {
1160: NState state = (NState)PetscRealPart(cpcol_2_state[cpid]);
1162: if (state == DELETED) {
1163: PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
1164: PetscInt sgid_old = (PetscInt)PetscRealPart(cpcol_2_par_orig[cpid]);
1166: if (sgid_old == -1 && sgid_new != -1) {
1167: PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
1169: PetscCall(PetscHMapISet(gid_cpid, gid, cpid));
1170: }
1171: }
1172: }
1174: /* look for deleted ghosts and see if they moved - remove it */
1175: for (lid = 0; lid < nloc; lid++) {
1176: NState state = lid_state[lid];
1178: if (IS_SELECTED(state)) {
1179: PetscCDIntNd *pos, *last = NULL;
1181: /* look for deleted ghosts and see if they moved */
1182: PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
1183: while (pos) {
1184: PetscInt gid;
1186: PetscCall(PetscCDIntNdGetID(pos, &gid));
1187: if (gid < my0 || gid >= Iend) {
1188: PetscCall(PetscHMapIGet(gid_cpid, gid, &cpid));
1189: if (cpid != -1) {
1190: /* a moved ghost - */
1191: /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove 'flid' from list *\/ */
1192: PetscCall(PetscCDRemoveNextNode(aggs_2, lid, last));
1193: } else last = pos;
1194: } else last = pos;
1196: PetscCall(PetscCDGetNextPos(aggs_2, lid, &pos));
1197: } /* loop over list of deleted */
1198: } /* selected */
1199: }
1200: PetscCall(PetscHMapIDestroy(&gid_cpid));
1202: /* look at ghosts, see if they changed - and it */
1203: for (cpid = 0; cpid < nghost_2; cpid++) {
1204: PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
1206: if (sgid_new >= my0 && sgid_new < Iend) { /* this is mine */
1207: PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
1208: PetscInt slid_new = sgid_new - my0, hav = 0;
1209: PetscCDIntNd *pos;
1211: /* search for this gid to see if I have it */
1212: PetscCall(PetscCDGetHeadPos(aggs_2, slid_new, &pos));
1213: while (pos) {
1214: PetscInt gidj;
1216: PetscCall(PetscCDIntNdGetID(pos, &gidj));
1217: PetscCall(PetscCDGetNextPos(aggs_2, slid_new, &pos));
1219: if (gidj == gid) {
1220: hav = 1;
1221: break;
1222: }
1223: }
1224: if (hav != 1) {
1225: /* insert 'flidj' into head of llist */
1226: PetscCall(PetscCDAppendID(aggs_2, slid_new, gid));
1227: }
1228: }
1229: }
1230: PetscCall(VecRestoreArray(mpimat_1->lvec, &cpcol_1_state));
1231: PetscCall(VecRestoreArray(mpimat_2->lvec, &cpcol_2_state));
1232: PetscCall(VecRestoreArray(ghostparents2, &cpcol_2_parent));
1233: PetscCall(VecRestoreArray(ghostgids2, &cpcol_2_gid));
1234: PetscCall(VecDestroy(&ghostgids2));
1235: PetscCall(VecDestroy(&ghostparents2));
1236: PetscCall(VecDestroy(&ghost_par_orig2));
1237: }
1238: PetscCall(PetscFree3(lid_state, lid_parent_gid, lid_cprowID_1));
1239: PetscFunctionReturn(PETSC_SUCCESS);
1240: }
1242: /*
1243: PCGAMGCoarsen_AGG - supports squaring the graph (deprecated) and new graph for
1244: communication of QR data used with HEM and MISk coarsening
1246: Input Parameter:
1247: . a_pc - this
1249: Input/Output Parameter:
1250: . a_Gmat1 - graph to coarsen (in), graph off processor edges for QR gather scatter (out)
1252: Output Parameter:
1253: . agg_lists - list of aggregates
1255: */
1256: static PetscErrorCode PCGAMGCoarsen_AGG(PC a_pc, Mat *a_Gmat1, PetscCoarsenData **agg_lists)
1257: {
1258: PC_MG *mg = (PC_MG *)a_pc->data;
1259: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1260: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
1261: Mat Gmat2, Gmat1 = *a_Gmat1; /* aggressive graph */
1262: IS perm;
1263: PetscInt Istart, Iend, Ii, nloc, bs, nn;
1264: PetscInt *permute, *degree;
1265: PetscBool *bIndexSet;
1266: PetscReal hashfact;
1267: PetscInt iSwapIndex;
1268: PetscRandom random;
1269: MPI_Comm comm;
1271: PetscFunctionBegin;
1272: PetscCall(PetscObjectGetComm((PetscObject)Gmat1, &comm));
1273: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
1274: PetscCall(MatGetLocalSize(Gmat1, &nn, NULL));
1275: PetscCall(MatGetBlockSize(Gmat1, &bs));
1276: PetscCheck(bs == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "bs %" PetscInt_FMT " must be 1", bs);
1277: nloc = nn / bs;
1278: /* get MIS aggs - randomize */
1279: PetscCall(PetscMalloc2(nloc, &permute, nloc, °ree));
1280: PetscCall(PetscCalloc1(nloc, &bIndexSet));
1281: for (Ii = 0; Ii < nloc; Ii++) permute[Ii] = Ii;
1282: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &random));
1283: PetscCall(MatGetOwnershipRange(Gmat1, &Istart, &Iend));
1284: for (Ii = 0; Ii < nloc; Ii++) {
1285: PetscInt nc;
1287: PetscCall(MatGetRow(Gmat1, Istart + Ii, &nc, NULL, NULL));
1288: degree[Ii] = nc;
1289: PetscCall(MatRestoreRow(Gmat1, Istart + Ii, &nc, NULL, NULL));
1290: }
1291: for (Ii = 0; Ii < nloc; Ii++) {
1292: PetscCall(PetscRandomGetValueReal(random, &hashfact));
1293: iSwapIndex = (PetscInt)(hashfact * nloc) % nloc;
1294: if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) {
1295: PetscInt iTemp = permute[iSwapIndex];
1297: permute[iSwapIndex] = permute[Ii];
1298: permute[Ii] = iTemp;
1299: iTemp = degree[iSwapIndex];
1300: degree[iSwapIndex] = degree[Ii];
1301: degree[Ii] = iTemp;
1302: bIndexSet[iSwapIndex] = PETSC_TRUE;
1303: }
1304: }
1305: // apply minimum degree ordering -- NEW
1306: if (pc_gamg_agg->use_minimum_degree_ordering) PetscCall(PetscSortIntWithArray(nloc, degree, permute));
1307: PetscCall(PetscFree(bIndexSet));
1308: PetscCall(PetscRandomDestroy(&random));
1309: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_USE_POINTER, &perm));
1310: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0));
1311: // square graph
1312: if (pc_gamg->current_level < pc_gamg_agg->aggressive_coarsening_levels && pc_gamg_agg->use_aggressive_square_graph) PetscCall(PCGAMGSquareGraph_GAMG(a_pc, Gmat1, &Gmat2));
1313: else Gmat2 = Gmat1;
1314: // switch to old MIS-1 for square graph
1315: if (pc_gamg->current_level < pc_gamg_agg->aggressive_coarsening_levels) {
1316: if (!pc_gamg_agg->use_aggressive_square_graph) PetscCall(MatCoarsenMISKSetDistance(pc_gamg_agg->crs, pc_gamg_agg->aggressive_mis_k)); // hardwire to MIS-2
1317: else PetscCall(MatCoarsenSetType(pc_gamg_agg->crs, MATCOARSENMIS)); // old MIS -- side effect
1318: } else if (pc_gamg_agg->use_aggressive_square_graph && pc_gamg_agg->aggressive_coarsening_levels > 0) { // we reset the MIS
1319: const char *prefix;
1321: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)a_pc, &prefix));
1322: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc_gamg_agg->crs, prefix));
1323: PetscCall(MatCoarsenSetFromOptions(pc_gamg_agg->crs)); // get the default back on non-aggressive levels when square graph switched to old MIS
1324: }
1325: PetscCall(MatCoarsenSetAdjacency(pc_gamg_agg->crs, Gmat2));
1326: PetscCall(MatCoarsenSetStrictAggs(pc_gamg_agg->crs, PETSC_TRUE));
1327: PetscCall(MatCoarsenSetGreedyOrdering(pc_gamg_agg->crs, perm));
1328: PetscCall(MatCoarsenApply(pc_gamg_agg->crs));
1329: PetscCall(MatCoarsenGetData(pc_gamg_agg->crs, agg_lists)); /* output */
1331: PetscCall(ISDestroy(&perm));
1332: PetscCall(PetscFree2(permute, degree));
1333: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0));
1335: if (Gmat2 != Gmat1) { // square graph, we need ghosts for selected
1336: PetscCoarsenData *llist = *agg_lists;
1338: PetscCall(fixAggregatesWithSquare(a_pc, Gmat2, Gmat1, *agg_lists));
1339: PetscCall(MatDestroy(&Gmat1));
1340: *a_Gmat1 = Gmat2; /* output */
1341: PetscCall(PetscCDSetMat(llist, *a_Gmat1)); /* Need a graph with ghosts here */
1342: }
1343: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
1344: PetscFunctionReturn(PETSC_SUCCESS);
1345: }
1347: /*
1348: PCGAMGConstructProlongator_AGG
1350: Input Parameter:
1351: . pc - this
1352: . Amat - matrix on this fine level
1353: . Graph - used to get ghost data for nodes in
1354: . agg_lists - list of aggregates
1355: Output Parameter:
1356: . a_P_out - prolongation operator to the next level
1357: */
1358: static PetscErrorCode PCGAMGConstructProlongator_AGG(PC pc, Mat Amat, PetscCoarsenData *agg_lists, Mat *a_P_out)
1359: {
1360: PC_MG *mg = (PC_MG *)pc->data;
1361: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1362: const PetscInt col_bs = pc_gamg->data_cell_cols;
1363: PetscInt Istart, Iend, nloc, ii, jj, kk, my0, nLocalSelected, bs;
1364: Mat Gmat, Prol;
1365: PetscMPIInt size;
1366: MPI_Comm comm;
1367: PetscReal *data_w_ghost;
1368: PetscInt myCrs0, nbnodes = 0, *flid_fgid;
1369: MatType mtype;
1371: PetscFunctionBegin;
1372: PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
1373: PetscCheck(col_bs >= 1, comm, PETSC_ERR_PLIB, "Column bs cannot be less than 1");
1374: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1375: PetscCallMPI(MPI_Comm_size(comm, &size));
1376: PetscCall(MatGetOwnershipRange(Amat, &Istart, &Iend));
1377: PetscCall(MatGetBlockSize(Amat, &bs));
1378: nloc = (Iend - Istart) / bs;
1379: my0 = Istart / bs;
1380: PetscCheck((Iend - Istart) % bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "(Iend %" PetscInt_FMT " - Istart %" PetscInt_FMT ") not divisible by bs %" PetscInt_FMT, Iend, Istart, bs);
1381: PetscCall(PetscCDGetMat(agg_lists, &Gmat)); // get auxiliary matrix for ghost edges for size > 1
1383: /* get 'nLocalSelected' */
1384: for (ii = 0, nLocalSelected = 0; ii < nloc; ii++) {
1385: PetscBool ise;
1387: /* filter out singletons 0 or 1? */
1388: PetscCall(PetscCDIsEmptyAt(agg_lists, ii, &ise));
1389: if (!ise) nLocalSelected++;
1390: }
1392: /* create prolongator, create P matrix */
1393: PetscCall(MatGetType(Amat, &mtype));
1394: PetscCall(MatCreate(comm, &Prol));
1395: PetscCall(MatSetSizes(Prol, nloc * bs, nLocalSelected * col_bs, PETSC_DETERMINE, PETSC_DETERMINE));
1396: PetscCall(MatSetBlockSizes(Prol, bs, col_bs)); // should this be before MatSetSizes?
1397: PetscCall(MatSetType(Prol, mtype));
1398: #if PetscDefined(HAVE_DEVICE)
1399: PetscBool flg;
1400: PetscCall(MatBoundToCPU(Amat, &flg));
1401: PetscCall(MatBindToCPU(Prol, flg));
1402: if (flg) PetscCall(MatSetBindingPropagates(Prol, PETSC_TRUE));
1403: #endif
1404: PetscCall(MatSeqAIJSetPreallocation(Prol, col_bs, NULL));
1405: PetscCall(MatMPIAIJSetPreallocation(Prol, col_bs, NULL, col_bs, NULL));
1407: /* can get all points "removed" */
1408: PetscCall(MatGetSize(Prol, &kk, &ii));
1409: if (!ii) {
1410: PetscCall(PetscInfo(pc, "%s: No selected points on coarse grid\n", ((PetscObject)pc)->prefix));
1411: PetscCall(MatDestroy(&Prol));
1412: *a_P_out = NULL; /* out */
1413: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1414: PetscFunctionReturn(PETSC_SUCCESS);
1415: }
1416: PetscCall(PetscInfo(pc, "%s: New grid %" PetscInt_FMT " nodes\n", ((PetscObject)pc)->prefix, ii / col_bs));
1417: PetscCall(MatGetOwnershipRangeColumn(Prol, &myCrs0, &kk));
1419: PetscCheck((kk - myCrs0) % col_bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "(kk %" PetscInt_FMT " -myCrs0 %" PetscInt_FMT ") not divisible by col_bs %" PetscInt_FMT, kk, myCrs0, col_bs);
1420: myCrs0 = myCrs0 / col_bs;
1421: PetscCheck((kk / col_bs - myCrs0) == nLocalSelected, PETSC_COMM_SELF, PETSC_ERR_PLIB, "(kk %" PetscInt_FMT "/col_bs %" PetscInt_FMT " - myCrs0 %" PetscInt_FMT ") != nLocalSelected %" PetscInt_FMT ")", kk, col_bs, myCrs0, nLocalSelected);
1423: /* create global vector of data in 'data_w_ghost' */
1424: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0));
1425: if (size > 1) { /* get ghost null space data */
1426: PetscReal *tmp_gdata, *tmp_ldata, *tp2;
1428: PetscCall(PetscMalloc1(nloc, &tmp_ldata));
1429: for (jj = 0; jj < col_bs; jj++) {
1430: for (kk = 0; kk < bs; kk++) {
1431: PetscInt stride;
1432: const PetscReal *tp = PetscSafePointerPlusOffset(pc_gamg->data, jj * bs * nloc + kk);
1434: for (PetscInt ii = 0; ii < nloc; ii++, tp += bs) tmp_ldata[ii] = *tp;
1436: PetscCall(PCGAMGGetDataWithGhosts(Gmat, 1, tmp_ldata, &stride, &tmp_gdata));
1438: if (!jj && !kk) { /* now I know how many total nodes - allocate TODO: move below and do in one 'col_bs' call */
1439: PetscCall(PetscMalloc1(stride * bs * col_bs, &data_w_ghost));
1440: nbnodes = bs * stride;
1441: }
1442: tp2 = PetscSafePointerPlusOffset(data_w_ghost, jj * bs * stride + kk);
1443: for (PetscInt ii = 0; ii < stride; ii++, tp2 += bs) *tp2 = tmp_gdata[ii];
1444: PetscCall(PetscFree(tmp_gdata));
1445: }
1446: }
1447: PetscCall(PetscFree(tmp_ldata));
1448: } else {
1449: nbnodes = bs * nloc;
1450: data_w_ghost = pc_gamg->data;
1451: }
1453: /* get 'flid_fgid' TODO - move up to get 'stride' and do get null space data above in one step (jj loop) */
1454: if (size > 1) {
1455: PetscReal *fid_glid_loc, *fiddata;
1456: PetscInt stride;
1458: PetscCall(PetscMalloc1(nloc, &fid_glid_loc));
1459: for (kk = 0; kk < nloc; kk++) fid_glid_loc[kk] = (PetscReal)(my0 + kk);
1460: PetscCall(PCGAMGGetDataWithGhosts(Gmat, 1, fid_glid_loc, &stride, &fiddata));
1461: PetscCall(PetscMalloc1(stride, &flid_fgid)); /* copy real data to in */
1462: for (kk = 0; kk < stride; kk++) flid_fgid[kk] = (PetscInt)fiddata[kk];
1463: PetscCall(PetscFree(fiddata));
1465: PetscCheck(stride == nbnodes / bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "stride %" PetscInt_FMT " != nbnodes %" PetscInt_FMT "/bs %" PetscInt_FMT, stride, nbnodes, bs);
1466: PetscCall(PetscFree(fid_glid_loc));
1467: } else {
1468: PetscCall(PetscMalloc1(nloc, &flid_fgid));
1469: for (kk = 0; kk < nloc; kk++) flid_fgid[kk] = my0 + kk;
1470: }
1471: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0));
1472: /* get P0 */
1473: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0));
1474: {
1475: PetscReal *data_out = NULL;
1477: PetscCall(formProl0(agg_lists, bs, col_bs, myCrs0, nbnodes, data_w_ghost, flid_fgid, &data_out, Prol));
1478: PetscCall(PetscFree(pc_gamg->data));
1480: pc_gamg->data = data_out;
1481: pc_gamg->data_cell_rows = col_bs;
1482: pc_gamg->data_sz = col_bs * col_bs * nLocalSelected;
1483: }
1484: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0));
1485: if (size > 1) PetscCall(PetscFree(data_w_ghost));
1486: PetscCall(PetscFree(flid_fgid));
1488: *a_P_out = Prol; /* out */
1489: PetscCall(MatViewFromOptions(Prol, NULL, "-pc_gamg_agg_view_initial_prolongation"));
1491: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1492: PetscFunctionReturn(PETSC_SUCCESS);
1493: }
1495: /*
1496: PCGAMGKernelPreservingFilter_AGG - filter the prolongator while preserving the near-null space constraint P*B_c = B
1498: Applies `MatFilter()` to drop small entries, then corrects each row so that
1499: P_filtered * B_c = B (the fine near-null space) is restored.
1501: For nSAvec == 1: rescale each row by B[i] / (P_filtered[i,:] * B_c[J_i]).
1502: For nSAvec > 1: solve a small nSAvec x nSAvec SPD system per row and add
1503: a rank-nSAvec correction to the row entries.
1504: */
1505: static PetscErrorCode PCGAMGKernelPreservingFilter_AGG(PC pc, Mat Prol, PetscReal threshold)
1506: {
1507: PC_MG *mg = (PC_MG *)pc->data;
1508: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1509: const PetscInt nSAvec = pc_gamg->data_cell_rows; /* == data_cell_cols after formProl0 */
1510: PetscInt cStart, cEnd, rStart, rEnd;
1511: const PetscReal *Bc_data = pc_gamg->data;
1512: Vec *Bc_vecs, *B_vecs;
1513: PetscScalar *Bc_arr;
1515: PetscFunctionBegin;
1516: PetscCall(PetscInfo(pc, "Kernel-preserving filter of prolongator with threshold %g, nSAvec=%" PetscInt_FMT "\n", (double)threshold, nSAvec));
1518: PetscCall(MatGetOwnershipRange(Prol, &rStart, &rEnd));
1519: PetscCall(MatGetOwnershipRangeColumn(Prol, &cStart, &cEnd));
1521: /* Step 1: build coarse null-space vectors and compute B = P_original * B_c */
1522: PetscCall(PetscMalloc1(nSAvec, &Bc_vecs));
1523: PetscCall(PetscMalloc1(nSAvec, &B_vecs));
1525: {
1526: PetscInt nloc = cEnd - cStart;
1527: for (PetscInt k = 0; k < nSAvec; k++) {
1528: PetscCall(MatCreateVecs(Prol, &Bc_vecs[k], &B_vecs[k]));
1529: /* fill local entries: Bc_data layout is Bc_data[k * nloc + c] (stride == nloc) */
1530: PetscCall(VecGetArray(Bc_vecs[k], &Bc_arr));
1531: for (PetscInt c = 0; c < nloc; c++) Bc_arr[c] = (PetscScalar)Bc_data[k * nloc + c];
1532: PetscCall(VecRestoreArray(Bc_vecs[k], &Bc_arr));
1533: PetscCall(MatMult(Prol, Bc_vecs[k], B_vecs[k]));
1534: }
1535: }
1537: /* Step 2: apply the threshold filter */
1538: {
1539: PetscBool info_active = PETSC_FALSE;
1540: MatInfo info0, info1;
1541: PetscCall(PetscInfoEnabled(((PetscObject)pc)->classid, &info_active));
1542: if (info_active) PetscCall(MatGetInfo(Prol, MAT_GLOBAL_SUM, &info0));
1543: PetscCall(MatFilter(Prol, threshold, PETSC_TRUE, PETSC_TRUE));
1544: if (info_active) {
1545: PetscCall(MatGetInfo(Prol, MAT_GLOBAL_SUM, &info1));
1546: PetscCall(PetscInfo(pc, "Prolongator filter: nnz before=%g after=%g reduction=%g%%\n", info0.nz_used, info1.nz_used, (info0.nz_used > 0) ? 100.0 * (info0.nz_used - info1.nz_used) / info0.nz_used : 0.0));
1547: }
1548: }
1550: /* Step 3: correct rows to restore P_filtered * B_c = B */
1551: if (nSAvec == 1) {
1552: /*
1553: Scalar case: use `MatMult()` + element-wise scaling + `MatDiagonalScale()`.
1554: scale_i = B_i / (P_filtered * Bc)_i, then P_new = diag(scale) * P_filtered.
1555: Guard against zero denominators (empty rows after filter).
1556: No ghost column access needed.
1557: */
1558: Vec d_vec, scale_vec;
1559: PetscInt n_local;
1560: PetscScalar *s_arr;
1561: const PetscScalar *b_arr, *d_arr;
1563: PetscCall(MatCreateVecs(Prol, NULL, &d_vec));
1564: PetscCall(MatMult(Prol, Bc_vecs[0], d_vec));
1565: PetscCall(VecDuplicate(d_vec, &scale_vec));
1566: PetscCall(VecGetLocalSize(d_vec, &n_local));
1567: PetscCall(VecGetArrayRead(B_vecs[0], &b_arr));
1568: PetscCall(VecGetArrayRead(d_vec, &d_arr));
1569: PetscCall(VecGetArray(scale_vec, &s_arr));
1570: for (PetscInt i = 0; i < n_local; i++) s_arr[i] = (PetscAbsScalar(d_arr[i]) > 0.0) ? b_arr[i] / d_arr[i] : 1.0;
1571: PetscCall(VecRestoreArray(scale_vec, &s_arr));
1572: PetscCall(VecRestoreArrayRead(d_vec, &d_arr));
1573: PetscCall(VecRestoreArrayRead(B_vecs[0], &b_arr));
1574: PetscCall(MatDiagonalScale(Prol, scale_vec, NULL));
1575: PetscCall(VecDestroy(&d_vec));
1576: PetscCall(VecDestroy(&scale_vec));
1577: } else {
1578: /*
1579: Vector case (nSAvec > 1): per-row least-squares correction.
1580: Scatter Bc_data to include ghost column values using Prol's Mvctx,
1581: then build a hash map from global ghost column index to local ghost index
1582: so that `MatGetRow()` global column indices can be mapped to the ghosted array.
1583: */
1584: PetscInt nloc = cEnd - cStart;
1585: PetscInt ghost_stride;
1586: PetscReal *Bc_ghosted = NULL;
1587: const PetscReal *Bc_ghosted_ro;
1588: PetscMPIInt comm_size;
1589: PetscHMapI ghost_gid_to_lid; /* global ghost col index -> local ghost index (0-based) */
1590: PetscInt num_ghosts = 0;
1592: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)Prol), &comm_size));
1593: if (comm_size > 1) {
1594: Mat_MPIAIJ *mpimat = (Mat_MPIAIJ *)Prol->data;
1595: Vec tmp_vec;
1596: PetscScalar *data_arr;
1597: PetscInt nnodes;
1599: PetscCall(VecGetLocalSize(mpimat->lvec, &num_ghosts));
1600: nnodes = nloc + num_ghosts;
1601: ghost_stride = nnodes;
1602: /*
1603: Scatter Bc_data to include ghost column values using Prol's Mvctx.
1604: Cannot use PCGAMGGetDataWithGhosts() because it assumes square matrix
1605: (uses MatGetOwnershipRange() for row indices, but Prol is rectangular).
1606: */
1607: PetscCall(MatCreateVecs(Prol, &tmp_vec, NULL));
1608: PetscCall(PetscMalloc1(nSAvec * nnodes, &Bc_ghosted));
1609: for (PetscInt dir = 0; dir < nSAvec; dir++) {
1610: PetscScalar *tmp_arr;
1611: PetscCall(VecGetArray(tmp_vec, &tmp_arr));
1612: for (PetscInt kk = 0; kk < nloc; kk++) {
1613: PetscReal val = Bc_data[dir * nloc + kk];
1614: Bc_ghosted[dir * nnodes + kk] = val;
1615: tmp_arr[kk] = (PetscScalar)val;
1616: }
1617: PetscCall(VecRestoreArray(tmp_vec, &tmp_arr));
1618: PetscCall(VecScatterBegin(mpimat->Mvctx, tmp_vec, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD));
1619: PetscCall(VecScatterEnd(mpimat->Mvctx, tmp_vec, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD));
1620: PetscCall(VecGetArray(mpimat->lvec, &data_arr));
1621: for (PetscInt g = 0; g < num_ghosts; g++) Bc_ghosted[dir * nnodes + nloc + g] = PetscRealPart(data_arr[g]);
1622: PetscCall(VecRestoreArray(mpimat->lvec, &data_arr));
1623: }
1624: PetscCall(VecDestroy(&tmp_vec));
1625: Bc_ghosted_ro = Bc_ghosted;
1626: /* build hash: global ghost col index -> local ghost index (0-based into ghost portion) */
1627: PetscCall(PetscHMapICreateWithSize(2 * num_ghosts + 1, &ghost_gid_to_lid));
1628: for (PetscInt g = 0; g < num_ghosts; g++) PetscCall(PetscHMapISet(ghost_gid_to_lid, mpimat->garray[g], g));
1629: } else {
1630: /* sequential: no ghosts, ghost_stride == nloc, use Bc_data directly (read-only) */
1631: ghost_stride = nloc;
1632: Bc_ghosted_ro = Bc_data;
1633: PetscCall(PetscHMapICreateWithSize(1, &ghost_gid_to_lid));
1634: }
1636: {
1637: PetscInt nrows = rEnd - rStart, max_ncols = 0;
1638: const PetscScalar **B_arrays;
1639: PetscScalar *work, *new_vals, *G, *rhs, *x, *bc_col;
1640: PetscInt *ghosted_idx, *col_buf;
1641: PetscBLASInt *ipiv;
1642: PetscBLASInt N_b;
1644: PetscCall(PetscMalloc1(nSAvec, &B_arrays));
1645: for (PetscInt k = 0; k < nSAvec; k++) PetscCall(VecGetArrayRead(B_vecs[k], &B_arrays[k]));
1646: /* work: nSAvec*nSAvec Gram + nSAvec rhs + nSAvec solution + nSAvec bc_col scratch */
1647: PetscCall(PetscMalloc1(nSAvec * nSAvec + 3 * nSAvec, &work));
1648: PetscCall(PetscMalloc1(nSAvec, &ipiv));
1649: PetscCall(PetscBLASIntCast(nSAvec, &N_b));
1650: G = work;
1651: rhs = work + nSAvec * nSAvec;
1652: x = rhs + nSAvec;
1653: bc_col = x + nSAvec;
1655: /* find max row width and total nnz for pre-allocation */
1656: {
1657: PetscInt total_nnz = 0;
1658: for (PetscInt row = 0; row < nrows; row++) {
1659: PetscInt ncols;
1660: PetscCall(MatGetRow(Prol, rStart + row, &ncols, NULL, NULL));
1661: if (ncols > max_ncols) max_ncols = ncols;
1662: total_nnz += ncols;
1663: PetscCall(MatRestoreRow(Prol, rStart + row, &ncols, NULL, NULL));
1664: }
1665: /* allocate flat CSR-like buffers to store all corrections before applying */
1666: PetscCall(PetscMalloc1(total_nnz, &new_vals));
1667: PetscCall(PetscMalloc1(total_nnz, &col_buf));
1668: }
1669: PetscCall(PetscMalloc1(max_ncols, &ghosted_idx));
1671: /* Pass 1: read rows, compute corrections, store in flat buffers */
1672: {
1673: PetscInt *row_offsets;
1674: PetscInt offset = 0, n_singular = 0, n_zero_rows = 0, n_corrected = 0, n_underdetermined = 0;
1675: PetscReal max_xnorm = 0.0;
1677: PetscCall(PetscMalloc1(nrows + 1, &row_offsets));
1678: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1680: for (PetscInt row = 0; row < nrows; row++) {
1681: PetscInt ncols;
1682: const PetscInt *cols;
1683: const PetscScalar *vals;
1684: PetscInt grow = rStart + row;
1685: PetscBLASInt NRHS = 1, LDA = N_b, LDB = N_b, info;
1687: row_offsets[row] = offset;
1688: PetscCall(MatGetRow(Prol, grow, &ncols, &cols, &vals));
1689: if (ncols == 0) {
1690: n_zero_rows++;
1691: PetscCall(MatRestoreRow(Prol, grow, &ncols, &cols, &vals));
1692: continue;
1693: }
1695: /* When ncols < nSAvec the Gram matrix G is rank-deficient by construction;
1696: skip correction for this row (keep filtered values as-is).
1697: Note: the near-null space constraint P*Bc = B is NOT enforced for these rows.
1698: This typically occurs at boundary or isolated nodes where few coarse neighbors
1699: remain after filtering; the impact on convergence is generally small. */
1700: if (ncols < nSAvec) {
1701: n_underdetermined++;
1702: for (PetscInt j = 0; j < ncols; j++) {
1703: col_buf[offset + j] = cols[j];
1704: new_vals[offset + j] = vals[j];
1705: }
1706: offset += ncols;
1707: PetscCall(MatRestoreRow(Prol, grow, &ncols, &cols, &vals));
1708: continue;
1709: }
1711: /* map global column indices to ghosted array indices and save cols */
1712: for (PetscInt j = 0; j < ncols; j++) {
1713: col_buf[offset + j] = cols[j];
1714: if (cols[j] >= cStart && cols[j] < cEnd) ghosted_idx[j] = cols[j] - cStart;
1715: else {
1716: PetscInt g = -1;
1717: PetscCall(PetscHMapIGet(ghost_gid_to_lid, cols[j], &g));
1718: PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Off-diagonal column %" PetscInt_FMT " not found in ghost map for prolongator filter", cols[j]);
1719: ghosted_idx[j] = nloc + g;
1720: }
1721: }
1723: for (PetscInt i = 0; i < nSAvec * nSAvec; i++) G[i] = 0.0;
1725: /* rhs[k] = B[row,k] - sum_j P[row,j] * Bc[ghosted_idx[j], k] */
1726: for (PetscInt k = 0; k < nSAvec; k++) {
1727: PetscScalar dot = 0.0;
1728: for (PetscInt j = 0; j < ncols; j++) dot += vals[j] * (PetscScalar)Bc_ghosted_ro[k * ghost_stride + ghosted_idx[j]];
1729: rhs[k] = B_arrays[k][row] - dot;
1730: }
1732: /* G[k1,k2] = sum_j Bc[j,k1] * Bc[j,k2] using pre-gathered bc_col */
1733: for (PetscInt j = 0; j < ncols; j++) {
1734: PetscInt gidx = ghosted_idx[j];
1735: for (PetscInt k = 0; k < nSAvec; k++) bc_col[k] = (PetscScalar)Bc_ghosted_ro[k * ghost_stride + gidx];
1736: for (PetscInt k1 = 0; k1 < nSAvec; k1++)
1737: for (PetscInt k2 = k1; k2 < nSAvec; k2++) G[k1 * nSAvec + k2] += bc_col[k1] * bc_col[k2];
1738: }
1739: /* fill lower triangle from upper (G is symmetric) */
1740: for (PetscInt k1 = 1; k1 < nSAvec; k1++)
1741: for (PetscInt k2 = 0; k2 < k1; k2++) G[k1 * nSAvec + k2] = G[k2 * nSAvec + k1];
1743: /* solve G * x = rhs */
1744: for (PetscInt i = 0; i < nSAvec; i++) x[i] = rhs[i];
1745: PetscCallBLAS("LAPACKgesv", LAPACKgesv_(&N_b, &NRHS, G, &LDA, ipiv, x, &LDB, &info));
1746: if (info != 0) {
1747: /* G is singular despite ncols >= nSAvec (Bc columns linearly dependent);
1748: keep filtered values as-is (near-null space constraint not enforced for this row) */
1749: n_singular++;
1750: for (PetscInt j = 0; j < ncols; j++) new_vals[offset + j] = vals[j];
1751: offset += ncols;
1752: PetscCall(MatRestoreRow(Prol, grow, &ncols, &cols, &vals));
1753: continue;
1754: }
1756: /* track ||x||^2 */
1757: {
1758: PetscReal xnorm2 = 0.0;
1759: for (PetscInt k = 0; k < nSAvec; k++) xnorm2 += PetscSqr(PetscAbsScalar(x[k]));
1760: if (xnorm2 > max_xnorm) max_xnorm = xnorm2;
1761: }
1762: n_corrected++;
1764: /* new_vals[j] = vals[j] + sum_k Bc[ghosted_idx[j],k] * x[k] */
1765: for (PetscInt j = 0; j < ncols; j++) {
1766: PetscScalar delta = 0.0;
1767: PetscInt gidx = ghosted_idx[j];
1768: for (PetscInt k = 0; k < nSAvec; k++) delta += (PetscScalar)Bc_ghosted_ro[k * ghost_stride + gidx] * x[k];
1769: new_vals[offset + j] = vals[j] + delta;
1770: }
1771: offset += ncols;
1772: PetscCall(MatRestoreRow(Prol, grow, &ncols, &cols, &vals));
1773: }
1774: row_offsets[nrows] = offset;
1775: PetscCall(PetscFPTrapPop());
1776: PetscCall(PetscInfo(pc, "PCGAMGKernelPreservingFilter_AGG: nrows=%" PetscInt_FMT " corrected=%" PetscInt_FMT " zero_rows=%" PetscInt_FMT " underdetermined(ncols<nSAvec)=%" PetscInt_FMT " singular_G=%" PetscInt_FMT " max_xnorm2=%g\n", nrows, n_corrected, n_zero_rows, n_underdetermined, n_singular, (double)max_xnorm));
1778: /* Pass 2: apply all corrections at once */
1779: for (PetscInt row = 0; row < nrows; row++) {
1780: PetscInt grow = rStart + row;
1781: PetscInt nc = row_offsets[row + 1] - row_offsets[row];
1782: if (nc > 0) PetscCall(MatSetValues(Prol, 1, &grow, nc, col_buf + row_offsets[row], new_vals + row_offsets[row], INSERT_VALUES));
1783: }
1784: PetscCall(PetscFree(row_offsets));
1785: }
1787: for (PetscInt k = 0; k < nSAvec; k++) PetscCall(VecRestoreArrayRead(B_vecs[k], &B_arrays[k]));
1788: PetscCall(PetscFree(B_arrays));
1789: PetscCall(PetscFree(work));
1790: PetscCall(PetscFree(ipiv));
1791: PetscCall(PetscFree(ghosted_idx));
1792: PetscCall(PetscFree(new_vals));
1793: PetscCall(PetscFree(col_buf));
1794: }
1796: PetscCall(PetscHMapIDestroy(&ghost_gid_to_lid));
1797: if (comm_size > 1) PetscCall(PetscFree(Bc_ghosted));
1798: }
1800: PetscCall(MatAssemblyBegin(Prol, MAT_FINAL_ASSEMBLY));
1801: PetscCall(MatAssemblyEnd(Prol, MAT_FINAL_ASSEMBLY));
1803: for (PetscInt k = 0; k < nSAvec; k++) {
1804: PetscCall(VecDestroy(&Bc_vecs[k]));
1805: PetscCall(VecDestroy(&B_vecs[k]));
1806: }
1807: PetscCall(PetscFree(Bc_vecs));
1808: PetscCall(PetscFree(B_vecs));
1809: PetscFunctionReturn(PETSC_SUCCESS);
1810: }
1812: /*
1813: PCGAMGOptimizeProlongator_AGG - given the initial prolongator optimizes it by smoothed aggregation pc_gamg_agg->nsmooths times
1815: Input Parameter:
1816: . pc - this
1817: . Amat - matrix on this fine level
1818: In/Output Parameter:
1819: . a_P - prolongation operator to the next level
1820: */
1821: static PetscErrorCode PCGAMGOptimizeProlongator_AGG(PC pc, Mat Amat, Mat *a_P)
1822: {
1823: PC_MG *mg = (PC_MG *)pc->data;
1824: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1825: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
1826: Mat Prol = *a_P;
1827: MPI_Comm comm;
1828: KSP eksp;
1829: Vec bb, xx;
1830: PC epc;
1831: PetscReal alpha, emax, emin;
1833: PetscFunctionBegin;
1834: PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
1835: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0));
1837: /* compute maximum singular value of operator to be used in smoother */
1838: if (0 < pc_gamg_agg->nsmooths) {
1839: /* get eigen estimates */
1840: if (pc_gamg->emax > 0) {
1841: emin = pc_gamg->emin;
1842: emax = pc_gamg->emax;
1843: } else {
1844: const char *prefix;
1846: PetscCall(MatCreateVecs(Amat, &bb, NULL));
1847: PetscCall(MatCreateVecs(Amat, &xx, NULL));
1848: PetscCall(KSPSetNoisy_Private(Amat, bb));
1850: PetscCall(KSPCreate(comm, &eksp));
1851: PetscCall(KSPSetNestLevel(eksp, pc->kspnestlevel));
1852: PetscCall(PCGetOptionsPrefix(pc, &prefix));
1853: PetscCall(KSPSetOptionsPrefix(eksp, prefix));
1854: PetscCall(KSPAppendOptionsPrefix(eksp, "pc_gamg_esteig_"));
1855: {
1856: PetscBool isset, sflg;
1858: PetscCall(MatIsSPDKnown(Amat, &isset, &sflg));
1859: if (isset && sflg) PetscCall(KSPSetType(eksp, KSPCG));
1860: }
1861: PetscCall(KSPSetErrorIfNotConverged(eksp, pc->erroriffailure));
1862: PetscCall(KSPSetNormType(eksp, KSP_NORM_NONE));
1864: PetscCall(KSPSetInitialGuessNonzero(eksp, PETSC_FALSE));
1865: PetscCall(KSPSetOperators(eksp, Amat, Amat));
1867: PetscCall(KSPGetPC(eksp, &epc));
1868: PetscCall(PCSetType(epc, PCJACOBI)); /* smoother in smoothed agg. */
1870: PetscCall(KSPSetTolerances(eksp, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, 10)); // 10 is safer, but 5 is often fine, can override with -pc_gamg_esteig_ksp_max_it -mg_levels_ksp_chebyshev_esteig 0,0.25,0,1.2
1872: PetscCall(KSPSetFromOptions(eksp));
1873: PetscCall(KSPSetComputeSingularValues(eksp, PETSC_TRUE));
1874: PetscCall(KSPSolve(eksp, bb, xx));
1875: PetscCall(KSPCheckSolve(eksp, pc, xx));
1877: PetscCall(KSPComputeExtremeSingularValues(eksp, &emax, &emin));
1878: PetscCall(PetscInfo(pc, "%s: Smooth P0: max eigen=%e min=%e PC=%s\n", ((PetscObject)pc)->prefix, (double)emax, (double)emin, PCJACOBI));
1879: PetscCall(VecDestroy(&xx));
1880: PetscCall(VecDestroy(&bb));
1881: PetscCall(KSPDestroy(&eksp));
1882: }
1883: if (pc_gamg->use_sa_esteig) {
1884: mg->min_eigen_DinvA[pc_gamg->current_level] = emin;
1885: mg->max_eigen_DinvA[pc_gamg->current_level] = emax;
1886: PetscCall(PetscInfo(pc, "%s: Smooth P0: level %" PetscInt_FMT ", cache spectra %g %g\n", ((PetscObject)pc)->prefix, pc_gamg->current_level, (double)emin, (double)emax));
1887: } else {
1888: mg->min_eigen_DinvA[pc_gamg->current_level] = 0;
1889: mg->max_eigen_DinvA[pc_gamg->current_level] = 0;
1890: }
1891: } else {
1892: mg->min_eigen_DinvA[pc_gamg->current_level] = 0;
1893: mg->max_eigen_DinvA[pc_gamg->current_level] = 0;
1894: }
1896: /* smooth P0 */
1897: if (pc_gamg_agg->nsmooths > 0) {
1898: Vec diag;
1900: /* TODO: Set a PCFailedReason and exit the building of the AMG preconditioner */
1901: PetscCheck(emax != 0.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Computed maximum singular value as zero");
1903: PetscCall(MatCreateVecs(Amat, &diag, NULL));
1904: PetscCall(MatGetDiagonal(Amat, diag)); /* effectively PCJACOBI */
1905: PetscCall(VecReciprocal(diag));
1907: for (PetscInt jj = 0; jj < pc_gamg_agg->nsmooths; jj++) {
1908: Mat tMat;
1910: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0));
1911: /*
1912: Smooth aggregation on the prolongator
1914: P_{i} := (I - 1.4/emax D^{-1}A) P_i\{i-1}
1915: */
1916: PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0));
1917: PetscCall(MatMatMult(Amat, Prol, MAT_INITIAL_MATRIX, PETSC_CURRENT, &tMat));
1918: PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0));
1919: PetscCall(MatProductClear(tMat));
1920: PetscCall(MatDiagonalScale(tMat, diag, NULL));
1922: /* TODO: Document the 1.4 and don't hardwire it in this routine */
1923: alpha = -1.4 / emax;
1924: PetscCall(MatAYPX(tMat, alpha, Prol, SUBSET_NONZERO_PATTERN));
1925: PetscCall(MatDestroy(&Prol));
1926: Prol = tMat;
1927: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0));
1928: }
1929: PetscCall(VecDestroy(&diag));
1930: }
1931: if (pc_gamg->prolongator_filter > 0.0) PetscCall(PCGAMGKernelPreservingFilter_AGG(pc, Prol, pc_gamg->prolongator_filter));
1932: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0));
1933: PetscCall(MatViewFromOptions(Prol, NULL, "-pc_gamg_agg_view_prolongation"));
1934: *a_P = Prol;
1935: PetscFunctionReturn(PETSC_SUCCESS);
1936: }
1938: /*MC
1939: PCGAMGAGG - Smooth aggregation, {cite}`vanek1996algebraic`, {cite}`vanek2001convergence`, variant of PETSc's algebraic multigrid (`PCGAMG`) preconditioner
1941: Options Database Keys:
1942: + -pc_gamg_agg_nsmooths nsmooth - number of smoothing steps to use with smooth aggregation to construct prolongation
1943: . -pc_gamg_prolongator_filter thr - filter small entries from the prolongator, preserving the near-null space (0=disabled, 0.0025=typical)
1944: . -pc_gamg_aggressive_coarsening n - number of aggressive coarsening (MIS-2 or square graph) levels from finest.
1945: . -pc_gamg_aggressive_square_graph (true|false) - use square graph ($A^T A$), alternative is MIS-k (k=2), for aggressive coarsening
1946: . -pc_gamg_mis_k_minimum_degree_ordering (true|false) - use minimum degree ordering in greedy MIS algorithm
1947: . -pc_gamg_asm_hem_aggs n - number of HEM aggregation steps for ASM smoother
1948: - -pc_gamg_aggressive_mis_k n - number (k) distance in MIS coarsening (>2 is 'aggressive')
1950: Level: intermediate
1952: Notes:
1953: To obtain good performance for `PCGAMG` for vector valued problems you must
1954: call `MatSetBlockSize()` to indicate the number of degrees of freedom per grid point.
1955: Call `MatSetNearNullSpace()` (or `PCSetCoordinates()` if solving the equations of elasticity) to indicate the near null space of the operator
1957: When `-pc_gamg_aggressive_square_graph` is used, the coarsening is obtained by first squaring the graph and then applying, by default, a
1958: MIS-1 coarsening with `MatCoarsenApply()` on the squared graph.
1960: The many options for `PCMG` and `PCGAMG` such as controlling the smoothers on each level etc. also work for `PCGAMGAGG`
1962: .seealso: `PCGAMG`, [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCCreate()`, `PCSetType()`,
1963: `MatSetBlockSize()`, `PCMGType`, `PCSetCoordinates()`, `MatSetNearNullSpace()`, `PCGAMGSetType()`,
1964: `PCGAMGAGG`, `PCGAMGGEO`, `PCGAMGCLASSICAL`, `PCGAMGSetProcEqLim()`, `PCGAMGSetCoarseEqLim()`, `PCGAMGSetRepartition()`, `PCGAMGRegister()`,
1965: `PCGAMGSetReuseInterpolation()`, `PCGAMGASMSetUseAggs()`, `PCGAMGSetParallelCoarseGridSolve()`, `PCGAMGSetNlevels()`, `PCGAMGSetThreshold()`,
1966: `PCGAMGGetType()`, `PCGAMGSetUseSAEstEig()`
1967: M*/
1968: PetscErrorCode PCCreateGAMG_AGG(PC pc)
1969: {
1970: PC_MG *mg = (PC_MG *)pc->data;
1971: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1972: PC_GAMG_AGG *pc_gamg_agg;
1974: PetscFunctionBegin;
1975: /* create sub context for SA */
1976: PetscCall(PetscNew(&pc_gamg_agg));
1977: pc_gamg->subctx = pc_gamg_agg;
1979: pc_gamg->ops->setfromoptions = PCSetFromOptions_GAMG_AGG;
1980: pc_gamg->ops->destroy = PCDestroy_GAMG_AGG;
1981: /* reset does not do anything; setup not virtual */
1983: /* set internal function pointers */
1984: pc_gamg->ops->creategraph = PCGAMGCreateGraph_AGG;
1985: pc_gamg->ops->coarsen = PCGAMGCoarsen_AGG;
1986: pc_gamg->ops->prolongator = PCGAMGConstructProlongator_AGG;
1987: pc_gamg->ops->optprolongator = PCGAMGOptimizeProlongator_AGG;
1988: pc_gamg->ops->createdefaultdata = PCSetData_AGG;
1989: pc_gamg->ops->view = PCView_GAMG_AGG;
1991: pc_gamg_agg->nsmooths = 1;
1992: pc_gamg_agg->aggressive_coarsening_levels = 1;
1993: pc_gamg_agg->use_aggressive_square_graph = PETSC_TRUE;
1994: pc_gamg_agg->use_minimum_degree_ordering = PETSC_FALSE;
1995: pc_gamg_agg->use_low_mem_filter = PETSC_FALSE;
1996: pc_gamg_agg->aggressive_mis_k = 2;
1997: pc_gamg_agg->graph_symmetrize = PETSC_TRUE;
1999: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", PCGAMGSetNSmooths_AGG));
2000: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", PCGAMGSetAggressiveLevels_AGG));
2001: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveSquareGraph_C", PCGAMGSetAggressiveSquareGraph_AGG));
2002: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetMinDegreeOrdering_C", PCGAMGMISkSetMinDegreeOrdering_AGG));
2003: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetLowMemoryFilter_C", PCGAMGSetLowMemoryFilter_AGG));
2004: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetAggressive_C", PCGAMGMISkSetAggressive_AGG));
2005: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetGraphSymmetrize_C", PCGAMGSetGraphSymmetrize_AGG));
2006: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetProlongatorFilter_C", PCGAMGSetProlongatorFilter_AGG));
2007: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGGetProlongatorFilter_C", PCGAMGGetProlongatorFilter_AGG));
2008: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_AGG));
2009: PetscFunctionReturn(PETSC_SUCCESS);
2010: }