Actual source code: classical.c
1: #include <../src/ksp/pc/impls/gamg/gamg.h>
2: #include <petscsf.h>
4: static PetscFunctionList PCGAMGClassicalProlongatorList = NULL;
5: static PetscBool PCGAMGClassicalPackageInitialized = PETSC_FALSE;
7: typedef struct {
8: PetscReal interp_threshold; /* interpolation threshold */
9: char prolongtype[256];
10: PetscInt nsmooths; /* number of jacobi smoothings on the prolongator */
11: } PC_GAMG_Classical;
13: /*@
14: PCGAMGClassicalSetType - Sets the type of classical interpolation to use with `PCGAMG`
16: Collective
18: Input Parameters:
19: + pc - the preconditioner context
20: - type - the interpolation to use, see `PCGAMGClassicalType()`
22: Options Database Key:
23: . -pc_gamg_classical_type <direct,standard> - set type of classical AMG prolongation
25: Level: intermediate
27: .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGClassicalType`, `PCGAMGClassicalGetType()`
28: @*/
29: PetscErrorCode PCGAMGClassicalSetType(PC pc, PCGAMGClassicalType type)
30: {
31: PetscFunctionBegin;
33: PetscTryMethod(pc, "PCGAMGClassicalSetType_C", (PC, PCGAMGClassicalType), (pc, type));
34: PetscFunctionReturn(PETSC_SUCCESS);
35: }
37: /*@
38: PCGAMGClassicalGetType - Gets the type of classical interpolation to use with `PCGAMG`
40: Collective
42: Input Parameter:
43: . pc - the preconditioner context
45: Output Parameter:
46: . type - the type used, see `PCGAMGClassicalType()`
48: Level: intermediate
50: .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGClassicalType`, `PCGAMGClassicalSetType()`
51: @*/
52: PetscErrorCode PCGAMGClassicalGetType(PC pc, PCGAMGClassicalType *type)
53: {
54: PetscFunctionBegin;
56: PetscUseMethod(pc, "PCGAMGClassicalGetType_C", (PC, PCGAMGClassicalType *), (pc, type));
57: PetscFunctionReturn(PETSC_SUCCESS);
58: }
60: static PetscErrorCode PCGAMGClassicalSetType_GAMG(PC pc, PCGAMGClassicalType type)
61: {
62: PC_MG *mg = (PC_MG *)pc->data;
63: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
64: PC_GAMG_Classical *cls = (PC_GAMG_Classical *)pc_gamg->subctx;
66: PetscFunctionBegin;
67: PetscCall(PetscStrncpy(cls->prolongtype, type, sizeof(cls->prolongtype)));
68: PetscFunctionReturn(PETSC_SUCCESS);
69: }
71: static PetscErrorCode PCGAMGClassicalGetType_GAMG(PC pc, PCGAMGClassicalType *type)
72: {
73: PC_MG *mg = (PC_MG *)pc->data;
74: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
75: PC_GAMG_Classical *cls = (PC_GAMG_Classical *)pc_gamg->subctx;
77: PetscFunctionBegin;
78: *type = cls->prolongtype;
79: PetscFunctionReturn(PETSC_SUCCESS);
80: }
82: static PetscErrorCode PCGAMGCreateGraph_Classical(PC pc, Mat A, Mat *G)
83: {
84: PetscInt s, f, n, idx, lidx, gidx;
85: PetscInt r, c, ncols;
86: const PetscInt *rcol;
87: const PetscScalar *rval;
88: PetscInt *gcol;
89: PetscScalar *gval;
90: PetscReal rmax;
91: PetscInt cmax = 0;
92: PC_MG *mg = (PC_MG *)pc->data;
93: PC_GAMG *gamg = (PC_GAMG *)mg->innerctx;
94: PetscInt *gsparse, *lsparse;
95: PetscScalar *Amax;
96: MatType mtype;
98: PetscFunctionBegin;
99: PetscCall(MatGetOwnershipRange(A, &s, &f));
100: n = f - s;
101: PetscCall(PetscMalloc3(n, &lsparse, n, &gsparse, n, &Amax));
103: for (r = 0; r < n; r++) {
104: lsparse[r] = 0;
105: gsparse[r] = 0;
106: }
108: for (r = s; r < f; r++) {
109: /* determine the maximum off-diagonal in each row */
110: rmax = 0.;
111: PetscCall(MatGetRow(A, r, &ncols, &rcol, &rval));
112: for (c = 0; c < ncols; c++) {
113: if (PetscRealPart(-rval[c]) > rmax && rcol[c] != r) rmax = PetscRealPart(-rval[c]);
114: }
115: Amax[r - s] = rmax;
116: if (ncols > cmax) cmax = ncols;
117: lidx = 0;
118: gidx = 0;
119: /* create the local and global sparsity patterns */
120: for (c = 0; c < ncols; c++) {
121: if (PetscRealPart(-rval[c]) > gamg->threshold[0] * PetscRealPart(Amax[r - s]) || rcol[c] == r) {
122: if (rcol[c] < f && rcol[c] >= s) {
123: lidx++;
124: } else {
125: gidx++;
126: }
127: }
128: }
129: PetscCall(MatRestoreRow(A, r, &ncols, &rcol, &rval));
130: lsparse[r - s] = lidx;
131: gsparse[r - s] = gidx;
132: }
133: PetscCall(PetscMalloc2(cmax, &gval, cmax, &gcol));
135: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), G));
136: PetscCall(MatGetType(A, &mtype));
137: PetscCall(MatSetType(*G, mtype));
138: PetscCall(MatSetSizes(*G, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
139: PetscCall(MatMPIAIJSetPreallocation(*G, 0, lsparse, 0, gsparse));
140: PetscCall(MatSeqAIJSetPreallocation(*G, 0, lsparse));
141: for (r = s; r < f; r++) {
142: PetscCall(MatGetRow(A, r, &ncols, &rcol, &rval));
143: idx = 0;
144: for (c = 0; c < ncols; c++) {
145: /* classical strength of connection */
146: if (PetscRealPart(-rval[c]) > gamg->threshold[0] * PetscRealPart(Amax[r - s]) || rcol[c] == r) {
147: gcol[idx] = rcol[c];
148: gval[idx] = rval[c];
149: idx++;
150: }
151: }
152: PetscCall(MatSetValues(*G, 1, &r, idx, gcol, gval, INSERT_VALUES));
153: PetscCall(MatRestoreRow(A, r, &ncols, &rcol, &rval));
154: }
155: PetscCall(MatAssemblyBegin(*G, MAT_FINAL_ASSEMBLY));
156: PetscCall(MatAssemblyEnd(*G, MAT_FINAL_ASSEMBLY));
158: PetscCall(PetscFree2(gval, gcol));
159: PetscCall(PetscFree3(lsparse, gsparse, Amax));
160: PetscFunctionReturn(PETSC_SUCCESS);
161: }
163: static PetscErrorCode PCGAMGCoarsen_Classical(PC pc, Mat *G, PetscCoarsenData **agg_lists)
164: {
165: MatCoarsen crs;
166: MPI_Comm fcomm = ((PetscObject)pc)->comm;
167: const char *prefix;
169: PetscFunctionBegin;
170: PetscCheck(G, fcomm, PETSC_ERR_ARG_WRONGSTATE, "Must set Graph in PC in PCGAMG before coarsening");
172: PetscCall(MatCoarsenCreate(fcomm, &crs));
173: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix));
174: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)crs, prefix));
175: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)crs, "pc_gamg_"));
176: PetscCall(MatCoarsenSetFromOptions(crs));
177: PetscCall(MatCoarsenSetAdjacency(crs, *G));
178: PetscCall(MatCoarsenSetStrictAggs(crs, PETSC_TRUE));
179: PetscCall(MatCoarsenApply(crs));
180: PetscCall(MatCoarsenGetData(crs, agg_lists));
181: PetscCall(MatCoarsenDestroy(&crs));
182: PetscFunctionReturn(PETSC_SUCCESS);
183: }
185: static PetscErrorCode PCGAMGProlongator_Classical_Direct(PC pc, Mat A, PetscCoarsenData *agg_lists, Mat *P)
186: {
187: PC_MG *mg = (PC_MG *)pc->data;
188: PC_GAMG *gamg = (PC_GAMG *)mg->innerctx;
189: PetscBool iscoarse, isMPIAIJ, isSEQAIJ;
190: PetscInt fn, cn, fs, fe, cs, ce, i, j, ncols, col, row_f, row_c, cmax = 0, idx, noff;
191: PetscInt *lcid, *gcid, *lsparse, *gsparse, *colmap, *pcols;
192: const PetscInt *rcol;
193: PetscReal *Amax_pos, *Amax_neg;
194: PetscScalar g_pos, g_neg, a_pos, a_neg, diag, invdiag, alpha, beta, pij;
195: PetscScalar *pvals;
196: const PetscScalar *rval;
197: Mat lA, gA = NULL;
198: MatType mtype;
199: Vec C, lvec;
200: PetscLayout clayout;
201: PetscSF sf;
202: Mat_MPIAIJ *mpiaij;
204: PetscFunctionBegin;
205: PetscCall(MatGetOwnershipRange(A, &fs, &fe));
206: fn = fe - fs;
207: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPIAIJ, &isMPIAIJ));
208: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQAIJ, &isSEQAIJ));
209: PetscCheck(isMPIAIJ || isSEQAIJ, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Classical AMG requires MPIAIJ matrix");
210: if (isMPIAIJ) {
211: mpiaij = (Mat_MPIAIJ *)A->data;
212: lA = mpiaij->A;
213: gA = mpiaij->B;
214: lvec = mpiaij->lvec;
215: PetscCall(VecGetSize(lvec, &noff));
216: colmap = mpiaij->garray;
217: PetscCall(MatGetLayouts(A, NULL, &clayout));
218: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)A), &sf));
219: PetscCall(PetscSFSetGraphLayout(sf, clayout, noff, NULL, PETSC_COPY_VALUES, colmap));
220: PetscCall(PetscMalloc1(noff, &gcid));
221: } else {
222: lA = A;
223: }
224: PetscCall(PetscMalloc5(fn, &lsparse, fn, &gsparse, fn, &lcid, fn, &Amax_pos, fn, &Amax_neg));
226: /* count the number of coarse unknowns */
227: cn = 0;
228: for (i = 0; i < fn; i++) {
229: /* filter out singletons */
230: PetscCall(PetscCDIsEmptyAt(agg_lists, i, &iscoarse));
231: lcid[i] = -1;
232: if (!iscoarse) cn++;
233: }
235: /* create the coarse vector */
236: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)A), cn, PETSC_DECIDE, &C));
237: PetscCall(VecGetOwnershipRange(C, &cs, &ce));
239: cn = 0;
240: for (i = 0; i < fn; i++) {
241: PetscCall(PetscCDIsEmptyAt(agg_lists, i, &iscoarse));
242: if (!iscoarse) {
243: lcid[i] = cs + cn;
244: cn++;
245: } else {
246: lcid[i] = -1;
247: }
248: }
250: if (gA) {
251: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, lcid, gcid, MPI_REPLACE));
252: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, lcid, gcid, MPI_REPLACE));
253: }
255: /* determine the largest off-diagonal entries in each row */
256: for (i = fs; i < fe; i++) {
257: Amax_pos[i - fs] = 0.;
258: Amax_neg[i - fs] = 0.;
259: PetscCall(MatGetRow(A, i, &ncols, &rcol, &rval));
260: for (j = 0; j < ncols; j++) {
261: if ((PetscRealPart(-rval[j]) > Amax_neg[i - fs]) && i != rcol[j]) Amax_neg[i - fs] = PetscAbsScalar(rval[j]);
262: if ((PetscRealPart(rval[j]) > Amax_pos[i - fs]) && i != rcol[j]) Amax_pos[i - fs] = PetscAbsScalar(rval[j]);
263: }
264: if (ncols > cmax) cmax = ncols;
265: PetscCall(MatRestoreRow(A, i, &ncols, &rcol, &rval));
266: }
267: PetscCall(PetscMalloc2(cmax, &pcols, cmax, &pvals));
268: PetscCall(VecDestroy(&C));
270: /* count the on and off processor sparsity patterns for the prolongator */
271: for (i = 0; i < fn; i++) {
272: /* on */
273: lsparse[i] = 0;
274: gsparse[i] = 0;
275: if (lcid[i] >= 0) {
276: lsparse[i] = 1;
277: gsparse[i] = 0;
278: } else {
279: PetscCall(MatGetRow(lA, i, &ncols, &rcol, &rval));
280: for (j = 0; j < ncols; j++) {
281: col = rcol[j];
282: if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0] * Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0] * Amax_neg[i])) lsparse[i] += 1;
283: }
284: PetscCall(MatRestoreRow(lA, i, &ncols, &rcol, &rval));
285: /* off */
286: if (gA) {
287: PetscCall(MatGetRow(gA, i, &ncols, &rcol, &rval));
288: for (j = 0; j < ncols; j++) {
289: col = rcol[j];
290: if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0] * Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0] * Amax_neg[i])) gsparse[i] += 1;
291: }
292: PetscCall(MatRestoreRow(gA, i, &ncols, &rcol, &rval));
293: }
294: }
295: }
297: /* preallocate and create the prolongator */
298: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), P));
299: PetscCall(MatGetType(A, &mtype));
300: PetscCall(MatSetType(*P, mtype));
301: PetscCall(MatSetSizes(*P, fn, cn, PETSC_DETERMINE, PETSC_DETERMINE));
302: PetscCall(MatMPIAIJSetPreallocation(*P, 0, lsparse, 0, gsparse));
303: PetscCall(MatSeqAIJSetPreallocation(*P, 0, lsparse));
305: /* loop over local fine nodes -- get the diagonal, the sum of positive and negative strong and weak weights, and set up the row */
306: for (i = 0; i < fn; i++) {
307: /* determine on or off */
308: row_f = i + fs;
309: row_c = lcid[i];
310: if (row_c >= 0) {
311: pij = 1.;
312: PetscCall(MatSetValues(*P, 1, &row_f, 1, &row_c, &pij, INSERT_VALUES));
313: } else {
314: g_pos = 0.;
315: g_neg = 0.;
316: a_pos = 0.;
317: a_neg = 0.;
318: diag = 0.;
320: /* local connections */
321: PetscCall(MatGetRow(lA, i, &ncols, &rcol, &rval));
322: for (j = 0; j < ncols; j++) {
323: col = rcol[j];
324: if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0] * Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0] * Amax_neg[i])) {
325: if (PetscRealPart(rval[j]) > 0.) {
326: g_pos += rval[j];
327: } else {
328: g_neg += rval[j];
329: }
330: }
331: if (col != i) {
332: if (PetscRealPart(rval[j]) > 0.) {
333: a_pos += rval[j];
334: } else {
335: a_neg += rval[j];
336: }
337: } else {
338: diag = rval[j];
339: }
340: }
341: PetscCall(MatRestoreRow(lA, i, &ncols, &rcol, &rval));
343: /* ghosted connections */
344: if (gA) {
345: PetscCall(MatGetRow(gA, i, &ncols, &rcol, &rval));
346: for (j = 0; j < ncols; j++) {
347: col = rcol[j];
348: if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0] * Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0] * Amax_neg[i])) {
349: if (PetscRealPart(rval[j]) > 0.) {
350: g_pos += rval[j];
351: } else {
352: g_neg += rval[j];
353: }
354: }
355: if (PetscRealPart(rval[j]) > 0.) {
356: a_pos += rval[j];
357: } else {
358: a_neg += rval[j];
359: }
360: }
361: PetscCall(MatRestoreRow(gA, i, &ncols, &rcol, &rval));
362: }
364: if (g_neg == 0.) {
365: alpha = 0.;
366: } else {
367: alpha = -a_neg / g_neg;
368: }
370: if (g_pos == 0.) {
371: diag += a_pos;
372: beta = 0.;
373: } else {
374: beta = -a_pos / g_pos;
375: }
376: if (diag == 0.) {
377: invdiag = 0.;
378: } else invdiag = 1. / diag;
379: /* on */
380: PetscCall(MatGetRow(lA, i, &ncols, &rcol, &rval));
381: idx = 0;
382: for (j = 0; j < ncols; j++) {
383: col = rcol[j];
384: if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0] * Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0] * Amax_neg[i])) {
385: row_f = i + fs;
386: row_c = lcid[col];
387: /* set the values for on-processor ones */
388: if (PetscRealPart(rval[j]) < 0.) {
389: pij = rval[j] * alpha * invdiag;
390: } else {
391: pij = rval[j] * beta * invdiag;
392: }
393: if (PetscAbsScalar(pij) != 0.) {
394: pvals[idx] = pij;
395: pcols[idx] = row_c;
396: idx++;
397: }
398: }
399: }
400: PetscCall(MatRestoreRow(lA, i, &ncols, &rcol, &rval));
401: /* off */
402: if (gA) {
403: PetscCall(MatGetRow(gA, i, &ncols, &rcol, &rval));
404: for (j = 0; j < ncols; j++) {
405: col = rcol[j];
406: if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0] * Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0] * Amax_neg[i])) {
407: row_f = i + fs;
408: row_c = gcid[col];
409: /* set the values for on-processor ones */
410: if (PetscRealPart(rval[j]) < 0.) {
411: pij = rval[j] * alpha * invdiag;
412: } else {
413: pij = rval[j] * beta * invdiag;
414: }
415: if (PetscAbsScalar(pij) != 0.) {
416: pvals[idx] = pij;
417: pcols[idx] = row_c;
418: idx++;
419: }
420: }
421: }
422: PetscCall(MatRestoreRow(gA, i, &ncols, &rcol, &rval));
423: }
424: PetscCall(MatSetValues(*P, 1, &row_f, idx, pcols, pvals, INSERT_VALUES));
425: }
426: }
428: PetscCall(MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY));
429: PetscCall(MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY));
431: PetscCall(PetscFree5(lsparse, gsparse, lcid, Amax_pos, Amax_neg));
433: PetscCall(PetscFree2(pcols, pvals));
434: if (gA) {
435: PetscCall(PetscSFDestroy(&sf));
436: PetscCall(PetscFree(gcid));
437: }
438: PetscFunctionReturn(PETSC_SUCCESS);
439: }
441: static PetscErrorCode PCGAMGTruncateProlongator_Private(PC pc, Mat *P)
442: {
443: PetscInt j, i, ps, pf, pn, pcs, pcf, pcn, idx, cmax;
444: const PetscScalar *pval;
445: const PetscInt *pcol;
446: PetscScalar *pnval;
447: PetscInt *pncol;
448: PetscInt ncols;
449: Mat Pnew;
450: PetscInt *lsparse, *gsparse;
451: PetscReal pmax_pos, pmax_neg, ptot_pos, ptot_neg, pthresh_pos, pthresh_neg;
452: PC_MG *mg = (PC_MG *)pc->data;
453: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
454: PC_GAMG_Classical *cls = (PC_GAMG_Classical *)pc_gamg->subctx;
455: MatType mtype;
457: PetscFunctionBegin;
458: /* trim and rescale with reallocation */
459: PetscCall(MatGetOwnershipRange(*P, &ps, &pf));
460: PetscCall(MatGetOwnershipRangeColumn(*P, &pcs, &pcf));
461: pn = pf - ps;
462: pcn = pcf - pcs;
463: PetscCall(PetscMalloc2(pn, &lsparse, pn, &gsparse));
464: /* allocate */
465: cmax = 0;
466: for (i = ps; i < pf; i++) {
467: lsparse[i - ps] = 0;
468: gsparse[i - ps] = 0;
469: PetscCall(MatGetRow(*P, i, &ncols, &pcol, &pval));
470: if (ncols > cmax) cmax = ncols;
471: pmax_pos = 0.;
472: pmax_neg = 0.;
473: for (j = 0; j < ncols; j++) {
474: if (PetscRealPart(pval[j]) > pmax_pos) {
475: pmax_pos = PetscRealPart(pval[j]);
476: } else if (PetscRealPart(pval[j]) < pmax_neg) {
477: pmax_neg = PetscRealPart(pval[j]);
478: }
479: }
480: for (j = 0; j < ncols; j++) {
481: if (PetscRealPart(pval[j]) >= pmax_pos * cls->interp_threshold || PetscRealPart(pval[j]) <= pmax_neg * cls->interp_threshold) {
482: if (pcol[j] >= pcs && pcol[j] < pcf) {
483: lsparse[i - ps]++;
484: } else {
485: gsparse[i - ps]++;
486: }
487: }
488: }
489: PetscCall(MatRestoreRow(*P, i, &ncols, &pcol, &pval));
490: }
492: PetscCall(PetscMalloc2(cmax, &pnval, cmax, &pncol));
494: PetscCall(MatGetType(*P, &mtype));
495: PetscCall(MatCreate(PetscObjectComm((PetscObject)*P), &Pnew));
496: PetscCall(MatSetType(Pnew, mtype));
497: PetscCall(MatSetSizes(Pnew, pn, pcn, PETSC_DETERMINE, PETSC_DETERMINE));
498: PetscCall(MatSeqAIJSetPreallocation(Pnew, 0, lsparse));
499: PetscCall(MatMPIAIJSetPreallocation(Pnew, 0, lsparse, 0, gsparse));
501: for (i = ps; i < pf; i++) {
502: PetscCall(MatGetRow(*P, i, &ncols, &pcol, &pval));
503: pmax_pos = 0.;
504: pmax_neg = 0.;
505: for (j = 0; j < ncols; j++) {
506: if (PetscRealPart(pval[j]) > pmax_pos) {
507: pmax_pos = PetscRealPart(pval[j]);
508: } else if (PetscRealPart(pval[j]) < pmax_neg) {
509: pmax_neg = PetscRealPart(pval[j]);
510: }
511: }
512: pthresh_pos = 0.;
513: pthresh_neg = 0.;
514: ptot_pos = 0.;
515: ptot_neg = 0.;
516: for (j = 0; j < ncols; j++) {
517: if (PetscRealPart(pval[j]) >= cls->interp_threshold * pmax_pos) {
518: pthresh_pos += PetscRealPart(pval[j]);
519: } else if (PetscRealPart(pval[j]) <= cls->interp_threshold * pmax_neg) {
520: pthresh_neg += PetscRealPart(pval[j]);
521: }
522: if (PetscRealPart(pval[j]) > 0.) {
523: ptot_pos += PetscRealPart(pval[j]);
524: } else {
525: ptot_neg += PetscRealPart(pval[j]);
526: }
527: }
528: if (PetscAbsReal(pthresh_pos) > 0.) ptot_pos /= pthresh_pos;
529: if (PetscAbsReal(pthresh_neg) > 0.) ptot_neg /= pthresh_neg;
530: idx = 0;
531: for (j = 0; j < ncols; j++) {
532: if (PetscRealPart(pval[j]) >= pmax_pos * cls->interp_threshold) {
533: pnval[idx] = ptot_pos * pval[j];
534: pncol[idx] = pcol[j];
535: idx++;
536: } else if (PetscRealPart(pval[j]) <= pmax_neg * cls->interp_threshold) {
537: pnval[idx] = ptot_neg * pval[j];
538: pncol[idx] = pcol[j];
539: idx++;
540: }
541: }
542: PetscCall(MatRestoreRow(*P, i, &ncols, &pcol, &pval));
543: PetscCall(MatSetValues(Pnew, 1, &i, idx, pncol, pnval, INSERT_VALUES));
544: }
546: PetscCall(MatAssemblyBegin(Pnew, MAT_FINAL_ASSEMBLY));
547: PetscCall(MatAssemblyEnd(Pnew, MAT_FINAL_ASSEMBLY));
548: PetscCall(MatDestroy(P));
550: *P = Pnew;
551: PetscCall(PetscFree2(lsparse, gsparse));
552: PetscCall(PetscFree2(pnval, pncol));
553: PetscFunctionReturn(PETSC_SUCCESS);
554: }
556: static PetscErrorCode PCGAMGProlongator_Classical_Standard(PC pc, Mat A, PetscCoarsenData *agg_lists, Mat *P)
557: {
558: Mat lA, *lAs;
559: MatType mtype;
560: Vec cv;
561: PetscInt *gcid, *lcid, *lsparse, *gsparse, *picol;
562: PetscInt fs, fe, cs, ce, nl, i, j, k, li, lni, ci, ncols, maxcols, fn, cn, cid;
563: PetscMPIInt size;
564: const PetscInt *lidx, *icol, *gidx;
565: PetscBool iscoarse;
566: PetscScalar vi, pentry, pjentry;
567: PetscScalar *pcontrib, *pvcol;
568: const PetscScalar *vcol;
569: PetscReal diag, jdiag, jwttotal;
570: PetscInt pncols;
571: PetscSF sf;
572: PetscLayout clayout;
573: IS lis;
575: PetscFunctionBegin;
576: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
577: PetscCall(MatGetOwnershipRange(A, &fs, &fe));
578: fn = fe - fs;
579: PetscCall(ISCreateStride(PETSC_COMM_SELF, fe - fs, fs, 1, &lis));
580: if (size > 1) {
581: PetscCall(MatGetLayouts(A, NULL, &clayout));
582: /* increase the overlap by two to get neighbors of neighbors */
583: PetscCall(MatIncreaseOverlap(A, 1, &lis, 2));
584: PetscCall(ISSort(lis));
585: /* get the local part of A */
586: PetscCall(MatCreateSubMatrices(A, 1, &lis, &lis, MAT_INITIAL_MATRIX, &lAs));
587: lA = lAs[0];
588: /* build an SF out of it */
589: PetscCall(ISGetLocalSize(lis, &nl));
590: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)A), &sf));
591: PetscCall(ISGetIndices(lis, &lidx));
592: PetscCall(PetscSFSetGraphLayout(sf, clayout, nl, NULL, PETSC_COPY_VALUES, lidx));
593: PetscCall(ISRestoreIndices(lis, &lidx));
594: } else {
595: lA = A;
596: nl = fn;
597: }
598: /* create a communication structure for the overlapped portion and transmit coarse indices */
599: PetscCall(PetscMalloc3(fn, &lsparse, fn, &gsparse, nl, &pcontrib));
600: /* create coarse vector */
601: cn = 0;
602: for (i = 0; i < fn; i++) {
603: PetscCall(PetscCDIsEmptyAt(agg_lists, i, &iscoarse));
604: if (!iscoarse) cn++;
605: }
606: PetscCall(PetscMalloc1(fn, &gcid));
607: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)A), cn, PETSC_DECIDE, &cv));
608: PetscCall(VecGetOwnershipRange(cv, &cs, &ce));
609: cn = 0;
610: for (i = 0; i < fn; i++) {
611: PetscCall(PetscCDIsEmptyAt(agg_lists, i, &iscoarse));
612: if (!iscoarse) {
613: gcid[i] = cs + cn;
614: cn++;
615: } else {
616: gcid[i] = -1;
617: }
618: }
619: if (size > 1) {
620: PetscCall(PetscMalloc1(nl, &lcid));
621: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, gcid, lcid, MPI_REPLACE));
622: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, gcid, lcid, MPI_REPLACE));
623: } else {
624: lcid = gcid;
625: }
626: /* count to preallocate the prolongator */
627: PetscCall(ISGetIndices(lis, &gidx));
628: maxcols = 0;
629: /* count the number of unique contributing coarse cells for each fine */
630: for (i = 0; i < nl; i++) {
631: pcontrib[i] = 0.;
632: PetscCall(MatGetRow(lA, i, &ncols, &icol, NULL));
633: if (gidx[i] >= fs && gidx[i] < fe) {
634: li = gidx[i] - fs;
635: lsparse[li] = 0;
636: gsparse[li] = 0;
637: cid = lcid[i];
638: if (cid >= 0) {
639: lsparse[li] = 1;
640: } else {
641: for (j = 0; j < ncols; j++) {
642: if (lcid[icol[j]] >= 0) {
643: pcontrib[icol[j]] = 1.;
644: } else {
645: ci = icol[j];
646: PetscCall(MatRestoreRow(lA, i, &ncols, &icol, NULL));
647: PetscCall(MatGetRow(lA, ci, &ncols, &icol, NULL));
648: for (k = 0; k < ncols; k++) {
649: if (lcid[icol[k]] >= 0) pcontrib[icol[k]] = 1.;
650: }
651: PetscCall(MatRestoreRow(lA, ci, &ncols, &icol, NULL));
652: PetscCall(MatGetRow(lA, i, &ncols, &icol, NULL));
653: }
654: }
655: for (j = 0; j < ncols; j++) {
656: if (lcid[icol[j]] >= 0 && pcontrib[icol[j]] != 0.) {
657: lni = lcid[icol[j]];
658: if (lni >= cs && lni < ce) {
659: lsparse[li]++;
660: } else {
661: gsparse[li]++;
662: }
663: pcontrib[icol[j]] = 0.;
664: } else {
665: ci = icol[j];
666: PetscCall(MatRestoreRow(lA, i, &ncols, &icol, NULL));
667: PetscCall(MatGetRow(lA, ci, &ncols, &icol, NULL));
668: for (k = 0; k < ncols; k++) {
669: if (lcid[icol[k]] >= 0 && pcontrib[icol[k]] != 0.) {
670: lni = lcid[icol[k]];
671: if (lni >= cs && lni < ce) {
672: lsparse[li]++;
673: } else {
674: gsparse[li]++;
675: }
676: pcontrib[icol[k]] = 0.;
677: }
678: }
679: PetscCall(MatRestoreRow(lA, ci, &ncols, &icol, NULL));
680: PetscCall(MatGetRow(lA, i, &ncols, &icol, NULL));
681: }
682: }
683: }
684: if (lsparse[li] + gsparse[li] > maxcols) maxcols = lsparse[li] + gsparse[li];
685: }
686: PetscCall(MatRestoreRow(lA, i, &ncols, &icol, &vcol));
687: }
688: PetscCall(PetscMalloc2(maxcols, &picol, maxcols, &pvcol));
689: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), P));
690: PetscCall(MatGetType(A, &mtype));
691: PetscCall(MatSetType(*P, mtype));
692: PetscCall(MatSetSizes(*P, fn, cn, PETSC_DETERMINE, PETSC_DETERMINE));
693: PetscCall(MatMPIAIJSetPreallocation(*P, 0, lsparse, 0, gsparse));
694: PetscCall(MatSeqAIJSetPreallocation(*P, 0, lsparse));
695: for (i = 0; i < nl; i++) {
696: diag = 0.;
697: if (gidx[i] >= fs && gidx[i] < fe) {
698: pncols = 0;
699: cid = lcid[i];
700: if (cid >= 0) {
701: pncols = 1;
702: picol[0] = cid;
703: pvcol[0] = 1.;
704: } else {
705: PetscCall(MatGetRow(lA, i, &ncols, &icol, &vcol));
706: for (j = 0; j < ncols; j++) {
707: pentry = vcol[j];
708: if (lcid[icol[j]] >= 0) {
709: /* coarse neighbor */
710: pcontrib[icol[j]] += pentry;
711: } else if (icol[j] != i) {
712: /* the neighbor is a strongly connected fine node */
713: ci = icol[j];
714: vi = vcol[j];
715: PetscCall(MatRestoreRow(lA, i, &ncols, &icol, &vcol));
716: PetscCall(MatGetRow(lA, ci, &ncols, &icol, &vcol));
717: jwttotal = 0.;
718: jdiag = 0.;
719: for (k = 0; k < ncols; k++) {
720: if (ci == icol[k]) jdiag = PetscRealPart(vcol[k]);
721: }
722: for (k = 0; k < ncols; k++) {
723: if (lcid[icol[k]] >= 0 && jdiag * PetscRealPart(vcol[k]) < 0.) {
724: pjentry = vcol[k];
725: jwttotal += PetscRealPart(pjentry);
726: }
727: }
728: if (jwttotal != 0.) {
729: jwttotal = PetscRealPart(vi) / jwttotal;
730: for (k = 0; k < ncols; k++) {
731: if (lcid[icol[k]] >= 0 && jdiag * PetscRealPart(vcol[k]) < 0.) {
732: pjentry = vcol[k] * jwttotal;
733: pcontrib[icol[k]] += pjentry;
734: }
735: }
736: } else {
737: diag += PetscRealPart(vi);
738: }
739: PetscCall(MatRestoreRow(lA, ci, &ncols, &icol, &vcol));
740: PetscCall(MatGetRow(lA, i, &ncols, &icol, &vcol));
741: } else {
742: diag += PetscRealPart(vcol[j]);
743: }
744: }
745: if (diag != 0.) {
746: diag = 1. / diag;
747: for (j = 0; j < ncols; j++) {
748: if (lcid[icol[j]] >= 0 && pcontrib[icol[j]] != 0.) {
749: /* the neighbor is a coarse node */
750: if (PetscAbsScalar(pcontrib[icol[j]]) > 0.0) {
751: lni = lcid[icol[j]];
752: pvcol[pncols] = -pcontrib[icol[j]] * diag;
753: picol[pncols] = lni;
754: pncols++;
755: }
756: pcontrib[icol[j]] = 0.;
757: } else {
758: /* the neighbor is a strongly connected fine node */
759: ci = icol[j];
760: PetscCall(MatRestoreRow(lA, i, &ncols, &icol, &vcol));
761: PetscCall(MatGetRow(lA, ci, &ncols, &icol, &vcol));
762: for (k = 0; k < ncols; k++) {
763: if (lcid[icol[k]] >= 0 && pcontrib[icol[k]] != 0.) {
764: if (PetscAbsScalar(pcontrib[icol[k]]) > 0.0) {
765: lni = lcid[icol[k]];
766: pvcol[pncols] = -pcontrib[icol[k]] * diag;
767: picol[pncols] = lni;
768: pncols++;
769: }
770: pcontrib[icol[k]] = 0.;
771: }
772: }
773: PetscCall(MatRestoreRow(lA, ci, &ncols, &icol, &vcol));
774: PetscCall(MatGetRow(lA, i, &ncols, &icol, &vcol));
775: }
776: pcontrib[icol[j]] = 0.;
777: }
778: PetscCall(MatRestoreRow(lA, i, &ncols, &icol, &vcol));
779: }
780: }
781: ci = gidx[i];
782: if (pncols > 0) PetscCall(MatSetValues(*P, 1, &ci, pncols, picol, pvcol, INSERT_VALUES));
783: }
784: }
785: PetscCall(ISRestoreIndices(lis, &gidx));
786: PetscCall(PetscFree2(picol, pvcol));
787: PetscCall(PetscFree3(lsparse, gsparse, pcontrib));
788: PetscCall(ISDestroy(&lis));
789: PetscCall(PetscFree(gcid));
790: if (size > 1) {
791: PetscCall(PetscFree(lcid));
792: PetscCall(MatDestroyMatrices(1, &lAs));
793: PetscCall(PetscSFDestroy(&sf));
794: }
795: PetscCall(VecDestroy(&cv));
796: PetscCall(MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY));
797: PetscCall(MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY));
798: PetscFunctionReturn(PETSC_SUCCESS);
799: }
801: static PetscErrorCode PCGAMGOptProlongator_Classical_Jacobi(PC pc, Mat A, Mat *P)
802: {
803: PetscInt f, s, n, cf, cs, i, idx;
804: PetscInt *coarserows;
805: PetscInt ncols;
806: const PetscInt *pcols;
807: const PetscScalar *pvals;
808: Mat Pnew;
809: Vec diag;
810: PC_MG *mg = (PC_MG *)pc->data;
811: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
812: PC_GAMG_Classical *cls = (PC_GAMG_Classical *)pc_gamg->subctx;
814: PetscFunctionBegin;
815: if (cls->nsmooths == 0) {
816: PetscCall(PCGAMGTruncateProlongator_Private(pc, P));
817: PetscFunctionReturn(PETSC_SUCCESS);
818: }
819: PetscCall(MatGetOwnershipRange(*P, &s, &f));
820: n = f - s;
821: PetscCall(MatGetOwnershipRangeColumn(*P, &cs, &cf));
822: PetscCall(PetscMalloc1(n, &coarserows));
823: /* identify the rows corresponding to coarse unknowns */
824: idx = 0;
825: for (i = s; i < f; i++) {
826: PetscCall(MatGetRow(*P, i, &ncols, &pcols, &pvals));
827: /* assume, for now, that it's a coarse unknown if it has a single unit entry */
828: if (ncols == 1) {
829: if (pvals[0] == 1.) {
830: coarserows[idx] = i;
831: idx++;
832: }
833: }
834: PetscCall(MatRestoreRow(*P, i, &ncols, &pcols, &pvals));
835: }
836: PetscCall(MatCreateVecs(A, &diag, NULL));
837: PetscCall(MatGetDiagonal(A, diag));
838: PetscCall(VecReciprocal(diag));
839: for (i = 0; i < cls->nsmooths; i++) {
840: PetscCall(MatMatMult(A, *P, MAT_INITIAL_MATRIX, PETSC_CURRENT, &Pnew));
841: PetscCall(MatZeroRows(Pnew, idx, coarserows, 0., NULL, NULL));
842: PetscCall(MatDiagonalScale(Pnew, diag, NULL));
843: PetscCall(MatAYPX(Pnew, -1.0, *P, DIFFERENT_NONZERO_PATTERN));
844: PetscCall(MatDestroy(P));
845: *P = Pnew;
846: Pnew = NULL;
847: }
848: PetscCall(VecDestroy(&diag));
849: PetscCall(PetscFree(coarserows));
850: PetscCall(PCGAMGTruncateProlongator_Private(pc, P));
851: PetscFunctionReturn(PETSC_SUCCESS);
852: }
854: static PetscErrorCode PCGAMGProlongator_Classical(PC pc, Mat A, PetscCoarsenData *agg_lists, Mat *P)
855: {
856: PetscErrorCode (*f)(PC, Mat, PetscCoarsenData *, Mat *);
857: PC_MG *mg = (PC_MG *)pc->data;
858: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
859: PC_GAMG_Classical *cls = (PC_GAMG_Classical *)pc_gamg->subctx;
861: PetscFunctionBegin;
862: PetscCall(PetscFunctionListFind(PCGAMGClassicalProlongatorList, cls->prolongtype, &f));
863: PetscCheck(f, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Cannot find PCGAMG Classical prolongator type");
864: PetscCall((*f)(pc, A, agg_lists, P));
865: PetscFunctionReturn(PETSC_SUCCESS);
866: }
868: static PetscErrorCode PCGAMGDestroy_Classical(PC pc)
869: {
870: PC_MG *mg = (PC_MG *)pc->data;
871: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
873: PetscFunctionBegin;
874: PetscCall(PetscFree(pc_gamg->subctx));
875: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGClassicalSetType_C", NULL));
876: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGClassicalGetType_C", NULL));
877: PetscFunctionReturn(PETSC_SUCCESS);
878: }
880: static PetscErrorCode PCGAMGSetFromOptions_Classical(PC pc, PetscOptionItems *PetscOptionsObject)
881: {
882: PC_MG *mg = (PC_MG *)pc->data;
883: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
884: PC_GAMG_Classical *cls = (PC_GAMG_Classical *)pc_gamg->subctx;
885: char tname[256];
886: PetscBool flg;
888: PetscFunctionBegin;
889: PetscOptionsHeadBegin(PetscOptionsObject, "GAMG-Classical options");
890: PetscCall(PetscOptionsFList("-pc_gamg_classical_type", "Type of Classical AMG prolongation", "PCGAMGClassicalSetType", PCGAMGClassicalProlongatorList, cls->prolongtype, tname, sizeof(tname), &flg));
891: if (flg) PetscCall(PCGAMGClassicalSetType(pc, tname));
892: PetscCall(PetscOptionsReal("-pc_gamg_classical_interp_threshold", "Threshold for classical interpolator entries", "", cls->interp_threshold, &cls->interp_threshold, NULL));
893: PetscCall(PetscOptionsInt("-pc_gamg_classical_nsmooths", "Threshold for classical interpolator entries", "", cls->nsmooths, &cls->nsmooths, NULL));
894: PetscOptionsHeadEnd();
895: PetscFunctionReturn(PETSC_SUCCESS);
896: }
898: static PetscErrorCode PCGAMGSetData_Classical(PC pc, Mat A)
899: {
900: PC_MG *mg = (PC_MG *)pc->data;
901: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
903: PetscFunctionBegin;
904: /* no data for classical AMG */
905: pc_gamg->data = NULL;
906: pc_gamg->data_cell_cols = 0;
907: pc_gamg->data_cell_rows = 0;
908: pc_gamg->data_sz = 0;
909: PetscFunctionReturn(PETSC_SUCCESS);
910: }
912: static PetscErrorCode PCGAMGClassicalFinalizePackage(void)
913: {
914: PetscFunctionBegin;
915: PCGAMGClassicalPackageInitialized = PETSC_FALSE;
916: PetscCall(PetscFunctionListDestroy(&PCGAMGClassicalProlongatorList));
917: PetscFunctionReturn(PETSC_SUCCESS);
918: }
920: static PetscErrorCode PCGAMGClassicalInitializePackage(void)
921: {
922: PetscFunctionBegin;
923: if (PCGAMGClassicalPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
924: PetscCall(PetscFunctionListAdd(&PCGAMGClassicalProlongatorList, PCGAMGCLASSICALDIRECT, PCGAMGProlongator_Classical_Direct));
925: PetscCall(PetscFunctionListAdd(&PCGAMGClassicalProlongatorList, PCGAMGCLASSICALSTANDARD, PCGAMGProlongator_Classical_Standard));
926: PetscCall(PetscRegisterFinalize(PCGAMGClassicalFinalizePackage));
927: PetscFunctionReturn(PETSC_SUCCESS);
928: }
930: PetscErrorCode PCCreateGAMG_Classical(PC pc)
931: {
932: PC_MG *mg = (PC_MG *)pc->data;
933: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
934: PC_GAMG_Classical *pc_gamg_classical;
936: PetscFunctionBegin;
937: PetscCall(PCGAMGClassicalInitializePackage());
938: if (pc_gamg->subctx) {
939: /* call base class */
940: PetscCall(PCDestroy_GAMG(pc));
941: }
943: /* create sub context for SA */
944: PetscCall(PetscNew(&pc_gamg_classical));
945: pc_gamg->subctx = pc_gamg_classical;
946: pc->ops->setfromoptions = PCGAMGSetFromOptions_Classical;
947: /* reset does not do anything; setup not virtual */
949: /* set internal function pointers */
950: pc_gamg->ops->destroy = PCGAMGDestroy_Classical;
951: pc_gamg->ops->creategraph = PCGAMGCreateGraph_Classical;
952: pc_gamg->ops->coarsen = PCGAMGCoarsen_Classical;
953: pc_gamg->ops->prolongator = PCGAMGProlongator_Classical;
954: pc_gamg->ops->optprolongator = PCGAMGOptProlongator_Classical_Jacobi;
955: pc_gamg->ops->setfromoptions = PCGAMGSetFromOptions_Classical;
957: pc_gamg->ops->createdefaultdata = PCGAMGSetData_Classical;
958: pc_gamg_classical->interp_threshold = 0.2;
959: pc_gamg_classical->nsmooths = 0;
960: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGClassicalSetType_C", PCGAMGClassicalSetType_GAMG));
961: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGClassicalGetType_C", PCGAMGClassicalGetType_GAMG));
962: PetscCall(PCGAMGClassicalSetType(pc, PCGAMGCLASSICALSTANDARD));
963: PetscFunctionReturn(PETSC_SUCCESS);
964: }