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: MatCoarsen crs;
14: } PC_GAMG_AGG;
16: /*@
17: PCGAMGSetNSmooths - Set number of smoothing steps (1 is typical) used for multigrid on all the levels
19: Logically Collective
21: Input Parameters:
22: + pc - the preconditioner context
23: - n - the number of smooths
25: Options Database Key:
26: . -pc_gamg_agg_nsmooths <nsmooth, default=1> - number of smoothing steps to use with smooth aggregation
28: Level: intermediate
30: .seealso: `PCMG`, `PCGAMG`
31: @*/
32: PetscErrorCode PCGAMGSetNSmooths(PC pc, PetscInt n)
33: {
34: PetscFunctionBegin;
37: PetscTryMethod(pc, "PCGAMGSetNSmooths_C", (PC, PetscInt), (pc, n));
38: PetscFunctionReturn(PETSC_SUCCESS);
39: }
41: static PetscErrorCode PCGAMGSetNSmooths_AGG(PC pc, PetscInt n)
42: {
43: PC_MG *mg = (PC_MG *)pc->data;
44: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
45: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
47: PetscFunctionBegin;
48: pc_gamg_agg->nsmooths = n;
49: PetscFunctionReturn(PETSC_SUCCESS);
50: }
52: /*@
53: PCGAMGSetAggressiveLevels - Use aggressive coarsening on first n levels
55: Logically Collective
57: Input Parameters:
58: + pc - the preconditioner context
59: - n - 0, 1 or more
61: Options Database Key:
62: . -pc_gamg_aggressive_coarsening <n,default = 1> - Number of levels to square the graph on before aggregating it
64: Level: intermediate
66: .seealso: `PCGAMG`, `PCGAMGSetThreshold()`
67: @*/
68: PetscErrorCode PCGAMGSetAggressiveLevels(PC pc, PetscInt n)
69: {
70: PetscFunctionBegin;
73: PetscTryMethod(pc, "PCGAMGSetAggressiveLevels_C", (PC, PetscInt), (pc, n));
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
77: static PetscErrorCode PCGAMGSetAggressiveLevels_AGG(PC pc, PetscInt n)
78: {
79: PC_MG *mg = (PC_MG *)pc->data;
80: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
81: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
83: PetscFunctionBegin;
84: pc_gamg_agg->aggressive_coarsening_levels = n;
85: PetscFunctionReturn(PETSC_SUCCESS);
86: }
88: static PetscErrorCode PCSetFromOptions_GAMG_AGG(PC pc, PetscOptionItems *PetscOptionsObject)
89: {
90: PC_MG *mg = (PC_MG *)pc->data;
91: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
92: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
94: PetscFunctionBegin;
95: PetscOptionsHeadBegin(PetscOptionsObject, "GAMG-AGG options");
96: {
97: PetscBool flg;
98: PetscCall(PetscOptionsInt("-pc_gamg_agg_nsmooths", "smoothing steps for smoothed aggregation, usually 1", "PCGAMGSetNSmooths", pc_gamg_agg->nsmooths, &pc_gamg_agg->nsmooths, NULL));
99: PetscCall(
100: PetscOptionsInt("-pc_gamg_square_graph", "Number of aggressive coarsening (MIS-2) levels from finest (alias for -pc_gamg_aggressive_coarsening, deprecated)", "PCGAMGSetAggressiveLevels", pc_gamg_agg->aggressive_coarsening_levels, &pc_gamg_agg->aggressive_coarsening_levels, &flg));
101: if (!flg) {
102: 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, NULL));
103: } else {
104: 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, &flg));
105: if (flg) 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));
106: }
107: }
108: PetscOptionsHeadEnd();
109: PetscFunctionReturn(PETSC_SUCCESS);
110: }
112: static PetscErrorCode PCDestroy_GAMG_AGG(PC pc)
113: {
114: PC_MG *mg = (PC_MG *)pc->data;
115: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
117: PetscFunctionBegin;
118: PetscCall(PetscFree(pc_gamg->subctx));
119: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", NULL));
120: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", NULL));
121: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
122: PetscFunctionReturn(PETSC_SUCCESS);
123: }
125: /*
126: PCSetCoordinates_AGG
128: Collective
130: Input Parameter:
131: . pc - the preconditioner context
132: . ndm - dimesion of data (used for dof/vertex for Stokes)
133: . a_nloc - number of vertices local
134: . coords - [a_nloc][ndm] - interleaved coordinate data: {x_0, y_0, z_0, x_1, y_1, ...}
135: */
137: static PetscErrorCode PCSetCoordinates_AGG(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords)
138: {
139: PC_MG *mg = (PC_MG *)pc->data;
140: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
141: PetscInt arrsz, kk, ii, jj, nloc, ndatarows, ndf;
142: Mat mat = pc->pmat;
144: PetscFunctionBegin;
147: nloc = a_nloc;
149: /* SA: null space vectors */
150: PetscCall(MatGetBlockSize(mat, &ndf)); /* this does not work for Stokes */
151: if (coords && ndf == 1) pc_gamg->data_cell_cols = 1; /* scalar w/ coords and SA (not needed) */
152: else if (coords) {
153: PetscCheck(ndm <= ndf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "degrees of motion %" PetscInt_FMT " > block size %" PetscInt_FMT, ndm, ndf);
154: pc_gamg->data_cell_cols = (ndm == 2 ? 3 : 6); /* displacement elasticity */
155: 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);
156: } else pc_gamg->data_cell_cols = ndf; /* no data, force SA with constant null space vectors */
157: pc_gamg->data_cell_rows = ndatarows = ndf;
158: 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);
159: arrsz = nloc * pc_gamg->data_cell_rows * pc_gamg->data_cell_cols;
161: if (!pc_gamg->data || (pc_gamg->data_sz != arrsz)) {
162: PetscCall(PetscFree(pc_gamg->data));
163: PetscCall(PetscMalloc1(arrsz + 1, &pc_gamg->data));
164: }
165: /* copy data in - column oriented */
166: for (kk = 0; kk < nloc; kk++) {
167: const PetscInt M = nloc * pc_gamg->data_cell_rows; /* stride into data */
168: PetscReal *data = &pc_gamg->data[kk * ndatarows]; /* start of cell */
169: if (pc_gamg->data_cell_cols == 1) *data = 1.0;
170: else {
171: /* translational modes */
172: for (ii = 0; ii < ndatarows; ii++) {
173: for (jj = 0; jj < ndatarows; jj++) {
174: if (ii == jj) data[ii * M + jj] = 1.0;
175: else data[ii * M + jj] = 0.0;
176: }
177: }
179: /* rotational modes */
180: if (coords) {
181: if (ndm == 2) {
182: data += 2 * M;
183: data[0] = -coords[2 * kk + 1];
184: data[1] = coords[2 * kk];
185: } else {
186: data += 3 * M;
187: data[0] = 0.0;
188: data[M + 0] = coords[3 * kk + 2];
189: data[2 * M + 0] = -coords[3 * kk + 1];
190: data[1] = -coords[3 * kk + 2];
191: data[M + 1] = 0.0;
192: data[2 * M + 1] = coords[3 * kk];
193: data[2] = coords[3 * kk + 1];
194: data[M + 2] = -coords[3 * kk];
195: data[2 * M + 2] = 0.0;
196: }
197: }
198: }
199: }
200: pc_gamg->data_sz = arrsz;
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: /*
205: PCSetData_AGG - called if data is not set with PCSetCoordinates.
206: Looks in Mat for near null space.
207: Does not work for Stokes
209: Input Parameter:
210: . pc -
211: . a_A - matrix to get (near) null space out of.
212: */
213: static PetscErrorCode PCSetData_AGG(PC pc, Mat a_A)
214: {
215: PC_MG *mg = (PC_MG *)pc->data;
216: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
217: MatNullSpace mnull;
219: PetscFunctionBegin;
220: PetscCall(MatGetNearNullSpace(a_A, &mnull));
221: if (!mnull) {
222: DM dm;
223: PetscCall(PCGetDM(pc, &dm));
224: if (!dm) PetscCall(MatGetDM(a_A, &dm));
225: if (dm) {
226: PetscObject deformation;
227: PetscInt Nf;
229: PetscCall(DMGetNumFields(dm, &Nf));
230: if (Nf) {
231: PetscCall(DMGetField(dm, 0, NULL, &deformation));
232: PetscCall(PetscObjectQuery((PetscObject)deformation, "nearnullspace", (PetscObject *)&mnull));
233: if (!mnull) PetscCall(PetscObjectQuery((PetscObject)deformation, "nullspace", (PetscObject *)&mnull));
234: }
235: }
236: }
238: if (!mnull) {
239: PetscInt bs, NN, MM;
240: PetscCall(MatGetBlockSize(a_A, &bs));
241: PetscCall(MatGetLocalSize(a_A, &MM, &NN));
242: PetscCheck(MM % bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MM %" PetscInt_FMT " must be divisible by bs %" PetscInt_FMT, MM, bs);
243: PetscCall(PCSetCoordinates_AGG(pc, bs, MM / bs, NULL));
244: } else {
245: PetscReal *nullvec;
246: PetscBool has_const;
247: PetscInt i, j, mlocal, nvec, bs;
248: const Vec *vecs;
249: const PetscScalar *v;
251: PetscCall(MatGetLocalSize(a_A, &mlocal, NULL));
252: PetscCall(MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs));
253: for (i = 0; i < nvec; i++) {
254: PetscCall(VecGetLocalSize(vecs[i], &j));
255: PetscCheck(j == mlocal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Attached null space vector size %" PetscInt_FMT " != matrix size %" PetscInt_FMT, j, mlocal);
256: }
257: pc_gamg->data_sz = (nvec + !!has_const) * mlocal;
258: PetscCall(PetscMalloc1((nvec + !!has_const) * mlocal, &nullvec));
259: if (has_const)
260: for (i = 0; i < mlocal; i++) nullvec[i] = 1.0;
261: for (i = 0; i < nvec; i++) {
262: PetscCall(VecGetArrayRead(vecs[i], &v));
263: for (j = 0; j < mlocal; j++) nullvec[(i + !!has_const) * mlocal + j] = PetscRealPart(v[j]);
264: PetscCall(VecRestoreArrayRead(vecs[i], &v));
265: }
266: pc_gamg->data = nullvec;
267: pc_gamg->data_cell_cols = (nvec + !!has_const);
268: PetscCall(MatGetBlockSize(a_A, &bs));
269: pc_gamg->data_cell_rows = bs;
270: }
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
274: /*
275: formProl0 - collect null space data for each aggregate, do QR, put R in coarse grid data and Q in P_0
277: Input Parameter:
278: . agg_llists - list of arrays with aggregates -- list from selected vertices of aggregate unselected vertices
279: . bs - row block size
280: . nSAvec - column bs of new P
281: . my0crs - global index of start of locals
282: . data_stride - bs*(nloc nodes + ghost nodes) [data_stride][nSAvec]
283: . data_in[data_stride*nSAvec] - local data on fine grid
284: . flid_fgid[data_stride/bs] - make local to global IDs, includes ghosts in 'locals_llist'
286: Output Parameter:
287: . a_data_out - in with fine grid data (w/ghosts), out with coarse grid data
288: . a_Prol - prolongation operator
289: */
290: 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)
291: {
292: PetscInt Istart, my0, Iend, nloc, clid, flid = 0, aggID, kk, jj, ii, mm, nSelected, minsz, nghosts, out_data_stride;
293: MPI_Comm comm;
294: PetscReal *out_data;
295: PetscCDIntNd *pos;
296: PCGAMGHashTable fgid_flid;
298: PetscFunctionBegin;
299: PetscCall(PetscObjectGetComm((PetscObject)a_Prol, &comm));
300: PetscCall(MatGetOwnershipRange(a_Prol, &Istart, &Iend));
301: nloc = (Iend - Istart) / bs;
302: my0 = Istart / bs;
303: 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);
304: Iend /= bs;
305: nghosts = data_stride / bs - nloc;
307: PetscCall(PCGAMGHashTableCreate(2 * nghosts + 1, &fgid_flid));
308: for (kk = 0; kk < nghosts; kk++) PetscCall(PCGAMGHashTableAdd(&fgid_flid, flid_fgid[nloc + kk], nloc + kk));
310: /* count selected -- same as number of cols of P */
311: for (nSelected = mm = 0; mm < nloc; mm++) {
312: PetscBool ise;
313: PetscCall(PetscCDEmptyAt(agg_llists, mm, &ise));
314: if (!ise) nSelected++;
315: }
316: PetscCall(MatGetOwnershipRangeColumn(a_Prol, &ii, &jj));
317: PetscCheck((ii / nSAvec) == my0crs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ii %" PetscInt_FMT " /nSAvec %" PetscInt_FMT " != my0crs %" PetscInt_FMT, ii, nSAvec, my0crs);
318: 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);
320: /* aloc space for coarse point data (output) */
321: out_data_stride = nSelected * nSAvec;
323: PetscCall(PetscMalloc1(out_data_stride * nSAvec, &out_data));
324: for (ii = 0; ii < out_data_stride * nSAvec; ii++) out_data[ii] = PETSC_MAX_REAL;
325: *a_data_out = out_data; /* output - stride nSelected*nSAvec */
327: /* find points and set prolongation */
328: minsz = 100;
329: for (mm = clid = 0; mm < nloc; mm++) {
330: PetscCall(PetscCDSizeAt(agg_llists, mm, &jj));
331: if (jj > 0) {
332: const PetscInt lid = mm, cgid = my0crs + clid;
333: PetscInt cids[100]; /* max bs */
334: PetscBLASInt asz = jj, M = asz * bs, N = nSAvec, INFO;
335: PetscBLASInt Mdata = M + ((N - M > 0) ? N - M : 0), LDA = Mdata, LWORK = N * bs;
336: PetscScalar *qqc, *qqr, *TAU, *WORK;
337: PetscInt *fids;
338: PetscReal *data;
340: /* count agg */
341: if (asz < minsz) minsz = asz;
343: /* get block */
344: PetscCall(PetscMalloc5(Mdata * N, &qqc, M * N, &qqr, N, &TAU, LWORK, &WORK, M, &fids));
346: aggID = 0;
347: PetscCall(PetscCDGetHeadPos(agg_llists, lid, &pos));
348: while (pos) {
349: PetscInt gid1;
350: PetscCall(PetscCDIntNdGetID(pos, &gid1));
351: PetscCall(PetscCDGetNextPos(agg_llists, lid, &pos));
353: if (gid1 >= my0 && gid1 < Iend) flid = gid1 - my0;
354: else {
355: PetscCall(PCGAMGHashTableFind(&fgid_flid, gid1, &flid));
356: PetscCheck(flid >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cannot find gid1 in table");
357: }
358: /* copy in B_i matrix - column oriented */
359: data = &data_in[flid * bs];
360: for (ii = 0; ii < bs; ii++) {
361: for (jj = 0; jj < N; jj++) {
362: PetscReal d = data[jj * data_stride + ii];
363: qqc[jj * Mdata + aggID * bs + ii] = d;
364: }
365: }
366: /* set fine IDs */
367: for (kk = 0; kk < bs; kk++) fids[aggID * bs + kk] = flid_fgid[flid] * bs + kk;
368: aggID++;
369: }
371: /* pad with zeros */
372: for (ii = asz * bs; ii < Mdata; ii++) {
373: for (jj = 0; jj < N; jj++, kk++) qqc[jj * Mdata + ii] = .0;
374: }
376: /* QR */
377: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
378: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&Mdata, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO));
379: PetscCall(PetscFPTrapPop());
380: PetscCheck(INFO == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xGEQRF error");
381: /* get R - column oriented - output B_{i+1} */
382: {
383: PetscReal *data = &out_data[clid * nSAvec];
384: for (jj = 0; jj < nSAvec; jj++) {
385: for (ii = 0; ii < nSAvec; ii++) {
386: 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);
387: if (ii <= jj) data[jj * out_data_stride + ii] = PetscRealPart(qqc[jj * Mdata + ii]);
388: else data[jj * out_data_stride + ii] = 0.;
389: }
390: }
391: }
393: /* get Q - row oriented */
394: PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&Mdata, &N, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO));
395: PetscCheck(INFO == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xORGQR error arg %" PetscBLASInt_FMT, -INFO);
397: for (ii = 0; ii < M; ii++) {
398: for (jj = 0; jj < N; jj++) qqr[N * ii + jj] = qqc[jj * Mdata + ii];
399: }
401: /* add diagonal block of P0 */
402: for (kk = 0; kk < N; kk++) { cids[kk] = N * cgid + kk; /* global col IDs in P0 */ }
403: PetscCall(MatSetValues(a_Prol, M, fids, N, cids, qqr, INSERT_VALUES));
404: PetscCall(PetscFree5(qqc, qqr, TAU, WORK, fids));
405: clid++;
406: } /* coarse agg */
407: } /* for all fine nodes */
408: PetscCall(MatAssemblyBegin(a_Prol, MAT_FINAL_ASSEMBLY));
409: PetscCall(MatAssemblyEnd(a_Prol, MAT_FINAL_ASSEMBLY));
410: PetscCall(PCGAMGHashTableDestroy(&fgid_flid));
411: PetscFunctionReturn(PETSC_SUCCESS);
412: }
414: static PetscErrorCode PCView_GAMG_AGG(PC pc, PetscViewer viewer)
415: {
416: PC_MG *mg = (PC_MG *)pc->data;
417: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
418: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
420: PetscFunctionBegin;
421: PetscCall(PetscViewerASCIIPrintf(viewer, " AGG specific options\n"));
422: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of levels to square graph %d\n", (int)pc_gamg_agg->aggressive_coarsening_levels));
423: PetscCall(PetscViewerASCIIPrintf(viewer, " Number smoothing steps %d\n", (int)pc_gamg_agg->nsmooths));
424: PetscFunctionReturn(PETSC_SUCCESS);
425: }
427: static PetscErrorCode PCGAMGCreateGraph_AGG(PC pc, Mat Amat, Mat *a_Gmat)
428: {
429: PC_MG *mg = (PC_MG *)pc->data;
430: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
431: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
432: const PetscReal vfilter = pc_gamg->threshold[pc_gamg->current_level];
433: PetscBool ishem;
434: const char *prefix;
436: PetscFunctionBegin;
437: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_GRAPH], 0, 0, 0, 0));
438: /* 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 */
440: /* MATCOARSENHEM requires numerical weights for edges so ensure they are computed */
441: PetscCall(MatCoarsenCreate(PetscObjectComm((PetscObject)pc), &pc_gamg_agg->crs));
442: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix));
443: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc_gamg_agg->crs, prefix));
444: PetscCall(MatCoarsenSetFromOptions(pc_gamg_agg->crs));
445: PetscCall(PetscObjectTypeCompare((PetscObject)pc_gamg_agg->crs, MATCOARSENHEM, &ishem));
447: PetscCall(MatCreateGraph(Amat, PETSC_TRUE, (vfilter >= 0 || ishem) ? PETSC_TRUE : PETSC_FALSE, vfilter, a_Gmat));
448: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_GRAPH], 0, 0, 0, 0));
449: PetscFunctionReturn(PETSC_SUCCESS);
450: }
452: /*
453: PCGAMGCoarsen_AGG - supports squaring the graph (deprecated) and new graph for
454: communication of QR data used with HEM and MISk coarsening
456: Input Parameter:
457: . a_pc - this
459: Input/Output Parameter:
460: . a_Gmat1 - graph to coarsen (in), graph off processor edges for QR gather scatter (out)
462: Output Parameter:
463: . agg_lists - list of aggregates
465: */
466: static PetscErrorCode PCGAMGCoarsen_AGG(PC a_pc, Mat *a_Gmat1, PetscCoarsenData **agg_lists)
467: {
468: PC_MG *mg = (PC_MG *)a_pc->data;
469: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
470: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
471: Mat mat, Gmat1 = *a_Gmat1; /* aggressive graph */
472: IS perm;
473: PetscInt Istart, Iend, Ii, nloc, bs, nn;
474: PetscInt *permute, *degree;
475: PetscBool *bIndexSet;
476: PetscReal hashfact;
477: PetscInt iSwapIndex;
478: PetscRandom random;
480: PetscFunctionBegin;
481: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
482: PetscCall(MatGetLocalSize(Gmat1, &nn, NULL));
483: PetscCall(MatGetBlockSize(Gmat1, &bs));
484: PetscCheck(bs == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "bs %" PetscInt_FMT " must be 1", bs);
485: nloc = nn / bs;
487: /* get MIS aggs - randomize */
488: PetscCall(PetscMalloc2(nloc, &permute, nloc, °ree));
489: PetscCall(PetscCalloc1(nloc, &bIndexSet));
490: for (Ii = 0; Ii < nloc; Ii++) permute[Ii] = Ii;
491: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &random));
492: PetscCall(MatGetOwnershipRange(Gmat1, &Istart, &Iend));
493: for (Ii = 0; Ii < nloc; Ii++) {
494: PetscInt nc;
495: PetscCall(MatGetRow(Gmat1, Istart + Ii, &nc, NULL, NULL));
496: degree[Ii] = nc;
497: PetscCall(MatRestoreRow(Gmat1, Istart + Ii, &nc, NULL, NULL));
498: }
499: for (Ii = 0; Ii < nloc; Ii++) {
500: PetscCall(PetscRandomGetValueReal(random, &hashfact));
501: iSwapIndex = (PetscInt)(hashfact * nloc) % nloc;
502: if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) {
503: PetscInt iTemp = permute[iSwapIndex];
504: permute[iSwapIndex] = permute[Ii];
505: permute[Ii] = iTemp;
506: iTemp = degree[iSwapIndex];
507: degree[iSwapIndex] = degree[Ii];
508: degree[Ii] = iTemp;
509: bIndexSet[iSwapIndex] = PETSC_TRUE;
510: }
511: }
512: // create minimum degree ordering
513: PetscCall(PetscSortIntWithArray(nloc, degree, permute));
515: PetscCall(PetscFree(bIndexSet));
516: PetscCall(PetscRandomDestroy(&random));
517: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_USE_POINTER, &perm));
518: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0));
519: PetscCall(MatCoarsenSetGreedyOrdering(pc_gamg_agg->crs, perm));
520: PetscCall(MatCoarsenSetAdjacency(pc_gamg_agg->crs, Gmat1));
521: PetscCall(MatCoarsenSetStrictAggs(pc_gamg_agg->crs, PETSC_TRUE));
522: if (pc_gamg->current_level < pc_gamg_agg->aggressive_coarsening_levels) PetscCall(MatCoarsenMISKSetDistance(pc_gamg_agg->crs, 2)); // hardwire to MIS-2
523: else PetscCall(MatCoarsenMISKSetDistance(pc_gamg_agg->crs, 1)); // MIS
524: PetscCall(MatCoarsenApply(pc_gamg_agg->crs));
525: PetscCall(MatCoarsenViewFromOptions(pc_gamg_agg->crs, NULL, "-mat_coarsen_view"));
526: PetscCall(MatCoarsenGetData(pc_gamg_agg->crs, agg_lists)); /* output */
527: PetscCall(MatCoarsenDestroy(&pc_gamg_agg->crs));
529: PetscCall(ISDestroy(&perm));
530: PetscCall(PetscFree2(permute, degree));
531: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0));
533: {
534: PetscCoarsenData *llist = *agg_lists;
535: /* see if we have a matrix that takes precedence (returned from MatCoarsenApply) */
536: PetscCall(PetscCDGetMat(llist, &mat));
537: if (mat) {
538: PetscCall(MatDestroy(&Gmat1));
539: *a_Gmat1 = mat; /* output */
540: }
541: }
542: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
543: PetscFunctionReturn(PETSC_SUCCESS);
544: }
546: /*
547: PCGAMGProlongator_AGG
549: Input Parameter:
550: . pc - this
551: . Amat - matrix on this fine level
552: . Graph - used to get ghost data for nodes in
553: . agg_lists - list of aggregates
554: Output Parameter:
555: . a_P_out - prolongation operator to the next level
556: */
557: static PetscErrorCode PCGAMGProlongator_AGG(PC pc, Mat Amat, Mat Gmat, PetscCoarsenData *agg_lists, Mat *a_P_out)
558: {
559: PC_MG *mg = (PC_MG *)pc->data;
560: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
561: const PetscInt col_bs = pc_gamg->data_cell_cols;
562: PetscInt Istart, Iend, nloc, ii, jj, kk, my0, nLocalSelected, bs;
563: Mat Prol;
564: PetscMPIInt size;
565: MPI_Comm comm;
566: PetscReal *data_w_ghost;
567: PetscInt myCrs0, nbnodes = 0, *flid_fgid;
568: MatType mtype;
570: PetscFunctionBegin;
571: PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
572: PetscCheck(col_bs >= 1, comm, PETSC_ERR_PLIB, "Column bs cannot be less than 1");
573: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
574: PetscCallMPI(MPI_Comm_size(comm, &size));
575: PetscCall(MatGetOwnershipRange(Amat, &Istart, &Iend));
576: PetscCall(MatGetBlockSize(Amat, &bs));
577: nloc = (Iend - Istart) / bs;
578: my0 = Istart / bs;
579: 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);
581: /* get 'nLocalSelected' */
582: for (ii = 0, nLocalSelected = 0; ii < nloc; ii++) {
583: PetscBool ise;
584: /* filter out singletons 0 or 1? */
585: PetscCall(PetscCDEmptyAt(agg_lists, ii, &ise));
586: if (!ise) nLocalSelected++;
587: }
589: /* create prolongator, create P matrix */
590: PetscCall(MatGetType(Amat, &mtype));
591: PetscCall(MatCreate(comm, &Prol));
592: PetscCall(MatSetSizes(Prol, nloc * bs, nLocalSelected * col_bs, PETSC_DETERMINE, PETSC_DETERMINE));
593: PetscCall(MatSetBlockSizes(Prol, bs, col_bs));
594: PetscCall(MatSetType(Prol, mtype));
595: PetscCall(MatSeqAIJSetPreallocation(Prol, col_bs, NULL));
596: PetscCall(MatMPIAIJSetPreallocation(Prol, col_bs, NULL, col_bs, NULL));
598: /* can get all points "removed" */
599: PetscCall(MatGetSize(Prol, &kk, &ii));
600: if (!ii) {
601: PetscCall(PetscInfo(pc, "%s: No selected points on coarse grid\n", ((PetscObject)pc)->prefix));
602: PetscCall(MatDestroy(&Prol));
603: *a_P_out = NULL; /* out */
604: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
605: PetscFunctionReturn(PETSC_SUCCESS);
606: }
607: PetscCall(PetscInfo(pc, "%s: New grid %" PetscInt_FMT " nodes\n", ((PetscObject)pc)->prefix, ii / col_bs));
608: PetscCall(MatGetOwnershipRangeColumn(Prol, &myCrs0, &kk));
610: 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);
611: myCrs0 = myCrs0 / col_bs;
612: 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);
614: /* create global vector of data in 'data_w_ghost' */
615: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0));
616: if (size > 1) { /* get ghost null space data */
617: PetscReal *tmp_gdata, *tmp_ldata, *tp2;
618: PetscCall(PetscMalloc1(nloc, &tmp_ldata));
619: for (jj = 0; jj < col_bs; jj++) {
620: for (kk = 0; kk < bs; kk++) {
621: PetscInt ii, stride;
622: const PetscReal *tp = pc_gamg->data + jj * bs * nloc + kk;
623: for (ii = 0; ii < nloc; ii++, tp += bs) tmp_ldata[ii] = *tp;
625: PetscCall(PCGAMGGetDataWithGhosts(Gmat, 1, tmp_ldata, &stride, &tmp_gdata));
627: if (!jj && !kk) { /* now I know how many total nodes - allocate TODO: move below and do in one 'col_bs' call */
628: PetscCall(PetscMalloc1(stride * bs * col_bs, &data_w_ghost));
629: nbnodes = bs * stride;
630: }
631: tp2 = data_w_ghost + jj * bs * stride + kk;
632: for (ii = 0; ii < stride; ii++, tp2 += bs) *tp2 = tmp_gdata[ii];
633: PetscCall(PetscFree(tmp_gdata));
634: }
635: }
636: PetscCall(PetscFree(tmp_ldata));
637: } else {
638: nbnodes = bs * nloc;
639: data_w_ghost = (PetscReal *)pc_gamg->data;
640: }
642: /* get 'flid_fgid' TODO - move up to get 'stride' and do get null space data above in one step (jj loop) */
643: if (size > 1) {
644: PetscReal *fid_glid_loc, *fiddata;
645: PetscInt stride;
647: PetscCall(PetscMalloc1(nloc, &fid_glid_loc));
648: for (kk = 0; kk < nloc; kk++) fid_glid_loc[kk] = (PetscReal)(my0 + kk);
649: PetscCall(PCGAMGGetDataWithGhosts(Gmat, 1, fid_glid_loc, &stride, &fiddata));
650: PetscCall(PetscMalloc1(stride, &flid_fgid)); /* copy real data to in */
651: for (kk = 0; kk < stride; kk++) flid_fgid[kk] = (PetscInt)fiddata[kk];
652: PetscCall(PetscFree(fiddata));
654: PetscCheck(stride == nbnodes / bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "stride %" PetscInt_FMT " != nbnodes %" PetscInt_FMT "/bs %" PetscInt_FMT, stride, nbnodes, bs);
655: PetscCall(PetscFree(fid_glid_loc));
656: } else {
657: PetscCall(PetscMalloc1(nloc, &flid_fgid));
658: for (kk = 0; kk < nloc; kk++) flid_fgid[kk] = my0 + kk;
659: }
660: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0));
661: /* get P0 */
662: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0));
663: {
664: PetscReal *data_out = NULL;
665: PetscCall(formProl0(agg_lists, bs, col_bs, myCrs0, nbnodes, data_w_ghost, flid_fgid, &data_out, Prol));
666: PetscCall(PetscFree(pc_gamg->data));
668: pc_gamg->data = data_out;
669: pc_gamg->data_cell_rows = col_bs;
670: pc_gamg->data_sz = col_bs * col_bs * nLocalSelected;
671: }
672: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0));
673: if (size > 1) PetscCall(PetscFree(data_w_ghost));
674: PetscCall(PetscFree(flid_fgid));
676: *a_P_out = Prol; /* out */
678: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
679: PetscFunctionReturn(PETSC_SUCCESS);
680: }
682: /*
683: PCGAMGOptProlongator_AGG
685: Input Parameter:
686: . pc - this
687: . Amat - matrix on this fine level
688: In/Output Parameter:
689: . a_P - prolongation operator to the next level
690: */
691: static PetscErrorCode PCGAMGOptProlongator_AGG(PC pc, Mat Amat, Mat *a_P)
692: {
693: PC_MG *mg = (PC_MG *)pc->data;
694: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
695: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
696: PetscInt jj;
697: Mat Prol = *a_P;
698: MPI_Comm comm;
699: KSP eksp;
700: Vec bb, xx;
701: PC epc;
702: PetscReal alpha, emax, emin;
704: PetscFunctionBegin;
705: PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
706: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0));
708: /* compute maximum singular value of operator to be used in smoother */
709: if (0 < pc_gamg_agg->nsmooths) {
710: /* get eigen estimates */
711: if (pc_gamg->emax > 0) {
712: emin = pc_gamg->emin;
713: emax = pc_gamg->emax;
714: } else {
715: const char *prefix;
717: PetscCall(MatCreateVecs(Amat, &bb, NULL));
718: PetscCall(MatCreateVecs(Amat, &xx, NULL));
719: PetscCall(KSPSetNoisy_Private(bb));
721: PetscCall(KSPCreate(comm, &eksp));
722: PetscCall(PCGetOptionsPrefix(pc, &prefix));
723: PetscCall(KSPSetOptionsPrefix(eksp, prefix));
724: PetscCall(KSPAppendOptionsPrefix(eksp, "pc_gamg_esteig_"));
725: {
726: PetscBool isset, sflg;
727: PetscCall(MatIsSPDKnown(Amat, &isset, &sflg));
728: if (isset && sflg) PetscCall(KSPSetType(eksp, KSPCG));
729: }
730: PetscCall(KSPSetErrorIfNotConverged(eksp, pc->erroriffailure));
731: PetscCall(KSPSetNormType(eksp, KSP_NORM_NONE));
733: PetscCall(KSPSetInitialGuessNonzero(eksp, PETSC_FALSE));
734: PetscCall(KSPSetOperators(eksp, Amat, Amat));
736: PetscCall(KSPGetPC(eksp, &epc));
737: PetscCall(PCSetType(epc, PCJACOBI)); /* smoother in smoothed agg. */
739: 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
741: PetscCall(KSPSetFromOptions(eksp));
742: PetscCall(KSPSetComputeSingularValues(eksp, PETSC_TRUE));
743: PetscCall(KSPSolve(eksp, bb, xx));
744: PetscCall(KSPCheckSolve(eksp, pc, xx));
746: PetscCall(KSPComputeExtremeSingularValues(eksp, &emax, &emin));
747: PetscCall(PetscInfo(pc, "%s: Smooth P0: max eigen=%e min=%e PC=%s\n", ((PetscObject)pc)->prefix, (double)emax, (double)emin, PCJACOBI));
748: PetscCall(VecDestroy(&xx));
749: PetscCall(VecDestroy(&bb));
750: PetscCall(KSPDestroy(&eksp));
751: }
752: if (pc_gamg->use_sa_esteig) {
753: mg->min_eigen_DinvA[pc_gamg->current_level] = emin;
754: mg->max_eigen_DinvA[pc_gamg->current_level] = emax;
755: 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));
756: } else {
757: mg->min_eigen_DinvA[pc_gamg->current_level] = 0;
758: mg->max_eigen_DinvA[pc_gamg->current_level] = 0;
759: }
760: } else {
761: mg->min_eigen_DinvA[pc_gamg->current_level] = 0;
762: mg->max_eigen_DinvA[pc_gamg->current_level] = 0;
763: }
765: /* smooth P0 */
766: for (jj = 0; jj < pc_gamg_agg->nsmooths; jj++) {
767: Mat tMat;
768: Vec diag;
770: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0));
772: /* smooth P1 := (I - omega/lam D^{-1}A)P0 */
773: PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0));
774: PetscCall(MatMatMult(Amat, Prol, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &tMat));
775: PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0));
776: PetscCall(MatProductClear(tMat));
777: PetscCall(MatCreateVecs(Amat, &diag, NULL));
778: PetscCall(MatGetDiagonal(Amat, diag)); /* effectively PCJACOBI */
779: PetscCall(VecReciprocal(diag));
780: PetscCall(MatDiagonalScale(tMat, diag, NULL));
781: PetscCall(VecDestroy(&diag));
783: /* TODO: Set a PCFailedReason and exit the building of the AMG preconditioner */
784: PetscCheck(emax != 0.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Computed maximum singular value as zero");
785: /* TODO: Document the 1.4 and don't hardwire it in this routine */
786: alpha = -1.4 / emax;
788: PetscCall(MatAYPX(tMat, alpha, Prol, SUBSET_NONZERO_PATTERN));
789: PetscCall(MatDestroy(&Prol));
790: Prol = tMat;
791: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0));
792: }
793: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0));
794: *a_P = Prol;
795: PetscFunctionReturn(PETSC_SUCCESS);
796: }
798: /*
799: PCCreateGAMG_AGG
801: Input Parameter:
802: . pc -
803: */
804: PetscErrorCode PCCreateGAMG_AGG(PC pc)
805: {
806: PC_MG *mg = (PC_MG *)pc->data;
807: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
808: PC_GAMG_AGG *pc_gamg_agg;
810: PetscFunctionBegin;
811: /* create sub context for SA */
812: PetscCall(PetscNew(&pc_gamg_agg));
813: pc_gamg->subctx = pc_gamg_agg;
815: pc_gamg->ops->setfromoptions = PCSetFromOptions_GAMG_AGG;
816: pc_gamg->ops->destroy = PCDestroy_GAMG_AGG;
817: /* reset does not do anything; setup not virtual */
819: /* set internal function pointers */
820: pc_gamg->ops->creategraph = PCGAMGCreateGraph_AGG;
821: pc_gamg->ops->coarsen = PCGAMGCoarsen_AGG;
822: pc_gamg->ops->prolongator = PCGAMGProlongator_AGG;
823: pc_gamg->ops->optprolongator = PCGAMGOptProlongator_AGG;
824: pc_gamg->ops->createdefaultdata = PCSetData_AGG;
825: pc_gamg->ops->view = PCView_GAMG_AGG;
827: pc_gamg_agg->aggressive_coarsening_levels = 1;
828: pc_gamg_agg->nsmooths = 1;
830: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", PCGAMGSetNSmooths_AGG));
831: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", PCGAMGSetAggressiveLevels_AGG));
832: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_AGG));
833: PetscFunctionReturn(PETSC_SUCCESS);
834: }