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>
9: #include <petscsf.h>
11: typedef struct {
12: PetscInt nsmooths; // number of smoothing steps to construct prolongation
13: PetscInt aggressive_coarsening_levels; // number of aggressive coarsening levels (square or MISk)
14: PetscInt aggressive_mis_k; // the k in MIS-k
15: PetscBool use_aggressive_square_graph;
16: PetscBool use_minimum_degree_ordering;
17: PetscBool use_low_mem_filter;
18: PetscBool graph_symmetrize;
19: MatCoarsen crs;
20: } PC_GAMG_AGG;
22: /*@
23: PCGAMGSetNSmooths - Set number of smoothing steps (1 is typical) used to construct the prolongation operator
25: Logically Collective
27: Input Parameters:
28: + pc - the preconditioner context
29: - n - the number of smooths, default is 1
31: Options Database Key:
32: . -pc_gamg_agg_nsmooths nsmooth - number of smoothing steps to use
34: Level: intermediate
36: Note:
37: This is a different concept from the number smoothing steps used during the linear solution process which
38: can be set with `-mg_levels_ksp_max_it`
40: Developer Note:
41: This should be named `PCGAMGAGGSetNSmooths()`.
43: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCMG`, `PCGAMG`
44: @*/
45: PetscErrorCode PCGAMGSetNSmooths(PC pc, PetscInt n)
46: {
47: PetscFunctionBegin;
50: PetscTryMethod(pc, "PCGAMGSetNSmooths_C", (PC, PetscInt), (pc, n));
51: PetscFunctionReturn(PETSC_SUCCESS);
52: }
54: static PetscErrorCode PCGAMGSetNSmooths_AGG(PC pc, PetscInt n)
55: {
56: PC_MG *mg = (PC_MG *)pc->data;
57: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
58: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
60: PetscFunctionBegin;
61: pc_gamg_agg->nsmooths = n;
62: PetscFunctionReturn(PETSC_SUCCESS);
63: }
65: /*@
66: PCGAMGSetAggressiveLevels - Use aggressive coarsening on first n levels
68: Logically Collective
70: Input Parameters:
71: + pc - the preconditioner context
72: - n - 0, 1 or more, the default is 1
74: Options Database Key:
75: . -pc_gamg_aggressive_coarsening n - the number of coarsenings to do aggressively
77: Level: intermediate
79: Note:
80: By default, aggressive coarsening squares the matrix (computes $A^T A$) before coarsening.
81: Calling `PCGAMGSetAggressiveSquareGraph()` with a value of `PETSC_FALSE` changes the aggressive coarsening strategy to use MIS-k, see `PCGAMGMISkSetAggressive()`.
83: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGMISkSetAggressive()`,
84: `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`, `PCGAMGSetLowMemoryFilter()`
85: @*/
86: PetscErrorCode PCGAMGSetAggressiveLevels(PC pc, PetscInt n)
87: {
88: PetscFunctionBegin;
91: PetscTryMethod(pc, "PCGAMGSetAggressiveLevels_C", (PC, PetscInt), (pc, n));
92: PetscFunctionReturn(PETSC_SUCCESS);
93: }
95: /*@
96: PCGAMGMISkSetAggressive - Number (k) distance in MIS coarsening (> 2 is aggressive)
98: Logically Collective
100: Input Parameters:
101: + pc - the preconditioner context
102: - n - 1 or more (default = 2)
104: Options Database Key:
105: . -pc_gamg_aggressive_mis_k n - the distance to use
107: Level: intermediate
109: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`,
110: `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`, `PCGAMGSetLowMemoryFilter()`
111: @*/
112: PetscErrorCode PCGAMGMISkSetAggressive(PC pc, PetscInt n)
113: {
114: PetscFunctionBegin;
117: PetscTryMethod(pc, "PCGAMGMISkSetAggressive_C", (PC, PetscInt), (pc, n));
118: PetscFunctionReturn(PETSC_SUCCESS);
119: }
121: /*@
122: 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
124: Logically Collective
126: Input Parameters:
127: + pc - the preconditioner context
128: - b - default true
130: Options Database Key:
131: . -pc_gamg_aggressive_square_graph (true|false) - whether to use the graph square to aggressively coarsen
133: Level: intermediate
135: Notes:
136: If `b` is `PETSC_FALSE` then MIS-k is used for aggressive coarsening, see `PCGAMGMISkSetAggressive()`
138: 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
139: that converges faster.
141: .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()`
142: @*/
143: PetscErrorCode PCGAMGSetAggressiveSquareGraph(PC pc, PetscBool b)
144: {
145: PetscFunctionBegin;
148: PetscTryMethod(pc, "PCGAMGSetAggressiveSquareGraph_C", (PC, PetscBool), (pc, b));
149: PetscFunctionReturn(PETSC_SUCCESS);
150: }
152: /*@
153: PCGAMGMISkSetMinDegreeOrdering - Use minimum degree ordering in greedy MIS algorithm
155: Logically Collective
157: Input Parameters:
158: + pc - the preconditioner context
159: - b - default false
161: Options Database Key:
162: . -pc_gamg_mis_k_minimum_degree_ordering (true|false) - use the minimum degree ordering
164: Level: intermediate
166: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`,
167: `PCGAMGSetAggressiveLevels()`, `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGSetLowMemoryFilter()`
168: @*/
169: PetscErrorCode PCGAMGMISkSetMinDegreeOrdering(PC pc, PetscBool b)
170: {
171: PetscFunctionBegin;
174: PetscTryMethod(pc, "PCGAMGMISkSetMinDegreeOrdering_C", (PC, PetscBool), (pc, b));
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: /*@
179: PCGAMGSetLowMemoryFilter - Use low memory graph/matrix filter
181: Logically Collective
183: Input Parameters:
184: + pc - the preconditioner context
185: - b - default false
187: Options Database Key:
188: . -pc_gamg_low_memory_threshold_filter (true|false) - use the low memory filter
190: Level: intermediate
192: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`,
193: `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`
194: @*/
195: PetscErrorCode PCGAMGSetLowMemoryFilter(PC pc, PetscBool b)
196: {
197: PetscFunctionBegin;
200: PetscTryMethod(pc, "PCGAMGSetLowMemoryFilter_C", (PC, PetscBool), (pc, b));
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: /*@
205: 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
207: Logically Collective
209: Input Parameters:
210: + pc - the preconditioner context
211: - b - default true
213: Options Database Key:
214: . -pc_gamg_graph_symmetrize (true|false) - symmetrize the graph
216: Level: intermediate
218: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`, `MatCreateGraph()`,
219: `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`
220: @*/
221: PetscErrorCode PCGAMGSetGraphSymmetrize(PC pc, PetscBool b)
222: {
223: PetscFunctionBegin;
226: PetscTryMethod(pc, "PCGAMGSetGraphSymmetrize_C", (PC, PetscBool), (pc, b));
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
230: /*@
231: PCGAMGSetProlongatorFilter - Set the relative threshold for block filtering of the prolongator (a kernel-preserving correction is applied afterward)
233: Logically Collective
235: Input Parameters:
236: + pc - the preconditioner context
237: - thr - relative threshold in [0,1); the block of prolongator entries coupling a fine node to a coarse node is dropped when its Frobenius norm is below `thr` times the
238: largest such block norm in that fine node's block row (0 disables filtering)
240: Options Database Key:
241: . -pc_gamg_prolongator_filter thr - relative threshold for block filtering of the prolongator (0=disabled, 0.01-0.1=typical)
243: Level: intermediate
245: Notes:
246: Each fine node corresponds to a block of rows (one per degree of freedom of the node, as given by the block size of the operator) and each coarse node to a
247: block of columns (one per near-null space vector), so the filtering drops small dense sub-blocks of the prolongator, not individual entries. The threshold is
248: relative to the largest block Frobenius norm in the same fine-node block row so the decision is invariant to the differing scales of the near-null space modes.
249: The comparison is strict, so the strongest block of a fine node always survives. On coarser levels the threshold is scaled by `PCGAMGSetProlongatorFilterScale()`.
250: Dropping whole blocks (rather than individual entries) keeps complete coarse-node blocks in every surviving fine row, so the near-null space correction below
251: remains full rank. The dropped entries are removed from the sparsity pattern with `MatEliminateZeros()`; on matrix types that do not implement it, they are
252: zeroed but remain in the pattern, so the coarse operators are unchanged in structure and the complexity and memory reduction is not realized (reported with `-info`).
254: After filtering, each row of the prolongator is corrected so that the filtered prolongator still reproduces the near-null space exactly, that is, P applied to the coarse
255: representation of the near-null space equals the fine near-null space. With a single near-null space vector each row is simply rescaled; with several, a small symmetric
256: positive-definite system (of size the number of near-null space vectors) is solved for each row and the resulting correction, a combination of the coarse near-null space
257: vectors, is added to the surviving entries of the row. Rows with fewer surviving entries than near-null space vectors are left uncorrected, as are, in the single-vector
258: case, rows whose near-null space entry is zero or whose required scale factor would be extremely large (an empty or nearly empty filtered row).
260: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGGetProlongatorFilter()`, `PCGAMGSetProlongatorFilterScale()`, `PCGAMGSetLowMemoryFilter()`
261: @*/
262: PetscErrorCode PCGAMGSetProlongatorFilter(PC pc, PetscReal thr)
263: {
264: PetscFunctionBegin;
267: PetscTryMethod(pc, "PCGAMGSetProlongatorFilter_C", (PC, PetscReal), (pc, thr));
268: PetscFunctionReturn(PETSC_SUCCESS);
269: }
271: /*@
272: PCGAMGGetProlongatorFilter - Get the relative threshold for block filtering of the prolongator
274: Not Collective
276: Input Parameter:
277: . pc - the preconditioner context
279: Output Parameter:
280: . thr - relative block-filtering threshold; the block of prolongator entries coupling a fine node to a coarse node is dropped when its Frobenius norm is below `thr` times
281: the largest such block norm in that fine node's block row (0 disables filtering, see `PCGAMGSetProlongatorFilter()`)
283: Level: intermediate
285: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetProlongatorFilter()`, `PCGAMGSetProlongatorFilterScale()`, `PCGAMGSetLowMemoryFilter()`
286: @*/
287: PetscErrorCode PCGAMGGetProlongatorFilter(PC pc, PetscReal *thr)
288: {
289: PetscFunctionBegin;
291: PetscAssertPointer(thr, 2);
292: PetscUseMethod(pc, "PCGAMGGetProlongatorFilter_C", (PC, PetscReal *), (pc, thr));
293: PetscFunctionReturn(PETSC_SUCCESS);
294: }
296: /*@
297: PCGAMGSetProlongatorFilterScale - Set the per-level scaling of the prolongator filter threshold (see `PCGAMGSetProlongatorFilter()`)
299: Logically Collective
301: Input Parameters:
302: + pc - the preconditioner context
303: - scale - per-level multiplier in [0,1]; the effective threshold on level l is `prolongator_filter` times `scale` raised to the power l, where level 0 is the finest
305: Options Database Key:
306: . -pc_gamg_prolongator_filter_scale scale - per-level scaling of the prolongator filter threshold (1.0=default)
308: Level: intermediate
310: Note:
311: A scale below 1 filters less aggressively on the coarser levels, where the prolongator is denser. Values above 1, which would make coarser levels filter more
312: aggressively than the finest, are not allowed. A scale of 0 disables filtering on all levels but the finest.
314: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetProlongatorFilter()`, `PCGAMGGetProlongatorFilterScale()`
315: @*/
316: PetscErrorCode PCGAMGSetProlongatorFilterScale(PC pc, PetscReal scale)
317: {
318: PetscFunctionBegin;
321: PetscTryMethod(pc, "PCGAMGSetProlongatorFilterScale_C", (PC, PetscReal), (pc, scale));
322: PetscFunctionReturn(PETSC_SUCCESS);
323: }
325: /*@
326: PCGAMGGetProlongatorFilterScale - Get the per-level scaling of the prolongator filter threshold
328: Not Collective
330: Input Parameter:
331: . pc - the preconditioner context
333: Output Parameter:
334: . scale - per-level multiplier in [0,1]; the effective threshold on level l is `prolongator_filter` times `scale` raised to the power l, where level 0 is the finest
336: Level: intermediate
338: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetProlongatorFilter()`, `PCGAMGSetProlongatorFilterScale()`
339: @*/
340: PetscErrorCode PCGAMGGetProlongatorFilterScale(PC pc, PetscReal *scale)
341: {
342: PetscFunctionBegin;
344: PetscAssertPointer(scale, 2);
345: PetscUseMethod(pc, "PCGAMGGetProlongatorFilterScale_C", (PC, PetscReal *), (pc, scale));
346: PetscFunctionReturn(PETSC_SUCCESS);
347: }
349: static PetscErrorCode PCGAMGSetAggressiveLevels_AGG(PC pc, PetscInt n)
350: {
351: PC_MG *mg = (PC_MG *)pc->data;
352: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
353: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
355: PetscFunctionBegin;
356: pc_gamg_agg->aggressive_coarsening_levels = n;
357: PetscFunctionReturn(PETSC_SUCCESS);
358: }
360: static PetscErrorCode PCGAMGMISkSetAggressive_AGG(PC pc, PetscInt n)
361: {
362: PC_MG *mg = (PC_MG *)pc->data;
363: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
364: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
366: PetscFunctionBegin;
367: pc_gamg_agg->aggressive_mis_k = n;
368: PetscFunctionReturn(PETSC_SUCCESS);
369: }
371: static PetscErrorCode PCGAMGSetAggressiveSquareGraph_AGG(PC pc, PetscBool b)
372: {
373: PC_MG *mg = (PC_MG *)pc->data;
374: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
375: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
377: PetscFunctionBegin;
378: pc_gamg_agg->use_aggressive_square_graph = b;
379: PetscFunctionReturn(PETSC_SUCCESS);
380: }
382: static PetscErrorCode PCGAMGSetLowMemoryFilter_AGG(PC pc, PetscBool b)
383: {
384: PC_MG *mg = (PC_MG *)pc->data;
385: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
386: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
388: PetscFunctionBegin;
389: pc_gamg_agg->use_low_mem_filter = b;
390: PetscFunctionReturn(PETSC_SUCCESS);
391: }
393: static PetscErrorCode PCGAMGSetGraphSymmetrize_AGG(PC pc, PetscBool b)
394: {
395: PC_MG *mg = (PC_MG *)pc->data;
396: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
397: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
399: PetscFunctionBegin;
400: pc_gamg_agg->graph_symmetrize = b;
401: PetscFunctionReturn(PETSC_SUCCESS);
402: }
404: static PetscErrorCode PCGAMGMISkSetMinDegreeOrdering_AGG(PC pc, PetscBool b)
405: {
406: PC_MG *mg = (PC_MG *)pc->data;
407: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
408: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
410: PetscFunctionBegin;
411: pc_gamg_agg->use_minimum_degree_ordering = b;
412: PetscFunctionReturn(PETSC_SUCCESS);
413: }
415: static PetscErrorCode PCGAMGSetProlongatorFilter_AGG(PC pc, PetscReal thr)
416: {
417: PC_MG *mg = (PC_MG *)pc->data;
418: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
420: PetscFunctionBegin;
421: PetscCheck(thr >= 0.0 && thr < 1.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Relative prolongator filter threshold %g must be in [0,1)", (double)thr);
422: if (thr > 0.2) PetscCall(PetscInfo(pc, "Warning: prolongator filter threshold %g is unusually large; typical values are 0.01 to 0.1\n", (double)thr));
423: pc_gamg->prolongator_filter = thr;
424: PetscFunctionReturn(PETSC_SUCCESS);
425: }
427: static PetscErrorCode PCGAMGGetProlongatorFilter_AGG(PC pc, PetscReal *thr)
428: {
429: PC_MG *mg = (PC_MG *)pc->data;
430: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
432: PetscFunctionBegin;
433: *thr = pc_gamg->prolongator_filter;
434: PetscFunctionReturn(PETSC_SUCCESS);
435: }
437: static PetscErrorCode PCGAMGSetProlongatorFilterScale_AGG(PC pc, PetscReal scale)
438: {
439: PC_MG *mg = (PC_MG *)pc->data;
440: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
442: PetscFunctionBegin;
443: PetscCheck(scale >= 0.0 && scale <= 1.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Prolongator filter scale %g must be in [0,1]", (double)scale);
444: pc_gamg->prolongator_filter_scale = scale;
445: PetscFunctionReturn(PETSC_SUCCESS);
446: }
448: static PetscErrorCode PCGAMGGetProlongatorFilterScale_AGG(PC pc, PetscReal *scale)
449: {
450: PC_MG *mg = (PC_MG *)pc->data;
451: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
453: PetscFunctionBegin;
454: *scale = pc_gamg->prolongator_filter_scale;
455: PetscFunctionReturn(PETSC_SUCCESS);
456: }
458: static PetscErrorCode PCSetFromOptions_GAMG_AGG(PC pc, PetscOptionItems PetscOptionsObject)
459: {
460: PC_MG *mg = (PC_MG *)pc->data;
461: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
462: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
463: PetscBool n_aggressive_flg, old_sq_provided = PETSC_FALSE, new_sq_provided = PETSC_FALSE, new_sqr_graph = pc_gamg_agg->use_aggressive_square_graph;
464: PetscInt nsq_graph_old = 0;
465: PetscReal thr = pc_gamg->prolongator_filter;
466: PetscReal scale = pc_gamg->prolongator_filter_scale;
467: PetscBool flg;
469: PetscFunctionBegin;
470: PetscOptionsHeadBegin(PetscOptionsObject, "GAMG-AGG options");
471: 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));
472: // aggressive coarsening logic with deprecated -pc_gamg_square_graph
473: 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));
474: if (!n_aggressive_flg)
475: 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));
476: 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));
477: if (!new_sq_provided && old_sq_provided) {
478: pc_gamg_agg->aggressive_coarsening_levels = nsq_graph_old; // could be zero
479: pc_gamg_agg->use_aggressive_square_graph = PETSC_TRUE;
480: }
481: if (new_sq_provided && old_sq_provided)
482: 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));
483: 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));
484: 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));
485: 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));
486: PetscCall(PetscOptionsBool("-pc_gamg_graph_symmetrize", "Symmetrize graph for coarsening", "PCGAMGSetGraphSymmetrize", pc_gamg_agg->graph_symmetrize, &pc_gamg_agg->graph_symmetrize, NULL));
487: PetscCall(PetscOptionsBoundedReal("-pc_gamg_prolongator_filter", "Relative Frobenius-norm threshold for block filtering of the prolongator (0=disabled)", "PCGAMGSetProlongatorFilter", thr, &thr, &flg, 0.0));
488: if (flg) PetscCall(PCGAMGSetProlongatorFilter(pc, thr));
489: PetscCall(PetscOptionsRangeReal("-pc_gamg_prolongator_filter_scale", "Per-level scaling of the prolongator filter threshold", "PCGAMGSetProlongatorFilterScale", scale, &scale, &flg, 0.0, 1.0));
490: if (flg) PetscCall(PCGAMGSetProlongatorFilterScale(pc, scale));
492: PetscOptionsHeadEnd();
493: PetscFunctionReturn(PETSC_SUCCESS);
494: }
496: static PetscErrorCode PCDestroy_GAMG_AGG(PC pc)
497: {
498: PC_MG *mg = (PC_MG *)pc->data;
499: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
500: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
502: PetscFunctionBegin;
503: PetscCall(MatCoarsenDestroy(&pc_gamg_agg->crs));
504: PetscCall(PetscFree(pc_gamg->subctx));
505: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", NULL));
506: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", NULL));
507: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetAggressive_C", NULL));
508: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetMinDegreeOrdering_C", NULL));
509: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetLowMemoryFilter_C", NULL));
510: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveSquareGraph_C", NULL));
511: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetGraphSymmetrize_C", NULL));
512: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetProlongatorFilter_C", NULL));
513: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGGetProlongatorFilter_C", NULL));
514: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetProlongatorFilterScale_C", NULL));
515: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGGetProlongatorFilterScale_C", NULL));
516: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
517: PetscFunctionReturn(PETSC_SUCCESS);
518: }
520: /*
521: PCSetCoordinates_AGG
523: Collective
525: Input Parameter:
526: . pc - the preconditioner context
527: . ndm - dimension of data (used for dof/vertex for Stokes)
528: . a_nloc - number of vertices local
529: . coords - [a_nloc][ndm] - interleaved coordinate data: {x_0, y_0, z_0, x_1, y_1, ...}
530: */
532: static PetscErrorCode PCSetCoordinates_AGG(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords)
533: {
534: PC_MG *mg = (PC_MG *)pc->data;
535: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
536: PetscInt arrsz, kk, ii, jj, nloc, ndatarows, ndf;
537: Mat mat = pc->pmat;
539: PetscFunctionBegin;
542: nloc = a_nloc;
544: /* SA: null space vectors */
545: PetscCall(MatGetBlockSize(mat, &ndf)); /* this does not work for Stokes */
546: if (coords && ndf == 1) pc_gamg->data_cell_cols = 1; /* scalar w/ coords and SA (not needed) */
547: else if (coords) {
548: PetscCheck(ndm <= ndf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "degrees of motion %" PetscInt_FMT " > block size %" PetscInt_FMT, ndm, ndf);
549: pc_gamg->data_cell_cols = (ndm == 2 ? 3 : 6); /* displacement elasticity */
550: 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);
551: } else pc_gamg->data_cell_cols = ndf; /* no data, force SA with constant null space vectors */
552: pc_gamg->data_cell_rows = ndatarows = ndf;
553: 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);
554: arrsz = nloc * pc_gamg->data_cell_rows * pc_gamg->data_cell_cols;
556: if (!pc_gamg->data || (pc_gamg->data_sz != arrsz)) {
557: PetscCall(PetscFree(pc_gamg->data));
558: PetscCall(PetscMalloc1(arrsz + 1, &pc_gamg->data));
559: }
560: /* copy data in - column-oriented */
561: for (kk = 0; kk < nloc; kk++) {
562: const PetscInt M = nloc * pc_gamg->data_cell_rows; /* stride into data */
563: PetscReal *data = &pc_gamg->data[kk * ndatarows]; /* start of cell */
565: if (pc_gamg->data_cell_cols == 1) *data = 1.0;
566: else {
567: /* translational modes */
568: for (ii = 0; ii < ndatarows; ii++) {
569: for (jj = 0; jj < ndatarows; jj++) {
570: if (ii == jj) data[ii * M + jj] = 1.0;
571: else data[ii * M + jj] = 0.0;
572: }
573: }
575: /* rotational modes */
576: if (coords) {
577: if (ndm == 2) {
578: data += 2 * M;
579: data[0] = -coords[2 * kk + 1];
580: data[1] = coords[2 * kk];
581: } else {
582: data += 3 * M;
583: data[0] = 0.0;
584: data[M + 0] = coords[3 * kk + 2];
585: data[2 * M + 0] = -coords[3 * kk + 1];
586: data[1] = -coords[3 * kk + 2];
587: data[M + 1] = 0.0;
588: data[2 * M + 1] = coords[3 * kk];
589: data[2] = coords[3 * kk + 1];
590: data[M + 2] = -coords[3 * kk];
591: data[2 * M + 2] = 0.0;
592: }
593: }
594: }
595: }
596: pc_gamg->data_sz = arrsz;
597: PetscFunctionReturn(PETSC_SUCCESS);
598: }
600: /*
601: PCSetData_AGG - called if data is not set with PCSetCoordinates.
602: Looks in Mat for near null space.
603: Does not work for Stokes
605: Input Parameter:
606: . pc -
607: . a_A - matrix to get (near) null space out of.
608: */
609: static PetscErrorCode PCSetData_AGG(PC pc, Mat a_A)
610: {
611: PC_MG *mg = (PC_MG *)pc->data;
612: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
613: MatNullSpace mnull;
615: PetscFunctionBegin;
616: PetscCall(MatGetNearNullSpace(a_A, &mnull));
617: if (!mnull) {
618: DM dm;
620: PetscCall(PCGetDM(pc, &dm));
621: if (!dm) PetscCall(MatGetDM(a_A, &dm));
622: if (dm) {
623: PetscObject deformation;
624: PetscInt Nf;
626: PetscCall(DMGetNumFields(dm, &Nf));
627: if (Nf) {
628: PetscCall(DMGetField(dm, 0, NULL, &deformation));
629: if (deformation) {
630: PetscCall(PetscObjectQuery(deformation, "nearnullspace", (PetscObject *)&mnull));
631: if (!mnull) PetscCall(PetscObjectQuery(deformation, "nullspace", (PetscObject *)&mnull));
632: }
633: }
634: }
635: }
637: if (!mnull) {
638: PetscInt bs, NN, MM;
640: PetscCall(MatGetBlockSize(a_A, &bs));
641: PetscCall(MatGetLocalSize(a_A, &MM, &NN));
642: PetscCheck(MM % bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MM %" PetscInt_FMT " must be divisible by bs %" PetscInt_FMT, MM, bs);
643: PetscCall(PCSetCoordinates_AGG(pc, bs, MM / bs, NULL));
644: } else {
645: PetscReal *nullvec;
646: PetscBool has_const;
647: PetscInt i, j, mlocal, nvec, bs;
648: const Vec *vecs;
649: const PetscScalar *v;
651: PetscCall(MatGetLocalSize(a_A, &mlocal, NULL));
652: PetscCall(MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs));
653: for (i = 0; i < nvec; i++) {
654: PetscCall(VecGetLocalSize(vecs[i], &j));
655: PetscCheck(j == mlocal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Attached null space vector size %" PetscInt_FMT " != matrix size %" PetscInt_FMT, j, mlocal);
656: }
657: pc_gamg->data_sz = (nvec + !!has_const) * mlocal;
658: PetscCall(PetscMalloc1((nvec + !!has_const) * mlocal, &nullvec));
659: if (has_const)
660: for (i = 0; i < mlocal; i++) nullvec[i] = 1.0;
661: for (i = 0; i < nvec; i++) {
662: PetscCall(VecGetArrayRead(vecs[i], &v));
663: for (j = 0; j < mlocal; j++) nullvec[(i + !!has_const) * mlocal + j] = PetscRealPart(v[j]);
664: PetscCall(VecRestoreArrayRead(vecs[i], &v));
665: }
666: pc_gamg->data = nullvec;
667: pc_gamg->data_cell_cols = (nvec + !!has_const);
668: PetscCall(MatGetBlockSize(a_A, &bs));
669: pc_gamg->data_cell_rows = bs;
670: }
671: PetscFunctionReturn(PETSC_SUCCESS);
672: }
674: /*
675: formProl0 - collect null space data for each aggregate, do QR, put R in coarse grid data and Q in P_0
677: Input Parameter:
678: . agg_llists - list of arrays with aggregates -- list from selected vertices of aggregate unselected vertices
679: . bs - row block size
680: . nSAvec - column bs of new P
681: . my0crs - global index of start of locals
682: . data_stride - bs*(nloc nodes + ghost nodes) [data_stride][nSAvec]
683: . data_in[data_stride*nSAvec] - local data on fine grid
684: . flid_fgid[data_stride/bs] - make local to global IDs, includes ghosts in 'locals_llist'
686: Output Parameter:
687: . a_data_out - in with fine grid data (w/ghosts), out with coarse grid data
688: . a_Prol - prolongation operator
689: */
690: 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)
691: {
692: PetscInt Istart, my0, Iend, nloc, clid, flid = 0, aggID, kk, jj, ii, mm, nSelected, minsz, nghosts, out_data_stride;
693: MPI_Comm comm;
694: PetscReal *out_data;
695: PetscCDIntNd *pos;
696: PetscHMapI fgid_flid;
698: PetscFunctionBegin;
699: PetscCall(PetscObjectGetComm((PetscObject)a_Prol, &comm));
700: PetscCall(MatGetOwnershipRange(a_Prol, &Istart, &Iend));
701: nloc = (Iend - Istart) / bs;
702: my0 = Istart / bs;
703: 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);
704: Iend /= bs;
705: nghosts = data_stride / bs - nloc;
707: PetscCall(PetscHMapICreateWithSize(2 * nghosts + 1, &fgid_flid));
709: for (kk = 0; kk < nghosts; kk++) PetscCall(PetscHMapISet(fgid_flid, flid_fgid[nloc + kk], nloc + kk));
711: /* count selected -- same as number of cols of P */
712: for (nSelected = mm = 0; mm < nloc; mm++) {
713: PetscBool ise;
715: PetscCall(PetscCDIsEmptyAt(agg_llists, mm, &ise));
716: if (!ise) nSelected++;
717: }
718: PetscCall(MatGetOwnershipRangeColumn(a_Prol, &ii, &jj));
719: PetscCheck((ii / nSAvec) == my0crs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ii %" PetscInt_FMT " /nSAvec %" PetscInt_FMT " != my0crs %" PetscInt_FMT, ii, nSAvec, my0crs);
720: 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);
722: /* aloc space for coarse point data (output) */
723: out_data_stride = nSelected * nSAvec;
725: PetscCall(PetscMalloc1(out_data_stride * nSAvec, &out_data));
726: for (ii = 0; ii < out_data_stride * nSAvec; ii++) out_data[ii] = PETSC_MAX_REAL;
727: *a_data_out = out_data; /* output - stride nSelected*nSAvec */
729: /* find points and set prolongation */
730: minsz = 100;
731: for (mm = clid = 0; mm < nloc; mm++) {
732: PetscCall(PetscCDCountAt(agg_llists, mm, &jj));
733: if (jj > 0) {
734: const PetscInt lid = mm, cgid = my0crs + clid;
735: PetscInt cids[100]; /* max bs */
736: PetscBLASInt asz, M, N;
737: PetscBLASInt Mdata, LDA, LWORK;
738: PetscScalar *qqc, *qqr, *TAU, *WORK;
739: PetscInt *fids;
740: PetscReal *data;
742: PetscCall(PetscBLASIntCast(jj, &asz));
743: PetscCall(PetscBLASIntCast(asz * bs, &M));
744: PetscCall(PetscBLASIntCast(nSAvec, &N));
745: PetscCall(PetscBLASIntCast(M + ((N - M > 0) ? N - M : 0), &Mdata));
746: PetscCall(PetscBLASIntCast(Mdata, &LDA));
747: PetscCall(PetscBLASIntCast(N * bs, &LWORK));
748: /* count agg */
749: if (asz < minsz) minsz = asz;
751: /* get block */
752: PetscCall(PetscMalloc5(Mdata * N, &qqc, M * N, &qqr, N, &TAU, LWORK, &WORK, M, &fids));
754: aggID = 0;
755: PetscCall(PetscCDGetHeadPos(agg_llists, lid, &pos));
756: while (pos) {
757: PetscInt gid1;
759: PetscCall(PetscCDIntNdGetID(pos, &gid1));
760: PetscCall(PetscCDGetNextPos(agg_llists, lid, &pos));
762: if (gid1 >= my0 && gid1 < Iend) flid = gid1 - my0;
763: else {
764: PetscCall(PetscHMapIGet(fgid_flid, gid1, &flid));
765: PetscCheck(flid >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cannot find gid1 in table");
766: }
767: /* copy in B_i matrix - column-oriented */
768: data = &data_in[flid * bs];
769: for (ii = 0; ii < bs; ii++) {
770: for (jj = 0; jj < N; jj++) {
771: PetscReal d = data[jj * data_stride + ii];
773: qqc[jj * Mdata + aggID * bs + ii] = d;
774: }
775: }
776: /* set fine IDs */
777: for (kk = 0; kk < bs; kk++) fids[aggID * bs + kk] = flid_fgid[flid] * bs + kk;
778: aggID++;
779: }
781: /* pad with zeros */
782: for (ii = asz * bs; ii < Mdata; ii++) {
783: for (jj = 0; jj < N; jj++, kk++) qqc[jj * Mdata + ii] = .0;
784: }
786: /* QR */
787: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
788: PetscCallLAPACKInfo("LAPACKgeqrf", LAPACKgeqrf_(&Mdata, &N, qqc, &LDA, TAU, WORK, &LWORK, &info));
789: PetscCall(PetscFPTrapPop());
790: /* get R - column-oriented - output B_{i+1} */
791: {
792: PetscReal *data = &out_data[clid * nSAvec];
794: for (jj = 0; jj < nSAvec; jj++) {
795: for (ii = 0; ii < nSAvec; ii++) {
796: 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);
797: if (ii <= jj) data[jj * out_data_stride + ii] = PetscRealPart(qqc[jj * Mdata + ii]);
798: else data[jj * out_data_stride + ii] = 0.;
799: }
800: }
801: }
803: /* get Q - row-oriented */
804: PetscCallLAPACKInfo("LAPACKorgqr", LAPACKorgqr_(&Mdata, &N, &N, qqc, &LDA, TAU, WORK, &LWORK, &info));
806: for (ii = 0; ii < M; ii++) {
807: for (jj = 0; jj < N; jj++) qqr[N * ii + jj] = qqc[jj * Mdata + ii];
808: }
810: /* add diagonal block of P0 */
811: for (kk = 0; kk < N; kk++) cids[kk] = N * cgid + kk; /* global col IDs in P0 */
812: PetscCall(MatSetValues(a_Prol, M, fids, N, cids, qqr, INSERT_VALUES));
813: PetscCall(PetscFree5(qqc, qqr, TAU, WORK, fids));
814: clid++;
815: } /* coarse agg */
816: } /* for all fine nodes */
817: PetscCall(MatAssemblyBegin(a_Prol, MAT_FINAL_ASSEMBLY));
818: PetscCall(MatAssemblyEnd(a_Prol, MAT_FINAL_ASSEMBLY));
819: PetscCall(PetscHMapIDestroy(&fgid_flid));
820: PetscFunctionReturn(PETSC_SUCCESS);
821: }
823: static PetscErrorCode PCView_GAMG_AGG(PC pc, PetscViewer viewer)
824: {
825: PC_MG *mg = (PC_MG *)pc->data;
826: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
827: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
829: PetscFunctionBegin;
830: PetscCall(PetscViewerASCIIPrintf(viewer, " AGG specific options\n"));
831: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of levels of aggressive coarsening %" PetscInt_FMT "\n", pc_gamg_agg->aggressive_coarsening_levels));
832: if (pc_gamg_agg->aggressive_coarsening_levels > 0) {
833: PetscCall(PetscViewerASCIIPrintf(viewer, " %s aggressive coarsening\n", !pc_gamg_agg->use_aggressive_square_graph ? "MIS-k" : "Square graph"));
834: 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));
835: }
836: PetscCall(PetscViewerASCIIPushTab(viewer));
837: PetscCall(PetscViewerASCIIPushTab(viewer));
838: PetscCall(PetscViewerASCIIPushTab(viewer));
839: PetscCall(PetscViewerASCIIPushTab(viewer));
840: if (pc_gamg_agg->crs) PetscCall(MatCoarsenView(pc_gamg_agg->crs, viewer));
841: else PetscCall(PetscViewerASCIIPrintf(viewer, "Coarsening algorithm not yet selected\n"));
842: PetscCall(PetscViewerASCIIPopTab(viewer));
843: PetscCall(PetscViewerASCIIPopTab(viewer));
844: PetscCall(PetscViewerASCIIPopTab(viewer));
845: PetscCall(PetscViewerASCIIPopTab(viewer));
846: PetscCall(PetscViewerASCIIPrintf(viewer, " Number smoothing steps to construct prolongation %" PetscInt_FMT "\n", pc_gamg_agg->nsmooths));
847: PetscFunctionReturn(PETSC_SUCCESS);
848: }
850: static PetscErrorCode PCGAMGCreateGraph_AGG(PC pc, Mat Amat, Mat *a_Gmat)
851: {
852: PC_MG *mg = (PC_MG *)pc->data;
853: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
854: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
855: const PetscReal vfilter = pc_gamg->threshold[pc_gamg->current_level];
856: PetscBool ishem, ismis;
857: const char *prefix;
858: MatInfo info0, info1;
859: PetscInt bs;
861: PetscFunctionBegin;
862: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
863: /* 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 */
864: /* MATCOARSENHEM requires numerical weights for edges so ensure they are computed */
865: PetscCall(MatCoarsenDestroy(&pc_gamg_agg->crs));
866: PetscCall(MatCoarsenCreate(PetscObjectComm((PetscObject)pc), &pc_gamg_agg->crs));
867: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix));
868: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc_gamg_agg->crs, prefix));
869: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)pc_gamg_agg->crs, "pc_gamg_"));
870: PetscCall(MatCoarsenSetFromOptions(pc_gamg_agg->crs));
871: PetscCall(MatGetBlockSize(Amat, &bs));
872: // check for valid indices wrt bs
873: for (int ii = 0; ii < pc_gamg_agg->crs->strength_index_size; ii++) {
874: 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",
875: pc_gamg_agg->crs->strength_index[ii], bs);
876: }
877: PetscCall(PetscObjectTypeCompare((PetscObject)pc_gamg_agg->crs, MATCOARSENHEM, &ishem));
878: if (ishem) {
879: 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));
880: pc_gamg_agg->aggressive_coarsening_levels = 0; // aggressive and HEM does not make sense
881: PetscCall(MatCoarsenSetMaximumIterations(pc_gamg_agg->crs, pc_gamg_agg->crs->max_it)); // for code coverage
882: PetscCall(MatCoarsenSetThreshold(pc_gamg_agg->crs, vfilter)); // for code coverage
883: } else {
884: PetscCall(PetscObjectTypeCompare((PetscObject)pc_gamg_agg->crs, MATCOARSENMIS, &ismis));
885: if (ismis && pc_gamg_agg->aggressive_coarsening_levels && !pc_gamg_agg->use_aggressive_square_graph) {
886: PetscCall(PetscInfo(pc, "MIS and aggressive coarsening and no square graph: force square graph\n"));
887: pc_gamg_agg->use_aggressive_square_graph = PETSC_TRUE;
888: }
889: }
890: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
891: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_GRAPH], 0, 0, 0, 0));
892: PetscCall(MatGetInfo(Amat, MAT_LOCAL, &info0)); /* global reduction */
894: if (ishem || pc_gamg_agg->use_low_mem_filter) {
895: 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));
896: } else {
897: // make scalar graph, symmetrize if not known to be symmetric, scale, but do not filter (expensive)
898: 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));
899: if (vfilter >= 0) {
900: PetscInt Istart, Iend, ncols, nnz0, nnz1, NN, MM, nloc;
901: Mat tGmat, Gmat = *a_Gmat;
902: MPI_Comm comm;
903: const PetscScalar *vals;
904: const PetscInt *idx;
905: PetscInt *d_nnz, *o_nnz, kk, *garray = NULL, *AJ, maxcols = 0;
906: MatScalar *AA; // this is checked in graph
907: PetscBool isseqaij;
908: Mat a, b, c;
909: MatType jtype;
911: PetscCall(PetscObjectGetComm((PetscObject)Gmat, &comm));
912: PetscCall(PetscObjectBaseTypeCompare((PetscObject)Gmat, MATSEQAIJ, &isseqaij));
913: PetscCall(MatGetType(Gmat, &jtype));
914: PetscCall(MatCreate(comm, &tGmat));
915: PetscCall(MatSetType(tGmat, jtype));
917: /* TODO GPU: this can be called when filter = 0 -> Probably provide MatAIJThresholdCompress that compresses the entries below a threshold?
918: Also, if the matrix is symmetric, can we skip this
919: operation? It can be very expensive on large matrices. */
921: // global sizes
922: PetscCall(MatGetSize(Gmat, &MM, &NN));
923: PetscCall(MatGetOwnershipRange(Gmat, &Istart, &Iend));
924: nloc = Iend - Istart;
925: PetscCall(PetscMalloc2(nloc, &d_nnz, nloc, &o_nnz));
926: if (isseqaij) {
927: a = Gmat;
928: b = NULL;
929: } else {
930: Mat_MPIAIJ *d = (Mat_MPIAIJ *)Gmat->data;
932: a = d->A;
933: b = d->B;
934: garray = d->garray;
935: }
936: /* Determine upper bound on non-zeros needed in new filtered matrix */
937: for (PetscInt row = 0; row < nloc; row++) {
938: PetscCall(MatGetRow(a, row, &ncols, NULL, NULL));
939: d_nnz[row] = ncols;
940: if (ncols > maxcols) maxcols = ncols;
941: PetscCall(MatRestoreRow(a, row, &ncols, NULL, NULL));
942: }
943: if (b) {
944: for (PetscInt row = 0; row < nloc; row++) {
945: PetscCall(MatGetRow(b, row, &ncols, NULL, NULL));
946: o_nnz[row] = ncols;
947: if (ncols > maxcols) maxcols = ncols;
948: PetscCall(MatRestoreRow(b, row, &ncols, NULL, NULL));
949: }
950: }
951: PetscCall(MatSetSizes(tGmat, nloc, nloc, MM, MM));
952: PetscCall(MatSetBlockSizes(tGmat, 1, 1));
953: PetscCall(MatSeqAIJSetPreallocation(tGmat, 0, d_nnz));
954: PetscCall(MatMPIAIJSetPreallocation(tGmat, 0, d_nnz, 0, o_nnz));
955: PetscCall(MatSetOption(tGmat, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
956: PetscCall(PetscFree2(d_nnz, o_nnz));
957: PetscCall(PetscMalloc2(maxcols, &AA, maxcols, &AJ));
958: nnz0 = nnz1 = 0;
959: for (c = a, kk = 0; c && kk < 2; c = b, kk++) {
960: for (PetscInt row = 0, grow = Istart, ncol_row, jj; row < nloc; row++, grow++) {
961: PetscCall(MatGetRow(c, row, &ncols, &idx, &vals));
962: for (ncol_row = jj = 0; jj < ncols; jj++, nnz0++) {
963: PetscScalar sv = PetscAbs(PetscRealPart(vals[jj]));
964: if (PetscRealPart(sv) > vfilter) {
965: PetscInt cid = idx[jj] + Istart; //diag
967: nnz1++;
968: if (c != a) cid = garray[idx[jj]];
969: AA[ncol_row] = vals[jj];
970: AJ[ncol_row] = cid;
971: ncol_row++;
972: }
973: }
974: PetscCall(MatRestoreRow(c, row, &ncols, &idx, &vals));
975: PetscCall(MatSetValues(tGmat, 1, &grow, ncol_row, AJ, AA, INSERT_VALUES));
976: }
977: }
978: PetscCall(PetscFree2(AA, AJ));
979: PetscCall(MatAssemblyBegin(tGmat, MAT_FINAL_ASSEMBLY));
980: PetscCall(MatAssemblyEnd(tGmat, MAT_FINAL_ASSEMBLY));
981: PetscCall(MatPropagateSymmetryOptions(Gmat, tGmat)); /* Normal Mat options are not relevant ? */
982: 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));
983: PetscCall(MatViewFromOptions(tGmat, NULL, "-mat_filter_graph_view"));
984: PetscCall(MatDestroy(&Gmat));
985: *a_Gmat = tGmat;
986: }
987: }
989: PetscCall(MatGetInfo(*a_Gmat, MAT_LOCAL, &info1)); /* global reduction */
990: 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));
991: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_GRAPH], 0, 0, 0, 0));
992: PetscFunctionReturn(PETSC_SUCCESS);
993: }
995: typedef PetscInt NState;
996: static const NState NOT_DONE = -2;
997: static const NState DELETED = -1;
998: static const NState REMOVED = -3;
999: #define IS_SELECTED(s) ((s) != DELETED && (s) != NOT_DONE && (s) != REMOVED)
1001: /*
1002: fixAggregatesWithSquare - greedy grab of with G1 (unsquared graph) -- AIJ specific -- change to fixAggregatesWithSquare -- TODD
1003: - AGG-MG specific: clears singletons out of 'selected_2'
1005: Input Parameter:
1006: . Gmat_2 - global matrix of squared graph (data not defined)
1007: . Gmat_1 - base graph to grab with base graph
1008: Input/Output Parameter:
1009: . aggs_2 - linked list of aggs with gids)
1010: */
1011: static PetscErrorCode fixAggregatesWithSquare(PC pc, Mat Gmat_2, Mat Gmat_1, PetscCoarsenData *aggs_2)
1012: {
1013: PetscBool isMPI;
1014: Mat_SeqAIJ *matA_1, *matB_1 = NULL;
1015: MPI_Comm comm;
1016: PetscInt lid, *ii, *idx, ix, Iend, my0, n, j;
1017: Mat_MPIAIJ *mpimat_2 = NULL, *mpimat_1 = NULL;
1018: const PetscInt nloc = Gmat_2->rmap->n;
1019: PetscInt *lid_cprowID_1 = NULL, *lid_parent_gid, *cpcol_2_par_orig = NULL;
1020: NState *lid_state, *cpcol_1_state = NULL, *cpcol_2_state = NULL;
1021: PetscSF sf_1 = NULL, sf_2 = NULL;
1022: PetscMPIInt rank;
1024: PetscFunctionBegin;
1025: PetscCall(PetscObjectGetComm((PetscObject)Gmat_2, &comm));
1026: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1027: PetscCall(MatGetOwnershipRange(Gmat_1, &my0, &Iend));
1029: /* get submatrices */
1030: PetscCall(PetscStrbeginswith(((PetscObject)Gmat_1)->type_name, MATMPIAIJ, &isMPI));
1031: PetscCall(PetscInfo(pc, "isMPI = %s\n", isMPI ? "yes" : "no"));
1032: PetscCall(PetscMalloc3(nloc, &lid_state, nloc, &lid_parent_gid, nloc, &lid_cprowID_1));
1033: for (lid = 0; lid < nloc; lid++) lid_cprowID_1[lid] = -1;
1034: if (isMPI) {
1035: /* grab matrix objects */
1036: mpimat_2 = (Mat_MPIAIJ *)Gmat_2->data;
1037: mpimat_1 = (Mat_MPIAIJ *)Gmat_1->data;
1038: matA_1 = (Mat_SeqAIJ *)mpimat_1->A->data;
1039: matB_1 = (Mat_SeqAIJ *)mpimat_1->B->data;
1041: /* force compressed row storage for B matrix in AuxMat */
1042: PetscCall(MatCheckCompressedRow(mpimat_1->B, matB_1->nonzerorowcnt, &matB_1->compressedrow, matB_1->i, Gmat_1->rmap->n, -1.0));
1043: for (ix = 0; ix < matB_1->compressedrow.nrows; ix++) {
1044: PetscInt lid = matB_1->compressedrow.rindex[ix];
1046: PetscCheck(lid <= nloc && lid >= -1, PETSC_COMM_SELF, PETSC_ERR_USER, "lid %" PetscInt_FMT " out of range. nloc = %" PetscInt_FMT, lid, nloc);
1047: if (lid != -1) lid_cprowID_1[lid] = ix;
1048: }
1049: } else {
1050: PetscBool isAIJ;
1052: PetscCall(PetscStrbeginswith(((PetscObject)Gmat_1)->type_name, MATSEQAIJ, &isAIJ));
1053: PetscCheck(isAIJ, PETSC_COMM_SELF, PETSC_ERR_USER, "Require AIJ matrix.");
1054: matA_1 = (Mat_SeqAIJ *)Gmat_1->data;
1055: }
1056: if (nloc > 0) PetscCheck(!matB_1 || matB_1->compressedrow.use, PETSC_COMM_SELF, PETSC_ERR_PLIB, "matB_1 && !matB_1->compressedrow.use: PETSc bug???");
1057: /* get state of locals and selected gid for deleted */
1058: for (lid = 0; lid < nloc; lid++) {
1059: lid_parent_gid[lid] = -1;
1060: lid_state[lid] = DELETED;
1061: }
1063: /* set lid_state */
1064: for (lid = 0; lid < nloc; lid++) {
1065: PetscCDIntNd *pos;
1067: PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
1068: if (pos) {
1069: PetscInt gid1;
1071: PetscCall(PetscCDIntNdGetID(pos, &gid1));
1072: PetscCheck(gid1 == lid + my0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "gid1 %" PetscInt_FMT " != lid %" PetscInt_FMT " + my0 %" PetscInt_FMT, gid1, lid, my0);
1073: lid_state[lid] = gid1;
1074: }
1075: }
1077: /* map local to selected local, DELETED means a ghost owns it */
1078: for (lid = 0; lid < nloc; lid++) {
1079: NState state = lid_state[lid];
1081: if (IS_SELECTED(state)) {
1082: PetscCDIntNd *pos;
1084: PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
1085: while (pos) {
1086: PetscInt gid1;
1088: PetscCall(PetscCDIntNdGetID(pos, &gid1));
1089: PetscCall(PetscCDGetNextPos(aggs_2, lid, &pos));
1090: if (gid1 >= my0 && gid1 < Iend) lid_parent_gid[gid1 - my0] = lid + my0;
1091: }
1092: }
1093: }
1094: /* get 'cpcol_1/2_state' & cpcol_2_par_orig */
1095: if (isMPI) {
1096: PetscInt nghost_1, nghost_2;
1098: PetscCall(VecGetLocalSize(mpimat_1->lvec, &nghost_1));
1099: PetscCall(VecGetLocalSize(mpimat_2->lvec, &nghost_2));
1100: PetscCall(PetscMalloc3(nghost_1, &cpcol_1_state, nghost_2, &cpcol_2_state, nghost_2, &cpcol_2_par_orig));
1101: PetscCall(MatGetMultPetscSF(Gmat_1, &sf_1));
1102: PetscCall(MatGetMultPetscSF(Gmat_2, &sf_2));
1103: /* get 'cpcol_1_state' */
1104: PetscCall(PetscSFBcastBegin(sf_1, MPIU_INT, lid_state, cpcol_1_state, MPI_REPLACE));
1105: PetscCall(PetscSFBcastEnd(sf_1, MPIU_INT, lid_state, cpcol_1_state, MPI_REPLACE));
1106: /* get 'cpcol_2_state' */
1107: PetscCall(PetscSFBcastBegin(sf_2, MPIU_INT, lid_state, cpcol_2_state, MPI_REPLACE));
1108: PetscCall(PetscSFBcastEnd(sf_2, MPIU_INT, lid_state, cpcol_2_state, MPI_REPLACE));
1109: /* get 'cpcol_2_par_orig' */
1110: PetscCall(PetscSFBcastBegin(sf_2, MPIU_INT, lid_parent_gid, cpcol_2_par_orig, MPI_REPLACE));
1111: PetscCall(PetscSFBcastEnd(sf_2, MPIU_INT, lid_parent_gid, cpcol_2_par_orig, MPI_REPLACE));
1112: } /* ismpi */
1113: for (lid = 0; lid < nloc; lid++) {
1114: NState state = lid_state[lid];
1116: if (IS_SELECTED(state)) {
1117: /* steal locals */
1118: ii = matA_1->i;
1119: n = ii[lid + 1] - ii[lid];
1120: idx = matA_1->j + ii[lid];
1121: for (j = 0; j < n; j++) {
1122: PetscInt lidj = idx[j], sgid;
1123: NState statej = lid_state[lidj];
1125: if (statej == DELETED && (sgid = lid_parent_gid[lidj]) != lid + my0) { /* steal local */
1126: lid_parent_gid[lidj] = lid + my0; /* send this if sgid is not local */
1127: if (sgid >= my0 && sgid < Iend) { /* I'm stealing this local from a local sgid */
1128: PetscInt hav = 0, slid = sgid - my0, gidj = lidj + my0;
1129: PetscCDIntNd *pos, *last = NULL;
1131: /* looking for local from local so id_llist_2 works */
1132: PetscCall(PetscCDGetHeadPos(aggs_2, slid, &pos));
1133: while (pos) {
1134: PetscInt gid;
1136: PetscCall(PetscCDIntNdGetID(pos, &gid));
1137: if (gid == gidj) {
1138: PetscCheck(last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "last cannot be null");
1139: PetscCall(PetscCDRemoveNextNode(aggs_2, slid, last));
1140: PetscCall(PetscCDAppendNode(aggs_2, lid, pos));
1141: hav = 1;
1142: break;
1143: } else last = pos;
1144: PetscCall(PetscCDGetNextPos(aggs_2, slid, &pos));
1145: }
1146: if (hav != 1) {
1147: PetscCheck(hav, PETSC_COMM_SELF, PETSC_ERR_PLIB, "failed to find adj in 'selected' lists - structurally unsymmetric matrix");
1148: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "found node %" PetscInt_FMT " times???", hav);
1149: }
1150: } else { /* I'm stealing this local, owned by a ghost */
1151: 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.",
1152: ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "");
1153: PetscCall(PetscCDAppendID(aggs_2, lid, lidj + my0));
1154: }
1155: }
1156: } /* local neighbors */
1157: } else if (state == DELETED /* && lid_cprowID_1 */) {
1158: PetscInt sgidold = lid_parent_gid[lid];
1160: /* see if I have a selected ghost neighbor that will steal me */
1161: if ((ix = lid_cprowID_1[lid]) != -1) {
1162: ii = matB_1->compressedrow.i;
1163: n = ii[ix + 1] - ii[ix];
1164: idx = matB_1->j + ii[ix];
1165: for (j = 0; j < n; j++) {
1166: PetscInt cpid = idx[j];
1167: NState statej = cpcol_1_state[cpid];
1169: if (IS_SELECTED(statej) && sgidold != statej) { /* ghost will steal this, remove from my list */
1170: lid_parent_gid[lid] = statej; /* send who selected */
1171: if (sgidold >= my0 && sgidold < Iend) { /* this was mine */
1172: PetscInt hav = 0, oldslidj = sgidold - my0;
1173: PetscCDIntNd *pos, *last = NULL;
1175: /* remove from 'oldslidj' list */
1176: PetscCall(PetscCDGetHeadPos(aggs_2, oldslidj, &pos));
1177: while (pos) {
1178: PetscInt gid;
1180: PetscCall(PetscCDIntNdGetID(pos, &gid));
1181: if (lid + my0 == gid) {
1182: /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove lid from oldslidj list *\/ */
1183: PetscCheck(last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "last cannot be null");
1184: PetscCall(PetscCDRemoveNextNode(aggs_2, oldslidj, last));
1185: /* ghost statej will add this later */
1186: hav = 1;
1187: break;
1188: } else last = pos;
1189: PetscCall(PetscCDGetNextPos(aggs_2, oldslidj, &pos));
1190: }
1191: if (hav != 1) {
1192: PetscCheck(hav, PETSC_COMM_SELF, PETSC_ERR_PLIB, "failed to find (hav=%" PetscInt_FMT ") adj in 'selected' lists - structurally unsymmetric matrix", hav);
1193: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "found node %" PetscInt_FMT " times???", hav);
1194: }
1195: } else {
1196: /* TODO: ghosts remove this later */
1197: }
1198: }
1199: }
1200: }
1201: } /* selected/deleted */
1202: } /* node loop */
1204: if (isMPI) {
1205: PetscInt *cpcol_2_parent;
1206: const PetscInt *cpcol_2_gid = mpimat_2->garray;
1207: PetscInt cpid, nghost_2;
1208: PetscHMapI gid_cpid;
1210: PetscCall(VecGetLocalSize(mpimat_2->lvec, &nghost_2));
1211: PetscCall(PetscMalloc1(nghost_2, &cpcol_2_parent));
1213: /* get 'cpcol_2_parent' */
1214: PetscCall(PetscSFBcastBegin(sf_2, MPIU_INT, lid_parent_gid, cpcol_2_parent, MPI_REPLACE));
1215: PetscCall(PetscSFBcastEnd(sf_2, MPIU_INT, lid_parent_gid, cpcol_2_parent, MPI_REPLACE));
1217: /* look for deleted ghosts and add to table */
1218: PetscCall(PetscHMapICreateWithSize(2 * nghost_2 + 1, &gid_cpid));
1219: for (cpid = 0; cpid < nghost_2; cpid++) {
1220: NState state = cpcol_2_state[cpid];
1222: if (state == DELETED) {
1223: PetscInt sgid_new = cpcol_2_parent[cpid];
1224: PetscInt sgid_old = cpcol_2_par_orig[cpid];
1226: if (sgid_old == -1 && sgid_new != -1) {
1227: PetscInt gid = cpcol_2_gid[cpid];
1229: PetscCall(PetscHMapISet(gid_cpid, gid, cpid));
1230: }
1231: }
1232: }
1234: /* look for deleted ghosts and see if they moved - remove it */
1235: for (lid = 0; lid < nloc; lid++) {
1236: NState state = lid_state[lid];
1238: if (IS_SELECTED(state)) {
1239: PetscCDIntNd *pos, *last = NULL;
1241: /* look for deleted ghosts and see if they moved */
1242: PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
1243: while (pos) {
1244: PetscInt gid;
1246: PetscCall(PetscCDIntNdGetID(pos, &gid));
1247: if (gid < my0 || gid >= Iend) {
1248: PetscCall(PetscHMapIGet(gid_cpid, gid, &cpid));
1249: if (cpid != -1) {
1250: /* a moved ghost - */
1251: /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove 'flid' from list *\/ */
1252: PetscCall(PetscCDRemoveNextNode(aggs_2, lid, last));
1253: } else last = pos;
1254: } else last = pos;
1256: PetscCall(PetscCDGetNextPos(aggs_2, lid, &pos));
1257: } /* loop over list of deleted */
1258: } /* selected */
1259: }
1260: PetscCall(PetscHMapIDestroy(&gid_cpid));
1262: /* look at ghosts, see if they changed - and it */
1263: for (cpid = 0; cpid < nghost_2; cpid++) {
1264: PetscInt sgid_new = cpcol_2_parent[cpid];
1266: if (sgid_new >= my0 && sgid_new < Iend) { /* this is mine */
1267: PetscInt gid = cpcol_2_gid[cpid];
1268: PetscInt slid_new = sgid_new - my0, hav = 0;
1269: PetscCDIntNd *pos;
1271: /* search for this gid to see if I have it */
1272: PetscCall(PetscCDGetHeadPos(aggs_2, slid_new, &pos));
1273: while (pos) {
1274: PetscInt gidj;
1276: PetscCall(PetscCDIntNdGetID(pos, &gidj));
1277: PetscCall(PetscCDGetNextPos(aggs_2, slid_new, &pos));
1279: if (gidj == gid) {
1280: hav = 1;
1281: break;
1282: }
1283: }
1284: if (hav != 1) {
1285: /* insert 'flidj' into head of llist */
1286: PetscCall(PetscCDAppendID(aggs_2, slid_new, gid));
1287: }
1288: }
1289: }
1290: PetscCall(PetscFree(cpcol_2_parent));
1291: PetscCall(PetscFree3(cpcol_1_state, cpcol_2_state, cpcol_2_par_orig));
1292: }
1293: PetscCall(PetscFree3(lid_state, lid_parent_gid, lid_cprowID_1));
1294: PetscFunctionReturn(PETSC_SUCCESS);
1295: }
1297: /*
1298: PCGAMGCoarsen_AGG - supports squaring the graph (deprecated) and new graph for
1299: communication of QR data used with HEM and MISk coarsening
1301: Input Parameter:
1302: . a_pc - this
1304: Input/Output Parameter:
1305: . a_Gmat1 - graph to coarsen (in), graph off processor edges for QR gather scatter (out)
1307: Output Parameter:
1308: . agg_lists - list of aggregates
1310: */
1311: static PetscErrorCode PCGAMGCoarsen_AGG(PC a_pc, Mat *a_Gmat1, PetscCoarsenData **agg_lists)
1312: {
1313: PC_MG *mg = (PC_MG *)a_pc->data;
1314: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1315: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
1316: Mat Gmat2, Gmat1 = *a_Gmat1; /* aggressive graph */
1317: IS perm;
1318: PetscInt Istart, Iend, Ii, nloc, bs, nn;
1319: PetscInt *permute, *degree;
1320: PetscBool *bIndexSet;
1321: PetscReal hashfact;
1322: PetscInt iSwapIndex;
1323: PetscRandom random;
1324: MPI_Comm comm;
1326: PetscFunctionBegin;
1327: PetscCall(PetscObjectGetComm((PetscObject)Gmat1, &comm));
1328: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
1329: PetscCall(MatGetLocalSize(Gmat1, &nn, NULL));
1330: PetscCall(MatGetBlockSize(Gmat1, &bs));
1331: PetscCheck(bs == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "bs %" PetscInt_FMT " must be 1", bs);
1332: nloc = nn / bs;
1333: /* get MIS aggs - randomize */
1334: PetscCall(PetscMalloc2(nloc, &permute, nloc, °ree));
1335: PetscCall(PetscCalloc1(nloc, &bIndexSet));
1336: for (Ii = 0; Ii < nloc; Ii++) permute[Ii] = Ii;
1337: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &random));
1338: PetscCall(MatGetOwnershipRange(Gmat1, &Istart, &Iend));
1339: for (Ii = 0; Ii < nloc; Ii++) {
1340: PetscInt nc;
1342: PetscCall(MatGetRow(Gmat1, Istart + Ii, &nc, NULL, NULL));
1343: degree[Ii] = nc;
1344: PetscCall(MatRestoreRow(Gmat1, Istart + Ii, &nc, NULL, NULL));
1345: }
1346: for (Ii = 0; Ii < nloc; Ii++) {
1347: PetscCall(PetscRandomGetValueReal(random, &hashfact));
1348: iSwapIndex = (PetscInt)(hashfact * nloc) % nloc;
1349: if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) {
1350: PetscInt iTemp = permute[iSwapIndex];
1352: permute[iSwapIndex] = permute[Ii];
1353: permute[Ii] = iTemp;
1354: iTemp = degree[iSwapIndex];
1355: degree[iSwapIndex] = degree[Ii];
1356: degree[Ii] = iTemp;
1357: bIndexSet[iSwapIndex] = PETSC_TRUE;
1358: }
1359: }
1360: // apply minimum degree ordering -- NEW
1361: if (pc_gamg_agg->use_minimum_degree_ordering) PetscCall(PetscSortIntWithArray(nloc, degree, permute));
1362: PetscCall(PetscFree(bIndexSet));
1363: PetscCall(PetscRandomDestroy(&random));
1364: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_USE_POINTER, &perm));
1365: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0));
1366: // square graph
1367: 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));
1368: else Gmat2 = Gmat1;
1369: // switch to old MIS-1 for square graph
1370: if (pc_gamg->current_level < pc_gamg_agg->aggressive_coarsening_levels) {
1371: if (!pc_gamg_agg->use_aggressive_square_graph) PetscCall(MatCoarsenMISKSetDistance(pc_gamg_agg->crs, pc_gamg_agg->aggressive_mis_k)); // hardwire to MIS-2
1372: else PetscCall(MatCoarsenSetType(pc_gamg_agg->crs, MATCOARSENMIS)); // old MIS -- side effect
1373: } else if (pc_gamg_agg->use_aggressive_square_graph && pc_gamg_agg->aggressive_coarsening_levels > 0) { // we reset the MIS
1374: const char *prefix;
1376: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)a_pc, &prefix));
1377: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc_gamg_agg->crs, prefix));
1378: PetscCall(MatCoarsenSetFromOptions(pc_gamg_agg->crs)); // get the default back on non-aggressive levels when square graph switched to old MIS
1379: }
1380: PetscCall(MatCoarsenSetAdjacency(pc_gamg_agg->crs, Gmat2));
1381: PetscCall(MatCoarsenSetStrictAggs(pc_gamg_agg->crs, PETSC_TRUE));
1382: PetscCall(MatCoarsenSetGreedyOrdering(pc_gamg_agg->crs, perm));
1383: PetscCall(MatCoarsenApply(pc_gamg_agg->crs));
1384: PetscCall(MatCoarsenGetData(pc_gamg_agg->crs, agg_lists)); /* output */
1386: PetscCall(ISDestroy(&perm));
1387: PetscCall(PetscFree2(permute, degree));
1388: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0));
1390: if (Gmat2 != Gmat1) { // square graph, we need ghosts for selected
1391: PetscCoarsenData *llist = *agg_lists;
1393: PetscCall(fixAggregatesWithSquare(a_pc, Gmat2, Gmat1, *agg_lists));
1394: PetscCall(MatDestroy(&Gmat1));
1395: *a_Gmat1 = Gmat2; /* output */
1396: PetscCall(PetscCDSetMat(llist, *a_Gmat1)); /* Need a graph with ghosts here */
1397: }
1398: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
1399: PetscFunctionReturn(PETSC_SUCCESS);
1400: }
1402: /*
1403: PCGAMGConstructProlongator_AGG
1405: Input Parameter:
1406: . pc - this
1407: . Amat - matrix on this fine level
1408: . Graph - used to get ghost data for nodes in
1409: . agg_lists - list of aggregates
1410: Output Parameter:
1411: . a_P_out - prolongation operator to the next level
1412: */
1413: static PetscErrorCode PCGAMGConstructProlongator_AGG(PC pc, Mat Amat, PetscCoarsenData *agg_lists, Mat *a_P_out)
1414: {
1415: PC_MG *mg = (PC_MG *)pc->data;
1416: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1417: const PetscInt col_bs = pc_gamg->data_cell_cols;
1418: PetscInt Istart, Iend, nloc, ii, jj, kk, my0, nLocalSelected, bs;
1419: Mat Gmat, Prol;
1420: PetscMPIInt size;
1421: MPI_Comm comm;
1422: PetscReal *data_w_ghost;
1423: PetscInt myCrs0, nbnodes = 0, *flid_fgid;
1424: MatType mtype;
1426: PetscFunctionBegin;
1427: PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
1428: PetscCheck(col_bs >= 1, comm, PETSC_ERR_PLIB, "Column bs cannot be less than 1");
1429: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1430: PetscCallMPI(MPI_Comm_size(comm, &size));
1431: PetscCall(MatGetOwnershipRange(Amat, &Istart, &Iend));
1432: PetscCall(MatGetBlockSize(Amat, &bs));
1433: nloc = (Iend - Istart) / bs;
1434: my0 = Istart / bs;
1435: 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);
1436: PetscCall(PetscCDGetMat(agg_lists, &Gmat)); // get auxiliary matrix for ghost edges for size > 1
1438: /* get 'nLocalSelected' */
1439: for (ii = 0, nLocalSelected = 0; ii < nloc; ii++) {
1440: PetscBool ise;
1442: /* filter out singletons 0 or 1? */
1443: PetscCall(PetscCDIsEmptyAt(agg_lists, ii, &ise));
1444: if (!ise) nLocalSelected++;
1445: }
1447: /* create prolongator, create P matrix */
1448: PetscCall(MatGetType(Amat, &mtype));
1449: PetscCall(MatCreate(comm, &Prol));
1450: PetscCall(MatSetSizes(Prol, nloc * bs, nLocalSelected * col_bs, PETSC_DETERMINE, PETSC_DETERMINE));
1451: PetscCall(MatSetBlockSizes(Prol, bs, col_bs)); // should this be before MatSetSizes?
1452: PetscCall(MatSetType(Prol, mtype));
1453: #if PetscDefined(HAVE_DEVICE)
1454: PetscBool flg;
1455: PetscCall(MatBoundToCPU(Amat, &flg));
1456: PetscCall(MatBindToCPU(Prol, flg));
1457: if (flg) PetscCall(MatSetBindingPropagates(Prol, PETSC_TRUE));
1458: #endif
1459: PetscCall(MatSeqAIJSetPreallocation(Prol, col_bs, NULL));
1460: PetscCall(MatMPIAIJSetPreallocation(Prol, col_bs, NULL, col_bs, NULL));
1462: /* can get all points "removed" */
1463: PetscCall(MatGetSize(Prol, &kk, &ii));
1464: if (!ii) {
1465: PetscCall(PetscInfo(pc, "%s: No selected points on coarse grid\n", ((PetscObject)pc)->prefix));
1466: PetscCall(MatDestroy(&Prol));
1467: *a_P_out = NULL; /* out */
1468: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1469: PetscFunctionReturn(PETSC_SUCCESS);
1470: }
1471: PetscCall(PetscInfo(pc, "%s: New grid %" PetscInt_FMT " nodes\n", ((PetscObject)pc)->prefix, ii / col_bs));
1472: PetscCall(MatGetOwnershipRangeColumn(Prol, &myCrs0, &kk));
1474: 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);
1475: myCrs0 = myCrs0 / col_bs;
1476: 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);
1478: /* create global vector of data in 'data_w_ghost' */
1479: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0));
1480: if (size > 1) { /* get ghost null space data */
1481: PetscReal *tmp_gdata, *tmp_ldata, *tp2;
1483: PetscCall(PetscMalloc1(nloc, &tmp_ldata));
1484: for (jj = 0; jj < col_bs; jj++) {
1485: for (kk = 0; kk < bs; kk++) {
1486: PetscInt stride;
1487: const PetscReal *tp = PetscSafePointerPlusOffset(pc_gamg->data, jj * bs * nloc + kk);
1489: for (PetscInt ii = 0; ii < nloc; ii++, tp += bs) tmp_ldata[ii] = *tp;
1491: PetscCall(PCGAMGGetDataWithGhosts(Gmat, 1, tmp_ldata, &stride, &tmp_gdata));
1493: if (!jj && !kk) { /* now I know how many total nodes - allocate TODO: move below and do in one 'col_bs' call */
1494: PetscCall(PetscMalloc1(stride * bs * col_bs, &data_w_ghost));
1495: nbnodes = bs * stride;
1496: }
1497: tp2 = PetscSafePointerPlusOffset(data_w_ghost, jj * bs * stride + kk);
1498: for (PetscInt ii = 0; ii < stride; ii++, tp2 += bs) *tp2 = tmp_gdata[ii];
1499: PetscCall(PetscFree(tmp_gdata));
1500: }
1501: }
1502: PetscCall(PetscFree(tmp_ldata));
1503: } else {
1504: nbnodes = bs * nloc;
1505: data_w_ghost = pc_gamg->data;
1506: }
1508: /* get 'flid_fgid' TODO - move up to get 'stride' and do get null space data above in one step (jj loop) */
1509: if (size > 1) {
1510: Mat_MPIAIJ *mpimat = (Mat_MPIAIJ *)Gmat->data;
1511: PetscInt stride, num_ghosts;
1513: /* global IDs of the local vertices, then of the ghosts, which are given by garray */
1514: PetscCall(VecGetLocalSize(mpimat->lvec, &num_ghosts));
1515: stride = nloc + num_ghosts;
1516: PetscCall(PetscMalloc1(stride, &flid_fgid));
1517: for (kk = 0; kk < nloc; kk++) flid_fgid[kk] = my0 + kk;
1518: if (num_ghosts > 0) PetscCall(PetscArraycpy(flid_fgid + nloc, mpimat->garray, num_ghosts));
1520: PetscCheck(stride == nbnodes / bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "stride %" PetscInt_FMT " != nbnodes %" PetscInt_FMT "/bs %" PetscInt_FMT, stride, nbnodes, bs);
1521: } else {
1522: PetscCall(PetscMalloc1(nloc, &flid_fgid));
1523: for (kk = 0; kk < nloc; kk++) flid_fgid[kk] = my0 + kk;
1524: }
1525: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0));
1526: /* get P0 */
1527: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0));
1528: {
1529: PetscReal *data_out = NULL;
1531: PetscCall(formProl0(agg_lists, bs, col_bs, myCrs0, nbnodes, data_w_ghost, flid_fgid, &data_out, Prol));
1532: PetscCall(PetscFree(pc_gamg->data));
1534: pc_gamg->data = data_out;
1535: pc_gamg->data_cell_rows = col_bs;
1536: pc_gamg->data_sz = col_bs * col_bs * nLocalSelected;
1537: }
1538: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0));
1539: if (size > 1) PetscCall(PetscFree(data_w_ghost));
1540: PetscCall(PetscFree(flid_fgid));
1542: *a_P_out = Prol; /* out */
1543: PetscCall(MatViewFromOptions(Prol, NULL, "-pc_gamg_agg_view_initial_prolongation"));
1545: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1546: PetscFunctionReturn(PETSC_SUCCESS);
1547: }
1549: // Drop small node-coupling blocks of the prolongator; see the PCGAMGSetProlongatorFilter() manual page
1550: static PetscErrorCode PCGAMGProlongatorBlockFilter_AGG(PC pc, Mat Prol, PetscInt col_bs, PetscReal thr)
1551: {
1552: PetscInt rbs, cbs, rStart, rEnd, nfn, local_nnz = 0, max_fn_cols = 0, max_row_cols = 0, ndrows = 0, zoff = 0;
1553: PetscInt *cn_gid, *fn_col, *fn_slot, *fn_rcnt, *zcols, *drow, *doff, *dcnt;
1554: PetscReal *cn_n2;
1555: PetscReal thr2 = thr * thr;
1556: PetscScalar *zeros;
1557: PetscBool no_off_proc;
1559: PetscFunctionBegin;
1560: PetscCall(MatGetBlockSizes(Prol, &rbs, &cbs));
1561: /* Prol is built internally with MatSetBlockSizes(..., col_bs); a mismatch here means smoothing corrupted the block structure */
1562: PetscCheck(cbs == col_bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Prolongator column block size %" PetscInt_FMT " != nSAvec %" PetscInt_FMT " (block structure lost during smoothing; should be unreachable with a user-facing config)", cbs, col_bs);
1563: PetscCall(MatGetOwnershipRange(Prol, &rStart, &rEnd));
1564: PetscCheck((rEnd - rStart) % rbs == 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Local rows %" PetscInt_FMT " not divisible by row block size %" PetscInt_FMT, rEnd - rStart, rbs);
1565: nfn = (rEnd - rStart) / rbs;
1567: /* Pre-pass over the row widths to size the scratch tightly: the total local nonzeros bound
1568: the dropped columns (zcols), the widest fine-node block row bounds the distinct coarse
1569: nodes of one fine node (cn_gid/cn_n2), and the widest single row bounds the entries
1570: zeroed by one MatSetValues() call (zeros) */
1571: for (PetscInt fn = 0; fn < nfn; fn++) {
1572: PetscInt fn_cols = 0;
1574: for (PetscInt rr = 0; rr < rbs; rr++) {
1575: PetscInt grow = rStart + fn * rbs + rr, ncols;
1577: PetscCall(MatGetRow(Prol, grow, &ncols, NULL, NULL));
1578: fn_cols += ncols;
1579: if (ncols > max_row_cols) max_row_cols = ncols;
1580: PetscCall(MatRestoreRow(Prol, grow, &ncols, NULL, NULL));
1581: }
1582: local_nnz += fn_cols;
1583: if (fn_cols > max_fn_cols) max_fn_cols = fn_cols;
1584: }
1586: PetscCall(PetscMalloc5(max_fn_cols, &cn_gid, max_fn_cols, &cn_n2, max_fn_cols, &fn_col, max_fn_cols, &fn_slot, rbs, &fn_rcnt));
1587: PetscCall(PetscCalloc1(max_row_cols, &zeros));
1588: PetscCall(PetscMalloc1(local_nnz, &zcols));
1589: PetscCall(PetscMalloc3(rEnd - rStart, &drow, rEnd - rStart, &doff, rEnd - rStart, &dcnt));
1591: for (PetscInt fn = 0; fn < nfn; fn++) {
1592: PetscInt ncn = 0; /* distinct coarse nodes touched by this fine node-block */
1593: PetscInt nent = 0; /* entries of this fine node-block cached in fn_col/fn_slot */
1594: PetscReal maxn2 = 0.0;
1596: /* Pass A: accumulate the per-coarse-node block Frobenius norm^2 over the rbs rows, caching each
1597: entry's column and resolved coarse-node slot so that Pass B needs no second MatGetRow() sweep */
1598: for (PetscInt rr = 0; rr < rbs; rr++) {
1599: PetscInt grow = rStart + fn * rbs + rr, ncols;
1600: const PetscInt *cols;
1601: const PetscScalar *vals;
1603: PetscCall(MatGetRow(Prol, grow, &ncols, &cols, &vals));
1604: fn_rcnt[rr] = ncols;
1605: for (PetscInt k = 0; k < ncols; k++) {
1606: PetscInt cn = cols[k] / col_bs, s = -1;
1607: PetscReal av = PetscAbsScalar(vals[k]);
1609: /* Linear dedup scan: ncn is bounded by the coarse-node degree of this
1610: fine node (typically ~5-20 for AMG-coarsened elasticity), so O(ncn)
1611: per entry is faster than a hash map at this scale */
1612: for (PetscInt t = 0; t < ncn; t++)
1613: if (cn_gid[t] == cn) {
1614: s = t;
1615: break;
1616: }
1617: if (s < 0) {
1618: s = ncn++;
1619: cn_gid[s] = cn;
1620: cn_n2[s] = 0.0;
1621: }
1622: cn_n2[s] += av * av;
1623: fn_col[nent] = cols[k];
1624: fn_slot[nent] = s;
1625: nent++;
1626: }
1627: PetscCall(MatRestoreRow(Prol, grow, &ncols, &cols, &vals));
1628: }
1629: for (PetscInt t = 0; t < ncn; t++)
1630: if (cn_n2[t] > maxn2) maxn2 = cn_n2[t];
1632: /* Pass B: collect the entries of blocks below thr*max for zeroing. The test is strict, so the
1633: strongest block of the node survives for any thr <= 1 */
1634: nent = 0;
1635: for (PetscInt rr = 0; rr < rbs; rr++) {
1636: PetscInt rec_off = zoff;
1638: for (PetscInt k = 0; k < fn_rcnt[rr]; k++, nent++)
1639: if (cn_n2[fn_slot[nent]] < thr2 * maxn2) zcols[zoff++] = fn_col[nent];
1640: if (zoff > rec_off) {
1641: drow[ndrows] = rStart + fn * rbs + rr;
1642: doff[ndrows] = rec_off;
1643: dcnt[ndrows] = zoff - rec_off;
1644: ndrows++;
1645: }
1646: }
1647: }
1649: /* all insertions are in local rows; skip the off-process assembly communication */
1650: PetscCall(MatGetOption(Prol, MAT_NO_OFF_PROC_ENTRIES, &no_off_proc));
1651: PetscCall(MatSetOption(Prol, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
1652: for (PetscInt i = 0; i < ndrows; i++) PetscCall(MatSetValues(Prol, 1, &drow[i], dcnt[i], &zcols[doff[i]], zeros, INSERT_VALUES));
1653: PetscCall(MatAssemblyBegin(Prol, MAT_FINAL_ASSEMBLY));
1654: PetscCall(MatAssemblyEnd(Prol, MAT_FINAL_ASSEMBLY));
1655: PetscCall(MatSetOption(Prol, MAT_NO_OFF_PROC_ENTRIES, no_off_proc));
1657: PetscCall(PetscFree5(cn_gid, cn_n2, fn_col, fn_slot, fn_rcnt));
1658: PetscCall(PetscFree(zeros));
1659: PetscCall(PetscFree(zcols));
1660: PetscCall(PetscFree3(drow, doff, dcnt));
1662: /* Compress out the explicit zeros just set, so the filter actually sparsifies. keep must be
1663: PETSC_FALSE: with keep, a zero whose local column index equals its local row index survives (an
1664: index-based diagonal test that is meaningless for the rectangular Prol). The compression is done
1665: in place, deliberately: a MatDuplicate()/MatHeaderReplace() copy as in MatFilter() would return
1666: the freed CSR tail to the allocator, at the cost of a peak-memory spike. MatEliminateZeros() is not
1667: implemented by all matrix types; in those cases the zeros are left in the sparsity pattern;
1668: the step-3 correction in PCGAMGKernelPreservingFilter_AGG() skips exactly-zero entries, so a dropped
1669: block stays dropped either way, it just still costs storage here. */
1670: if (Prol->ops->eliminatezeros) PetscCall(MatEliminateZeros(Prol, PETSC_FALSE));
1671: else PetscCall(PetscInfo(pc, "PCGAMGProlongatorBlockFilter_AGG: skipping zero elimination for %s; filtered entries are zeroed but not removed\n", ((PetscObject)Prol)->type_name));
1672: PetscFunctionReturn(PETSC_SUCCESS);
1673: }
1675: // Filter the prolongator by fine-node/coarse-node coupling blocks, then restore the near-null space constraint P*B_c = B; see the PCGAMGSetProlongatorFilter() manual page
1676: static PetscErrorCode PCGAMGKernelPreservingFilter_AGG(PC pc, Mat Prol, PetscReal threshold)
1677: {
1678: PC_MG *mg = (PC_MG *)pc->data;
1679: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1680: const PetscInt nSAvec = pc_gamg->data_cell_rows; /* == data_cell_cols after formProl0 */
1681: PetscInt cStart, cEnd, rStart, rEnd;
1682: const PetscReal *Bc_data = pc_gamg->data;
1683: Vec *Bc_vecs, *B_vecs;
1684: PetscScalar *Bc_arr;
1686: PetscFunctionBegin;
1687: PetscCall(PetscInfo(pc, "Kernel-preserving filter of prolongator with threshold %g, nSAvec=%" PetscInt_FMT "\n", (double)threshold, nSAvec));
1689: PetscCall(MatGetOwnershipRange(Prol, &rStart, &rEnd));
1690: PetscCall(MatGetOwnershipRangeColumn(Prol, &cStart, &cEnd));
1692: /* Step 1: build coarse null-space vectors and compute B = P_original * B_c */
1693: PetscCall(PetscMalloc1(nSAvec, &Bc_vecs));
1694: PetscCall(PetscMalloc1(nSAvec, &B_vecs));
1696: {
1697: PetscInt nloc = cEnd - cStart;
1698: for (PetscInt k = 0; k < nSAvec; k++) {
1699: PetscCall(MatCreateVecs(Prol, &Bc_vecs[k], &B_vecs[k]));
1700: /* fill local entries: Bc_data layout is Bc_data[k * nloc + c] (stride == nloc) */
1701: PetscCall(VecGetArray(Bc_vecs[k], &Bc_arr));
1702: for (PetscInt c = 0; c < nloc; c++) Bc_arr[c] = (PetscScalar)Bc_data[k * nloc + c];
1703: PetscCall(VecRestoreArray(Bc_vecs[k], &Bc_arr));
1704: PetscCall(MatMult(Prol, Bc_vecs[k], B_vecs[k]));
1705: }
1706: }
1708: /* Step 2: apply the block-aware threshold filter (drops whole node-coupling
1709: blocks, preserving the coarse-node block structure; see
1710: PCGAMGProlongatorBlockFilter_AGG()) */
1711: {
1712: PetscBool info_active = PETSC_FALSE;
1713: MatInfo info0, info1;
1714: PetscCall(PetscInfoEnabled(((PetscObject)pc)->classid, &info_active));
1715: if (info_active) PetscCall(MatGetInfo(Prol, MAT_GLOBAL_SUM, &info0));
1716: PetscCall(PCGAMGProlongatorBlockFilter_AGG(pc, Prol, nSAvec, threshold));
1717: if (info_active) {
1718: PetscCall(MatGetInfo(Prol, MAT_GLOBAL_SUM, &info1));
1719: 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));
1720: }
1721: }
1723: /* Step 3: correct rows to restore P_filtered * B_c = B */
1724: if (nSAvec == 1) {
1725: /*
1726: Scalar case: use `MatMult()` + element-wise scaling + `MatDiagonalScale()`.
1727: scale_i = B_i / (P_filtered * Bc)_i, then P_new = diag(scale) * P_filtered.
1728: A row is left unscaled, that is, the near-null space constraint is not enforced
1729: for it, in the two cases where scaling would do more harm than the constraint is
1730: worth: when |scale_i| would exceed smax, which covers a denominator that is zero
1731: (an empty row after the filter) or tiny relative to B_i, and when B_i is zero,
1732: where enforcing the constraint would zero the whole row of the prolongator.
1733: Zeros left in the sparsity pattern stay zero under `MatDiagonalScale()`.
1734: No ghost column access needed.
1735: */
1736: Vec d_vec, scale_vec;
1737: PetscInt n_local, n_unscaled = 0;
1738: PetscReal smax = 1.0e4; /* cap on the row scale; the filter drops small blocks, so a healthy row has scale ~1 */
1739: PetscScalar *s_arr;
1740: const PetscScalar *b_arr, *d_arr;
1742: PetscCall(MatCreateVecs(Prol, NULL, &d_vec));
1743: PetscCall(MatMult(Prol, Bc_vecs[0], d_vec));
1744: PetscCall(VecDuplicate(d_vec, &scale_vec));
1745: PetscCall(VecGetLocalSize(d_vec, &n_local));
1746: PetscCall(VecGetArrayRead(B_vecs[0], &b_arr));
1747: PetscCall(VecGetArrayRead(d_vec, &d_arr));
1748: PetscCall(VecGetArray(scale_vec, &s_arr));
1749: for (PetscInt i = 0; i < n_local; i++) {
1750: PetscReal b = PetscAbsScalar(b_arr[i]), d = PetscAbsScalar(d_arr[i]);
1752: if (b > 0.0 && smax * d > b) s_arr[i] = b_arr[i] / d_arr[i];
1753: else {
1754: s_arr[i] = 1.0;
1755: n_unscaled++;
1756: }
1757: }
1758: if (n_unscaled > 0) PetscCall(PetscInfo(pc, "PCGAMGKernelPreservingFilter_AGG: %" PetscInt_FMT " rows left unscaled (zero target or row scale above %g)\n", n_unscaled, (double)smax));
1759: PetscCall(VecRestoreArray(scale_vec, &s_arr));
1760: PetscCall(VecRestoreArrayRead(d_vec, &d_arr));
1761: PetscCall(VecRestoreArrayRead(B_vecs[0], &b_arr));
1762: PetscCall(MatDiagonalScale(Prol, scale_vec, NULL));
1763: PetscCall(VecDestroy(&d_vec));
1764: PetscCall(VecDestroy(&scale_vec));
1765: } else {
1766: /*
1767: Vector case (nSAvec > 1): per-row minimum-norm correction (Gram-matrix solve).
1768: Scatter Bc_data to include ghost column values using Prol's Mvctx,
1769: then build a hash map from global ghost column index to local ghost index
1770: so that `MatGetRow()` global column indices can be mapped to the ghosted array.
1771: */
1772: PetscInt nloc = cEnd - cStart;
1773: PetscInt ghost_stride;
1774: PetscReal *Bc_ghosted = NULL;
1775: const PetscReal *Bc_ghosted_ro;
1776: PetscMPIInt comm_size;
1777: PetscHMapI ghost_gid_to_lid; /* global ghost col index -> local ghost index (0-based) */
1778: PetscInt num_ghosts = 0;
1779: PetscBool no_off_proc;
1781: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)Prol), &comm_size));
1782: if (comm_size > 1) {
1783: Mat_MPIAIJ *mpimat = (Mat_MPIAIJ *)Prol->data;
1784: Vec tmp_vec;
1785: PetscScalar *data_arr;
1786: PetscInt nnodes;
1787: PetscBool isaij;
1789: PetscCall(PetscObjectBaseTypeCompare((PetscObject)Prol, MATMPIAIJ, &isaij));
1790: PetscCheck(isaij, PetscObjectComm((PetscObject)Prol), PETSC_ERR_SUP, "Prolongator filter requires an MPIAIJ-based prolongator, not %s", ((PetscObject)Prol)->type_name);
1791: PetscCall(VecGetLocalSize(mpimat->lvec, &num_ghosts));
1792: nnodes = nloc + num_ghosts;
1793: ghost_stride = nnodes;
1794: /*
1795: Scatter Bc_data to include ghost column values using Prol's Mvctx.
1796: Cannot use PCGAMGGetDataWithGhosts() because it assumes square matrix
1797: (uses MatGetOwnershipRange() for row indices, but Prol is rectangular).
1798: */
1799: PetscCall(MatCreateVecs(Prol, &tmp_vec, NULL));
1800: PetscCall(PetscMalloc1(nSAvec * nnodes, &Bc_ghosted));
1801: for (PetscInt dir = 0; dir < nSAvec; dir++) {
1802: PetscScalar *tmp_arr;
1803: PetscCall(VecGetArray(tmp_vec, &tmp_arr));
1804: for (PetscInt kk = 0; kk < nloc; kk++) {
1805: PetscReal val = Bc_data[dir * nloc + kk];
1806: Bc_ghosted[dir * nnodes + kk] = val;
1807: tmp_arr[kk] = (PetscScalar)val;
1808: }
1809: PetscCall(VecRestoreArray(tmp_vec, &tmp_arr));
1810: PetscCall(VecScatterBegin(mpimat->Mvctx, tmp_vec, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD));
1811: PetscCall(VecScatterEnd(mpimat->Mvctx, tmp_vec, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD));
1812: PetscCall(VecGetArray(mpimat->lvec, &data_arr));
1813: for (PetscInt g = 0; g < num_ghosts; g++) Bc_ghosted[dir * nnodes + nloc + g] = PetscRealPart(data_arr[g]);
1814: PetscCall(VecRestoreArray(mpimat->lvec, &data_arr));
1815: }
1816: PetscCall(VecDestroy(&tmp_vec));
1817: Bc_ghosted_ro = Bc_ghosted;
1818: /* build hash: global ghost col index -> local ghost index (0-based into ghost portion) */
1819: PetscCall(PetscHMapICreateWithSize(2 * num_ghosts + 1, &ghost_gid_to_lid));
1820: for (PetscInt g = 0; g < num_ghosts; g++) PetscCall(PetscHMapISet(ghost_gid_to_lid, mpimat->garray[g], g));
1821: } else {
1822: /* sequential: no ghosts, ghost_stride == nloc, use Bc_data directly (read-only) */
1823: ghost_stride = nloc;
1824: Bc_ghosted_ro = Bc_data;
1825: PetscCall(PetscHMapICreateWithSize(1, &ghost_gid_to_lid));
1826: }
1828: {
1829: PetscInt nrows = rEnd - rStart, max_ncols = 0;
1830: const PetscScalar **B_arrays;
1831: PetscScalar *work, *new_vals, *G, *rhs, *x, *bc_col;
1832: PetscReal *dscale;
1833: PetscInt *ghosted_idx, *act, *col_buf;
1834: PetscBLASInt *ipiv;
1835: PetscBLASInt N_b;
1837: PetscCall(PetscMalloc1(nSAvec, &B_arrays));
1838: for (PetscInt k = 0; k < nSAvec; k++) PetscCall(VecGetArrayRead(B_vecs[k], &B_arrays[k]));
1839: /* work: nSAvec*nSAvec Gram + nSAvec rhs + nSAvec solution + nSAvec bc_col scratch */
1840: PetscCall(PetscMalloc1(nSAvec * nSAvec + 3 * nSAvec, &work));
1841: PetscCall(PetscMalloc1(nSAvec, &dscale));
1842: PetscCall(PetscMalloc1(nSAvec, &ipiv));
1843: PetscCall(PetscBLASIntCast(nSAvec, &N_b));
1844: G = work;
1845: rhs = work + nSAvec * nSAvec;
1846: x = rhs + nSAvec;
1847: bc_col = x + nSAvec;
1849: /* find max row width and total nnz for pre-allocation */
1850: {
1851: PetscInt total_nnz = 0;
1852: for (PetscInt row = 0; row < nrows; row++) {
1853: PetscInt ncols;
1854: PetscCall(MatGetRow(Prol, rStart + row, &ncols, NULL, NULL));
1855: if (ncols > max_ncols) max_ncols = ncols;
1856: total_nnz += ncols;
1857: PetscCall(MatRestoreRow(Prol, rStart + row, &ncols, NULL, NULL));
1858: }
1859: /* allocate flat CSR-like buffers to store all corrections before applying */
1860: PetscCall(PetscMalloc1(total_nnz, &new_vals));
1861: PetscCall(PetscMalloc1(total_nnz, &col_buf));
1862: }
1863: PetscCall(PetscMalloc2(max_ncols, &ghosted_idx, max_ncols, &act));
1865: /* Pass 1: read rows, compute corrections, store in flat buffers */
1866: {
1867: PetscInt *row_offsets;
1868: PetscInt offset = 0, n_singular = 0, n_zero_rows = 0, n_corrected = 0, n_underdetermined = 0;
1869: PetscReal max_ynorm = 0.0;
1871: PetscCall(PetscMalloc1(nrows + 1, &row_offsets));
1872: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1874: for (PetscInt row = 0; row < nrows; row++) {
1875: PetscInt ncols, nact = 0, grow = rStart + row, roff = offset;
1876: PetscBLASInt NRHS = 1, LDA = N_b, LDB = N_b, info;
1877: const PetscInt *cols;
1878: const PetscScalar *vals;
1880: row_offsets[row] = roff;
1881: PetscCall(MatGetRow(Prol, grow, &ncols, &cols, &vals));
1882: /* Save the row unchanged, then correct only its nonzero entries: the block filter zeroes
1883: dropped entries but they are not removed from the sparsity pattern for every matrix type
1884: (see PCGAMGProlongatorBlockFilter_AGG()), and adding the correction to them would turn a
1885: dropped block back into a nonzero one */
1886: for (PetscInt j = 0; j < ncols; j++) {
1887: col_buf[roff + j] = cols[j];
1888: new_vals[roff + j] = vals[j];
1889: if (vals[j] != 0.0) act[nact++] = j;
1890: }
1891: offset = roff + ncols;
1892: if (nact == 0) {
1893: n_zero_rows++;
1894: PetscCall(MatRestoreRow(Prol, grow, &ncols, &cols, &vals));
1895: continue;
1896: }
1898: /* When nact < nSAvec the Gram matrix G is rank-deficient by construction;
1899: skip correction for this row (keep filtered values as-is).
1900: Note: the near-null space constraint P*Bc = B is NOT enforced for these rows.
1901: This typically occurs at boundary or isolated nodes where few coarse neighbors
1902: remain after filtering; the impact on convergence is generally small. */
1903: if (nact < nSAvec) {
1904: n_underdetermined++;
1905: PetscCall(MatRestoreRow(Prol, grow, &ncols, &cols, &vals));
1906: continue;
1907: }
1909: /* map the global column indices of the surviving entries to ghosted array indices */
1910: for (PetscInt a = 0; a < nact; a++) {
1911: PetscInt col = cols[act[a]];
1912: if (col >= cStart && col < cEnd) ghosted_idx[a] = col - cStart;
1913: else {
1914: PetscInt g = -1;
1915: PetscCall(PetscHMapIGet(ghost_gid_to_lid, col, &g));
1916: PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Off-diagonal column %" PetscInt_FMT " not found in ghost map for prolongator filter", col);
1917: ghosted_idx[a] = nloc + g;
1918: }
1919: }
1921: for (PetscInt i = 0; i < nSAvec * nSAvec; i++) G[i] = 0.0;
1923: /* rhs[k] = B[row,k] - sum_a P[row,act[a]] * Bc[ghosted_idx[a], k] */
1924: for (PetscInt k = 0; k < nSAvec; k++) {
1925: PetscScalar dot = 0.0;
1926: for (PetscInt a = 0; a < nact; a++) dot += vals[act[a]] * (PetscScalar)Bc_ghosted_ro[k * ghost_stride + ghosted_idx[a]];
1927: rhs[k] = B_arrays[k][row] - dot;
1928: }
1930: /* G[k1,k2] = sum_a Bc[a,k1] * Bc[a,k2] using pre-gathered bc_col */
1931: for (PetscInt a = 0; a < nact; a++) {
1932: PetscInt gidx = ghosted_idx[a];
1933: for (PetscInt k = 0; k < nSAvec; k++) bc_col[k] = (PetscScalar)Bc_ghosted_ro[k * ghost_stride + gidx];
1934: for (PetscInt k1 = 0; k1 < nSAvec; k1++)
1935: for (PetscInt k2 = k1; k2 < nSAvec; k2++) G[k1 * nSAvec + k2] += bc_col[k1] * bc_col[k2];
1936: }
1937: /* fill lower triangle from upper (G is symmetric) */
1938: for (PetscInt k1 = 1; k1 < nSAvec; k1++)
1939: for (PetscInt k2 = 0; k2 < k1; k2++) G[k1 * nSAvec + k2] = G[k2 * nSAvec + k1];
1941: /* Symmetric (Jacobi) equilibration: G is severely ill-conditioned for
1942: elasticity because the near-null modes have disparate scales (O(1)
1943: translations vs rotations that scale with the coordinates). Scale by
1944: dscale[k] = 1/sqrt(G[k,k]) so the rescaled Gram has a unit diagonal,
1945: then solve (D G D) y = D rhs and recover x = D y. This is identical to
1946: solving G x = rhs in exact arithmetic but removes the mode-scale
1947: ill-conditioning, making the correction accurate and FP-robust. */
1948: for (PetscInt k = 0; k < nSAvec; k++) {
1949: PetscReal gkk = PetscRealPart(G[k * nSAvec + k]);
1950: dscale[k] = gkk > 0.0 ? 1.0 / PetscSqrtReal(gkk) : 1.0;
1951: }
1952: for (PetscInt k1 = 0; k1 < nSAvec; k1++)
1953: for (PetscInt k2 = 0; k2 < nSAvec; k2++) G[k1 * nSAvec + k2] *= dscale[k1] * dscale[k2];
1955: /* solve (D G D) y = D rhs, then x = D y */
1956: for (PetscInt k = 0; k < nSAvec; k++) x[k] = dscale[k] * rhs[k];
1957: PetscCallBLAS("LAPACKgesv", LAPACKgesv_(&N_b, &NRHS, G, &LDA, ipiv, x, &LDB, &info));
1958: PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "LAPACKgesv: %" PetscBLASInt_FMT "-th argument had an illegal value", -info);
1959: if (info > 0) {
1960: /* G is singular despite nact >= nSAvec (Bc columns linearly dependent);
1961: keep filtered values as-is (near-null space constraint not enforced for this row) */
1962: n_singular++;
1963: PetscCall(MatRestoreRow(Prol, grow, &ncols, &cols, &vals));
1964: continue;
1965: }
1966: /* track the equilibrated solution norm ||y||^2 (x still holds y here): a
1967: basis-independent health metric, large values flag a problematic row */
1968: {
1969: PetscReal ynorm2 = 0.0;
1970: for (PetscInt k = 0; k < nSAvec; k++) ynorm2 += PetscSqr(PetscAbsScalar(x[k]));
1971: if (ynorm2 > max_ynorm) max_ynorm = ynorm2;
1972: }
1973: for (PetscInt k = 0; k < nSAvec; k++) x[k] *= dscale[k]; /* recover x = D y */
1974: n_corrected++;
1976: /* new_vals[act[a]] = vals[act[a]] + sum_k Bc[ghosted_idx[a],k] * x[k] */
1977: for (PetscInt a = 0; a < nact; a++) {
1978: PetscScalar delta = 0.0;
1979: PetscInt gidx = ghosted_idx[a];
1980: for (PetscInt k = 0; k < nSAvec; k++) delta += (PetscScalar)Bc_ghosted_ro[k * ghost_stride + gidx] * x[k];
1981: new_vals[roff + act[a]] += delta;
1982: }
1983: PetscCall(MatRestoreRow(Prol, grow, &ncols, &cols, &vals));
1984: }
1985: row_offsets[nrows] = offset;
1986: PetscCall(PetscFPTrapPop());
1987: PetscCall(PetscInfo(pc, "PCGAMGKernelPreservingFilter_AGG: corrected %" PetscInt_FMT "/%" PetscInt_FMT " rows, max equilibrated correction ||y||^2=%g\n", n_corrected, nrows, (double)max_ynorm));
1988: if (n_zero_rows + n_underdetermined + n_singular > 0)
1989: PetscCall(PetscInfo(pc, "PCGAMGKernelPreservingFilter_AGG: %" PetscInt_FMT " rows left uncorrected (zero=%" PetscInt_FMT " underdetermined=%" PetscInt_FMT " singular_G=%" PetscInt_FMT ")\n", n_zero_rows + n_underdetermined + n_singular, n_zero_rows, n_underdetermined, n_singular));
1991: /* Pass 2: apply all corrections at once */
1992: for (PetscInt row = 0; row < nrows; row++) {
1993: PetscInt grow = rStart + row;
1994: PetscInt nc = row_offsets[row + 1] - row_offsets[row];
1995: if (nc > 0) PetscCall(MatSetValues(Prol, 1, &grow, nc, col_buf + row_offsets[row], new_vals + row_offsets[row], INSERT_VALUES));
1996: }
1997: PetscCall(PetscFree(row_offsets));
1998: }
2000: for (PetscInt k = 0; k < nSAvec; k++) PetscCall(VecRestoreArrayRead(B_vecs[k], &B_arrays[k]));
2001: PetscCall(PetscFree(B_arrays));
2002: PetscCall(PetscFree(work));
2003: PetscCall(PetscFree(dscale));
2004: PetscCall(PetscFree(ipiv));
2005: PetscCall(PetscFree2(ghosted_idx, act));
2006: PetscCall(PetscFree(new_vals));
2007: PetscCall(PetscFree(col_buf));
2008: }
2010: PetscCall(PetscHMapIDestroy(&ghost_gid_to_lid));
2011: if (comm_size > 1) PetscCall(PetscFree(Bc_ghosted));
2013: /* all insertions are in local rows; skip the off-process assembly communication. The scalar
2014: branch needs no assembly at all: MatDiagonalScale() requires and preserves assembly */
2015: PetscCall(MatGetOption(Prol, MAT_NO_OFF_PROC_ENTRIES, &no_off_proc));
2016: PetscCall(MatSetOption(Prol, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
2017: PetscCall(MatAssemblyBegin(Prol, MAT_FINAL_ASSEMBLY));
2018: PetscCall(MatAssemblyEnd(Prol, MAT_FINAL_ASSEMBLY));
2019: PetscCall(MatSetOption(Prol, MAT_NO_OFF_PROC_ENTRIES, no_off_proc));
2020: }
2022: for (PetscInt k = 0; k < nSAvec; k++) {
2023: PetscCall(VecDestroy(&Bc_vecs[k]));
2024: PetscCall(VecDestroy(&B_vecs[k]));
2025: }
2026: PetscCall(PetscFree(Bc_vecs));
2027: PetscCall(PetscFree(B_vecs));
2028: PetscFunctionReturn(PETSC_SUCCESS);
2029: }
2031: /*
2032: PCGAMGOptimizeProlongator_AGG - given the initial prolongator optimizes it by smoothed aggregation pc_gamg_agg->nsmooths times
2034: Input Parameter:
2035: . pc - this
2036: . Amat - matrix on this fine level
2037: In/Output Parameter:
2038: . a_P - prolongation operator to the next level
2039: */
2040: static PetscErrorCode PCGAMGOptimizeProlongator_AGG(PC pc, Mat Amat, Mat *a_P)
2041: {
2042: PC_MG *mg = (PC_MG *)pc->data;
2043: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
2044: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
2045: Mat Prol = *a_P;
2046: MPI_Comm comm;
2047: KSP eksp;
2048: Vec bb, xx;
2049: PC epc;
2050: PetscReal alpha, emax, emin;
2051: PetscReal pfilter = pc_gamg->prolongator_filter * PetscPowRealInt(pc_gamg->prolongator_filter_scale, pc_gamg->current_level);
2053: PetscFunctionBegin;
2054: PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
2055: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0));
2057: /* compute maximum singular value of operator to be used in smoother */
2058: if (0 < pc_gamg_agg->nsmooths) {
2059: /* get eigen estimates */
2060: if (pc_gamg->emax > 0) {
2061: emin = pc_gamg->emin;
2062: emax = pc_gamg->emax;
2063: } else {
2064: const char *prefix;
2066: PetscCall(MatCreateVecs(Amat, &bb, NULL));
2067: PetscCall(MatCreateVecs(Amat, &xx, NULL));
2068: PetscCall(KSPSetNoisy_Private(Amat, bb));
2070: PetscCall(KSPCreate(comm, &eksp));
2071: PetscCall(KSPSetNestLevel(eksp, pc->kspnestlevel));
2072: PetscCall(PCGetOptionsPrefix(pc, &prefix));
2073: PetscCall(KSPSetOptionsPrefix(eksp, prefix));
2074: PetscCall(KSPAppendOptionsPrefix(eksp, "pc_gamg_esteig_"));
2075: {
2076: PetscBool isset, sflg;
2078: PetscCall(MatIsSPDKnown(Amat, &isset, &sflg));
2079: if (isset && sflg) PetscCall(KSPSetType(eksp, KSPCG));
2080: }
2081: PetscCall(KSPSetErrorIfNotConverged(eksp, pc->erroriffailure));
2082: PetscCall(KSPSetNormType(eksp, KSP_NORM_NONE));
2084: PetscCall(KSPSetInitialGuessNonzero(eksp, PETSC_FALSE));
2085: PetscCall(KSPSetOperators(eksp, Amat, Amat));
2087: PetscCall(KSPGetPC(eksp, &epc));
2088: PetscCall(PCSetType(epc, PCJACOBI)); /* smoother in smoothed agg. */
2090: 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
2092: PetscCall(KSPSetFromOptions(eksp));
2093: PetscCall(KSPSetComputeSingularValues(eksp, PETSC_TRUE));
2094: PetscCall(KSPSolve(eksp, bb, xx));
2095: PetscCall(KSPCheckSolve(eksp, pc, xx));
2097: PetscCall(KSPComputeExtremeSingularValues(eksp, &emax, &emin));
2098: PetscCall(PetscInfo(pc, "%s: Smooth P0: max eigen=%e min=%e PC=%s\n", ((PetscObject)pc)->prefix, (double)emax, (double)emin, PCJACOBI));
2099: PetscCall(VecDestroy(&xx));
2100: PetscCall(VecDestroy(&bb));
2101: PetscCall(KSPDestroy(&eksp));
2102: }
2103: if (pc_gamg->use_sa_esteig) {
2104: mg->min_eigen_DinvA[pc_gamg->current_level] = emin;
2105: mg->max_eigen_DinvA[pc_gamg->current_level] = emax;
2106: 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));
2107: } else {
2108: mg->min_eigen_DinvA[pc_gamg->current_level] = 0;
2109: mg->max_eigen_DinvA[pc_gamg->current_level] = 0;
2110: }
2111: } else {
2112: mg->min_eigen_DinvA[pc_gamg->current_level] = 0;
2113: mg->max_eigen_DinvA[pc_gamg->current_level] = 0;
2114: }
2116: /* smooth P0 */
2117: if (pc_gamg_agg->nsmooths > 0) {
2118: Vec diag;
2120: /* TODO: Set a PCFailedReason and exit the building of the AMG preconditioner */
2121: PetscCheck(emax != 0.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Computed maximum singular value as zero");
2123: PetscCall(MatCreateVecs(Amat, &diag, NULL));
2124: PetscCall(MatGetDiagonal(Amat, diag)); /* effectively PCJACOBI */
2125: PetscCall(VecReciprocal(diag));
2127: for (PetscInt jj = 0; jj < pc_gamg_agg->nsmooths; jj++) {
2128: Mat tMat;
2130: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0));
2131: /*
2132: Smooth aggregation on the prolongator
2134: P_{i} := (I - 1.4/emax D^{-1}A) P_i\{i-1}
2135: */
2136: PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0));
2137: PetscCall(MatMatMult(Amat, Prol, MAT_INITIAL_MATRIX, PETSC_CURRENT, &tMat));
2138: PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0));
2139: PetscCall(MatProductClear(tMat));
2140: PetscCall(MatDiagonalScale(tMat, diag, NULL));
2142: /* TODO: Document the 1.4 and don't hardwire it in this routine */
2143: alpha = -1.4 / emax;
2144: PetscCall(MatAYPX(tMat, alpha, Prol, SUBSET_NONZERO_PATTERN));
2145: PetscCall(MatDestroy(&Prol));
2146: Prol = tMat;
2147: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0));
2148: }
2149: PetscCall(VecDestroy(&diag));
2150: }
2151: /* a per-level threshold of 0 (from prolongator_filter_scale == 0 on the coarser levels) drops
2152: nothing, so skip the whole filter rather than pay for its passes and per-row solves */
2153: if (pfilter > 0.0) {
2154: PetscCall(PetscInfo(pc, "%s: level %" PetscInt_FMT " prolongator filter threshold %g (base %g, scale %g^%" PetscInt_FMT ")\n", ((PetscObject)pc)->prefix, pc_gamg->current_level, (double)pfilter, (double)pc_gamg->prolongator_filter,
2155: (double)pc_gamg->prolongator_filter_scale, pc_gamg->current_level));
2156: PetscCall(PCGAMGKernelPreservingFilter_AGG(pc, Prol, pfilter));
2157: }
2158: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0));
2159: PetscCall(MatViewFromOptions(Prol, NULL, "-pc_gamg_agg_view_prolongation"));
2160: *a_P = Prol;
2161: PetscFunctionReturn(PETSC_SUCCESS);
2162: }
2164: /*MC
2165: PCGAMGAGG - Smooth aggregation, {cite}`vanek1996algebraic`, {cite}`vanek2001convergence`, variant of PETSc's algebraic multigrid (`PCGAMG`) preconditioner
2167: Options Database Keys:
2168: + -pc_gamg_agg_nsmooths nsmooth - number of smoothing steps to use with smooth aggregation to construct prolongation
2169: . -pc_gamg_prolongator_filter thr - relative threshold for block filtering of the prolongator, preserving the near-null space (0=disabled, 0.01-0.1=typical)
2170: . -pc_gamg_prolongator_filter_scale scale - per-level scaling of the prolongator filter threshold (1.0=default)
2171: . -pc_gamg_aggressive_coarsening n - number of aggressive coarsening (MIS-2 or square graph) levels from finest.
2172: . -pc_gamg_aggressive_square_graph (true|false) - use square graph ($A^T A$), alternative is MIS-k (k=2), for aggressive coarsening
2173: . -pc_gamg_mis_k_minimum_degree_ordering (true|false) - use minimum degree ordering in greedy MIS algorithm
2174: . -pc_gamg_asm_hem_aggs n - number of HEM aggregation steps for ASM smoother
2175: - -pc_gamg_aggressive_mis_k n - number (k) distance in MIS coarsening (>2 is 'aggressive')
2177: Level: intermediate
2179: Notes:
2180: To obtain good performance for `PCGAMG` for vector valued problems you must
2181: call `MatSetBlockSize()` to indicate the number of degrees of freedom per grid point.
2182: Call `MatSetNearNullSpace()` (or `PCSetCoordinates()` if solving the equations of elasticity) to indicate the near null space of the operator
2184: When `-pc_gamg_aggressive_square_graph` is used, the coarsening is obtained by first squaring the graph and then applying, by default, a
2185: MIS-1 coarsening with `MatCoarsenApply()` on the squared graph.
2187: The many options for `PCMG` and `PCGAMG` such as controlling the smoothers on each level etc. also work for `PCGAMGAGG`
2189: .seealso: `PCGAMG`, [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCCreate()`, `PCSetType()`,
2190: `MatSetBlockSize()`, `PCMGType`, `PCSetCoordinates()`, `MatSetNearNullSpace()`, `PCGAMGSetType()`,
2191: `PCGAMGAGG`, `PCGAMGGEO`, `PCGAMGCLASSICAL`, `PCGAMGSetProcEqLim()`, `PCGAMGSetCoarseEqLim()`, `PCGAMGSetRepartition()`, `PCGAMGRegister()`,
2192: `PCGAMGSetReuseInterpolation()`, `PCGAMGASMSetUseAggs()`, `PCGAMGSetParallelCoarseGridSolve()`, `PCGAMGSetNlevels()`, `PCGAMGSetThreshold()`,
2193: `PCGAMGGetType()`, `PCGAMGSetUseSAEstEig()`
2194: M*/
2195: PetscErrorCode PCCreateGAMG_AGG(PC pc)
2196: {
2197: PC_MG *mg = (PC_MG *)pc->data;
2198: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
2199: PC_GAMG_AGG *pc_gamg_agg;
2201: PetscFunctionBegin;
2202: /* create sub context for SA */
2203: PetscCall(PetscNew(&pc_gamg_agg));
2204: pc_gamg->subctx = pc_gamg_agg;
2206: pc_gamg->ops->setfromoptions = PCSetFromOptions_GAMG_AGG;
2207: pc_gamg->ops->destroy = PCDestroy_GAMG_AGG;
2208: /* reset does not do anything; setup not virtual */
2210: /* set internal function pointers */
2211: pc_gamg->ops->creategraph = PCGAMGCreateGraph_AGG;
2212: pc_gamg->ops->coarsen = PCGAMGCoarsen_AGG;
2213: pc_gamg->ops->prolongator = PCGAMGConstructProlongator_AGG;
2214: pc_gamg->ops->optprolongator = PCGAMGOptimizeProlongator_AGG;
2215: pc_gamg->ops->createdefaultdata = PCSetData_AGG;
2216: pc_gamg->ops->view = PCView_GAMG_AGG;
2218: pc_gamg_agg->nsmooths = 1;
2219: pc_gamg_agg->aggressive_coarsening_levels = 1;
2220: pc_gamg_agg->use_aggressive_square_graph = PETSC_TRUE;
2221: pc_gamg_agg->use_minimum_degree_ordering = PETSC_FALSE;
2222: pc_gamg_agg->use_low_mem_filter = PETSC_FALSE;
2223: pc_gamg_agg->aggressive_mis_k = 2;
2224: pc_gamg_agg->graph_symmetrize = PETSC_TRUE;
2226: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", PCGAMGSetNSmooths_AGG));
2227: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", PCGAMGSetAggressiveLevels_AGG));
2228: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveSquareGraph_C", PCGAMGSetAggressiveSquareGraph_AGG));
2229: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetMinDegreeOrdering_C", PCGAMGMISkSetMinDegreeOrdering_AGG));
2230: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetLowMemoryFilter_C", PCGAMGSetLowMemoryFilter_AGG));
2231: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetAggressive_C", PCGAMGMISkSetAggressive_AGG));
2232: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetGraphSymmetrize_C", PCGAMGSetGraphSymmetrize_AGG));
2233: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetProlongatorFilter_C", PCGAMGSetProlongatorFilter_AGG));
2234: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGGetProlongatorFilter_C", PCGAMGGetProlongatorFilter_AGG));
2235: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetProlongatorFilterScale_C", PCGAMGSetProlongatorFilterScale_AGG));
2236: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGGetProlongatorFilterScale_C", PCGAMGGetProlongatorFilterScale_AGG));
2237: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_AGG));
2238: PetscFunctionReturn(PETSC_SUCCESS);
2239: }