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;
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: MatCoarsen crs;
18: } PC_GAMG_AGG;
20: /*@
21: PCGAMGSetNSmooths - Set number of smoothing steps (1 is typical) used for multigrid on all the levels
23: Logically Collective
25: Input Parameters:
26: + pc - the preconditioner context
27: - n - the number of smooths
29: Options Database Key:
30: . -pc_gamg_agg_nsmooths <nsmooth, default=1> - number of smoothing steps to use with smooth aggregation
32: Level: intermediate
34: .seealso: [](ch_ksp), `PCMG`, `PCGAMG`
35: @*/
36: PetscErrorCode PCGAMGSetNSmooths(PC pc, PetscInt n)
37: {
38: PetscFunctionBegin;
41: PetscTryMethod(pc, "PCGAMGSetNSmooths_C", (PC, PetscInt), (pc, n));
42: PetscFunctionReturn(PETSC_SUCCESS);
43: }
45: static PetscErrorCode PCGAMGSetNSmooths_AGG(PC pc, PetscInt n)
46: {
47: PC_MG *mg = (PC_MG *)pc->data;
48: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
49: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
51: PetscFunctionBegin;
52: pc_gamg_agg->nsmooths = n;
53: PetscFunctionReturn(PETSC_SUCCESS);
54: }
56: /*@
57: PCGAMGSetAggressiveLevels - Use aggressive coarsening on first n levels
59: Logically Collective
61: Input Parameters:
62: + pc - the preconditioner context
63: - n - 0, 1 or more
65: Options Database Key:
66: . -pc_gamg_aggressive_coarsening <n,default = 1> - Number of levels to square the graph on before aggregating it
68: Level: intermediate
70: .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`, `PCGAMGSetLowMemoryFilter()`
71: @*/
72: PetscErrorCode PCGAMGSetAggressiveLevels(PC pc, PetscInt n)
73: {
74: PetscFunctionBegin;
77: PetscTryMethod(pc, "PCGAMGSetAggressiveLevels_C", (PC, PetscInt), (pc, n));
78: PetscFunctionReturn(PETSC_SUCCESS);
79: }
81: /*@
82: PCGAMGMISkSetAggressive - Number (k) distance in MIS coarsening (>2 is 'aggressive')
84: Logically Collective
86: Input Parameters:
87: + pc - the preconditioner context
88: - n - 1 or more (default = 2)
90: Options Database Key:
91: . -pc_gamg_aggressive_mis_k <n,default=2> - Number (k) distance in MIS coarsening (>2 is 'aggressive')
93: Level: intermediate
95: .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`, `PCGAMGSetLowMemoryFilter()`
96: @*/
97: PetscErrorCode PCGAMGMISkSetAggressive(PC pc, PetscInt n)
98: {
99: PetscFunctionBegin;
102: PetscTryMethod(pc, "PCGAMGMISkSetAggressive_C", (PC, PetscInt), (pc, n));
103: PetscFunctionReturn(PETSC_SUCCESS);
104: }
106: /*@
107: PCGAMGSetAggressiveSquareGraph - Use graph square A'A for aggressive coarsening, old method
109: Logically Collective
111: Input Parameters:
112: + pc - the preconditioner context
113: - b - default false - MIS-k is faster
115: Options Database Key:
116: . -pc_gamg_aggressive_square_graph <bool,default=false> - Use square graph (A'A) or MIS-k (k=2) for aggressive coarsening
118: Level: intermediate
120: .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`, `PCGAMGMISkSetAggressive()`, `PCGAMGMISkSetMinDegreeOrdering()`, `PCGAMGSetLowMemoryFilter()`
121: @*/
122: PetscErrorCode PCGAMGSetAggressiveSquareGraph(PC pc, PetscBool b)
123: {
124: PetscFunctionBegin;
127: PetscTryMethod(pc, "PCGAMGSetAggressiveSquareGraph_C", (PC, PetscBool), (pc, b));
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: /*@
132: PCGAMGMISkSetMinDegreeOrdering - Use minimum degree ordering in greedy MIS algorithm
134: Logically Collective
136: Input Parameters:
137: + pc - the preconditioner context
138: - b - default true
140: Options Database Key:
141: . -pc_gamg_mis_k_minimum_degree_ordering <bool,default=true> - Use minimum degree ordering in greedy MIS algorithm
143: Level: intermediate
145: .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`, `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGSetLowMemoryFilter()`
146: @*/
147: PetscErrorCode PCGAMGMISkSetMinDegreeOrdering(PC pc, PetscBool b)
148: {
149: PetscFunctionBegin;
152: PetscTryMethod(pc, "PCGAMGMISkSetMinDegreeOrdering_C", (PC, PetscBool), (pc, b));
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
156: /*@
157: PCGAMGSetLowMemoryFilter - Use low memory graph/matrix filter
159: Logically Collective
161: Input Parameters:
162: + pc - the preconditioner context
163: - b - default false
165: Options Database Key:
166: . -pc_gamg_low_memory_threshold_filter <bool,default=false> - Use low memory graph/matrix filter
168: Level: intermediate
170: .seealso: `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`, `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`
171: @*/
172: PetscErrorCode PCGAMGSetLowMemoryFilter(PC pc, PetscBool b)
173: {
174: PetscFunctionBegin;
177: PetscTryMethod(pc, "PCGAMGSetLowMemoryFilter_C", (PC, PetscBool), (pc, b));
178: PetscFunctionReturn(PETSC_SUCCESS);
179: }
181: static PetscErrorCode PCGAMGSetAggressiveLevels_AGG(PC pc, PetscInt n)
182: {
183: PC_MG *mg = (PC_MG *)pc->data;
184: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
185: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
187: PetscFunctionBegin;
188: pc_gamg_agg->aggressive_coarsening_levels = n;
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: static PetscErrorCode PCGAMGMISkSetAggressive_AGG(PC pc, PetscInt n)
193: {
194: PC_MG *mg = (PC_MG *)pc->data;
195: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
196: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
198: PetscFunctionBegin;
199: pc_gamg_agg->aggressive_mis_k = n;
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: static PetscErrorCode PCGAMGSetAggressiveSquareGraph_AGG(PC pc, PetscBool b)
204: {
205: PC_MG *mg = (PC_MG *)pc->data;
206: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
207: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
209: PetscFunctionBegin;
210: pc_gamg_agg->use_aggressive_square_graph = b;
211: PetscFunctionReturn(PETSC_SUCCESS);
212: }
214: static PetscErrorCode PCGAMGSetLowMemoryFilter_AGG(PC pc, PetscBool b)
215: {
216: PC_MG *mg = (PC_MG *)pc->data;
217: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
218: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
220: PetscFunctionBegin;
221: pc_gamg_agg->use_low_mem_filter = b;
222: PetscFunctionReturn(PETSC_SUCCESS);
223: }
225: static PetscErrorCode PCGAMGMISkSetMinDegreeOrdering_AGG(PC pc, PetscBool b)
226: {
227: PC_MG *mg = (PC_MG *)pc->data;
228: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
229: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
231: PetscFunctionBegin;
232: pc_gamg_agg->use_minimum_degree_ordering = b;
233: PetscFunctionReturn(PETSC_SUCCESS);
234: }
236: static PetscErrorCode PCSetFromOptions_GAMG_AGG(PC pc, PetscOptionItems *PetscOptionsObject)
237: {
238: PC_MG *mg = (PC_MG *)pc->data;
239: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
240: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
241: PetscBool n_aggressive_flg, old_sq_provided = PETSC_FALSE, new_sq_provided = PETSC_FALSE, new_sqr_graph = pc_gamg_agg->use_aggressive_square_graph;
242: PetscInt nsq_graph_old = 0;
244: PetscFunctionBegin;
245: PetscOptionsHeadBegin(PetscOptionsObject, "GAMG-AGG options");
246: PetscCall(PetscOptionsInt("-pc_gamg_agg_nsmooths", "smoothing steps for smoothed aggregation, usually 1", "PCGAMGSetNSmooths", pc_gamg_agg->nsmooths, &pc_gamg_agg->nsmooths, NULL));
247: // aggressive coarsening logic with deprecated -pc_gamg_square_graph
248: 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));
249: if (!n_aggressive_flg)
250: 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));
251: PetscCall(PetscOptionsBool("-pc_gamg_aggressive_square_graph", "Use square graph (A'A) or MIS-k (k=2) for aggressive coarsening", "PCGAMGSetAggressiveSquareGraph", new_sqr_graph, &pc_gamg_agg->use_aggressive_square_graph, &new_sq_provided));
252: if (!new_sq_provided && old_sq_provided) {
253: pc_gamg_agg->aggressive_coarsening_levels = nsq_graph_old; // could be zero
254: pc_gamg_agg->use_aggressive_square_graph = PETSC_TRUE;
255: }
256: if (new_sq_provided && old_sq_provided)
257: 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 %d\n", (int)pc_gamg_agg->aggressive_coarsening_levels));
258: 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));
259: 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));
260: 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));
261: PetscOptionsHeadEnd();
262: PetscFunctionReturn(PETSC_SUCCESS);
263: }
265: static PetscErrorCode PCDestroy_GAMG_AGG(PC pc)
266: {
267: PC_MG *mg = (PC_MG *)pc->data;
268: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
270: PetscFunctionBegin;
271: PetscCall(PetscFree(pc_gamg->subctx));
272: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", NULL));
273: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", NULL));
274: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetAggressive_C", NULL));
275: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetMinDegreeOrdering_C", NULL));
276: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetLowMemoryFilter_C", NULL));
277: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveSquareGraph_C", NULL));
278: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
279: PetscFunctionReturn(PETSC_SUCCESS);
280: }
282: /*
283: PCSetCoordinates_AGG
285: Collective
287: Input Parameter:
288: . pc - the preconditioner context
289: . ndm - dimension of data (used for dof/vertex for Stokes)
290: . a_nloc - number of vertices local
291: . coords - [a_nloc][ndm] - interleaved coordinate data: {x_0, y_0, z_0, x_1, y_1, ...}
292: */
294: static PetscErrorCode PCSetCoordinates_AGG(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords)
295: {
296: PC_MG *mg = (PC_MG *)pc->data;
297: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
298: PetscInt arrsz, kk, ii, jj, nloc, ndatarows, ndf;
299: Mat mat = pc->pmat;
301: PetscFunctionBegin;
304: nloc = a_nloc;
306: /* SA: null space vectors */
307: PetscCall(MatGetBlockSize(mat, &ndf)); /* this does not work for Stokes */
308: if (coords && ndf == 1) pc_gamg->data_cell_cols = 1; /* scalar w/ coords and SA (not needed) */
309: else if (coords) {
310: PetscCheck(ndm <= ndf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "degrees of motion %" PetscInt_FMT " > block size %" PetscInt_FMT, ndm, ndf);
311: pc_gamg->data_cell_cols = (ndm == 2 ? 3 : 6); /* displacement elasticity */
312: 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);
313: } else pc_gamg->data_cell_cols = ndf; /* no data, force SA with constant null space vectors */
314: pc_gamg->data_cell_rows = ndatarows = ndf;
315: 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);
316: arrsz = nloc * pc_gamg->data_cell_rows * pc_gamg->data_cell_cols;
318: if (!pc_gamg->data || (pc_gamg->data_sz != arrsz)) {
319: PetscCall(PetscFree(pc_gamg->data));
320: PetscCall(PetscMalloc1(arrsz + 1, &pc_gamg->data));
321: }
322: /* copy data in - column-oriented */
323: for (kk = 0; kk < nloc; kk++) {
324: const PetscInt M = nloc * pc_gamg->data_cell_rows; /* stride into data */
325: PetscReal *data = &pc_gamg->data[kk * ndatarows]; /* start of cell */
326: if (pc_gamg->data_cell_cols == 1) *data = 1.0;
327: else {
328: /* translational modes */
329: for (ii = 0; ii < ndatarows; ii++) {
330: for (jj = 0; jj < ndatarows; jj++) {
331: if (ii == jj) data[ii * M + jj] = 1.0;
332: else data[ii * M + jj] = 0.0;
333: }
334: }
336: /* rotational modes */
337: if (coords) {
338: if (ndm == 2) {
339: data += 2 * M;
340: data[0] = -coords[2 * kk + 1];
341: data[1] = coords[2 * kk];
342: } else {
343: data += 3 * M;
344: data[0] = 0.0;
345: data[M + 0] = coords[3 * kk + 2];
346: data[2 * M + 0] = -coords[3 * kk + 1];
347: data[1] = -coords[3 * kk + 2];
348: data[M + 1] = 0.0;
349: data[2 * M + 1] = coords[3 * kk];
350: data[2] = coords[3 * kk + 1];
351: data[M + 2] = -coords[3 * kk];
352: data[2 * M + 2] = 0.0;
353: }
354: }
355: }
356: }
357: pc_gamg->data_sz = arrsz;
358: PetscFunctionReturn(PETSC_SUCCESS);
359: }
361: /*
362: PCSetData_AGG - called if data is not set with PCSetCoordinates.
363: Looks in Mat for near null space.
364: Does not work for Stokes
366: Input Parameter:
367: . pc -
368: . a_A - matrix to get (near) null space out of.
369: */
370: static PetscErrorCode PCSetData_AGG(PC pc, Mat a_A)
371: {
372: PC_MG *mg = (PC_MG *)pc->data;
373: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
374: MatNullSpace mnull;
376: PetscFunctionBegin;
377: PetscCall(MatGetNearNullSpace(a_A, &mnull));
378: if (!mnull) {
379: DM dm;
380: PetscCall(PCGetDM(pc, &dm));
381: if (!dm) PetscCall(MatGetDM(a_A, &dm));
382: if (dm) {
383: PetscObject deformation;
384: PetscInt Nf;
386: PetscCall(DMGetNumFields(dm, &Nf));
387: if (Nf) {
388: PetscCall(DMGetField(dm, 0, NULL, &deformation));
389: PetscCall(PetscObjectQuery((PetscObject)deformation, "nearnullspace", (PetscObject *)&mnull));
390: if (!mnull) PetscCall(PetscObjectQuery((PetscObject)deformation, "nullspace", (PetscObject *)&mnull));
391: }
392: }
393: }
395: if (!mnull) {
396: PetscInt bs, NN, MM;
397: PetscCall(MatGetBlockSize(a_A, &bs));
398: PetscCall(MatGetLocalSize(a_A, &MM, &NN));
399: PetscCheck(MM % bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MM %" PetscInt_FMT " must be divisible by bs %" PetscInt_FMT, MM, bs);
400: PetscCall(PCSetCoordinates_AGG(pc, bs, MM / bs, NULL));
401: } else {
402: PetscReal *nullvec;
403: PetscBool has_const;
404: PetscInt i, j, mlocal, nvec, bs;
405: const Vec *vecs;
406: const PetscScalar *v;
408: PetscCall(MatGetLocalSize(a_A, &mlocal, NULL));
409: PetscCall(MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs));
410: for (i = 0; i < nvec; i++) {
411: PetscCall(VecGetLocalSize(vecs[i], &j));
412: PetscCheck(j == mlocal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Attached null space vector size %" PetscInt_FMT " != matrix size %" PetscInt_FMT, j, mlocal);
413: }
414: pc_gamg->data_sz = (nvec + !!has_const) * mlocal;
415: PetscCall(PetscMalloc1((nvec + !!has_const) * mlocal, &nullvec));
416: if (has_const)
417: for (i = 0; i < mlocal; i++) nullvec[i] = 1.0;
418: for (i = 0; i < nvec; i++) {
419: PetscCall(VecGetArrayRead(vecs[i], &v));
420: for (j = 0; j < mlocal; j++) nullvec[(i + !!has_const) * mlocal + j] = PetscRealPart(v[j]);
421: PetscCall(VecRestoreArrayRead(vecs[i], &v));
422: }
423: pc_gamg->data = nullvec;
424: pc_gamg->data_cell_cols = (nvec + !!has_const);
425: PetscCall(MatGetBlockSize(a_A, &bs));
426: pc_gamg->data_cell_rows = bs;
427: }
428: PetscFunctionReturn(PETSC_SUCCESS);
429: }
431: /*
432: formProl0 - collect null space data for each aggregate, do QR, put R in coarse grid data and Q in P_0
434: Input Parameter:
435: . agg_llists - list of arrays with aggregates -- list from selected vertices of aggregate unselected vertices
436: . bs - row block size
437: . nSAvec - column bs of new P
438: . my0crs - global index of start of locals
439: . data_stride - bs*(nloc nodes + ghost nodes) [data_stride][nSAvec]
440: . data_in[data_stride*nSAvec] - local data on fine grid
441: . flid_fgid[data_stride/bs] - make local to global IDs, includes ghosts in 'locals_llist'
443: Output Parameter:
444: . a_data_out - in with fine grid data (w/ghosts), out with coarse grid data
445: . a_Prol - prolongation operator
446: */
447: 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)
448: {
449: PetscInt Istart, my0, Iend, nloc, clid, flid = 0, aggID, kk, jj, ii, mm, nSelected, minsz, nghosts, out_data_stride;
450: MPI_Comm comm;
451: PetscReal *out_data;
452: PetscCDIntNd *pos;
453: PCGAMGHashTable fgid_flid;
455: PetscFunctionBegin;
456: PetscCall(PetscObjectGetComm((PetscObject)a_Prol, &comm));
457: PetscCall(MatGetOwnershipRange(a_Prol, &Istart, &Iend));
458: nloc = (Iend - Istart) / bs;
459: my0 = Istart / bs;
460: 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);
461: Iend /= bs;
462: nghosts = data_stride / bs - nloc;
464: PetscCall(PCGAMGHashTableCreate(2 * nghosts + 1, &fgid_flid));
465: for (kk = 0; kk < nghosts; kk++) PetscCall(PCGAMGHashTableAdd(&fgid_flid, flid_fgid[nloc + kk], nloc + kk));
467: /* count selected -- same as number of cols of P */
468: for (nSelected = mm = 0; mm < nloc; mm++) {
469: PetscBool ise;
470: PetscCall(PetscCDIsEmptyAt(agg_llists, mm, &ise));
471: if (!ise) nSelected++;
472: }
473: PetscCall(MatGetOwnershipRangeColumn(a_Prol, &ii, &jj));
474: PetscCheck((ii / nSAvec) == my0crs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ii %" PetscInt_FMT " /nSAvec %" PetscInt_FMT " != my0crs %" PetscInt_FMT, ii, nSAvec, my0crs);
475: 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);
477: /* aloc space for coarse point data (output) */
478: out_data_stride = nSelected * nSAvec;
480: PetscCall(PetscMalloc1(out_data_stride * nSAvec, &out_data));
481: for (ii = 0; ii < out_data_stride * nSAvec; ii++) out_data[ii] = PETSC_MAX_REAL;
482: *a_data_out = out_data; /* output - stride nSelected*nSAvec */
484: /* find points and set prolongation */
485: minsz = 100;
486: for (mm = clid = 0; mm < nloc; mm++) {
487: PetscCall(PetscCDCountAt(agg_llists, mm, &jj));
488: if (jj > 0) {
489: const PetscInt lid = mm, cgid = my0crs + clid;
490: PetscInt cids[100]; /* max bs */
491: PetscBLASInt asz = jj, M = asz * bs, N = nSAvec, INFO;
492: PetscBLASInt Mdata = M + ((N - M > 0) ? N - M : 0), LDA = Mdata, LWORK = N * bs;
493: PetscScalar *qqc, *qqr, *TAU, *WORK;
494: PetscInt *fids;
495: PetscReal *data;
497: /* count agg */
498: if (asz < minsz) minsz = asz;
500: /* get block */
501: PetscCall(PetscMalloc5(Mdata * N, &qqc, M * N, &qqr, N, &TAU, LWORK, &WORK, M, &fids));
503: aggID = 0;
504: PetscCall(PetscCDGetHeadPos(agg_llists, lid, &pos));
505: while (pos) {
506: PetscInt gid1;
507: PetscCall(PetscCDIntNdGetID(pos, &gid1));
508: PetscCall(PetscCDGetNextPos(agg_llists, lid, &pos));
510: if (gid1 >= my0 && gid1 < Iend) flid = gid1 - my0;
511: else {
512: PetscCall(PCGAMGHashTableFind(&fgid_flid, gid1, &flid));
513: PetscCheck(flid >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cannot find gid1 in table");
514: }
515: /* copy in B_i matrix - column-oriented */
516: data = &data_in[flid * bs];
517: for (ii = 0; ii < bs; ii++) {
518: for (jj = 0; jj < N; jj++) {
519: PetscReal d = data[jj * data_stride + ii];
520: qqc[jj * Mdata + aggID * bs + ii] = d;
521: }
522: }
523: /* set fine IDs */
524: for (kk = 0; kk < bs; kk++) fids[aggID * bs + kk] = flid_fgid[flid] * bs + kk;
525: aggID++;
526: }
528: /* pad with zeros */
529: for (ii = asz * bs; ii < Mdata; ii++) {
530: for (jj = 0; jj < N; jj++, kk++) qqc[jj * Mdata + ii] = .0;
531: }
533: /* QR */
534: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
535: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&Mdata, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO));
536: PetscCall(PetscFPTrapPop());
537: PetscCheck(INFO == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xGEQRF error");
538: /* get R - column-oriented - output B_{i+1} */
539: {
540: PetscReal *data = &out_data[clid * nSAvec];
541: for (jj = 0; jj < nSAvec; jj++) {
542: for (ii = 0; ii < nSAvec; ii++) {
543: 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);
544: if (ii <= jj) data[jj * out_data_stride + ii] = PetscRealPart(qqc[jj * Mdata + ii]);
545: else data[jj * out_data_stride + ii] = 0.;
546: }
547: }
548: }
550: /* get Q - row-oriented */
551: PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&Mdata, &N, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO));
552: PetscCheck(INFO == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xORGQR error arg %" PetscBLASInt_FMT, -INFO);
554: for (ii = 0; ii < M; ii++) {
555: for (jj = 0; jj < N; jj++) qqr[N * ii + jj] = qqc[jj * Mdata + ii];
556: }
558: /* add diagonal block of P0 */
559: for (kk = 0; kk < N; kk++) { cids[kk] = N * cgid + kk; /* global col IDs in P0 */ }
560: PetscCall(MatSetValues(a_Prol, M, fids, N, cids, qqr, INSERT_VALUES));
561: PetscCall(PetscFree5(qqc, qqr, TAU, WORK, fids));
562: clid++;
563: } /* coarse agg */
564: } /* for all fine nodes */
565: PetscCall(MatAssemblyBegin(a_Prol, MAT_FINAL_ASSEMBLY));
566: PetscCall(MatAssemblyEnd(a_Prol, MAT_FINAL_ASSEMBLY));
567: PetscCall(PCGAMGHashTableDestroy(&fgid_flid));
568: PetscFunctionReturn(PETSC_SUCCESS);
569: }
571: static PetscErrorCode PCView_GAMG_AGG(PC pc, PetscViewer viewer)
572: {
573: PC_MG *mg = (PC_MG *)pc->data;
574: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
575: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
577: PetscFunctionBegin;
578: PetscCall(PetscViewerASCIIPrintf(viewer, " AGG specific options\n"));
579: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of levels of aggressive coarsening %d\n", (int)pc_gamg_agg->aggressive_coarsening_levels));
580: if (pc_gamg_agg->aggressive_coarsening_levels > 0) {
581: PetscCall(PetscViewerASCIIPrintf(viewer, " %s aggressive coarsening\n", !pc_gamg_agg->use_aggressive_square_graph ? "MIS-k" : "Square graph"));
582: if (!pc_gamg_agg->use_aggressive_square_graph) PetscCall(PetscViewerASCIIPrintf(viewer, " MIS-%d coarsening on aggressive levels\n", (int)pc_gamg_agg->aggressive_mis_k));
583: }
584: PetscCall(PetscViewerASCIIPrintf(viewer, " Number smoothing steps %d\n", (int)pc_gamg_agg->nsmooths));
585: PetscFunctionReturn(PETSC_SUCCESS);
586: }
588: static PetscErrorCode PCGAMGCreateGraph_AGG(PC pc, Mat Amat, Mat *a_Gmat)
589: {
590: PC_MG *mg = (PC_MG *)pc->data;
591: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
592: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
593: const PetscReal vfilter = pc_gamg->threshold[pc_gamg->current_level];
594: PetscBool ishem, ismis;
595: const char *prefix;
596: MatInfo info0, info1;
597: PetscInt bs;
599: PetscFunctionBegin;
600: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
601: /* 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 */
602: /* MATCOARSENHEM requires numerical weights for edges so ensure they are computed */
603: PetscCall(MatCoarsenCreate(PetscObjectComm((PetscObject)pc), &pc_gamg_agg->crs));
604: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix));
605: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc_gamg_agg->crs, prefix));
606: PetscCall(MatCoarsenSetFromOptions(pc_gamg_agg->crs));
607: PetscCall(MatGetBlockSize(Amat, &bs));
608: // check for valid indices wrt bs
609: for (int ii = 0; ii < pc_gamg_agg->crs->strength_index_size; ii++) {
610: 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 (%d) must be non-negative and < block size (%d), NB, can not use -mat_coarsen_strength_index with -mat_coarsen_strength_index",
611: (int)pc_gamg_agg->crs->strength_index[ii], (int)bs);
612: }
613: PetscCall(PetscObjectTypeCompare((PetscObject)pc_gamg_agg->crs, MATCOARSENHEM, &ishem));
614: if (ishem) {
615: if (pc_gamg_agg->aggressive_coarsening_levels) PetscCall(PetscInfo(pc, "HEM and aggressive coarsening ignored: HEM using %d iterations\n", (int)pc_gamg_agg->crs->max_it));
616: pc_gamg_agg->aggressive_coarsening_levels = 0; // aggressive and HEM does not make sense
617: PetscCall(MatCoarsenSetMaximumIterations(pc_gamg_agg->crs, pc_gamg_agg->crs->max_it)); // for code coverage
618: PetscCall(MatCoarsenSetThreshold(pc_gamg_agg->crs, vfilter)); // for code coverage
619: } else {
620: PetscCall(PetscObjectTypeCompare((PetscObject)pc_gamg_agg->crs, MATCOARSENMIS, &ismis));
621: if (ismis && pc_gamg_agg->aggressive_coarsening_levels && !pc_gamg_agg->use_aggressive_square_graph) {
622: PetscCall(PetscInfo(pc, "MIS and aggressive coarsening and no square graph: force square graph\n"));
623: pc_gamg_agg->use_aggressive_square_graph = PETSC_TRUE;
624: }
625: }
626: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
627: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_GRAPH], 0, 0, 0, 0));
628: PetscCall(MatGetInfo(Amat, MAT_LOCAL, &info0)); /* global reduction */
630: if (ishem || pc_gamg_agg->use_low_mem_filter) {
631: PetscCall(MatCreateGraph(Amat, PETSC_TRUE, (vfilter >= 0 || ishem) ? PETSC_TRUE : PETSC_FALSE, vfilter, pc_gamg_agg->crs->strength_index_size, pc_gamg_agg->crs->strength_index, a_Gmat));
632: } else {
633: // make scalar graph, symetrize if not know to be symmetric, scale, but do not filter (expensive)
634: PetscCall(MatCreateGraph(Amat, PETSC_TRUE, PETSC_TRUE, -1, pc_gamg_agg->crs->strength_index_size, pc_gamg_agg->crs->strength_index, a_Gmat));
635: if (vfilter >= 0) {
636: PetscInt Istart, Iend, ncols, nnz0, nnz1, NN, MM, nloc;
637: Mat tGmat, Gmat = *a_Gmat;
638: MPI_Comm comm;
639: const PetscScalar *vals;
640: const PetscInt *idx;
641: PetscInt *d_nnz, *o_nnz, kk, *garray = NULL, *AJ, maxcols = 0;
642: MatScalar *AA; // this is checked in graph
643: PetscBool isseqaij;
644: Mat a, b, c;
645: MatType jtype;
647: PetscCall(PetscObjectGetComm((PetscObject)Gmat, &comm));
648: PetscCall(PetscObjectBaseTypeCompare((PetscObject)Gmat, MATSEQAIJ, &isseqaij));
649: PetscCall(MatGetType(Gmat, &jtype));
650: PetscCall(MatCreate(comm, &tGmat));
651: PetscCall(MatSetType(tGmat, jtype));
653: /* TODO GPU: this can be called when filter = 0 -> Probably provide MatAIJThresholdCompress that compresses the entries below a threshold?
654: Also, if the matrix is symmetric, can we skip this
655: operation? It can be very expensive on large matrices. */
657: // global sizes
658: PetscCall(MatGetSize(Gmat, &MM, &NN));
659: PetscCall(MatGetOwnershipRange(Gmat, &Istart, &Iend));
660: nloc = Iend - Istart;
661: PetscCall(PetscMalloc2(nloc, &d_nnz, nloc, &o_nnz));
662: if (isseqaij) {
663: a = Gmat;
664: b = NULL;
665: } else {
666: Mat_MPIAIJ *d = (Mat_MPIAIJ *)Gmat->data;
667: a = d->A;
668: b = d->B;
669: garray = d->garray;
670: }
671: /* Determine upper bound on non-zeros needed in new filtered matrix */
672: for (PetscInt row = 0; row < nloc; row++) {
673: PetscCall(MatGetRow(a, row, &ncols, NULL, NULL));
674: d_nnz[row] = ncols;
675: if (ncols > maxcols) maxcols = ncols;
676: PetscCall(MatRestoreRow(a, row, &ncols, NULL, NULL));
677: }
678: if (b) {
679: for (PetscInt row = 0; row < nloc; row++) {
680: PetscCall(MatGetRow(b, row, &ncols, NULL, NULL));
681: o_nnz[row] = ncols;
682: if (ncols > maxcols) maxcols = ncols;
683: PetscCall(MatRestoreRow(b, row, &ncols, NULL, NULL));
684: }
685: }
686: PetscCall(MatSetSizes(tGmat, nloc, nloc, MM, MM));
687: PetscCall(MatSetBlockSizes(tGmat, 1, 1));
688: PetscCall(MatSeqAIJSetPreallocation(tGmat, 0, d_nnz));
689: PetscCall(MatMPIAIJSetPreallocation(tGmat, 0, d_nnz, 0, o_nnz));
690: PetscCall(MatSetOption(tGmat, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
691: PetscCall(PetscFree2(d_nnz, o_nnz));
692: PetscCall(PetscMalloc2(maxcols, &AA, maxcols, &AJ));
693: nnz0 = nnz1 = 0;
694: for (c = a, kk = 0; c && kk < 2; c = b, kk++) {
695: for (PetscInt row = 0, grow = Istart, ncol_row, jj; row < nloc; row++, grow++) {
696: PetscCall(MatGetRow(c, row, &ncols, &idx, &vals));
697: for (ncol_row = jj = 0; jj < ncols; jj++, nnz0++) {
698: PetscScalar sv = PetscAbs(PetscRealPart(vals[jj]));
699: if (PetscRealPart(sv) > vfilter) {
700: PetscInt cid = idx[jj] + Istart; //diag
701: nnz1++;
702: if (c != a) cid = garray[idx[jj]];
703: AA[ncol_row] = vals[jj];
704: AJ[ncol_row] = cid;
705: ncol_row++;
706: }
707: }
708: PetscCall(MatRestoreRow(c, row, &ncols, &idx, &vals));
709: PetscCall(MatSetValues(tGmat, 1, &grow, ncol_row, AJ, AA, INSERT_VALUES));
710: }
711: }
712: PetscCall(PetscFree2(AA, AJ));
713: PetscCall(MatAssemblyBegin(tGmat, MAT_FINAL_ASSEMBLY));
714: PetscCall(MatAssemblyEnd(tGmat, MAT_FINAL_ASSEMBLY));
715: PetscCall(MatPropagateSymmetryOptions(Gmat, tGmat)); /* Normal Mat options are not relevant ? */
716: 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));
717: PetscCall(MatViewFromOptions(tGmat, NULL, "-mat_filter_graph_view"));
718: PetscCall(MatDestroy(&Gmat));
719: *a_Gmat = tGmat;
720: }
721: }
723: PetscCall(MatGetInfo(*a_Gmat, MAT_LOCAL, &info1)); /* global reduction */
724: 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));
725: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_GRAPH], 0, 0, 0, 0));
726: PetscFunctionReturn(PETSC_SUCCESS);
727: }
729: typedef PetscInt NState;
730: static const NState NOT_DONE = -2;
731: static const NState DELETED = -1;
732: static const NState REMOVED = -3;
733: #define IS_SELECTED(s) (s != DELETED && s != NOT_DONE && s != REMOVED)
735: /*
736: fixAggregatesWithSquare - greedy grab of with G1 (unsquared graph) -- AIJ specific -- change to fixAggregatesWithSquare -- TODD
737: - AGG-MG specific: clears singletons out of 'selected_2'
739: Input Parameter:
740: . Gmat_2 - global matrix of squared graph (data not defined)
741: . Gmat_1 - base graph to grab with base graph
742: Input/Output Parameter:
743: . aggs_2 - linked list of aggs with gids)
744: */
745: static PetscErrorCode fixAggregatesWithSquare(PC pc, Mat Gmat_2, Mat Gmat_1, PetscCoarsenData *aggs_2)
746: {
747: PetscBool isMPI;
748: Mat_SeqAIJ *matA_1, *matB_1 = NULL;
749: MPI_Comm comm;
750: PetscInt lid, *ii, *idx, ix, Iend, my0, kk, n, j;
751: Mat_MPIAIJ *mpimat_2 = NULL, *mpimat_1 = NULL;
752: const PetscInt nloc = Gmat_2->rmap->n;
753: PetscScalar *cpcol_1_state, *cpcol_2_state, *cpcol_2_par_orig, *lid_parent_gid;
754: PetscInt *lid_cprowID_1 = NULL;
755: NState *lid_state;
756: Vec ghost_par_orig2;
757: PetscMPIInt rank;
759: PetscFunctionBegin;
760: PetscCall(PetscObjectGetComm((PetscObject)Gmat_2, &comm));
761: PetscCallMPI(MPI_Comm_rank(comm, &rank));
762: PetscCall(MatGetOwnershipRange(Gmat_1, &my0, &Iend));
764: /* get submatrices */
765: PetscCall(PetscStrbeginswith(((PetscObject)Gmat_1)->type_name, MATMPIAIJ, &isMPI));
766: PetscCall(PetscInfo(pc, "isMPI = %s\n", isMPI ? "yes" : "no"));
767: PetscCall(PetscMalloc3(nloc, &lid_state, nloc, &lid_parent_gid, nloc, &lid_cprowID_1));
768: for (lid = 0; lid < nloc; lid++) lid_cprowID_1[lid] = -1;
769: if (isMPI) {
770: /* grab matrix objects */
771: mpimat_2 = (Mat_MPIAIJ *)Gmat_2->data;
772: mpimat_1 = (Mat_MPIAIJ *)Gmat_1->data;
773: matA_1 = (Mat_SeqAIJ *)mpimat_1->A->data;
774: matB_1 = (Mat_SeqAIJ *)mpimat_1->B->data;
776: /* force compressed row storage for B matrix in AuxMat */
777: PetscCall(MatCheckCompressedRow(mpimat_1->B, matB_1->nonzerorowcnt, &matB_1->compressedrow, matB_1->i, Gmat_1->rmap->n, -1.0));
778: for (ix = 0; ix < matB_1->compressedrow.nrows; ix++) {
779: PetscInt lid = matB_1->compressedrow.rindex[ix];
780: PetscCheck(lid <= nloc && lid >= -1, PETSC_COMM_SELF, PETSC_ERR_USER, "lid %d out of range. nloc = %d", (int)lid, (int)nloc);
781: if (lid != -1) lid_cprowID_1[lid] = ix;
782: }
783: } else {
784: PetscBool isAIJ;
785: PetscCall(PetscStrbeginswith(((PetscObject)Gmat_1)->type_name, MATSEQAIJ, &isAIJ));
786: PetscCheck(isAIJ, PETSC_COMM_SELF, PETSC_ERR_USER, "Require AIJ matrix.");
787: matA_1 = (Mat_SeqAIJ *)Gmat_1->data;
788: }
789: if (nloc > 0) { PetscCheck(!matB_1 || matB_1->compressedrow.use, PETSC_COMM_SELF, PETSC_ERR_PLIB, "matB_1 && !matB_1->compressedrow.use: PETSc bug???"); }
790: /* get state of locals and selected gid for deleted */
791: for (lid = 0; lid < nloc; lid++) {
792: lid_parent_gid[lid] = -1.0;
793: lid_state[lid] = DELETED;
794: }
796: /* set lid_state */
797: for (lid = 0; lid < nloc; lid++) {
798: PetscCDIntNd *pos;
799: PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
800: if (pos) {
801: PetscInt gid1;
803: PetscCall(PetscCDIntNdGetID(pos, &gid1));
804: PetscCheck(gid1 == lid + my0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "gid1 %d != lid %d + my0 %d", (int)gid1, (int)lid, (int)my0);
805: lid_state[lid] = gid1;
806: }
807: }
809: /* map local to selected local, DELETED means a ghost owns it */
810: for (lid = kk = 0; lid < nloc; lid++) {
811: NState state = lid_state[lid];
812: if (IS_SELECTED(state)) {
813: PetscCDIntNd *pos;
814: PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
815: while (pos) {
816: PetscInt gid1;
817: PetscCall(PetscCDIntNdGetID(pos, &gid1));
818: PetscCall(PetscCDGetNextPos(aggs_2, lid, &pos));
819: if (gid1 >= my0 && gid1 < Iend) lid_parent_gid[gid1 - my0] = (PetscScalar)(lid + my0);
820: }
821: }
822: }
823: /* get 'cpcol_1/2_state' & cpcol_2_par_orig - uses mpimat_1/2->lvec for temp space */
824: if (isMPI) {
825: Vec tempVec;
826: /* get 'cpcol_1_state' */
827: PetscCall(MatCreateVecs(Gmat_1, &tempVec, NULL));
828: for (kk = 0, j = my0; kk < nloc; kk++, j++) {
829: PetscScalar v = (PetscScalar)lid_state[kk];
830: PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
831: }
832: PetscCall(VecAssemblyBegin(tempVec));
833: PetscCall(VecAssemblyEnd(tempVec));
834: PetscCall(VecScatterBegin(mpimat_1->Mvctx, tempVec, mpimat_1->lvec, INSERT_VALUES, SCATTER_FORWARD));
835: PetscCall(VecScatterEnd(mpimat_1->Mvctx, tempVec, mpimat_1->lvec, INSERT_VALUES, SCATTER_FORWARD));
836: PetscCall(VecGetArray(mpimat_1->lvec, &cpcol_1_state));
837: /* get 'cpcol_2_state' */
838: PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, mpimat_2->lvec, INSERT_VALUES, SCATTER_FORWARD));
839: PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, mpimat_2->lvec, INSERT_VALUES, SCATTER_FORWARD));
840: PetscCall(VecGetArray(mpimat_2->lvec, &cpcol_2_state));
841: /* get 'cpcol_2_par_orig' */
842: for (kk = 0, j = my0; kk < nloc; kk++, j++) {
843: PetscScalar v = (PetscScalar)lid_parent_gid[kk];
844: PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
845: }
846: PetscCall(VecAssemblyBegin(tempVec));
847: PetscCall(VecAssemblyEnd(tempVec));
848: PetscCall(VecDuplicate(mpimat_2->lvec, &ghost_par_orig2));
849: PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghost_par_orig2, INSERT_VALUES, SCATTER_FORWARD));
850: PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghost_par_orig2, INSERT_VALUES, SCATTER_FORWARD));
851: PetscCall(VecGetArray(ghost_par_orig2, &cpcol_2_par_orig));
853: PetscCall(VecDestroy(&tempVec));
854: } /* ismpi */
855: for (lid = 0; lid < nloc; lid++) {
856: NState state = lid_state[lid];
857: if (IS_SELECTED(state)) {
858: /* steal locals */
859: ii = matA_1->i;
860: n = ii[lid + 1] - ii[lid];
861: idx = matA_1->j + ii[lid];
862: for (j = 0; j < n; j++) {
863: PetscInt lidj = idx[j], sgid;
864: NState statej = lid_state[lidj];
865: if (statej == DELETED && (sgid = (PetscInt)PetscRealPart(lid_parent_gid[lidj])) != lid + my0) { /* steal local */
866: lid_parent_gid[lidj] = (PetscScalar)(lid + my0); /* send this if sgid is not local */
867: if (sgid >= my0 && sgid < Iend) { /* I'm stealing this local from a local sgid */
868: PetscInt hav = 0, slid = sgid - my0, gidj = lidj + my0;
869: PetscCDIntNd *pos, *last = NULL;
870: /* looking for local from local so id_llist_2 works */
871: PetscCall(PetscCDGetHeadPos(aggs_2, slid, &pos));
872: while (pos) {
873: PetscInt gid;
874: PetscCall(PetscCDIntNdGetID(pos, &gid));
875: if (gid == gidj) {
876: PetscCheck(last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "last cannot be null");
877: PetscCall(PetscCDRemoveNextNode(aggs_2, slid, last));
878: PetscCall(PetscCDAppendNode(aggs_2, lid, pos));
879: hav = 1;
880: break;
881: } else last = pos;
882: PetscCall(PetscCDGetNextPos(aggs_2, slid, &pos));
883: }
884: if (hav != 1) {
885: PetscCheck(hav, PETSC_COMM_SELF, PETSC_ERR_PLIB, "failed to find adj in 'selected' lists - structurally unsymmetric matrix");
886: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "found node %d times???", (int)hav);
887: }
888: } else { /* I'm stealing this local, owned by a ghost */
889: 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.",
890: ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "");
891: PetscCall(PetscCDAppendID(aggs_2, lid, lidj + my0));
892: }
893: }
894: } /* local neighbors */
895: } else if (state == DELETED /* && lid_cprowID_1 */) {
896: PetscInt sgidold = (PetscInt)PetscRealPart(lid_parent_gid[lid]);
897: /* see if I have a selected ghost neighbor that will steal me */
898: if ((ix = lid_cprowID_1[lid]) != -1) {
899: ii = matB_1->compressedrow.i;
900: n = ii[ix + 1] - ii[ix];
901: idx = matB_1->j + ii[ix];
902: for (j = 0; j < n; j++) {
903: PetscInt cpid = idx[j];
904: NState statej = (NState)PetscRealPart(cpcol_1_state[cpid]);
905: if (IS_SELECTED(statej) && sgidold != (PetscInt)statej) { /* ghost will steal this, remove from my list */
906: lid_parent_gid[lid] = (PetscScalar)statej; /* send who selected */
907: if (sgidold >= my0 && sgidold < Iend) { /* this was mine */
908: PetscInt hav = 0, oldslidj = sgidold - my0;
909: PetscCDIntNd *pos, *last = NULL;
910: /* remove from 'oldslidj' list */
911: PetscCall(PetscCDGetHeadPos(aggs_2, oldslidj, &pos));
912: while (pos) {
913: PetscInt gid;
914: PetscCall(PetscCDIntNdGetID(pos, &gid));
915: if (lid + my0 == gid) {
916: /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove lid from oldslidj list *\/ */
917: PetscCheck(last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "last cannot be null");
918: PetscCall(PetscCDRemoveNextNode(aggs_2, oldslidj, last));
919: /* ghost (PetscScalar)statej will add this later */
920: hav = 1;
921: break;
922: } else last = pos;
923: PetscCall(PetscCDGetNextPos(aggs_2, oldslidj, &pos));
924: }
925: if (hav != 1) {
926: PetscCheck(hav, PETSC_COMM_SELF, PETSC_ERR_PLIB, "failed to find (hav=%d) adj in 'selected' lists - structurally unsymmetric matrix", (int)hav);
927: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "found node %d times???", (int)hav);
928: }
929: } else {
930: /* TODO: ghosts remove this later */
931: }
932: }
933: }
934: }
935: } /* selected/deleted */
936: } /* node loop */
938: if (isMPI) {
939: PetscScalar *cpcol_2_parent, *cpcol_2_gid;
940: Vec tempVec, ghostgids2, ghostparents2;
941: PetscInt cpid, nghost_2;
942: PCGAMGHashTable gid_cpid;
944: PetscCall(VecGetSize(mpimat_2->lvec, &nghost_2));
945: PetscCall(MatCreateVecs(Gmat_2, &tempVec, NULL));
947: /* get 'cpcol_2_parent' */
948: for (kk = 0, j = my0; kk < nloc; kk++, j++) { PetscCall(VecSetValues(tempVec, 1, &j, &lid_parent_gid[kk], INSERT_VALUES)); }
949: PetscCall(VecAssemblyBegin(tempVec));
950: PetscCall(VecAssemblyEnd(tempVec));
951: PetscCall(VecDuplicate(mpimat_2->lvec, &ghostparents2));
952: PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghostparents2, INSERT_VALUES, SCATTER_FORWARD));
953: PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghostparents2, INSERT_VALUES, SCATTER_FORWARD));
954: PetscCall(VecGetArray(ghostparents2, &cpcol_2_parent));
956: /* get 'cpcol_2_gid' */
957: for (kk = 0, j = my0; kk < nloc; kk++, j++) {
958: PetscScalar v = (PetscScalar)j;
959: PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
960: }
961: PetscCall(VecAssemblyBegin(tempVec));
962: PetscCall(VecAssemblyEnd(tempVec));
963: PetscCall(VecDuplicate(mpimat_2->lvec, &ghostgids2));
964: PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghostgids2, INSERT_VALUES, SCATTER_FORWARD));
965: PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghostgids2, INSERT_VALUES, SCATTER_FORWARD));
966: PetscCall(VecGetArray(ghostgids2, &cpcol_2_gid));
967: PetscCall(VecDestroy(&tempVec));
969: /* look for deleted ghosts and add to table */
970: PetscCall(PCGAMGHashTableCreate(2 * nghost_2 + 1, &gid_cpid));
971: for (cpid = 0; cpid < nghost_2; cpid++) {
972: NState state = (NState)PetscRealPart(cpcol_2_state[cpid]);
973: if (state == DELETED) {
974: PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
975: PetscInt sgid_old = (PetscInt)PetscRealPart(cpcol_2_par_orig[cpid]);
976: if (sgid_old == -1 && sgid_new != -1) {
977: PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
978: PetscCall(PCGAMGHashTableAdd(&gid_cpid, gid, cpid));
979: }
980: }
981: }
983: /* look for deleted ghosts and see if they moved - remove it */
984: for (lid = 0; lid < nloc; lid++) {
985: NState state = lid_state[lid];
986: if (IS_SELECTED(state)) {
987: PetscCDIntNd *pos, *last = NULL;
988: /* look for deleted ghosts and see if they moved */
989: PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
990: while (pos) {
991: PetscInt gid;
992: PetscCall(PetscCDIntNdGetID(pos, &gid));
994: if (gid < my0 || gid >= Iend) {
995: PetscCall(PCGAMGHashTableFind(&gid_cpid, gid, &cpid));
996: if (cpid != -1) {
997: /* a moved ghost - */
998: /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove 'flid' from list *\/ */
999: PetscCall(PetscCDRemoveNextNode(aggs_2, lid, last));
1000: } else last = pos;
1001: } else last = pos;
1003: PetscCall(PetscCDGetNextPos(aggs_2, lid, &pos));
1004: } /* loop over list of deleted */
1005: } /* selected */
1006: }
1007: PetscCall(PCGAMGHashTableDestroy(&gid_cpid));
1009: /* look at ghosts, see if they changed - and it */
1010: for (cpid = 0; cpid < nghost_2; cpid++) {
1011: PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
1012: if (sgid_new >= my0 && sgid_new < Iend) { /* this is mine */
1013: PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
1014: PetscInt slid_new = sgid_new - my0, hav = 0;
1015: PetscCDIntNd *pos;
1017: /* search for this gid to see if I have it */
1018: PetscCall(PetscCDGetHeadPos(aggs_2, slid_new, &pos));
1019: while (pos) {
1020: PetscInt gidj;
1021: PetscCall(PetscCDIntNdGetID(pos, &gidj));
1022: PetscCall(PetscCDGetNextPos(aggs_2, slid_new, &pos));
1024: if (gidj == gid) {
1025: hav = 1;
1026: break;
1027: }
1028: }
1029: if (hav != 1) {
1030: /* insert 'flidj' into head of llist */
1031: PetscCall(PetscCDAppendID(aggs_2, slid_new, gid));
1032: }
1033: }
1034: }
1035: PetscCall(VecRestoreArray(mpimat_1->lvec, &cpcol_1_state));
1036: PetscCall(VecRestoreArray(mpimat_2->lvec, &cpcol_2_state));
1037: PetscCall(VecRestoreArray(ghostparents2, &cpcol_2_parent));
1038: PetscCall(VecRestoreArray(ghostgids2, &cpcol_2_gid));
1039: PetscCall(VecDestroy(&ghostgids2));
1040: PetscCall(VecDestroy(&ghostparents2));
1041: PetscCall(VecDestroy(&ghost_par_orig2));
1042: }
1043: PetscCall(PetscFree3(lid_state, lid_parent_gid, lid_cprowID_1));
1044: PetscFunctionReturn(PETSC_SUCCESS);
1045: }
1047: /*
1048: PCGAMGCoarsen_AGG - supports squaring the graph (deprecated) and new graph for
1049: communication of QR data used with HEM and MISk coarsening
1051: Input Parameter:
1052: . a_pc - this
1054: Input/Output Parameter:
1055: . a_Gmat1 - graph to coarsen (in), graph off processor edges for QR gather scatter (out)
1057: Output Parameter:
1058: . agg_lists - list of aggregates
1060: */
1061: static PetscErrorCode PCGAMGCoarsen_AGG(PC a_pc, Mat *a_Gmat1, PetscCoarsenData **agg_lists)
1062: {
1063: PC_MG *mg = (PC_MG *)a_pc->data;
1064: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1065: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
1066: Mat Gmat2, Gmat1 = *a_Gmat1; /* aggressive graph */
1067: IS perm;
1068: PetscInt Istart, Iend, Ii, nloc, bs, nn;
1069: PetscInt *permute, *degree;
1070: PetscBool *bIndexSet;
1071: PetscReal hashfact;
1072: PetscInt iSwapIndex;
1073: PetscRandom random;
1074: MPI_Comm comm;
1076: PetscFunctionBegin;
1077: PetscCall(PetscObjectGetComm((PetscObject)Gmat1, &comm));
1078: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
1079: PetscCall(MatGetLocalSize(Gmat1, &nn, NULL));
1080: PetscCall(MatGetBlockSize(Gmat1, &bs));
1081: PetscCheck(bs == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "bs %" PetscInt_FMT " must be 1", bs);
1082: nloc = nn / bs;
1083: /* get MIS aggs - randomize */
1084: PetscCall(PetscMalloc2(nloc, &permute, nloc, °ree));
1085: PetscCall(PetscCalloc1(nloc, &bIndexSet));
1086: for (Ii = 0; Ii < nloc; Ii++) permute[Ii] = Ii;
1087: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &random));
1088: PetscCall(MatGetOwnershipRange(Gmat1, &Istart, &Iend));
1089: for (Ii = 0; Ii < nloc; Ii++) {
1090: PetscInt nc;
1091: PetscCall(MatGetRow(Gmat1, Istart + Ii, &nc, NULL, NULL));
1092: degree[Ii] = nc;
1093: PetscCall(MatRestoreRow(Gmat1, Istart + Ii, &nc, NULL, NULL));
1094: }
1095: for (Ii = 0; Ii < nloc; Ii++) {
1096: PetscCall(PetscRandomGetValueReal(random, &hashfact));
1097: iSwapIndex = (PetscInt)(hashfact * nloc) % nloc;
1098: if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) {
1099: PetscInt iTemp = permute[iSwapIndex];
1100: permute[iSwapIndex] = permute[Ii];
1101: permute[Ii] = iTemp;
1102: iTemp = degree[iSwapIndex];
1103: degree[iSwapIndex] = degree[Ii];
1104: degree[Ii] = iTemp;
1105: bIndexSet[iSwapIndex] = PETSC_TRUE;
1106: }
1107: }
1108: // apply minimum degree ordering -- NEW
1109: if (pc_gamg_agg->use_minimum_degree_ordering) { PetscCall(PetscSortIntWithArray(nloc, degree, permute)); }
1110: PetscCall(PetscFree(bIndexSet));
1111: PetscCall(PetscRandomDestroy(&random));
1112: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_USE_POINTER, &perm));
1113: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0));
1114: // square graph
1115: if (pc_gamg->current_level < pc_gamg_agg->aggressive_coarsening_levels && pc_gamg_agg->use_aggressive_square_graph) {
1116: PetscCall(PCGAMGSquareGraph_GAMG(a_pc, Gmat1, &Gmat2));
1117: } else Gmat2 = Gmat1;
1118: // switch to old MIS-1 for square graph
1119: if (pc_gamg->current_level < pc_gamg_agg->aggressive_coarsening_levels) {
1120: if (!pc_gamg_agg->use_aggressive_square_graph) PetscCall(MatCoarsenMISKSetDistance(pc_gamg_agg->crs, pc_gamg_agg->aggressive_mis_k)); // hardwire to MIS-2
1121: else PetscCall(MatCoarsenSetType(pc_gamg_agg->crs, MATCOARSENMIS)); // old MIS -- side effect
1122: } else if (pc_gamg_agg->use_aggressive_square_graph && pc_gamg_agg->aggressive_coarsening_levels > 0) { // we reset the MIS
1123: const char *prefix;
1124: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)a_pc, &prefix));
1125: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc_gamg_agg->crs, prefix));
1126: PetscCall(MatCoarsenSetFromOptions(pc_gamg_agg->crs)); // get the default back on non-aggressive levels when square graph switched to old MIS
1127: }
1128: PetscCall(MatCoarsenSetAdjacency(pc_gamg_agg->crs, Gmat2));
1129: PetscCall(MatCoarsenSetStrictAggs(pc_gamg_agg->crs, PETSC_TRUE));
1130: PetscCall(MatCoarsenSetGreedyOrdering(pc_gamg_agg->crs, perm));
1131: PetscCall(MatCoarsenApply(pc_gamg_agg->crs));
1132: PetscCall(MatCoarsenViewFromOptions(pc_gamg_agg->crs, NULL, "-mat_coarsen_view"));
1133: PetscCall(MatCoarsenGetData(pc_gamg_agg->crs, agg_lists)); /* output */
1134: PetscCall(MatCoarsenDestroy(&pc_gamg_agg->crs));
1136: PetscCall(ISDestroy(&perm));
1137: PetscCall(PetscFree2(permute, degree));
1138: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0));
1140: if (Gmat2 != Gmat1) { // square graph, we need ghosts for selected
1141: PetscCoarsenData *llist = *agg_lists;
1142: PetscCall(fixAggregatesWithSquare(a_pc, Gmat2, Gmat1, *agg_lists));
1143: PetscCall(MatDestroy(&Gmat1));
1144: *a_Gmat1 = Gmat2; /* output */
1145: PetscCall(PetscCDSetMat(llist, *a_Gmat1)); /* Need a graph with ghosts here */
1146: }
1147: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
1148: PetscFunctionReturn(PETSC_SUCCESS);
1149: }
1151: /*
1152: PCGAMGProlongator_AGG
1154: Input Parameter:
1155: . pc - this
1156: . Amat - matrix on this fine level
1157: . Graph - used to get ghost data for nodes in
1158: . agg_lists - list of aggregates
1159: Output Parameter:
1160: . a_P_out - prolongation operator to the next level
1161: */
1162: static PetscErrorCode PCGAMGProlongator_AGG(PC pc, Mat Amat, PetscCoarsenData *agg_lists, Mat *a_P_out)
1163: {
1164: PC_MG *mg = (PC_MG *)pc->data;
1165: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1166: const PetscInt col_bs = pc_gamg->data_cell_cols;
1167: PetscInt Istart, Iend, nloc, ii, jj, kk, my0, nLocalSelected, bs;
1168: Mat Gmat, Prol;
1169: PetscMPIInt size;
1170: MPI_Comm comm;
1171: PetscReal *data_w_ghost;
1172: PetscInt myCrs0, nbnodes = 0, *flid_fgid;
1173: MatType mtype;
1175: PetscFunctionBegin;
1176: PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
1177: PetscCheck(col_bs >= 1, comm, PETSC_ERR_PLIB, "Column bs cannot be less than 1");
1178: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1179: PetscCallMPI(MPI_Comm_size(comm, &size));
1180: PetscCall(MatGetOwnershipRange(Amat, &Istart, &Iend));
1181: PetscCall(MatGetBlockSize(Amat, &bs));
1182: nloc = (Iend - Istart) / bs;
1183: my0 = Istart / bs;
1184: 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);
1185: PetscCall(PetscCDGetMat(agg_lists, &Gmat)); // get auxiliary matrix for ghost edges for size > 1
1187: /* get 'nLocalSelected' */
1188: for (ii = 0, nLocalSelected = 0; ii < nloc; ii++) {
1189: PetscBool ise;
1190: /* filter out singletons 0 or 1? */
1191: PetscCall(PetscCDIsEmptyAt(agg_lists, ii, &ise));
1192: if (!ise) nLocalSelected++;
1193: }
1195: /* create prolongator, create P matrix */
1196: PetscCall(MatGetType(Amat, &mtype));
1197: PetscCall(MatCreate(comm, &Prol));
1198: PetscCall(MatSetSizes(Prol, nloc * bs, nLocalSelected * col_bs, PETSC_DETERMINE, PETSC_DETERMINE));
1199: PetscCall(MatSetBlockSizes(Prol, bs, col_bs)); // should this be before MatSetSizes?
1200: PetscCall(MatSetType(Prol, mtype));
1201: #if PetscDefined(HAVE_DEVICE)
1202: PetscBool flg;
1203: PetscCall(MatBoundToCPU(Amat, &flg));
1204: PetscCall(MatBindToCPU(Prol, flg));
1205: if (flg) PetscCall(MatSetBindingPropagates(Prol, PETSC_TRUE));
1206: #endif
1207: PetscCall(MatSeqAIJSetPreallocation(Prol, col_bs, NULL));
1208: PetscCall(MatMPIAIJSetPreallocation(Prol, col_bs, NULL, col_bs, NULL));
1210: /* can get all points "removed" */
1211: PetscCall(MatGetSize(Prol, &kk, &ii));
1212: if (!ii) {
1213: PetscCall(PetscInfo(pc, "%s: No selected points on coarse grid\n", ((PetscObject)pc)->prefix));
1214: PetscCall(MatDestroy(&Prol));
1215: *a_P_out = NULL; /* out */
1216: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1217: PetscFunctionReturn(PETSC_SUCCESS);
1218: }
1219: PetscCall(PetscInfo(pc, "%s: New grid %" PetscInt_FMT " nodes\n", ((PetscObject)pc)->prefix, ii / col_bs));
1220: PetscCall(MatGetOwnershipRangeColumn(Prol, &myCrs0, &kk));
1222: 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);
1223: myCrs0 = myCrs0 / col_bs;
1224: 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);
1226: /* create global vector of data in 'data_w_ghost' */
1227: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0));
1228: if (size > 1) { /* get ghost null space data */
1229: PetscReal *tmp_gdata, *tmp_ldata, *tp2;
1230: PetscCall(PetscMalloc1(nloc, &tmp_ldata));
1231: for (jj = 0; jj < col_bs; jj++) {
1232: for (kk = 0; kk < bs; kk++) {
1233: PetscInt ii, stride;
1234: const PetscReal *tp = PetscSafePointerPlusOffset(pc_gamg->data, jj * bs * nloc + kk);
1235: for (ii = 0; ii < nloc; ii++, tp += bs) tmp_ldata[ii] = *tp;
1237: PetscCall(PCGAMGGetDataWithGhosts(Gmat, 1, tmp_ldata, &stride, &tmp_gdata));
1239: if (!jj && !kk) { /* now I know how many total nodes - allocate TODO: move below and do in one 'col_bs' call */
1240: PetscCall(PetscMalloc1(stride * bs * col_bs, &data_w_ghost));
1241: nbnodes = bs * stride;
1242: }
1243: tp2 = PetscSafePointerPlusOffset(data_w_ghost, jj * bs * stride + kk);
1244: for (ii = 0; ii < stride; ii++, tp2 += bs) *tp2 = tmp_gdata[ii];
1245: PetscCall(PetscFree(tmp_gdata));
1246: }
1247: }
1248: PetscCall(PetscFree(tmp_ldata));
1249: } else {
1250: nbnodes = bs * nloc;
1251: data_w_ghost = (PetscReal *)pc_gamg->data;
1252: }
1254: /* get 'flid_fgid' TODO - move up to get 'stride' and do get null space data above in one step (jj loop) */
1255: if (size > 1) {
1256: PetscReal *fid_glid_loc, *fiddata;
1257: PetscInt stride;
1259: PetscCall(PetscMalloc1(nloc, &fid_glid_loc));
1260: for (kk = 0; kk < nloc; kk++) fid_glid_loc[kk] = (PetscReal)(my0 + kk);
1261: PetscCall(PCGAMGGetDataWithGhosts(Gmat, 1, fid_glid_loc, &stride, &fiddata));
1262: PetscCall(PetscMalloc1(stride, &flid_fgid)); /* copy real data to in */
1263: for (kk = 0; kk < stride; kk++) flid_fgid[kk] = (PetscInt)fiddata[kk];
1264: PetscCall(PetscFree(fiddata));
1266: PetscCheck(stride == nbnodes / bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "stride %" PetscInt_FMT " != nbnodes %" PetscInt_FMT "/bs %" PetscInt_FMT, stride, nbnodes, bs);
1267: PetscCall(PetscFree(fid_glid_loc));
1268: } else {
1269: PetscCall(PetscMalloc1(nloc, &flid_fgid));
1270: for (kk = 0; kk < nloc; kk++) flid_fgid[kk] = my0 + kk;
1271: }
1272: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0));
1273: /* get P0 */
1274: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0));
1275: {
1276: PetscReal *data_out = NULL;
1277: PetscCall(formProl0(agg_lists, bs, col_bs, myCrs0, nbnodes, data_w_ghost, flid_fgid, &data_out, Prol));
1278: PetscCall(PetscFree(pc_gamg->data));
1280: pc_gamg->data = data_out;
1281: pc_gamg->data_cell_rows = col_bs;
1282: pc_gamg->data_sz = col_bs * col_bs * nLocalSelected;
1283: }
1284: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0));
1285: if (size > 1) PetscCall(PetscFree(data_w_ghost));
1286: PetscCall(PetscFree(flid_fgid));
1288: *a_P_out = Prol; /* out */
1289: PetscCall(MatViewFromOptions(Prol, NULL, "-view_P"));
1291: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1292: PetscFunctionReturn(PETSC_SUCCESS);
1293: }
1295: /*
1296: PCGAMGOptProlongator_AGG
1298: Input Parameter:
1299: . pc - this
1300: . Amat - matrix on this fine level
1301: In/Output Parameter:
1302: . a_P - prolongation operator to the next level
1303: */
1304: static PetscErrorCode PCGAMGOptProlongator_AGG(PC pc, Mat Amat, Mat *a_P)
1305: {
1306: PC_MG *mg = (PC_MG *)pc->data;
1307: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1308: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
1309: PetscInt jj;
1310: Mat Prol = *a_P;
1311: MPI_Comm comm;
1312: KSP eksp;
1313: Vec bb, xx;
1314: PC epc;
1315: PetscReal alpha, emax, emin;
1317: PetscFunctionBegin;
1318: PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
1319: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0));
1321: /* compute maximum singular value of operator to be used in smoother */
1322: if (0 < pc_gamg_agg->nsmooths) {
1323: /* get eigen estimates */
1324: if (pc_gamg->emax > 0) {
1325: emin = pc_gamg->emin;
1326: emax = pc_gamg->emax;
1327: } else {
1328: const char *prefix;
1330: PetscCall(MatCreateVecs(Amat, &bb, NULL));
1331: PetscCall(MatCreateVecs(Amat, &xx, NULL));
1332: PetscCall(KSPSetNoisy_Private(bb));
1334: PetscCall(KSPCreate(comm, &eksp));
1335: PetscCall(KSPSetNestLevel(eksp, pc->kspnestlevel));
1336: PetscCall(PCGetOptionsPrefix(pc, &prefix));
1337: PetscCall(KSPSetOptionsPrefix(eksp, prefix));
1338: PetscCall(KSPAppendOptionsPrefix(eksp, "pc_gamg_esteig_"));
1339: {
1340: PetscBool isset, sflg;
1341: PetscCall(MatIsSPDKnown(Amat, &isset, &sflg));
1342: if (isset && sflg) PetscCall(KSPSetType(eksp, KSPCG));
1343: }
1344: PetscCall(KSPSetErrorIfNotConverged(eksp, pc->erroriffailure));
1345: PetscCall(KSPSetNormType(eksp, KSP_NORM_NONE));
1347: PetscCall(KSPSetInitialGuessNonzero(eksp, PETSC_FALSE));
1348: PetscCall(KSPSetOperators(eksp, Amat, Amat));
1350: PetscCall(KSPGetPC(eksp, &epc));
1351: PetscCall(PCSetType(epc, PCJACOBI)); /* smoother in smoothed agg. */
1353: PetscCall(KSPSetTolerances(eksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, 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
1355: PetscCall(KSPSetFromOptions(eksp));
1356: PetscCall(KSPSetComputeSingularValues(eksp, PETSC_TRUE));
1357: PetscCall(KSPSolve(eksp, bb, xx));
1358: PetscCall(KSPCheckSolve(eksp, pc, xx));
1360: PetscCall(KSPComputeExtremeSingularValues(eksp, &emax, &emin));
1361: PetscCall(PetscInfo(pc, "%s: Smooth P0: max eigen=%e min=%e PC=%s\n", ((PetscObject)pc)->prefix, (double)emax, (double)emin, PCJACOBI));
1362: PetscCall(VecDestroy(&xx));
1363: PetscCall(VecDestroy(&bb));
1364: PetscCall(KSPDestroy(&eksp));
1365: }
1366: if (pc_gamg->use_sa_esteig) {
1367: mg->min_eigen_DinvA[pc_gamg->current_level] = emin;
1368: mg->max_eigen_DinvA[pc_gamg->current_level] = emax;
1369: 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));
1370: } else {
1371: mg->min_eigen_DinvA[pc_gamg->current_level] = 0;
1372: mg->max_eigen_DinvA[pc_gamg->current_level] = 0;
1373: }
1374: } else {
1375: mg->min_eigen_DinvA[pc_gamg->current_level] = 0;
1376: mg->max_eigen_DinvA[pc_gamg->current_level] = 0;
1377: }
1379: /* smooth P0 */
1380: for (jj = 0; jj < pc_gamg_agg->nsmooths; jj++) {
1381: Mat tMat;
1382: Vec diag;
1384: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0));
1386: /* smooth P1 := (I - omega/lam D^{-1}A)P0 */
1387: PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0));
1388: PetscCall(MatMatMult(Amat, Prol, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &tMat));
1389: PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0));
1390: PetscCall(MatProductClear(tMat));
1391: PetscCall(MatCreateVecs(Amat, &diag, NULL));
1392: PetscCall(MatGetDiagonal(Amat, diag)); /* effectively PCJACOBI */
1393: PetscCall(VecReciprocal(diag));
1394: PetscCall(MatDiagonalScale(tMat, diag, NULL));
1395: PetscCall(VecDestroy(&diag));
1397: /* TODO: Set a PCFailedReason and exit the building of the AMG preconditioner */
1398: PetscCheck(emax != 0.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Computed maximum singular value as zero");
1399: /* TODO: Document the 1.4 and don't hardwire it in this routine */
1400: alpha = -1.4 / emax;
1402: PetscCall(MatAYPX(tMat, alpha, Prol, SUBSET_NONZERO_PATTERN));
1403: PetscCall(MatDestroy(&Prol));
1404: Prol = tMat;
1405: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0));
1406: }
1407: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0));
1408: *a_P = Prol;
1409: PetscFunctionReturn(PETSC_SUCCESS);
1410: }
1412: /*
1413: PCCreateGAMG_AGG
1415: Input Parameter:
1416: . pc -
1417: */
1418: PetscErrorCode PCCreateGAMG_AGG(PC pc)
1419: {
1420: PC_MG *mg = (PC_MG *)pc->data;
1421: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1422: PC_GAMG_AGG *pc_gamg_agg;
1424: PetscFunctionBegin;
1425: /* create sub context for SA */
1426: PetscCall(PetscNew(&pc_gamg_agg));
1427: pc_gamg->subctx = pc_gamg_agg;
1429: pc_gamg->ops->setfromoptions = PCSetFromOptions_GAMG_AGG;
1430: pc_gamg->ops->destroy = PCDestroy_GAMG_AGG;
1431: /* reset does not do anything; setup not virtual */
1433: /* set internal function pointers */
1434: pc_gamg->ops->creategraph = PCGAMGCreateGraph_AGG;
1435: pc_gamg->ops->coarsen = PCGAMGCoarsen_AGG;
1436: pc_gamg->ops->prolongator = PCGAMGProlongator_AGG;
1437: pc_gamg->ops->optprolongator = PCGAMGOptProlongator_AGG;
1438: pc_gamg->ops->createdefaultdata = PCSetData_AGG;
1439: pc_gamg->ops->view = PCView_GAMG_AGG;
1441: pc_gamg_agg->nsmooths = 1;
1442: pc_gamg_agg->aggressive_coarsening_levels = 1;
1443: pc_gamg_agg->use_aggressive_square_graph = PETSC_TRUE;
1444: pc_gamg_agg->use_minimum_degree_ordering = PETSC_FALSE;
1445: pc_gamg_agg->use_low_mem_filter = PETSC_FALSE;
1446: pc_gamg_agg->aggressive_mis_k = 2;
1448: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", PCGAMGSetNSmooths_AGG));
1449: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", PCGAMGSetAggressiveLevels_AGG));
1450: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveSquareGraph_C", PCGAMGSetAggressiveSquareGraph_AGG));
1451: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetMinDegreeOrdering_C", PCGAMGMISkSetMinDegreeOrdering_AGG));
1452: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetLowMemoryFilter_C", PCGAMGSetLowMemoryFilter_AGG));
1453: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetAggressive_C", PCGAMGMISkSetAggressive_AGG));
1454: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_AGG));
1455: PetscFunctionReturn(PETSC_SUCCESS);
1456: }