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