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;
168: PetscFunctionBegin;
169: PetscCheck(G, fcomm, PETSC_ERR_ARG_WRONGSTATE, "Must set Graph in PC in PCGAMG before coarsening");
171: PetscCall(MatCoarsenCreate(fcomm, &crs));
172: PetscCall(MatCoarsenSetFromOptions(crs));
173: PetscCall(MatCoarsenSetAdjacency(crs, *G));
174: PetscCall(MatCoarsenSetStrictAggs(crs, PETSC_TRUE));
175: PetscCall(MatCoarsenApply(crs));
176: PetscCall(MatCoarsenGetData(crs, agg_lists));
177: PetscCall(MatCoarsenDestroy(&crs));
178: PetscFunctionReturn(PETSC_SUCCESS);
179: }
181: static PetscErrorCode PCGAMGProlongator_Classical_Direct(PC pc, Mat A, PetscCoarsenData *agg_lists, Mat *P)
182: {
183: PC_MG *mg = (PC_MG *)pc->data;
184: PC_GAMG *gamg = (PC_GAMG *)mg->innerctx;
185: PetscBool iscoarse, isMPIAIJ, isSEQAIJ;
186: PetscInt fn, cn, fs, fe, cs, ce, i, j, ncols, col, row_f, row_c, cmax = 0, idx, noff;
187: PetscInt *lcid, *gcid, *lsparse, *gsparse, *colmap, *pcols;
188: const PetscInt *rcol;
189: PetscReal *Amax_pos, *Amax_neg;
190: PetscScalar g_pos, g_neg, a_pos, a_neg, diag, invdiag, alpha, beta, pij;
191: PetscScalar *pvals;
192: const PetscScalar *rval;
193: Mat lA, gA = NULL;
194: MatType mtype;
195: Vec C, lvec;
196: PetscLayout clayout;
197: PetscSF sf;
198: Mat_MPIAIJ *mpiaij;
200: PetscFunctionBegin;
201: PetscCall(MatGetOwnershipRange(A, &fs, &fe));
202: fn = fe - fs;
203: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPIAIJ, &isMPIAIJ));
204: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQAIJ, &isSEQAIJ));
205: PetscCheck(isMPIAIJ || isSEQAIJ, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Classical AMG requires MPIAIJ matrix");
206: if (isMPIAIJ) {
207: mpiaij = (Mat_MPIAIJ *)A->data;
208: lA = mpiaij->A;
209: gA = mpiaij->B;
210: lvec = mpiaij->lvec;
211: PetscCall(VecGetSize(lvec, &noff));
212: colmap = mpiaij->garray;
213: PetscCall(MatGetLayouts(A, NULL, &clayout));
214: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)A), &sf));
215: PetscCall(PetscSFSetGraphLayout(sf, clayout, noff, NULL, PETSC_COPY_VALUES, colmap));
216: PetscCall(PetscMalloc1(noff, &gcid));
217: } else {
218: lA = A;
219: }
220: PetscCall(PetscMalloc5(fn, &lsparse, fn, &gsparse, fn, &lcid, fn, &Amax_pos, fn, &Amax_neg));
222: /* count the number of coarse unknowns */
223: cn = 0;
224: for (i = 0; i < fn; i++) {
225: /* filter out singletons */
226: PetscCall(PetscCDIsEmptyAt(agg_lists, i, &iscoarse));
227: lcid[i] = -1;
228: if (!iscoarse) cn++;
229: }
231: /* create the coarse vector */
232: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)A), cn, PETSC_DECIDE, &C));
233: PetscCall(VecGetOwnershipRange(C, &cs, &ce));
235: cn = 0;
236: for (i = 0; i < fn; i++) {
237: PetscCall(PetscCDIsEmptyAt(agg_lists, i, &iscoarse));
238: if (!iscoarse) {
239: lcid[i] = cs + cn;
240: cn++;
241: } else {
242: lcid[i] = -1;
243: }
244: }
246: if (gA) {
247: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, lcid, gcid, MPI_REPLACE));
248: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, lcid, gcid, MPI_REPLACE));
249: }
251: /* determine the largest off-diagonal entries in each row */
252: for (i = fs; i < fe; i++) {
253: Amax_pos[i - fs] = 0.;
254: Amax_neg[i - fs] = 0.;
255: PetscCall(MatGetRow(A, i, &ncols, &rcol, &rval));
256: for (j = 0; j < ncols; j++) {
257: if ((PetscRealPart(-rval[j]) > Amax_neg[i - fs]) && i != rcol[j]) Amax_neg[i - fs] = PetscAbsScalar(rval[j]);
258: if ((PetscRealPart(rval[j]) > Amax_pos[i - fs]) && i != rcol[j]) Amax_pos[i - fs] = PetscAbsScalar(rval[j]);
259: }
260: if (ncols > cmax) cmax = ncols;
261: PetscCall(MatRestoreRow(A, i, &ncols, &rcol, &rval));
262: }
263: PetscCall(PetscMalloc2(cmax, &pcols, cmax, &pvals));
264: PetscCall(VecDestroy(&C));
266: /* count the on and off processor sparsity patterns for the prolongator */
267: for (i = 0; i < fn; i++) {
268: /* on */
269: lsparse[i] = 0;
270: gsparse[i] = 0;
271: if (lcid[i] >= 0) {
272: lsparse[i] = 1;
273: gsparse[i] = 0;
274: } else {
275: PetscCall(MatGetRow(lA, i, &ncols, &rcol, &rval));
276: for (j = 0; j < ncols; j++) {
277: col = rcol[j];
278: 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;
279: }
280: PetscCall(MatRestoreRow(lA, i, &ncols, &rcol, &rval));
281: /* off */
282: if (gA) {
283: PetscCall(MatGetRow(gA, i, &ncols, &rcol, &rval));
284: for (j = 0; j < ncols; j++) {
285: col = rcol[j];
286: 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;
287: }
288: PetscCall(MatRestoreRow(gA, i, &ncols, &rcol, &rval));
289: }
290: }
291: }
293: /* preallocate and create the prolongator */
294: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), P));
295: PetscCall(MatGetType(A, &mtype));
296: PetscCall(MatSetType(*P, mtype));
297: PetscCall(MatSetSizes(*P, fn, cn, PETSC_DETERMINE, PETSC_DETERMINE));
298: PetscCall(MatMPIAIJSetPreallocation(*P, 0, lsparse, 0, gsparse));
299: PetscCall(MatSeqAIJSetPreallocation(*P, 0, lsparse));
301: /* loop over local fine nodes -- get the diagonal, the sum of positive and negative strong and weak weights, and set up the row */
302: for (i = 0; i < fn; i++) {
303: /* determine on or off */
304: row_f = i + fs;
305: row_c = lcid[i];
306: if (row_c >= 0) {
307: pij = 1.;
308: PetscCall(MatSetValues(*P, 1, &row_f, 1, &row_c, &pij, INSERT_VALUES));
309: } else {
310: g_pos = 0.;
311: g_neg = 0.;
312: a_pos = 0.;
313: a_neg = 0.;
314: diag = 0.;
316: /* local connections */
317: PetscCall(MatGetRow(lA, i, &ncols, &rcol, &rval));
318: for (j = 0; j < ncols; j++) {
319: col = rcol[j];
320: if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0] * Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0] * Amax_neg[i])) {
321: if (PetscRealPart(rval[j]) > 0.) {
322: g_pos += rval[j];
323: } else {
324: g_neg += rval[j];
325: }
326: }
327: if (col != i) {
328: if (PetscRealPart(rval[j]) > 0.) {
329: a_pos += rval[j];
330: } else {
331: a_neg += rval[j];
332: }
333: } else {
334: diag = rval[j];
335: }
336: }
337: PetscCall(MatRestoreRow(lA, i, &ncols, &rcol, &rval));
339: /* ghosted connections */
340: if (gA) {
341: PetscCall(MatGetRow(gA, i, &ncols, &rcol, &rval));
342: for (j = 0; j < ncols; j++) {
343: col = rcol[j];
344: if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0] * Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0] * Amax_neg[i])) {
345: if (PetscRealPart(rval[j]) > 0.) {
346: g_pos += rval[j];
347: } else {
348: g_neg += rval[j];
349: }
350: }
351: if (PetscRealPart(rval[j]) > 0.) {
352: a_pos += rval[j];
353: } else {
354: a_neg += rval[j];
355: }
356: }
357: PetscCall(MatRestoreRow(gA, i, &ncols, &rcol, &rval));
358: }
360: if (g_neg == 0.) {
361: alpha = 0.;
362: } else {
363: alpha = -a_neg / g_neg;
364: }
366: if (g_pos == 0.) {
367: diag += a_pos;
368: beta = 0.;
369: } else {
370: beta = -a_pos / g_pos;
371: }
372: if (diag == 0.) {
373: invdiag = 0.;
374: } else invdiag = 1. / diag;
375: /* on */
376: PetscCall(MatGetRow(lA, i, &ncols, &rcol, &rval));
377: idx = 0;
378: for (j = 0; j < ncols; j++) {
379: col = rcol[j];
380: if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0] * Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0] * Amax_neg[i])) {
381: row_f = i + fs;
382: row_c = lcid[col];
383: /* set the values for on-processor ones */
384: if (PetscRealPart(rval[j]) < 0.) {
385: pij = rval[j] * alpha * invdiag;
386: } else {
387: pij = rval[j] * beta * invdiag;
388: }
389: if (PetscAbsScalar(pij) != 0.) {
390: pvals[idx] = pij;
391: pcols[idx] = row_c;
392: idx++;
393: }
394: }
395: }
396: PetscCall(MatRestoreRow(lA, i, &ncols, &rcol, &rval));
397: /* off */
398: if (gA) {
399: PetscCall(MatGetRow(gA, i, &ncols, &rcol, &rval));
400: for (j = 0; j < ncols; j++) {
401: col = rcol[j];
402: if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0] * Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0] * Amax_neg[i])) {
403: row_f = i + fs;
404: row_c = gcid[col];
405: /* set the values for on-processor ones */
406: if (PetscRealPart(rval[j]) < 0.) {
407: pij = rval[j] * alpha * invdiag;
408: } else {
409: pij = rval[j] * beta * invdiag;
410: }
411: if (PetscAbsScalar(pij) != 0.) {
412: pvals[idx] = pij;
413: pcols[idx] = row_c;
414: idx++;
415: }
416: }
417: }
418: PetscCall(MatRestoreRow(gA, i, &ncols, &rcol, &rval));
419: }
420: PetscCall(MatSetValues(*P, 1, &row_f, idx, pcols, pvals, INSERT_VALUES));
421: }
422: }
424: PetscCall(MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY));
425: PetscCall(MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY));
427: PetscCall(PetscFree5(lsparse, gsparse, lcid, Amax_pos, Amax_neg));
429: PetscCall(PetscFree2(pcols, pvals));
430: if (gA) {
431: PetscCall(PetscSFDestroy(&sf));
432: PetscCall(PetscFree(gcid));
433: }
434: PetscFunctionReturn(PETSC_SUCCESS);
435: }
437: static PetscErrorCode PCGAMGTruncateProlongator_Private(PC pc, Mat *P)
438: {
439: PetscInt j, i, ps, pf, pn, pcs, pcf, pcn, idx, cmax;
440: const PetscScalar *pval;
441: const PetscInt *pcol;
442: PetscScalar *pnval;
443: PetscInt *pncol;
444: PetscInt ncols;
445: Mat Pnew;
446: PetscInt *lsparse, *gsparse;
447: PetscReal pmax_pos, pmax_neg, ptot_pos, ptot_neg, pthresh_pos, pthresh_neg;
448: PC_MG *mg = (PC_MG *)pc->data;
449: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
450: PC_GAMG_Classical *cls = (PC_GAMG_Classical *)pc_gamg->subctx;
451: MatType mtype;
453: PetscFunctionBegin;
454: /* trim and rescale with reallocation */
455: PetscCall(MatGetOwnershipRange(*P, &ps, &pf));
456: PetscCall(MatGetOwnershipRangeColumn(*P, &pcs, &pcf));
457: pn = pf - ps;
458: pcn = pcf - pcs;
459: PetscCall(PetscMalloc2(pn, &lsparse, pn, &gsparse));
460: /* allocate */
461: cmax = 0;
462: for (i = ps; i < pf; i++) {
463: lsparse[i - ps] = 0;
464: gsparse[i - ps] = 0;
465: PetscCall(MatGetRow(*P, i, &ncols, &pcol, &pval));
466: if (ncols > cmax) cmax = ncols;
467: pmax_pos = 0.;
468: pmax_neg = 0.;
469: for (j = 0; j < ncols; j++) {
470: if (PetscRealPart(pval[j]) > pmax_pos) {
471: pmax_pos = PetscRealPart(pval[j]);
472: } else if (PetscRealPart(pval[j]) < pmax_neg) {
473: pmax_neg = PetscRealPart(pval[j]);
474: }
475: }
476: for (j = 0; j < ncols; j++) {
477: if (PetscRealPart(pval[j]) >= pmax_pos * cls->interp_threshold || PetscRealPart(pval[j]) <= pmax_neg * cls->interp_threshold) {
478: if (pcol[j] >= pcs && pcol[j] < pcf) {
479: lsparse[i - ps]++;
480: } else {
481: gsparse[i - ps]++;
482: }
483: }
484: }
485: PetscCall(MatRestoreRow(*P, i, &ncols, &pcol, &pval));
486: }
488: PetscCall(PetscMalloc2(cmax, &pnval, cmax, &pncol));
490: PetscCall(MatGetType(*P, &mtype));
491: PetscCall(MatCreate(PetscObjectComm((PetscObject)*P), &Pnew));
492: PetscCall(MatSetType(Pnew, mtype));
493: PetscCall(MatSetSizes(Pnew, pn, pcn, PETSC_DETERMINE, PETSC_DETERMINE));
494: PetscCall(MatSeqAIJSetPreallocation(Pnew, 0, lsparse));
495: PetscCall(MatMPIAIJSetPreallocation(Pnew, 0, lsparse, 0, gsparse));
497: for (i = ps; i < pf; i++) {
498: PetscCall(MatGetRow(*P, i, &ncols, &pcol, &pval));
499: pmax_pos = 0.;
500: pmax_neg = 0.;
501: for (j = 0; j < ncols; j++) {
502: if (PetscRealPart(pval[j]) > pmax_pos) {
503: pmax_pos = PetscRealPart(pval[j]);
504: } else if (PetscRealPart(pval[j]) < pmax_neg) {
505: pmax_neg = PetscRealPart(pval[j]);
506: }
507: }
508: pthresh_pos = 0.;
509: pthresh_neg = 0.;
510: ptot_pos = 0.;
511: ptot_neg = 0.;
512: for (j = 0; j < ncols; j++) {
513: if (PetscRealPart(pval[j]) >= cls->interp_threshold * pmax_pos) {
514: pthresh_pos += PetscRealPart(pval[j]);
515: } else if (PetscRealPart(pval[j]) <= cls->interp_threshold * pmax_neg) {
516: pthresh_neg += PetscRealPart(pval[j]);
517: }
518: if (PetscRealPart(pval[j]) > 0.) {
519: ptot_pos += PetscRealPart(pval[j]);
520: } else {
521: ptot_neg += PetscRealPart(pval[j]);
522: }
523: }
524: if (PetscAbsReal(pthresh_pos) > 0.) ptot_pos /= pthresh_pos;
525: if (PetscAbsReal(pthresh_neg) > 0.) ptot_neg /= pthresh_neg;
526: idx = 0;
527: for (j = 0; j < ncols; j++) {
528: if (PetscRealPart(pval[j]) >= pmax_pos * cls->interp_threshold) {
529: pnval[idx] = ptot_pos * pval[j];
530: pncol[idx] = pcol[j];
531: idx++;
532: } else if (PetscRealPart(pval[j]) <= pmax_neg * cls->interp_threshold) {
533: pnval[idx] = ptot_neg * pval[j];
534: pncol[idx] = pcol[j];
535: idx++;
536: }
537: }
538: PetscCall(MatRestoreRow(*P, i, &ncols, &pcol, &pval));
539: PetscCall(MatSetValues(Pnew, 1, &i, idx, pncol, pnval, INSERT_VALUES));
540: }
542: PetscCall(MatAssemblyBegin(Pnew, MAT_FINAL_ASSEMBLY));
543: PetscCall(MatAssemblyEnd(Pnew, MAT_FINAL_ASSEMBLY));
544: PetscCall(MatDestroy(P));
546: *P = Pnew;
547: PetscCall(PetscFree2(lsparse, gsparse));
548: PetscCall(PetscFree2(pnval, pncol));
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }
552: static PetscErrorCode PCGAMGProlongator_Classical_Standard(PC pc, Mat A, PetscCoarsenData *agg_lists, Mat *P)
553: {
554: Mat lA, *lAs;
555: MatType mtype;
556: Vec cv;
557: PetscInt *gcid, *lcid, *lsparse, *gsparse, *picol;
558: PetscInt fs, fe, cs, ce, nl, i, j, k, li, lni, ci, ncols, maxcols, fn, cn, cid;
559: PetscMPIInt size;
560: const PetscInt *lidx, *icol, *gidx;
561: PetscBool iscoarse;
562: PetscScalar vi, pentry, pjentry;
563: PetscScalar *pcontrib, *pvcol;
564: const PetscScalar *vcol;
565: PetscReal diag, jdiag, jwttotal;
566: PetscInt pncols;
567: PetscSF sf;
568: PetscLayout clayout;
569: IS lis;
571: PetscFunctionBegin;
572: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
573: PetscCall(MatGetOwnershipRange(A, &fs, &fe));
574: fn = fe - fs;
575: PetscCall(ISCreateStride(PETSC_COMM_SELF, fe - fs, fs, 1, &lis));
576: if (size > 1) {
577: PetscCall(MatGetLayouts(A, NULL, &clayout));
578: /* increase the overlap by two to get neighbors of neighbors */
579: PetscCall(MatIncreaseOverlap(A, 1, &lis, 2));
580: PetscCall(ISSort(lis));
581: /* get the local part of A */
582: PetscCall(MatCreateSubMatrices(A, 1, &lis, &lis, MAT_INITIAL_MATRIX, &lAs));
583: lA = lAs[0];
584: /* build an SF out of it */
585: PetscCall(ISGetLocalSize(lis, &nl));
586: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)A), &sf));
587: PetscCall(ISGetIndices(lis, &lidx));
588: PetscCall(PetscSFSetGraphLayout(sf, clayout, nl, NULL, PETSC_COPY_VALUES, lidx));
589: PetscCall(ISRestoreIndices(lis, &lidx));
590: } else {
591: lA = A;
592: nl = fn;
593: }
594: /* create a communication structure for the overlapped portion and transmit coarse indices */
595: PetscCall(PetscMalloc3(fn, &lsparse, fn, &gsparse, nl, &pcontrib));
596: /* create coarse vector */
597: cn = 0;
598: for (i = 0; i < fn; i++) {
599: PetscCall(PetscCDIsEmptyAt(agg_lists, i, &iscoarse));
600: if (!iscoarse) cn++;
601: }
602: PetscCall(PetscMalloc1(fn, &gcid));
603: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)A), cn, PETSC_DECIDE, &cv));
604: PetscCall(VecGetOwnershipRange(cv, &cs, &ce));
605: cn = 0;
606: for (i = 0; i < fn; i++) {
607: PetscCall(PetscCDIsEmptyAt(agg_lists, i, &iscoarse));
608: if (!iscoarse) {
609: gcid[i] = cs + cn;
610: cn++;
611: } else {
612: gcid[i] = -1;
613: }
614: }
615: if (size > 1) {
616: PetscCall(PetscMalloc1(nl, &lcid));
617: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, gcid, lcid, MPI_REPLACE));
618: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, gcid, lcid, MPI_REPLACE));
619: } else {
620: lcid = gcid;
621: }
622: /* count to preallocate the prolongator */
623: PetscCall(ISGetIndices(lis, &gidx));
624: maxcols = 0;
625: /* count the number of unique contributing coarse cells for each fine */
626: for (i = 0; i < nl; i++) {
627: pcontrib[i] = 0.;
628: PetscCall(MatGetRow(lA, i, &ncols, &icol, NULL));
629: if (gidx[i] >= fs && gidx[i] < fe) {
630: li = gidx[i] - fs;
631: lsparse[li] = 0;
632: gsparse[li] = 0;
633: cid = lcid[i];
634: if (cid >= 0) {
635: lsparse[li] = 1;
636: } else {
637: for (j = 0; j < ncols; j++) {
638: if (lcid[icol[j]] >= 0) {
639: pcontrib[icol[j]] = 1.;
640: } else {
641: ci = icol[j];
642: PetscCall(MatRestoreRow(lA, i, &ncols, &icol, NULL));
643: PetscCall(MatGetRow(lA, ci, &ncols, &icol, NULL));
644: for (k = 0; k < ncols; k++) {
645: if (lcid[icol[k]] >= 0) pcontrib[icol[k]] = 1.;
646: }
647: PetscCall(MatRestoreRow(lA, ci, &ncols, &icol, NULL));
648: PetscCall(MatGetRow(lA, i, &ncols, &icol, NULL));
649: }
650: }
651: for (j = 0; j < ncols; j++) {
652: if (lcid[icol[j]] >= 0 && pcontrib[icol[j]] != 0.) {
653: lni = lcid[icol[j]];
654: if (lni >= cs && lni < ce) {
655: lsparse[li]++;
656: } else {
657: gsparse[li]++;
658: }
659: pcontrib[icol[j]] = 0.;
660: } else {
661: ci = icol[j];
662: PetscCall(MatRestoreRow(lA, i, &ncols, &icol, NULL));
663: PetscCall(MatGetRow(lA, ci, &ncols, &icol, NULL));
664: for (k = 0; k < ncols; k++) {
665: if (lcid[icol[k]] >= 0 && pcontrib[icol[k]] != 0.) {
666: lni = lcid[icol[k]];
667: if (lni >= cs && lni < ce) {
668: lsparse[li]++;
669: } else {
670: gsparse[li]++;
671: }
672: pcontrib[icol[k]] = 0.;
673: }
674: }
675: PetscCall(MatRestoreRow(lA, ci, &ncols, &icol, NULL));
676: PetscCall(MatGetRow(lA, i, &ncols, &icol, NULL));
677: }
678: }
679: }
680: if (lsparse[li] + gsparse[li] > maxcols) maxcols = lsparse[li] + gsparse[li];
681: }
682: PetscCall(MatRestoreRow(lA, i, &ncols, &icol, &vcol));
683: }
684: PetscCall(PetscMalloc2(maxcols, &picol, maxcols, &pvcol));
685: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), P));
686: PetscCall(MatGetType(A, &mtype));
687: PetscCall(MatSetType(*P, mtype));
688: PetscCall(MatSetSizes(*P, fn, cn, PETSC_DETERMINE, PETSC_DETERMINE));
689: PetscCall(MatMPIAIJSetPreallocation(*P, 0, lsparse, 0, gsparse));
690: PetscCall(MatSeqAIJSetPreallocation(*P, 0, lsparse));
691: for (i = 0; i < nl; i++) {
692: diag = 0.;
693: if (gidx[i] >= fs && gidx[i] < fe) {
694: pncols = 0;
695: cid = lcid[i];
696: if (cid >= 0) {
697: pncols = 1;
698: picol[0] = cid;
699: pvcol[0] = 1.;
700: } else {
701: PetscCall(MatGetRow(lA, i, &ncols, &icol, &vcol));
702: for (j = 0; j < ncols; j++) {
703: pentry = vcol[j];
704: if (lcid[icol[j]] >= 0) {
705: /* coarse neighbor */
706: pcontrib[icol[j]] += pentry;
707: } else if (icol[j] != i) {
708: /* the neighbor is a strongly connected fine node */
709: ci = icol[j];
710: vi = vcol[j];
711: PetscCall(MatRestoreRow(lA, i, &ncols, &icol, &vcol));
712: PetscCall(MatGetRow(lA, ci, &ncols, &icol, &vcol));
713: jwttotal = 0.;
714: jdiag = 0.;
715: for (k = 0; k < ncols; k++) {
716: if (ci == icol[k]) jdiag = PetscRealPart(vcol[k]);
717: }
718: for (k = 0; k < ncols; k++) {
719: if (lcid[icol[k]] >= 0 && jdiag * PetscRealPart(vcol[k]) < 0.) {
720: pjentry = vcol[k];
721: jwttotal += PetscRealPart(pjentry);
722: }
723: }
724: if (jwttotal != 0.) {
725: jwttotal = PetscRealPart(vi) / jwttotal;
726: for (k = 0; k < ncols; k++) {
727: if (lcid[icol[k]] >= 0 && jdiag * PetscRealPart(vcol[k]) < 0.) {
728: pjentry = vcol[k] * jwttotal;
729: pcontrib[icol[k]] += pjentry;
730: }
731: }
732: } else {
733: diag += PetscRealPart(vi);
734: }
735: PetscCall(MatRestoreRow(lA, ci, &ncols, &icol, &vcol));
736: PetscCall(MatGetRow(lA, i, &ncols, &icol, &vcol));
737: } else {
738: diag += PetscRealPart(vcol[j]);
739: }
740: }
741: if (diag != 0.) {
742: diag = 1. / diag;
743: for (j = 0; j < ncols; j++) {
744: if (lcid[icol[j]] >= 0 && pcontrib[icol[j]] != 0.) {
745: /* the neighbor is a coarse node */
746: if (PetscAbsScalar(pcontrib[icol[j]]) > 0.0) {
747: lni = lcid[icol[j]];
748: pvcol[pncols] = -pcontrib[icol[j]] * diag;
749: picol[pncols] = lni;
750: pncols++;
751: }
752: pcontrib[icol[j]] = 0.;
753: } else {
754: /* the neighbor is a strongly connected fine node */
755: ci = icol[j];
756: PetscCall(MatRestoreRow(lA, i, &ncols, &icol, &vcol));
757: PetscCall(MatGetRow(lA, ci, &ncols, &icol, &vcol));
758: for (k = 0; k < ncols; k++) {
759: if (lcid[icol[k]] >= 0 && pcontrib[icol[k]] != 0.) {
760: if (PetscAbsScalar(pcontrib[icol[k]]) > 0.0) {
761: lni = lcid[icol[k]];
762: pvcol[pncols] = -pcontrib[icol[k]] * diag;
763: picol[pncols] = lni;
764: pncols++;
765: }
766: pcontrib[icol[k]] = 0.;
767: }
768: }
769: PetscCall(MatRestoreRow(lA, ci, &ncols, &icol, &vcol));
770: PetscCall(MatGetRow(lA, i, &ncols, &icol, &vcol));
771: }
772: pcontrib[icol[j]] = 0.;
773: }
774: PetscCall(MatRestoreRow(lA, i, &ncols, &icol, &vcol));
775: }
776: }
777: ci = gidx[i];
778: if (pncols > 0) PetscCall(MatSetValues(*P, 1, &ci, pncols, picol, pvcol, INSERT_VALUES));
779: }
780: }
781: PetscCall(ISRestoreIndices(lis, &gidx));
782: PetscCall(PetscFree2(picol, pvcol));
783: PetscCall(PetscFree3(lsparse, gsparse, pcontrib));
784: PetscCall(ISDestroy(&lis));
785: PetscCall(PetscFree(gcid));
786: if (size > 1) {
787: PetscCall(PetscFree(lcid));
788: PetscCall(MatDestroyMatrices(1, &lAs));
789: PetscCall(PetscSFDestroy(&sf));
790: }
791: PetscCall(VecDestroy(&cv));
792: PetscCall(MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY));
793: PetscCall(MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY));
794: PetscFunctionReturn(PETSC_SUCCESS);
795: }
797: static PetscErrorCode PCGAMGOptProlongator_Classical_Jacobi(PC pc, Mat A, Mat *P)
798: {
799: PetscInt f, s, n, cf, cs, i, idx;
800: PetscInt *coarserows;
801: PetscInt ncols;
802: const PetscInt *pcols;
803: const PetscScalar *pvals;
804: Mat Pnew;
805: Vec diag;
806: PC_MG *mg = (PC_MG *)pc->data;
807: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
808: PC_GAMG_Classical *cls = (PC_GAMG_Classical *)pc_gamg->subctx;
810: PetscFunctionBegin;
811: if (cls->nsmooths == 0) {
812: PetscCall(PCGAMGTruncateProlongator_Private(pc, P));
813: PetscFunctionReturn(PETSC_SUCCESS);
814: }
815: PetscCall(MatGetOwnershipRange(*P, &s, &f));
816: n = f - s;
817: PetscCall(MatGetOwnershipRangeColumn(*P, &cs, &cf));
818: PetscCall(PetscMalloc1(n, &coarserows));
819: /* identify the rows corresponding to coarse unknowns */
820: idx = 0;
821: for (i = s; i < f; i++) {
822: PetscCall(MatGetRow(*P, i, &ncols, &pcols, &pvals));
823: /* assume, for now, that it's a coarse unknown if it has a single unit entry */
824: if (ncols == 1) {
825: if (pvals[0] == 1.) {
826: coarserows[idx] = i;
827: idx++;
828: }
829: }
830: PetscCall(MatRestoreRow(*P, i, &ncols, &pcols, &pvals));
831: }
832: PetscCall(MatCreateVecs(A, &diag, NULL));
833: PetscCall(MatGetDiagonal(A, diag));
834: PetscCall(VecReciprocal(diag));
835: for (i = 0; i < cls->nsmooths; i++) {
836: PetscCall(MatMatMult(A, *P, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Pnew));
837: PetscCall(MatZeroRows(Pnew, idx, coarserows, 0., NULL, NULL));
838: PetscCall(MatDiagonalScale(Pnew, diag, NULL));
839: PetscCall(MatAYPX(Pnew, -1.0, *P, DIFFERENT_NONZERO_PATTERN));
840: PetscCall(MatDestroy(P));
841: *P = Pnew;
842: Pnew = NULL;
843: }
844: PetscCall(VecDestroy(&diag));
845: PetscCall(PetscFree(coarserows));
846: PetscCall(PCGAMGTruncateProlongator_Private(pc, P));
847: PetscFunctionReturn(PETSC_SUCCESS);
848: }
850: static PetscErrorCode PCGAMGProlongator_Classical(PC pc, Mat A, PetscCoarsenData *agg_lists, Mat *P)
851: {
852: PetscErrorCode (*f)(PC, Mat, PetscCoarsenData *, Mat *);
853: PC_MG *mg = (PC_MG *)pc->data;
854: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
855: PC_GAMG_Classical *cls = (PC_GAMG_Classical *)pc_gamg->subctx;
857: PetscFunctionBegin;
858: PetscCall(PetscFunctionListFind(PCGAMGClassicalProlongatorList, cls->prolongtype, &f));
859: PetscCheck(f, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Cannot find PCGAMG Classical prolongator type");
860: PetscCall((*f)(pc, A, agg_lists, P));
861: PetscFunctionReturn(PETSC_SUCCESS);
862: }
864: static PetscErrorCode PCGAMGDestroy_Classical(PC pc)
865: {
866: PC_MG *mg = (PC_MG *)pc->data;
867: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
869: PetscFunctionBegin;
870: PetscCall(PetscFree(pc_gamg->subctx));
871: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGClassicalSetType_C", NULL));
872: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGClassicalGetType_C", NULL));
873: PetscFunctionReturn(PETSC_SUCCESS);
874: }
876: static PetscErrorCode PCGAMGSetFromOptions_Classical(PC pc, PetscOptionItems *PetscOptionsObject)
877: {
878: PC_MG *mg = (PC_MG *)pc->data;
879: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
880: PC_GAMG_Classical *cls = (PC_GAMG_Classical *)pc_gamg->subctx;
881: char tname[256];
882: PetscBool flg;
884: PetscFunctionBegin;
885: PetscOptionsHeadBegin(PetscOptionsObject, "GAMG-Classical options");
886: PetscCall(PetscOptionsFList("-pc_gamg_classical_type", "Type of Classical AMG prolongation", "PCGAMGClassicalSetType", PCGAMGClassicalProlongatorList, cls->prolongtype, tname, sizeof(tname), &flg));
887: if (flg) PetscCall(PCGAMGClassicalSetType(pc, tname));
888: PetscCall(PetscOptionsReal("-pc_gamg_classical_interp_threshold", "Threshold for classical interpolator entries", "", cls->interp_threshold, &cls->interp_threshold, NULL));
889: PetscCall(PetscOptionsInt("-pc_gamg_classical_nsmooths", "Threshold for classical interpolator entries", "", cls->nsmooths, &cls->nsmooths, NULL));
890: PetscOptionsHeadEnd();
891: PetscFunctionReturn(PETSC_SUCCESS);
892: }
894: static PetscErrorCode PCGAMGSetData_Classical(PC pc, Mat A)
895: {
896: PC_MG *mg = (PC_MG *)pc->data;
897: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
899: PetscFunctionBegin;
900: /* no data for classical AMG */
901: pc_gamg->data = NULL;
902: pc_gamg->data_cell_cols = 0;
903: pc_gamg->data_cell_rows = 0;
904: pc_gamg->data_sz = 0;
905: PetscFunctionReturn(PETSC_SUCCESS);
906: }
908: static PetscErrorCode PCGAMGClassicalFinalizePackage(void)
909: {
910: PetscFunctionBegin;
911: PCGAMGClassicalPackageInitialized = PETSC_FALSE;
912: PetscCall(PetscFunctionListDestroy(&PCGAMGClassicalProlongatorList));
913: PetscFunctionReturn(PETSC_SUCCESS);
914: }
916: static PetscErrorCode PCGAMGClassicalInitializePackage(void)
917: {
918: PetscFunctionBegin;
919: if (PCGAMGClassicalPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
920: PetscCall(PetscFunctionListAdd(&PCGAMGClassicalProlongatorList, PCGAMGCLASSICALDIRECT, PCGAMGProlongator_Classical_Direct));
921: PetscCall(PetscFunctionListAdd(&PCGAMGClassicalProlongatorList, PCGAMGCLASSICALSTANDARD, PCGAMGProlongator_Classical_Standard));
922: PetscCall(PetscRegisterFinalize(PCGAMGClassicalFinalizePackage));
923: PetscFunctionReturn(PETSC_SUCCESS);
924: }
926: PetscErrorCode PCCreateGAMG_Classical(PC pc)
927: {
928: PC_MG *mg = (PC_MG *)pc->data;
929: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
930: PC_GAMG_Classical *pc_gamg_classical;
932: PetscFunctionBegin;
933: PetscCall(PCGAMGClassicalInitializePackage());
934: if (pc_gamg->subctx) {
935: /* call base class */
936: PetscCall(PCDestroy_GAMG(pc));
937: }
939: /* create sub context for SA */
940: PetscCall(PetscNew(&pc_gamg_classical));
941: pc_gamg->subctx = pc_gamg_classical;
942: pc->ops->setfromoptions = PCGAMGSetFromOptions_Classical;
943: /* reset does not do anything; setup not virtual */
945: /* set internal function pointers */
946: pc_gamg->ops->destroy = PCGAMGDestroy_Classical;
947: pc_gamg->ops->creategraph = PCGAMGCreateGraph_Classical;
948: pc_gamg->ops->coarsen = PCGAMGCoarsen_Classical;
949: pc_gamg->ops->prolongator = PCGAMGProlongator_Classical;
950: pc_gamg->ops->optprolongator = PCGAMGOptProlongator_Classical_Jacobi;
951: pc_gamg->ops->setfromoptions = PCGAMGSetFromOptions_Classical;
953: pc_gamg->ops->createdefaultdata = PCGAMGSetData_Classical;
954: pc_gamg_classical->interp_threshold = 0.2;
955: pc_gamg_classical->nsmooths = 0;
956: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGClassicalSetType_C", PCGAMGClassicalSetType_GAMG));
957: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGClassicalGetType_C", PCGAMGClassicalGetType_GAMG));
958: PetscCall(PCGAMGClassicalSetType(pc, PCGAMGCLASSICALSTANDARD));
959: PetscFunctionReturn(PETSC_SUCCESS);
960: }