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: MatCoarsen crs;
17: } PC_GAMG_AGG;
19: /*@
20: PCGAMGSetNSmooths - Set number of smoothing steps (1 is typical) used for multigrid on all the levels
22: Logically Collective
24: Input Parameters:
25: + pc - the preconditioner context
26: - n - the number of smooths
28: Options Database Key:
29: . -pc_gamg_agg_nsmooths <nsmooth, default=1> - number of smoothing steps to use with smooth aggregation
31: Level: intermediate
33: .seealso: [](ch_ksp), `PCMG`, `PCGAMG`
34: @*/
35: PetscErrorCode PCGAMGSetNSmooths(PC pc, PetscInt n)
36: {
37: PetscFunctionBegin;
40: PetscTryMethod(pc, "PCGAMGSetNSmooths_C", (PC, PetscInt), (pc, n));
41: PetscFunctionReturn(PETSC_SUCCESS);
42: }
44: static PetscErrorCode PCGAMGSetNSmooths_AGG(PC pc, PetscInt n)
45: {
46: PC_MG *mg = (PC_MG *)pc->data;
47: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
48: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
50: PetscFunctionBegin;
51: pc_gamg_agg->nsmooths = n;
52: PetscFunctionReturn(PETSC_SUCCESS);
53: }
55: /*@
56: PCGAMGSetAggressiveLevels - Use aggressive coarsening on first n levels
58: Logically Collective
60: Input Parameters:
61: + pc - the preconditioner context
62: - n - 0, 1 or more
64: Options Database Key:
65: . -pc_gamg_aggressive_coarsening <n,default = 1> - Number of levels to square the graph on before aggregating it
67: Level: intermediate
69: .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`
70: @*/
71: PetscErrorCode PCGAMGSetAggressiveLevels(PC pc, PetscInt n)
72: {
73: PetscFunctionBegin;
76: PetscTryMethod(pc, "PCGAMGSetAggressiveLevels_C", (PC, PetscInt), (pc, n));
77: PetscFunctionReturn(PETSC_SUCCESS);
78: }
80: /*@
81: PCGAMGMISkSetAggressive - Number (k) distance in MIS coarsening (>2 is 'aggressive')
83: Logically Collective
85: Input Parameters:
86: + pc - the preconditioner context
87: - n - 1 or more (default = 2)
89: Options Database Key:
90: . -pc_gamg_aggressive_mis_k <n,default=2> - Number (k) distance in MIS coarsening (>2 is 'aggressive')
92: Level: intermediate
94: .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`
95: @*/
96: PetscErrorCode PCGAMGMISkSetAggressive(PC pc, PetscInt n)
97: {
98: PetscFunctionBegin;
101: PetscTryMethod(pc, "PCGAMGMISkSetAggressive_C", (PC, PetscInt), (pc, n));
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: /*@
106: PCGAMGSetAggressiveSquareGraph - Use graph square A'A for aggressive coarsening, old method
108: Logically Collective
110: Input Parameters:
111: + pc - the preconditioner context
112: - b - default false - MIS-k is faster
114: Options Database Key:
115: . -pc_gamg_aggressive_square_graph <bool,default=false> - Use square graph (A'A) or MIS-k (k=2) for aggressive coarsening
117: Level: intermediate
119: .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`, `PCGAMGMISkSetAggressive()`, `PCGAMGMISkSetMinDegreeOrdering()`
120: @*/
121: PetscErrorCode PCGAMGSetAggressiveSquareGraph(PC pc, PetscBool b)
122: {
123: PetscFunctionBegin;
126: PetscTryMethod(pc, "PCGAMGSetAggressiveSquareGraph_C", (PC, PetscBool), (pc, b));
127: PetscFunctionReturn(PETSC_SUCCESS);
128: }
130: /*@
131: PCGAMGMISkSetMinDegreeOrdering - Use minimum degree ordering in greedy MIS algorithm
133: Logically Collective
135: Input Parameters:
136: + pc - the preconditioner context
137: - b - default true
139: Options Database Key:
140: . -pc_gamg_mis_k_minimum_degree_ordering <bool,default=true> - Use minimum degree ordering in greedy MIS algorithm
142: Level: intermediate
144: .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`, `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`
145: @*/
146: PetscErrorCode PCGAMGMISkSetMinDegreeOrdering(PC pc, PetscBool b)
147: {
148: PetscFunctionBegin;
151: PetscTryMethod(pc, "PCGAMGMISkSetMinDegreeOrdering_C", (PC, PetscBool), (pc, b));
152: PetscFunctionReturn(PETSC_SUCCESS);
153: }
155: static PetscErrorCode PCGAMGSetAggressiveLevels_AGG(PC pc, PetscInt n)
156: {
157: PC_MG *mg = (PC_MG *)pc->data;
158: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
159: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
161: PetscFunctionBegin;
162: pc_gamg_agg->aggressive_coarsening_levels = n;
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: static PetscErrorCode PCGAMGMISkSetAggressive_AGG(PC pc, PetscInt n)
167: {
168: PC_MG *mg = (PC_MG *)pc->data;
169: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
170: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
172: PetscFunctionBegin;
173: pc_gamg_agg->aggressive_mis_k = n;
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: static PetscErrorCode PCGAMGSetAggressiveSquareGraph_AGG(PC pc, PetscBool b)
178: {
179: PC_MG *mg = (PC_MG *)pc->data;
180: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
181: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
183: PetscFunctionBegin;
184: pc_gamg_agg->use_aggressive_square_graph = b;
185: PetscFunctionReturn(PETSC_SUCCESS);
186: }
188: static PetscErrorCode PCGAMGMISkSetMinDegreeOrdering_AGG(PC pc, PetscBool b)
189: {
190: PC_MG *mg = (PC_MG *)pc->data;
191: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
192: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
194: PetscFunctionBegin;
195: pc_gamg_agg->use_minimum_degree_ordering = b;
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
199: static PetscErrorCode PCSetFromOptions_GAMG_AGG(PC pc, PetscOptionItems *PetscOptionsObject)
200: {
201: PC_MG *mg = (PC_MG *)pc->data;
202: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
203: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
205: PetscFunctionBegin;
206: PetscOptionsHeadBegin(PetscOptionsObject, "GAMG-AGG options");
207: {
208: PetscBool flg;
209: PetscCall(PetscOptionsInt("-pc_gamg_agg_nsmooths", "smoothing steps for smoothed aggregation, usually 1", "PCGAMGSetNSmooths", pc_gamg_agg->nsmooths, &pc_gamg_agg->nsmooths, NULL));
210: 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));
211: if (!flg) {
212: PetscCall(
213: PetscOptionsInt("-pc_gamg_square_graph", "Number of aggressive coarsening (MIS-2) levels from finest (deprecated alias for -pc_gamg_aggressive_coarsening)", "PCGAMGSetAggressiveLevels", pc_gamg_agg->aggressive_coarsening_levels, &pc_gamg_agg->aggressive_coarsening_levels, NULL));
214: } else {
215: PetscCall(
216: 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));
217: 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));
218: }
219: if (pc_gamg_agg->aggressive_coarsening_levels > 0) {
220: PetscCall(PetscOptionsBool("-pc_gamg_aggressive_square_graph", "Use square graph (A'A) or MIS-k (k=2) for aggressive coarsening", "PCGAMGSetAggressiveSquareGraph", pc_gamg_agg->use_aggressive_square_graph, &pc_gamg_agg->use_aggressive_square_graph, NULL));
221: }
222: 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));
223: 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));
224: }
225: PetscOptionsHeadEnd();
226: PetscFunctionReturn(PETSC_SUCCESS);
227: }
229: static PetscErrorCode PCDestroy_GAMG_AGG(PC pc)
230: {
231: PC_MG *mg = (PC_MG *)pc->data;
232: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
234: PetscFunctionBegin;
235: PetscCall(PetscFree(pc_gamg->subctx));
236: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", NULL));
237: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", NULL));
238: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetAggressive_C", NULL));
239: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetMinDegreeOrdering_C", NULL));
240: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveSquareGraph_C", NULL));
241: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }
245: /*
246: PCSetCoordinates_AGG
248: Collective
250: Input Parameter:
251: . pc - the preconditioner context
252: . ndm - dimension of data (used for dof/vertex for Stokes)
253: . a_nloc - number of vertices local
254: . coords - [a_nloc][ndm] - interleaved coordinate data: {x_0, y_0, z_0, x_1, y_1, ...}
255: */
257: static PetscErrorCode PCSetCoordinates_AGG(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords)
258: {
259: PC_MG *mg = (PC_MG *)pc->data;
260: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
261: PetscInt arrsz, kk, ii, jj, nloc, ndatarows, ndf;
262: Mat mat = pc->pmat;
264: PetscFunctionBegin;
267: nloc = a_nloc;
269: /* SA: null space vectors */
270: PetscCall(MatGetBlockSize(mat, &ndf)); /* this does not work for Stokes */
271: if (coords && ndf == 1) pc_gamg->data_cell_cols = 1; /* scalar w/ coords and SA (not needed) */
272: else if (coords) {
273: PetscCheck(ndm <= ndf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "degrees of motion %" PetscInt_FMT " > block size %" PetscInt_FMT, ndm, ndf);
274: pc_gamg->data_cell_cols = (ndm == 2 ? 3 : 6); /* displacement elasticity */
275: 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);
276: } else pc_gamg->data_cell_cols = ndf; /* no data, force SA with constant null space vectors */
277: pc_gamg->data_cell_rows = ndatarows = ndf;
278: 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);
279: arrsz = nloc * pc_gamg->data_cell_rows * pc_gamg->data_cell_cols;
281: if (!pc_gamg->data || (pc_gamg->data_sz != arrsz)) {
282: PetscCall(PetscFree(pc_gamg->data));
283: PetscCall(PetscMalloc1(arrsz + 1, &pc_gamg->data));
284: }
285: /* copy data in - column oriented */
286: for (kk = 0; kk < nloc; kk++) {
287: const PetscInt M = nloc * pc_gamg->data_cell_rows; /* stride into data */
288: PetscReal *data = &pc_gamg->data[kk * ndatarows]; /* start of cell */
289: if (pc_gamg->data_cell_cols == 1) *data = 1.0;
290: else {
291: /* translational modes */
292: for (ii = 0; ii < ndatarows; ii++) {
293: for (jj = 0; jj < ndatarows; jj++) {
294: if (ii == jj) data[ii * M + jj] = 1.0;
295: else data[ii * M + jj] = 0.0;
296: }
297: }
299: /* rotational modes */
300: if (coords) {
301: if (ndm == 2) {
302: data += 2 * M;
303: data[0] = -coords[2 * kk + 1];
304: data[1] = coords[2 * kk];
305: } else {
306: data += 3 * M;
307: data[0] = 0.0;
308: data[M + 0] = coords[3 * kk + 2];
309: data[2 * M + 0] = -coords[3 * kk + 1];
310: data[1] = -coords[3 * kk + 2];
311: data[M + 1] = 0.0;
312: data[2 * M + 1] = coords[3 * kk];
313: data[2] = coords[3 * kk + 1];
314: data[M + 2] = -coords[3 * kk];
315: data[2 * M + 2] = 0.0;
316: }
317: }
318: }
319: }
320: pc_gamg->data_sz = arrsz;
321: PetscFunctionReturn(PETSC_SUCCESS);
322: }
324: /*
325: PCSetData_AGG - called if data is not set with PCSetCoordinates.
326: Looks in Mat for near null space.
327: Does not work for Stokes
329: Input Parameter:
330: . pc -
331: . a_A - matrix to get (near) null space out of.
332: */
333: static PetscErrorCode PCSetData_AGG(PC pc, Mat a_A)
334: {
335: PC_MG *mg = (PC_MG *)pc->data;
336: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
337: MatNullSpace mnull;
339: PetscFunctionBegin;
340: PetscCall(MatGetNearNullSpace(a_A, &mnull));
341: if (!mnull) {
342: DM dm;
343: PetscCall(PCGetDM(pc, &dm));
344: if (!dm) PetscCall(MatGetDM(a_A, &dm));
345: if (dm) {
346: PetscObject deformation;
347: PetscInt Nf;
349: PetscCall(DMGetNumFields(dm, &Nf));
350: if (Nf) {
351: PetscCall(DMGetField(dm, 0, NULL, &deformation));
352: PetscCall(PetscObjectQuery((PetscObject)deformation, "nearnullspace", (PetscObject *)&mnull));
353: if (!mnull) PetscCall(PetscObjectQuery((PetscObject)deformation, "nullspace", (PetscObject *)&mnull));
354: }
355: }
356: }
358: if (!mnull) {
359: PetscInt bs, NN, MM;
360: PetscCall(MatGetBlockSize(a_A, &bs));
361: PetscCall(MatGetLocalSize(a_A, &MM, &NN));
362: PetscCheck(MM % bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MM %" PetscInt_FMT " must be divisible by bs %" PetscInt_FMT, MM, bs);
363: PetscCall(PCSetCoordinates_AGG(pc, bs, MM / bs, NULL));
364: } else {
365: PetscReal *nullvec;
366: PetscBool has_const;
367: PetscInt i, j, mlocal, nvec, bs;
368: const Vec *vecs;
369: const PetscScalar *v;
371: PetscCall(MatGetLocalSize(a_A, &mlocal, NULL));
372: PetscCall(MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs));
373: for (i = 0; i < nvec; i++) {
374: PetscCall(VecGetLocalSize(vecs[i], &j));
375: PetscCheck(j == mlocal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Attached null space vector size %" PetscInt_FMT " != matrix size %" PetscInt_FMT, j, mlocal);
376: }
377: pc_gamg->data_sz = (nvec + !!has_const) * mlocal;
378: PetscCall(PetscMalloc1((nvec + !!has_const) * mlocal, &nullvec));
379: if (has_const)
380: for (i = 0; i < mlocal; i++) nullvec[i] = 1.0;
381: for (i = 0; i < nvec; i++) {
382: PetscCall(VecGetArrayRead(vecs[i], &v));
383: for (j = 0; j < mlocal; j++) nullvec[(i + !!has_const) * mlocal + j] = PetscRealPart(v[j]);
384: PetscCall(VecRestoreArrayRead(vecs[i], &v));
385: }
386: pc_gamg->data = nullvec;
387: pc_gamg->data_cell_cols = (nvec + !!has_const);
388: PetscCall(MatGetBlockSize(a_A, &bs));
389: pc_gamg->data_cell_rows = bs;
390: }
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
394: /*
395: formProl0 - collect null space data for each aggregate, do QR, put R in coarse grid data and Q in P_0
397: Input Parameter:
398: . agg_llists - list of arrays with aggregates -- list from selected vertices of aggregate unselected vertices
399: . bs - row block size
400: . nSAvec - column bs of new P
401: . my0crs - global index of start of locals
402: . data_stride - bs*(nloc nodes + ghost nodes) [data_stride][nSAvec]
403: . data_in[data_stride*nSAvec] - local data on fine grid
404: . flid_fgid[data_stride/bs] - make local to global IDs, includes ghosts in 'locals_llist'
406: Output Parameter:
407: . a_data_out - in with fine grid data (w/ghosts), out with coarse grid data
408: . a_Prol - prolongation operator
409: */
410: 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)
411: {
412: PetscInt Istart, my0, Iend, nloc, clid, flid = 0, aggID, kk, jj, ii, mm, nSelected, minsz, nghosts, out_data_stride;
413: MPI_Comm comm;
414: PetscReal *out_data;
415: PetscCDIntNd *pos;
416: PCGAMGHashTable fgid_flid;
418: PetscFunctionBegin;
419: PetscCall(PetscObjectGetComm((PetscObject)a_Prol, &comm));
420: PetscCall(MatGetOwnershipRange(a_Prol, &Istart, &Iend));
421: nloc = (Iend - Istart) / bs;
422: my0 = Istart / bs;
423: 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);
424: Iend /= bs;
425: nghosts = data_stride / bs - nloc;
427: PetscCall(PCGAMGHashTableCreate(2 * nghosts + 1, &fgid_flid));
428: for (kk = 0; kk < nghosts; kk++) PetscCall(PCGAMGHashTableAdd(&fgid_flid, flid_fgid[nloc + kk], nloc + kk));
430: /* count selected -- same as number of cols of P */
431: for (nSelected = mm = 0; mm < nloc; mm++) {
432: PetscBool ise;
433: PetscCall(PetscCDEmptyAt(agg_llists, mm, &ise));
434: if (!ise) nSelected++;
435: }
436: PetscCall(MatGetOwnershipRangeColumn(a_Prol, &ii, &jj));
437: PetscCheck((ii / nSAvec) == my0crs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ii %" PetscInt_FMT " /nSAvec %" PetscInt_FMT " != my0crs %" PetscInt_FMT, ii, nSAvec, my0crs);
438: 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);
440: /* aloc space for coarse point data (output) */
441: out_data_stride = nSelected * nSAvec;
443: PetscCall(PetscMalloc1(out_data_stride * nSAvec, &out_data));
444: for (ii = 0; ii < out_data_stride * nSAvec; ii++) out_data[ii] = PETSC_MAX_REAL;
445: *a_data_out = out_data; /* output - stride nSelected*nSAvec */
447: /* find points and set prolongation */
448: minsz = 100;
449: for (mm = clid = 0; mm < nloc; mm++) {
450: PetscCall(PetscCDSizeAt(agg_llists, mm, &jj));
451: if (jj > 0) {
452: const PetscInt lid = mm, cgid = my0crs + clid;
453: PetscInt cids[100]; /* max bs */
454: PetscBLASInt asz = jj, M = asz * bs, N = nSAvec, INFO;
455: PetscBLASInt Mdata = M + ((N - M > 0) ? N - M : 0), LDA = Mdata, LWORK = N * bs;
456: PetscScalar *qqc, *qqr, *TAU, *WORK;
457: PetscInt *fids;
458: PetscReal *data;
460: /* count agg */
461: if (asz < minsz) minsz = asz;
463: /* get block */
464: PetscCall(PetscMalloc5(Mdata * N, &qqc, M * N, &qqr, N, &TAU, LWORK, &WORK, M, &fids));
466: aggID = 0;
467: PetscCall(PetscCDGetHeadPos(agg_llists, lid, &pos));
468: while (pos) {
469: PetscInt gid1;
470: PetscCall(PetscCDIntNdGetID(pos, &gid1));
471: PetscCall(PetscCDGetNextPos(agg_llists, lid, &pos));
473: if (gid1 >= my0 && gid1 < Iend) flid = gid1 - my0;
474: else {
475: PetscCall(PCGAMGHashTableFind(&fgid_flid, gid1, &flid));
476: PetscCheck(flid >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cannot find gid1 in table");
477: }
478: /* copy in B_i matrix - column oriented */
479: data = &data_in[flid * bs];
480: for (ii = 0; ii < bs; ii++) {
481: for (jj = 0; jj < N; jj++) {
482: PetscReal d = data[jj * data_stride + ii];
483: qqc[jj * Mdata + aggID * bs + ii] = d;
484: }
485: }
486: /* set fine IDs */
487: for (kk = 0; kk < bs; kk++) fids[aggID * bs + kk] = flid_fgid[flid] * bs + kk;
488: aggID++;
489: }
491: /* pad with zeros */
492: for (ii = asz * bs; ii < Mdata; ii++) {
493: for (jj = 0; jj < N; jj++, kk++) qqc[jj * Mdata + ii] = .0;
494: }
496: /* QR */
497: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
498: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&Mdata, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO));
499: PetscCall(PetscFPTrapPop());
500: PetscCheck(INFO == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xGEQRF error");
501: /* get R - column oriented - output B_{i+1} */
502: {
503: PetscReal *data = &out_data[clid * nSAvec];
504: for (jj = 0; jj < nSAvec; jj++) {
505: for (ii = 0; ii < nSAvec; ii++) {
506: 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);
507: if (ii <= jj) data[jj * out_data_stride + ii] = PetscRealPart(qqc[jj * Mdata + ii]);
508: else data[jj * out_data_stride + ii] = 0.;
509: }
510: }
511: }
513: /* get Q - row oriented */
514: PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&Mdata, &N, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO));
515: PetscCheck(INFO == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xORGQR error arg %" PetscBLASInt_FMT, -INFO);
517: for (ii = 0; ii < M; ii++) {
518: for (jj = 0; jj < N; jj++) qqr[N * ii + jj] = qqc[jj * Mdata + ii];
519: }
521: /* add diagonal block of P0 */
522: for (kk = 0; kk < N; kk++) { cids[kk] = N * cgid + kk; /* global col IDs in P0 */ }
523: PetscCall(MatSetValues(a_Prol, M, fids, N, cids, qqr, INSERT_VALUES));
524: PetscCall(PetscFree5(qqc, qqr, TAU, WORK, fids));
525: clid++;
526: } /* coarse agg */
527: } /* for all fine nodes */
528: PetscCall(MatAssemblyBegin(a_Prol, MAT_FINAL_ASSEMBLY));
529: PetscCall(MatAssemblyEnd(a_Prol, MAT_FINAL_ASSEMBLY));
530: PetscCall(PCGAMGHashTableDestroy(&fgid_flid));
531: PetscFunctionReturn(PETSC_SUCCESS);
532: }
534: static PetscErrorCode PCView_GAMG_AGG(PC pc, PetscViewer viewer)
535: {
536: PC_MG *mg = (PC_MG *)pc->data;
537: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
538: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
540: PetscFunctionBegin;
541: PetscCall(PetscViewerASCIIPrintf(viewer, " AGG specific options\n"));
542: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of levels of aggressive coarsening %d\n", (int)pc_gamg_agg->aggressive_coarsening_levels));
543: if (pc_gamg_agg->aggressive_coarsening_levels > 0) {
544: PetscCall(PetscViewerASCIIPrintf(viewer, " %s aggressive coarsening\n", !pc_gamg_agg->use_aggressive_square_graph ? "MIS-k" : "Square graph"));
545: 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));
546: }
547: PetscCall(PetscViewerASCIIPrintf(viewer, " Number smoothing steps %d\n", (int)pc_gamg_agg->nsmooths));
548: PetscFunctionReturn(PETSC_SUCCESS);
549: }
551: static PetscErrorCode PCGAMGCreateGraph_AGG(PC pc, Mat Amat, Mat *a_Gmat)
552: {
553: PC_MG *mg = (PC_MG *)pc->data;
554: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
555: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
556: const PetscReal vfilter = pc_gamg->threshold[pc_gamg->current_level];
557: PetscBool ishem;
558: const char *prefix;
559: MatInfo info0, info1;
560: PetscInt bs;
562: PetscFunctionBegin;
563: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_GRAPH], 0, 0, 0, 0));
564: /* 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 */
566: /* MATCOARSENHEM requires numerical weights for edges so ensure they are computed */
567: PetscCall(MatCoarsenCreate(PetscObjectComm((PetscObject)pc), &pc_gamg_agg->crs));
568: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix));
569: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc_gamg_agg->crs, prefix));
570: PetscCall(MatCoarsenSetFromOptions(pc_gamg_agg->crs));
571: PetscCall(PetscObjectTypeCompare((PetscObject)pc_gamg_agg->crs, MATCOARSENHEM, &ishem));
572: if (ishem) pc_gamg_agg->aggressive_coarsening_levels = 0; // aggressive and HEM does not make sense
573: PetscCall(MatGetInfo(Amat, MAT_LOCAL, &info0)); /* global reduction */
574: PetscCall(MatCreateGraph(Amat, PETSC_TRUE, (vfilter >= 0 || ishem) ? PETSC_TRUE : PETSC_FALSE, vfilter, a_Gmat));
575: PetscCall(MatGetInfo(*a_Gmat, MAT_LOCAL, &info1)); /* global reduction */
576: PetscCall(MatGetBlockSize(Amat, &bs));
577: 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));
578: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_GRAPH], 0, 0, 0, 0));
579: PetscFunctionReturn(PETSC_SUCCESS);
580: }
582: typedef PetscInt NState;
583: static const NState NOT_DONE = -2;
584: static const NState DELETED = -1;
585: static const NState REMOVED = -3;
586: #define IS_SELECTED(s) (s != DELETED && s != NOT_DONE && s != REMOVED)
588: /*
589: fixAggregatesWithSquare - greedy grab of with G1 (unsquared graph) -- AIJ specific -- change to fixAggregatesWithSquare -- TODD
590: - AGG-MG specific: clears singletons out of 'selected_2'
592: Input Parameter:
593: . Gmat_2 - global matrix of squared graph (data not defined)
594: . Gmat_1 - base graph to grab with base graph
595: Input/Output Parameter:
596: . aggs_2 - linked list of aggs with gids)
597: */
598: static PetscErrorCode fixAggregatesWithSquare(PC pc, Mat Gmat_2, Mat Gmat_1, PetscCoarsenData *aggs_2)
599: {
600: PetscBool isMPI;
601: Mat_SeqAIJ *matA_1, *matB_1 = NULL;
602: MPI_Comm comm;
603: PetscInt lid, *ii, *idx, ix, Iend, my0, kk, n, j;
604: Mat_MPIAIJ *mpimat_2 = NULL, *mpimat_1 = NULL;
605: const PetscInt nloc = Gmat_2->rmap->n;
606: PetscScalar *cpcol_1_state, *cpcol_2_state, *cpcol_2_par_orig, *lid_parent_gid;
607: PetscInt *lid_cprowID_1 = NULL;
608: NState *lid_state;
609: Vec ghost_par_orig2;
610: PetscMPIInt rank;
612: PetscFunctionBegin;
613: PetscCall(PetscObjectGetComm((PetscObject)Gmat_2, &comm));
614: PetscCallMPI(MPI_Comm_rank(comm, &rank));
615: PetscCall(MatGetOwnershipRange(Gmat_1, &my0, &Iend));
617: /* get submatrices */
618: PetscCall(PetscStrbeginswith(((PetscObject)Gmat_1)->type_name, MATMPIAIJ, &isMPI));
619: PetscCall(PetscInfo(pc, "isMPI = %s\n", isMPI ? "yes" : "no"));
620: PetscCall(PetscMalloc3(nloc, &lid_state, nloc, &lid_parent_gid, nloc, &lid_cprowID_1));
621: for (lid = 0; lid < nloc; lid++) lid_cprowID_1[lid] = -1;
622: if (isMPI) {
623: /* grab matrix objects */
624: mpimat_2 = (Mat_MPIAIJ *)Gmat_2->data;
625: mpimat_1 = (Mat_MPIAIJ *)Gmat_1->data;
626: matA_1 = (Mat_SeqAIJ *)mpimat_1->A->data;
627: matB_1 = (Mat_SeqAIJ *)mpimat_1->B->data;
629: /* force compressed row storage for B matrix in AuxMat */
630: PetscCall(MatCheckCompressedRow(mpimat_1->B, matB_1->nonzerorowcnt, &matB_1->compressedrow, matB_1->i, Gmat_1->rmap->n, -1.0));
631: for (ix = 0; ix < matB_1->compressedrow.nrows; ix++) {
632: PetscInt lid = matB_1->compressedrow.rindex[ix];
633: PetscCheck(lid <= nloc && lid >= -1, PETSC_COMM_SELF, PETSC_ERR_USER, "lid %d out of range. nloc = %d", (int)lid, (int)nloc);
634: if (lid != -1) lid_cprowID_1[lid] = ix;
635: }
636: } else {
637: PetscBool isAIJ;
638: PetscCall(PetscStrbeginswith(((PetscObject)Gmat_1)->type_name, MATSEQAIJ, &isAIJ));
639: PetscCheck(isAIJ, PETSC_COMM_SELF, PETSC_ERR_USER, "Require AIJ matrix.");
640: matA_1 = (Mat_SeqAIJ *)Gmat_1->data;
641: }
642: if (nloc > 0) { PetscCheck(!matB_1 || matB_1->compressedrow.use, PETSC_COMM_SELF, PETSC_ERR_PLIB, "matB_1 && !matB_1->compressedrow.use: PETSc bug???"); }
643: /* get state of locals and selected gid for deleted */
644: for (lid = 0; lid < nloc; lid++) {
645: lid_parent_gid[lid] = -1.0;
646: lid_state[lid] = DELETED;
647: }
649: /* set lid_state */
650: for (lid = 0; lid < nloc; lid++) {
651: PetscCDIntNd *pos;
652: PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
653: if (pos) {
654: PetscInt gid1;
656: PetscCall(PetscCDIntNdGetID(pos, &gid1));
657: PetscCheck(gid1 == lid + my0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "gid1 %d != lid %d + my0 %d", (int)gid1, (int)lid, (int)my0);
658: lid_state[lid] = gid1;
659: }
660: }
662: /* map local to selected local, DELETED means a ghost owns it */
663: for (lid = kk = 0; lid < nloc; lid++) {
664: NState state = lid_state[lid];
665: if (IS_SELECTED(state)) {
666: PetscCDIntNd *pos;
667: PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
668: while (pos) {
669: PetscInt gid1;
670: PetscCall(PetscCDIntNdGetID(pos, &gid1));
671: PetscCall(PetscCDGetNextPos(aggs_2, lid, &pos));
672: if (gid1 >= my0 && gid1 < Iend) lid_parent_gid[gid1 - my0] = (PetscScalar)(lid + my0);
673: }
674: }
675: }
676: /* get 'cpcol_1/2_state' & cpcol_2_par_orig - uses mpimat_1/2->lvec for temp space */
677: if (isMPI) {
678: Vec tempVec;
679: /* get 'cpcol_1_state' */
680: PetscCall(MatCreateVecs(Gmat_1, &tempVec, NULL));
681: for (kk = 0, j = my0; kk < nloc; kk++, j++) {
682: PetscScalar v = (PetscScalar)lid_state[kk];
683: PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
684: }
685: PetscCall(VecAssemblyBegin(tempVec));
686: PetscCall(VecAssemblyEnd(tempVec));
687: PetscCall(VecScatterBegin(mpimat_1->Mvctx, tempVec, mpimat_1->lvec, INSERT_VALUES, SCATTER_FORWARD));
688: PetscCall(VecScatterEnd(mpimat_1->Mvctx, tempVec, mpimat_1->lvec, INSERT_VALUES, SCATTER_FORWARD));
689: PetscCall(VecGetArray(mpimat_1->lvec, &cpcol_1_state));
690: /* get 'cpcol_2_state' */
691: PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, mpimat_2->lvec, INSERT_VALUES, SCATTER_FORWARD));
692: PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, mpimat_2->lvec, INSERT_VALUES, SCATTER_FORWARD));
693: PetscCall(VecGetArray(mpimat_2->lvec, &cpcol_2_state));
694: /* get 'cpcol_2_par_orig' */
695: for (kk = 0, j = my0; kk < nloc; kk++, j++) {
696: PetscScalar v = (PetscScalar)lid_parent_gid[kk];
697: PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
698: }
699: PetscCall(VecAssemblyBegin(tempVec));
700: PetscCall(VecAssemblyEnd(tempVec));
701: PetscCall(VecDuplicate(mpimat_2->lvec, &ghost_par_orig2));
702: PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghost_par_orig2, INSERT_VALUES, SCATTER_FORWARD));
703: PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghost_par_orig2, INSERT_VALUES, SCATTER_FORWARD));
704: PetscCall(VecGetArray(ghost_par_orig2, &cpcol_2_par_orig));
706: PetscCall(VecDestroy(&tempVec));
707: } /* ismpi */
708: for (lid = 0; lid < nloc; lid++) {
709: NState state = lid_state[lid];
710: if (IS_SELECTED(state)) {
711: /* steal locals */
712: ii = matA_1->i;
713: n = ii[lid + 1] - ii[lid];
714: idx = matA_1->j + ii[lid];
715: for (j = 0; j < n; j++) {
716: PetscInt lidj = idx[j], sgid;
717: NState statej = lid_state[lidj];
718: if (statej == DELETED && (sgid = (PetscInt)PetscRealPart(lid_parent_gid[lidj])) != lid + my0) { /* steal local */
719: lid_parent_gid[lidj] = (PetscScalar)(lid + my0); /* send this if sgid is not local */
720: if (sgid >= my0 && sgid < Iend) { /* I'm stealing this local from a local sgid */
721: PetscInt hav = 0, slid = sgid - my0, gidj = lidj + my0;
722: PetscCDIntNd *pos, *last = NULL;
723: /* looking for local from local so id_llist_2 works */
724: PetscCall(PetscCDGetHeadPos(aggs_2, slid, &pos));
725: while (pos) {
726: PetscInt gid;
727: PetscCall(PetscCDIntNdGetID(pos, &gid));
728: if (gid == gidj) {
729: PetscCheck(last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "last cannot be null");
730: PetscCall(PetscCDRemoveNextNode(aggs_2, slid, last));
731: PetscCall(PetscCDAppendNode(aggs_2, lid, pos));
732: hav = 1;
733: break;
734: } else last = pos;
735: PetscCall(PetscCDGetNextPos(aggs_2, slid, &pos));
736: }
737: if (hav != 1) {
738: PetscCheck(hav, PETSC_COMM_SELF, PETSC_ERR_PLIB, "failed to find adj in 'selected' lists - structurally unsymmetric matrix");
739: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "found node %d times???", (int)hav);
740: }
741: } else { /* I'm stealing this local, owned by a ghost */
742: 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.",
743: ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "");
744: PetscCall(PetscCDAppendID(aggs_2, lid, lidj + my0));
745: }
746: }
747: } /* local neighbors */
748: } else if (state == DELETED /* && lid_cprowID_1 */) {
749: PetscInt sgidold = (PetscInt)PetscRealPart(lid_parent_gid[lid]);
750: /* see if I have a selected ghost neighbor that will steal me */
751: if ((ix = lid_cprowID_1[lid]) != -1) {
752: ii = matB_1->compressedrow.i;
753: n = ii[ix + 1] - ii[ix];
754: idx = matB_1->j + ii[ix];
755: for (j = 0; j < n; j++) {
756: PetscInt cpid = idx[j];
757: NState statej = (NState)PetscRealPart(cpcol_1_state[cpid]);
758: if (IS_SELECTED(statej) && sgidold != (PetscInt)statej) { /* ghost will steal this, remove from my list */
759: lid_parent_gid[lid] = (PetscScalar)statej; /* send who selected */
760: if (sgidold >= my0 && sgidold < Iend) { /* this was mine */
761: PetscInt hav = 0, oldslidj = sgidold - my0;
762: PetscCDIntNd *pos, *last = NULL;
763: /* remove from 'oldslidj' list */
764: PetscCall(PetscCDGetHeadPos(aggs_2, oldslidj, &pos));
765: while (pos) {
766: PetscInt gid;
767: PetscCall(PetscCDIntNdGetID(pos, &gid));
768: if (lid + my0 == gid) {
769: /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove lid from oldslidj list *\/ */
770: PetscCheck(last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "last cannot be null");
771: PetscCall(PetscCDRemoveNextNode(aggs_2, oldslidj, last));
772: /* ghost (PetscScalar)statej will add this later */
773: hav = 1;
774: break;
775: } else last = pos;
776: PetscCall(PetscCDGetNextPos(aggs_2, oldslidj, &pos));
777: }
778: if (hav != 1) {
779: PetscCheck(hav, PETSC_COMM_SELF, PETSC_ERR_PLIB, "failed to find (hav=%d) adj in 'selected' lists - structurally unsymmetric matrix", (int)hav);
780: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "found node %d times???", (int)hav);
781: }
782: } else {
783: /* TODO: ghosts remove this later */
784: }
785: }
786: }
787: }
788: } /* selected/deleted */
789: } /* node loop */
791: if (isMPI) {
792: PetscScalar *cpcol_2_parent, *cpcol_2_gid;
793: Vec tempVec, ghostgids2, ghostparents2;
794: PetscInt cpid, nghost_2;
795: PCGAMGHashTable gid_cpid;
797: PetscCall(VecGetSize(mpimat_2->lvec, &nghost_2));
798: PetscCall(MatCreateVecs(Gmat_2, &tempVec, NULL));
800: /* get 'cpcol_2_parent' */
801: for (kk = 0, j = my0; kk < nloc; kk++, j++) { PetscCall(VecSetValues(tempVec, 1, &j, &lid_parent_gid[kk], INSERT_VALUES)); }
802: PetscCall(VecAssemblyBegin(tempVec));
803: PetscCall(VecAssemblyEnd(tempVec));
804: PetscCall(VecDuplicate(mpimat_2->lvec, &ghostparents2));
805: PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghostparents2, INSERT_VALUES, SCATTER_FORWARD));
806: PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghostparents2, INSERT_VALUES, SCATTER_FORWARD));
807: PetscCall(VecGetArray(ghostparents2, &cpcol_2_parent));
809: /* get 'cpcol_2_gid' */
810: for (kk = 0, j = my0; kk < nloc; kk++, j++) {
811: PetscScalar v = (PetscScalar)j;
812: PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
813: }
814: PetscCall(VecAssemblyBegin(tempVec));
815: PetscCall(VecAssemblyEnd(tempVec));
816: PetscCall(VecDuplicate(mpimat_2->lvec, &ghostgids2));
817: PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghostgids2, INSERT_VALUES, SCATTER_FORWARD));
818: PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghostgids2, INSERT_VALUES, SCATTER_FORWARD));
819: PetscCall(VecGetArray(ghostgids2, &cpcol_2_gid));
820: PetscCall(VecDestroy(&tempVec));
822: /* look for deleted ghosts and add to table */
823: PetscCall(PCGAMGHashTableCreate(2 * nghost_2 + 1, &gid_cpid));
824: for (cpid = 0; cpid < nghost_2; cpid++) {
825: NState state = (NState)PetscRealPart(cpcol_2_state[cpid]);
826: if (state == DELETED) {
827: PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
828: PetscInt sgid_old = (PetscInt)PetscRealPart(cpcol_2_par_orig[cpid]);
829: if (sgid_old == -1 && sgid_new != -1) {
830: PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
831: PetscCall(PCGAMGHashTableAdd(&gid_cpid, gid, cpid));
832: }
833: }
834: }
836: /* look for deleted ghosts and see if they moved - remove it */
837: for (lid = 0; lid < nloc; lid++) {
838: NState state = lid_state[lid];
839: if (IS_SELECTED(state)) {
840: PetscCDIntNd *pos, *last = NULL;
841: /* look for deleted ghosts and see if they moved */
842: PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
843: while (pos) {
844: PetscInt gid;
845: PetscCall(PetscCDIntNdGetID(pos, &gid));
847: if (gid < my0 || gid >= Iend) {
848: PetscCall(PCGAMGHashTableFind(&gid_cpid, gid, &cpid));
849: if (cpid != -1) {
850: /* a moved ghost - */
851: /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove 'flid' from list *\/ */
852: PetscCall(PetscCDRemoveNextNode(aggs_2, lid, last));
853: } else last = pos;
854: } else last = pos;
856: PetscCall(PetscCDGetNextPos(aggs_2, lid, &pos));
857: } /* loop over list of deleted */
858: } /* selected */
859: }
860: PetscCall(PCGAMGHashTableDestroy(&gid_cpid));
862: /* look at ghosts, see if they changed - and it */
863: for (cpid = 0; cpid < nghost_2; cpid++) {
864: PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
865: if (sgid_new >= my0 && sgid_new < Iend) { /* this is mine */
866: PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
867: PetscInt slid_new = sgid_new - my0, hav = 0;
868: PetscCDIntNd *pos;
870: /* search for this gid to see if I have it */
871: PetscCall(PetscCDGetHeadPos(aggs_2, slid_new, &pos));
872: while (pos) {
873: PetscInt gidj;
874: PetscCall(PetscCDIntNdGetID(pos, &gidj));
875: PetscCall(PetscCDGetNextPos(aggs_2, slid_new, &pos));
877: if (gidj == gid) {
878: hav = 1;
879: break;
880: }
881: }
882: if (hav != 1) {
883: /* insert 'flidj' into head of llist */
884: PetscCall(PetscCDAppendID(aggs_2, slid_new, gid));
885: }
886: }
887: }
888: PetscCall(VecRestoreArray(mpimat_1->lvec, &cpcol_1_state));
889: PetscCall(VecRestoreArray(mpimat_2->lvec, &cpcol_2_state));
890: PetscCall(VecRestoreArray(ghostparents2, &cpcol_2_parent));
891: PetscCall(VecRestoreArray(ghostgids2, &cpcol_2_gid));
892: PetscCall(VecDestroy(&ghostgids2));
893: PetscCall(VecDestroy(&ghostparents2));
894: PetscCall(VecDestroy(&ghost_par_orig2));
895: }
896: PetscCall(PetscFree3(lid_state, lid_parent_gid, lid_cprowID_1));
897: PetscFunctionReturn(PETSC_SUCCESS);
898: }
900: /*
901: PCGAMGCoarsen_AGG - supports squaring the graph (deprecated) and new graph for
902: communication of QR data used with HEM and MISk coarsening
904: Input Parameter:
905: . a_pc - this
907: Input/Output Parameter:
908: . a_Gmat1 - graph to coarsen (in), graph off processor edges for QR gather scatter (out)
910: Output Parameter:
911: . agg_lists - list of aggregates
913: */
914: static PetscErrorCode PCGAMGCoarsen_AGG(PC a_pc, Mat *a_Gmat1, PetscCoarsenData **agg_lists)
915: {
916: PC_MG *mg = (PC_MG *)a_pc->data;
917: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
918: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
919: Mat mat, Gmat2, Gmat1 = *a_Gmat1; /* aggressive graph */
920: IS perm;
921: PetscInt Istart, Iend, Ii, nloc, bs, nn;
922: PetscInt *permute, *degree;
923: PetscBool *bIndexSet;
924: PetscReal hashfact;
925: PetscInt iSwapIndex;
926: PetscRandom random;
927: MPI_Comm comm;
929: PetscFunctionBegin;
930: PetscCall(PetscObjectGetComm((PetscObject)Gmat1, &comm));
931: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
932: PetscCall(MatGetLocalSize(Gmat1, &nn, NULL));
933: PetscCall(MatGetBlockSize(Gmat1, &bs));
934: PetscCheck(bs == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "bs %" PetscInt_FMT " must be 1", bs);
935: nloc = nn / bs;
936: /* get MIS aggs - randomize */
937: PetscCall(PetscMalloc2(nloc, &permute, nloc, °ree));
938: PetscCall(PetscCalloc1(nloc, &bIndexSet));
939: for (Ii = 0; Ii < nloc; Ii++) permute[Ii] = Ii;
940: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &random));
941: PetscCall(MatGetOwnershipRange(Gmat1, &Istart, &Iend));
942: for (Ii = 0; Ii < nloc; Ii++) {
943: PetscInt nc;
944: PetscCall(MatGetRow(Gmat1, Istart + Ii, &nc, NULL, NULL));
945: degree[Ii] = nc;
946: PetscCall(MatRestoreRow(Gmat1, Istart + Ii, &nc, NULL, NULL));
947: }
948: for (Ii = 0; Ii < nloc; Ii++) {
949: PetscCall(PetscRandomGetValueReal(random, &hashfact));
950: iSwapIndex = (PetscInt)(hashfact * nloc) % nloc;
951: if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) {
952: PetscInt iTemp = permute[iSwapIndex];
953: permute[iSwapIndex] = permute[Ii];
954: permute[Ii] = iTemp;
955: iTemp = degree[iSwapIndex];
956: degree[iSwapIndex] = degree[Ii];
957: degree[Ii] = iTemp;
958: bIndexSet[iSwapIndex] = PETSC_TRUE;
959: }
960: }
961: // apply minimum degree ordering -- NEW
962: if (pc_gamg_agg->use_minimum_degree_ordering) { PetscCall(PetscSortIntWithArray(nloc, degree, permute)); }
963: PetscCall(PetscFree(bIndexSet));
964: PetscCall(PetscRandomDestroy(&random));
965: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_USE_POINTER, &perm));
966: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0));
967: // square graph
968: if (pc_gamg->current_level < pc_gamg_agg->aggressive_coarsening_levels && pc_gamg_agg->use_aggressive_square_graph) {
969: PetscCall(PCGAMGSquareGraph_GAMG(a_pc, Gmat1, &Gmat2));
970: } else Gmat2 = Gmat1;
971: // switch to old MIS-1 for square graph
972: if (pc_gamg->current_level < pc_gamg_agg->aggressive_coarsening_levels) {
973: if (!pc_gamg_agg->use_aggressive_square_graph) PetscCall(MatCoarsenMISKSetDistance(pc_gamg_agg->crs, pc_gamg_agg->aggressive_mis_k)); // hardwire to MIS-2
974: else PetscCall(MatCoarsenSetType(pc_gamg_agg->crs, MATCOARSENMIS)); // old MIS -- side effect
975: } else if (pc_gamg_agg->use_aggressive_square_graph && pc_gamg_agg->aggressive_coarsening_levels > 0) { // we reset the MIS
976: const char *prefix;
977: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)a_pc, &prefix));
978: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc_gamg_agg->crs, prefix));
979: PetscCall(MatCoarsenSetFromOptions(pc_gamg_agg->crs)); // get the default back on non-aggressive levels when square graph switched to old MIS
980: }
981: PetscCall(MatCoarsenSetAdjacency(pc_gamg_agg->crs, Gmat2));
982: PetscCall(MatCoarsenSetStrictAggs(pc_gamg_agg->crs, PETSC_TRUE));
983: PetscCall(MatCoarsenSetGreedyOrdering(pc_gamg_agg->crs, perm));
984: PetscCall(MatCoarsenApply(pc_gamg_agg->crs));
985: PetscCall(MatCoarsenViewFromOptions(pc_gamg_agg->crs, NULL, "-mat_coarsen_view"));
986: PetscCall(MatCoarsenGetData(pc_gamg_agg->crs, agg_lists)); /* output */
987: PetscCall(MatCoarsenDestroy(&pc_gamg_agg->crs));
989: PetscCall(ISDestroy(&perm));
990: PetscCall(PetscFree2(permute, degree));
991: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0));
993: if (Gmat2 != Gmat1) {
994: PetscCoarsenData *llist = *agg_lists;
995: PetscCall(fixAggregatesWithSquare(a_pc, Gmat2, Gmat1, *agg_lists));
996: PetscCall(MatDestroy(&Gmat1));
997: *a_Gmat1 = Gmat2; /* output */
998: PetscCall(PetscCDGetMat(llist, &mat));
999: PetscCheck(!mat, comm, PETSC_ERR_ARG_WRONG, "Unexpected auxiliary matrix with squared graph");
1000: } else {
1001: PetscCoarsenData *llist = *agg_lists;
1002: /* see if we have a matrix that takes precedence (returned from MatCoarsenApply) */
1003: PetscCall(PetscCDGetMat(llist, &mat));
1004: if (mat) {
1005: PetscCall(MatDestroy(a_Gmat1));
1006: *a_Gmat1 = mat; /* output */
1007: }
1008: }
1009: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
1010: PetscFunctionReturn(PETSC_SUCCESS);
1011: }
1013: /*
1014: PCGAMGProlongator_AGG
1016: Input Parameter:
1017: . pc - this
1018: . Amat - matrix on this fine level
1019: . Graph - used to get ghost data for nodes in
1020: . agg_lists - list of aggregates
1021: Output Parameter:
1022: . a_P_out - prolongation operator to the next level
1023: */
1024: static PetscErrorCode PCGAMGProlongator_AGG(PC pc, Mat Amat, Mat Gmat, PetscCoarsenData *agg_lists, Mat *a_P_out)
1025: {
1026: PC_MG *mg = (PC_MG *)pc->data;
1027: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1028: const PetscInt col_bs = pc_gamg->data_cell_cols;
1029: PetscInt Istart, Iend, nloc, ii, jj, kk, my0, nLocalSelected, bs;
1030: Mat Prol;
1031: PetscMPIInt size;
1032: MPI_Comm comm;
1033: PetscReal *data_w_ghost;
1034: PetscInt myCrs0, nbnodes = 0, *flid_fgid;
1035: MatType mtype;
1037: PetscFunctionBegin;
1038: PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
1039: PetscCheck(col_bs >= 1, comm, PETSC_ERR_PLIB, "Column bs cannot be less than 1");
1040: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1041: PetscCallMPI(MPI_Comm_size(comm, &size));
1042: PetscCall(MatGetOwnershipRange(Amat, &Istart, &Iend));
1043: PetscCall(MatGetBlockSize(Amat, &bs));
1044: nloc = (Iend - Istart) / bs;
1045: my0 = Istart / bs;
1046: 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);
1048: /* get 'nLocalSelected' */
1049: for (ii = 0, nLocalSelected = 0; ii < nloc; ii++) {
1050: PetscBool ise;
1051: /* filter out singletons 0 or 1? */
1052: PetscCall(PetscCDEmptyAt(agg_lists, ii, &ise));
1053: if (!ise) nLocalSelected++;
1054: }
1056: /* create prolongator, create P matrix */
1057: PetscCall(MatGetType(Amat, &mtype));
1058: PetscCall(MatCreate(comm, &Prol));
1059: PetscCall(MatSetSizes(Prol, nloc * bs, nLocalSelected * col_bs, PETSC_DETERMINE, PETSC_DETERMINE));
1060: PetscCall(MatSetBlockSizes(Prol, bs, col_bs));
1061: PetscCall(MatSetType(Prol, mtype));
1062: #if PetscDefined(HAVE_DEVICE)
1063: PetscBool flg;
1064: PetscCall(MatBoundToCPU(Amat, &flg));
1065: PetscCall(MatBindToCPU(Prol, flg));
1066: if (flg) PetscCall(MatSetBindingPropagates(Prol, PETSC_TRUE));
1067: #endif
1068: PetscCall(MatSeqAIJSetPreallocation(Prol, col_bs, NULL));
1069: PetscCall(MatMPIAIJSetPreallocation(Prol, col_bs, NULL, col_bs, NULL));
1071: /* can get all points "removed" */
1072: PetscCall(MatGetSize(Prol, &kk, &ii));
1073: if (!ii) {
1074: PetscCall(PetscInfo(pc, "%s: No selected points on coarse grid\n", ((PetscObject)pc)->prefix));
1075: PetscCall(MatDestroy(&Prol));
1076: *a_P_out = NULL; /* out */
1077: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1078: PetscFunctionReturn(PETSC_SUCCESS);
1079: }
1080: PetscCall(PetscInfo(pc, "%s: New grid %" PetscInt_FMT " nodes\n", ((PetscObject)pc)->prefix, ii / col_bs));
1081: PetscCall(MatGetOwnershipRangeColumn(Prol, &myCrs0, &kk));
1083: 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);
1084: myCrs0 = myCrs0 / col_bs;
1085: 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);
1087: /* create global vector of data in 'data_w_ghost' */
1088: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0));
1089: if (size > 1) { /* get ghost null space data */
1090: PetscReal *tmp_gdata, *tmp_ldata, *tp2;
1091: PetscCall(PetscMalloc1(nloc, &tmp_ldata));
1092: for (jj = 0; jj < col_bs; jj++) {
1093: for (kk = 0; kk < bs; kk++) {
1094: PetscInt ii, stride;
1095: const PetscReal *tp = pc_gamg->data + jj * bs * nloc + kk;
1096: for (ii = 0; ii < nloc; ii++, tp += bs) tmp_ldata[ii] = *tp;
1098: PetscCall(PCGAMGGetDataWithGhosts(Gmat, 1, tmp_ldata, &stride, &tmp_gdata));
1100: if (!jj && !kk) { /* now I know how many total nodes - allocate TODO: move below and do in one 'col_bs' call */
1101: PetscCall(PetscMalloc1(stride * bs * col_bs, &data_w_ghost));
1102: nbnodes = bs * stride;
1103: }
1104: tp2 = data_w_ghost + jj * bs * stride + kk;
1105: for (ii = 0; ii < stride; ii++, tp2 += bs) *tp2 = tmp_gdata[ii];
1106: PetscCall(PetscFree(tmp_gdata));
1107: }
1108: }
1109: PetscCall(PetscFree(tmp_ldata));
1110: } else {
1111: nbnodes = bs * nloc;
1112: data_w_ghost = (PetscReal *)pc_gamg->data;
1113: }
1115: /* get 'flid_fgid' TODO - move up to get 'stride' and do get null space data above in one step (jj loop) */
1116: if (size > 1) {
1117: PetscReal *fid_glid_loc, *fiddata;
1118: PetscInt stride;
1120: PetscCall(PetscMalloc1(nloc, &fid_glid_loc));
1121: for (kk = 0; kk < nloc; kk++) fid_glid_loc[kk] = (PetscReal)(my0 + kk);
1122: PetscCall(PCGAMGGetDataWithGhosts(Gmat, 1, fid_glid_loc, &stride, &fiddata));
1123: PetscCall(PetscMalloc1(stride, &flid_fgid)); /* copy real data to in */
1124: for (kk = 0; kk < stride; kk++) flid_fgid[kk] = (PetscInt)fiddata[kk];
1125: PetscCall(PetscFree(fiddata));
1127: PetscCheck(stride == nbnodes / bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "stride %" PetscInt_FMT " != nbnodes %" PetscInt_FMT "/bs %" PetscInt_FMT, stride, nbnodes, bs);
1128: PetscCall(PetscFree(fid_glid_loc));
1129: } else {
1130: PetscCall(PetscMalloc1(nloc, &flid_fgid));
1131: for (kk = 0; kk < nloc; kk++) flid_fgid[kk] = my0 + kk;
1132: }
1133: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0));
1134: /* get P0 */
1135: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0));
1136: {
1137: PetscReal *data_out = NULL;
1138: PetscCall(formProl0(agg_lists, bs, col_bs, myCrs0, nbnodes, data_w_ghost, flid_fgid, &data_out, Prol));
1139: PetscCall(PetscFree(pc_gamg->data));
1141: pc_gamg->data = data_out;
1142: pc_gamg->data_cell_rows = col_bs;
1143: pc_gamg->data_sz = col_bs * col_bs * nLocalSelected;
1144: }
1145: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0));
1146: if (size > 1) PetscCall(PetscFree(data_w_ghost));
1147: PetscCall(PetscFree(flid_fgid));
1149: *a_P_out = Prol; /* out */
1151: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1152: PetscFunctionReturn(PETSC_SUCCESS);
1153: }
1155: /*
1156: PCGAMGOptProlongator_AGG
1158: Input Parameter:
1159: . pc - this
1160: . Amat - matrix on this fine level
1161: In/Output Parameter:
1162: . a_P - prolongation operator to the next level
1163: */
1164: static PetscErrorCode PCGAMGOptProlongator_AGG(PC pc, Mat Amat, Mat *a_P)
1165: {
1166: PC_MG *mg = (PC_MG *)pc->data;
1167: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1168: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
1169: PetscInt jj;
1170: Mat Prol = *a_P;
1171: MPI_Comm comm;
1172: KSP eksp;
1173: Vec bb, xx;
1174: PC epc;
1175: PetscReal alpha, emax, emin;
1177: PetscFunctionBegin;
1178: PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
1179: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0));
1181: /* compute maximum singular value of operator to be used in smoother */
1182: if (0 < pc_gamg_agg->nsmooths) {
1183: /* get eigen estimates */
1184: if (pc_gamg->emax > 0) {
1185: emin = pc_gamg->emin;
1186: emax = pc_gamg->emax;
1187: } else {
1188: const char *prefix;
1190: PetscCall(MatCreateVecs(Amat, &bb, NULL));
1191: PetscCall(MatCreateVecs(Amat, &xx, NULL));
1192: PetscCall(KSPSetNoisy_Private(bb));
1194: PetscCall(KSPCreate(comm, &eksp));
1195: PetscCall(KSPSetNestLevel(eksp, pc->kspnestlevel));
1196: PetscCall(PCGetOptionsPrefix(pc, &prefix));
1197: PetscCall(KSPSetOptionsPrefix(eksp, prefix));
1198: PetscCall(KSPAppendOptionsPrefix(eksp, "pc_gamg_esteig_"));
1199: {
1200: PetscBool isset, sflg;
1201: PetscCall(MatIsSPDKnown(Amat, &isset, &sflg));
1202: if (isset && sflg) PetscCall(KSPSetType(eksp, KSPCG));
1203: }
1204: PetscCall(KSPSetErrorIfNotConverged(eksp, pc->erroriffailure));
1205: PetscCall(KSPSetNormType(eksp, KSP_NORM_NONE));
1207: PetscCall(KSPSetInitialGuessNonzero(eksp, PETSC_FALSE));
1208: PetscCall(KSPSetOperators(eksp, Amat, Amat));
1210: PetscCall(KSPGetPC(eksp, &epc));
1211: PetscCall(PCSetType(epc, PCJACOBI)); /* smoother in smoothed agg. */
1213: 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
1215: PetscCall(KSPSetFromOptions(eksp));
1216: PetscCall(KSPSetComputeSingularValues(eksp, PETSC_TRUE));
1217: PetscCall(KSPSolve(eksp, bb, xx));
1218: PetscCall(KSPCheckSolve(eksp, pc, xx));
1220: PetscCall(KSPComputeExtremeSingularValues(eksp, &emax, &emin));
1221: PetscCall(PetscInfo(pc, "%s: Smooth P0: max eigen=%e min=%e PC=%s\n", ((PetscObject)pc)->prefix, (double)emax, (double)emin, PCJACOBI));
1222: PetscCall(VecDestroy(&xx));
1223: PetscCall(VecDestroy(&bb));
1224: PetscCall(KSPDestroy(&eksp));
1225: }
1226: if (pc_gamg->use_sa_esteig) {
1227: mg->min_eigen_DinvA[pc_gamg->current_level] = emin;
1228: mg->max_eigen_DinvA[pc_gamg->current_level] = emax;
1229: 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));
1230: } else {
1231: mg->min_eigen_DinvA[pc_gamg->current_level] = 0;
1232: mg->max_eigen_DinvA[pc_gamg->current_level] = 0;
1233: }
1234: } else {
1235: mg->min_eigen_DinvA[pc_gamg->current_level] = 0;
1236: mg->max_eigen_DinvA[pc_gamg->current_level] = 0;
1237: }
1239: /* smooth P0 */
1240: for (jj = 0; jj < pc_gamg_agg->nsmooths; jj++) {
1241: Mat tMat;
1242: Vec diag;
1244: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0));
1246: /* smooth P1 := (I - omega/lam D^{-1}A)P0 */
1247: PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0));
1248: PetscCall(MatMatMult(Amat, Prol, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &tMat));
1249: PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0));
1250: PetscCall(MatProductClear(tMat));
1251: PetscCall(MatCreateVecs(Amat, &diag, NULL));
1252: PetscCall(MatGetDiagonal(Amat, diag)); /* effectively PCJACOBI */
1253: PetscCall(VecReciprocal(diag));
1254: PetscCall(MatDiagonalScale(tMat, diag, NULL));
1255: PetscCall(VecDestroy(&diag));
1257: /* TODO: Set a PCFailedReason and exit the building of the AMG preconditioner */
1258: PetscCheck(emax != 0.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Computed maximum singular value as zero");
1259: /* TODO: Document the 1.4 and don't hardwire it in this routine */
1260: alpha = -1.4 / emax;
1262: PetscCall(MatAYPX(tMat, alpha, Prol, SUBSET_NONZERO_PATTERN));
1263: PetscCall(MatDestroy(&Prol));
1264: Prol = tMat;
1265: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0));
1266: }
1267: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0));
1268: *a_P = Prol;
1269: PetscFunctionReturn(PETSC_SUCCESS);
1270: }
1272: /*
1273: PCCreateGAMG_AGG
1275: Input Parameter:
1276: . pc -
1277: */
1278: PetscErrorCode PCCreateGAMG_AGG(PC pc)
1279: {
1280: PC_MG *mg = (PC_MG *)pc->data;
1281: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1282: PC_GAMG_AGG *pc_gamg_agg;
1284: PetscFunctionBegin;
1285: /* create sub context for SA */
1286: PetscCall(PetscNew(&pc_gamg_agg));
1287: pc_gamg->subctx = pc_gamg_agg;
1289: pc_gamg->ops->setfromoptions = PCSetFromOptions_GAMG_AGG;
1290: pc_gamg->ops->destroy = PCDestroy_GAMG_AGG;
1291: /* reset does not do anything; setup not virtual */
1293: /* set internal function pointers */
1294: pc_gamg->ops->creategraph = PCGAMGCreateGraph_AGG;
1295: pc_gamg->ops->coarsen = PCGAMGCoarsen_AGG;
1296: pc_gamg->ops->prolongator = PCGAMGProlongator_AGG;
1297: pc_gamg->ops->optprolongator = PCGAMGOptProlongator_AGG;
1298: pc_gamg->ops->createdefaultdata = PCSetData_AGG;
1299: pc_gamg->ops->view = PCView_GAMG_AGG;
1301: pc_gamg_agg->nsmooths = 1;
1302: pc_gamg_agg->aggressive_coarsening_levels = 1;
1303: pc_gamg_agg->use_aggressive_square_graph = PETSC_FALSE;
1304: pc_gamg_agg->use_minimum_degree_ordering = PETSC_FALSE;
1305: pc_gamg_agg->aggressive_mis_k = 2;
1307: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", PCGAMGSetNSmooths_AGG));
1308: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", PCGAMGSetAggressiveLevels_AGG));
1309: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveSquareGraph_C", PCGAMGSetAggressiveSquareGraph_AGG));
1310: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetMinDegreeOrdering_C", PCGAMGMISkSetMinDegreeOrdering_AGG));
1311: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetAggressive_C", PCGAMGMISkSetAggressive_AGG));
1312: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_AGG));
1313: PetscFunctionReturn(PETSC_SUCCESS);
1314: }