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