Actual source code: agg.c
1: /*
2: GAMG geometric-algebric multigrid PC - Mark Adams 2011
3: */
5: #include <../src/ksp/pc/impls/gamg/gamg.h>
6: #include <petscblaslapack.h>
7: #include <petscdm.h>
8: #include <petsc/private/kspimpl.h>
10: typedef struct {
11: PetscInt nsmooths; // number of smoothing steps to construct prolongation
12: PetscInt aggressive_coarsening_levels; // number of aggressive coarsening levels (square or MISk)
13: PetscInt aggressive_mis_k; // the k in MIS-k
14: PetscBool use_aggressive_square_graph;
15: PetscBool use_minimum_degree_ordering;
16: PetscBool use_low_mem_filter;
17: PetscBool graph_symmetrize;
18: MatCoarsen crs;
19: } PC_GAMG_AGG;
21: /*@
22: PCGAMGSetNSmooths - Set number of smoothing steps (1 is typical) used to construct the prolongation operator
24: Logically Collective
26: Input Parameters:
27: + pc - the preconditioner context
28: - n - the number of smooths, default is 1
30: Options Database Key:
31: . -pc_gamg_agg_nsmooths nsmooth - number of smoothing steps to use
33: Level: intermediate
35: Note:
36: This is a different concept from the number smoothing steps used during the linear solution process which
37: can be set with `-mg_levels_ksp_max_it`
39: Developer Note:
40: This should be named `PCGAMGAGGSetNSmooths()`.
42: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCMG`, `PCGAMG`
43: @*/
44: PetscErrorCode PCGAMGSetNSmooths(PC pc, PetscInt n)
45: {
46: PetscFunctionBegin;
49: PetscTryMethod(pc, "PCGAMGSetNSmooths_C", (PC, PetscInt), (pc, n));
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: static PetscErrorCode PCGAMGSetNSmooths_AGG(PC pc, PetscInt n)
54: {
55: PC_MG *mg = (PC_MG *)pc->data;
56: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
57: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
59: PetscFunctionBegin;
60: pc_gamg_agg->nsmooths = n;
61: PetscFunctionReturn(PETSC_SUCCESS);
62: }
64: /*@
65: PCGAMGSetAggressiveLevels - Use aggressive coarsening on first n levels
67: Logically Collective
69: Input Parameters:
70: + pc - the preconditioner context
71: - n - 0, 1 or more, the default is 1
73: Options Database Key:
74: . -pc_gamg_aggressive_coarsening n - the number of coarsenings to do aggressively
76: Level: intermediate
78: Note:
79: By default, aggressive coarsening squares the matrix (computes $A^T A$) before coarsening.
80: Calling `PCGAMGSetAggressiveSquareGraph()` with a value of `PETSC_FALSE` changes the aggressive coarsening strategy to use MIS-k, see `PCGAMGMISkSetAggressive()`.
82: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGMISkSetAggressive()`,
83: `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`, `PCGAMGSetLowMemoryFilter()`
84: @*/
85: PetscErrorCode PCGAMGSetAggressiveLevels(PC pc, PetscInt n)
86: {
87: PetscFunctionBegin;
90: PetscTryMethod(pc, "PCGAMGSetAggressiveLevels_C", (PC, PetscInt), (pc, n));
91: PetscFunctionReturn(PETSC_SUCCESS);
92: }
94: /*@
95: PCGAMGMISkSetAggressive - Number (k) distance in MIS coarsening (> 2 is aggressive)
97: Logically Collective
99: Input Parameters:
100: + pc - the preconditioner context
101: - n - 1 or more (default = 2)
103: Options Database Key:
104: . -pc_gamg_aggressive_mis_k n - the distance to use
106: Level: intermediate
108: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`,
109: `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`, `PCGAMGSetLowMemoryFilter()`
110: @*/
111: PetscErrorCode PCGAMGMISkSetAggressive(PC pc, PetscInt n)
112: {
113: PetscFunctionBegin;
116: PetscTryMethod(pc, "PCGAMGMISkSetAggressive_C", (PC, PetscInt), (pc, n));
117: PetscFunctionReturn(PETSC_SUCCESS);
118: }
120: /*@
121: PCGAMGSetAggressiveSquareGraph - Use graph square ($A^T A$) for aggressive coarsening. Coarsening is slower than the alternative (MIS-2), which is faster and uses less memory
123: Logically Collective
125: Input Parameters:
126: + pc - the preconditioner context
127: - b - default true
129: Options Database Key:
130: . -pc_gamg_aggressive_square_graph (true|false) - whether to use the graph square to aggressively coarsen
132: Level: intermediate
134: Notes:
135: If `b` is `PETSC_FALSE` then MIS-k is used for aggressive coarsening, see `PCGAMGMISkSetAggressive()`
137: Squaring the matrix to perform the aggressive coarsening is slower and requires more memory than using MIS-k, but may result in a better preconditioner
138: that converges faster.
140: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`, `PCGAMGMISkSetAggressive()`, `PCGAMGMISkSetMinDegreeOrdering()`, `PCGAMGSetLowMemoryFilter()`
141: @*/
142: PetscErrorCode PCGAMGSetAggressiveSquareGraph(PC pc, PetscBool b)
143: {
144: PetscFunctionBegin;
147: PetscTryMethod(pc, "PCGAMGSetAggressiveSquareGraph_C", (PC, PetscBool), (pc, b));
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: /*@
152: PCGAMGMISkSetMinDegreeOrdering - Use minimum degree ordering in greedy MIS algorithm
154: Logically Collective
156: Input Parameters:
157: + pc - the preconditioner context
158: - b - default false
160: Options Database Key:
161: . -pc_gamg_mis_k_minimum_degree_ordering (true|false) - use the minimum degree ordering
163: Level: intermediate
165: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`,
166: `PCGAMGSetAggressiveLevels()`, `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGSetLowMemoryFilter()`
167: @*/
168: PetscErrorCode PCGAMGMISkSetMinDegreeOrdering(PC pc, PetscBool b)
169: {
170: PetscFunctionBegin;
173: PetscTryMethod(pc, "PCGAMGMISkSetMinDegreeOrdering_C", (PC, PetscBool), (pc, b));
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: /*@
178: PCGAMGSetLowMemoryFilter - Use low memory graph/matrix filter
180: Logically Collective
182: Input Parameters:
183: + pc - the preconditioner context
184: - b - default false
186: Options Database Key:
187: . -pc_gamg_low_memory_threshold_filter (true|false) - use the low memory filter
189: Level: intermediate
191: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`,
192: `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`
193: @*/
194: PetscErrorCode PCGAMGSetLowMemoryFilter(PC pc, PetscBool b)
195: {
196: PetscFunctionBegin;
199: PetscTryMethod(pc, "PCGAMGSetLowMemoryFilter_C", (PC, PetscBool), (pc, b));
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: /*@
204: PCGAMGSetGraphSymmetrize - Symmetrize graph used for coarsening. Defaults to true, but if matrix has symmetric attribute, then not needed since the graph is already known to be symmetric
206: Logically Collective
208: Input Parameters:
209: + pc - the preconditioner context
210: - b - default true
212: Options Database Key:
213: . -pc_gamg_graph_symmetrize (true|false) - symmetrize the graph
215: Level: intermediate
217: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`, `MatCreateGraph()`,
218: `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`
219: @*/
220: PetscErrorCode PCGAMGSetGraphSymmetrize(PC pc, PetscBool b)
221: {
222: PetscFunctionBegin;
225: PetscTryMethod(pc, "PCGAMGSetGraphSymmetrize_C", (PC, PetscBool), (pc, b));
226: PetscFunctionReturn(PETSC_SUCCESS);
227: }
229: /*@
230: PCGAMGSetProlongatorFilter - Set the relative threshold for block filtering of the prolongator (a kernel-preserving correction is applied afterward)
232: Logically Collective
234: Input Parameters:
235: + pc - the preconditioner context
236: - 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
237: largest such block norm in that fine node's block row (0 disables filtering)
239: Options Database Key:
240: . -pc_gamg_prolongator_filter thr - relative threshold for block filtering of the prolongator (0=disabled, 0.01-0.1=typical)
242: Level: intermediate
244: Notes:
245: 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
246: 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
247: 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.
248: The comparison is strict, so the strongest block of a fine node always survives. On coarser levels the threshold is scaled by `PCGAMGSetProlongatorFilterScale()`.
249: Dropping whole blocks (rather than individual entries) keeps complete coarse-node blocks in every surviving fine row, so the near-null space correction below
250: remains full rank. The dropped entries are removed from the sparsity pattern with `MatEliminateZeros()`; on matrix types that do not implement it, and on
251: HIPSPARSE where it is bypassed due to a known issue, they are zeroed but remain in the pattern, so the coarse operators are unchanged in structure and the
252: 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, kk, n, j;
1017: Mat_MPIAIJ *mpimat_2 = NULL, *mpimat_1 = NULL;
1018: const PetscInt nloc = Gmat_2->rmap->n;
1019: PetscScalar *cpcol_1_state, *cpcol_2_state, *cpcol_2_par_orig, *lid_parent_gid;
1020: PetscInt *lid_cprowID_1 = NULL;
1021: NState *lid_state;
1022: Vec ghost_par_orig2;
1023: PetscMPIInt rank;
1025: PetscFunctionBegin;
1026: PetscCall(PetscObjectGetComm((PetscObject)Gmat_2, &comm));
1027: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1028: PetscCall(MatGetOwnershipRange(Gmat_1, &my0, &Iend));
1030: /* get submatrices */
1031: PetscCall(PetscStrbeginswith(((PetscObject)Gmat_1)->type_name, MATMPIAIJ, &isMPI));
1032: PetscCall(PetscInfo(pc, "isMPI = %s\n", isMPI ? "yes" : "no"));
1033: PetscCall(PetscMalloc3(nloc, &lid_state, nloc, &lid_parent_gid, nloc, &lid_cprowID_1));
1034: for (lid = 0; lid < nloc; lid++) lid_cprowID_1[lid] = -1;
1035: if (isMPI) {
1036: /* grab matrix objects */
1037: mpimat_2 = (Mat_MPIAIJ *)Gmat_2->data;
1038: mpimat_1 = (Mat_MPIAIJ *)Gmat_1->data;
1039: matA_1 = (Mat_SeqAIJ *)mpimat_1->A->data;
1040: matB_1 = (Mat_SeqAIJ *)mpimat_1->B->data;
1042: /* force compressed row storage for B matrix in AuxMat */
1043: PetscCall(MatCheckCompressedRow(mpimat_1->B, matB_1->nonzerorowcnt, &matB_1->compressedrow, matB_1->i, Gmat_1->rmap->n, -1.0));
1044: for (ix = 0; ix < matB_1->compressedrow.nrows; ix++) {
1045: PetscInt lid = matB_1->compressedrow.rindex[ix];
1047: PetscCheck(lid <= nloc && lid >= -1, PETSC_COMM_SELF, PETSC_ERR_USER, "lid %" PetscInt_FMT " out of range. nloc = %" PetscInt_FMT, lid, nloc);
1048: if (lid != -1) lid_cprowID_1[lid] = ix;
1049: }
1050: } else {
1051: PetscBool isAIJ;
1053: PetscCall(PetscStrbeginswith(((PetscObject)Gmat_1)->type_name, MATSEQAIJ, &isAIJ));
1054: PetscCheck(isAIJ, PETSC_COMM_SELF, PETSC_ERR_USER, "Require AIJ matrix.");
1055: matA_1 = (Mat_SeqAIJ *)Gmat_1->data;
1056: }
1057: if (nloc > 0) PetscCheck(!matB_1 || matB_1->compressedrow.use, PETSC_COMM_SELF, PETSC_ERR_PLIB, "matB_1 && !matB_1->compressedrow.use: PETSc bug???");
1058: /* get state of locals and selected gid for deleted */
1059: for (lid = 0; lid < nloc; lid++) {
1060: lid_parent_gid[lid] = -1.0;
1061: lid_state[lid] = DELETED;
1062: }
1064: /* set lid_state */
1065: for (lid = 0; lid < nloc; lid++) {
1066: PetscCDIntNd *pos;
1068: PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
1069: if (pos) {
1070: PetscInt gid1;
1072: PetscCall(PetscCDIntNdGetID(pos, &gid1));
1073: PetscCheck(gid1 == lid + my0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "gid1 %" PetscInt_FMT " != lid %" PetscInt_FMT " + my0 %" PetscInt_FMT, gid1, lid, my0);
1074: lid_state[lid] = gid1;
1075: }
1076: }
1078: /* map local to selected local, DELETED means a ghost owns it */
1079: for (lid = 0; lid < nloc; lid++) {
1080: NState state = lid_state[lid];
1082: if (IS_SELECTED(state)) {
1083: PetscCDIntNd *pos;
1085: PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
1086: while (pos) {
1087: PetscInt gid1;
1089: PetscCall(PetscCDIntNdGetID(pos, &gid1));
1090: PetscCall(PetscCDGetNextPos(aggs_2, lid, &pos));
1091: if (gid1 >= my0 && gid1 < Iend) lid_parent_gid[gid1 - my0] = (PetscScalar)(lid + my0);
1092: }
1093: }
1094: }
1095: /* get 'cpcol_1/2_state' & cpcol_2_par_orig - uses mpimat_1/2->lvec for temp space */
1096: if (isMPI) {
1097: Vec tempVec;
1099: /* get 'cpcol_1_state' */
1100: PetscCall(MatCreateVecs(Gmat_1, &tempVec, NULL));
1101: for (kk = 0, j = my0; kk < nloc; kk++, j++) {
1102: PetscScalar v = (PetscScalar)lid_state[kk];
1104: PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
1105: }
1106: PetscCall(VecAssemblyBegin(tempVec));
1107: PetscCall(VecAssemblyEnd(tempVec));
1108: PetscCall(VecScatterBegin(mpimat_1->Mvctx, tempVec, mpimat_1->lvec, INSERT_VALUES, SCATTER_FORWARD));
1109: PetscCall(VecScatterEnd(mpimat_1->Mvctx, tempVec, mpimat_1->lvec, INSERT_VALUES, SCATTER_FORWARD));
1110: PetscCall(VecGetArray(mpimat_1->lvec, &cpcol_1_state));
1111: /* get 'cpcol_2_state' */
1112: PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, mpimat_2->lvec, INSERT_VALUES, SCATTER_FORWARD));
1113: PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, mpimat_2->lvec, INSERT_VALUES, SCATTER_FORWARD));
1114: PetscCall(VecGetArray(mpimat_2->lvec, &cpcol_2_state));
1115: /* get 'cpcol_2_par_orig' */
1116: for (kk = 0, j = my0; kk < nloc; kk++, j++) {
1117: PetscScalar v = lid_parent_gid[kk];
1119: PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
1120: }
1121: PetscCall(VecAssemblyBegin(tempVec));
1122: PetscCall(VecAssemblyEnd(tempVec));
1123: PetscCall(VecDuplicate(mpimat_2->lvec, &ghost_par_orig2));
1124: PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghost_par_orig2, INSERT_VALUES, SCATTER_FORWARD));
1125: PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghost_par_orig2, INSERT_VALUES, SCATTER_FORWARD));
1126: PetscCall(VecGetArray(ghost_par_orig2, &cpcol_2_par_orig));
1128: PetscCall(VecDestroy(&tempVec));
1129: } /* ismpi */
1130: for (lid = 0; lid < nloc; lid++) {
1131: NState state = lid_state[lid];
1133: if (IS_SELECTED(state)) {
1134: /* steal locals */
1135: ii = matA_1->i;
1136: n = ii[lid + 1] - ii[lid];
1137: idx = matA_1->j + ii[lid];
1138: for (j = 0; j < n; j++) {
1139: PetscInt lidj = idx[j], sgid;
1140: NState statej = lid_state[lidj];
1142: if (statej == DELETED && (sgid = (PetscInt)PetscRealPart(lid_parent_gid[lidj])) != lid + my0) { /* steal local */
1143: lid_parent_gid[lidj] = (PetscScalar)(lid + my0); /* send this if sgid is not local */
1144: if (sgid >= my0 && sgid < Iend) { /* I'm stealing this local from a local sgid */
1145: PetscInt hav = 0, slid = sgid - my0, gidj = lidj + my0;
1146: PetscCDIntNd *pos, *last = NULL;
1148: /* looking for local from local so id_llist_2 works */
1149: PetscCall(PetscCDGetHeadPos(aggs_2, slid, &pos));
1150: while (pos) {
1151: PetscInt gid;
1153: PetscCall(PetscCDIntNdGetID(pos, &gid));
1154: if (gid == gidj) {
1155: PetscCheck(last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "last cannot be null");
1156: PetscCall(PetscCDRemoveNextNode(aggs_2, slid, last));
1157: PetscCall(PetscCDAppendNode(aggs_2, lid, pos));
1158: hav = 1;
1159: break;
1160: } else last = pos;
1161: PetscCall(PetscCDGetNextPos(aggs_2, slid, &pos));
1162: }
1163: if (hav != 1) {
1164: PetscCheck(hav, PETSC_COMM_SELF, PETSC_ERR_PLIB, "failed to find adj in 'selected' lists - structurally unsymmetric matrix");
1165: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "found node %" PetscInt_FMT " times???", hav);
1166: }
1167: } else { /* I'm stealing this local, owned by a ghost */
1168: 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.",
1169: ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "");
1170: PetscCall(PetscCDAppendID(aggs_2, lid, lidj + my0));
1171: }
1172: }
1173: } /* local neighbors */
1174: } else if (state == DELETED /* && lid_cprowID_1 */) {
1175: PetscInt sgidold = (PetscInt)PetscRealPart(lid_parent_gid[lid]);
1177: /* see if I have a selected ghost neighbor that will steal me */
1178: if ((ix = lid_cprowID_1[lid]) != -1) {
1179: ii = matB_1->compressedrow.i;
1180: n = ii[ix + 1] - ii[ix];
1181: idx = matB_1->j + ii[ix];
1182: for (j = 0; j < n; j++) {
1183: PetscInt cpid = idx[j];
1184: NState statej = (NState)PetscRealPart(cpcol_1_state[cpid]);
1186: if (IS_SELECTED(statej) && sgidold != statej) { /* ghost will steal this, remove from my list */
1187: lid_parent_gid[lid] = (PetscScalar)statej; /* send who selected */
1188: if (sgidold >= my0 && sgidold < Iend) { /* this was mine */
1189: PetscInt hav = 0, oldslidj = sgidold - my0;
1190: PetscCDIntNd *pos, *last = NULL;
1192: /* remove from 'oldslidj' list */
1193: PetscCall(PetscCDGetHeadPos(aggs_2, oldslidj, &pos));
1194: while (pos) {
1195: PetscInt gid;
1197: PetscCall(PetscCDIntNdGetID(pos, &gid));
1198: if (lid + my0 == gid) {
1199: /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove lid from oldslidj list *\/ */
1200: PetscCheck(last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "last cannot be null");
1201: PetscCall(PetscCDRemoveNextNode(aggs_2, oldslidj, last));
1202: /* ghost (PetscScalar)statej will add this later */
1203: hav = 1;
1204: break;
1205: } else last = pos;
1206: PetscCall(PetscCDGetNextPos(aggs_2, oldslidj, &pos));
1207: }
1208: if (hav != 1) {
1209: PetscCheck(hav, PETSC_COMM_SELF, PETSC_ERR_PLIB, "failed to find (hav=%" PetscInt_FMT ") adj in 'selected' lists - structurally unsymmetric matrix", hav);
1210: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "found node %" PetscInt_FMT " times???", hav);
1211: }
1212: } else {
1213: /* TODO: ghosts remove this later */
1214: }
1215: }
1216: }
1217: }
1218: } /* selected/deleted */
1219: } /* node loop */
1221: if (isMPI) {
1222: PetscScalar *cpcol_2_parent, *cpcol_2_gid;
1223: Vec tempVec, ghostgids2, ghostparents2;
1224: PetscInt cpid, nghost_2;
1225: PetscHMapI gid_cpid;
1227: PetscCall(VecGetSize(mpimat_2->lvec, &nghost_2));
1228: PetscCall(MatCreateVecs(Gmat_2, &tempVec, NULL));
1230: /* get 'cpcol_2_parent' */
1231: for (kk = 0, j = my0; kk < nloc; kk++, j++) PetscCall(VecSetValues(tempVec, 1, &j, &lid_parent_gid[kk], INSERT_VALUES));
1232: PetscCall(VecAssemblyBegin(tempVec));
1233: PetscCall(VecAssemblyEnd(tempVec));
1234: PetscCall(VecDuplicate(mpimat_2->lvec, &ghostparents2));
1235: PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghostparents2, INSERT_VALUES, SCATTER_FORWARD));
1236: PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghostparents2, INSERT_VALUES, SCATTER_FORWARD));
1237: PetscCall(VecGetArray(ghostparents2, &cpcol_2_parent));
1239: /* get 'cpcol_2_gid' */
1240: for (kk = 0, j = my0; kk < nloc; kk++, j++) {
1241: PetscScalar v = (PetscScalar)j;
1243: PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
1244: }
1245: PetscCall(VecAssemblyBegin(tempVec));
1246: PetscCall(VecAssemblyEnd(tempVec));
1247: PetscCall(VecDuplicate(mpimat_2->lvec, &ghostgids2));
1248: PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghostgids2, INSERT_VALUES, SCATTER_FORWARD));
1249: PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghostgids2, INSERT_VALUES, SCATTER_FORWARD));
1250: PetscCall(VecGetArray(ghostgids2, &cpcol_2_gid));
1251: PetscCall(VecDestroy(&tempVec));
1253: /* look for deleted ghosts and add to table */
1254: PetscCall(PetscHMapICreateWithSize(2 * nghost_2 + 1, &gid_cpid));
1255: for (cpid = 0; cpid < nghost_2; cpid++) {
1256: NState state = (NState)PetscRealPart(cpcol_2_state[cpid]);
1258: if (state == DELETED) {
1259: PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
1260: PetscInt sgid_old = (PetscInt)PetscRealPart(cpcol_2_par_orig[cpid]);
1262: if (sgid_old == -1 && sgid_new != -1) {
1263: PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
1265: PetscCall(PetscHMapISet(gid_cpid, gid, cpid));
1266: }
1267: }
1268: }
1270: /* look for deleted ghosts and see if they moved - remove it */
1271: for (lid = 0; lid < nloc; lid++) {
1272: NState state = lid_state[lid];
1274: if (IS_SELECTED(state)) {
1275: PetscCDIntNd *pos, *last = NULL;
1277: /* look for deleted ghosts and see if they moved */
1278: PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
1279: while (pos) {
1280: PetscInt gid;
1282: PetscCall(PetscCDIntNdGetID(pos, &gid));
1283: if (gid < my0 || gid >= Iend) {
1284: PetscCall(PetscHMapIGet(gid_cpid, gid, &cpid));
1285: if (cpid != -1) {
1286: /* a moved ghost - */
1287: /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove 'flid' from list *\/ */
1288: PetscCall(PetscCDRemoveNextNode(aggs_2, lid, last));
1289: } else last = pos;
1290: } else last = pos;
1292: PetscCall(PetscCDGetNextPos(aggs_2, lid, &pos));
1293: } /* loop over list of deleted */
1294: } /* selected */
1295: }
1296: PetscCall(PetscHMapIDestroy(&gid_cpid));
1298: /* look at ghosts, see if they changed - and it */
1299: for (cpid = 0; cpid < nghost_2; cpid++) {
1300: PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
1302: if (sgid_new >= my0 && sgid_new < Iend) { /* this is mine */
1303: PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
1304: PetscInt slid_new = sgid_new - my0, hav = 0;
1305: PetscCDIntNd *pos;
1307: /* search for this gid to see if I have it */
1308: PetscCall(PetscCDGetHeadPos(aggs_2, slid_new, &pos));
1309: while (pos) {
1310: PetscInt gidj;
1312: PetscCall(PetscCDIntNdGetID(pos, &gidj));
1313: PetscCall(PetscCDGetNextPos(aggs_2, slid_new, &pos));
1315: if (gidj == gid) {
1316: hav = 1;
1317: break;
1318: }
1319: }
1320: if (hav != 1) {
1321: /* insert 'flidj' into head of llist */
1322: PetscCall(PetscCDAppendID(aggs_2, slid_new, gid));
1323: }
1324: }
1325: }
1326: PetscCall(VecRestoreArray(mpimat_1->lvec, &cpcol_1_state));
1327: PetscCall(VecRestoreArray(mpimat_2->lvec, &cpcol_2_state));
1328: PetscCall(VecRestoreArray(ghostparents2, &cpcol_2_parent));
1329: PetscCall(VecRestoreArray(ghostgids2, &cpcol_2_gid));
1330: PetscCall(VecDestroy(&ghostgids2));
1331: PetscCall(VecDestroy(&ghostparents2));
1332: PetscCall(VecDestroy(&ghost_par_orig2));
1333: }
1334: PetscCall(PetscFree3(lid_state, lid_parent_gid, lid_cprowID_1));
1335: PetscFunctionReturn(PETSC_SUCCESS);
1336: }
1338: /*
1339: PCGAMGCoarsen_AGG - supports squaring the graph (deprecated) and new graph for
1340: communication of QR data used with HEM and MISk coarsening
1342: Input Parameter:
1343: . a_pc - this
1345: Input/Output Parameter:
1346: . a_Gmat1 - graph to coarsen (in), graph off processor edges for QR gather scatter (out)
1348: Output Parameter:
1349: . agg_lists - list of aggregates
1351: */
1352: static PetscErrorCode PCGAMGCoarsen_AGG(PC a_pc, Mat *a_Gmat1, PetscCoarsenData **agg_lists)
1353: {
1354: PC_MG *mg = (PC_MG *)a_pc->data;
1355: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1356: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
1357: Mat Gmat2, Gmat1 = *a_Gmat1; /* aggressive graph */
1358: IS perm;
1359: PetscInt Istart, Iend, Ii, nloc, bs, nn;
1360: PetscInt *permute, *degree;
1361: PetscBool *bIndexSet;
1362: PetscReal hashfact;
1363: PetscInt iSwapIndex;
1364: PetscRandom random;
1365: MPI_Comm comm;
1367: PetscFunctionBegin;
1368: PetscCall(PetscObjectGetComm((PetscObject)Gmat1, &comm));
1369: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
1370: PetscCall(MatGetLocalSize(Gmat1, &nn, NULL));
1371: PetscCall(MatGetBlockSize(Gmat1, &bs));
1372: PetscCheck(bs == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "bs %" PetscInt_FMT " must be 1", bs);
1373: nloc = nn / bs;
1374: /* get MIS aggs - randomize */
1375: PetscCall(PetscMalloc2(nloc, &permute, nloc, °ree));
1376: PetscCall(PetscCalloc1(nloc, &bIndexSet));
1377: for (Ii = 0; Ii < nloc; Ii++) permute[Ii] = Ii;
1378: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &random));
1379: PetscCall(MatGetOwnershipRange(Gmat1, &Istart, &Iend));
1380: for (Ii = 0; Ii < nloc; Ii++) {
1381: PetscInt nc;
1383: PetscCall(MatGetRow(Gmat1, Istart + Ii, &nc, NULL, NULL));
1384: degree[Ii] = nc;
1385: PetscCall(MatRestoreRow(Gmat1, Istart + Ii, &nc, NULL, NULL));
1386: }
1387: for (Ii = 0; Ii < nloc; Ii++) {
1388: PetscCall(PetscRandomGetValueReal(random, &hashfact));
1389: iSwapIndex = (PetscInt)(hashfact * nloc) % nloc;
1390: if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) {
1391: PetscInt iTemp = permute[iSwapIndex];
1393: permute[iSwapIndex] = permute[Ii];
1394: permute[Ii] = iTemp;
1395: iTemp = degree[iSwapIndex];
1396: degree[iSwapIndex] = degree[Ii];
1397: degree[Ii] = iTemp;
1398: bIndexSet[iSwapIndex] = PETSC_TRUE;
1399: }
1400: }
1401: // apply minimum degree ordering -- NEW
1402: if (pc_gamg_agg->use_minimum_degree_ordering) PetscCall(PetscSortIntWithArray(nloc, degree, permute));
1403: PetscCall(PetscFree(bIndexSet));
1404: PetscCall(PetscRandomDestroy(&random));
1405: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_USE_POINTER, &perm));
1406: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0));
1407: // square graph
1408: 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));
1409: else Gmat2 = Gmat1;
1410: // switch to old MIS-1 for square graph
1411: if (pc_gamg->current_level < pc_gamg_agg->aggressive_coarsening_levels) {
1412: if (!pc_gamg_agg->use_aggressive_square_graph) PetscCall(MatCoarsenMISKSetDistance(pc_gamg_agg->crs, pc_gamg_agg->aggressive_mis_k)); // hardwire to MIS-2
1413: else PetscCall(MatCoarsenSetType(pc_gamg_agg->crs, MATCOARSENMIS)); // old MIS -- side effect
1414: } else if (pc_gamg_agg->use_aggressive_square_graph && pc_gamg_agg->aggressive_coarsening_levels > 0) { // we reset the MIS
1415: const char *prefix;
1417: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)a_pc, &prefix));
1418: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc_gamg_agg->crs, prefix));
1419: PetscCall(MatCoarsenSetFromOptions(pc_gamg_agg->crs)); // get the default back on non-aggressive levels when square graph switched to old MIS
1420: }
1421: PetscCall(MatCoarsenSetAdjacency(pc_gamg_agg->crs, Gmat2));
1422: PetscCall(MatCoarsenSetStrictAggs(pc_gamg_agg->crs, PETSC_TRUE));
1423: PetscCall(MatCoarsenSetGreedyOrdering(pc_gamg_agg->crs, perm));
1424: PetscCall(MatCoarsenApply(pc_gamg_agg->crs));
1425: PetscCall(MatCoarsenGetData(pc_gamg_agg->crs, agg_lists)); /* output */
1427: PetscCall(ISDestroy(&perm));
1428: PetscCall(PetscFree2(permute, degree));
1429: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0));
1431: if (Gmat2 != Gmat1) { // square graph, we need ghosts for selected
1432: PetscCoarsenData *llist = *agg_lists;
1434: PetscCall(fixAggregatesWithSquare(a_pc, Gmat2, Gmat1, *agg_lists));
1435: PetscCall(MatDestroy(&Gmat1));
1436: *a_Gmat1 = Gmat2; /* output */
1437: PetscCall(PetscCDSetMat(llist, *a_Gmat1)); /* Need a graph with ghosts here */
1438: }
1439: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
1440: PetscFunctionReturn(PETSC_SUCCESS);
1441: }
1443: /*
1444: PCGAMGConstructProlongator_AGG
1446: Input Parameter:
1447: . pc - this
1448: . Amat - matrix on this fine level
1449: . Graph - used to get ghost data for nodes in
1450: . agg_lists - list of aggregates
1451: Output Parameter:
1452: . a_P_out - prolongation operator to the next level
1453: */
1454: static PetscErrorCode PCGAMGConstructProlongator_AGG(PC pc, Mat Amat, PetscCoarsenData *agg_lists, Mat *a_P_out)
1455: {
1456: PC_MG *mg = (PC_MG *)pc->data;
1457: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1458: const PetscInt col_bs = pc_gamg->data_cell_cols;
1459: PetscInt Istart, Iend, nloc, ii, jj, kk, my0, nLocalSelected, bs;
1460: Mat Gmat, Prol;
1461: PetscMPIInt size;
1462: MPI_Comm comm;
1463: PetscReal *data_w_ghost;
1464: PetscInt myCrs0, nbnodes = 0, *flid_fgid;
1465: MatType mtype;
1467: PetscFunctionBegin;
1468: PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
1469: PetscCheck(col_bs >= 1, comm, PETSC_ERR_PLIB, "Column bs cannot be less than 1");
1470: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1471: PetscCallMPI(MPI_Comm_size(comm, &size));
1472: PetscCall(MatGetOwnershipRange(Amat, &Istart, &Iend));
1473: PetscCall(MatGetBlockSize(Amat, &bs));
1474: nloc = (Iend - Istart) / bs;
1475: my0 = Istart / bs;
1476: 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);
1477: PetscCall(PetscCDGetMat(agg_lists, &Gmat)); // get auxiliary matrix for ghost edges for size > 1
1479: /* get 'nLocalSelected' */
1480: for (ii = 0, nLocalSelected = 0; ii < nloc; ii++) {
1481: PetscBool ise;
1483: /* filter out singletons 0 or 1? */
1484: PetscCall(PetscCDIsEmptyAt(agg_lists, ii, &ise));
1485: if (!ise) nLocalSelected++;
1486: }
1488: /* create prolongator, create P matrix */
1489: PetscCall(MatGetType(Amat, &mtype));
1490: PetscCall(MatCreate(comm, &Prol));
1491: PetscCall(MatSetSizes(Prol, nloc * bs, nLocalSelected * col_bs, PETSC_DETERMINE, PETSC_DETERMINE));
1492: PetscCall(MatSetBlockSizes(Prol, bs, col_bs)); // should this be before MatSetSizes?
1493: PetscCall(MatSetType(Prol, mtype));
1494: #if PetscDefined(HAVE_DEVICE)
1495: PetscBool flg;
1496: PetscCall(MatBoundToCPU(Amat, &flg));
1497: PetscCall(MatBindToCPU(Prol, flg));
1498: if (flg) PetscCall(MatSetBindingPropagates(Prol, PETSC_TRUE));
1499: #endif
1500: PetscCall(MatSeqAIJSetPreallocation(Prol, col_bs, NULL));
1501: PetscCall(MatMPIAIJSetPreallocation(Prol, col_bs, NULL, col_bs, NULL));
1503: /* can get all points "removed" */
1504: PetscCall(MatGetSize(Prol, &kk, &ii));
1505: if (!ii) {
1506: PetscCall(PetscInfo(pc, "%s: No selected points on coarse grid\n", ((PetscObject)pc)->prefix));
1507: PetscCall(MatDestroy(&Prol));
1508: *a_P_out = NULL; /* out */
1509: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1510: PetscFunctionReturn(PETSC_SUCCESS);
1511: }
1512: PetscCall(PetscInfo(pc, "%s: New grid %" PetscInt_FMT " nodes\n", ((PetscObject)pc)->prefix, ii / col_bs));
1513: PetscCall(MatGetOwnershipRangeColumn(Prol, &myCrs0, &kk));
1515: 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);
1516: myCrs0 = myCrs0 / col_bs;
1517: 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);
1519: /* create global vector of data in 'data_w_ghost' */
1520: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0));
1521: if (size > 1) { /* get ghost null space data */
1522: PetscReal *tmp_gdata, *tmp_ldata, *tp2;
1524: PetscCall(PetscMalloc1(nloc, &tmp_ldata));
1525: for (jj = 0; jj < col_bs; jj++) {
1526: for (kk = 0; kk < bs; kk++) {
1527: PetscInt stride;
1528: const PetscReal *tp = PetscSafePointerPlusOffset(pc_gamg->data, jj * bs * nloc + kk);
1530: for (PetscInt ii = 0; ii < nloc; ii++, tp += bs) tmp_ldata[ii] = *tp;
1532: PetscCall(PCGAMGGetDataWithGhosts(Gmat, 1, tmp_ldata, &stride, &tmp_gdata));
1534: if (!jj && !kk) { /* now I know how many total nodes - allocate TODO: move below and do in one 'col_bs' call */
1535: PetscCall(PetscMalloc1(stride * bs * col_bs, &data_w_ghost));
1536: nbnodes = bs * stride;
1537: }
1538: tp2 = PetscSafePointerPlusOffset(data_w_ghost, jj * bs * stride + kk);
1539: for (PetscInt ii = 0; ii < stride; ii++, tp2 += bs) *tp2 = tmp_gdata[ii];
1540: PetscCall(PetscFree(tmp_gdata));
1541: }
1542: }
1543: PetscCall(PetscFree(tmp_ldata));
1544: } else {
1545: nbnodes = bs * nloc;
1546: data_w_ghost = pc_gamg->data;
1547: }
1549: /* get 'flid_fgid' TODO - move up to get 'stride' and do get null space data above in one step (jj loop) */
1550: if (size > 1) {
1551: PetscReal *fid_glid_loc, *fiddata;
1552: PetscInt stride;
1554: PetscCall(PetscMalloc1(nloc, &fid_glid_loc));
1555: for (kk = 0; kk < nloc; kk++) fid_glid_loc[kk] = (PetscReal)(my0 + kk);
1556: PetscCall(PCGAMGGetDataWithGhosts(Gmat, 1, fid_glid_loc, &stride, &fiddata));
1557: PetscCall(PetscMalloc1(stride, &flid_fgid)); /* copy real data to in */
1558: for (kk = 0; kk < stride; kk++) flid_fgid[kk] = (PetscInt)fiddata[kk];
1559: PetscCall(PetscFree(fiddata));
1561: PetscCheck(stride == nbnodes / bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "stride %" PetscInt_FMT " != nbnodes %" PetscInt_FMT "/bs %" PetscInt_FMT, stride, nbnodes, bs);
1562: PetscCall(PetscFree(fid_glid_loc));
1563: } else {
1564: PetscCall(PetscMalloc1(nloc, &flid_fgid));
1565: for (kk = 0; kk < nloc; kk++) flid_fgid[kk] = my0 + kk;
1566: }
1567: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0));
1568: /* get P0 */
1569: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0));
1570: {
1571: PetscReal *data_out = NULL;
1573: PetscCall(formProl0(agg_lists, bs, col_bs, myCrs0, nbnodes, data_w_ghost, flid_fgid, &data_out, Prol));
1574: PetscCall(PetscFree(pc_gamg->data));
1576: pc_gamg->data = data_out;
1577: pc_gamg->data_cell_rows = col_bs;
1578: pc_gamg->data_sz = col_bs * col_bs * nLocalSelected;
1579: }
1580: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0));
1581: if (size > 1) PetscCall(PetscFree(data_w_ghost));
1582: PetscCall(PetscFree(flid_fgid));
1584: *a_P_out = Prol; /* out */
1585: PetscCall(MatViewFromOptions(Prol, NULL, "-pc_gamg_agg_view_initial_prolongation"));
1587: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1588: PetscFunctionReturn(PETSC_SUCCESS);
1589: }
1591: // Drop small node-coupling blocks of the prolongator; see the PCGAMGSetProlongatorFilter() manual page
1592: static PetscErrorCode PCGAMGProlongatorBlockFilter_AGG(PC pc, Mat Prol, PetscInt col_bs, PetscReal thr)
1593: {
1594: PetscInt rbs, cbs, rStart, rEnd, nfn, local_nnz = 0, max_fn_cols = 0, max_row_cols = 0, ndrows = 0, zoff = 0;
1595: PetscInt *cn_gid, *fn_col, *fn_slot, *fn_rcnt, *zcols, *drow, *doff, *dcnt;
1596: PetscReal *cn_n2;
1597: PetscReal thr2 = thr * thr;
1598: PetscScalar *zeros;
1599: PetscBool no_off_proc, ishipsparse;
1601: PetscFunctionBegin;
1602: PetscCall(MatGetBlockSizes(Prol, &rbs, &cbs));
1603: /* Prol is built internally with MatSetBlockSizes(..., col_bs); a mismatch here means smoothing corrupted the block structure */
1604: 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);
1605: PetscCall(MatGetOwnershipRange(Prol, &rStart, &rEnd));
1606: 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);
1607: nfn = (rEnd - rStart) / rbs;
1609: /* Pre-pass over the row widths to size the scratch tightly: the total local nonzeros bound
1610: the dropped columns (zcols), the widest fine-node block row bounds the distinct coarse
1611: nodes of one fine node (cn_gid/cn_n2), and the widest single row bounds the entries
1612: zeroed by one MatSetValues() call (zeros) */
1613: for (PetscInt fn = 0; fn < nfn; fn++) {
1614: PetscInt fn_cols = 0;
1616: for (PetscInt rr = 0; rr < rbs; rr++) {
1617: PetscInt grow = rStart + fn * rbs + rr, ncols;
1619: PetscCall(MatGetRow(Prol, grow, &ncols, NULL, NULL));
1620: fn_cols += ncols;
1621: if (ncols > max_row_cols) max_row_cols = ncols;
1622: PetscCall(MatRestoreRow(Prol, grow, &ncols, NULL, NULL));
1623: }
1624: local_nnz += fn_cols;
1625: if (fn_cols > max_fn_cols) max_fn_cols = fn_cols;
1626: }
1628: 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));
1629: PetscCall(PetscCalloc1(max_row_cols, &zeros));
1630: PetscCall(PetscMalloc1(local_nnz, &zcols));
1631: PetscCall(PetscMalloc3(rEnd - rStart, &drow, rEnd - rStart, &doff, rEnd - rStart, &dcnt));
1633: for (PetscInt fn = 0; fn < nfn; fn++) {
1634: PetscInt ncn = 0; /* distinct coarse nodes touched by this fine node-block */
1635: PetscInt nent = 0; /* entries of this fine node-block cached in fn_col/fn_slot */
1636: PetscReal maxn2 = 0.0;
1638: /* Pass A: accumulate the per-coarse-node block Frobenius norm^2 over the rbs rows, caching each
1639: entry's column and resolved coarse-node slot so that Pass B needs no second MatGetRow() sweep */
1640: for (PetscInt rr = 0; rr < rbs; rr++) {
1641: PetscInt grow = rStart + fn * rbs + rr, ncols;
1642: const PetscInt *cols;
1643: const PetscScalar *vals;
1645: PetscCall(MatGetRow(Prol, grow, &ncols, &cols, &vals));
1646: fn_rcnt[rr] = ncols;
1647: for (PetscInt k = 0; k < ncols; k++) {
1648: PetscInt cn = cols[k] / col_bs, s = -1;
1649: PetscReal av = PetscAbsScalar(vals[k]);
1651: /* Linear dedup scan: ncn is bounded by the coarse-node degree of this
1652: fine node (typically ~5-20 for AMG-coarsened elasticity), so O(ncn)
1653: per entry is faster than a hash map at this scale */
1654: for (PetscInt t = 0; t < ncn; t++)
1655: if (cn_gid[t] == cn) {
1656: s = t;
1657: break;
1658: }
1659: if (s < 0) {
1660: s = ncn++;
1661: cn_gid[s] = cn;
1662: cn_n2[s] = 0.0;
1663: }
1664: cn_n2[s] += av * av;
1665: fn_col[nent] = cols[k];
1666: fn_slot[nent] = s;
1667: nent++;
1668: }
1669: PetscCall(MatRestoreRow(Prol, grow, &ncols, &cols, &vals));
1670: }
1671: for (PetscInt t = 0; t < ncn; t++)
1672: if (cn_n2[t] > maxn2) maxn2 = cn_n2[t];
1674: /* Pass B: collect the entries of blocks below thr*max for zeroing. The test is strict, so the
1675: strongest block of the node survives for any thr <= 1 */
1676: nent = 0;
1677: for (PetscInt rr = 0; rr < rbs; rr++) {
1678: PetscInt rec_off = zoff;
1680: for (PetscInt k = 0; k < fn_rcnt[rr]; k++, nent++)
1681: if (cn_n2[fn_slot[nent]] < thr2 * maxn2) zcols[zoff++] = fn_col[nent];
1682: if (zoff > rec_off) {
1683: drow[ndrows] = rStart + fn * rbs + rr;
1684: doff[ndrows] = rec_off;
1685: dcnt[ndrows] = zoff - rec_off;
1686: ndrows++;
1687: }
1688: }
1689: }
1691: /* all insertions are in local rows; skip the off-process assembly communication */
1692: PetscCall(MatGetOption(Prol, MAT_NO_OFF_PROC_ENTRIES, &no_off_proc));
1693: PetscCall(MatSetOption(Prol, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
1694: for (PetscInt i = 0; i < ndrows; i++) PetscCall(MatSetValues(Prol, 1, &drow[i], dcnt[i], &zcols[doff[i]], zeros, INSERT_VALUES));
1695: PetscCall(MatAssemblyBegin(Prol, MAT_FINAL_ASSEMBLY));
1696: PetscCall(MatAssemblyEnd(Prol, MAT_FINAL_ASSEMBLY));
1697: PetscCall(MatSetOption(Prol, MAT_NO_OFF_PROC_ENTRIES, no_off_proc));
1699: PetscCall(PetscFree5(cn_gid, cn_n2, fn_col, fn_slot, fn_rcnt));
1700: PetscCall(PetscFree(zeros));
1701: PetscCall(PetscFree(zcols));
1702: PetscCall(PetscFree3(drow, doff, dcnt));
1704: /* Compress out the explicit zeros just set, so the filter actually sparsifies. keep must be
1705: PETSC_FALSE: with keep, a zero whose local column index equals its local row index survives (an
1706: index-based diagonal test that is meaningless for the rectangular Prol). The compression is done
1707: in place, deliberately: a MatDuplicate()/MatHeaderReplace() copy as in MatFilter() would return
1708: the freed CSR tail to the allocator, at the cost of a peak-memory spike. MatEliminateZeros() has
1709: a known issue with HIPSPARSE (see the bypass in MatFilter()) and is not implemented by all matrix
1710: types; in those cases the zeros are left in the sparsity pattern; the step-3 correction in
1711: PCGAMGKernelPreservingFilter_AGG() skips exactly-zero entries, so a dropped block stays dropped
1712: either way, it just still costs storage here. */
1713: PetscCall(PetscObjectTypeCompareAny((PetscObject)Prol, &ishipsparse, MATSEQAIJHIPSPARSE, MATMPIAIJHIPSPARSE, ""));
1714: if (!ishipsparse && Prol->ops->eliminatezeros) PetscCall(MatEliminateZeros(Prol, PETSC_FALSE));
1715: else PetscCall(PetscInfo(pc, "PCGAMGProlongatorBlockFilter_AGG: skipping zero elimination for %s; filtered entries are zeroed but not removed\n", ((PetscObject)Prol)->type_name));
1716: PetscFunctionReturn(PETSC_SUCCESS);
1717: }
1719: // 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
1720: static PetscErrorCode PCGAMGKernelPreservingFilter_AGG(PC pc, Mat Prol, PetscReal threshold)
1721: {
1722: PC_MG *mg = (PC_MG *)pc->data;
1723: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1724: const PetscInt nSAvec = pc_gamg->data_cell_rows; /* == data_cell_cols after formProl0 */
1725: PetscInt cStart, cEnd, rStart, rEnd;
1726: const PetscReal *Bc_data = pc_gamg->data;
1727: Vec *Bc_vecs, *B_vecs;
1728: PetscScalar *Bc_arr;
1730: PetscFunctionBegin;
1731: PetscCall(PetscInfo(pc, "Kernel-preserving filter of prolongator with threshold %g, nSAvec=%" PetscInt_FMT "\n", (double)threshold, nSAvec));
1733: PetscCall(MatGetOwnershipRange(Prol, &rStart, &rEnd));
1734: PetscCall(MatGetOwnershipRangeColumn(Prol, &cStart, &cEnd));
1736: /* Step 1: build coarse null-space vectors and compute B = P_original * B_c */
1737: PetscCall(PetscMalloc1(nSAvec, &Bc_vecs));
1738: PetscCall(PetscMalloc1(nSAvec, &B_vecs));
1740: {
1741: PetscInt nloc = cEnd - cStart;
1742: for (PetscInt k = 0; k < nSAvec; k++) {
1743: PetscCall(MatCreateVecs(Prol, &Bc_vecs[k], &B_vecs[k]));
1744: /* fill local entries: Bc_data layout is Bc_data[k * nloc + c] (stride == nloc) */
1745: PetscCall(VecGetArray(Bc_vecs[k], &Bc_arr));
1746: for (PetscInt c = 0; c < nloc; c++) Bc_arr[c] = (PetscScalar)Bc_data[k * nloc + c];
1747: PetscCall(VecRestoreArray(Bc_vecs[k], &Bc_arr));
1748: PetscCall(MatMult(Prol, Bc_vecs[k], B_vecs[k]));
1749: }
1750: }
1752: /* Step 2: apply the block-aware threshold filter (drops whole node-coupling
1753: blocks, preserving the coarse-node block structure; see
1754: PCGAMGProlongatorBlockFilter_AGG()) */
1755: {
1756: PetscBool info_active = PETSC_FALSE;
1757: MatInfo info0, info1;
1758: PetscCall(PetscInfoEnabled(((PetscObject)pc)->classid, &info_active));
1759: if (info_active) PetscCall(MatGetInfo(Prol, MAT_GLOBAL_SUM, &info0));
1760: PetscCall(PCGAMGProlongatorBlockFilter_AGG(pc, Prol, nSAvec, threshold));
1761: if (info_active) {
1762: PetscCall(MatGetInfo(Prol, MAT_GLOBAL_SUM, &info1));
1763: 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));
1764: }
1765: }
1767: /* Step 3: correct rows to restore P_filtered * B_c = B */
1768: if (nSAvec == 1) {
1769: /*
1770: Scalar case: use `MatMult()` + element-wise scaling + `MatDiagonalScale()`.
1771: scale_i = B_i / (P_filtered * Bc)_i, then P_new = diag(scale) * P_filtered.
1772: A row is left unscaled, that is, the near-null space constraint is not enforced
1773: for it, in the two cases where scaling would do more harm than the constraint is
1774: worth: when |scale_i| would exceed smax, which covers a denominator that is zero
1775: (an empty row after the filter) or tiny relative to B_i, and when B_i is zero,
1776: where enforcing the constraint would zero the whole row of the prolongator.
1777: Zeros left in the sparsity pattern stay zero under `MatDiagonalScale()`.
1778: No ghost column access needed.
1779: */
1780: Vec d_vec, scale_vec;
1781: PetscInt n_local, n_unscaled = 0;
1782: PetscReal smax = 1.0e4; /* cap on the row scale; the filter drops small blocks, so a healthy row has scale ~1 */
1783: PetscScalar *s_arr;
1784: const PetscScalar *b_arr, *d_arr;
1786: PetscCall(MatCreateVecs(Prol, NULL, &d_vec));
1787: PetscCall(MatMult(Prol, Bc_vecs[0], d_vec));
1788: PetscCall(VecDuplicate(d_vec, &scale_vec));
1789: PetscCall(VecGetLocalSize(d_vec, &n_local));
1790: PetscCall(VecGetArrayRead(B_vecs[0], &b_arr));
1791: PetscCall(VecGetArrayRead(d_vec, &d_arr));
1792: PetscCall(VecGetArray(scale_vec, &s_arr));
1793: for (PetscInt i = 0; i < n_local; i++) {
1794: PetscReal b = PetscAbsScalar(b_arr[i]), d = PetscAbsScalar(d_arr[i]);
1796: if (b > 0.0 && smax * d > b) s_arr[i] = b_arr[i] / d_arr[i];
1797: else {
1798: s_arr[i] = 1.0;
1799: n_unscaled++;
1800: }
1801: }
1802: 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));
1803: PetscCall(VecRestoreArray(scale_vec, &s_arr));
1804: PetscCall(VecRestoreArrayRead(d_vec, &d_arr));
1805: PetscCall(VecRestoreArrayRead(B_vecs[0], &b_arr));
1806: PetscCall(MatDiagonalScale(Prol, scale_vec, NULL));
1807: PetscCall(VecDestroy(&d_vec));
1808: PetscCall(VecDestroy(&scale_vec));
1809: } else {
1810: /*
1811: Vector case (nSAvec > 1): per-row minimum-norm correction (Gram-matrix solve).
1812: Scatter Bc_data to include ghost column values using Prol's Mvctx,
1813: then build a hash map from global ghost column index to local ghost index
1814: so that `MatGetRow()` global column indices can be mapped to the ghosted array.
1815: */
1816: PetscInt nloc = cEnd - cStart;
1817: PetscInt ghost_stride;
1818: PetscReal *Bc_ghosted = NULL;
1819: const PetscReal *Bc_ghosted_ro;
1820: PetscMPIInt comm_size;
1821: PetscHMapI ghost_gid_to_lid; /* global ghost col index -> local ghost index (0-based) */
1822: PetscInt num_ghosts = 0;
1823: PetscBool no_off_proc;
1825: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)Prol), &comm_size));
1826: if (comm_size > 1) {
1827: Mat_MPIAIJ *mpimat = (Mat_MPIAIJ *)Prol->data;
1828: Vec tmp_vec;
1829: PetscScalar *data_arr;
1830: PetscInt nnodes;
1831: PetscBool isaij;
1833: PetscCall(PetscObjectBaseTypeCompare((PetscObject)Prol, MATMPIAIJ, &isaij));
1834: PetscCheck(isaij, PetscObjectComm((PetscObject)Prol), PETSC_ERR_SUP, "Prolongator filter requires an MPIAIJ-based prolongator, not %s", ((PetscObject)Prol)->type_name);
1835: PetscCall(VecGetLocalSize(mpimat->lvec, &num_ghosts));
1836: nnodes = nloc + num_ghosts;
1837: ghost_stride = nnodes;
1838: /*
1839: Scatter Bc_data to include ghost column values using Prol's Mvctx.
1840: Cannot use PCGAMGGetDataWithGhosts() because it assumes square matrix
1841: (uses MatGetOwnershipRange() for row indices, but Prol is rectangular).
1842: */
1843: PetscCall(MatCreateVecs(Prol, &tmp_vec, NULL));
1844: PetscCall(PetscMalloc1(nSAvec * nnodes, &Bc_ghosted));
1845: for (PetscInt dir = 0; dir < nSAvec; dir++) {
1846: PetscScalar *tmp_arr;
1847: PetscCall(VecGetArray(tmp_vec, &tmp_arr));
1848: for (PetscInt kk = 0; kk < nloc; kk++) {
1849: PetscReal val = Bc_data[dir * nloc + kk];
1850: Bc_ghosted[dir * nnodes + kk] = val;
1851: tmp_arr[kk] = (PetscScalar)val;
1852: }
1853: PetscCall(VecRestoreArray(tmp_vec, &tmp_arr));
1854: PetscCall(VecScatterBegin(mpimat->Mvctx, tmp_vec, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD));
1855: PetscCall(VecScatterEnd(mpimat->Mvctx, tmp_vec, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD));
1856: PetscCall(VecGetArray(mpimat->lvec, &data_arr));
1857: for (PetscInt g = 0; g < num_ghosts; g++) Bc_ghosted[dir * nnodes + nloc + g] = PetscRealPart(data_arr[g]);
1858: PetscCall(VecRestoreArray(mpimat->lvec, &data_arr));
1859: }
1860: PetscCall(VecDestroy(&tmp_vec));
1861: Bc_ghosted_ro = Bc_ghosted;
1862: /* build hash: global ghost col index -> local ghost index (0-based into ghost portion) */
1863: PetscCall(PetscHMapICreateWithSize(2 * num_ghosts + 1, &ghost_gid_to_lid));
1864: for (PetscInt g = 0; g < num_ghosts; g++) PetscCall(PetscHMapISet(ghost_gid_to_lid, mpimat->garray[g], g));
1865: } else {
1866: /* sequential: no ghosts, ghost_stride == nloc, use Bc_data directly (read-only) */
1867: ghost_stride = nloc;
1868: Bc_ghosted_ro = Bc_data;
1869: PetscCall(PetscHMapICreateWithSize(1, &ghost_gid_to_lid));
1870: }
1872: {
1873: PetscInt nrows = rEnd - rStart, max_ncols = 0;
1874: const PetscScalar **B_arrays;
1875: PetscScalar *work, *new_vals, *G, *rhs, *x, *bc_col;
1876: PetscReal *dscale;
1877: PetscInt *ghosted_idx, *act, *col_buf;
1878: PetscBLASInt *ipiv;
1879: PetscBLASInt N_b;
1881: PetscCall(PetscMalloc1(nSAvec, &B_arrays));
1882: for (PetscInt k = 0; k < nSAvec; k++) PetscCall(VecGetArrayRead(B_vecs[k], &B_arrays[k]));
1883: /* work: nSAvec*nSAvec Gram + nSAvec rhs + nSAvec solution + nSAvec bc_col scratch */
1884: PetscCall(PetscMalloc1(nSAvec * nSAvec + 3 * nSAvec, &work));
1885: PetscCall(PetscMalloc1(nSAvec, &dscale));
1886: PetscCall(PetscMalloc1(nSAvec, &ipiv));
1887: PetscCall(PetscBLASIntCast(nSAvec, &N_b));
1888: G = work;
1889: rhs = work + nSAvec * nSAvec;
1890: x = rhs + nSAvec;
1891: bc_col = x + nSAvec;
1893: /* find max row width and total nnz for pre-allocation */
1894: {
1895: PetscInt total_nnz = 0;
1896: for (PetscInt row = 0; row < nrows; row++) {
1897: PetscInt ncols;
1898: PetscCall(MatGetRow(Prol, rStart + row, &ncols, NULL, NULL));
1899: if (ncols > max_ncols) max_ncols = ncols;
1900: total_nnz += ncols;
1901: PetscCall(MatRestoreRow(Prol, rStart + row, &ncols, NULL, NULL));
1902: }
1903: /* allocate flat CSR-like buffers to store all corrections before applying */
1904: PetscCall(PetscMalloc1(total_nnz, &new_vals));
1905: PetscCall(PetscMalloc1(total_nnz, &col_buf));
1906: }
1907: PetscCall(PetscMalloc2(max_ncols, &ghosted_idx, max_ncols, &act));
1909: /* Pass 1: read rows, compute corrections, store in flat buffers */
1910: {
1911: PetscInt *row_offsets;
1912: PetscInt offset = 0, n_singular = 0, n_zero_rows = 0, n_corrected = 0, n_underdetermined = 0;
1913: PetscReal max_ynorm = 0.0;
1915: PetscCall(PetscMalloc1(nrows + 1, &row_offsets));
1916: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1918: for (PetscInt row = 0; row < nrows; row++) {
1919: PetscInt ncols, nact = 0, grow = rStart + row, roff = offset;
1920: PetscBLASInt NRHS = 1, LDA = N_b, LDB = N_b, info;
1921: const PetscInt *cols;
1922: const PetscScalar *vals;
1924: row_offsets[row] = roff;
1925: PetscCall(MatGetRow(Prol, grow, &ncols, &cols, &vals));
1926: /* Save the row unchanged, then correct only its nonzero entries: the block filter zeroes
1927: dropped entries but they are not removed from the sparsity pattern for every matrix type
1928: (see PCGAMGProlongatorBlockFilter_AGG()), and adding the correction to them would turn a
1929: dropped block back into a nonzero one */
1930: for (PetscInt j = 0; j < ncols; j++) {
1931: col_buf[roff + j] = cols[j];
1932: new_vals[roff + j] = vals[j];
1933: if (vals[j] != 0.0) act[nact++] = j;
1934: }
1935: offset = roff + ncols;
1936: if (nact == 0) {
1937: n_zero_rows++;
1938: PetscCall(MatRestoreRow(Prol, grow, &ncols, &cols, &vals));
1939: continue;
1940: }
1942: /* When nact < nSAvec the Gram matrix G is rank-deficient by construction;
1943: skip correction for this row (keep filtered values as-is).
1944: Note: the near-null space constraint P*Bc = B is NOT enforced for these rows.
1945: This typically occurs at boundary or isolated nodes where few coarse neighbors
1946: remain after filtering; the impact on convergence is generally small. */
1947: if (nact < nSAvec) {
1948: n_underdetermined++;
1949: PetscCall(MatRestoreRow(Prol, grow, &ncols, &cols, &vals));
1950: continue;
1951: }
1953: /* map the global column indices of the surviving entries to ghosted array indices */
1954: for (PetscInt a = 0; a < nact; a++) {
1955: PetscInt col = cols[act[a]];
1956: if (col >= cStart && col < cEnd) ghosted_idx[a] = col - cStart;
1957: else {
1958: PetscInt g = -1;
1959: PetscCall(PetscHMapIGet(ghost_gid_to_lid, col, &g));
1960: PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Off-diagonal column %" PetscInt_FMT " not found in ghost map for prolongator filter", col);
1961: ghosted_idx[a] = nloc + g;
1962: }
1963: }
1965: for (PetscInt i = 0; i < nSAvec * nSAvec; i++) G[i] = 0.0;
1967: /* rhs[k] = B[row,k] - sum_a P[row,act[a]] * Bc[ghosted_idx[a], k] */
1968: for (PetscInt k = 0; k < nSAvec; k++) {
1969: PetscScalar dot = 0.0;
1970: for (PetscInt a = 0; a < nact; a++) dot += vals[act[a]] * (PetscScalar)Bc_ghosted_ro[k * ghost_stride + ghosted_idx[a]];
1971: rhs[k] = B_arrays[k][row] - dot;
1972: }
1974: /* G[k1,k2] = sum_a Bc[a,k1] * Bc[a,k2] using pre-gathered bc_col */
1975: for (PetscInt a = 0; a < nact; a++) {
1976: PetscInt gidx = ghosted_idx[a];
1977: for (PetscInt k = 0; k < nSAvec; k++) bc_col[k] = (PetscScalar)Bc_ghosted_ro[k * ghost_stride + gidx];
1978: for (PetscInt k1 = 0; k1 < nSAvec; k1++)
1979: for (PetscInt k2 = k1; k2 < nSAvec; k2++) G[k1 * nSAvec + k2] += bc_col[k1] * bc_col[k2];
1980: }
1981: /* fill lower triangle from upper (G is symmetric) */
1982: for (PetscInt k1 = 1; k1 < nSAvec; k1++)
1983: for (PetscInt k2 = 0; k2 < k1; k2++) G[k1 * nSAvec + k2] = G[k2 * nSAvec + k1];
1985: /* Symmetric (Jacobi) equilibration: G is severely ill-conditioned for
1986: elasticity because the near-null modes have disparate scales (O(1)
1987: translations vs rotations that scale with the coordinates). Scale by
1988: dscale[k] = 1/sqrt(G[k,k]) so the rescaled Gram has a unit diagonal,
1989: then solve (D G D) y = D rhs and recover x = D y. This is identical to
1990: solving G x = rhs in exact arithmetic but removes the mode-scale
1991: ill-conditioning, making the correction accurate and FP-robust. */
1992: for (PetscInt k = 0; k < nSAvec; k++) {
1993: PetscReal gkk = PetscRealPart(G[k * nSAvec + k]);
1994: dscale[k] = gkk > 0.0 ? 1.0 / PetscSqrtReal(gkk) : 1.0;
1995: }
1996: for (PetscInt k1 = 0; k1 < nSAvec; k1++)
1997: for (PetscInt k2 = 0; k2 < nSAvec; k2++) G[k1 * nSAvec + k2] *= dscale[k1] * dscale[k2];
1999: /* solve (D G D) y = D rhs, then x = D y */
2000: for (PetscInt k = 0; k < nSAvec; k++) x[k] = dscale[k] * rhs[k];
2001: PetscCallBLAS("LAPACKgesv", LAPACKgesv_(&N_b, &NRHS, G, &LDA, ipiv, x, &LDB, &info));
2002: PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "LAPACKgesv: %" PetscBLASInt_FMT "-th argument had an illegal value", -info);
2003: if (info > 0) {
2004: /* G is singular despite nact >= nSAvec (Bc columns linearly dependent);
2005: keep filtered values as-is (near-null space constraint not enforced for this row) */
2006: n_singular++;
2007: PetscCall(MatRestoreRow(Prol, grow, &ncols, &cols, &vals));
2008: continue;
2009: }
2010: /* track the equilibrated solution norm ||y||^2 (x still holds y here): a
2011: basis-independent health metric, large values flag a problematic row */
2012: {
2013: PetscReal ynorm2 = 0.0;
2014: for (PetscInt k = 0; k < nSAvec; k++) ynorm2 += PetscSqr(PetscAbsScalar(x[k]));
2015: if (ynorm2 > max_ynorm) max_ynorm = ynorm2;
2016: }
2017: for (PetscInt k = 0; k < nSAvec; k++) x[k] *= dscale[k]; /* recover x = D y */
2018: n_corrected++;
2020: /* new_vals[act[a]] = vals[act[a]] + sum_k Bc[ghosted_idx[a],k] * x[k] */
2021: for (PetscInt a = 0; a < nact; a++) {
2022: PetscScalar delta = 0.0;
2023: PetscInt gidx = ghosted_idx[a];
2024: for (PetscInt k = 0; k < nSAvec; k++) delta += (PetscScalar)Bc_ghosted_ro[k * ghost_stride + gidx] * x[k];
2025: new_vals[roff + act[a]] += delta;
2026: }
2027: PetscCall(MatRestoreRow(Prol, grow, &ncols, &cols, &vals));
2028: }
2029: row_offsets[nrows] = offset;
2030: PetscCall(PetscFPTrapPop());
2031: PetscCall(PetscInfo(pc, "PCGAMGKernelPreservingFilter_AGG: corrected %" PetscInt_FMT "/%" PetscInt_FMT " rows, max equilibrated correction ||y||^2=%g\n", n_corrected, nrows, (double)max_ynorm));
2032: if (n_zero_rows + n_underdetermined + n_singular > 0)
2033: 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));
2035: /* Pass 2: apply all corrections at once */
2036: for (PetscInt row = 0; row < nrows; row++) {
2037: PetscInt grow = rStart + row;
2038: PetscInt nc = row_offsets[row + 1] - row_offsets[row];
2039: if (nc > 0) PetscCall(MatSetValues(Prol, 1, &grow, nc, col_buf + row_offsets[row], new_vals + row_offsets[row], INSERT_VALUES));
2040: }
2041: PetscCall(PetscFree(row_offsets));
2042: }
2044: for (PetscInt k = 0; k < nSAvec; k++) PetscCall(VecRestoreArrayRead(B_vecs[k], &B_arrays[k]));
2045: PetscCall(PetscFree(B_arrays));
2046: PetscCall(PetscFree(work));
2047: PetscCall(PetscFree(dscale));
2048: PetscCall(PetscFree(ipiv));
2049: PetscCall(PetscFree2(ghosted_idx, act));
2050: PetscCall(PetscFree(new_vals));
2051: PetscCall(PetscFree(col_buf));
2052: }
2054: PetscCall(PetscHMapIDestroy(&ghost_gid_to_lid));
2055: if (comm_size > 1) PetscCall(PetscFree(Bc_ghosted));
2057: /* all insertions are in local rows; skip the off-process assembly communication. The scalar
2058: branch needs no assembly at all: MatDiagonalScale() requires and preserves assembly */
2059: PetscCall(MatGetOption(Prol, MAT_NO_OFF_PROC_ENTRIES, &no_off_proc));
2060: PetscCall(MatSetOption(Prol, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
2061: PetscCall(MatAssemblyBegin(Prol, MAT_FINAL_ASSEMBLY));
2062: PetscCall(MatAssemblyEnd(Prol, MAT_FINAL_ASSEMBLY));
2063: PetscCall(MatSetOption(Prol, MAT_NO_OFF_PROC_ENTRIES, no_off_proc));
2064: }
2066: for (PetscInt k = 0; k < nSAvec; k++) {
2067: PetscCall(VecDestroy(&Bc_vecs[k]));
2068: PetscCall(VecDestroy(&B_vecs[k]));
2069: }
2070: PetscCall(PetscFree(Bc_vecs));
2071: PetscCall(PetscFree(B_vecs));
2072: PetscFunctionReturn(PETSC_SUCCESS);
2073: }
2075: /*
2076: PCGAMGOptimizeProlongator_AGG - given the initial prolongator optimizes it by smoothed aggregation pc_gamg_agg->nsmooths times
2078: Input Parameter:
2079: . pc - this
2080: . Amat - matrix on this fine level
2081: In/Output Parameter:
2082: . a_P - prolongation operator to the next level
2083: */
2084: static PetscErrorCode PCGAMGOptimizeProlongator_AGG(PC pc, Mat Amat, Mat *a_P)
2085: {
2086: PC_MG *mg = (PC_MG *)pc->data;
2087: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
2088: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
2089: Mat Prol = *a_P;
2090: MPI_Comm comm;
2091: KSP eksp;
2092: Vec bb, xx;
2093: PC epc;
2094: PetscReal alpha, emax, emin;
2095: PetscReal pfilter = pc_gamg->prolongator_filter * PetscPowRealInt(pc_gamg->prolongator_filter_scale, pc_gamg->current_level);
2097: PetscFunctionBegin;
2098: PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
2099: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0));
2101: /* compute maximum singular value of operator to be used in smoother */
2102: if (0 < pc_gamg_agg->nsmooths) {
2103: /* get eigen estimates */
2104: if (pc_gamg->emax > 0) {
2105: emin = pc_gamg->emin;
2106: emax = pc_gamg->emax;
2107: } else {
2108: const char *prefix;
2110: PetscCall(MatCreateVecs(Amat, &bb, NULL));
2111: PetscCall(MatCreateVecs(Amat, &xx, NULL));
2112: PetscCall(KSPSetNoisy_Private(Amat, bb));
2114: PetscCall(KSPCreate(comm, &eksp));
2115: PetscCall(KSPSetNestLevel(eksp, pc->kspnestlevel));
2116: PetscCall(PCGetOptionsPrefix(pc, &prefix));
2117: PetscCall(KSPSetOptionsPrefix(eksp, prefix));
2118: PetscCall(KSPAppendOptionsPrefix(eksp, "pc_gamg_esteig_"));
2119: {
2120: PetscBool isset, sflg;
2122: PetscCall(MatIsSPDKnown(Amat, &isset, &sflg));
2123: if (isset && sflg) PetscCall(KSPSetType(eksp, KSPCG));
2124: }
2125: PetscCall(KSPSetErrorIfNotConverged(eksp, pc->erroriffailure));
2126: PetscCall(KSPSetNormType(eksp, KSP_NORM_NONE));
2128: PetscCall(KSPSetInitialGuessNonzero(eksp, PETSC_FALSE));
2129: PetscCall(KSPSetOperators(eksp, Amat, Amat));
2131: PetscCall(KSPGetPC(eksp, &epc));
2132: PetscCall(PCSetType(epc, PCJACOBI)); /* smoother in smoothed agg. */
2134: 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
2136: PetscCall(KSPSetFromOptions(eksp));
2137: PetscCall(KSPSetComputeSingularValues(eksp, PETSC_TRUE));
2138: PetscCall(KSPSolve(eksp, bb, xx));
2139: PetscCall(KSPCheckSolve(eksp, pc, xx));
2141: PetscCall(KSPComputeExtremeSingularValues(eksp, &emax, &emin));
2142: PetscCall(PetscInfo(pc, "%s: Smooth P0: max eigen=%e min=%e PC=%s\n", ((PetscObject)pc)->prefix, (double)emax, (double)emin, PCJACOBI));
2143: PetscCall(VecDestroy(&xx));
2144: PetscCall(VecDestroy(&bb));
2145: PetscCall(KSPDestroy(&eksp));
2146: }
2147: if (pc_gamg->use_sa_esteig) {
2148: mg->min_eigen_DinvA[pc_gamg->current_level] = emin;
2149: mg->max_eigen_DinvA[pc_gamg->current_level] = emax;
2150: 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));
2151: } else {
2152: mg->min_eigen_DinvA[pc_gamg->current_level] = 0;
2153: mg->max_eigen_DinvA[pc_gamg->current_level] = 0;
2154: }
2155: } else {
2156: mg->min_eigen_DinvA[pc_gamg->current_level] = 0;
2157: mg->max_eigen_DinvA[pc_gamg->current_level] = 0;
2158: }
2160: /* smooth P0 */
2161: if (pc_gamg_agg->nsmooths > 0) {
2162: Vec diag;
2164: /* TODO: Set a PCFailedReason and exit the building of the AMG preconditioner */
2165: PetscCheck(emax != 0.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Computed maximum singular value as zero");
2167: PetscCall(MatCreateVecs(Amat, &diag, NULL));
2168: PetscCall(MatGetDiagonal(Amat, diag)); /* effectively PCJACOBI */
2169: PetscCall(VecReciprocal(diag));
2171: for (PetscInt jj = 0; jj < pc_gamg_agg->nsmooths; jj++) {
2172: Mat tMat;
2174: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0));
2175: /*
2176: Smooth aggregation on the prolongator
2178: P_{i} := (I - 1.4/emax D^{-1}A) P_i\{i-1}
2179: */
2180: PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0));
2181: PetscCall(MatMatMult(Amat, Prol, MAT_INITIAL_MATRIX, PETSC_CURRENT, &tMat));
2182: PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0));
2183: PetscCall(MatProductClear(tMat));
2184: PetscCall(MatDiagonalScale(tMat, diag, NULL));
2186: /* TODO: Document the 1.4 and don't hardwire it in this routine */
2187: alpha = -1.4 / emax;
2188: PetscCall(MatAYPX(tMat, alpha, Prol, SUBSET_NONZERO_PATTERN));
2189: PetscCall(MatDestroy(&Prol));
2190: Prol = tMat;
2191: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0));
2192: }
2193: PetscCall(VecDestroy(&diag));
2194: }
2195: /* a per-level threshold of 0 (from prolongator_filter_scale == 0 on the coarser levels) drops
2196: nothing, so skip the whole filter rather than pay for its passes and per-row solves */
2197: if (pfilter > 0.0) {
2198: 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,
2199: (double)pc_gamg->prolongator_filter_scale, pc_gamg->current_level));
2200: PetscCall(PCGAMGKernelPreservingFilter_AGG(pc, Prol, pfilter));
2201: }
2202: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0));
2203: PetscCall(MatViewFromOptions(Prol, NULL, "-pc_gamg_agg_view_prolongation"));
2204: *a_P = Prol;
2205: PetscFunctionReturn(PETSC_SUCCESS);
2206: }
2208: /*MC
2209: PCGAMGAGG - Smooth aggregation, {cite}`vanek1996algebraic`, {cite}`vanek2001convergence`, variant of PETSc's algebraic multigrid (`PCGAMG`) preconditioner
2211: Options Database Keys:
2212: + -pc_gamg_agg_nsmooths nsmooth - number of smoothing steps to use with smooth aggregation to construct prolongation
2213: . -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)
2214: . -pc_gamg_prolongator_filter_scale scale - per-level scaling of the prolongator filter threshold (1.0=default)
2215: . -pc_gamg_aggressive_coarsening n - number of aggressive coarsening (MIS-2 or square graph) levels from finest.
2216: . -pc_gamg_aggressive_square_graph (true|false) - use square graph ($A^T A$), alternative is MIS-k (k=2), for aggressive coarsening
2217: . -pc_gamg_mis_k_minimum_degree_ordering (true|false) - use minimum degree ordering in greedy MIS algorithm
2218: . -pc_gamg_asm_hem_aggs n - number of HEM aggregation steps for ASM smoother
2219: - -pc_gamg_aggressive_mis_k n - number (k) distance in MIS coarsening (>2 is 'aggressive')
2221: Level: intermediate
2223: Notes:
2224: To obtain good performance for `PCGAMG` for vector valued problems you must
2225: call `MatSetBlockSize()` to indicate the number of degrees of freedom per grid point.
2226: Call `MatSetNearNullSpace()` (or `PCSetCoordinates()` if solving the equations of elasticity) to indicate the near null space of the operator
2228: When `-pc_gamg_aggressive_square_graph` is used, the coarsening is obtained by first squaring the graph and then applying, by default, a
2229: MIS-1 coarsening with `MatCoarsenApply()` on the squared graph.
2231: The many options for `PCMG` and `PCGAMG` such as controlling the smoothers on each level etc. also work for `PCGAMGAGG`
2233: .seealso: `PCGAMG`, [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCCreate()`, `PCSetType()`,
2234: `MatSetBlockSize()`, `PCMGType`, `PCSetCoordinates()`, `MatSetNearNullSpace()`, `PCGAMGSetType()`,
2235: `PCGAMGAGG`, `PCGAMGGEO`, `PCGAMGCLASSICAL`, `PCGAMGSetProcEqLim()`, `PCGAMGSetCoarseEqLim()`, `PCGAMGSetRepartition()`, `PCGAMGRegister()`,
2236: `PCGAMGSetReuseInterpolation()`, `PCGAMGASMSetUseAggs()`, `PCGAMGSetParallelCoarseGridSolve()`, `PCGAMGSetNlevels()`, `PCGAMGSetThreshold()`,
2237: `PCGAMGGetType()`, `PCGAMGSetUseSAEstEig()`
2238: M*/
2239: PetscErrorCode PCCreateGAMG_AGG(PC pc)
2240: {
2241: PC_MG *mg = (PC_MG *)pc->data;
2242: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
2243: PC_GAMG_AGG *pc_gamg_agg;
2245: PetscFunctionBegin;
2246: /* create sub context for SA */
2247: PetscCall(PetscNew(&pc_gamg_agg));
2248: pc_gamg->subctx = pc_gamg_agg;
2250: pc_gamg->ops->setfromoptions = PCSetFromOptions_GAMG_AGG;
2251: pc_gamg->ops->destroy = PCDestroy_GAMG_AGG;
2252: /* reset does not do anything; setup not virtual */
2254: /* set internal function pointers */
2255: pc_gamg->ops->creategraph = PCGAMGCreateGraph_AGG;
2256: pc_gamg->ops->coarsen = PCGAMGCoarsen_AGG;
2257: pc_gamg->ops->prolongator = PCGAMGConstructProlongator_AGG;
2258: pc_gamg->ops->optprolongator = PCGAMGOptimizeProlongator_AGG;
2259: pc_gamg->ops->createdefaultdata = PCSetData_AGG;
2260: pc_gamg->ops->view = PCView_GAMG_AGG;
2262: pc_gamg_agg->nsmooths = 1;
2263: pc_gamg_agg->aggressive_coarsening_levels = 1;
2264: pc_gamg_agg->use_aggressive_square_graph = PETSC_TRUE;
2265: pc_gamg_agg->use_minimum_degree_ordering = PETSC_FALSE;
2266: pc_gamg_agg->use_low_mem_filter = PETSC_FALSE;
2267: pc_gamg_agg->aggressive_mis_k = 2;
2268: pc_gamg_agg->graph_symmetrize = PETSC_TRUE;
2270: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", PCGAMGSetNSmooths_AGG));
2271: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", PCGAMGSetAggressiveLevels_AGG));
2272: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveSquareGraph_C", PCGAMGSetAggressiveSquareGraph_AGG));
2273: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetMinDegreeOrdering_C", PCGAMGMISkSetMinDegreeOrdering_AGG));
2274: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetLowMemoryFilter_C", PCGAMGSetLowMemoryFilter_AGG));
2275: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetAggressive_C", PCGAMGMISkSetAggressive_AGG));
2276: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetGraphSymmetrize_C", PCGAMGSetGraphSymmetrize_AGG));
2277: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetProlongatorFilter_C", PCGAMGSetProlongatorFilter_AGG));
2278: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGGetProlongatorFilter_C", PCGAMGGetProlongatorFilter_AGG));
2279: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetProlongatorFilterScale_C", PCGAMGSetProlongatorFilterScale_AGG));
2280: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGGetProlongatorFilterScale_C", PCGAMGGetProlongatorFilterScale_AGG));
2281: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_AGG));
2282: PetscFunctionReturn(PETSC_SUCCESS);
2283: }