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: static PetscErrorCode PCGAMGSetAggressiveLevels_AGG(PC pc, PetscInt n)
230: {
231: PC_MG *mg = (PC_MG *)pc->data;
232: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
233: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
235: PetscFunctionBegin;
236: pc_gamg_agg->aggressive_coarsening_levels = n;
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: static PetscErrorCode PCGAMGMISkSetAggressive_AGG(PC pc, PetscInt n)
241: {
242: PC_MG *mg = (PC_MG *)pc->data;
243: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
244: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
246: PetscFunctionBegin;
247: pc_gamg_agg->aggressive_mis_k = n;
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }
251: static PetscErrorCode PCGAMGSetAggressiveSquareGraph_AGG(PC pc, PetscBool b)
252: {
253: PC_MG *mg = (PC_MG *)pc->data;
254: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
255: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
257: PetscFunctionBegin;
258: pc_gamg_agg->use_aggressive_square_graph = b;
259: PetscFunctionReturn(PETSC_SUCCESS);
260: }
262: static PetscErrorCode PCGAMGSetLowMemoryFilter_AGG(PC pc, PetscBool b)
263: {
264: PC_MG *mg = (PC_MG *)pc->data;
265: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
266: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
268: PetscFunctionBegin;
269: pc_gamg_agg->use_low_mem_filter = b;
270: PetscFunctionReturn(PETSC_SUCCESS);
271: }
273: static PetscErrorCode PCGAMGSetGraphSymmetrize_AGG(PC pc, PetscBool b)
274: {
275: PC_MG *mg = (PC_MG *)pc->data;
276: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
277: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
279: PetscFunctionBegin;
280: pc_gamg_agg->graph_symmetrize = b;
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }
284: static PetscErrorCode PCGAMGMISkSetMinDegreeOrdering_AGG(PC pc, PetscBool b)
285: {
286: PC_MG *mg = (PC_MG *)pc->data;
287: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
288: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
290: PetscFunctionBegin;
291: pc_gamg_agg->use_minimum_degree_ordering = b;
292: PetscFunctionReturn(PETSC_SUCCESS);
293: }
295: static PetscErrorCode PCSetFromOptions_GAMG_AGG(PC pc, PetscOptionItems PetscOptionsObject)
296: {
297: PC_MG *mg = (PC_MG *)pc->data;
298: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
299: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
300: PetscBool n_aggressive_flg, old_sq_provided = PETSC_FALSE, new_sq_provided = PETSC_FALSE, new_sqr_graph = pc_gamg_agg->use_aggressive_square_graph;
301: PetscInt nsq_graph_old = 0;
303: PetscFunctionBegin;
304: PetscOptionsHeadBegin(PetscOptionsObject, "GAMG-AGG options");
305: 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));
306: // aggressive coarsening logic with deprecated -pc_gamg_square_graph
307: 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));
308: if (!n_aggressive_flg)
309: 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));
310: 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));
311: if (!new_sq_provided && old_sq_provided) {
312: pc_gamg_agg->aggressive_coarsening_levels = nsq_graph_old; // could be zero
313: pc_gamg_agg->use_aggressive_square_graph = PETSC_TRUE;
314: }
315: if (new_sq_provided && old_sq_provided)
316: 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));
317: 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));
318: 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));
319: 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));
320: PetscCall(PetscOptionsBool("-pc_gamg_graph_symmetrize", "Symmetrize graph for coarsening", "PCGAMGSetGraphSymmetrize", pc_gamg_agg->graph_symmetrize, &pc_gamg_agg->graph_symmetrize, NULL));
321: PetscOptionsHeadEnd();
322: PetscFunctionReturn(PETSC_SUCCESS);
323: }
325: static PetscErrorCode PCDestroy_GAMG_AGG(PC pc)
326: {
327: PC_MG *mg = (PC_MG *)pc->data;
328: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
329: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
331: PetscFunctionBegin;
332: PetscCall(MatCoarsenDestroy(&pc_gamg_agg->crs));
333: PetscCall(PetscFree(pc_gamg->subctx));
334: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", NULL));
335: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", NULL));
336: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetAggressive_C", NULL));
337: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetMinDegreeOrdering_C", NULL));
338: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetLowMemoryFilter_C", NULL));
339: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveSquareGraph_C", NULL));
340: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetGraphSymmetrize_C", NULL));
341: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
342: PetscFunctionReturn(PETSC_SUCCESS);
343: }
345: /*
346: PCSetCoordinates_AGG
348: Collective
350: Input Parameter:
351: . pc - the preconditioner context
352: . ndm - dimension of data (used for dof/vertex for Stokes)
353: . a_nloc - number of vertices local
354: . coords - [a_nloc][ndm] - interleaved coordinate data: {x_0, y_0, z_0, x_1, y_1, ...}
355: */
357: static PetscErrorCode PCSetCoordinates_AGG(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords)
358: {
359: PC_MG *mg = (PC_MG *)pc->data;
360: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
361: PetscInt arrsz, kk, ii, jj, nloc, ndatarows, ndf;
362: Mat mat = pc->pmat;
364: PetscFunctionBegin;
367: nloc = a_nloc;
369: /* SA: null space vectors */
370: PetscCall(MatGetBlockSize(mat, &ndf)); /* this does not work for Stokes */
371: if (coords && ndf == 1) pc_gamg->data_cell_cols = 1; /* scalar w/ coords and SA (not needed) */
372: else if (coords) {
373: PetscCheck(ndm <= ndf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "degrees of motion %" PetscInt_FMT " > block size %" PetscInt_FMT, ndm, ndf);
374: pc_gamg->data_cell_cols = (ndm == 2 ? 3 : 6); /* displacement elasticity */
375: 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);
376: } else pc_gamg->data_cell_cols = ndf; /* no data, force SA with constant null space vectors */
377: pc_gamg->data_cell_rows = ndatarows = ndf;
378: 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);
379: arrsz = nloc * pc_gamg->data_cell_rows * pc_gamg->data_cell_cols;
381: if (!pc_gamg->data || (pc_gamg->data_sz != arrsz)) {
382: PetscCall(PetscFree(pc_gamg->data));
383: PetscCall(PetscMalloc1(arrsz + 1, &pc_gamg->data));
384: }
385: /* copy data in - column-oriented */
386: for (kk = 0; kk < nloc; kk++) {
387: const PetscInt M = nloc * pc_gamg->data_cell_rows; /* stride into data */
388: PetscReal *data = &pc_gamg->data[kk * ndatarows]; /* start of cell */
390: if (pc_gamg->data_cell_cols == 1) *data = 1.0;
391: else {
392: /* translational modes */
393: for (ii = 0; ii < ndatarows; ii++) {
394: for (jj = 0; jj < ndatarows; jj++) {
395: if (ii == jj) data[ii * M + jj] = 1.0;
396: else data[ii * M + jj] = 0.0;
397: }
398: }
400: /* rotational modes */
401: if (coords) {
402: if (ndm == 2) {
403: data += 2 * M;
404: data[0] = -coords[2 * kk + 1];
405: data[1] = coords[2 * kk];
406: } else {
407: data += 3 * M;
408: data[0] = 0.0;
409: data[M + 0] = coords[3 * kk + 2];
410: data[2 * M + 0] = -coords[3 * kk + 1];
411: data[1] = -coords[3 * kk + 2];
412: data[M + 1] = 0.0;
413: data[2 * M + 1] = coords[3 * kk];
414: data[2] = coords[3 * kk + 1];
415: data[M + 2] = -coords[3 * kk];
416: data[2 * M + 2] = 0.0;
417: }
418: }
419: }
420: }
421: pc_gamg->data_sz = arrsz;
422: PetscFunctionReturn(PETSC_SUCCESS);
423: }
425: /*
426: PCSetData_AGG - called if data is not set with PCSetCoordinates.
427: Looks in Mat for near null space.
428: Does not work for Stokes
430: Input Parameter:
431: . pc -
432: . a_A - matrix to get (near) null space out of.
433: */
434: static PetscErrorCode PCSetData_AGG(PC pc, Mat a_A)
435: {
436: PC_MG *mg = (PC_MG *)pc->data;
437: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
438: MatNullSpace mnull;
440: PetscFunctionBegin;
441: PetscCall(MatGetNearNullSpace(a_A, &mnull));
442: if (!mnull) {
443: DM dm;
445: PetscCall(PCGetDM(pc, &dm));
446: if (!dm) PetscCall(MatGetDM(a_A, &dm));
447: if (dm) {
448: PetscObject deformation;
449: PetscInt Nf;
451: PetscCall(DMGetNumFields(dm, &Nf));
452: if (Nf) {
453: PetscCall(DMGetField(dm, 0, NULL, &deformation));
454: if (deformation) {
455: PetscCall(PetscObjectQuery(deformation, "nearnullspace", (PetscObject *)&mnull));
456: if (!mnull) PetscCall(PetscObjectQuery(deformation, "nullspace", (PetscObject *)&mnull));
457: }
458: }
459: }
460: }
462: if (!mnull) {
463: PetscInt bs, NN, MM;
465: PetscCall(MatGetBlockSize(a_A, &bs));
466: PetscCall(MatGetLocalSize(a_A, &MM, &NN));
467: PetscCheck(MM % bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MM %" PetscInt_FMT " must be divisible by bs %" PetscInt_FMT, MM, bs);
468: PetscCall(PCSetCoordinates_AGG(pc, bs, MM / bs, NULL));
469: } else {
470: PetscReal *nullvec;
471: PetscBool has_const;
472: PetscInt i, j, mlocal, nvec, bs;
473: const Vec *vecs;
474: const PetscScalar *v;
476: PetscCall(MatGetLocalSize(a_A, &mlocal, NULL));
477: PetscCall(MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs));
478: for (i = 0; i < nvec; i++) {
479: PetscCall(VecGetLocalSize(vecs[i], &j));
480: PetscCheck(j == mlocal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Attached null space vector size %" PetscInt_FMT " != matrix size %" PetscInt_FMT, j, mlocal);
481: }
482: pc_gamg->data_sz = (nvec + !!has_const) * mlocal;
483: PetscCall(PetscMalloc1((nvec + !!has_const) * mlocal, &nullvec));
484: if (has_const)
485: for (i = 0; i < mlocal; i++) nullvec[i] = 1.0;
486: for (i = 0; i < nvec; i++) {
487: PetscCall(VecGetArrayRead(vecs[i], &v));
488: for (j = 0; j < mlocal; j++) nullvec[(i + !!has_const) * mlocal + j] = PetscRealPart(v[j]);
489: PetscCall(VecRestoreArrayRead(vecs[i], &v));
490: }
491: pc_gamg->data = nullvec;
492: pc_gamg->data_cell_cols = (nvec + !!has_const);
493: PetscCall(MatGetBlockSize(a_A, &bs));
494: pc_gamg->data_cell_rows = bs;
495: }
496: PetscFunctionReturn(PETSC_SUCCESS);
497: }
499: /*
500: formProl0 - collect null space data for each aggregate, do QR, put R in coarse grid data and Q in P_0
502: Input Parameter:
503: . agg_llists - list of arrays with aggregates -- list from selected vertices of aggregate unselected vertices
504: . bs - row block size
505: . nSAvec - column bs of new P
506: . my0crs - global index of start of locals
507: . data_stride - bs*(nloc nodes + ghost nodes) [data_stride][nSAvec]
508: . data_in[data_stride*nSAvec] - local data on fine grid
509: . flid_fgid[data_stride/bs] - make local to global IDs, includes ghosts in 'locals_llist'
511: Output Parameter:
512: . a_data_out - in with fine grid data (w/ghosts), out with coarse grid data
513: . a_Prol - prolongation operator
514: */
515: 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)
516: {
517: PetscInt Istart, my0, Iend, nloc, clid, flid = 0, aggID, kk, jj, ii, mm, nSelected, minsz, nghosts, out_data_stride;
518: MPI_Comm comm;
519: PetscReal *out_data;
520: PetscCDIntNd *pos;
521: PetscHMapI fgid_flid;
523: PetscFunctionBegin;
524: PetscCall(PetscObjectGetComm((PetscObject)a_Prol, &comm));
525: PetscCall(MatGetOwnershipRange(a_Prol, &Istart, &Iend));
526: nloc = (Iend - Istart) / bs;
527: my0 = Istart / bs;
528: 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);
529: Iend /= bs;
530: nghosts = data_stride / bs - nloc;
532: PetscCall(PetscHMapICreateWithSize(2 * nghosts + 1, &fgid_flid));
534: for (kk = 0; kk < nghosts; kk++) PetscCall(PetscHMapISet(fgid_flid, flid_fgid[nloc + kk], nloc + kk));
536: /* count selected -- same as number of cols of P */
537: for (nSelected = mm = 0; mm < nloc; mm++) {
538: PetscBool ise;
540: PetscCall(PetscCDIsEmptyAt(agg_llists, mm, &ise));
541: if (!ise) nSelected++;
542: }
543: PetscCall(MatGetOwnershipRangeColumn(a_Prol, &ii, &jj));
544: PetscCheck((ii / nSAvec) == my0crs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ii %" PetscInt_FMT " /nSAvec %" PetscInt_FMT " != my0crs %" PetscInt_FMT, ii, nSAvec, my0crs);
545: 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);
547: /* aloc space for coarse point data (output) */
548: out_data_stride = nSelected * nSAvec;
550: PetscCall(PetscMalloc1(out_data_stride * nSAvec, &out_data));
551: for (ii = 0; ii < out_data_stride * nSAvec; ii++) out_data[ii] = PETSC_MAX_REAL;
552: *a_data_out = out_data; /* output - stride nSelected*nSAvec */
554: /* find points and set prolongation */
555: minsz = 100;
556: for (mm = clid = 0; mm < nloc; mm++) {
557: PetscCall(PetscCDCountAt(agg_llists, mm, &jj));
558: if (jj > 0) {
559: const PetscInt lid = mm, cgid = my0crs + clid;
560: PetscInt cids[100]; /* max bs */
561: PetscBLASInt asz, M, N, INFO;
562: PetscBLASInt Mdata, LDA, LWORK;
563: PetscScalar *qqc, *qqr, *TAU, *WORK;
564: PetscInt *fids;
565: PetscReal *data;
567: PetscCall(PetscBLASIntCast(jj, &asz));
568: PetscCall(PetscBLASIntCast(asz * bs, &M));
569: PetscCall(PetscBLASIntCast(nSAvec, &N));
570: PetscCall(PetscBLASIntCast(M + ((N - M > 0) ? N - M : 0), &Mdata));
571: PetscCall(PetscBLASIntCast(Mdata, &LDA));
572: PetscCall(PetscBLASIntCast(N * bs, &LWORK));
573: /* count agg */
574: if (asz < minsz) minsz = asz;
576: /* get block */
577: PetscCall(PetscMalloc5(Mdata * N, &qqc, M * N, &qqr, N, &TAU, LWORK, &WORK, M, &fids));
579: aggID = 0;
580: PetscCall(PetscCDGetHeadPos(agg_llists, lid, &pos));
581: while (pos) {
582: PetscInt gid1;
584: PetscCall(PetscCDIntNdGetID(pos, &gid1));
585: PetscCall(PetscCDGetNextPos(agg_llists, lid, &pos));
587: if (gid1 >= my0 && gid1 < Iend) flid = gid1 - my0;
588: else {
589: PetscCall(PetscHMapIGet(fgid_flid, gid1, &flid));
590: PetscCheck(flid >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cannot find gid1 in table");
591: }
592: /* copy in B_i matrix - column-oriented */
593: data = &data_in[flid * bs];
594: for (ii = 0; ii < bs; ii++) {
595: for (jj = 0; jj < N; jj++) {
596: PetscReal d = data[jj * data_stride + ii];
598: qqc[jj * Mdata + aggID * bs + ii] = d;
599: }
600: }
601: /* set fine IDs */
602: for (kk = 0; kk < bs; kk++) fids[aggID * bs + kk] = flid_fgid[flid] * bs + kk;
603: aggID++;
604: }
606: /* pad with zeros */
607: for (ii = asz * bs; ii < Mdata; ii++) {
608: for (jj = 0; jj < N; jj++, kk++) qqc[jj * Mdata + ii] = .0;
609: }
611: /* QR */
612: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
613: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&Mdata, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO));
614: PetscCall(PetscFPTrapPop());
615: PetscCheck(INFO == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xGEQRF error");
616: /* get R - column-oriented - output B_{i+1} */
617: {
618: PetscReal *data = &out_data[clid * nSAvec];
620: for (jj = 0; jj < nSAvec; jj++) {
621: for (ii = 0; ii < nSAvec; ii++) {
622: 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);
623: if (ii <= jj) data[jj * out_data_stride + ii] = PetscRealPart(qqc[jj * Mdata + ii]);
624: else data[jj * out_data_stride + ii] = 0.;
625: }
626: }
627: }
629: /* get Q - row-oriented */
630: PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&Mdata, &N, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO));
631: PetscCheck(INFO == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xORGQR error arg %" PetscBLASInt_FMT, -INFO);
633: for (ii = 0; ii < M; ii++) {
634: for (jj = 0; jj < N; jj++) qqr[N * ii + jj] = qqc[jj * Mdata + ii];
635: }
637: /* add diagonal block of P0 */
638: for (kk = 0; kk < N; kk++) cids[kk] = N * cgid + kk; /* global col IDs in P0 */
639: PetscCall(MatSetValues(a_Prol, M, fids, N, cids, qqr, INSERT_VALUES));
640: PetscCall(PetscFree5(qqc, qqr, TAU, WORK, fids));
641: clid++;
642: } /* coarse agg */
643: } /* for all fine nodes */
644: PetscCall(MatAssemblyBegin(a_Prol, MAT_FINAL_ASSEMBLY));
645: PetscCall(MatAssemblyEnd(a_Prol, MAT_FINAL_ASSEMBLY));
646: PetscCall(PetscHMapIDestroy(&fgid_flid));
647: PetscFunctionReturn(PETSC_SUCCESS);
648: }
650: static PetscErrorCode PCView_GAMG_AGG(PC pc, PetscViewer viewer)
651: {
652: PC_MG *mg = (PC_MG *)pc->data;
653: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
654: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
656: PetscFunctionBegin;
657: PetscCall(PetscViewerASCIIPrintf(viewer, " AGG specific options\n"));
658: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of levels of aggressive coarsening %" PetscInt_FMT "\n", pc_gamg_agg->aggressive_coarsening_levels));
659: if (pc_gamg_agg->aggressive_coarsening_levels > 0) {
660: PetscCall(PetscViewerASCIIPrintf(viewer, " %s aggressive coarsening\n", !pc_gamg_agg->use_aggressive_square_graph ? "MIS-k" : "Square graph"));
661: 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));
662: }
663: PetscCall(PetscViewerASCIIPushTab(viewer));
664: PetscCall(PetscViewerASCIIPushTab(viewer));
665: PetscCall(PetscViewerASCIIPushTab(viewer));
666: PetscCall(PetscViewerASCIIPushTab(viewer));
667: if (pc_gamg_agg->crs) PetscCall(MatCoarsenView(pc_gamg_agg->crs, viewer));
668: else PetscCall(PetscViewerASCIIPrintf(viewer, "Coarsening algorithm not yet selected\n"));
669: PetscCall(PetscViewerASCIIPopTab(viewer));
670: PetscCall(PetscViewerASCIIPopTab(viewer));
671: PetscCall(PetscViewerASCIIPopTab(viewer));
672: PetscCall(PetscViewerASCIIPopTab(viewer));
673: PetscCall(PetscViewerASCIIPrintf(viewer, " Number smoothing steps to construct prolongation %" PetscInt_FMT "\n", pc_gamg_agg->nsmooths));
674: PetscFunctionReturn(PETSC_SUCCESS);
675: }
677: static PetscErrorCode PCGAMGCreateGraph_AGG(PC pc, Mat Amat, Mat *a_Gmat)
678: {
679: PC_MG *mg = (PC_MG *)pc->data;
680: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
681: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
682: const PetscReal vfilter = pc_gamg->threshold[pc_gamg->current_level];
683: PetscBool ishem, ismis;
684: const char *prefix;
685: MatInfo info0, info1;
686: PetscInt bs;
688: PetscFunctionBegin;
689: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
690: /* 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 */
691: /* MATCOARSENHEM requires numerical weights for edges so ensure they are computed */
692: PetscCall(MatCoarsenDestroy(&pc_gamg_agg->crs));
693: PetscCall(MatCoarsenCreate(PetscObjectComm((PetscObject)pc), &pc_gamg_agg->crs));
694: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix));
695: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc_gamg_agg->crs, prefix));
696: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)pc_gamg_agg->crs, "pc_gamg_"));
697: PetscCall(MatCoarsenSetFromOptions(pc_gamg_agg->crs));
698: PetscCall(MatGetBlockSize(Amat, &bs));
699: // check for valid indices wrt bs
700: for (int ii = 0; ii < pc_gamg_agg->crs->strength_index_size; ii++) {
701: 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",
702: pc_gamg_agg->crs->strength_index[ii], bs);
703: }
704: PetscCall(PetscObjectTypeCompare((PetscObject)pc_gamg_agg->crs, MATCOARSENHEM, &ishem));
705: if (ishem) {
706: 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));
707: pc_gamg_agg->aggressive_coarsening_levels = 0; // aggressive and HEM does not make sense
708: PetscCall(MatCoarsenSetMaximumIterations(pc_gamg_agg->crs, pc_gamg_agg->crs->max_it)); // for code coverage
709: PetscCall(MatCoarsenSetThreshold(pc_gamg_agg->crs, vfilter)); // for code coverage
710: } else {
711: PetscCall(PetscObjectTypeCompare((PetscObject)pc_gamg_agg->crs, MATCOARSENMIS, &ismis));
712: if (ismis && pc_gamg_agg->aggressive_coarsening_levels && !pc_gamg_agg->use_aggressive_square_graph) {
713: PetscCall(PetscInfo(pc, "MIS and aggressive coarsening and no square graph: force square graph\n"));
714: pc_gamg_agg->use_aggressive_square_graph = PETSC_TRUE;
715: }
716: }
717: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
718: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_GRAPH], 0, 0, 0, 0));
719: PetscCall(MatGetInfo(Amat, MAT_LOCAL, &info0)); /* global reduction */
721: if (ishem || pc_gamg_agg->use_low_mem_filter) {
722: 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));
723: } else {
724: // make scalar graph, symmetrize if not known to be symmetric, scale, but do not filter (expensive)
725: 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));
726: if (vfilter >= 0) {
727: PetscInt Istart, Iend, ncols, nnz0, nnz1, NN, MM, nloc;
728: Mat tGmat, Gmat = *a_Gmat;
729: MPI_Comm comm;
730: const PetscScalar *vals;
731: const PetscInt *idx;
732: PetscInt *d_nnz, *o_nnz, kk, *garray = NULL, *AJ, maxcols = 0;
733: MatScalar *AA; // this is checked in graph
734: PetscBool isseqaij;
735: Mat a, b, c;
736: MatType jtype;
738: PetscCall(PetscObjectGetComm((PetscObject)Gmat, &comm));
739: PetscCall(PetscObjectBaseTypeCompare((PetscObject)Gmat, MATSEQAIJ, &isseqaij));
740: PetscCall(MatGetType(Gmat, &jtype));
741: PetscCall(MatCreate(comm, &tGmat));
742: PetscCall(MatSetType(tGmat, jtype));
744: /* TODO GPU: this can be called when filter = 0 -> Probably provide MatAIJThresholdCompress that compresses the entries below a threshold?
745: Also, if the matrix is symmetric, can we skip this
746: operation? It can be very expensive on large matrices. */
748: // global sizes
749: PetscCall(MatGetSize(Gmat, &MM, &NN));
750: PetscCall(MatGetOwnershipRange(Gmat, &Istart, &Iend));
751: nloc = Iend - Istart;
752: PetscCall(PetscMalloc2(nloc, &d_nnz, nloc, &o_nnz));
753: if (isseqaij) {
754: a = Gmat;
755: b = NULL;
756: } else {
757: Mat_MPIAIJ *d = (Mat_MPIAIJ *)Gmat->data;
759: a = d->A;
760: b = d->B;
761: garray = d->garray;
762: }
763: /* Determine upper bound on non-zeros needed in new filtered matrix */
764: for (PetscInt row = 0; row < nloc; row++) {
765: PetscCall(MatGetRow(a, row, &ncols, NULL, NULL));
766: d_nnz[row] = ncols;
767: if (ncols > maxcols) maxcols = ncols;
768: PetscCall(MatRestoreRow(a, row, &ncols, NULL, NULL));
769: }
770: if (b) {
771: for (PetscInt row = 0; row < nloc; row++) {
772: PetscCall(MatGetRow(b, row, &ncols, NULL, NULL));
773: o_nnz[row] = ncols;
774: if (ncols > maxcols) maxcols = ncols;
775: PetscCall(MatRestoreRow(b, row, &ncols, NULL, NULL));
776: }
777: }
778: PetscCall(MatSetSizes(tGmat, nloc, nloc, MM, MM));
779: PetscCall(MatSetBlockSizes(tGmat, 1, 1));
780: PetscCall(MatSeqAIJSetPreallocation(tGmat, 0, d_nnz));
781: PetscCall(MatMPIAIJSetPreallocation(tGmat, 0, d_nnz, 0, o_nnz));
782: PetscCall(MatSetOption(tGmat, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
783: PetscCall(PetscFree2(d_nnz, o_nnz));
784: PetscCall(PetscMalloc2(maxcols, &AA, maxcols, &AJ));
785: nnz0 = nnz1 = 0;
786: for (c = a, kk = 0; c && kk < 2; c = b, kk++) {
787: for (PetscInt row = 0, grow = Istart, ncol_row, jj; row < nloc; row++, grow++) {
788: PetscCall(MatGetRow(c, row, &ncols, &idx, &vals));
789: for (ncol_row = jj = 0; jj < ncols; jj++, nnz0++) {
790: PetscScalar sv = PetscAbs(PetscRealPart(vals[jj]));
791: if (PetscRealPart(sv) > vfilter) {
792: PetscInt cid = idx[jj] + Istart; //diag
794: nnz1++;
795: if (c != a) cid = garray[idx[jj]];
796: AA[ncol_row] = vals[jj];
797: AJ[ncol_row] = cid;
798: ncol_row++;
799: }
800: }
801: PetscCall(MatRestoreRow(c, row, &ncols, &idx, &vals));
802: PetscCall(MatSetValues(tGmat, 1, &grow, ncol_row, AJ, AA, INSERT_VALUES));
803: }
804: }
805: PetscCall(PetscFree2(AA, AJ));
806: PetscCall(MatAssemblyBegin(tGmat, MAT_FINAL_ASSEMBLY));
807: PetscCall(MatAssemblyEnd(tGmat, MAT_FINAL_ASSEMBLY));
808: PetscCall(MatPropagateSymmetryOptions(Gmat, tGmat)); /* Normal Mat options are not relevant ? */
809: 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));
810: PetscCall(MatViewFromOptions(tGmat, NULL, "-mat_filter_graph_view"));
811: PetscCall(MatDestroy(&Gmat));
812: *a_Gmat = tGmat;
813: }
814: }
816: PetscCall(MatGetInfo(*a_Gmat, MAT_LOCAL, &info1)); /* global reduction */
817: 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));
818: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_GRAPH], 0, 0, 0, 0));
819: PetscFunctionReturn(PETSC_SUCCESS);
820: }
822: typedef PetscInt NState;
823: static const NState NOT_DONE = -2;
824: static const NState DELETED = -1;
825: static const NState REMOVED = -3;
826: #define IS_SELECTED(s) (s != DELETED && s != NOT_DONE && s != REMOVED)
828: /*
829: fixAggregatesWithSquare - greedy grab of with G1 (unsquared graph) -- AIJ specific -- change to fixAggregatesWithSquare -- TODD
830: - AGG-MG specific: clears singletons out of 'selected_2'
832: Input Parameter:
833: . Gmat_2 - global matrix of squared graph (data not defined)
834: . Gmat_1 - base graph to grab with base graph
835: Input/Output Parameter:
836: . aggs_2 - linked list of aggs with gids)
837: */
838: static PetscErrorCode fixAggregatesWithSquare(PC pc, Mat Gmat_2, Mat Gmat_1, PetscCoarsenData *aggs_2)
839: {
840: PetscBool isMPI;
841: Mat_SeqAIJ *matA_1, *matB_1 = NULL;
842: MPI_Comm comm;
843: PetscInt lid, *ii, *idx, ix, Iend, my0, kk, n, j;
844: Mat_MPIAIJ *mpimat_2 = NULL, *mpimat_1 = NULL;
845: const PetscInt nloc = Gmat_2->rmap->n;
846: PetscScalar *cpcol_1_state, *cpcol_2_state, *cpcol_2_par_orig, *lid_parent_gid;
847: PetscInt *lid_cprowID_1 = NULL;
848: NState *lid_state;
849: Vec ghost_par_orig2;
850: PetscMPIInt rank;
852: PetscFunctionBegin;
853: PetscCall(PetscObjectGetComm((PetscObject)Gmat_2, &comm));
854: PetscCallMPI(MPI_Comm_rank(comm, &rank));
855: PetscCall(MatGetOwnershipRange(Gmat_1, &my0, &Iend));
857: /* get submatrices */
858: PetscCall(PetscStrbeginswith(((PetscObject)Gmat_1)->type_name, MATMPIAIJ, &isMPI));
859: PetscCall(PetscInfo(pc, "isMPI = %s\n", isMPI ? "yes" : "no"));
860: PetscCall(PetscMalloc3(nloc, &lid_state, nloc, &lid_parent_gid, nloc, &lid_cprowID_1));
861: for (lid = 0; lid < nloc; lid++) lid_cprowID_1[lid] = -1;
862: if (isMPI) {
863: /* grab matrix objects */
864: mpimat_2 = (Mat_MPIAIJ *)Gmat_2->data;
865: mpimat_1 = (Mat_MPIAIJ *)Gmat_1->data;
866: matA_1 = (Mat_SeqAIJ *)mpimat_1->A->data;
867: matB_1 = (Mat_SeqAIJ *)mpimat_1->B->data;
869: /* force compressed row storage for B matrix in AuxMat */
870: PetscCall(MatCheckCompressedRow(mpimat_1->B, matB_1->nonzerorowcnt, &matB_1->compressedrow, matB_1->i, Gmat_1->rmap->n, -1.0));
871: for (ix = 0; ix < matB_1->compressedrow.nrows; ix++) {
872: PetscInt lid = matB_1->compressedrow.rindex[ix];
874: PetscCheck(lid <= nloc && lid >= -1, PETSC_COMM_SELF, PETSC_ERR_USER, "lid %" PetscInt_FMT " out of range. nloc = %" PetscInt_FMT, lid, nloc);
875: if (lid != -1) lid_cprowID_1[lid] = ix;
876: }
877: } else {
878: PetscBool isAIJ;
880: PetscCall(PetscStrbeginswith(((PetscObject)Gmat_1)->type_name, MATSEQAIJ, &isAIJ));
881: PetscCheck(isAIJ, PETSC_COMM_SELF, PETSC_ERR_USER, "Require AIJ matrix.");
882: matA_1 = (Mat_SeqAIJ *)Gmat_1->data;
883: }
884: if (nloc > 0) PetscCheck(!matB_1 || matB_1->compressedrow.use, PETSC_COMM_SELF, PETSC_ERR_PLIB, "matB_1 && !matB_1->compressedrow.use: PETSc bug???");
885: /* get state of locals and selected gid for deleted */
886: for (lid = 0; lid < nloc; lid++) {
887: lid_parent_gid[lid] = -1.0;
888: lid_state[lid] = DELETED;
889: }
891: /* set lid_state */
892: for (lid = 0; lid < nloc; lid++) {
893: PetscCDIntNd *pos;
895: PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
896: if (pos) {
897: PetscInt gid1;
899: PetscCall(PetscCDIntNdGetID(pos, &gid1));
900: PetscCheck(gid1 == lid + my0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "gid1 %" PetscInt_FMT " != lid %" PetscInt_FMT " + my0 %" PetscInt_FMT, gid1, lid, my0);
901: lid_state[lid] = gid1;
902: }
903: }
905: /* map local to selected local, DELETED means a ghost owns it */
906: for (lid = 0; lid < nloc; lid++) {
907: NState state = lid_state[lid];
909: if (IS_SELECTED(state)) {
910: PetscCDIntNd *pos;
912: PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
913: while (pos) {
914: PetscInt gid1;
916: PetscCall(PetscCDIntNdGetID(pos, &gid1));
917: PetscCall(PetscCDGetNextPos(aggs_2, lid, &pos));
918: if (gid1 >= my0 && gid1 < Iend) lid_parent_gid[gid1 - my0] = (PetscScalar)(lid + my0);
919: }
920: }
921: }
922: /* get 'cpcol_1/2_state' & cpcol_2_par_orig - uses mpimat_1/2->lvec for temp space */
923: if (isMPI) {
924: Vec tempVec;
926: /* get 'cpcol_1_state' */
927: PetscCall(MatCreateVecs(Gmat_1, &tempVec, NULL));
928: for (kk = 0, j = my0; kk < nloc; kk++, j++) {
929: PetscScalar v = (PetscScalar)lid_state[kk];
931: PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
932: }
933: PetscCall(VecAssemblyBegin(tempVec));
934: PetscCall(VecAssemblyEnd(tempVec));
935: PetscCall(VecScatterBegin(mpimat_1->Mvctx, tempVec, mpimat_1->lvec, INSERT_VALUES, SCATTER_FORWARD));
936: PetscCall(VecScatterEnd(mpimat_1->Mvctx, tempVec, mpimat_1->lvec, INSERT_VALUES, SCATTER_FORWARD));
937: PetscCall(VecGetArray(mpimat_1->lvec, &cpcol_1_state));
938: /* get 'cpcol_2_state' */
939: PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, mpimat_2->lvec, INSERT_VALUES, SCATTER_FORWARD));
940: PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, mpimat_2->lvec, INSERT_VALUES, SCATTER_FORWARD));
941: PetscCall(VecGetArray(mpimat_2->lvec, &cpcol_2_state));
942: /* get 'cpcol_2_par_orig' */
943: for (kk = 0, j = my0; kk < nloc; kk++, j++) {
944: PetscScalar v = lid_parent_gid[kk];
946: PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
947: }
948: PetscCall(VecAssemblyBegin(tempVec));
949: PetscCall(VecAssemblyEnd(tempVec));
950: PetscCall(VecDuplicate(mpimat_2->lvec, &ghost_par_orig2));
951: PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghost_par_orig2, INSERT_VALUES, SCATTER_FORWARD));
952: PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghost_par_orig2, INSERT_VALUES, SCATTER_FORWARD));
953: PetscCall(VecGetArray(ghost_par_orig2, &cpcol_2_par_orig));
955: PetscCall(VecDestroy(&tempVec));
956: } /* ismpi */
957: for (lid = 0; lid < nloc; lid++) {
958: NState state = lid_state[lid];
960: if (IS_SELECTED(state)) {
961: /* steal locals */
962: ii = matA_1->i;
963: n = ii[lid + 1] - ii[lid];
964: idx = matA_1->j + ii[lid];
965: for (j = 0; j < n; j++) {
966: PetscInt lidj = idx[j], sgid;
967: NState statej = lid_state[lidj];
969: if (statej == DELETED && (sgid = (PetscInt)PetscRealPart(lid_parent_gid[lidj])) != lid + my0) { /* steal local */
970: lid_parent_gid[lidj] = (PetscScalar)(lid + my0); /* send this if sgid is not local */
971: if (sgid >= my0 && sgid < Iend) { /* I'm stealing this local from a local sgid */
972: PetscInt hav = 0, slid = sgid - my0, gidj = lidj + my0;
973: PetscCDIntNd *pos, *last = NULL;
975: /* looking for local from local so id_llist_2 works */
976: PetscCall(PetscCDGetHeadPos(aggs_2, slid, &pos));
977: while (pos) {
978: PetscInt gid;
980: PetscCall(PetscCDIntNdGetID(pos, &gid));
981: if (gid == gidj) {
982: PetscCheck(last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "last cannot be null");
983: PetscCall(PetscCDRemoveNextNode(aggs_2, slid, last));
984: PetscCall(PetscCDAppendNode(aggs_2, lid, pos));
985: hav = 1;
986: break;
987: } else last = pos;
988: PetscCall(PetscCDGetNextPos(aggs_2, slid, &pos));
989: }
990: if (hav != 1) {
991: PetscCheck(hav, PETSC_COMM_SELF, PETSC_ERR_PLIB, "failed to find adj in 'selected' lists - structurally unsymmetric matrix");
992: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "found node %" PetscInt_FMT " times???", hav);
993: }
994: } else { /* I'm stealing this local, owned by a ghost */
995: 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.",
996: ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "");
997: PetscCall(PetscCDAppendID(aggs_2, lid, lidj + my0));
998: }
999: }
1000: } /* local neighbors */
1001: } else if (state == DELETED /* && lid_cprowID_1 */) {
1002: PetscInt sgidold = (PetscInt)PetscRealPart(lid_parent_gid[lid]);
1004: /* see if I have a selected ghost neighbor that will steal me */
1005: if ((ix = lid_cprowID_1[lid]) != -1) {
1006: ii = matB_1->compressedrow.i;
1007: n = ii[ix + 1] - ii[ix];
1008: idx = matB_1->j + ii[ix];
1009: for (j = 0; j < n; j++) {
1010: PetscInt cpid = idx[j];
1011: NState statej = (NState)PetscRealPart(cpcol_1_state[cpid]);
1013: if (IS_SELECTED(statej) && sgidold != statej) { /* ghost will steal this, remove from my list */
1014: lid_parent_gid[lid] = (PetscScalar)statej; /* send who selected */
1015: if (sgidold >= my0 && sgidold < Iend) { /* this was mine */
1016: PetscInt hav = 0, oldslidj = sgidold - my0;
1017: PetscCDIntNd *pos, *last = NULL;
1019: /* remove from 'oldslidj' list */
1020: PetscCall(PetscCDGetHeadPos(aggs_2, oldslidj, &pos));
1021: while (pos) {
1022: PetscInt gid;
1024: PetscCall(PetscCDIntNdGetID(pos, &gid));
1025: if (lid + my0 == gid) {
1026: /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove lid from oldslidj list *\/ */
1027: PetscCheck(last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "last cannot be null");
1028: PetscCall(PetscCDRemoveNextNode(aggs_2, oldslidj, last));
1029: /* ghost (PetscScalar)statej will add this later */
1030: hav = 1;
1031: break;
1032: } else last = pos;
1033: PetscCall(PetscCDGetNextPos(aggs_2, oldslidj, &pos));
1034: }
1035: if (hav != 1) {
1036: PetscCheck(hav, PETSC_COMM_SELF, PETSC_ERR_PLIB, "failed to find (hav=%" PetscInt_FMT ") adj in 'selected' lists - structurally unsymmetric matrix", hav);
1037: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "found node %" PetscInt_FMT " times???", hav);
1038: }
1039: } else {
1040: /* TODO: ghosts remove this later */
1041: }
1042: }
1043: }
1044: }
1045: } /* selected/deleted */
1046: } /* node loop */
1048: if (isMPI) {
1049: PetscScalar *cpcol_2_parent, *cpcol_2_gid;
1050: Vec tempVec, ghostgids2, ghostparents2;
1051: PetscInt cpid, nghost_2;
1052: PetscHMapI gid_cpid;
1054: PetscCall(VecGetSize(mpimat_2->lvec, &nghost_2));
1055: PetscCall(MatCreateVecs(Gmat_2, &tempVec, NULL));
1057: /* get 'cpcol_2_parent' */
1058: for (kk = 0, j = my0; kk < nloc; kk++, j++) PetscCall(VecSetValues(tempVec, 1, &j, &lid_parent_gid[kk], INSERT_VALUES));
1059: PetscCall(VecAssemblyBegin(tempVec));
1060: PetscCall(VecAssemblyEnd(tempVec));
1061: PetscCall(VecDuplicate(mpimat_2->lvec, &ghostparents2));
1062: PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghostparents2, INSERT_VALUES, SCATTER_FORWARD));
1063: PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghostparents2, INSERT_VALUES, SCATTER_FORWARD));
1064: PetscCall(VecGetArray(ghostparents2, &cpcol_2_parent));
1066: /* get 'cpcol_2_gid' */
1067: for (kk = 0, j = my0; kk < nloc; kk++, j++) {
1068: PetscScalar v = (PetscScalar)j;
1070: PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
1071: }
1072: PetscCall(VecAssemblyBegin(tempVec));
1073: PetscCall(VecAssemblyEnd(tempVec));
1074: PetscCall(VecDuplicate(mpimat_2->lvec, &ghostgids2));
1075: PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghostgids2, INSERT_VALUES, SCATTER_FORWARD));
1076: PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghostgids2, INSERT_VALUES, SCATTER_FORWARD));
1077: PetscCall(VecGetArray(ghostgids2, &cpcol_2_gid));
1078: PetscCall(VecDestroy(&tempVec));
1080: /* look for deleted ghosts and add to table */
1081: PetscCall(PetscHMapICreateWithSize(2 * nghost_2 + 1, &gid_cpid));
1082: for (cpid = 0; cpid < nghost_2; cpid++) {
1083: NState state = (NState)PetscRealPart(cpcol_2_state[cpid]);
1085: if (state == DELETED) {
1086: PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
1087: PetscInt sgid_old = (PetscInt)PetscRealPart(cpcol_2_par_orig[cpid]);
1089: if (sgid_old == -1 && sgid_new != -1) {
1090: PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
1092: PetscCall(PetscHMapISet(gid_cpid, gid, cpid));
1093: }
1094: }
1095: }
1097: /* look for deleted ghosts and see if they moved - remove it */
1098: for (lid = 0; lid < nloc; lid++) {
1099: NState state = lid_state[lid];
1101: if (IS_SELECTED(state)) {
1102: PetscCDIntNd *pos, *last = NULL;
1104: /* look for deleted ghosts and see if they moved */
1105: PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
1106: while (pos) {
1107: PetscInt gid;
1109: PetscCall(PetscCDIntNdGetID(pos, &gid));
1110: if (gid < my0 || gid >= Iend) {
1111: PetscCall(PetscHMapIGet(gid_cpid, gid, &cpid));
1112: if (cpid != -1) {
1113: /* a moved ghost - */
1114: /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove 'flid' from list *\/ */
1115: PetscCall(PetscCDRemoveNextNode(aggs_2, lid, last));
1116: } else last = pos;
1117: } else last = pos;
1119: PetscCall(PetscCDGetNextPos(aggs_2, lid, &pos));
1120: } /* loop over list of deleted */
1121: } /* selected */
1122: }
1123: PetscCall(PetscHMapIDestroy(&gid_cpid));
1125: /* look at ghosts, see if they changed - and it */
1126: for (cpid = 0; cpid < nghost_2; cpid++) {
1127: PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
1129: if (sgid_new >= my0 && sgid_new < Iend) { /* this is mine */
1130: PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
1131: PetscInt slid_new = sgid_new - my0, hav = 0;
1132: PetscCDIntNd *pos;
1134: /* search for this gid to see if I have it */
1135: PetscCall(PetscCDGetHeadPos(aggs_2, slid_new, &pos));
1136: while (pos) {
1137: PetscInt gidj;
1139: PetscCall(PetscCDIntNdGetID(pos, &gidj));
1140: PetscCall(PetscCDGetNextPos(aggs_2, slid_new, &pos));
1142: if (gidj == gid) {
1143: hav = 1;
1144: break;
1145: }
1146: }
1147: if (hav != 1) {
1148: /* insert 'flidj' into head of llist */
1149: PetscCall(PetscCDAppendID(aggs_2, slid_new, gid));
1150: }
1151: }
1152: }
1153: PetscCall(VecRestoreArray(mpimat_1->lvec, &cpcol_1_state));
1154: PetscCall(VecRestoreArray(mpimat_2->lvec, &cpcol_2_state));
1155: PetscCall(VecRestoreArray(ghostparents2, &cpcol_2_parent));
1156: PetscCall(VecRestoreArray(ghostgids2, &cpcol_2_gid));
1157: PetscCall(VecDestroy(&ghostgids2));
1158: PetscCall(VecDestroy(&ghostparents2));
1159: PetscCall(VecDestroy(&ghost_par_orig2));
1160: }
1161: PetscCall(PetscFree3(lid_state, lid_parent_gid, lid_cprowID_1));
1162: PetscFunctionReturn(PETSC_SUCCESS);
1163: }
1165: /*
1166: PCGAMGCoarsen_AGG - supports squaring the graph (deprecated) and new graph for
1167: communication of QR data used with HEM and MISk coarsening
1169: Input Parameter:
1170: . a_pc - this
1172: Input/Output Parameter:
1173: . a_Gmat1 - graph to coarsen (in), graph off processor edges for QR gather scatter (out)
1175: Output Parameter:
1176: . agg_lists - list of aggregates
1178: */
1179: static PetscErrorCode PCGAMGCoarsen_AGG(PC a_pc, Mat *a_Gmat1, PetscCoarsenData **agg_lists)
1180: {
1181: PC_MG *mg = (PC_MG *)a_pc->data;
1182: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1183: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
1184: Mat Gmat2, Gmat1 = *a_Gmat1; /* aggressive graph */
1185: IS perm;
1186: PetscInt Istart, Iend, Ii, nloc, bs, nn;
1187: PetscInt *permute, *degree;
1188: PetscBool *bIndexSet;
1189: PetscReal hashfact;
1190: PetscInt iSwapIndex;
1191: PetscRandom random;
1192: MPI_Comm comm;
1194: PetscFunctionBegin;
1195: PetscCall(PetscObjectGetComm((PetscObject)Gmat1, &comm));
1196: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
1197: PetscCall(MatGetLocalSize(Gmat1, &nn, NULL));
1198: PetscCall(MatGetBlockSize(Gmat1, &bs));
1199: PetscCheck(bs == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "bs %" PetscInt_FMT " must be 1", bs);
1200: nloc = nn / bs;
1201: /* get MIS aggs - randomize */
1202: PetscCall(PetscMalloc2(nloc, &permute, nloc, °ree));
1203: PetscCall(PetscCalloc1(nloc, &bIndexSet));
1204: for (Ii = 0; Ii < nloc; Ii++) permute[Ii] = Ii;
1205: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &random));
1206: PetscCall(MatGetOwnershipRange(Gmat1, &Istart, &Iend));
1207: for (Ii = 0; Ii < nloc; Ii++) {
1208: PetscInt nc;
1210: PetscCall(MatGetRow(Gmat1, Istart + Ii, &nc, NULL, NULL));
1211: degree[Ii] = nc;
1212: PetscCall(MatRestoreRow(Gmat1, Istart + Ii, &nc, NULL, NULL));
1213: }
1214: for (Ii = 0; Ii < nloc; Ii++) {
1215: PetscCall(PetscRandomGetValueReal(random, &hashfact));
1216: iSwapIndex = (PetscInt)(hashfact * nloc) % nloc;
1217: if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) {
1218: PetscInt iTemp = permute[iSwapIndex];
1220: permute[iSwapIndex] = permute[Ii];
1221: permute[Ii] = iTemp;
1222: iTemp = degree[iSwapIndex];
1223: degree[iSwapIndex] = degree[Ii];
1224: degree[Ii] = iTemp;
1225: bIndexSet[iSwapIndex] = PETSC_TRUE;
1226: }
1227: }
1228: // apply minimum degree ordering -- NEW
1229: if (pc_gamg_agg->use_minimum_degree_ordering) PetscCall(PetscSortIntWithArray(nloc, degree, permute));
1230: PetscCall(PetscFree(bIndexSet));
1231: PetscCall(PetscRandomDestroy(&random));
1232: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_USE_POINTER, &perm));
1233: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0));
1234: // square graph
1235: 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));
1236: else Gmat2 = Gmat1;
1237: // switch to old MIS-1 for square graph
1238: if (pc_gamg->current_level < pc_gamg_agg->aggressive_coarsening_levels) {
1239: if (!pc_gamg_agg->use_aggressive_square_graph) PetscCall(MatCoarsenMISKSetDistance(pc_gamg_agg->crs, pc_gamg_agg->aggressive_mis_k)); // hardwire to MIS-2
1240: else PetscCall(MatCoarsenSetType(pc_gamg_agg->crs, MATCOARSENMIS)); // old MIS -- side effect
1241: } else if (pc_gamg_agg->use_aggressive_square_graph && pc_gamg_agg->aggressive_coarsening_levels > 0) { // we reset the MIS
1242: const char *prefix;
1244: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)a_pc, &prefix));
1245: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc_gamg_agg->crs, prefix));
1246: PetscCall(MatCoarsenSetFromOptions(pc_gamg_agg->crs)); // get the default back on non-aggressive levels when square graph switched to old MIS
1247: }
1248: PetscCall(MatCoarsenSetAdjacency(pc_gamg_agg->crs, Gmat2));
1249: PetscCall(MatCoarsenSetStrictAggs(pc_gamg_agg->crs, PETSC_TRUE));
1250: PetscCall(MatCoarsenSetGreedyOrdering(pc_gamg_agg->crs, perm));
1251: PetscCall(MatCoarsenApply(pc_gamg_agg->crs));
1252: PetscCall(MatCoarsenGetData(pc_gamg_agg->crs, agg_lists)); /* output */
1254: PetscCall(ISDestroy(&perm));
1255: PetscCall(PetscFree2(permute, degree));
1256: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0));
1258: if (Gmat2 != Gmat1) { // square graph, we need ghosts for selected
1259: PetscCoarsenData *llist = *agg_lists;
1261: PetscCall(fixAggregatesWithSquare(a_pc, Gmat2, Gmat1, *agg_lists));
1262: PetscCall(MatDestroy(&Gmat1));
1263: *a_Gmat1 = Gmat2; /* output */
1264: PetscCall(PetscCDSetMat(llist, *a_Gmat1)); /* Need a graph with ghosts here */
1265: }
1266: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
1267: PetscFunctionReturn(PETSC_SUCCESS);
1268: }
1270: /*
1271: PCGAMGConstructProlongator_AGG
1273: Input Parameter:
1274: . pc - this
1275: . Amat - matrix on this fine level
1276: . Graph - used to get ghost data for nodes in
1277: . agg_lists - list of aggregates
1278: Output Parameter:
1279: . a_P_out - prolongation operator to the next level
1280: */
1281: static PetscErrorCode PCGAMGConstructProlongator_AGG(PC pc, Mat Amat, PetscCoarsenData *agg_lists, Mat *a_P_out)
1282: {
1283: PC_MG *mg = (PC_MG *)pc->data;
1284: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1285: const PetscInt col_bs = pc_gamg->data_cell_cols;
1286: PetscInt Istart, Iend, nloc, ii, jj, kk, my0, nLocalSelected, bs;
1287: Mat Gmat, Prol;
1288: PetscMPIInt size;
1289: MPI_Comm comm;
1290: PetscReal *data_w_ghost;
1291: PetscInt myCrs0, nbnodes = 0, *flid_fgid;
1292: MatType mtype;
1294: PetscFunctionBegin;
1295: PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
1296: PetscCheck(col_bs >= 1, comm, PETSC_ERR_PLIB, "Column bs cannot be less than 1");
1297: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1298: PetscCallMPI(MPI_Comm_size(comm, &size));
1299: PetscCall(MatGetOwnershipRange(Amat, &Istart, &Iend));
1300: PetscCall(MatGetBlockSize(Amat, &bs));
1301: nloc = (Iend - Istart) / bs;
1302: my0 = Istart / bs;
1303: 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);
1304: PetscCall(PetscCDGetMat(agg_lists, &Gmat)); // get auxiliary matrix for ghost edges for size > 1
1306: /* get 'nLocalSelected' */
1307: for (ii = 0, nLocalSelected = 0; ii < nloc; ii++) {
1308: PetscBool ise;
1310: /* filter out singletons 0 or 1? */
1311: PetscCall(PetscCDIsEmptyAt(agg_lists, ii, &ise));
1312: if (!ise) nLocalSelected++;
1313: }
1315: /* create prolongator, create P matrix */
1316: PetscCall(MatGetType(Amat, &mtype));
1317: PetscCall(MatCreate(comm, &Prol));
1318: PetscCall(MatSetSizes(Prol, nloc * bs, nLocalSelected * col_bs, PETSC_DETERMINE, PETSC_DETERMINE));
1319: PetscCall(MatSetBlockSizes(Prol, bs, col_bs)); // should this be before MatSetSizes?
1320: PetscCall(MatSetType(Prol, mtype));
1321: #if PetscDefined(HAVE_DEVICE)
1322: PetscBool flg;
1323: PetscCall(MatBoundToCPU(Amat, &flg));
1324: PetscCall(MatBindToCPU(Prol, flg));
1325: if (flg) PetscCall(MatSetBindingPropagates(Prol, PETSC_TRUE));
1326: #endif
1327: PetscCall(MatSeqAIJSetPreallocation(Prol, col_bs, NULL));
1328: PetscCall(MatMPIAIJSetPreallocation(Prol, col_bs, NULL, col_bs, NULL));
1330: /* can get all points "removed" */
1331: PetscCall(MatGetSize(Prol, &kk, &ii));
1332: if (!ii) {
1333: PetscCall(PetscInfo(pc, "%s: No selected points on coarse grid\n", ((PetscObject)pc)->prefix));
1334: PetscCall(MatDestroy(&Prol));
1335: *a_P_out = NULL; /* out */
1336: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1337: PetscFunctionReturn(PETSC_SUCCESS);
1338: }
1339: PetscCall(PetscInfo(pc, "%s: New grid %" PetscInt_FMT " nodes\n", ((PetscObject)pc)->prefix, ii / col_bs));
1340: PetscCall(MatGetOwnershipRangeColumn(Prol, &myCrs0, &kk));
1342: 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);
1343: myCrs0 = myCrs0 / col_bs;
1344: 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);
1346: /* create global vector of data in 'data_w_ghost' */
1347: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0));
1348: if (size > 1) { /* get ghost null space data */
1349: PetscReal *tmp_gdata, *tmp_ldata, *tp2;
1351: PetscCall(PetscMalloc1(nloc, &tmp_ldata));
1352: for (jj = 0; jj < col_bs; jj++) {
1353: for (kk = 0; kk < bs; kk++) {
1354: PetscInt ii, stride;
1355: const PetscReal *tp = PetscSafePointerPlusOffset(pc_gamg->data, jj * bs * nloc + kk);
1357: for (ii = 0; ii < nloc; ii++, tp += bs) tmp_ldata[ii] = *tp;
1359: PetscCall(PCGAMGGetDataWithGhosts(Gmat, 1, tmp_ldata, &stride, &tmp_gdata));
1361: if (!jj && !kk) { /* now I know how many total nodes - allocate TODO: move below and do in one 'col_bs' call */
1362: PetscCall(PetscMalloc1(stride * bs * col_bs, &data_w_ghost));
1363: nbnodes = bs * stride;
1364: }
1365: tp2 = PetscSafePointerPlusOffset(data_w_ghost, jj * bs * stride + kk);
1366: for (ii = 0; ii < stride; ii++, tp2 += bs) *tp2 = tmp_gdata[ii];
1367: PetscCall(PetscFree(tmp_gdata));
1368: }
1369: }
1370: PetscCall(PetscFree(tmp_ldata));
1371: } else {
1372: nbnodes = bs * nloc;
1373: data_w_ghost = pc_gamg->data;
1374: }
1376: /* get 'flid_fgid' TODO - move up to get 'stride' and do get null space data above in one step (jj loop) */
1377: if (size > 1) {
1378: PetscReal *fid_glid_loc, *fiddata;
1379: PetscInt stride;
1381: PetscCall(PetscMalloc1(nloc, &fid_glid_loc));
1382: for (kk = 0; kk < nloc; kk++) fid_glid_loc[kk] = (PetscReal)(my0 + kk);
1383: PetscCall(PCGAMGGetDataWithGhosts(Gmat, 1, fid_glid_loc, &stride, &fiddata));
1384: PetscCall(PetscMalloc1(stride, &flid_fgid)); /* copy real data to in */
1385: for (kk = 0; kk < stride; kk++) flid_fgid[kk] = (PetscInt)fiddata[kk];
1386: PetscCall(PetscFree(fiddata));
1388: PetscCheck(stride == nbnodes / bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "stride %" PetscInt_FMT " != nbnodes %" PetscInt_FMT "/bs %" PetscInt_FMT, stride, nbnodes, bs);
1389: PetscCall(PetscFree(fid_glid_loc));
1390: } else {
1391: PetscCall(PetscMalloc1(nloc, &flid_fgid));
1392: for (kk = 0; kk < nloc; kk++) flid_fgid[kk] = my0 + kk;
1393: }
1394: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0));
1395: /* get P0 */
1396: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0));
1397: {
1398: PetscReal *data_out = NULL;
1400: PetscCall(formProl0(agg_lists, bs, col_bs, myCrs0, nbnodes, data_w_ghost, flid_fgid, &data_out, Prol));
1401: PetscCall(PetscFree(pc_gamg->data));
1403: pc_gamg->data = data_out;
1404: pc_gamg->data_cell_rows = col_bs;
1405: pc_gamg->data_sz = col_bs * col_bs * nLocalSelected;
1406: }
1407: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0));
1408: if (size > 1) PetscCall(PetscFree(data_w_ghost));
1409: PetscCall(PetscFree(flid_fgid));
1411: *a_P_out = Prol; /* out */
1412: PetscCall(MatViewFromOptions(Prol, NULL, "-pc_gamg_agg_view_initial_prolongation"));
1414: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1415: PetscFunctionReturn(PETSC_SUCCESS);
1416: }
1418: /*
1419: PCGAMGOptimizeProlongator_AGG - given the initial prolongator optimizes it by smoothed aggregation pc_gamg_agg->nsmooths times
1421: Input Parameter:
1422: . pc - this
1423: . Amat - matrix on this fine level
1424: In/Output Parameter:
1425: . a_P - prolongation operator to the next level
1426: */
1427: static PetscErrorCode PCGAMGOptimizeProlongator_AGG(PC pc, Mat Amat, Mat *a_P)
1428: {
1429: PC_MG *mg = (PC_MG *)pc->data;
1430: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1431: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
1432: PetscInt jj;
1433: Mat Prol = *a_P;
1434: MPI_Comm comm;
1435: KSP eksp;
1436: Vec bb, xx;
1437: PC epc;
1438: PetscReal alpha, emax, emin;
1440: PetscFunctionBegin;
1441: PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
1442: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0));
1444: /* compute maximum singular value of operator to be used in smoother */
1445: if (0 < pc_gamg_agg->nsmooths) {
1446: /* get eigen estimates */
1447: if (pc_gamg->emax > 0) {
1448: emin = pc_gamg->emin;
1449: emax = pc_gamg->emax;
1450: } else {
1451: const char *prefix;
1453: PetscCall(MatCreateVecs(Amat, &bb, NULL));
1454: PetscCall(MatCreateVecs(Amat, &xx, NULL));
1455: PetscCall(KSPSetNoisy_Private(Amat, bb));
1457: PetscCall(KSPCreate(comm, &eksp));
1458: PetscCall(KSPSetNestLevel(eksp, pc->kspnestlevel));
1459: PetscCall(PCGetOptionsPrefix(pc, &prefix));
1460: PetscCall(KSPSetOptionsPrefix(eksp, prefix));
1461: PetscCall(KSPAppendOptionsPrefix(eksp, "pc_gamg_esteig_"));
1462: {
1463: PetscBool isset, sflg;
1465: PetscCall(MatIsSPDKnown(Amat, &isset, &sflg));
1466: if (isset && sflg) PetscCall(KSPSetType(eksp, KSPCG));
1467: }
1468: PetscCall(KSPSetErrorIfNotConverged(eksp, pc->erroriffailure));
1469: PetscCall(KSPSetNormType(eksp, KSP_NORM_NONE));
1471: PetscCall(KSPSetInitialGuessNonzero(eksp, PETSC_FALSE));
1472: PetscCall(KSPSetOperators(eksp, Amat, Amat));
1474: PetscCall(KSPGetPC(eksp, &epc));
1475: PetscCall(PCSetType(epc, PCJACOBI)); /* smoother in smoothed agg. */
1477: 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
1479: PetscCall(KSPSetFromOptions(eksp));
1480: PetscCall(KSPSetComputeSingularValues(eksp, PETSC_TRUE));
1481: PetscCall(KSPSolve(eksp, bb, xx));
1482: PetscCall(KSPCheckSolve(eksp, pc, xx));
1484: PetscCall(KSPComputeExtremeSingularValues(eksp, &emax, &emin));
1485: PetscCall(PetscInfo(pc, "%s: Smooth P0: max eigen=%e min=%e PC=%s\n", ((PetscObject)pc)->prefix, (double)emax, (double)emin, PCJACOBI));
1486: PetscCall(VecDestroy(&xx));
1487: PetscCall(VecDestroy(&bb));
1488: PetscCall(KSPDestroy(&eksp));
1489: }
1490: if (pc_gamg->use_sa_esteig) {
1491: mg->min_eigen_DinvA[pc_gamg->current_level] = emin;
1492: mg->max_eigen_DinvA[pc_gamg->current_level] = emax;
1493: 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));
1494: } else {
1495: mg->min_eigen_DinvA[pc_gamg->current_level] = 0;
1496: mg->max_eigen_DinvA[pc_gamg->current_level] = 0;
1497: }
1498: } else {
1499: mg->min_eigen_DinvA[pc_gamg->current_level] = 0;
1500: mg->max_eigen_DinvA[pc_gamg->current_level] = 0;
1501: }
1503: /* smooth P0 */
1504: if (pc_gamg_agg->nsmooths > 0) {
1505: Vec diag;
1507: /* TODO: Set a PCFailedReason and exit the building of the AMG preconditioner */
1508: PetscCheck(emax != 0.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Computed maximum singular value as zero");
1510: PetscCall(MatCreateVecs(Amat, &diag, NULL));
1511: PetscCall(MatGetDiagonal(Amat, diag)); /* effectively PCJACOBI */
1512: PetscCall(VecReciprocal(diag));
1514: for (jj = 0; jj < pc_gamg_agg->nsmooths; jj++) {
1515: Mat tMat;
1517: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0));
1518: /*
1519: Smooth aggregation on the prolongator
1521: P_{i} := (I - 1.4/emax D^{-1}A) P_i\{i-1}
1522: */
1523: PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0));
1524: PetscCall(MatMatMult(Amat, Prol, MAT_INITIAL_MATRIX, PETSC_CURRENT, &tMat));
1525: PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0));
1526: PetscCall(MatProductClear(tMat));
1527: PetscCall(MatDiagonalScale(tMat, diag, NULL));
1529: /* TODO: Document the 1.4 and don't hardwire it in this routine */
1530: alpha = -1.4 / emax;
1531: PetscCall(MatAYPX(tMat, alpha, Prol, SUBSET_NONZERO_PATTERN));
1532: PetscCall(MatDestroy(&Prol));
1533: Prol = tMat;
1534: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0));
1535: }
1536: PetscCall(VecDestroy(&diag));
1537: }
1538: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0));
1539: PetscCall(MatViewFromOptions(Prol, NULL, "-pc_gamg_agg_view_prolongation"));
1540: *a_P = Prol;
1541: PetscFunctionReturn(PETSC_SUCCESS);
1542: }
1544: /*MC
1545: PCGAMGAGG - Smooth aggregation, {cite}`vanek1996algebraic`, {cite}`vanek2001convergence`, variant of PETSc's algebraic multigrid (`PCGAMG`) preconditioner
1547: Options Database Keys:
1548: + -pc_gamg_agg_nsmooths nsmooth - number of smoothing steps to use with smooth aggregation to construct prolongation
1549: . -pc_gamg_aggressive_coarsening n - number of aggressive coarsening (MIS-2 or square graph) levels from finest.
1550: . -pc_gamg_aggressive_square_graph (true|false) - use square graph ($A^T A$), alternative is MIS-k (k=2), for aggressive coarsening
1551: . -pc_gamg_mis_k_minimum_degree_ordering (true|false) - use minimum degree ordering in greedy MIS algorithm
1552: . -pc_gamg_asm_hem_aggs n - number of HEM aggregation steps for ASM smoother
1553: - -pc_gamg_aggressive_mis_k n - number (k) distance in MIS coarsening (>2 is 'aggressive')
1555: Level: intermediate
1557: Notes:
1558: To obtain good performance for `PCGAMG` for vector valued problems you must
1559: call `MatSetBlockSize()` to indicate the number of degrees of freedom per grid point.
1560: Call `MatSetNearNullSpace()` (or `PCSetCoordinates()` if solving the equations of elasticity) to indicate the near null space of the operator
1562: When `-pc_gamg_aggressive_square_graph` is used, the coarsening is obtained by first squaring the graph and then applying, by default, a
1563: MIS-1 coarsening with `MatCoarsenApply()` on the squared graph.
1565: The many options for `PCMG` and `PCGAMG` such as controlling the smoothers on each level etc. also work for `PCGAMGAGG`
1567: .seealso: `PCGAMG`, [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCCreate()`, `PCSetType()`,
1568: `MatSetBlockSize()`, `PCMGType`, `PCSetCoordinates()`, `MatSetNearNullSpace()`, `PCGAMGSetType()`,
1569: `PCGAMGAGG`, `PCGAMGGEO`, `PCGAMGCLASSICAL`, `PCGAMGSetProcEqLim()`, `PCGAMGSetCoarseEqLim()`, `PCGAMGSetRepartition()`, `PCGAMGRegister()`,
1570: `PCGAMGSetReuseInterpolation()`, `PCGAMGASMSetUseAggs()`, `PCGAMGSetParallelCoarseGridSolve()`, `PCGAMGSetNlevels()`, `PCGAMGSetThreshold()`,
1571: `PCGAMGGetType()`, `PCGAMGSetUseSAEstEig()`
1572: M*/
1573: PetscErrorCode PCCreateGAMG_AGG(PC pc)
1574: {
1575: PC_MG *mg = (PC_MG *)pc->data;
1576: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1577: PC_GAMG_AGG *pc_gamg_agg;
1579: PetscFunctionBegin;
1580: /* create sub context for SA */
1581: PetscCall(PetscNew(&pc_gamg_agg));
1582: pc_gamg->subctx = pc_gamg_agg;
1584: pc_gamg->ops->setfromoptions = PCSetFromOptions_GAMG_AGG;
1585: pc_gamg->ops->destroy = PCDestroy_GAMG_AGG;
1586: /* reset does not do anything; setup not virtual */
1588: /* set internal function pointers */
1589: pc_gamg->ops->creategraph = PCGAMGCreateGraph_AGG;
1590: pc_gamg->ops->coarsen = PCGAMGCoarsen_AGG;
1591: pc_gamg->ops->prolongator = PCGAMGConstructProlongator_AGG;
1592: pc_gamg->ops->optprolongator = PCGAMGOptimizeProlongator_AGG;
1593: pc_gamg->ops->createdefaultdata = PCSetData_AGG;
1594: pc_gamg->ops->view = PCView_GAMG_AGG;
1596: pc_gamg_agg->nsmooths = 1;
1597: pc_gamg_agg->aggressive_coarsening_levels = 1;
1598: pc_gamg_agg->use_aggressive_square_graph = PETSC_TRUE;
1599: pc_gamg_agg->use_minimum_degree_ordering = PETSC_FALSE;
1600: pc_gamg_agg->use_low_mem_filter = PETSC_FALSE;
1601: pc_gamg_agg->aggressive_mis_k = 2;
1602: pc_gamg_agg->graph_symmetrize = PETSC_TRUE;
1604: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", PCGAMGSetNSmooths_AGG));
1605: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", PCGAMGSetAggressiveLevels_AGG));
1606: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveSquareGraph_C", PCGAMGSetAggressiveSquareGraph_AGG));
1607: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetMinDegreeOrdering_C", PCGAMGMISkSetMinDegreeOrdering_AGG));
1608: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetLowMemoryFilter_C", PCGAMGSetLowMemoryFilter_AGG));
1609: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetAggressive_C", PCGAMGMISkSetAggressive_AGG));
1610: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetGraphSymmetrize_C", PCGAMGSetGraphSymmetrize_AGG));
1611: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_AGG));
1612: PetscFunctionReturn(PETSC_SUCCESS);
1613: }