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, they are
251: zeroed but remain in the pattern, so the coarse operators are unchanged in structure and the complexity and memory reduction is not realized (reported with `-info`).
253: 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
254: 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
255: 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
256: 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
257: 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).
259: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGGetProlongatorFilter()`, `PCGAMGSetProlongatorFilterScale()`, `PCGAMGSetLowMemoryFilter()`
260: @*/
261: PetscErrorCode PCGAMGSetProlongatorFilter(PC pc, PetscReal thr)
262: {
263: PetscFunctionBegin;
266: PetscTryMethod(pc, "PCGAMGSetProlongatorFilter_C", (PC, PetscReal), (pc, thr));
267: PetscFunctionReturn(PETSC_SUCCESS);
268: }
270: /*@
271: PCGAMGGetProlongatorFilter - Get the relative threshold for block filtering of the prolongator
273: Not Collective
275: Input Parameter:
276: . pc - the preconditioner context
278: Output Parameter:
279: . 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
280: the largest such block norm in that fine node's block row (0 disables filtering, see `PCGAMGSetProlongatorFilter()`)
282: Level: intermediate
284: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetProlongatorFilter()`, `PCGAMGSetProlongatorFilterScale()`, `PCGAMGSetLowMemoryFilter()`
285: @*/
286: PetscErrorCode PCGAMGGetProlongatorFilter(PC pc, PetscReal *thr)
287: {
288: PetscFunctionBegin;
290: PetscAssertPointer(thr, 2);
291: PetscUseMethod(pc, "PCGAMGGetProlongatorFilter_C", (PC, PetscReal *), (pc, thr));
292: PetscFunctionReturn(PETSC_SUCCESS);
293: }
295: /*@
296: PCGAMGSetProlongatorFilterScale - Set the per-level scaling of the prolongator filter threshold (see `PCGAMGSetProlongatorFilter()`)
298: Logically Collective
300: Input Parameters:
301: + pc - the preconditioner context
302: - 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
304: Options Database Key:
305: . -pc_gamg_prolongator_filter_scale scale - per-level scaling of the prolongator filter threshold (1.0=default)
307: Level: intermediate
309: Note:
310: 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
311: aggressively than the finest, are not allowed. A scale of 0 disables filtering on all levels but the finest.
313: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetProlongatorFilter()`, `PCGAMGGetProlongatorFilterScale()`
314: @*/
315: PetscErrorCode PCGAMGSetProlongatorFilterScale(PC pc, PetscReal scale)
316: {
317: PetscFunctionBegin;
320: PetscTryMethod(pc, "PCGAMGSetProlongatorFilterScale_C", (PC, PetscReal), (pc, scale));
321: PetscFunctionReturn(PETSC_SUCCESS);
322: }
324: /*@
325: PCGAMGGetProlongatorFilterScale - Get the per-level scaling of the prolongator filter threshold
327: Not Collective
329: Input Parameter:
330: . pc - the preconditioner context
332: Output Parameter:
333: . 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
335: Level: intermediate
337: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetProlongatorFilter()`, `PCGAMGSetProlongatorFilterScale()`
338: @*/
339: PetscErrorCode PCGAMGGetProlongatorFilterScale(PC pc, PetscReal *scale)
340: {
341: PetscFunctionBegin;
343: PetscAssertPointer(scale, 2);
344: PetscUseMethod(pc, "PCGAMGGetProlongatorFilterScale_C", (PC, PetscReal *), (pc, scale));
345: PetscFunctionReturn(PETSC_SUCCESS);
346: }
348: static PetscErrorCode PCGAMGSetAggressiveLevels_AGG(PC pc, PetscInt n)
349: {
350: PC_MG *mg = (PC_MG *)pc->data;
351: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
352: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
354: PetscFunctionBegin;
355: pc_gamg_agg->aggressive_coarsening_levels = n;
356: PetscFunctionReturn(PETSC_SUCCESS);
357: }
359: static PetscErrorCode PCGAMGMISkSetAggressive_AGG(PC pc, PetscInt n)
360: {
361: PC_MG *mg = (PC_MG *)pc->data;
362: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
363: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
365: PetscFunctionBegin;
366: pc_gamg_agg->aggressive_mis_k = n;
367: PetscFunctionReturn(PETSC_SUCCESS);
368: }
370: static PetscErrorCode PCGAMGSetAggressiveSquareGraph_AGG(PC pc, PetscBool b)
371: {
372: PC_MG *mg = (PC_MG *)pc->data;
373: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
374: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
376: PetscFunctionBegin;
377: pc_gamg_agg->use_aggressive_square_graph = b;
378: PetscFunctionReturn(PETSC_SUCCESS);
379: }
381: static PetscErrorCode PCGAMGSetLowMemoryFilter_AGG(PC pc, PetscBool b)
382: {
383: PC_MG *mg = (PC_MG *)pc->data;
384: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
385: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
387: PetscFunctionBegin;
388: pc_gamg_agg->use_low_mem_filter = b;
389: PetscFunctionReturn(PETSC_SUCCESS);
390: }
392: static PetscErrorCode PCGAMGSetGraphSymmetrize_AGG(PC pc, PetscBool b)
393: {
394: PC_MG *mg = (PC_MG *)pc->data;
395: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
396: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
398: PetscFunctionBegin;
399: pc_gamg_agg->graph_symmetrize = b;
400: PetscFunctionReturn(PETSC_SUCCESS);
401: }
403: static PetscErrorCode PCGAMGMISkSetMinDegreeOrdering_AGG(PC pc, PetscBool b)
404: {
405: PC_MG *mg = (PC_MG *)pc->data;
406: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
407: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
409: PetscFunctionBegin;
410: pc_gamg_agg->use_minimum_degree_ordering = b;
411: PetscFunctionReturn(PETSC_SUCCESS);
412: }
414: static PetscErrorCode PCGAMGSetProlongatorFilter_AGG(PC pc, PetscReal thr)
415: {
416: PC_MG *mg = (PC_MG *)pc->data;
417: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
419: PetscFunctionBegin;
420: 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);
421: 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));
422: pc_gamg->prolongator_filter = thr;
423: PetscFunctionReturn(PETSC_SUCCESS);
424: }
426: static PetscErrorCode PCGAMGGetProlongatorFilter_AGG(PC pc, PetscReal *thr)
427: {
428: PC_MG *mg = (PC_MG *)pc->data;
429: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
431: PetscFunctionBegin;
432: *thr = pc_gamg->prolongator_filter;
433: PetscFunctionReturn(PETSC_SUCCESS);
434: }
436: static PetscErrorCode PCGAMGSetProlongatorFilterScale_AGG(PC pc, PetscReal scale)
437: {
438: PC_MG *mg = (PC_MG *)pc->data;
439: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
441: PetscFunctionBegin;
442: 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);
443: pc_gamg->prolongator_filter_scale = scale;
444: PetscFunctionReturn(PETSC_SUCCESS);
445: }
447: static PetscErrorCode PCGAMGGetProlongatorFilterScale_AGG(PC pc, PetscReal *scale)
448: {
449: PC_MG *mg = (PC_MG *)pc->data;
450: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
452: PetscFunctionBegin;
453: *scale = pc_gamg->prolongator_filter_scale;
454: PetscFunctionReturn(PETSC_SUCCESS);
455: }
457: static PetscErrorCode PCSetFromOptions_GAMG_AGG(PC pc, PetscOptionItems PetscOptionsObject)
458: {
459: PC_MG *mg = (PC_MG *)pc->data;
460: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
461: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
462: PetscBool n_aggressive_flg, old_sq_provided = PETSC_FALSE, new_sq_provided = PETSC_FALSE, new_sqr_graph = pc_gamg_agg->use_aggressive_square_graph;
463: PetscInt nsq_graph_old = 0;
464: PetscReal thr = pc_gamg->prolongator_filter;
465: PetscReal scale = pc_gamg->prolongator_filter_scale;
466: PetscBool flg;
468: PetscFunctionBegin;
469: PetscOptionsHeadBegin(PetscOptionsObject, "GAMG-AGG options");
470: 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));
471: // aggressive coarsening logic with deprecated -pc_gamg_square_graph
472: 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));
473: if (!n_aggressive_flg)
474: 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));
475: 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));
476: if (!new_sq_provided && old_sq_provided) {
477: pc_gamg_agg->aggressive_coarsening_levels = nsq_graph_old; // could be zero
478: pc_gamg_agg->use_aggressive_square_graph = PETSC_TRUE;
479: }
480: if (new_sq_provided && old_sq_provided)
481: 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));
482: 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));
483: 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));
484: 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));
485: PetscCall(PetscOptionsBool("-pc_gamg_graph_symmetrize", "Symmetrize graph for coarsening", "PCGAMGSetGraphSymmetrize", pc_gamg_agg->graph_symmetrize, &pc_gamg_agg->graph_symmetrize, NULL));
486: PetscCall(PetscOptionsBoundedReal("-pc_gamg_prolongator_filter", "Relative Frobenius-norm threshold for block filtering of the prolongator (0=disabled)", "PCGAMGSetProlongatorFilter", thr, &thr, &flg, 0.0));
487: if (flg) PetscCall(PCGAMGSetProlongatorFilter(pc, thr));
488: PetscCall(PetscOptionsRangeReal("-pc_gamg_prolongator_filter_scale", "Per-level scaling of the prolongator filter threshold", "PCGAMGSetProlongatorFilterScale", scale, &scale, &flg, 0.0, 1.0));
489: if (flg) PetscCall(PCGAMGSetProlongatorFilterScale(pc, scale));
491: PetscOptionsHeadEnd();
492: PetscFunctionReturn(PETSC_SUCCESS);
493: }
495: static PetscErrorCode PCDestroy_GAMG_AGG(PC pc)
496: {
497: PC_MG *mg = (PC_MG *)pc->data;
498: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
499: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
501: PetscFunctionBegin;
502: PetscCall(MatCoarsenDestroy(&pc_gamg_agg->crs));
503: PetscCall(PetscFree(pc_gamg->subctx));
504: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", NULL));
505: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", NULL));
506: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetAggressive_C", NULL));
507: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetMinDegreeOrdering_C", NULL));
508: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetLowMemoryFilter_C", NULL));
509: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveSquareGraph_C", NULL));
510: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetGraphSymmetrize_C", NULL));
511: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetProlongatorFilter_C", NULL));
512: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGGetProlongatorFilter_C", NULL));
513: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetProlongatorFilterScale_C", NULL));
514: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGGetProlongatorFilterScale_C", NULL));
515: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
516: PetscFunctionReturn(PETSC_SUCCESS);
517: }
519: /*
520: PCSetCoordinates_AGG
522: Collective
524: Input Parameter:
525: . pc - the preconditioner context
526: . ndm - dimension of data (used for dof/vertex for Stokes)
527: . a_nloc - number of vertices local
528: . coords - [a_nloc][ndm] - interleaved coordinate data: {x_0, y_0, z_0, x_1, y_1, ...}
529: */
531: static PetscErrorCode PCSetCoordinates_AGG(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords)
532: {
533: PC_MG *mg = (PC_MG *)pc->data;
534: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
535: PetscInt arrsz, kk, ii, jj, nloc, ndatarows, ndf;
536: Mat mat = pc->pmat;
538: PetscFunctionBegin;
541: nloc = a_nloc;
543: /* SA: null space vectors */
544: PetscCall(MatGetBlockSize(mat, &ndf)); /* this does not work for Stokes */
545: if (coords && ndf == 1) pc_gamg->data_cell_cols = 1; /* scalar w/ coords and SA (not needed) */
546: else if (coords) {
547: PetscCheck(ndm <= ndf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "degrees of motion %" PetscInt_FMT " > block size %" PetscInt_FMT, ndm, ndf);
548: pc_gamg->data_cell_cols = (ndm == 2 ? 3 : 6); /* displacement elasticity */
549: 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);
550: } else pc_gamg->data_cell_cols = ndf; /* no data, force SA with constant null space vectors */
551: pc_gamg->data_cell_rows = ndatarows = ndf;
552: 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);
553: arrsz = nloc * pc_gamg->data_cell_rows * pc_gamg->data_cell_cols;
555: if (!pc_gamg->data || (pc_gamg->data_sz != arrsz)) {
556: PetscCall(PetscFree(pc_gamg->data));
557: PetscCall(PetscMalloc1(arrsz + 1, &pc_gamg->data));
558: }
559: /* copy data in - column-oriented */
560: for (kk = 0; kk < nloc; kk++) {
561: const PetscInt M = nloc * pc_gamg->data_cell_rows; /* stride into data */
562: PetscReal *data = &pc_gamg->data[kk * ndatarows]; /* start of cell */
564: if (pc_gamg->data_cell_cols == 1) *data = 1.0;
565: else {
566: /* translational modes */
567: for (ii = 0; ii < ndatarows; ii++) {
568: for (jj = 0; jj < ndatarows; jj++) {
569: if (ii == jj) data[ii * M + jj] = 1.0;
570: else data[ii * M + jj] = 0.0;
571: }
572: }
574: /* rotational modes */
575: if (coords) {
576: if (ndm == 2) {
577: data += 2 * M;
578: data[0] = -coords[2 * kk + 1];
579: data[1] = coords[2 * kk];
580: } else {
581: data += 3 * M;
582: data[0] = 0.0;
583: data[M + 0] = coords[3 * kk + 2];
584: data[2 * M + 0] = -coords[3 * kk + 1];
585: data[1] = -coords[3 * kk + 2];
586: data[M + 1] = 0.0;
587: data[2 * M + 1] = coords[3 * kk];
588: data[2] = coords[3 * kk + 1];
589: data[M + 2] = -coords[3 * kk];
590: data[2 * M + 2] = 0.0;
591: }
592: }
593: }
594: }
595: pc_gamg->data_sz = arrsz;
596: PetscFunctionReturn(PETSC_SUCCESS);
597: }
599: /*
600: PCSetData_AGG - called if data is not set with PCSetCoordinates.
601: Looks in Mat for near null space.
602: Does not work for Stokes
604: Input Parameter:
605: . pc -
606: . a_A - matrix to get (near) null space out of.
607: */
608: static PetscErrorCode PCSetData_AGG(PC pc, Mat a_A)
609: {
610: PC_MG *mg = (PC_MG *)pc->data;
611: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
612: MatNullSpace mnull;
614: PetscFunctionBegin;
615: PetscCall(MatGetNearNullSpace(a_A, &mnull));
616: if (!mnull) {
617: DM dm;
619: PetscCall(PCGetDM(pc, &dm));
620: if (!dm) PetscCall(MatGetDM(a_A, &dm));
621: if (dm) {
622: PetscObject deformation;
623: PetscInt Nf;
625: PetscCall(DMGetNumFields(dm, &Nf));
626: if (Nf) {
627: PetscCall(DMGetField(dm, 0, NULL, &deformation));
628: if (deformation) {
629: PetscCall(PetscObjectQuery(deformation, "nearnullspace", (PetscObject *)&mnull));
630: if (!mnull) PetscCall(PetscObjectQuery(deformation, "nullspace", (PetscObject *)&mnull));
631: }
632: }
633: }
634: }
636: if (!mnull) {
637: PetscInt bs, NN, MM;
639: PetscCall(MatGetBlockSize(a_A, &bs));
640: PetscCall(MatGetLocalSize(a_A, &MM, &NN));
641: PetscCheck(MM % bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MM %" PetscInt_FMT " must be divisible by bs %" PetscInt_FMT, MM, bs);
642: PetscCall(PCSetCoordinates_AGG(pc, bs, MM / bs, NULL));
643: } else {
644: PetscReal *nullvec;
645: PetscBool has_const;
646: PetscInt i, j, mlocal, nvec, bs;
647: const Vec *vecs;
648: const PetscScalar *v;
650: PetscCall(MatGetLocalSize(a_A, &mlocal, NULL));
651: PetscCall(MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs));
652: for (i = 0; i < nvec; i++) {
653: PetscCall(VecGetLocalSize(vecs[i], &j));
654: PetscCheck(j == mlocal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Attached null space vector size %" PetscInt_FMT " != matrix size %" PetscInt_FMT, j, mlocal);
655: }
656: pc_gamg->data_sz = (nvec + !!has_const) * mlocal;
657: PetscCall(PetscMalloc1((nvec + !!has_const) * mlocal, &nullvec));
658: if (has_const)
659: for (i = 0; i < mlocal; i++) nullvec[i] = 1.0;
660: for (i = 0; i < nvec; i++) {
661: PetscCall(VecGetArrayRead(vecs[i], &v));
662: for (j = 0; j < mlocal; j++) nullvec[(i + !!has_const) * mlocal + j] = PetscRealPart(v[j]);
663: PetscCall(VecRestoreArrayRead(vecs[i], &v));
664: }
665: pc_gamg->data = nullvec;
666: pc_gamg->data_cell_cols = (nvec + !!has_const);
667: PetscCall(MatGetBlockSize(a_A, &bs));
668: pc_gamg->data_cell_rows = bs;
669: }
670: PetscFunctionReturn(PETSC_SUCCESS);
671: }
673: /*
674: formProl0 - collect null space data for each aggregate, do QR, put R in coarse grid data and Q in P_0
676: Input Parameter:
677: . agg_llists - list of arrays with aggregates -- list from selected vertices of aggregate unselected vertices
678: . bs - row block size
679: . nSAvec - column bs of new P
680: . my0crs - global index of start of locals
681: . data_stride - bs*(nloc nodes + ghost nodes) [data_stride][nSAvec]
682: . data_in[data_stride*nSAvec] - local data on fine grid
683: . flid_fgid[data_stride/bs] - make local to global IDs, includes ghosts in 'locals_llist'
685: Output Parameter:
686: . a_data_out - in with fine grid data (w/ghosts), out with coarse grid data
687: . a_Prol - prolongation operator
688: */
689: 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)
690: {
691: PetscInt Istart, my0, Iend, nloc, clid, flid = 0, aggID, kk, jj, ii, mm, nSelected, minsz, nghosts, out_data_stride;
692: MPI_Comm comm;
693: PetscReal *out_data;
694: PetscCDIntNd *pos;
695: PetscHMapI fgid_flid;
697: PetscFunctionBegin;
698: PetscCall(PetscObjectGetComm((PetscObject)a_Prol, &comm));
699: PetscCall(MatGetOwnershipRange(a_Prol, &Istart, &Iend));
700: nloc = (Iend - Istart) / bs;
701: my0 = Istart / bs;
702: 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);
703: Iend /= bs;
704: nghosts = data_stride / bs - nloc;
706: PetscCall(PetscHMapICreateWithSize(2 * nghosts + 1, &fgid_flid));
708: for (kk = 0; kk < nghosts; kk++) PetscCall(PetscHMapISet(fgid_flid, flid_fgid[nloc + kk], nloc + kk));
710: /* count selected -- same as number of cols of P */
711: for (nSelected = mm = 0; mm < nloc; mm++) {
712: PetscBool ise;
714: PetscCall(PetscCDIsEmptyAt(agg_llists, mm, &ise));
715: if (!ise) nSelected++;
716: }
717: PetscCall(MatGetOwnershipRangeColumn(a_Prol, &ii, &jj));
718: PetscCheck((ii / nSAvec) == my0crs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ii %" PetscInt_FMT " /nSAvec %" PetscInt_FMT " != my0crs %" PetscInt_FMT, ii, nSAvec, my0crs);
719: 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);
721: /* aloc space for coarse point data (output) */
722: out_data_stride = nSelected * nSAvec;
724: PetscCall(PetscMalloc1(out_data_stride * nSAvec, &out_data));
725: for (ii = 0; ii < out_data_stride * nSAvec; ii++) out_data[ii] = PETSC_MAX_REAL;
726: *a_data_out = out_data; /* output - stride nSelected*nSAvec */
728: /* find points and set prolongation */
729: minsz = 100;
730: for (mm = clid = 0; mm < nloc; mm++) {
731: PetscCall(PetscCDCountAt(agg_llists, mm, &jj));
732: if (jj > 0) {
733: const PetscInt lid = mm, cgid = my0crs + clid;
734: PetscInt cids[100]; /* max bs */
735: PetscBLASInt asz, M, N;
736: PetscBLASInt Mdata, LDA, LWORK;
737: PetscScalar *qqc, *qqr, *TAU, *WORK;
738: PetscInt *fids;
739: PetscReal *data;
741: PetscCall(PetscBLASIntCast(jj, &asz));
742: PetscCall(PetscBLASIntCast(asz * bs, &M));
743: PetscCall(PetscBLASIntCast(nSAvec, &N));
744: PetscCall(PetscBLASIntCast(M + ((N - M > 0) ? N - M : 0), &Mdata));
745: PetscCall(PetscBLASIntCast(Mdata, &LDA));
746: PetscCall(PetscBLASIntCast(N * bs, &LWORK));
747: /* count agg */
748: if (asz < minsz) minsz = asz;
750: /* get block */
751: PetscCall(PetscMalloc5(Mdata * N, &qqc, M * N, &qqr, N, &TAU, LWORK, &WORK, M, &fids));
753: aggID = 0;
754: PetscCall(PetscCDGetHeadPos(agg_llists, lid, &pos));
755: while (pos) {
756: PetscInt gid1;
758: PetscCall(PetscCDIntNdGetID(pos, &gid1));
759: PetscCall(PetscCDGetNextPos(agg_llists, lid, &pos));
761: if (gid1 >= my0 && gid1 < Iend) flid = gid1 - my0;
762: else {
763: PetscCall(PetscHMapIGet(fgid_flid, gid1, &flid));
764: PetscCheck(flid >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cannot find gid1 in table");
765: }
766: /* copy in B_i matrix - column-oriented */
767: data = &data_in[flid * bs];
768: for (ii = 0; ii < bs; ii++) {
769: for (jj = 0; jj < N; jj++) {
770: PetscReal d = data[jj * data_stride + ii];
772: qqc[jj * Mdata + aggID * bs + ii] = d;
773: }
774: }
775: /* set fine IDs */
776: for (kk = 0; kk < bs; kk++) fids[aggID * bs + kk] = flid_fgid[flid] * bs + kk;
777: aggID++;
778: }
780: /* pad with zeros */
781: for (ii = asz * bs; ii < Mdata; ii++) {
782: for (jj = 0; jj < N; jj++, kk++) qqc[jj * Mdata + ii] = .0;
783: }
785: /* QR */
786: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
787: PetscCallLAPACKInfo("LAPACKgeqrf", LAPACKgeqrf_(&Mdata, &N, qqc, &LDA, TAU, WORK, &LWORK, &info));
788: PetscCall(PetscFPTrapPop());
789: /* get R - column-oriented - output B_{i+1} */
790: {
791: PetscReal *data = &out_data[clid * nSAvec];
793: for (jj = 0; jj < nSAvec; jj++) {
794: for (ii = 0; ii < nSAvec; ii++) {
795: 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);
796: if (ii <= jj) data[jj * out_data_stride + ii] = PetscRealPart(qqc[jj * Mdata + ii]);
797: else data[jj * out_data_stride + ii] = 0.;
798: }
799: }
800: }
802: /* get Q - row-oriented */
803: PetscCallLAPACKInfo("LAPACKorgqr", LAPACKorgqr_(&Mdata, &N, &N, qqc, &LDA, TAU, WORK, &LWORK, &info));
805: for (ii = 0; ii < M; ii++) {
806: for (jj = 0; jj < N; jj++) qqr[N * ii + jj] = qqc[jj * Mdata + ii];
807: }
809: /* add diagonal block of P0 */
810: for (kk = 0; kk < N; kk++) cids[kk] = N * cgid + kk; /* global col IDs in P0 */
811: PetscCall(MatSetValues(a_Prol, M, fids, N, cids, qqr, INSERT_VALUES));
812: PetscCall(PetscFree5(qqc, qqr, TAU, WORK, fids));
813: clid++;
814: } /* coarse agg */
815: } /* for all fine nodes */
816: PetscCall(MatAssemblyBegin(a_Prol, MAT_FINAL_ASSEMBLY));
817: PetscCall(MatAssemblyEnd(a_Prol, MAT_FINAL_ASSEMBLY));
818: PetscCall(PetscHMapIDestroy(&fgid_flid));
819: PetscFunctionReturn(PETSC_SUCCESS);
820: }
822: static PetscErrorCode PCView_GAMG_AGG(PC pc, PetscViewer viewer)
823: {
824: PC_MG *mg = (PC_MG *)pc->data;
825: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
826: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
828: PetscFunctionBegin;
829: PetscCall(PetscViewerASCIIPrintf(viewer, " AGG specific options\n"));
830: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of levels of aggressive coarsening %" PetscInt_FMT "\n", pc_gamg_agg->aggressive_coarsening_levels));
831: if (pc_gamg_agg->aggressive_coarsening_levels > 0) {
832: PetscCall(PetscViewerASCIIPrintf(viewer, " %s aggressive coarsening\n", !pc_gamg_agg->use_aggressive_square_graph ? "MIS-k" : "Square graph"));
833: 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));
834: }
835: PetscCall(PetscViewerASCIIPushTab(viewer));
836: PetscCall(PetscViewerASCIIPushTab(viewer));
837: PetscCall(PetscViewerASCIIPushTab(viewer));
838: PetscCall(PetscViewerASCIIPushTab(viewer));
839: if (pc_gamg_agg->crs) PetscCall(MatCoarsenView(pc_gamg_agg->crs, viewer));
840: else PetscCall(PetscViewerASCIIPrintf(viewer, "Coarsening algorithm not yet selected\n"));
841: PetscCall(PetscViewerASCIIPopTab(viewer));
842: PetscCall(PetscViewerASCIIPopTab(viewer));
843: PetscCall(PetscViewerASCIIPopTab(viewer));
844: PetscCall(PetscViewerASCIIPopTab(viewer));
845: PetscCall(PetscViewerASCIIPrintf(viewer, " Number smoothing steps to construct prolongation %" PetscInt_FMT "\n", pc_gamg_agg->nsmooths));
846: PetscFunctionReturn(PETSC_SUCCESS);
847: }
849: static PetscErrorCode PCGAMGCreateGraph_AGG(PC pc, Mat Amat, Mat *a_Gmat)
850: {
851: PC_MG *mg = (PC_MG *)pc->data;
852: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
853: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
854: const PetscReal vfilter = pc_gamg->threshold[pc_gamg->current_level];
855: PetscBool ishem, ismis;
856: const char *prefix;
857: MatInfo info0, info1;
858: PetscInt bs;
860: PetscFunctionBegin;
861: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
862: /* 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 */
863: /* MATCOARSENHEM requires numerical weights for edges so ensure they are computed */
864: PetscCall(MatCoarsenDestroy(&pc_gamg_agg->crs));
865: PetscCall(MatCoarsenCreate(PetscObjectComm((PetscObject)pc), &pc_gamg_agg->crs));
866: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix));
867: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc_gamg_agg->crs, prefix));
868: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)pc_gamg_agg->crs, "pc_gamg_"));
869: PetscCall(MatCoarsenSetFromOptions(pc_gamg_agg->crs));
870: PetscCall(MatGetBlockSize(Amat, &bs));
871: // check for valid indices wrt bs
872: for (int ii = 0; ii < pc_gamg_agg->crs->strength_index_size; ii++) {
873: 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",
874: pc_gamg_agg->crs->strength_index[ii], bs);
875: }
876: PetscCall(PetscObjectTypeCompare((PetscObject)pc_gamg_agg->crs, MATCOARSENHEM, &ishem));
877: if (ishem) {
878: 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));
879: pc_gamg_agg->aggressive_coarsening_levels = 0; // aggressive and HEM does not make sense
880: PetscCall(MatCoarsenSetMaximumIterations(pc_gamg_agg->crs, pc_gamg_agg->crs->max_it)); // for code coverage
881: PetscCall(MatCoarsenSetThreshold(pc_gamg_agg->crs, vfilter)); // for code coverage
882: } else {
883: PetscCall(PetscObjectTypeCompare((PetscObject)pc_gamg_agg->crs, MATCOARSENMIS, &ismis));
884: if (ismis && pc_gamg_agg->aggressive_coarsening_levels && !pc_gamg_agg->use_aggressive_square_graph) {
885: PetscCall(PetscInfo(pc, "MIS and aggressive coarsening and no square graph: force square graph\n"));
886: pc_gamg_agg->use_aggressive_square_graph = PETSC_TRUE;
887: }
888: }
889: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
890: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_GRAPH], 0, 0, 0, 0));
891: PetscCall(MatGetInfo(Amat, MAT_LOCAL, &info0)); /* global reduction */
893: if (ishem || pc_gamg_agg->use_low_mem_filter) {
894: 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));
895: } else {
896: // make scalar graph, symmetrize if not known to be symmetric, scale, but do not filter (expensive)
897: 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));
898: if (vfilter >= 0) {
899: PetscInt Istart, Iend, ncols, nnz0, nnz1, NN, MM, nloc;
900: Mat tGmat, Gmat = *a_Gmat;
901: MPI_Comm comm;
902: const PetscScalar *vals;
903: const PetscInt *idx;
904: PetscInt *d_nnz, *o_nnz, kk, *garray = NULL, *AJ, maxcols = 0;
905: MatScalar *AA; // this is checked in graph
906: PetscBool isseqaij;
907: Mat a, b, c;
908: MatType jtype;
910: PetscCall(PetscObjectGetComm((PetscObject)Gmat, &comm));
911: PetscCall(PetscObjectBaseTypeCompare((PetscObject)Gmat, MATSEQAIJ, &isseqaij));
912: PetscCall(MatGetType(Gmat, &jtype));
913: PetscCall(MatCreate(comm, &tGmat));
914: PetscCall(MatSetType(tGmat, jtype));
916: /* TODO GPU: this can be called when filter = 0 -> Probably provide MatAIJThresholdCompress that compresses the entries below a threshold?
917: Also, if the matrix is symmetric, can we skip this
918: operation? It can be very expensive on large matrices. */
920: // global sizes
921: PetscCall(MatGetSize(Gmat, &MM, &NN));
922: PetscCall(MatGetOwnershipRange(Gmat, &Istart, &Iend));
923: nloc = Iend - Istart;
924: PetscCall(PetscMalloc2(nloc, &d_nnz, nloc, &o_nnz));
925: if (isseqaij) {
926: a = Gmat;
927: b = NULL;
928: } else {
929: Mat_MPIAIJ *d = (Mat_MPIAIJ *)Gmat->data;
931: a = d->A;
932: b = d->B;
933: garray = d->garray;
934: }
935: /* Determine upper bound on non-zeros needed in new filtered matrix */
936: for (PetscInt row = 0; row < nloc; row++) {
937: PetscCall(MatGetRow(a, row, &ncols, NULL, NULL));
938: d_nnz[row] = ncols;
939: if (ncols > maxcols) maxcols = ncols;
940: PetscCall(MatRestoreRow(a, row, &ncols, NULL, NULL));
941: }
942: if (b) {
943: for (PetscInt row = 0; row < nloc; row++) {
944: PetscCall(MatGetRow(b, row, &ncols, NULL, NULL));
945: o_nnz[row] = ncols;
946: if (ncols > maxcols) maxcols = ncols;
947: PetscCall(MatRestoreRow(b, row, &ncols, NULL, NULL));
948: }
949: }
950: PetscCall(MatSetSizes(tGmat, nloc, nloc, MM, MM));
951: PetscCall(MatSetBlockSizes(tGmat, 1, 1));
952: PetscCall(MatSeqAIJSetPreallocation(tGmat, 0, d_nnz));
953: PetscCall(MatMPIAIJSetPreallocation(tGmat, 0, d_nnz, 0, o_nnz));
954: PetscCall(MatSetOption(tGmat, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
955: PetscCall(PetscFree2(d_nnz, o_nnz));
956: PetscCall(PetscMalloc2(maxcols, &AA, maxcols, &AJ));
957: nnz0 = nnz1 = 0;
958: for (c = a, kk = 0; c && kk < 2; c = b, kk++) {
959: for (PetscInt row = 0, grow = Istart, ncol_row, jj; row < nloc; row++, grow++) {
960: PetscCall(MatGetRow(c, row, &ncols, &idx, &vals));
961: for (ncol_row = jj = 0; jj < ncols; jj++, nnz0++) {
962: PetscScalar sv = PetscAbs(PetscRealPart(vals[jj]));
963: if (PetscRealPart(sv) > vfilter) {
964: PetscInt cid = idx[jj] + Istart; //diag
966: nnz1++;
967: if (c != a) cid = garray[idx[jj]];
968: AA[ncol_row] = vals[jj];
969: AJ[ncol_row] = cid;
970: ncol_row++;
971: }
972: }
973: PetscCall(MatRestoreRow(c, row, &ncols, &idx, &vals));
974: PetscCall(MatSetValues(tGmat, 1, &grow, ncol_row, AJ, AA, INSERT_VALUES));
975: }
976: }
977: PetscCall(PetscFree2(AA, AJ));
978: PetscCall(MatAssemblyBegin(tGmat, MAT_FINAL_ASSEMBLY));
979: PetscCall(MatAssemblyEnd(tGmat, MAT_FINAL_ASSEMBLY));
980: PetscCall(MatPropagateSymmetryOptions(Gmat, tGmat)); /* Normal Mat options are not relevant ? */
981: 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));
982: PetscCall(MatViewFromOptions(tGmat, NULL, "-mat_filter_graph_view"));
983: PetscCall(MatDestroy(&Gmat));
984: *a_Gmat = tGmat;
985: }
986: }
988: PetscCall(MatGetInfo(*a_Gmat, MAT_LOCAL, &info1)); /* global reduction */
989: 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));
990: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_GRAPH], 0, 0, 0, 0));
991: PetscFunctionReturn(PETSC_SUCCESS);
992: }
994: typedef PetscInt NState;
995: static const NState NOT_DONE = -2;
996: static const NState DELETED = -1;
997: static const NState REMOVED = -3;
998: #define IS_SELECTED(s) ((s) != DELETED && (s) != NOT_DONE && (s) != REMOVED)
1000: /*
1001: fixAggregatesWithSquare - greedy grab of with G1 (unsquared graph) -- AIJ specific -- change to fixAggregatesWithSquare -- TODD
1002: - AGG-MG specific: clears singletons out of 'selected_2'
1004: Input Parameter:
1005: . Gmat_2 - global matrix of squared graph (data not defined)
1006: . Gmat_1 - base graph to grab with base graph
1007: Input/Output Parameter:
1008: . aggs_2 - linked list of aggs with gids)
1009: */
1010: static PetscErrorCode fixAggregatesWithSquare(PC pc, Mat Gmat_2, Mat Gmat_1, PetscCoarsenData *aggs_2)
1011: {
1012: PetscBool isMPI;
1013: Mat_SeqAIJ *matA_1, *matB_1 = NULL;
1014: MPI_Comm comm;
1015: PetscInt lid, *ii, *idx, ix, Iend, my0, kk, n, j;
1016: Mat_MPIAIJ *mpimat_2 = NULL, *mpimat_1 = NULL;
1017: const PetscInt nloc = Gmat_2->rmap->n;
1018: PetscScalar *cpcol_1_state, *cpcol_2_state, *cpcol_2_par_orig, *lid_parent_gid;
1019: PetscInt *lid_cprowID_1 = NULL;
1020: NState *lid_state;
1021: Vec ghost_par_orig2;
1022: PetscMPIInt rank;
1024: PetscFunctionBegin;
1025: PetscCall(PetscObjectGetComm((PetscObject)Gmat_2, &comm));
1026: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1027: PetscCall(MatGetOwnershipRange(Gmat_1, &my0, &Iend));
1029: /* get submatrices */
1030: PetscCall(PetscStrbeginswith(((PetscObject)Gmat_1)->type_name, MATMPIAIJ, &isMPI));
1031: PetscCall(PetscInfo(pc, "isMPI = %s\n", isMPI ? "yes" : "no"));
1032: PetscCall(PetscMalloc3(nloc, &lid_state, nloc, &lid_parent_gid, nloc, &lid_cprowID_1));
1033: for (lid = 0; lid < nloc; lid++) lid_cprowID_1[lid] = -1;
1034: if (isMPI) {
1035: /* grab matrix objects */
1036: mpimat_2 = (Mat_MPIAIJ *)Gmat_2->data;
1037: mpimat_1 = (Mat_MPIAIJ *)Gmat_1->data;
1038: matA_1 = (Mat_SeqAIJ *)mpimat_1->A->data;
1039: matB_1 = (Mat_SeqAIJ *)mpimat_1->B->data;
1041: /* force compressed row storage for B matrix in AuxMat */
1042: PetscCall(MatCheckCompressedRow(mpimat_1->B, matB_1->nonzerorowcnt, &matB_1->compressedrow, matB_1->i, Gmat_1->rmap->n, -1.0));
1043: for (ix = 0; ix < matB_1->compressedrow.nrows; ix++) {
1044: PetscInt lid = matB_1->compressedrow.rindex[ix];
1046: PetscCheck(lid <= nloc && lid >= -1, PETSC_COMM_SELF, PETSC_ERR_USER, "lid %" PetscInt_FMT " out of range. nloc = %" PetscInt_FMT, lid, nloc);
1047: if (lid != -1) lid_cprowID_1[lid] = ix;
1048: }
1049: } else {
1050: PetscBool isAIJ;
1052: PetscCall(PetscStrbeginswith(((PetscObject)Gmat_1)->type_name, MATSEQAIJ, &isAIJ));
1053: PetscCheck(isAIJ, PETSC_COMM_SELF, PETSC_ERR_USER, "Require AIJ matrix.");
1054: matA_1 = (Mat_SeqAIJ *)Gmat_1->data;
1055: }
1056: if (nloc > 0) PetscCheck(!matB_1 || matB_1->compressedrow.use, PETSC_COMM_SELF, PETSC_ERR_PLIB, "matB_1 && !matB_1->compressedrow.use: PETSc bug???");
1057: /* get state of locals and selected gid for deleted */
1058: for (lid = 0; lid < nloc; lid++) {
1059: lid_parent_gid[lid] = -1.0;
1060: lid_state[lid] = DELETED;
1061: }
1063: /* set lid_state */
1064: for (lid = 0; lid < nloc; lid++) {
1065: PetscCDIntNd *pos;
1067: PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
1068: if (pos) {
1069: PetscInt gid1;
1071: PetscCall(PetscCDIntNdGetID(pos, &gid1));
1072: PetscCheck(gid1 == lid + my0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "gid1 %" PetscInt_FMT " != lid %" PetscInt_FMT " + my0 %" PetscInt_FMT, gid1, lid, my0);
1073: lid_state[lid] = gid1;
1074: }
1075: }
1077: /* map local to selected local, DELETED means a ghost owns it */
1078: for (lid = 0; lid < nloc; lid++) {
1079: NState state = lid_state[lid];
1081: if (IS_SELECTED(state)) {
1082: PetscCDIntNd *pos;
1084: PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
1085: while (pos) {
1086: PetscInt gid1;
1088: PetscCall(PetscCDIntNdGetID(pos, &gid1));
1089: PetscCall(PetscCDGetNextPos(aggs_2, lid, &pos));
1090: if (gid1 >= my0 && gid1 < Iend) lid_parent_gid[gid1 - my0] = (PetscScalar)(lid + my0);
1091: }
1092: }
1093: }
1094: /* get 'cpcol_1/2_state' & cpcol_2_par_orig - uses mpimat_1/2->lvec for temp space */
1095: if (isMPI) {
1096: Vec tempVec;
1098: /* get 'cpcol_1_state' */
1099: PetscCall(MatCreateVecs(Gmat_1, &tempVec, NULL));
1100: for (kk = 0, j = my0; kk < nloc; kk++, j++) {
1101: PetscScalar v = (PetscScalar)lid_state[kk];
1103: PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
1104: }
1105: PetscCall(VecAssemblyBegin(tempVec));
1106: PetscCall(VecAssemblyEnd(tempVec));
1107: PetscCall(VecScatterBegin(mpimat_1->Mvctx, tempVec, mpimat_1->lvec, INSERT_VALUES, SCATTER_FORWARD));
1108: PetscCall(VecScatterEnd(mpimat_1->Mvctx, tempVec, mpimat_1->lvec, INSERT_VALUES, SCATTER_FORWARD));
1109: PetscCall(VecGetArray(mpimat_1->lvec, &cpcol_1_state));
1110: /* get 'cpcol_2_state' */
1111: PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, mpimat_2->lvec, INSERT_VALUES, SCATTER_FORWARD));
1112: PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, mpimat_2->lvec, INSERT_VALUES, SCATTER_FORWARD));
1113: PetscCall(VecGetArray(mpimat_2->lvec, &cpcol_2_state));
1114: /* get 'cpcol_2_par_orig' */
1115: for (kk = 0, j = my0; kk < nloc; kk++, j++) {
1116: PetscScalar v = lid_parent_gid[kk];
1118: PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
1119: }
1120: PetscCall(VecAssemblyBegin(tempVec));
1121: PetscCall(VecAssemblyEnd(tempVec));
1122: PetscCall(VecDuplicate(mpimat_2->lvec, &ghost_par_orig2));
1123: PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghost_par_orig2, INSERT_VALUES, SCATTER_FORWARD));
1124: PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghost_par_orig2, INSERT_VALUES, SCATTER_FORWARD));
1125: PetscCall(VecGetArray(ghost_par_orig2, &cpcol_2_par_orig));
1127: PetscCall(VecDestroy(&tempVec));
1128: } /* ismpi */
1129: for (lid = 0; lid < nloc; lid++) {
1130: NState state = lid_state[lid];
1132: if (IS_SELECTED(state)) {
1133: /* steal locals */
1134: ii = matA_1->i;
1135: n = ii[lid + 1] - ii[lid];
1136: idx = matA_1->j + ii[lid];
1137: for (j = 0; j < n; j++) {
1138: PetscInt lidj = idx[j], sgid;
1139: NState statej = lid_state[lidj];
1141: if (statej == DELETED && (sgid = (PetscInt)PetscRealPart(lid_parent_gid[lidj])) != lid + my0) { /* steal local */
1142: lid_parent_gid[lidj] = (PetscScalar)(lid + my0); /* send this if sgid is not local */
1143: if (sgid >= my0 && sgid < Iend) { /* I'm stealing this local from a local sgid */
1144: PetscInt hav = 0, slid = sgid - my0, gidj = lidj + my0;
1145: PetscCDIntNd *pos, *last = NULL;
1147: /* looking for local from local so id_llist_2 works */
1148: PetscCall(PetscCDGetHeadPos(aggs_2, slid, &pos));
1149: while (pos) {
1150: PetscInt gid;
1152: PetscCall(PetscCDIntNdGetID(pos, &gid));
1153: if (gid == gidj) {
1154: PetscCheck(last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "last cannot be null");
1155: PetscCall(PetscCDRemoveNextNode(aggs_2, slid, last));
1156: PetscCall(PetscCDAppendNode(aggs_2, lid, pos));
1157: hav = 1;
1158: break;
1159: } else last = pos;
1160: PetscCall(PetscCDGetNextPos(aggs_2, slid, &pos));
1161: }
1162: if (hav != 1) {
1163: PetscCheck(hav, PETSC_COMM_SELF, PETSC_ERR_PLIB, "failed to find adj in 'selected' lists - structurally unsymmetric matrix");
1164: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "found node %" PetscInt_FMT " times???", hav);
1165: }
1166: } else { /* I'm stealing this local, owned by a ghost */
1167: 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.",
1168: ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "");
1169: PetscCall(PetscCDAppendID(aggs_2, lid, lidj + my0));
1170: }
1171: }
1172: } /* local neighbors */
1173: } else if (state == DELETED /* && lid_cprowID_1 */) {
1174: PetscInt sgidold = (PetscInt)PetscRealPart(lid_parent_gid[lid]);
1176: /* see if I have a selected ghost neighbor that will steal me */
1177: if ((ix = lid_cprowID_1[lid]) != -1) {
1178: ii = matB_1->compressedrow.i;
1179: n = ii[ix + 1] - ii[ix];
1180: idx = matB_1->j + ii[ix];
1181: for (j = 0; j < n; j++) {
1182: PetscInt cpid = idx[j];
1183: NState statej = (NState)PetscRealPart(cpcol_1_state[cpid]);
1185: if (IS_SELECTED(statej) && sgidold != statej) { /* ghost will steal this, remove from my list */
1186: lid_parent_gid[lid] = (PetscScalar)statej; /* send who selected */
1187: if (sgidold >= my0 && sgidold < Iend) { /* this was mine */
1188: PetscInt hav = 0, oldslidj = sgidold - my0;
1189: PetscCDIntNd *pos, *last = NULL;
1191: /* remove from 'oldslidj' list */
1192: PetscCall(PetscCDGetHeadPos(aggs_2, oldslidj, &pos));
1193: while (pos) {
1194: PetscInt gid;
1196: PetscCall(PetscCDIntNdGetID(pos, &gid));
1197: if (lid + my0 == gid) {
1198: /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove lid from oldslidj list *\/ */
1199: PetscCheck(last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "last cannot be null");
1200: PetscCall(PetscCDRemoveNextNode(aggs_2, oldslidj, last));
1201: /* ghost (PetscScalar)statej will add this later */
1202: hav = 1;
1203: break;
1204: } else last = pos;
1205: PetscCall(PetscCDGetNextPos(aggs_2, oldslidj, &pos));
1206: }
1207: if (hav != 1) {
1208: PetscCheck(hav, PETSC_COMM_SELF, PETSC_ERR_PLIB, "failed to find (hav=%" PetscInt_FMT ") adj in 'selected' lists - structurally unsymmetric matrix", hav);
1209: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "found node %" PetscInt_FMT " times???", hav);
1210: }
1211: } else {
1212: /* TODO: ghosts remove this later */
1213: }
1214: }
1215: }
1216: }
1217: } /* selected/deleted */
1218: } /* node loop */
1220: if (isMPI) {
1221: PetscScalar *cpcol_2_parent, *cpcol_2_gid;
1222: Vec tempVec, ghostgids2, ghostparents2;
1223: PetscInt cpid, nghost_2;
1224: PetscHMapI gid_cpid;
1226: PetscCall(VecGetSize(mpimat_2->lvec, &nghost_2));
1227: PetscCall(MatCreateVecs(Gmat_2, &tempVec, NULL));
1229: /* get 'cpcol_2_parent' */
1230: for (kk = 0, j = my0; kk < nloc; kk++, j++) PetscCall(VecSetValues(tempVec, 1, &j, &lid_parent_gid[kk], INSERT_VALUES));
1231: PetscCall(VecAssemblyBegin(tempVec));
1232: PetscCall(VecAssemblyEnd(tempVec));
1233: PetscCall(VecDuplicate(mpimat_2->lvec, &ghostparents2));
1234: PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghostparents2, INSERT_VALUES, SCATTER_FORWARD));
1235: PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghostparents2, INSERT_VALUES, SCATTER_FORWARD));
1236: PetscCall(VecGetArray(ghostparents2, &cpcol_2_parent));
1238: /* get 'cpcol_2_gid' */
1239: for (kk = 0, j = my0; kk < nloc; kk++, j++) {
1240: PetscScalar v = (PetscScalar)j;
1242: PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
1243: }
1244: PetscCall(VecAssemblyBegin(tempVec));
1245: PetscCall(VecAssemblyEnd(tempVec));
1246: PetscCall(VecDuplicate(mpimat_2->lvec, &ghostgids2));
1247: PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghostgids2, INSERT_VALUES, SCATTER_FORWARD));
1248: PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghostgids2, INSERT_VALUES, SCATTER_FORWARD));
1249: PetscCall(VecGetArray(ghostgids2, &cpcol_2_gid));
1250: PetscCall(VecDestroy(&tempVec));
1252: /* look for deleted ghosts and add to table */
1253: PetscCall(PetscHMapICreateWithSize(2 * nghost_2 + 1, &gid_cpid));
1254: for (cpid = 0; cpid < nghost_2; cpid++) {
1255: NState state = (NState)PetscRealPart(cpcol_2_state[cpid]);
1257: if (state == DELETED) {
1258: PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
1259: PetscInt sgid_old = (PetscInt)PetscRealPart(cpcol_2_par_orig[cpid]);
1261: if (sgid_old == -1 && sgid_new != -1) {
1262: PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
1264: PetscCall(PetscHMapISet(gid_cpid, gid, cpid));
1265: }
1266: }
1267: }
1269: /* look for deleted ghosts and see if they moved - remove it */
1270: for (lid = 0; lid < nloc; lid++) {
1271: NState state = lid_state[lid];
1273: if (IS_SELECTED(state)) {
1274: PetscCDIntNd *pos, *last = NULL;
1276: /* look for deleted ghosts and see if they moved */
1277: PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
1278: while (pos) {
1279: PetscInt gid;
1281: PetscCall(PetscCDIntNdGetID(pos, &gid));
1282: if (gid < my0 || gid >= Iend) {
1283: PetscCall(PetscHMapIGet(gid_cpid, gid, &cpid));
1284: if (cpid != -1) {
1285: /* a moved ghost - */
1286: /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove 'flid' from list *\/ */
1287: PetscCall(PetscCDRemoveNextNode(aggs_2, lid, last));
1288: } else last = pos;
1289: } else last = pos;
1291: PetscCall(PetscCDGetNextPos(aggs_2, lid, &pos));
1292: } /* loop over list of deleted */
1293: } /* selected */
1294: }
1295: PetscCall(PetscHMapIDestroy(&gid_cpid));
1297: /* look at ghosts, see if they changed - and it */
1298: for (cpid = 0; cpid < nghost_2; cpid++) {
1299: PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
1301: if (sgid_new >= my0 && sgid_new < Iend) { /* this is mine */
1302: PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
1303: PetscInt slid_new = sgid_new - my0, hav = 0;
1304: PetscCDIntNd *pos;
1306: /* search for this gid to see if I have it */
1307: PetscCall(PetscCDGetHeadPos(aggs_2, slid_new, &pos));
1308: while (pos) {
1309: PetscInt gidj;
1311: PetscCall(PetscCDIntNdGetID(pos, &gidj));
1312: PetscCall(PetscCDGetNextPos(aggs_2, slid_new, &pos));
1314: if (gidj == gid) {
1315: hav = 1;
1316: break;
1317: }
1318: }
1319: if (hav != 1) {
1320: /* insert 'flidj' into head of llist */
1321: PetscCall(PetscCDAppendID(aggs_2, slid_new, gid));
1322: }
1323: }
1324: }
1325: PetscCall(VecRestoreArray(mpimat_1->lvec, &cpcol_1_state));
1326: PetscCall(VecRestoreArray(mpimat_2->lvec, &cpcol_2_state));
1327: PetscCall(VecRestoreArray(ghostparents2, &cpcol_2_parent));
1328: PetscCall(VecRestoreArray(ghostgids2, &cpcol_2_gid));
1329: PetscCall(VecDestroy(&ghostgids2));
1330: PetscCall(VecDestroy(&ghostparents2));
1331: PetscCall(VecDestroy(&ghost_par_orig2));
1332: }
1333: PetscCall(PetscFree3(lid_state, lid_parent_gid, lid_cprowID_1));
1334: PetscFunctionReturn(PETSC_SUCCESS);
1335: }
1337: /*
1338: PCGAMGCoarsen_AGG - supports squaring the graph (deprecated) and new graph for
1339: communication of QR data used with HEM and MISk coarsening
1341: Input Parameter:
1342: . a_pc - this
1344: Input/Output Parameter:
1345: . a_Gmat1 - graph to coarsen (in), graph off processor edges for QR gather scatter (out)
1347: Output Parameter:
1348: . agg_lists - list of aggregates
1350: */
1351: static PetscErrorCode PCGAMGCoarsen_AGG(PC a_pc, Mat *a_Gmat1, PetscCoarsenData **agg_lists)
1352: {
1353: PC_MG *mg = (PC_MG *)a_pc->data;
1354: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1355: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
1356: Mat Gmat2, Gmat1 = *a_Gmat1; /* aggressive graph */
1357: IS perm;
1358: PetscInt Istart, Iend, Ii, nloc, bs, nn;
1359: PetscInt *permute, *degree;
1360: PetscBool *bIndexSet;
1361: PetscReal hashfact;
1362: PetscInt iSwapIndex;
1363: PetscRandom random;
1364: MPI_Comm comm;
1366: PetscFunctionBegin;
1367: PetscCall(PetscObjectGetComm((PetscObject)Gmat1, &comm));
1368: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
1369: PetscCall(MatGetLocalSize(Gmat1, &nn, NULL));
1370: PetscCall(MatGetBlockSize(Gmat1, &bs));
1371: PetscCheck(bs == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "bs %" PetscInt_FMT " must be 1", bs);
1372: nloc = nn / bs;
1373: /* get MIS aggs - randomize */
1374: PetscCall(PetscMalloc2(nloc, &permute, nloc, °ree));
1375: PetscCall(PetscCalloc1(nloc, &bIndexSet));
1376: for (Ii = 0; Ii < nloc; Ii++) permute[Ii] = Ii;
1377: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &random));
1378: PetscCall(MatGetOwnershipRange(Gmat1, &Istart, &Iend));
1379: for (Ii = 0; Ii < nloc; Ii++) {
1380: PetscInt nc;
1382: PetscCall(MatGetRow(Gmat1, Istart + Ii, &nc, NULL, NULL));
1383: degree[Ii] = nc;
1384: PetscCall(MatRestoreRow(Gmat1, Istart + Ii, &nc, NULL, NULL));
1385: }
1386: for (Ii = 0; Ii < nloc; Ii++) {
1387: PetscCall(PetscRandomGetValueReal(random, &hashfact));
1388: iSwapIndex = (PetscInt)(hashfact * nloc) % nloc;
1389: if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) {
1390: PetscInt iTemp = permute[iSwapIndex];
1392: permute[iSwapIndex] = permute[Ii];
1393: permute[Ii] = iTemp;
1394: iTemp = degree[iSwapIndex];
1395: degree[iSwapIndex] = degree[Ii];
1396: degree[Ii] = iTemp;
1397: bIndexSet[iSwapIndex] = PETSC_TRUE;
1398: }
1399: }
1400: // apply minimum degree ordering -- NEW
1401: if (pc_gamg_agg->use_minimum_degree_ordering) PetscCall(PetscSortIntWithArray(nloc, degree, permute));
1402: PetscCall(PetscFree(bIndexSet));
1403: PetscCall(PetscRandomDestroy(&random));
1404: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_USE_POINTER, &perm));
1405: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0));
1406: // square graph
1407: 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));
1408: else Gmat2 = Gmat1;
1409: // switch to old MIS-1 for square graph
1410: if (pc_gamg->current_level < pc_gamg_agg->aggressive_coarsening_levels) {
1411: if (!pc_gamg_agg->use_aggressive_square_graph) PetscCall(MatCoarsenMISKSetDistance(pc_gamg_agg->crs, pc_gamg_agg->aggressive_mis_k)); // hardwire to MIS-2
1412: else PetscCall(MatCoarsenSetType(pc_gamg_agg->crs, MATCOARSENMIS)); // old MIS -- side effect
1413: } else if (pc_gamg_agg->use_aggressive_square_graph && pc_gamg_agg->aggressive_coarsening_levels > 0) { // we reset the MIS
1414: const char *prefix;
1416: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)a_pc, &prefix));
1417: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc_gamg_agg->crs, prefix));
1418: PetscCall(MatCoarsenSetFromOptions(pc_gamg_agg->crs)); // get the default back on non-aggressive levels when square graph switched to old MIS
1419: }
1420: PetscCall(MatCoarsenSetAdjacency(pc_gamg_agg->crs, Gmat2));
1421: PetscCall(MatCoarsenSetStrictAggs(pc_gamg_agg->crs, PETSC_TRUE));
1422: PetscCall(MatCoarsenSetGreedyOrdering(pc_gamg_agg->crs, perm));
1423: PetscCall(MatCoarsenApply(pc_gamg_agg->crs));
1424: PetscCall(MatCoarsenGetData(pc_gamg_agg->crs, agg_lists)); /* output */
1426: PetscCall(ISDestroy(&perm));
1427: PetscCall(PetscFree2(permute, degree));
1428: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0));
1430: if (Gmat2 != Gmat1) { // square graph, we need ghosts for selected
1431: PetscCoarsenData *llist = *agg_lists;
1433: PetscCall(fixAggregatesWithSquare(a_pc, Gmat2, Gmat1, *agg_lists));
1434: PetscCall(MatDestroy(&Gmat1));
1435: *a_Gmat1 = Gmat2; /* output */
1436: PetscCall(PetscCDSetMat(llist, *a_Gmat1)); /* Need a graph with ghosts here */
1437: }
1438: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
1439: PetscFunctionReturn(PETSC_SUCCESS);
1440: }
1442: /*
1443: PCGAMGConstructProlongator_AGG
1445: Input Parameter:
1446: . pc - this
1447: . Amat - matrix on this fine level
1448: . Graph - used to get ghost data for nodes in
1449: . agg_lists - list of aggregates
1450: Output Parameter:
1451: . a_P_out - prolongation operator to the next level
1452: */
1453: static PetscErrorCode PCGAMGConstructProlongator_AGG(PC pc, Mat Amat, PetscCoarsenData *agg_lists, Mat *a_P_out)
1454: {
1455: PC_MG *mg = (PC_MG *)pc->data;
1456: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1457: const PetscInt col_bs = pc_gamg->data_cell_cols;
1458: PetscInt Istart, Iend, nloc, ii, jj, kk, my0, nLocalSelected, bs;
1459: Mat Gmat, Prol;
1460: PetscMPIInt size;
1461: MPI_Comm comm;
1462: PetscReal *data_w_ghost;
1463: PetscInt myCrs0, nbnodes = 0, *flid_fgid;
1464: MatType mtype;
1466: PetscFunctionBegin;
1467: PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
1468: PetscCheck(col_bs >= 1, comm, PETSC_ERR_PLIB, "Column bs cannot be less than 1");
1469: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1470: PetscCallMPI(MPI_Comm_size(comm, &size));
1471: PetscCall(MatGetOwnershipRange(Amat, &Istart, &Iend));
1472: PetscCall(MatGetBlockSize(Amat, &bs));
1473: nloc = (Iend - Istart) / bs;
1474: my0 = Istart / bs;
1475: 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);
1476: PetscCall(PetscCDGetMat(agg_lists, &Gmat)); // get auxiliary matrix for ghost edges for size > 1
1478: /* get 'nLocalSelected' */
1479: for (ii = 0, nLocalSelected = 0; ii < nloc; ii++) {
1480: PetscBool ise;
1482: /* filter out singletons 0 or 1? */
1483: PetscCall(PetscCDIsEmptyAt(agg_lists, ii, &ise));
1484: if (!ise) nLocalSelected++;
1485: }
1487: /* create prolongator, create P matrix */
1488: PetscCall(MatGetType(Amat, &mtype));
1489: PetscCall(MatCreate(comm, &Prol));
1490: PetscCall(MatSetSizes(Prol, nloc * bs, nLocalSelected * col_bs, PETSC_DETERMINE, PETSC_DETERMINE));
1491: PetscCall(MatSetBlockSizes(Prol, bs, col_bs)); // should this be before MatSetSizes?
1492: PetscCall(MatSetType(Prol, mtype));
1493: #if PetscDefined(HAVE_DEVICE)
1494: PetscBool flg;
1495: PetscCall(MatBoundToCPU(Amat, &flg));
1496: PetscCall(MatBindToCPU(Prol, flg));
1497: if (flg) PetscCall(MatSetBindingPropagates(Prol, PETSC_TRUE));
1498: #endif
1499: PetscCall(MatSeqAIJSetPreallocation(Prol, col_bs, NULL));
1500: PetscCall(MatMPIAIJSetPreallocation(Prol, col_bs, NULL, col_bs, NULL));
1502: /* can get all points "removed" */
1503: PetscCall(MatGetSize(Prol, &kk, &ii));
1504: if (!ii) {
1505: PetscCall(PetscInfo(pc, "%s: No selected points on coarse grid\n", ((PetscObject)pc)->prefix));
1506: PetscCall(MatDestroy(&Prol));
1507: *a_P_out = NULL; /* out */
1508: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1509: PetscFunctionReturn(PETSC_SUCCESS);
1510: }
1511: PetscCall(PetscInfo(pc, "%s: New grid %" PetscInt_FMT " nodes\n", ((PetscObject)pc)->prefix, ii / col_bs));
1512: PetscCall(MatGetOwnershipRangeColumn(Prol, &myCrs0, &kk));
1514: 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);
1515: myCrs0 = myCrs0 / col_bs;
1516: 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);
1518: /* create global vector of data in 'data_w_ghost' */
1519: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0));
1520: if (size > 1) { /* get ghost null space data */
1521: PetscReal *tmp_gdata, *tmp_ldata, *tp2;
1523: PetscCall(PetscMalloc1(nloc, &tmp_ldata));
1524: for (jj = 0; jj < col_bs; jj++) {
1525: for (kk = 0; kk < bs; kk++) {
1526: PetscInt stride;
1527: const PetscReal *tp = PetscSafePointerPlusOffset(pc_gamg->data, jj * bs * nloc + kk);
1529: for (PetscInt ii = 0; ii < nloc; ii++, tp += bs) tmp_ldata[ii] = *tp;
1531: PetscCall(PCGAMGGetDataWithGhosts(Gmat, 1, tmp_ldata, &stride, &tmp_gdata));
1533: if (!jj && !kk) { /* now I know how many total nodes - allocate TODO: move below and do in one 'col_bs' call */
1534: PetscCall(PetscMalloc1(stride * bs * col_bs, &data_w_ghost));
1535: nbnodes = bs * stride;
1536: }
1537: tp2 = PetscSafePointerPlusOffset(data_w_ghost, jj * bs * stride + kk);
1538: for (PetscInt ii = 0; ii < stride; ii++, tp2 += bs) *tp2 = tmp_gdata[ii];
1539: PetscCall(PetscFree(tmp_gdata));
1540: }
1541: }
1542: PetscCall(PetscFree(tmp_ldata));
1543: } else {
1544: nbnodes = bs * nloc;
1545: data_w_ghost = pc_gamg->data;
1546: }
1548: /* get 'flid_fgid' TODO - move up to get 'stride' and do get null space data above in one step (jj loop) */
1549: if (size > 1) {
1550: PetscReal *fid_glid_loc, *fiddata;
1551: PetscInt stride;
1553: PetscCall(PetscMalloc1(nloc, &fid_glid_loc));
1554: for (kk = 0; kk < nloc; kk++) fid_glid_loc[kk] = (PetscReal)(my0 + kk);
1555: PetscCall(PCGAMGGetDataWithGhosts(Gmat, 1, fid_glid_loc, &stride, &fiddata));
1556: PetscCall(PetscMalloc1(stride, &flid_fgid)); /* copy real data to in */
1557: for (kk = 0; kk < stride; kk++) flid_fgid[kk] = (PetscInt)fiddata[kk];
1558: PetscCall(PetscFree(fiddata));
1560: PetscCheck(stride == nbnodes / bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "stride %" PetscInt_FMT " != nbnodes %" PetscInt_FMT "/bs %" PetscInt_FMT, stride, nbnodes, bs);
1561: PetscCall(PetscFree(fid_glid_loc));
1562: } else {
1563: PetscCall(PetscMalloc1(nloc, &flid_fgid));
1564: for (kk = 0; kk < nloc; kk++) flid_fgid[kk] = my0 + kk;
1565: }
1566: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0));
1567: /* get P0 */
1568: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0));
1569: {
1570: PetscReal *data_out = NULL;
1572: PetscCall(formProl0(agg_lists, bs, col_bs, myCrs0, nbnodes, data_w_ghost, flid_fgid, &data_out, Prol));
1573: PetscCall(PetscFree(pc_gamg->data));
1575: pc_gamg->data = data_out;
1576: pc_gamg->data_cell_rows = col_bs;
1577: pc_gamg->data_sz = col_bs * col_bs * nLocalSelected;
1578: }
1579: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0));
1580: if (size > 1) PetscCall(PetscFree(data_w_ghost));
1581: PetscCall(PetscFree(flid_fgid));
1583: *a_P_out = Prol; /* out */
1584: PetscCall(MatViewFromOptions(Prol, NULL, "-pc_gamg_agg_view_initial_prolongation"));
1586: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1587: PetscFunctionReturn(PETSC_SUCCESS);
1588: }
1590: // Drop small node-coupling blocks of the prolongator; see the PCGAMGSetProlongatorFilter() manual page
1591: static PetscErrorCode PCGAMGProlongatorBlockFilter_AGG(PC pc, Mat Prol, PetscInt col_bs, PetscReal thr)
1592: {
1593: PetscInt rbs, cbs, rStart, rEnd, nfn, local_nnz = 0, max_fn_cols = 0, max_row_cols = 0, ndrows = 0, zoff = 0;
1594: PetscInt *cn_gid, *fn_col, *fn_slot, *fn_rcnt, *zcols, *drow, *doff, *dcnt;
1595: PetscReal *cn_n2;
1596: PetscReal thr2 = thr * thr;
1597: PetscScalar *zeros;
1598: PetscBool no_off_proc;
1600: PetscFunctionBegin;
1601: PetscCall(MatGetBlockSizes(Prol, &rbs, &cbs));
1602: /* Prol is built internally with MatSetBlockSizes(..., col_bs); a mismatch here means smoothing corrupted the block structure */
1603: 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);
1604: PetscCall(MatGetOwnershipRange(Prol, &rStart, &rEnd));
1605: 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);
1606: nfn = (rEnd - rStart) / rbs;
1608: /* Pre-pass over the row widths to size the scratch tightly: the total local nonzeros bound
1609: the dropped columns (zcols), the widest fine-node block row bounds the distinct coarse
1610: nodes of one fine node (cn_gid/cn_n2), and the widest single row bounds the entries
1611: zeroed by one MatSetValues() call (zeros) */
1612: for (PetscInt fn = 0; fn < nfn; fn++) {
1613: PetscInt fn_cols = 0;
1615: for (PetscInt rr = 0; rr < rbs; rr++) {
1616: PetscInt grow = rStart + fn * rbs + rr, ncols;
1618: PetscCall(MatGetRow(Prol, grow, &ncols, NULL, NULL));
1619: fn_cols += ncols;
1620: if (ncols > max_row_cols) max_row_cols = ncols;
1621: PetscCall(MatRestoreRow(Prol, grow, &ncols, NULL, NULL));
1622: }
1623: local_nnz += fn_cols;
1624: if (fn_cols > max_fn_cols) max_fn_cols = fn_cols;
1625: }
1627: 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));
1628: PetscCall(PetscCalloc1(max_row_cols, &zeros));
1629: PetscCall(PetscMalloc1(local_nnz, &zcols));
1630: PetscCall(PetscMalloc3(rEnd - rStart, &drow, rEnd - rStart, &doff, rEnd - rStart, &dcnt));
1632: for (PetscInt fn = 0; fn < nfn; fn++) {
1633: PetscInt ncn = 0; /* distinct coarse nodes touched by this fine node-block */
1634: PetscInt nent = 0; /* entries of this fine node-block cached in fn_col/fn_slot */
1635: PetscReal maxn2 = 0.0;
1637: /* Pass A: accumulate the per-coarse-node block Frobenius norm^2 over the rbs rows, caching each
1638: entry's column and resolved coarse-node slot so that Pass B needs no second MatGetRow() sweep */
1639: for (PetscInt rr = 0; rr < rbs; rr++) {
1640: PetscInt grow = rStart + fn * rbs + rr, ncols;
1641: const PetscInt *cols;
1642: const PetscScalar *vals;
1644: PetscCall(MatGetRow(Prol, grow, &ncols, &cols, &vals));
1645: fn_rcnt[rr] = ncols;
1646: for (PetscInt k = 0; k < ncols; k++) {
1647: PetscInt cn = cols[k] / col_bs, s = -1;
1648: PetscReal av = PetscAbsScalar(vals[k]);
1650: /* Linear dedup scan: ncn is bounded by the coarse-node degree of this
1651: fine node (typically ~5-20 for AMG-coarsened elasticity), so O(ncn)
1652: per entry is faster than a hash map at this scale */
1653: for (PetscInt t = 0; t < ncn; t++)
1654: if (cn_gid[t] == cn) {
1655: s = t;
1656: break;
1657: }
1658: if (s < 0) {
1659: s = ncn++;
1660: cn_gid[s] = cn;
1661: cn_n2[s] = 0.0;
1662: }
1663: cn_n2[s] += av * av;
1664: fn_col[nent] = cols[k];
1665: fn_slot[nent] = s;
1666: nent++;
1667: }
1668: PetscCall(MatRestoreRow(Prol, grow, &ncols, &cols, &vals));
1669: }
1670: for (PetscInt t = 0; t < ncn; t++)
1671: if (cn_n2[t] > maxn2) maxn2 = cn_n2[t];
1673: /* Pass B: collect the entries of blocks below thr*max for zeroing. The test is strict, so the
1674: strongest block of the node survives for any thr <= 1 */
1675: nent = 0;
1676: for (PetscInt rr = 0; rr < rbs; rr++) {
1677: PetscInt rec_off = zoff;
1679: for (PetscInt k = 0; k < fn_rcnt[rr]; k++, nent++)
1680: if (cn_n2[fn_slot[nent]] < thr2 * maxn2) zcols[zoff++] = fn_col[nent];
1681: if (zoff > rec_off) {
1682: drow[ndrows] = rStart + fn * rbs + rr;
1683: doff[ndrows] = rec_off;
1684: dcnt[ndrows] = zoff - rec_off;
1685: ndrows++;
1686: }
1687: }
1688: }
1690: /* all insertions are in local rows; skip the off-process assembly communication */
1691: PetscCall(MatGetOption(Prol, MAT_NO_OFF_PROC_ENTRIES, &no_off_proc));
1692: PetscCall(MatSetOption(Prol, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
1693: for (PetscInt i = 0; i < ndrows; i++) PetscCall(MatSetValues(Prol, 1, &drow[i], dcnt[i], &zcols[doff[i]], zeros, INSERT_VALUES));
1694: PetscCall(MatAssemblyBegin(Prol, MAT_FINAL_ASSEMBLY));
1695: PetscCall(MatAssemblyEnd(Prol, MAT_FINAL_ASSEMBLY));
1696: PetscCall(MatSetOption(Prol, MAT_NO_OFF_PROC_ENTRIES, no_off_proc));
1698: PetscCall(PetscFree5(cn_gid, cn_n2, fn_col, fn_slot, fn_rcnt));
1699: PetscCall(PetscFree(zeros));
1700: PetscCall(PetscFree(zcols));
1701: PetscCall(PetscFree3(drow, doff, dcnt));
1703: /* Compress out the explicit zeros just set, so the filter actually sparsifies. keep must be
1704: PETSC_FALSE: with keep, a zero whose local column index equals its local row index survives (an
1705: index-based diagonal test that is meaningless for the rectangular Prol). The compression is done
1706: in place, deliberately: a MatDuplicate()/MatHeaderReplace() copy as in MatFilter() would return
1707: the freed CSR tail to the allocator, at the cost of a peak-memory spike. MatEliminateZeros() is not
1708: implemented by all matrix types; in those cases the zeros are left in the sparsity pattern;
1709: the step-3 correction in PCGAMGKernelPreservingFilter_AGG() skips exactly-zero entries, so a dropped
1710: block stays dropped either way, it just still costs storage here. */
1711: if (Prol->ops->eliminatezeros) PetscCall(MatEliminateZeros(Prol, PETSC_FALSE));
1712: else PetscCall(PetscInfo(pc, "PCGAMGProlongatorBlockFilter_AGG: skipping zero elimination for %s; filtered entries are zeroed but not removed\n", ((PetscObject)Prol)->type_name));
1713: PetscFunctionReturn(PETSC_SUCCESS);
1714: }
1716: // 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
1717: static PetscErrorCode PCGAMGKernelPreservingFilter_AGG(PC pc, Mat Prol, PetscReal threshold)
1718: {
1719: PC_MG *mg = (PC_MG *)pc->data;
1720: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1721: const PetscInt nSAvec = pc_gamg->data_cell_rows; /* == data_cell_cols after formProl0 */
1722: PetscInt cStart, cEnd, rStart, rEnd;
1723: const PetscReal *Bc_data = pc_gamg->data;
1724: Vec *Bc_vecs, *B_vecs;
1725: PetscScalar *Bc_arr;
1727: PetscFunctionBegin;
1728: PetscCall(PetscInfo(pc, "Kernel-preserving filter of prolongator with threshold %g, nSAvec=%" PetscInt_FMT "\n", (double)threshold, nSAvec));
1730: PetscCall(MatGetOwnershipRange(Prol, &rStart, &rEnd));
1731: PetscCall(MatGetOwnershipRangeColumn(Prol, &cStart, &cEnd));
1733: /* Step 1: build coarse null-space vectors and compute B = P_original * B_c */
1734: PetscCall(PetscMalloc1(nSAvec, &Bc_vecs));
1735: PetscCall(PetscMalloc1(nSAvec, &B_vecs));
1737: {
1738: PetscInt nloc = cEnd - cStart;
1739: for (PetscInt k = 0; k < nSAvec; k++) {
1740: PetscCall(MatCreateVecs(Prol, &Bc_vecs[k], &B_vecs[k]));
1741: /* fill local entries: Bc_data layout is Bc_data[k * nloc + c] (stride == nloc) */
1742: PetscCall(VecGetArray(Bc_vecs[k], &Bc_arr));
1743: for (PetscInt c = 0; c < nloc; c++) Bc_arr[c] = (PetscScalar)Bc_data[k * nloc + c];
1744: PetscCall(VecRestoreArray(Bc_vecs[k], &Bc_arr));
1745: PetscCall(MatMult(Prol, Bc_vecs[k], B_vecs[k]));
1746: }
1747: }
1749: /* Step 2: apply the block-aware threshold filter (drops whole node-coupling
1750: blocks, preserving the coarse-node block structure; see
1751: PCGAMGProlongatorBlockFilter_AGG()) */
1752: {
1753: PetscBool info_active = PETSC_FALSE;
1754: MatInfo info0, info1;
1755: PetscCall(PetscInfoEnabled(((PetscObject)pc)->classid, &info_active));
1756: if (info_active) PetscCall(MatGetInfo(Prol, MAT_GLOBAL_SUM, &info0));
1757: PetscCall(PCGAMGProlongatorBlockFilter_AGG(pc, Prol, nSAvec, threshold));
1758: if (info_active) {
1759: PetscCall(MatGetInfo(Prol, MAT_GLOBAL_SUM, &info1));
1760: 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));
1761: }
1762: }
1764: /* Step 3: correct rows to restore P_filtered * B_c = B */
1765: if (nSAvec == 1) {
1766: /*
1767: Scalar case: use `MatMult()` + element-wise scaling + `MatDiagonalScale()`.
1768: scale_i = B_i / (P_filtered * Bc)_i, then P_new = diag(scale) * P_filtered.
1769: A row is left unscaled, that is, the near-null space constraint is not enforced
1770: for it, in the two cases where scaling would do more harm than the constraint is
1771: worth: when |scale_i| would exceed smax, which covers a denominator that is zero
1772: (an empty row after the filter) or tiny relative to B_i, and when B_i is zero,
1773: where enforcing the constraint would zero the whole row of the prolongator.
1774: Zeros left in the sparsity pattern stay zero under `MatDiagonalScale()`.
1775: No ghost column access needed.
1776: */
1777: Vec d_vec, scale_vec;
1778: PetscInt n_local, n_unscaled = 0;
1779: PetscReal smax = 1.0e4; /* cap on the row scale; the filter drops small blocks, so a healthy row has scale ~1 */
1780: PetscScalar *s_arr;
1781: const PetscScalar *b_arr, *d_arr;
1783: PetscCall(MatCreateVecs(Prol, NULL, &d_vec));
1784: PetscCall(MatMult(Prol, Bc_vecs[0], d_vec));
1785: PetscCall(VecDuplicate(d_vec, &scale_vec));
1786: PetscCall(VecGetLocalSize(d_vec, &n_local));
1787: PetscCall(VecGetArrayRead(B_vecs[0], &b_arr));
1788: PetscCall(VecGetArrayRead(d_vec, &d_arr));
1789: PetscCall(VecGetArray(scale_vec, &s_arr));
1790: for (PetscInt i = 0; i < n_local; i++) {
1791: PetscReal b = PetscAbsScalar(b_arr[i]), d = PetscAbsScalar(d_arr[i]);
1793: if (b > 0.0 && smax * d > b) s_arr[i] = b_arr[i] / d_arr[i];
1794: else {
1795: s_arr[i] = 1.0;
1796: n_unscaled++;
1797: }
1798: }
1799: 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));
1800: PetscCall(VecRestoreArray(scale_vec, &s_arr));
1801: PetscCall(VecRestoreArrayRead(d_vec, &d_arr));
1802: PetscCall(VecRestoreArrayRead(B_vecs[0], &b_arr));
1803: PetscCall(MatDiagonalScale(Prol, scale_vec, NULL));
1804: PetscCall(VecDestroy(&d_vec));
1805: PetscCall(VecDestroy(&scale_vec));
1806: } else {
1807: /*
1808: Vector case (nSAvec > 1): per-row minimum-norm correction (Gram-matrix solve).
1809: Scatter Bc_data to include ghost column values using Prol's Mvctx,
1810: then build a hash map from global ghost column index to local ghost index
1811: so that `MatGetRow()` global column indices can be mapped to the ghosted array.
1812: */
1813: PetscInt nloc = cEnd - cStart;
1814: PetscInt ghost_stride;
1815: PetscReal *Bc_ghosted = NULL;
1816: const PetscReal *Bc_ghosted_ro;
1817: PetscMPIInt comm_size;
1818: PetscHMapI ghost_gid_to_lid; /* global ghost col index -> local ghost index (0-based) */
1819: PetscInt num_ghosts = 0;
1820: PetscBool no_off_proc;
1822: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)Prol), &comm_size));
1823: if (comm_size > 1) {
1824: Mat_MPIAIJ *mpimat = (Mat_MPIAIJ *)Prol->data;
1825: Vec tmp_vec;
1826: PetscScalar *data_arr;
1827: PetscInt nnodes;
1828: PetscBool isaij;
1830: PetscCall(PetscObjectBaseTypeCompare((PetscObject)Prol, MATMPIAIJ, &isaij));
1831: PetscCheck(isaij, PetscObjectComm((PetscObject)Prol), PETSC_ERR_SUP, "Prolongator filter requires an MPIAIJ-based prolongator, not %s", ((PetscObject)Prol)->type_name);
1832: PetscCall(VecGetLocalSize(mpimat->lvec, &num_ghosts));
1833: nnodes = nloc + num_ghosts;
1834: ghost_stride = nnodes;
1835: /*
1836: Scatter Bc_data to include ghost column values using Prol's Mvctx.
1837: Cannot use PCGAMGGetDataWithGhosts() because it assumes square matrix
1838: (uses MatGetOwnershipRange() for row indices, but Prol is rectangular).
1839: */
1840: PetscCall(MatCreateVecs(Prol, &tmp_vec, NULL));
1841: PetscCall(PetscMalloc1(nSAvec * nnodes, &Bc_ghosted));
1842: for (PetscInt dir = 0; dir < nSAvec; dir++) {
1843: PetscScalar *tmp_arr;
1844: PetscCall(VecGetArray(tmp_vec, &tmp_arr));
1845: for (PetscInt kk = 0; kk < nloc; kk++) {
1846: PetscReal val = Bc_data[dir * nloc + kk];
1847: Bc_ghosted[dir * nnodes + kk] = val;
1848: tmp_arr[kk] = (PetscScalar)val;
1849: }
1850: PetscCall(VecRestoreArray(tmp_vec, &tmp_arr));
1851: PetscCall(VecScatterBegin(mpimat->Mvctx, tmp_vec, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD));
1852: PetscCall(VecScatterEnd(mpimat->Mvctx, tmp_vec, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD));
1853: PetscCall(VecGetArray(mpimat->lvec, &data_arr));
1854: for (PetscInt g = 0; g < num_ghosts; g++) Bc_ghosted[dir * nnodes + nloc + g] = PetscRealPart(data_arr[g]);
1855: PetscCall(VecRestoreArray(mpimat->lvec, &data_arr));
1856: }
1857: PetscCall(VecDestroy(&tmp_vec));
1858: Bc_ghosted_ro = Bc_ghosted;
1859: /* build hash: global ghost col index -> local ghost index (0-based into ghost portion) */
1860: PetscCall(PetscHMapICreateWithSize(2 * num_ghosts + 1, &ghost_gid_to_lid));
1861: for (PetscInt g = 0; g < num_ghosts; g++) PetscCall(PetscHMapISet(ghost_gid_to_lid, mpimat->garray[g], g));
1862: } else {
1863: /* sequential: no ghosts, ghost_stride == nloc, use Bc_data directly (read-only) */
1864: ghost_stride = nloc;
1865: Bc_ghosted_ro = Bc_data;
1866: PetscCall(PetscHMapICreateWithSize(1, &ghost_gid_to_lid));
1867: }
1869: {
1870: PetscInt nrows = rEnd - rStart, max_ncols = 0;
1871: const PetscScalar **B_arrays;
1872: PetscScalar *work, *new_vals, *G, *rhs, *x, *bc_col;
1873: PetscReal *dscale;
1874: PetscInt *ghosted_idx, *act, *col_buf;
1875: PetscBLASInt *ipiv;
1876: PetscBLASInt N_b;
1878: PetscCall(PetscMalloc1(nSAvec, &B_arrays));
1879: for (PetscInt k = 0; k < nSAvec; k++) PetscCall(VecGetArrayRead(B_vecs[k], &B_arrays[k]));
1880: /* work: nSAvec*nSAvec Gram + nSAvec rhs + nSAvec solution + nSAvec bc_col scratch */
1881: PetscCall(PetscMalloc1(nSAvec * nSAvec + 3 * nSAvec, &work));
1882: PetscCall(PetscMalloc1(nSAvec, &dscale));
1883: PetscCall(PetscMalloc1(nSAvec, &ipiv));
1884: PetscCall(PetscBLASIntCast(nSAvec, &N_b));
1885: G = work;
1886: rhs = work + nSAvec * nSAvec;
1887: x = rhs + nSAvec;
1888: bc_col = x + nSAvec;
1890: /* find max row width and total nnz for pre-allocation */
1891: {
1892: PetscInt total_nnz = 0;
1893: for (PetscInt row = 0; row < nrows; row++) {
1894: PetscInt ncols;
1895: PetscCall(MatGetRow(Prol, rStart + row, &ncols, NULL, NULL));
1896: if (ncols > max_ncols) max_ncols = ncols;
1897: total_nnz += ncols;
1898: PetscCall(MatRestoreRow(Prol, rStart + row, &ncols, NULL, NULL));
1899: }
1900: /* allocate flat CSR-like buffers to store all corrections before applying */
1901: PetscCall(PetscMalloc1(total_nnz, &new_vals));
1902: PetscCall(PetscMalloc1(total_nnz, &col_buf));
1903: }
1904: PetscCall(PetscMalloc2(max_ncols, &ghosted_idx, max_ncols, &act));
1906: /* Pass 1: read rows, compute corrections, store in flat buffers */
1907: {
1908: PetscInt *row_offsets;
1909: PetscInt offset = 0, n_singular = 0, n_zero_rows = 0, n_corrected = 0, n_underdetermined = 0;
1910: PetscReal max_ynorm = 0.0;
1912: PetscCall(PetscMalloc1(nrows + 1, &row_offsets));
1913: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1915: for (PetscInt row = 0; row < nrows; row++) {
1916: PetscInt ncols, nact = 0, grow = rStart + row, roff = offset;
1917: PetscBLASInt NRHS = 1, LDA = N_b, LDB = N_b, info;
1918: const PetscInt *cols;
1919: const PetscScalar *vals;
1921: row_offsets[row] = roff;
1922: PetscCall(MatGetRow(Prol, grow, &ncols, &cols, &vals));
1923: /* Save the row unchanged, then correct only its nonzero entries: the block filter zeroes
1924: dropped entries but they are not removed from the sparsity pattern for every matrix type
1925: (see PCGAMGProlongatorBlockFilter_AGG()), and adding the correction to them would turn a
1926: dropped block back into a nonzero one */
1927: for (PetscInt j = 0; j < ncols; j++) {
1928: col_buf[roff + j] = cols[j];
1929: new_vals[roff + j] = vals[j];
1930: if (vals[j] != 0.0) act[nact++] = j;
1931: }
1932: offset = roff + ncols;
1933: if (nact == 0) {
1934: n_zero_rows++;
1935: PetscCall(MatRestoreRow(Prol, grow, &ncols, &cols, &vals));
1936: continue;
1937: }
1939: /* When nact < nSAvec the Gram matrix G is rank-deficient by construction;
1940: skip correction for this row (keep filtered values as-is).
1941: Note: the near-null space constraint P*Bc = B is NOT enforced for these rows.
1942: This typically occurs at boundary or isolated nodes where few coarse neighbors
1943: remain after filtering; the impact on convergence is generally small. */
1944: if (nact < nSAvec) {
1945: n_underdetermined++;
1946: PetscCall(MatRestoreRow(Prol, grow, &ncols, &cols, &vals));
1947: continue;
1948: }
1950: /* map the global column indices of the surviving entries to ghosted array indices */
1951: for (PetscInt a = 0; a < nact; a++) {
1952: PetscInt col = cols[act[a]];
1953: if (col >= cStart && col < cEnd) ghosted_idx[a] = col - cStart;
1954: else {
1955: PetscInt g = -1;
1956: PetscCall(PetscHMapIGet(ghost_gid_to_lid, col, &g));
1957: PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Off-diagonal column %" PetscInt_FMT " not found in ghost map for prolongator filter", col);
1958: ghosted_idx[a] = nloc + g;
1959: }
1960: }
1962: for (PetscInt i = 0; i < nSAvec * nSAvec; i++) G[i] = 0.0;
1964: /* rhs[k] = B[row,k] - sum_a P[row,act[a]] * Bc[ghosted_idx[a], k] */
1965: for (PetscInt k = 0; k < nSAvec; k++) {
1966: PetscScalar dot = 0.0;
1967: for (PetscInt a = 0; a < nact; a++) dot += vals[act[a]] * (PetscScalar)Bc_ghosted_ro[k * ghost_stride + ghosted_idx[a]];
1968: rhs[k] = B_arrays[k][row] - dot;
1969: }
1971: /* G[k1,k2] = sum_a Bc[a,k1] * Bc[a,k2] using pre-gathered bc_col */
1972: for (PetscInt a = 0; a < nact; a++) {
1973: PetscInt gidx = ghosted_idx[a];
1974: for (PetscInt k = 0; k < nSAvec; k++) bc_col[k] = (PetscScalar)Bc_ghosted_ro[k * ghost_stride + gidx];
1975: for (PetscInt k1 = 0; k1 < nSAvec; k1++)
1976: for (PetscInt k2 = k1; k2 < nSAvec; k2++) G[k1 * nSAvec + k2] += bc_col[k1] * bc_col[k2];
1977: }
1978: /* fill lower triangle from upper (G is symmetric) */
1979: for (PetscInt k1 = 1; k1 < nSAvec; k1++)
1980: for (PetscInt k2 = 0; k2 < k1; k2++) G[k1 * nSAvec + k2] = G[k2 * nSAvec + k1];
1982: /* Symmetric (Jacobi) equilibration: G is severely ill-conditioned for
1983: elasticity because the near-null modes have disparate scales (O(1)
1984: translations vs rotations that scale with the coordinates). Scale by
1985: dscale[k] = 1/sqrt(G[k,k]) so the rescaled Gram has a unit diagonal,
1986: then solve (D G D) y = D rhs and recover x = D y. This is identical to
1987: solving G x = rhs in exact arithmetic but removes the mode-scale
1988: ill-conditioning, making the correction accurate and FP-robust. */
1989: for (PetscInt k = 0; k < nSAvec; k++) {
1990: PetscReal gkk = PetscRealPart(G[k * nSAvec + k]);
1991: dscale[k] = gkk > 0.0 ? 1.0 / PetscSqrtReal(gkk) : 1.0;
1992: }
1993: for (PetscInt k1 = 0; k1 < nSAvec; k1++)
1994: for (PetscInt k2 = 0; k2 < nSAvec; k2++) G[k1 * nSAvec + k2] *= dscale[k1] * dscale[k2];
1996: /* solve (D G D) y = D rhs, then x = D y */
1997: for (PetscInt k = 0; k < nSAvec; k++) x[k] = dscale[k] * rhs[k];
1998: PetscCallBLAS("LAPACKgesv", LAPACKgesv_(&N_b, &NRHS, G, &LDA, ipiv, x, &LDB, &info));
1999: PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "LAPACKgesv: %" PetscBLASInt_FMT "-th argument had an illegal value", -info);
2000: if (info > 0) {
2001: /* G is singular despite nact >= nSAvec (Bc columns linearly dependent);
2002: keep filtered values as-is (near-null space constraint not enforced for this row) */
2003: n_singular++;
2004: PetscCall(MatRestoreRow(Prol, grow, &ncols, &cols, &vals));
2005: continue;
2006: }
2007: /* track the equilibrated solution norm ||y||^2 (x still holds y here): a
2008: basis-independent health metric, large values flag a problematic row */
2009: {
2010: PetscReal ynorm2 = 0.0;
2011: for (PetscInt k = 0; k < nSAvec; k++) ynorm2 += PetscSqr(PetscAbsScalar(x[k]));
2012: if (ynorm2 > max_ynorm) max_ynorm = ynorm2;
2013: }
2014: for (PetscInt k = 0; k < nSAvec; k++) x[k] *= dscale[k]; /* recover x = D y */
2015: n_corrected++;
2017: /* new_vals[act[a]] = vals[act[a]] + sum_k Bc[ghosted_idx[a],k] * x[k] */
2018: for (PetscInt a = 0; a < nact; a++) {
2019: PetscScalar delta = 0.0;
2020: PetscInt gidx = ghosted_idx[a];
2021: for (PetscInt k = 0; k < nSAvec; k++) delta += (PetscScalar)Bc_ghosted_ro[k * ghost_stride + gidx] * x[k];
2022: new_vals[roff + act[a]] += delta;
2023: }
2024: PetscCall(MatRestoreRow(Prol, grow, &ncols, &cols, &vals));
2025: }
2026: row_offsets[nrows] = offset;
2027: PetscCall(PetscFPTrapPop());
2028: PetscCall(PetscInfo(pc, "PCGAMGKernelPreservingFilter_AGG: corrected %" PetscInt_FMT "/%" PetscInt_FMT " rows, max equilibrated correction ||y||^2=%g\n", n_corrected, nrows, (double)max_ynorm));
2029: if (n_zero_rows + n_underdetermined + n_singular > 0)
2030: 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));
2032: /* Pass 2: apply all corrections at once */
2033: for (PetscInt row = 0; row < nrows; row++) {
2034: PetscInt grow = rStart + row;
2035: PetscInt nc = row_offsets[row + 1] - row_offsets[row];
2036: if (nc > 0) PetscCall(MatSetValues(Prol, 1, &grow, nc, col_buf + row_offsets[row], new_vals + row_offsets[row], INSERT_VALUES));
2037: }
2038: PetscCall(PetscFree(row_offsets));
2039: }
2041: for (PetscInt k = 0; k < nSAvec; k++) PetscCall(VecRestoreArrayRead(B_vecs[k], &B_arrays[k]));
2042: PetscCall(PetscFree(B_arrays));
2043: PetscCall(PetscFree(work));
2044: PetscCall(PetscFree(dscale));
2045: PetscCall(PetscFree(ipiv));
2046: PetscCall(PetscFree2(ghosted_idx, act));
2047: PetscCall(PetscFree(new_vals));
2048: PetscCall(PetscFree(col_buf));
2049: }
2051: PetscCall(PetscHMapIDestroy(&ghost_gid_to_lid));
2052: if (comm_size > 1) PetscCall(PetscFree(Bc_ghosted));
2054: /* all insertions are in local rows; skip the off-process assembly communication. The scalar
2055: branch needs no assembly at all: MatDiagonalScale() requires and preserves assembly */
2056: PetscCall(MatGetOption(Prol, MAT_NO_OFF_PROC_ENTRIES, &no_off_proc));
2057: PetscCall(MatSetOption(Prol, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
2058: PetscCall(MatAssemblyBegin(Prol, MAT_FINAL_ASSEMBLY));
2059: PetscCall(MatAssemblyEnd(Prol, MAT_FINAL_ASSEMBLY));
2060: PetscCall(MatSetOption(Prol, MAT_NO_OFF_PROC_ENTRIES, no_off_proc));
2061: }
2063: for (PetscInt k = 0; k < nSAvec; k++) {
2064: PetscCall(VecDestroy(&Bc_vecs[k]));
2065: PetscCall(VecDestroy(&B_vecs[k]));
2066: }
2067: PetscCall(PetscFree(Bc_vecs));
2068: PetscCall(PetscFree(B_vecs));
2069: PetscFunctionReturn(PETSC_SUCCESS);
2070: }
2072: /*
2073: PCGAMGOptimizeProlongator_AGG - given the initial prolongator optimizes it by smoothed aggregation pc_gamg_agg->nsmooths times
2075: Input Parameter:
2076: . pc - this
2077: . Amat - matrix on this fine level
2078: In/Output Parameter:
2079: . a_P - prolongation operator to the next level
2080: */
2081: static PetscErrorCode PCGAMGOptimizeProlongator_AGG(PC pc, Mat Amat, Mat *a_P)
2082: {
2083: PC_MG *mg = (PC_MG *)pc->data;
2084: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
2085: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
2086: Mat Prol = *a_P;
2087: MPI_Comm comm;
2088: KSP eksp;
2089: Vec bb, xx;
2090: PC epc;
2091: PetscReal alpha, emax, emin;
2092: PetscReal pfilter = pc_gamg->prolongator_filter * PetscPowRealInt(pc_gamg->prolongator_filter_scale, pc_gamg->current_level);
2094: PetscFunctionBegin;
2095: PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
2096: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0));
2098: /* compute maximum singular value of operator to be used in smoother */
2099: if (0 < pc_gamg_agg->nsmooths) {
2100: /* get eigen estimates */
2101: if (pc_gamg->emax > 0) {
2102: emin = pc_gamg->emin;
2103: emax = pc_gamg->emax;
2104: } else {
2105: const char *prefix;
2107: PetscCall(MatCreateVecs(Amat, &bb, NULL));
2108: PetscCall(MatCreateVecs(Amat, &xx, NULL));
2109: PetscCall(KSPSetNoisy_Private(Amat, bb));
2111: PetscCall(KSPCreate(comm, &eksp));
2112: PetscCall(KSPSetNestLevel(eksp, pc->kspnestlevel));
2113: PetscCall(PCGetOptionsPrefix(pc, &prefix));
2114: PetscCall(KSPSetOptionsPrefix(eksp, prefix));
2115: PetscCall(KSPAppendOptionsPrefix(eksp, "pc_gamg_esteig_"));
2116: {
2117: PetscBool isset, sflg;
2119: PetscCall(MatIsSPDKnown(Amat, &isset, &sflg));
2120: if (isset && sflg) PetscCall(KSPSetType(eksp, KSPCG));
2121: }
2122: PetscCall(KSPSetErrorIfNotConverged(eksp, pc->erroriffailure));
2123: PetscCall(KSPSetNormType(eksp, KSP_NORM_NONE));
2125: PetscCall(KSPSetInitialGuessNonzero(eksp, PETSC_FALSE));
2126: PetscCall(KSPSetOperators(eksp, Amat, Amat));
2128: PetscCall(KSPGetPC(eksp, &epc));
2129: PetscCall(PCSetType(epc, PCJACOBI)); /* smoother in smoothed agg. */
2131: 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
2133: PetscCall(KSPSetFromOptions(eksp));
2134: PetscCall(KSPSetComputeSingularValues(eksp, PETSC_TRUE));
2135: PetscCall(KSPSolve(eksp, bb, xx));
2136: PetscCall(KSPCheckSolve(eksp, pc, xx));
2138: PetscCall(KSPComputeExtremeSingularValues(eksp, &emax, &emin));
2139: PetscCall(PetscInfo(pc, "%s: Smooth P0: max eigen=%e min=%e PC=%s\n", ((PetscObject)pc)->prefix, (double)emax, (double)emin, PCJACOBI));
2140: PetscCall(VecDestroy(&xx));
2141: PetscCall(VecDestroy(&bb));
2142: PetscCall(KSPDestroy(&eksp));
2143: }
2144: if (pc_gamg->use_sa_esteig) {
2145: mg->min_eigen_DinvA[pc_gamg->current_level] = emin;
2146: mg->max_eigen_DinvA[pc_gamg->current_level] = emax;
2147: 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));
2148: } else {
2149: mg->min_eigen_DinvA[pc_gamg->current_level] = 0;
2150: mg->max_eigen_DinvA[pc_gamg->current_level] = 0;
2151: }
2152: } else {
2153: mg->min_eigen_DinvA[pc_gamg->current_level] = 0;
2154: mg->max_eigen_DinvA[pc_gamg->current_level] = 0;
2155: }
2157: /* smooth P0 */
2158: if (pc_gamg_agg->nsmooths > 0) {
2159: Vec diag;
2161: /* TODO: Set a PCFailedReason and exit the building of the AMG preconditioner */
2162: PetscCheck(emax != 0.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Computed maximum singular value as zero");
2164: PetscCall(MatCreateVecs(Amat, &diag, NULL));
2165: PetscCall(MatGetDiagonal(Amat, diag)); /* effectively PCJACOBI */
2166: PetscCall(VecReciprocal(diag));
2168: for (PetscInt jj = 0; jj < pc_gamg_agg->nsmooths; jj++) {
2169: Mat tMat;
2171: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0));
2172: /*
2173: Smooth aggregation on the prolongator
2175: P_{i} := (I - 1.4/emax D^{-1}A) P_i\{i-1}
2176: */
2177: PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0));
2178: PetscCall(MatMatMult(Amat, Prol, MAT_INITIAL_MATRIX, PETSC_CURRENT, &tMat));
2179: PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0));
2180: PetscCall(MatProductClear(tMat));
2181: PetscCall(MatDiagonalScale(tMat, diag, NULL));
2183: /* TODO: Document the 1.4 and don't hardwire it in this routine */
2184: alpha = -1.4 / emax;
2185: PetscCall(MatAYPX(tMat, alpha, Prol, SUBSET_NONZERO_PATTERN));
2186: PetscCall(MatDestroy(&Prol));
2187: Prol = tMat;
2188: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0));
2189: }
2190: PetscCall(VecDestroy(&diag));
2191: }
2192: /* a per-level threshold of 0 (from prolongator_filter_scale == 0 on the coarser levels) drops
2193: nothing, so skip the whole filter rather than pay for its passes and per-row solves */
2194: if (pfilter > 0.0) {
2195: 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,
2196: (double)pc_gamg->prolongator_filter_scale, pc_gamg->current_level));
2197: PetscCall(PCGAMGKernelPreservingFilter_AGG(pc, Prol, pfilter));
2198: }
2199: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0));
2200: PetscCall(MatViewFromOptions(Prol, NULL, "-pc_gamg_agg_view_prolongation"));
2201: *a_P = Prol;
2202: PetscFunctionReturn(PETSC_SUCCESS);
2203: }
2205: /*MC
2206: PCGAMGAGG - Smooth aggregation, {cite}`vanek1996algebraic`, {cite}`vanek2001convergence`, variant of PETSc's algebraic multigrid (`PCGAMG`) preconditioner
2208: Options Database Keys:
2209: + -pc_gamg_agg_nsmooths nsmooth - number of smoothing steps to use with smooth aggregation to construct prolongation
2210: . -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)
2211: . -pc_gamg_prolongator_filter_scale scale - per-level scaling of the prolongator filter threshold (1.0=default)
2212: . -pc_gamg_aggressive_coarsening n - number of aggressive coarsening (MIS-2 or square graph) levels from finest.
2213: . -pc_gamg_aggressive_square_graph (true|false) - use square graph ($A^T A$), alternative is MIS-k (k=2), for aggressive coarsening
2214: . -pc_gamg_mis_k_minimum_degree_ordering (true|false) - use minimum degree ordering in greedy MIS algorithm
2215: . -pc_gamg_asm_hem_aggs n - number of HEM aggregation steps for ASM smoother
2216: - -pc_gamg_aggressive_mis_k n - number (k) distance in MIS coarsening (>2 is 'aggressive')
2218: Level: intermediate
2220: Notes:
2221: To obtain good performance for `PCGAMG` for vector valued problems you must
2222: call `MatSetBlockSize()` to indicate the number of degrees of freedom per grid point.
2223: Call `MatSetNearNullSpace()` (or `PCSetCoordinates()` if solving the equations of elasticity) to indicate the near null space of the operator
2225: When `-pc_gamg_aggressive_square_graph` is used, the coarsening is obtained by first squaring the graph and then applying, by default, a
2226: MIS-1 coarsening with `MatCoarsenApply()` on the squared graph.
2228: The many options for `PCMG` and `PCGAMG` such as controlling the smoothers on each level etc. also work for `PCGAMGAGG`
2230: .seealso: `PCGAMG`, [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCCreate()`, `PCSetType()`,
2231: `MatSetBlockSize()`, `PCMGType`, `PCSetCoordinates()`, `MatSetNearNullSpace()`, `PCGAMGSetType()`,
2232: `PCGAMGAGG`, `PCGAMGGEO`, `PCGAMGCLASSICAL`, `PCGAMGSetProcEqLim()`, `PCGAMGSetCoarseEqLim()`, `PCGAMGSetRepartition()`, `PCGAMGRegister()`,
2233: `PCGAMGSetReuseInterpolation()`, `PCGAMGASMSetUseAggs()`, `PCGAMGSetParallelCoarseGridSolve()`, `PCGAMGSetNlevels()`, `PCGAMGSetThreshold()`,
2234: `PCGAMGGetType()`, `PCGAMGSetUseSAEstEig()`
2235: M*/
2236: PetscErrorCode PCCreateGAMG_AGG(PC pc)
2237: {
2238: PC_MG *mg = (PC_MG *)pc->data;
2239: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
2240: PC_GAMG_AGG *pc_gamg_agg;
2242: PetscFunctionBegin;
2243: /* create sub context for SA */
2244: PetscCall(PetscNew(&pc_gamg_agg));
2245: pc_gamg->subctx = pc_gamg_agg;
2247: pc_gamg->ops->setfromoptions = PCSetFromOptions_GAMG_AGG;
2248: pc_gamg->ops->destroy = PCDestroy_GAMG_AGG;
2249: /* reset does not do anything; setup not virtual */
2251: /* set internal function pointers */
2252: pc_gamg->ops->creategraph = PCGAMGCreateGraph_AGG;
2253: pc_gamg->ops->coarsen = PCGAMGCoarsen_AGG;
2254: pc_gamg->ops->prolongator = PCGAMGConstructProlongator_AGG;
2255: pc_gamg->ops->optprolongator = PCGAMGOptimizeProlongator_AGG;
2256: pc_gamg->ops->createdefaultdata = PCSetData_AGG;
2257: pc_gamg->ops->view = PCView_GAMG_AGG;
2259: pc_gamg_agg->nsmooths = 1;
2260: pc_gamg_agg->aggressive_coarsening_levels = 1;
2261: pc_gamg_agg->use_aggressive_square_graph = PETSC_TRUE;
2262: pc_gamg_agg->use_minimum_degree_ordering = PETSC_FALSE;
2263: pc_gamg_agg->use_low_mem_filter = PETSC_FALSE;
2264: pc_gamg_agg->aggressive_mis_k = 2;
2265: pc_gamg_agg->graph_symmetrize = PETSC_TRUE;
2267: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", PCGAMGSetNSmooths_AGG));
2268: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", PCGAMGSetAggressiveLevels_AGG));
2269: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveSquareGraph_C", PCGAMGSetAggressiveSquareGraph_AGG));
2270: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetMinDegreeOrdering_C", PCGAMGMISkSetMinDegreeOrdering_AGG));
2271: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetLowMemoryFilter_C", PCGAMGSetLowMemoryFilter_AGG));
2272: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetAggressive_C", PCGAMGMISkSetAggressive_AGG));
2273: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetGraphSymmetrize_C", PCGAMGSetGraphSymmetrize_AGG));
2274: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetProlongatorFilter_C", PCGAMGSetProlongatorFilter_AGG));
2275: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGGetProlongatorFilter_C", PCGAMGGetProlongatorFilter_AGG));
2276: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetProlongatorFilterScale_C", PCGAMGSetProlongatorFilterScale_AGG));
2277: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGGetProlongatorFilterScale_C", PCGAMGGetProlongatorFilterScale_AGG));
2278: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_AGG));
2279: PetscFunctionReturn(PETSC_SUCCESS);
2280: }