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: }