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