Actual source code: agg.c

  1: /*
  2:  GAMG geometric-algebric multigrid PC - Mark Adams 2011
  3:  */

  5: #include <../src/ksp/pc/impls/gamg/gamg.h>
  6: #include <petscblaslapack.h>
  7: #include <petscdm.h>
  8: #include <petsc/private/kspimpl.h>

 10: typedef struct {
 11:   PetscInt   nsmooths;
 12:   PetscInt   aggressive_coarsening_levels; // number of aggressive coarsening levels (square or MISk)
 13:   MatCoarsen crs;
 14: } PC_GAMG_AGG;

 16: /*@
 17:    PCGAMGSetNSmooths - Set number of smoothing steps (1 is typical) used for multigrid on all the levels

 19:    Logically Collective

 21:    Input Parameters:
 22: .  pc - the preconditioner context

 24:    Options Database Key:
 25: .  -pc_gamg_agg_nsmooths <nsmooth, default=1> - number of smoothing steps to use with smooth aggregation

 27:    Level: intermediate

 29: .seealso: `PCMG`, `PCGAMG`
 30: @*/
 31: PetscErrorCode PCGAMGSetNSmooths(PC pc, PetscInt n)
 32: {
 35:   PetscTryMethod(pc, "PCGAMGSetNSmooths_C", (PC, PetscInt), (pc, n));
 36:   return 0;
 37: }

 39: static PetscErrorCode PCGAMGSetNSmooths_AGG(PC pc, PetscInt n)
 40: {
 41:   PC_MG       *mg          = (PC_MG *)pc->data;
 42:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
 43:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;

 45:   pc_gamg_agg->nsmooths = n;
 46:   return 0;
 47: }

 49: /*@
 50:    PCGAMGSetAggressiveLevels -  Use aggressive coarsening on first n levels

 52:    Logically Collective

 54:    Input Parameters:
 55: +  pc - the preconditioner context
 56: -  n - 0, 1 or more

 58:    Options Database Key:
 59: .  -pc_gamg_aggressive_coarsening <n,default = 1> - Number of levels to square the graph on before aggregating it

 61:    Level: intermediate

 63: .seealso: `PCGAMG`, `PCGAMGSetThreshold()`
 64: @*/
 65: PetscErrorCode PCGAMGSetAggressiveLevels(PC pc, PetscInt n)
 66: {
 69:   PetscTryMethod(pc, "PCGAMGSetAggressiveLevels_C", (PC, PetscInt), (pc, n));
 70:   return 0;
 71: }

 73: static PetscErrorCode PCGAMGSetAggressiveLevels_AGG(PC pc, PetscInt n)
 74: {
 75:   PC_MG       *mg          = (PC_MG *)pc->data;
 76:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
 77:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;

 79:   pc_gamg_agg->aggressive_coarsening_levels = n;
 80:   return 0;
 81: }

 83: static PetscErrorCode PCSetFromOptions_GAMG_AGG(PC pc, PetscOptionItems *PetscOptionsObject)
 84: {
 85:   PC_MG       *mg          = (PC_MG *)pc->data;
 86:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
 87:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;

 89:   PetscOptionsHeadBegin(PetscOptionsObject, "GAMG-AGG options");
 90:   {
 91:     PetscBool flg;
 92:     PetscOptionsInt("-pc_gamg_agg_nsmooths", "smoothing steps for smoothed aggregation, usually 1", "PCGAMGSetNSmooths", pc_gamg_agg->nsmooths, &pc_gamg_agg->nsmooths, NULL);
 93:     pc_gamg_agg->aggressive_coarsening_levels = 1;
 94:     PetscCall(
 95:       PetscOptionsInt("-pc_gamg_square_graph", "Number of aggressive coarsening (MIS-2) levels from finest (alias for -pc_gamg_aggressive_coarsening, deprecated)", "PCGAMGSetAggressiveLevels", pc_gamg_agg->aggressive_coarsening_levels, &pc_gamg_agg->aggressive_coarsening_levels, &flg));
 96:     if (!flg) {
 97:       PetscOptionsInt("-pc_gamg_aggressive_coarsening", "Number of aggressive coarsening (MIS-2) levels from finest", "PCGAMGSetAggressiveLevels", pc_gamg_agg->aggressive_coarsening_levels, &pc_gamg_agg->aggressive_coarsening_levels, NULL);
 98:     } else {
 99:       PetscOptionsInt("-pc_gamg_aggressive_coarsening", "Number of aggressive coarsening (MIS-2) levels from finest", "PCGAMGSetAggressiveLevels", pc_gamg_agg->aggressive_coarsening_levels, &pc_gamg_agg->aggressive_coarsening_levels, &flg);
100:       if (flg) PetscInfo(pc, "Warning: both -pc_gamg_square_graph and -pc_gamg_aggressive_coarsening are used. -pc_gamg_square_graph is deprecated, Number of aggressive levels is %d\n", (int)pc_gamg_agg->aggressive_coarsening_levels);
101:     }
102:   }
103:   PetscOptionsHeadEnd();
104:   return 0;
105: }

107: static PetscErrorCode PCDestroy_GAMG_AGG(PC pc)
108: {
109:   PC_MG   *mg      = (PC_MG *)pc->data;
110:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;

112:   PetscFree(pc_gamg->subctx);
113:   PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", NULL);
114:   PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", NULL);
115:   PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL);
116:   return 0;
117: }

119: /*
120:    PCSetCoordinates_AGG
121:      - collective

123:    Input Parameter:
124:    . pc - the preconditioner context
125:    . ndm - dimesion of data (used for dof/vertex for Stokes)
126:    . a_nloc - number of vertices local
127:    . coords - [a_nloc][ndm] - interleaved coordinate data: {x_0, y_0, z_0, x_1, y_1, ...}
128: */

130: static PetscErrorCode PCSetCoordinates_AGG(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords)
131: {
132:   PC_MG   *mg      = (PC_MG *)pc->data;
133:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
134:   PetscInt arrsz, kk, ii, jj, nloc, ndatarows, ndf;
135:   Mat      mat = pc->pmat;

139:   nloc = a_nloc;

141:   /* SA: null space vectors */
142:   MatGetBlockSize(mat, &ndf);               /* this does not work for Stokes */
143:   if (coords && ndf == 1) pc_gamg->data_cell_cols = 1; /* scalar w/ coords and SA (not needed) */
144:   else if (coords) {
146:     pc_gamg->data_cell_cols = (ndm == 2 ? 3 : 6); /* displacement elasticity */
148:   } else pc_gamg->data_cell_cols = ndf; /* no data, force SA with constant null space vectors */
149:   pc_gamg->data_cell_rows = ndatarows = ndf;
151:   arrsz = nloc * pc_gamg->data_cell_rows * pc_gamg->data_cell_cols;

153:   if (!pc_gamg->data || (pc_gamg->data_sz != arrsz)) {
154:     PetscFree(pc_gamg->data);
155:     PetscMalloc1(arrsz + 1, &pc_gamg->data);
156:   }
157:   /* copy data in - column oriented */
158:   for (kk = 0; kk < nloc; kk++) {
159:     const PetscInt M    = nloc * pc_gamg->data_cell_rows; /* stride into data */
160:     PetscReal     *data = &pc_gamg->data[kk * ndatarows]; /* start of cell */
161:     if (pc_gamg->data_cell_cols == 1) *data = 1.0;
162:     else {
163:       /* translational modes */
164:       for (ii = 0; ii < ndatarows; ii++) {
165:         for (jj = 0; jj < ndatarows; jj++) {
166:           if (ii == jj) data[ii * M + jj] = 1.0;
167:           else data[ii * M + jj] = 0.0;
168:         }
169:       }

171:       /* rotational modes */
172:       if (coords) {
173:         if (ndm == 2) {
174:           data += 2 * M;
175:           data[0] = -coords[2 * kk + 1];
176:           data[1] = coords[2 * kk];
177:         } else {
178:           data += 3 * M;
179:           data[0]         = 0.0;
180:           data[M + 0]     = coords[3 * kk + 2];
181:           data[2 * M + 0] = -coords[3 * kk + 1];
182:           data[1]         = -coords[3 * kk + 2];
183:           data[M + 1]     = 0.0;
184:           data[2 * M + 1] = coords[3 * kk];
185:           data[2]         = coords[3 * kk + 1];
186:           data[M + 2]     = -coords[3 * kk];
187:           data[2 * M + 2] = 0.0;
188:         }
189:       }
190:     }
191:   }
192:   pc_gamg->data_sz = arrsz;
193:   return 0;
194: }

196: /*
197:    PCSetData_AGG - called if data is not set with PCSetCoordinates.
198:       Looks in Mat for near null space.
199:       Does not work for Stokes

201:   Input Parameter:
202:    . pc -
203:    . a_A - matrix to get (near) null space out of.
204: */
205: static PetscErrorCode PCSetData_AGG(PC pc, Mat a_A)
206: {
207:   PC_MG       *mg      = (PC_MG *)pc->data;
208:   PC_GAMG     *pc_gamg = (PC_GAMG *)mg->innerctx;
209:   MatNullSpace mnull;

211:   MatGetNearNullSpace(a_A, &mnull);
212:   if (!mnull) {
213:     DM dm;
214:     PCGetDM(pc, &dm);
215:     if (!dm) MatGetDM(a_A, &dm);
216:     if (dm) {
217:       PetscObject deformation;
218:       PetscInt    Nf;

220:       DMGetNumFields(dm, &Nf);
221:       if (Nf) {
222:         DMGetField(dm, 0, NULL, &deformation);
223:         PetscObjectQuery((PetscObject)deformation, "nearnullspace", (PetscObject *)&mnull);
224:         if (!mnull) PetscObjectQuery((PetscObject)deformation, "nullspace", (PetscObject *)&mnull);
225:       }
226:     }
227:   }

229:   if (!mnull) {
230:     PetscInt bs, NN, MM;
231:     MatGetBlockSize(a_A, &bs);
232:     MatGetLocalSize(a_A, &MM, &NN);
234:     PCSetCoordinates_AGG(pc, bs, MM / bs, NULL);
235:   } else {
236:     PetscReal         *nullvec;
237:     PetscBool          has_const;
238:     PetscInt           i, j, mlocal, nvec, bs;
239:     const Vec         *vecs;
240:     const PetscScalar *v;

242:     MatGetLocalSize(a_A, &mlocal, NULL);
243:     MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs);
244:     pc_gamg->data_sz = (nvec + !!has_const) * mlocal;
245:     PetscMalloc1((nvec + !!has_const) * mlocal, &nullvec);
246:     if (has_const)
247:       for (i = 0; i < mlocal; i++) nullvec[i] = 1.0;
248:     for (i = 0; i < nvec; i++) {
249:       VecGetArrayRead(vecs[i], &v);
250:       for (j = 0; j < mlocal; j++) nullvec[(i + !!has_const) * mlocal + j] = PetscRealPart(v[j]);
251:       VecRestoreArrayRead(vecs[i], &v);
252:     }
253:     pc_gamg->data           = nullvec;
254:     pc_gamg->data_cell_cols = (nvec + !!has_const);
255:     MatGetBlockSize(a_A, &bs);
256:     pc_gamg->data_cell_rows = bs;
257:   }
258:   return 0;
259: }

261: /*
262:   formProl0 - collect null space data for each aggregate, do QR, put R in coarse grid data and Q in P_0

264:   Input Parameter:
265:    . agg_llists - list of arrays with aggregates -- list from selected vertices of aggregate unselected vertices
266:    . bs - row block size
267:    . nSAvec - column bs of new P
268:    . my0crs - global index of start of locals
269:    . data_stride - bs*(nloc nodes + ghost nodes) [data_stride][nSAvec]
270:    . data_in[data_stride*nSAvec] - local data on fine grid
271:    . flid_fgid[data_stride/bs] - make local to global IDs, includes ghosts in 'locals_llist'

273:   Output Parameter:
274:    . a_data_out - in with fine grid data (w/ghosts), out with coarse grid data
275:    . a_Prol - prolongation operator
276: */
277: static PetscErrorCode formProl0(PetscCoarsenData *agg_llists, PetscInt bs, PetscInt nSAvec, PetscInt my0crs, PetscInt data_stride, PetscReal data_in[], const PetscInt flid_fgid[], PetscReal **a_data_out, Mat a_Prol)
278: {
279:   PetscInt        Istart, my0, Iend, nloc, clid, flid = 0, aggID, kk, jj, ii, mm, nSelected, minsz, nghosts, out_data_stride;
280:   MPI_Comm        comm;
281:   PetscReal      *out_data;
282:   PetscCDIntNd   *pos;
283:   PCGAMGHashTable fgid_flid;

285:   PetscObjectGetComm((PetscObject)a_Prol, &comm);
286:   MatGetOwnershipRange(a_Prol, &Istart, &Iend);
287:   nloc = (Iend - Istart) / bs;
288:   my0  = Istart / bs;
290:   Iend /= bs;
291:   nghosts = data_stride / bs - nloc;

293:   PCGAMGHashTableCreate(2 * nghosts + 1, &fgid_flid);
294:   for (kk = 0; kk < nghosts; kk++) PCGAMGHashTableAdd(&fgid_flid, flid_fgid[nloc + kk], nloc + kk);

296:   /* count selected -- same as number of cols of P */
297:   for (nSelected = mm = 0; mm < nloc; mm++) {
298:     PetscBool ise;
299:     PetscCDEmptyAt(agg_llists, mm, &ise);
300:     if (!ise) nSelected++;
301:   }
302:   MatGetOwnershipRangeColumn(a_Prol, &ii, &jj);

306:   /* aloc space for coarse point data (output) */
307:   out_data_stride = nSelected * nSAvec;

309:   PetscMalloc1(out_data_stride * nSAvec, &out_data);
310:   for (ii = 0; ii < out_data_stride * nSAvec; ii++) out_data[ii] = PETSC_MAX_REAL;
311:   *a_data_out = out_data; /* output - stride nSelected*nSAvec */

313:   /* find points and set prolongation */
314:   minsz = 100;
315:   for (mm = clid = 0; mm < nloc; mm++) {
316:     PetscCDSizeAt(agg_llists, mm, &jj);
317:     if (jj > 0) {
318:       const PetscInt lid = mm, cgid = my0crs + clid;
319:       PetscInt       cids[100]; /* max bs */
320:       PetscBLASInt   asz = jj, M = asz * bs, N = nSAvec, INFO;
321:       PetscBLASInt   Mdata = M + ((N - M > 0) ? N - M : 0), LDA = Mdata, LWORK = N * bs;
322:       PetscScalar   *qqc, *qqr, *TAU, *WORK;
323:       PetscInt      *fids;
324:       PetscReal     *data;

326:       /* count agg */
327:       if (asz < minsz) minsz = asz;

329:       /* get block */
330:       PetscMalloc5(Mdata * N, &qqc, M * N, &qqr, N, &TAU, LWORK, &WORK, M, &fids);

332:       aggID = 0;
333:       PetscCDGetHeadPos(agg_llists, lid, &pos);
334:       while (pos) {
335:         PetscInt gid1;
336:         PetscCDIntNdGetID(pos, &gid1);
337:         PetscCDGetNextPos(agg_llists, lid, &pos);

339:         if (gid1 >= my0 && gid1 < Iend) flid = gid1 - my0;
340:         else {
341:           PCGAMGHashTableFind(&fgid_flid, gid1, &flid);
343:         }
344:         /* copy in B_i matrix - column oriented */
345:         data = &data_in[flid * bs];
346:         for (ii = 0; ii < bs; ii++) {
347:           for (jj = 0; jj < N; jj++) {
348:             PetscReal d                       = data[jj * data_stride + ii];
349:             qqc[jj * Mdata + aggID * bs + ii] = d;
350:           }
351:         }
352:         /* set fine IDs */
353:         for (kk = 0; kk < bs; kk++) fids[aggID * bs + kk] = flid_fgid[flid] * bs + kk;
354:         aggID++;
355:       }

357:       /* pad with zeros */
358:       for (ii = asz * bs; ii < Mdata; ii++) {
359:         for (jj = 0; jj < N; jj++, kk++) qqc[jj * Mdata + ii] = .0;
360:       }

362:       /* QR */
363:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
364:       PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&Mdata, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO));
365:       PetscFPTrapPop();
367:       /* get R - column oriented - output B_{i+1} */
368:       {
369:         PetscReal *data = &out_data[clid * nSAvec];
370:         for (jj = 0; jj < nSAvec; jj++) {
371:           for (ii = 0; ii < nSAvec; ii++) {
373:             if (ii <= jj) data[jj * out_data_stride + ii] = PetscRealPart(qqc[jj * Mdata + ii]);
374:             else data[jj * out_data_stride + ii] = 0.;
375:           }
376:         }
377:       }

379:       /* get Q - row oriented */
380:       PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&Mdata, &N, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO));

383:       for (ii = 0; ii < M; ii++) {
384:         for (jj = 0; jj < N; jj++) qqr[N * ii + jj] = qqc[jj * Mdata + ii];
385:       }

387:       /* add diagonal block of P0 */
388:       for (kk = 0; kk < N; kk++) { cids[kk] = N * cgid + kk; /* global col IDs in P0 */ }
389:       MatSetValues(a_Prol, M, fids, N, cids, qqr, INSERT_VALUES);
390:       PetscFree5(qqc, qqr, TAU, WORK, fids);
391:       clid++;
392:     } /* coarse agg */
393:   }   /* for all fine nodes */
394:   MatAssemblyBegin(a_Prol, MAT_FINAL_ASSEMBLY);
395:   MatAssemblyEnd(a_Prol, MAT_FINAL_ASSEMBLY);
396:   PCGAMGHashTableDestroy(&fgid_flid);
397:   return 0;
398: }

400: static PetscErrorCode PCView_GAMG_AGG(PC pc, PetscViewer viewer)
401: {
402:   PC_MG       *mg          = (PC_MG *)pc->data;
403:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
404:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;

406:   PetscViewerASCIIPrintf(viewer, "      AGG specific options\n");
407:   PetscViewerASCIIPrintf(viewer, "        Number of levels to square graph %d\n", (int)pc_gamg_agg->aggressive_coarsening_levels);
408:   PetscViewerASCIIPrintf(viewer, "        Number smoothing steps %d\n", (int)pc_gamg_agg->nsmooths);
409:   return 0;
410: }

412: static PetscErrorCode PCGAMGCreateGraph_AGG(PC pc, Mat Amat, Mat *a_Gmat)
413: {
414:   PC_MG          *mg          = (PC_MG *)pc->data;
415:   PC_GAMG        *pc_gamg     = (PC_GAMG *)mg->innerctx;
416:   PC_GAMG_AGG    *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
417:   const PetscReal vfilter     = pc_gamg->threshold[pc_gamg->current_level];
418:   PetscBool       ishem;
419:   const char     *prefix;

421:   PetscLogEventBegin(petsc_gamg_setup_events[GAMG_GRAPH], 0, 0, 0, 0);
422:   /* Note: depending on the algorithm that will be used for computing the coarse grid points this should pass PETSC_TRUE or PETSC_FALSE as the first argument */

424:   /* MATCOARSENHEM requires numerical weights for edges so ensure they are computed */
425:   MatCoarsenCreate(PetscObjectComm((PetscObject)pc), &pc_gamg_agg->crs);
426:   PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix);
427:   PetscObjectSetOptionsPrefix((PetscObject)pc_gamg_agg->crs, prefix);
428:   MatCoarsenSetFromOptions(pc_gamg_agg->crs);
429:   PetscObjectTypeCompare((PetscObject)pc_gamg_agg->crs, MATCOARSENHEM, &ishem);

431:   MatCreateGraph(Amat, PETSC_TRUE, (vfilter >= 0 || ishem) ? PETSC_TRUE : PETSC_FALSE, vfilter, a_Gmat);
432:   PetscLogEventEnd(petsc_gamg_setup_events[GAMG_GRAPH], 0, 0, 0, 0);
433:   return 0;
434: }

436: /*
437:    PCGAMGCoarsen_AGG - supports squaring the graph (deprecated) and new graph for
438:      communication of QR data used with HEM and MISk coarsening

440:   Input Parameter:
441:    . a_pc - this

443:   Input/Output Parameter:
444:    . a_Gmat1 - graph to coarsen (in), graph off processor edges for QR gather scatter (out)

446:   Output Parameter:
447:    . agg_lists - list of aggregates

449: */
450: static PetscErrorCode PCGAMGCoarsen_AGG(PC a_pc, Mat *a_Gmat1, PetscCoarsenData **agg_lists)
451: {
452:   PC_MG       *mg          = (PC_MG *)a_pc->data;
453:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
454:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
455:   Mat          mat, Gmat1 = *a_Gmat1; /* aggressive graph */
456:   IS           perm;
457:   PetscInt     Istart, Iend, Ii, nloc, bs, nn;
458:   PetscInt    *permute, *degree;
459:   PetscBool   *bIndexSet;
460:   PetscReal    hashfact;
461:   PetscInt     iSwapIndex;
462:   PetscRandom  random;

464:   PetscLogEventBegin(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0);
465:   MatGetLocalSize(Gmat1, &nn, NULL);
466:   MatGetBlockSize(Gmat1, &bs);
468:   nloc = nn / bs;

470:   /* get MIS aggs - randomize */
471:   PetscMalloc2(nloc, &permute, nloc, &degree);
472:   PetscCalloc1(nloc, &bIndexSet);
473:   for (Ii = 0; Ii < nloc; Ii++) permute[Ii] = Ii;
474:   PetscRandomCreate(PETSC_COMM_SELF, &random);
475:   MatGetOwnershipRange(Gmat1, &Istart, &Iend);
476:   for (Ii = 0; Ii < nloc; Ii++) {
477:     PetscInt nc;
478:     MatGetRow(Gmat1, Istart + Ii, &nc, NULL, NULL);
479:     degree[Ii] = nc;
480:     MatRestoreRow(Gmat1, Istart + Ii, &nc, NULL, NULL);
481:   }
482:   for (Ii = 0; Ii < nloc; Ii++) {
483:     PetscRandomGetValueReal(random, &hashfact);
484:     iSwapIndex = (PetscInt)(hashfact * nloc) % nloc;
485:     if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) {
486:       PetscInt iTemp        = permute[iSwapIndex];
487:       permute[iSwapIndex]   = permute[Ii];
488:       permute[Ii]           = iTemp;
489:       iTemp                 = degree[iSwapIndex];
490:       degree[iSwapIndex]    = degree[Ii];
491:       degree[Ii]            = iTemp;
492:       bIndexSet[iSwapIndex] = PETSC_TRUE;
493:     }
494:   }
495:   // create minimum degree ordering
496:   PetscSortIntWithArray(nloc, degree, permute);

498:   PetscFree(bIndexSet);
499:   PetscRandomDestroy(&random);
500:   ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_USE_POINTER, &perm);
501:   PetscLogEventBegin(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0);
502:   MatCoarsenSetGreedyOrdering(pc_gamg_agg->crs, perm);
503:   MatCoarsenSetAdjacency(pc_gamg_agg->crs, Gmat1);
504:   MatCoarsenSetStrictAggs(pc_gamg_agg->crs, PETSC_TRUE);
505:   if (pc_gamg->current_level < pc_gamg_agg->aggressive_coarsening_levels) MatCoarsenMISKSetDistance(pc_gamg_agg->crs, 2); // hardwire to MIS-2
506:   else MatCoarsenMISKSetDistance(pc_gamg_agg->crs, 1);                                                                    // MIS
507:   MatCoarsenApply(pc_gamg_agg->crs);
508:   MatCoarsenViewFromOptions(pc_gamg_agg->crs, NULL, "-mat_coarsen_view");
509:   MatCoarsenGetData(pc_gamg_agg->crs, agg_lists); /* output */
510:   MatCoarsenDestroy(&pc_gamg_agg->crs);

512:   ISDestroy(&perm);
513:   PetscFree2(permute, degree);
514:   PetscLogEventEnd(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0);

516:   {
517:     PetscCoarsenData *llist = *agg_lists;
518:     /* see if we have a matrix that takes precedence (returned from MatCoarsenApply) */
519:     PetscCDGetMat(llist, &mat);
520:     if (mat) {
521:       MatDestroy(&Gmat1);
522:       *a_Gmat1 = mat; /* output */
523:     }
524:   }
525:   PetscLogEventEnd(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0);
526:   return 0;
527: }

529: /*
530:  PCGAMGProlongator_AGG

532:  Input Parameter:
533:  . pc - this
534:  . Amat - matrix on this fine level
535:  . Graph - used to get ghost data for nodes in
536:  . agg_lists - list of aggregates
537:  Output Parameter:
538:  . a_P_out - prolongation operator to the next level
539:  */
540: static PetscErrorCode PCGAMGProlongator_AGG(PC pc, Mat Amat, Mat Gmat, PetscCoarsenData *agg_lists, Mat *a_P_out)
541: {
542:   PC_MG         *mg      = (PC_MG *)pc->data;
543:   PC_GAMG       *pc_gamg = (PC_GAMG *)mg->innerctx;
544:   const PetscInt col_bs  = pc_gamg->data_cell_cols;
545:   PetscInt       Istart, Iend, nloc, ii, jj, kk, my0, nLocalSelected, bs;
546:   Mat            Prol;
547:   PetscMPIInt    size;
548:   MPI_Comm       comm;
549:   PetscReal     *data_w_ghost;
550:   PetscInt       myCrs0, nbnodes = 0, *flid_fgid;
551:   MatType        mtype;

553:   PetscObjectGetComm((PetscObject)Amat, &comm);
555:   PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0);
556:   MPI_Comm_size(comm, &size);
557:   MatGetOwnershipRange(Amat, &Istart, &Iend);
558:   MatGetBlockSize(Amat, &bs);
559:   nloc = (Iend - Istart) / bs;
560:   my0  = Istart / bs;

563:   /* get 'nLocalSelected' */
564:   for (ii = 0, nLocalSelected = 0; ii < nloc; ii++) {
565:     PetscBool ise;
566:     /* filter out singletons 0 or 1? */
567:     PetscCDEmptyAt(agg_lists, ii, &ise);
568:     if (!ise) nLocalSelected++;
569:   }

571:   /* create prolongator, create P matrix */
572:   MatGetType(Amat, &mtype);
573:   MatCreate(comm, &Prol);
574:   MatSetSizes(Prol, nloc * bs, nLocalSelected * col_bs, PETSC_DETERMINE, PETSC_DETERMINE);
575:   MatSetBlockSizes(Prol, bs, col_bs);
576:   MatSetType(Prol, mtype);
577:   MatSeqAIJSetPreallocation(Prol, col_bs, NULL);
578:   MatMPIAIJSetPreallocation(Prol, col_bs, NULL, col_bs, NULL);

580:   /* can get all points "removed" */
581:   MatGetSize(Prol, &kk, &ii);
582:   if (!ii) {
583:     PetscInfo(pc, "%s: No selected points on coarse grid\n", ((PetscObject)pc)->prefix);
584:     MatDestroy(&Prol);
585:     *a_P_out = NULL; /* out */
586:     PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0);
587:     return 0;
588:   }
589:   PetscInfo(pc, "%s: New grid %" PetscInt_FMT " nodes\n", ((PetscObject)pc)->prefix, ii / col_bs);
590:   MatGetOwnershipRangeColumn(Prol, &myCrs0, &kk);

593:   myCrs0 = myCrs0 / col_bs;

596:   /* create global vector of data in 'data_w_ghost' */
597:   PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0);
598:   if (size > 1) { /* get ghost null space data */
599:     PetscReal *tmp_gdata, *tmp_ldata, *tp2;
600:     PetscMalloc1(nloc, &tmp_ldata);
601:     for (jj = 0; jj < col_bs; jj++) {
602:       for (kk = 0; kk < bs; kk++) {
603:         PetscInt         ii, stride;
604:         const PetscReal *tp = pc_gamg->data + jj * bs * nloc + kk;
605:         for (ii = 0; ii < nloc; ii++, tp += bs) tmp_ldata[ii] = *tp;

607:         PCGAMGGetDataWithGhosts(Gmat, 1, tmp_ldata, &stride, &tmp_gdata);

609:         if (!jj && !kk) { /* now I know how many total nodes - allocate TODO: move below and do in one 'col_bs' call */
610:           PetscMalloc1(stride * bs * col_bs, &data_w_ghost);
611:           nbnodes = bs * stride;
612:         }
613:         tp2 = data_w_ghost + jj * bs * stride + kk;
614:         for (ii = 0; ii < stride; ii++, tp2 += bs) *tp2 = tmp_gdata[ii];
615:         PetscFree(tmp_gdata);
616:       }
617:     }
618:     PetscFree(tmp_ldata);
619:   } else {
620:     nbnodes      = bs * nloc;
621:     data_w_ghost = (PetscReal *)pc_gamg->data;
622:   }

624:   /* get 'flid_fgid' TODO - move up to get 'stride' and do get null space data above in one step (jj loop) */
625:   if (size > 1) {
626:     PetscReal *fid_glid_loc, *fiddata;
627:     PetscInt   stride;

629:     PetscMalloc1(nloc, &fid_glid_loc);
630:     for (kk = 0; kk < nloc; kk++) fid_glid_loc[kk] = (PetscReal)(my0 + kk);
631:     PCGAMGGetDataWithGhosts(Gmat, 1, fid_glid_loc, &stride, &fiddata);
632:     PetscMalloc1(stride, &flid_fgid); /* copy real data to in */
633:     for (kk = 0; kk < stride; kk++) flid_fgid[kk] = (PetscInt)fiddata[kk];
634:     PetscFree(fiddata);

637:     PetscFree(fid_glid_loc);
638:   } else {
639:     PetscMalloc1(nloc, &flid_fgid);
640:     for (kk = 0; kk < nloc; kk++) flid_fgid[kk] = my0 + kk;
641:   }
642:   PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0);
643:   /* get P0 */
644:   PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0);
645:   {
646:     PetscReal *data_out = NULL;
647:     formProl0(agg_lists, bs, col_bs, myCrs0, nbnodes, data_w_ghost, flid_fgid, &data_out, Prol);
648:     PetscFree(pc_gamg->data);

650:     pc_gamg->data           = data_out;
651:     pc_gamg->data_cell_rows = col_bs;
652:     pc_gamg->data_sz        = col_bs * col_bs * nLocalSelected;
653:   }
654:   PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0);
655:   if (size > 1) PetscFree(data_w_ghost);
656:   PetscFree(flid_fgid);

658:   *a_P_out = Prol; /* out */

660:   PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0);
661:   return 0;
662: }

664: /*
665:    PCGAMGOptProlongator_AGG

667:   Input Parameter:
668:    . pc - this
669:    . Amat - matrix on this fine level
670:  In/Output Parameter:
671:    . a_P - prolongation operator to the next level
672: */
673: static PetscErrorCode PCGAMGOptProlongator_AGG(PC pc, Mat Amat, Mat *a_P)
674: {
675:   PC_MG       *mg          = (PC_MG *)pc->data;
676:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
677:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
678:   PetscInt     jj;
679:   Mat          Prol = *a_P;
680:   MPI_Comm     comm;
681:   KSP          eksp;
682:   Vec          bb, xx;
683:   PC           epc;
684:   PetscReal    alpha, emax, emin;

686:   PetscObjectGetComm((PetscObject)Amat, &comm);
687:   PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0);

689:   /* compute maximum singular value of operator to be used in smoother */
690:   if (0 < pc_gamg_agg->nsmooths) {
691:     /* get eigen estimates */
692:     if (pc_gamg->emax > 0) {
693:       emin = pc_gamg->emin;
694:       emax = pc_gamg->emax;
695:     } else {
696:       const char *prefix;

698:       MatCreateVecs(Amat, &bb, NULL);
699:       MatCreateVecs(Amat, &xx, NULL);
700:       KSPSetNoisy_Private(bb);

702:       KSPCreate(comm, &eksp);
703:       PCGetOptionsPrefix(pc, &prefix);
704:       KSPSetOptionsPrefix(eksp, prefix);
705:       KSPAppendOptionsPrefix(eksp, "pc_gamg_esteig_");
706:       {
707:         PetscBool isset, sflg;
708:         MatIsSPDKnown(Amat, &isset, &sflg);
709:         if (isset && sflg) KSPSetType(eksp, KSPCG);
710:       }
711:       KSPSetErrorIfNotConverged(eksp, pc->erroriffailure);
712:       KSPSetNormType(eksp, KSP_NORM_NONE);

714:       KSPSetInitialGuessNonzero(eksp, PETSC_FALSE);
715:       KSPSetOperators(eksp, Amat, Amat);

717:       KSPGetPC(eksp, &epc);
718:       PCSetType(epc, PCJACOBI); /* smoother in smoothed agg. */

720:       KSPSetTolerances(eksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, 10); // 10 is safer, but 5 is often fine, can override with -pc_gamg_esteig_ksp_max_it -mg_levels_ksp_chebyshev_esteig 0,0.25,0,1.2

722:       KSPSetFromOptions(eksp);
723:       KSPSetComputeSingularValues(eksp, PETSC_TRUE);
724:       KSPSolve(eksp, bb, xx);
725:       KSPCheckSolve(eksp, pc, xx);

727:       KSPComputeExtremeSingularValues(eksp, &emax, &emin);
728:       PetscInfo(pc, "%s: Smooth P0: max eigen=%e min=%e PC=%s\n", ((PetscObject)pc)->prefix, (double)emax, (double)emin, PCJACOBI);
729:       VecDestroy(&xx);
730:       VecDestroy(&bb);
731:       KSPDestroy(&eksp);
732:     }
733:     if (pc_gamg->use_sa_esteig) {
734:       mg->min_eigen_DinvA[pc_gamg->current_level] = emin;
735:       mg->max_eigen_DinvA[pc_gamg->current_level] = emax;
736:       PetscInfo(pc, "%s: Smooth P0: level %" PetscInt_FMT ", cache spectra %g %g\n", ((PetscObject)pc)->prefix, pc_gamg->current_level, (double)emin, (double)emax);
737:     } else {
738:       mg->min_eigen_DinvA[pc_gamg->current_level] = 0;
739:       mg->max_eigen_DinvA[pc_gamg->current_level] = 0;
740:     }
741:   } else {
742:     mg->min_eigen_DinvA[pc_gamg->current_level] = 0;
743:     mg->max_eigen_DinvA[pc_gamg->current_level] = 0;
744:   }

746:   /* smooth P0 */
747:   for (jj = 0; jj < pc_gamg_agg->nsmooths; jj++) {
748:     Mat tMat;
749:     Vec diag;

751:     PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0);

753:     /* smooth P1 := (I - omega/lam D^{-1}A)P0 */
754:     PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0);
755:     MatMatMult(Amat, Prol, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &tMat);
756:     PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0);
757:     MatProductClear(tMat);
758:     MatCreateVecs(Amat, &diag, NULL);
759:     MatGetDiagonal(Amat, diag); /* effectively PCJACOBI */
760:     VecReciprocal(diag);
761:     MatDiagonalScale(tMat, diag, NULL);
762:     VecDestroy(&diag);

764:     /* TODO: Set a PCFailedReason and exit the building of the AMG preconditioner */
766:     /* TODO: Document the 1.4 and don't hardwire it in this routine */
767:     alpha = -1.4 / emax;

769:     MatAYPX(tMat, alpha, Prol, SUBSET_NONZERO_PATTERN);
770:     MatDestroy(&Prol);
771:     Prol = tMat;
772:     PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0);
773:   }
774:   PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0);
775:   *a_P = Prol;
776:   return 0;
777: }

779: /*
780:    PCCreateGAMG_AGG

782:   Input Parameter:
783:    . pc -
784: */
785: PetscErrorCode PCCreateGAMG_AGG(PC pc)
786: {
787:   PC_MG       *mg      = (PC_MG *)pc->data;
788:   PC_GAMG     *pc_gamg = (PC_GAMG *)mg->innerctx;
789:   PC_GAMG_AGG *pc_gamg_agg;

791:   /* create sub context for SA */
792:   PetscNew(&pc_gamg_agg);
793:   pc_gamg->subctx = pc_gamg_agg;

795:   pc_gamg->ops->setfromoptions = PCSetFromOptions_GAMG_AGG;
796:   pc_gamg->ops->destroy        = PCDestroy_GAMG_AGG;
797:   /* reset does not do anything; setup not virtual */

799:   /* set internal function pointers */
800:   pc_gamg->ops->creategraph       = PCGAMGCreateGraph_AGG;
801:   pc_gamg->ops->coarsen           = PCGAMGCoarsen_AGG;
802:   pc_gamg->ops->prolongator       = PCGAMGProlongator_AGG;
803:   pc_gamg->ops->optprolongator    = PCGAMGOptProlongator_AGG;
804:   pc_gamg->ops->createdefaultdata = PCSetData_AGG;
805:   pc_gamg->ops->view              = PCView_GAMG_AGG;

807:   pc_gamg_agg->aggressive_coarsening_levels = 0;
808:   pc_gamg_agg->nsmooths                     = 1;

810:   PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", PCGAMGSetNSmooths_AGG);
811:   PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", PCGAMGSetAggressiveLevels_AGG);
812:   PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_AGG);
813:   return 0;
814: }