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