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: /*@
 14:   PCGAMGClassicalSetType - Sets the type of classical interpolation to use with `PCGAMG`

 16:   Collective

 18:   Input Parameters:
 19: + pc   - the preconditioner context
 20: - type - the interpolation to use, see `PCGAMGClassicalType()`

 22:   Options Database Key:
 23: . -pc_gamg_classical_type <direct,standard> - set type of classical AMG prolongation

 25:   Level: intermediate

 27: .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGClassicalType`, `PCGAMGClassicalGetType()`
 28: @*/
 29: PetscErrorCode PCGAMGClassicalSetType(PC pc, PCGAMGClassicalType type)
 30: {
 31:   PetscFunctionBegin;
 33:   PetscTryMethod(pc, "PCGAMGClassicalSetType_C", (PC, PCGAMGClassicalType), (pc, type));
 34:   PetscFunctionReturn(PETSC_SUCCESS);
 35: }

 37: /*@
 38:   PCGAMGClassicalGetType - Gets the type of classical interpolation to use with `PCGAMG`

 40:   Collective

 42:   Input Parameter:
 43: . pc - the preconditioner context

 45:   Output Parameter:
 46: . type - the type used, see `PCGAMGClassicalType()`

 48:   Level: intermediate

 50: .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGClassicalType`, `PCGAMGClassicalSetType()`
 51: @*/
 52: PetscErrorCode PCGAMGClassicalGetType(PC pc, PCGAMGClassicalType *type)
 53: {
 54:   PetscFunctionBegin;
 56:   PetscUseMethod(pc, "PCGAMGClassicalGetType_C", (PC, PCGAMGClassicalType *), (pc, type));
 57:   PetscFunctionReturn(PETSC_SUCCESS);
 58: }

 60: static PetscErrorCode PCGAMGClassicalSetType_GAMG(PC pc, PCGAMGClassicalType type)
 61: {
 62:   PC_MG             *mg      = (PC_MG *)pc->data;
 63:   PC_GAMG           *pc_gamg = (PC_GAMG *)mg->innerctx;
 64:   PC_GAMG_Classical *cls     = (PC_GAMG_Classical *)pc_gamg->subctx;

 66:   PetscFunctionBegin;
 67:   PetscCall(PetscStrncpy(cls->prolongtype, type, sizeof(cls->prolongtype)));
 68:   PetscFunctionReturn(PETSC_SUCCESS);
 69: }

 71: static PetscErrorCode PCGAMGClassicalGetType_GAMG(PC pc, PCGAMGClassicalType *type)
 72: {
 73:   PC_MG             *mg      = (PC_MG *)pc->data;
 74:   PC_GAMG           *pc_gamg = (PC_GAMG *)mg->innerctx;
 75:   PC_GAMG_Classical *cls     = (PC_GAMG_Classical *)pc_gamg->subctx;

 77:   PetscFunctionBegin;
 78:   *type = cls->prolongtype;
 79:   PetscFunctionReturn(PETSC_SUCCESS);
 80: }

 82: static PetscErrorCode PCGAMGCreateGraph_Classical(PC pc, Mat A, Mat *G)
 83: {
 84:   PetscInt           s, f, n, idx, lidx, gidx;
 85:   PetscInt           r, c, ncols;
 86:   const PetscInt    *rcol;
 87:   const PetscScalar *rval;
 88:   PetscInt          *gcol;
 89:   PetscScalar       *gval;
 90:   PetscReal          rmax;
 91:   PetscInt           cmax = 0;
 92:   PC_MG             *mg   = (PC_MG *)pc->data;
 93:   PC_GAMG           *gamg = (PC_GAMG *)mg->innerctx;
 94:   PetscInt          *gsparse, *lsparse;
 95:   PetscScalar       *Amax;
 96:   MatType            mtype;

 98:   PetscFunctionBegin;
 99:   PetscCall(MatGetOwnershipRange(A, &s, &f));
100:   n = f - s;
101:   PetscCall(PetscMalloc3(n, &lsparse, n, &gsparse, n, &Amax));

103:   for (r = 0; r < n; r++) {
104:     lsparse[r] = 0;
105:     gsparse[r] = 0;
106:   }

108:   for (r = s; r < f; r++) {
109:     /* determine the maximum off-diagonal in each row */
110:     rmax = 0.;
111:     PetscCall(MatGetRow(A, r, &ncols, &rcol, &rval));
112:     for (c = 0; c < ncols; c++) {
113:       if (PetscRealPart(-rval[c]) > rmax && rcol[c] != r) rmax = PetscRealPart(-rval[c]);
114:     }
115:     Amax[r - s] = rmax;
116:     if (ncols > cmax) cmax = ncols;
117:     lidx = 0;
118:     gidx = 0;
119:     /* create the local and global sparsity patterns */
120:     for (c = 0; c < ncols; c++) {
121:       if (PetscRealPart(-rval[c]) > gamg->threshold[0] * PetscRealPart(Amax[r - s]) || rcol[c] == r) {
122:         if (rcol[c] < f && rcol[c] >= s) {
123:           lidx++;
124:         } else {
125:           gidx++;
126:         }
127:       }
128:     }
129:     PetscCall(MatRestoreRow(A, r, &ncols, &rcol, &rval));
130:     lsparse[r - s] = lidx;
131:     gsparse[r - s] = gidx;
132:   }
133:   PetscCall(PetscMalloc2(cmax, &gval, cmax, &gcol));

135:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), G));
136:   PetscCall(MatGetType(A, &mtype));
137:   PetscCall(MatSetType(*G, mtype));
138:   PetscCall(MatSetSizes(*G, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
139:   PetscCall(MatMPIAIJSetPreallocation(*G, 0, lsparse, 0, gsparse));
140:   PetscCall(MatSeqAIJSetPreallocation(*G, 0, lsparse));
141:   for (r = s; r < f; r++) {
142:     PetscCall(MatGetRow(A, r, &ncols, &rcol, &rval));
143:     idx = 0;
144:     for (c = 0; c < ncols; c++) {
145:       /* classical strength of connection */
146:       if (PetscRealPart(-rval[c]) > gamg->threshold[0] * PetscRealPart(Amax[r - s]) || rcol[c] == r) {
147:         gcol[idx] = rcol[c];
148:         gval[idx] = rval[c];
149:         idx++;
150:       }
151:     }
152:     PetscCall(MatSetValues(*G, 1, &r, idx, gcol, gval, INSERT_VALUES));
153:     PetscCall(MatRestoreRow(A, r, &ncols, &rcol, &rval));
154:   }
155:   PetscCall(MatAssemblyBegin(*G, MAT_FINAL_ASSEMBLY));
156:   PetscCall(MatAssemblyEnd(*G, MAT_FINAL_ASSEMBLY));

158:   PetscCall(PetscFree2(gval, gcol));
159:   PetscCall(PetscFree3(lsparse, gsparse, Amax));
160:   PetscFunctionReturn(PETSC_SUCCESS);
161: }

163: static PetscErrorCode PCGAMGCoarsen_Classical(PC pc, Mat *G, PetscCoarsenData **agg_lists)
164: {
165:   MatCoarsen  crs;
166:   MPI_Comm    fcomm = ((PetscObject)pc)->comm;
167:   const char *prefix;

169:   PetscFunctionBegin;
170:   PetscCheck(G, fcomm, PETSC_ERR_ARG_WRONGSTATE, "Must set Graph in PC in PCGAMG before coarsening");

172:   PetscCall(MatCoarsenCreate(fcomm, &crs));
173:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix));
174:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)crs, prefix));
175:   PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)crs, "pc_gamg_"));
176:   PetscCall(MatCoarsenSetFromOptions(crs));
177:   PetscCall(MatCoarsenSetAdjacency(crs, *G));
178:   PetscCall(MatCoarsenSetStrictAggs(crs, PETSC_TRUE));
179:   PetscCall(MatCoarsenApply(crs));
180:   PetscCall(MatCoarsenGetData(crs, agg_lists));
181:   PetscCall(MatCoarsenDestroy(&crs));
182:   PetscFunctionReturn(PETSC_SUCCESS);
183: }

185: static PetscErrorCode PCGAMGProlongator_Classical_Direct(PC pc, Mat A, PetscCoarsenData *agg_lists, Mat *P)
186: {
187:   PC_MG             *mg   = (PC_MG *)pc->data;
188:   PC_GAMG           *gamg = (PC_GAMG *)mg->innerctx;
189:   PetscBool          iscoarse, isMPIAIJ, isSEQAIJ;
190:   PetscInt           fn, cn, fs, fe, cs, ce, i, j, ncols, col, row_f, row_c, cmax = 0, idx, noff;
191:   PetscInt          *lcid, *gcid, *lsparse, *gsparse, *colmap, *pcols;
192:   const PetscInt    *rcol;
193:   PetscReal         *Amax_pos, *Amax_neg;
194:   PetscScalar        g_pos, g_neg, a_pos, a_neg, diag, invdiag, alpha, beta, pij;
195:   PetscScalar       *pvals;
196:   const PetscScalar *rval;
197:   Mat                lA, gA = NULL;
198:   MatType            mtype;
199:   Vec                C, lvec;
200:   PetscLayout        clayout;
201:   PetscSF            sf;
202:   Mat_MPIAIJ        *mpiaij;

204:   PetscFunctionBegin;
205:   PetscCall(MatGetOwnershipRange(A, &fs, &fe));
206:   fn = fe - fs;
207:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPIAIJ, &isMPIAIJ));
208:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQAIJ, &isSEQAIJ));
209:   PetscCheck(isMPIAIJ || isSEQAIJ, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Classical AMG requires MPIAIJ matrix");
210:   if (isMPIAIJ) {
211:     mpiaij = (Mat_MPIAIJ *)A->data;
212:     lA     = mpiaij->A;
213:     gA     = mpiaij->B;
214:     lvec   = mpiaij->lvec;
215:     PetscCall(VecGetSize(lvec, &noff));
216:     colmap = mpiaij->garray;
217:     PetscCall(MatGetLayouts(A, NULL, &clayout));
218:     PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)A), &sf));
219:     PetscCall(PetscSFSetGraphLayout(sf, clayout, noff, NULL, PETSC_COPY_VALUES, colmap));
220:     PetscCall(PetscMalloc1(noff, &gcid));
221:   } else {
222:     lA = A;
223:   }
224:   PetscCall(PetscMalloc5(fn, &lsparse, fn, &gsparse, fn, &lcid, fn, &Amax_pos, fn, &Amax_neg));

226:   /* count the number of coarse unknowns */
227:   cn = 0;
228:   for (i = 0; i < fn; i++) {
229:     /* filter out singletons */
230:     PetscCall(PetscCDIsEmptyAt(agg_lists, i, &iscoarse));
231:     lcid[i] = -1;
232:     if (!iscoarse) cn++;
233:   }

235:   /* create the coarse vector */
236:   PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)A), cn, PETSC_DECIDE, &C));
237:   PetscCall(VecGetOwnershipRange(C, &cs, &ce));

239:   cn = 0;
240:   for (i = 0; i < fn; i++) {
241:     PetscCall(PetscCDIsEmptyAt(agg_lists, i, &iscoarse));
242:     if (!iscoarse) {
243:       lcid[i] = cs + cn;
244:       cn++;
245:     } else {
246:       lcid[i] = -1;
247:     }
248:   }

250:   if (gA) {
251:     PetscCall(PetscSFBcastBegin(sf, MPIU_INT, lcid, gcid, MPI_REPLACE));
252:     PetscCall(PetscSFBcastEnd(sf, MPIU_INT, lcid, gcid, MPI_REPLACE));
253:   }

255:   /* determine the largest off-diagonal entries in each row */
256:   for (i = fs; i < fe; i++) {
257:     Amax_pos[i - fs] = 0.;
258:     Amax_neg[i - fs] = 0.;
259:     PetscCall(MatGetRow(A, i, &ncols, &rcol, &rval));
260:     for (j = 0; j < ncols; j++) {
261:       if ((PetscRealPart(-rval[j]) > Amax_neg[i - fs]) && i != rcol[j]) Amax_neg[i - fs] = PetscAbsScalar(rval[j]);
262:       if ((PetscRealPart(rval[j]) > Amax_pos[i - fs]) && i != rcol[j]) Amax_pos[i - fs] = PetscAbsScalar(rval[j]);
263:     }
264:     if (ncols > cmax) cmax = ncols;
265:     PetscCall(MatRestoreRow(A, i, &ncols, &rcol, &rval));
266:   }
267:   PetscCall(PetscMalloc2(cmax, &pcols, cmax, &pvals));
268:   PetscCall(VecDestroy(&C));

270:   /* count the on and off processor sparsity patterns for the prolongator */
271:   for (i = 0; i < fn; i++) {
272:     /* on */
273:     lsparse[i] = 0;
274:     gsparse[i] = 0;
275:     if (lcid[i] >= 0) {
276:       lsparse[i] = 1;
277:       gsparse[i] = 0;
278:     } else {
279:       PetscCall(MatGetRow(lA, i, &ncols, &rcol, &rval));
280:       for (j = 0; j < ncols; j++) {
281:         col = rcol[j];
282:         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;
283:       }
284:       PetscCall(MatRestoreRow(lA, i, &ncols, &rcol, &rval));
285:       /* off */
286:       if (gA) {
287:         PetscCall(MatGetRow(gA, i, &ncols, &rcol, &rval));
288:         for (j = 0; j < ncols; j++) {
289:           col = rcol[j];
290:           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;
291:         }
292:         PetscCall(MatRestoreRow(gA, i, &ncols, &rcol, &rval));
293:       }
294:     }
295:   }

297:   /* preallocate and create the prolongator */
298:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), P));
299:   PetscCall(MatGetType(A, &mtype));
300:   PetscCall(MatSetType(*P, mtype));
301:   PetscCall(MatSetSizes(*P, fn, cn, PETSC_DETERMINE, PETSC_DETERMINE));
302:   PetscCall(MatMPIAIJSetPreallocation(*P, 0, lsparse, 0, gsparse));
303:   PetscCall(MatSeqAIJSetPreallocation(*P, 0, lsparse));

305:   /* loop over local fine nodes -- get the diagonal, the sum of positive and negative strong and weak weights, and set up the row */
306:   for (i = 0; i < fn; i++) {
307:     /* determine on or off */
308:     row_f = i + fs;
309:     row_c = lcid[i];
310:     if (row_c >= 0) {
311:       pij = 1.;
312:       PetscCall(MatSetValues(*P, 1, &row_f, 1, &row_c, &pij, INSERT_VALUES));
313:     } else {
314:       g_pos = 0.;
315:       g_neg = 0.;
316:       a_pos = 0.;
317:       a_neg = 0.;
318:       diag  = 0.;

320:       /* local connections */
321:       PetscCall(MatGetRow(lA, i, &ncols, &rcol, &rval));
322:       for (j = 0; j < ncols; j++) {
323:         col = rcol[j];
324:         if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0] * Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0] * Amax_neg[i])) {
325:           if (PetscRealPart(rval[j]) > 0.) {
326:             g_pos += rval[j];
327:           } else {
328:             g_neg += rval[j];
329:           }
330:         }
331:         if (col != i) {
332:           if (PetscRealPart(rval[j]) > 0.) {
333:             a_pos += rval[j];
334:           } else {
335:             a_neg += rval[j];
336:           }
337:         } else {
338:           diag = rval[j];
339:         }
340:       }
341:       PetscCall(MatRestoreRow(lA, i, &ncols, &rcol, &rval));

343:       /* ghosted connections */
344:       if (gA) {
345:         PetscCall(MatGetRow(gA, i, &ncols, &rcol, &rval));
346:         for (j = 0; j < ncols; j++) {
347:           col = rcol[j];
348:           if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0] * Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0] * Amax_neg[i])) {
349:             if (PetscRealPart(rval[j]) > 0.) {
350:               g_pos += rval[j];
351:             } else {
352:               g_neg += rval[j];
353:             }
354:           }
355:           if (PetscRealPart(rval[j]) > 0.) {
356:             a_pos += rval[j];
357:           } else {
358:             a_neg += rval[j];
359:           }
360:         }
361:         PetscCall(MatRestoreRow(gA, i, &ncols, &rcol, &rval));
362:       }

364:       if (g_neg == 0.) {
365:         alpha = 0.;
366:       } else {
367:         alpha = -a_neg / g_neg;
368:       }

370:       if (g_pos == 0.) {
371:         diag += a_pos;
372:         beta = 0.;
373:       } else {
374:         beta = -a_pos / g_pos;
375:       }
376:       if (diag == 0.) {
377:         invdiag = 0.;
378:       } else invdiag = 1. / diag;
379:       /* on */
380:       PetscCall(MatGetRow(lA, i, &ncols, &rcol, &rval));
381:       idx = 0;
382:       for (j = 0; j < ncols; j++) {
383:         col = rcol[j];
384:         if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0] * Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0] * Amax_neg[i])) {
385:           row_f = i + fs;
386:           row_c = lcid[col];
387:           /* set the values for on-processor ones */
388:           if (PetscRealPart(rval[j]) < 0.) {
389:             pij = rval[j] * alpha * invdiag;
390:           } else {
391:             pij = rval[j] * beta * invdiag;
392:           }
393:           if (PetscAbsScalar(pij) != 0.) {
394:             pvals[idx] = pij;
395:             pcols[idx] = row_c;
396:             idx++;
397:           }
398:         }
399:       }
400:       PetscCall(MatRestoreRow(lA, i, &ncols, &rcol, &rval));
401:       /* off */
402:       if (gA) {
403:         PetscCall(MatGetRow(gA, i, &ncols, &rcol, &rval));
404:         for (j = 0; j < ncols; j++) {
405:           col = rcol[j];
406:           if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0] * Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0] * Amax_neg[i])) {
407:             row_f = i + fs;
408:             row_c = gcid[col];
409:             /* set the values for on-processor ones */
410:             if (PetscRealPart(rval[j]) < 0.) {
411:               pij = rval[j] * alpha * invdiag;
412:             } else {
413:               pij = rval[j] * beta * invdiag;
414:             }
415:             if (PetscAbsScalar(pij) != 0.) {
416:               pvals[idx] = pij;
417:               pcols[idx] = row_c;
418:               idx++;
419:             }
420:           }
421:         }
422:         PetscCall(MatRestoreRow(gA, i, &ncols, &rcol, &rval));
423:       }
424:       PetscCall(MatSetValues(*P, 1, &row_f, idx, pcols, pvals, INSERT_VALUES));
425:     }
426:   }

428:   PetscCall(MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY));
429:   PetscCall(MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY));

431:   PetscCall(PetscFree5(lsparse, gsparse, lcid, Amax_pos, Amax_neg));

433:   PetscCall(PetscFree2(pcols, pvals));
434:   if (gA) {
435:     PetscCall(PetscSFDestroy(&sf));
436:     PetscCall(PetscFree(gcid));
437:   }
438:   PetscFunctionReturn(PETSC_SUCCESS);
439: }

441: static PetscErrorCode PCGAMGTruncateProlongator_Private(PC pc, Mat *P)
442: {
443:   PetscInt           j, i, ps, pf, pn, pcs, pcf, pcn, idx, cmax;
444:   const PetscScalar *pval;
445:   const PetscInt    *pcol;
446:   PetscScalar       *pnval;
447:   PetscInt          *pncol;
448:   PetscInt           ncols;
449:   Mat                Pnew;
450:   PetscInt          *lsparse, *gsparse;
451:   PetscReal          pmax_pos, pmax_neg, ptot_pos, ptot_neg, pthresh_pos, pthresh_neg;
452:   PC_MG             *mg      = (PC_MG *)pc->data;
453:   PC_GAMG           *pc_gamg = (PC_GAMG *)mg->innerctx;
454:   PC_GAMG_Classical *cls     = (PC_GAMG_Classical *)pc_gamg->subctx;
455:   MatType            mtype;

457:   PetscFunctionBegin;
458:   /* trim and rescale with reallocation */
459:   PetscCall(MatGetOwnershipRange(*P, &ps, &pf));
460:   PetscCall(MatGetOwnershipRangeColumn(*P, &pcs, &pcf));
461:   pn  = pf - ps;
462:   pcn = pcf - pcs;
463:   PetscCall(PetscMalloc2(pn, &lsparse, pn, &gsparse));
464:   /* allocate */
465:   cmax = 0;
466:   for (i = ps; i < pf; i++) {
467:     lsparse[i - ps] = 0;
468:     gsparse[i - ps] = 0;
469:     PetscCall(MatGetRow(*P, i, &ncols, &pcol, &pval));
470:     if (ncols > cmax) cmax = ncols;
471:     pmax_pos = 0.;
472:     pmax_neg = 0.;
473:     for (j = 0; j < ncols; j++) {
474:       if (PetscRealPart(pval[j]) > pmax_pos) {
475:         pmax_pos = PetscRealPart(pval[j]);
476:       } else if (PetscRealPart(pval[j]) < pmax_neg) {
477:         pmax_neg = PetscRealPart(pval[j]);
478:       }
479:     }
480:     for (j = 0; j < ncols; j++) {
481:       if (PetscRealPart(pval[j]) >= pmax_pos * cls->interp_threshold || PetscRealPart(pval[j]) <= pmax_neg * cls->interp_threshold) {
482:         if (pcol[j] >= pcs && pcol[j] < pcf) {
483:           lsparse[i - ps]++;
484:         } else {
485:           gsparse[i - ps]++;
486:         }
487:       }
488:     }
489:     PetscCall(MatRestoreRow(*P, i, &ncols, &pcol, &pval));
490:   }

492:   PetscCall(PetscMalloc2(cmax, &pnval, cmax, &pncol));

494:   PetscCall(MatGetType(*P, &mtype));
495:   PetscCall(MatCreate(PetscObjectComm((PetscObject)*P), &Pnew));
496:   PetscCall(MatSetType(Pnew, mtype));
497:   PetscCall(MatSetSizes(Pnew, pn, pcn, PETSC_DETERMINE, PETSC_DETERMINE));
498:   PetscCall(MatSeqAIJSetPreallocation(Pnew, 0, lsparse));
499:   PetscCall(MatMPIAIJSetPreallocation(Pnew, 0, lsparse, 0, gsparse));

501:   for (i = ps; i < pf; i++) {
502:     PetscCall(MatGetRow(*P, i, &ncols, &pcol, &pval));
503:     pmax_pos = 0.;
504:     pmax_neg = 0.;
505:     for (j = 0; j < ncols; j++) {
506:       if (PetscRealPart(pval[j]) > pmax_pos) {
507:         pmax_pos = PetscRealPart(pval[j]);
508:       } else if (PetscRealPart(pval[j]) < pmax_neg) {
509:         pmax_neg = PetscRealPart(pval[j]);
510:       }
511:     }
512:     pthresh_pos = 0.;
513:     pthresh_neg = 0.;
514:     ptot_pos    = 0.;
515:     ptot_neg    = 0.;
516:     for (j = 0; j < ncols; j++) {
517:       if (PetscRealPart(pval[j]) >= cls->interp_threshold * pmax_pos) {
518:         pthresh_pos += PetscRealPart(pval[j]);
519:       } else if (PetscRealPart(pval[j]) <= cls->interp_threshold * pmax_neg) {
520:         pthresh_neg += PetscRealPart(pval[j]);
521:       }
522:       if (PetscRealPart(pval[j]) > 0.) {
523:         ptot_pos += PetscRealPart(pval[j]);
524:       } else {
525:         ptot_neg += PetscRealPart(pval[j]);
526:       }
527:     }
528:     if (PetscAbsReal(pthresh_pos) > 0.) ptot_pos /= pthresh_pos;
529:     if (PetscAbsReal(pthresh_neg) > 0.) ptot_neg /= pthresh_neg;
530:     idx = 0;
531:     for (j = 0; j < ncols; j++) {
532:       if (PetscRealPart(pval[j]) >= pmax_pos * cls->interp_threshold) {
533:         pnval[idx] = ptot_pos * pval[j];
534:         pncol[idx] = pcol[j];
535:         idx++;
536:       } else if (PetscRealPart(pval[j]) <= pmax_neg * cls->interp_threshold) {
537:         pnval[idx] = ptot_neg * pval[j];
538:         pncol[idx] = pcol[j];
539:         idx++;
540:       }
541:     }
542:     PetscCall(MatRestoreRow(*P, i, &ncols, &pcol, &pval));
543:     PetscCall(MatSetValues(Pnew, 1, &i, idx, pncol, pnval, INSERT_VALUES));
544:   }

546:   PetscCall(MatAssemblyBegin(Pnew, MAT_FINAL_ASSEMBLY));
547:   PetscCall(MatAssemblyEnd(Pnew, MAT_FINAL_ASSEMBLY));
548:   PetscCall(MatDestroy(P));

550:   *P = Pnew;
551:   PetscCall(PetscFree2(lsparse, gsparse));
552:   PetscCall(PetscFree2(pnval, pncol));
553:   PetscFunctionReturn(PETSC_SUCCESS);
554: }

556: static PetscErrorCode PCGAMGProlongator_Classical_Standard(PC pc, Mat A, PetscCoarsenData *agg_lists, Mat *P)
557: {
558:   Mat                lA, *lAs;
559:   MatType            mtype;
560:   Vec                cv;
561:   PetscInt          *gcid, *lcid, *lsparse, *gsparse, *picol;
562:   PetscInt           fs, fe, cs, ce, nl, i, j, k, li, lni, ci, ncols, maxcols, fn, cn, cid;
563:   PetscMPIInt        size;
564:   const PetscInt    *lidx, *icol, *gidx;
565:   PetscBool          iscoarse;
566:   PetscScalar        vi, pentry, pjentry;
567:   PetscScalar       *pcontrib, *pvcol;
568:   const PetscScalar *vcol;
569:   PetscReal          diag, jdiag, jwttotal;
570:   PetscInt           pncols;
571:   PetscSF            sf;
572:   PetscLayout        clayout;
573:   IS                 lis;

575:   PetscFunctionBegin;
576:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
577:   PetscCall(MatGetOwnershipRange(A, &fs, &fe));
578:   fn = fe - fs;
579:   PetscCall(ISCreateStride(PETSC_COMM_SELF, fe - fs, fs, 1, &lis));
580:   if (size > 1) {
581:     PetscCall(MatGetLayouts(A, NULL, &clayout));
582:     /* increase the overlap by two to get neighbors of neighbors */
583:     PetscCall(MatIncreaseOverlap(A, 1, &lis, 2));
584:     PetscCall(ISSort(lis));
585:     /* get the local part of A */
586:     PetscCall(MatCreateSubMatrices(A, 1, &lis, &lis, MAT_INITIAL_MATRIX, &lAs));
587:     lA = lAs[0];
588:     /* build an SF out of it */
589:     PetscCall(ISGetLocalSize(lis, &nl));
590:     PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)A), &sf));
591:     PetscCall(ISGetIndices(lis, &lidx));
592:     PetscCall(PetscSFSetGraphLayout(sf, clayout, nl, NULL, PETSC_COPY_VALUES, lidx));
593:     PetscCall(ISRestoreIndices(lis, &lidx));
594:   } else {
595:     lA = A;
596:     nl = fn;
597:   }
598:   /* create a communication structure for the overlapped portion and transmit coarse indices */
599:   PetscCall(PetscMalloc3(fn, &lsparse, fn, &gsparse, nl, &pcontrib));
600:   /* create coarse vector */
601:   cn = 0;
602:   for (i = 0; i < fn; i++) {
603:     PetscCall(PetscCDIsEmptyAt(agg_lists, i, &iscoarse));
604:     if (!iscoarse) cn++;
605:   }
606:   PetscCall(PetscMalloc1(fn, &gcid));
607:   PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)A), cn, PETSC_DECIDE, &cv));
608:   PetscCall(VecGetOwnershipRange(cv, &cs, &ce));
609:   cn = 0;
610:   for (i = 0; i < fn; i++) {
611:     PetscCall(PetscCDIsEmptyAt(agg_lists, i, &iscoarse));
612:     if (!iscoarse) {
613:       gcid[i] = cs + cn;
614:       cn++;
615:     } else {
616:       gcid[i] = -1;
617:     }
618:   }
619:   if (size > 1) {
620:     PetscCall(PetscMalloc1(nl, &lcid));
621:     PetscCall(PetscSFBcastBegin(sf, MPIU_INT, gcid, lcid, MPI_REPLACE));
622:     PetscCall(PetscSFBcastEnd(sf, MPIU_INT, gcid, lcid, MPI_REPLACE));
623:   } else {
624:     lcid = gcid;
625:   }
626:   /* count to preallocate the prolongator */
627:   PetscCall(ISGetIndices(lis, &gidx));
628:   maxcols = 0;
629:   /* count the number of unique contributing coarse cells for each fine */
630:   for (i = 0; i < nl; i++) {
631:     pcontrib[i] = 0.;
632:     PetscCall(MatGetRow(lA, i, &ncols, &icol, NULL));
633:     if (gidx[i] >= fs && gidx[i] < fe) {
634:       li          = gidx[i] - fs;
635:       lsparse[li] = 0;
636:       gsparse[li] = 0;
637:       cid         = lcid[i];
638:       if (cid >= 0) {
639:         lsparse[li] = 1;
640:       } else {
641:         for (j = 0; j < ncols; j++) {
642:           if (lcid[icol[j]] >= 0) {
643:             pcontrib[icol[j]] = 1.;
644:           } else {
645:             ci = icol[j];
646:             PetscCall(MatRestoreRow(lA, i, &ncols, &icol, NULL));
647:             PetscCall(MatGetRow(lA, ci, &ncols, &icol, NULL));
648:             for (k = 0; k < ncols; k++) {
649:               if (lcid[icol[k]] >= 0) pcontrib[icol[k]] = 1.;
650:             }
651:             PetscCall(MatRestoreRow(lA, ci, &ncols, &icol, NULL));
652:             PetscCall(MatGetRow(lA, i, &ncols, &icol, NULL));
653:           }
654:         }
655:         for (j = 0; j < ncols; j++) {
656:           if (lcid[icol[j]] >= 0 && pcontrib[icol[j]] != 0.) {
657:             lni = lcid[icol[j]];
658:             if (lni >= cs && lni < ce) {
659:               lsparse[li]++;
660:             } else {
661:               gsparse[li]++;
662:             }
663:             pcontrib[icol[j]] = 0.;
664:           } else {
665:             ci = icol[j];
666:             PetscCall(MatRestoreRow(lA, i, &ncols, &icol, NULL));
667:             PetscCall(MatGetRow(lA, ci, &ncols, &icol, NULL));
668:             for (k = 0; k < ncols; k++) {
669:               if (lcid[icol[k]] >= 0 && pcontrib[icol[k]] != 0.) {
670:                 lni = lcid[icol[k]];
671:                 if (lni >= cs && lni < ce) {
672:                   lsparse[li]++;
673:                 } else {
674:                   gsparse[li]++;
675:                 }
676:                 pcontrib[icol[k]] = 0.;
677:               }
678:             }
679:             PetscCall(MatRestoreRow(lA, ci, &ncols, &icol, NULL));
680:             PetscCall(MatGetRow(lA, i, &ncols, &icol, NULL));
681:           }
682:         }
683:       }
684:       if (lsparse[li] + gsparse[li] > maxcols) maxcols = lsparse[li] + gsparse[li];
685:     }
686:     PetscCall(MatRestoreRow(lA, i, &ncols, &icol, &vcol));
687:   }
688:   PetscCall(PetscMalloc2(maxcols, &picol, maxcols, &pvcol));
689:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), P));
690:   PetscCall(MatGetType(A, &mtype));
691:   PetscCall(MatSetType(*P, mtype));
692:   PetscCall(MatSetSizes(*P, fn, cn, PETSC_DETERMINE, PETSC_DETERMINE));
693:   PetscCall(MatMPIAIJSetPreallocation(*P, 0, lsparse, 0, gsparse));
694:   PetscCall(MatSeqAIJSetPreallocation(*P, 0, lsparse));
695:   for (i = 0; i < nl; i++) {
696:     diag = 0.;
697:     if (gidx[i] >= fs && gidx[i] < fe) {
698:       pncols = 0;
699:       cid    = lcid[i];
700:       if (cid >= 0) {
701:         pncols   = 1;
702:         picol[0] = cid;
703:         pvcol[0] = 1.;
704:       } else {
705:         PetscCall(MatGetRow(lA, i, &ncols, &icol, &vcol));
706:         for (j = 0; j < ncols; j++) {
707:           pentry = vcol[j];
708:           if (lcid[icol[j]] >= 0) {
709:             /* coarse neighbor */
710:             pcontrib[icol[j]] += pentry;
711:           } else if (icol[j] != i) {
712:             /* the neighbor is a strongly connected fine node */
713:             ci = icol[j];
714:             vi = vcol[j];
715:             PetscCall(MatRestoreRow(lA, i, &ncols, &icol, &vcol));
716:             PetscCall(MatGetRow(lA, ci, &ncols, &icol, &vcol));
717:             jwttotal = 0.;
718:             jdiag    = 0.;
719:             for (k = 0; k < ncols; k++) {
720:               if (ci == icol[k]) jdiag = PetscRealPart(vcol[k]);
721:             }
722:             for (k = 0; k < ncols; k++) {
723:               if (lcid[icol[k]] >= 0 && jdiag * PetscRealPart(vcol[k]) < 0.) {
724:                 pjentry = vcol[k];
725:                 jwttotal += PetscRealPart(pjentry);
726:               }
727:             }
728:             if (jwttotal != 0.) {
729:               jwttotal = PetscRealPart(vi) / jwttotal;
730:               for (k = 0; k < ncols; k++) {
731:                 if (lcid[icol[k]] >= 0 && jdiag * PetscRealPart(vcol[k]) < 0.) {
732:                   pjentry = vcol[k] * jwttotal;
733:                   pcontrib[icol[k]] += pjentry;
734:                 }
735:               }
736:             } else {
737:               diag += PetscRealPart(vi);
738:             }
739:             PetscCall(MatRestoreRow(lA, ci, &ncols, &icol, &vcol));
740:             PetscCall(MatGetRow(lA, i, &ncols, &icol, &vcol));
741:           } else {
742:             diag += PetscRealPart(vcol[j]);
743:           }
744:         }
745:         if (diag != 0.) {
746:           diag = 1. / diag;
747:           for (j = 0; j < ncols; j++) {
748:             if (lcid[icol[j]] >= 0 && pcontrib[icol[j]] != 0.) {
749:               /* the neighbor is a coarse node */
750:               if (PetscAbsScalar(pcontrib[icol[j]]) > 0.0) {
751:                 lni           = lcid[icol[j]];
752:                 pvcol[pncols] = -pcontrib[icol[j]] * diag;
753:                 picol[pncols] = lni;
754:                 pncols++;
755:               }
756:               pcontrib[icol[j]] = 0.;
757:             } else {
758:               /* the neighbor is a strongly connected fine node */
759:               ci = icol[j];
760:               PetscCall(MatRestoreRow(lA, i, &ncols, &icol, &vcol));
761:               PetscCall(MatGetRow(lA, ci, &ncols, &icol, &vcol));
762:               for (k = 0; k < ncols; k++) {
763:                 if (lcid[icol[k]] >= 0 && pcontrib[icol[k]] != 0.) {
764:                   if (PetscAbsScalar(pcontrib[icol[k]]) > 0.0) {
765:                     lni           = lcid[icol[k]];
766:                     pvcol[pncols] = -pcontrib[icol[k]] * diag;
767:                     picol[pncols] = lni;
768:                     pncols++;
769:                   }
770:                   pcontrib[icol[k]] = 0.;
771:                 }
772:               }
773:               PetscCall(MatRestoreRow(lA, ci, &ncols, &icol, &vcol));
774:               PetscCall(MatGetRow(lA, i, &ncols, &icol, &vcol));
775:             }
776:             pcontrib[icol[j]] = 0.;
777:           }
778:           PetscCall(MatRestoreRow(lA, i, &ncols, &icol, &vcol));
779:         }
780:       }
781:       ci = gidx[i];
782:       if (pncols > 0) PetscCall(MatSetValues(*P, 1, &ci, pncols, picol, pvcol, INSERT_VALUES));
783:     }
784:   }
785:   PetscCall(ISRestoreIndices(lis, &gidx));
786:   PetscCall(PetscFree2(picol, pvcol));
787:   PetscCall(PetscFree3(lsparse, gsparse, pcontrib));
788:   PetscCall(ISDestroy(&lis));
789:   PetscCall(PetscFree(gcid));
790:   if (size > 1) {
791:     PetscCall(PetscFree(lcid));
792:     PetscCall(MatDestroyMatrices(1, &lAs));
793:     PetscCall(PetscSFDestroy(&sf));
794:   }
795:   PetscCall(VecDestroy(&cv));
796:   PetscCall(MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY));
797:   PetscCall(MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY));
798:   PetscFunctionReturn(PETSC_SUCCESS);
799: }

801: static PetscErrorCode PCGAMGOptProlongator_Classical_Jacobi(PC pc, Mat A, Mat *P)
802: {
803:   PetscInt           f, s, n, cf, cs, i, idx;
804:   PetscInt          *coarserows;
805:   PetscInt           ncols;
806:   const PetscInt    *pcols;
807:   const PetscScalar *pvals;
808:   Mat                Pnew;
809:   Vec                diag;
810:   PC_MG             *mg      = (PC_MG *)pc->data;
811:   PC_GAMG           *pc_gamg = (PC_GAMG *)mg->innerctx;
812:   PC_GAMG_Classical *cls     = (PC_GAMG_Classical *)pc_gamg->subctx;

814:   PetscFunctionBegin;
815:   if (cls->nsmooths == 0) {
816:     PetscCall(PCGAMGTruncateProlongator_Private(pc, P));
817:     PetscFunctionReturn(PETSC_SUCCESS);
818:   }
819:   PetscCall(MatGetOwnershipRange(*P, &s, &f));
820:   n = f - s;
821:   PetscCall(MatGetOwnershipRangeColumn(*P, &cs, &cf));
822:   PetscCall(PetscMalloc1(n, &coarserows));
823:   /* identify the rows corresponding to coarse unknowns */
824:   idx = 0;
825:   for (i = s; i < f; i++) {
826:     PetscCall(MatGetRow(*P, i, &ncols, &pcols, &pvals));
827:     /* assume, for now, that it's a coarse unknown if it has a single unit entry */
828:     if (ncols == 1) {
829:       if (pvals[0] == 1.) {
830:         coarserows[idx] = i;
831:         idx++;
832:       }
833:     }
834:     PetscCall(MatRestoreRow(*P, i, &ncols, &pcols, &pvals));
835:   }
836:   PetscCall(MatCreateVecs(A, &diag, NULL));
837:   PetscCall(MatGetDiagonal(A, diag));
838:   PetscCall(VecReciprocal(diag));
839:   for (i = 0; i < cls->nsmooths; i++) {
840:     PetscCall(MatMatMult(A, *P, MAT_INITIAL_MATRIX, PETSC_CURRENT, &Pnew));
841:     PetscCall(MatZeroRows(Pnew, idx, coarserows, 0., NULL, NULL));
842:     PetscCall(MatDiagonalScale(Pnew, diag, NULL));
843:     PetscCall(MatAYPX(Pnew, -1.0, *P, DIFFERENT_NONZERO_PATTERN));
844:     PetscCall(MatDestroy(P));
845:     *P   = Pnew;
846:     Pnew = NULL;
847:   }
848:   PetscCall(VecDestroy(&diag));
849:   PetscCall(PetscFree(coarserows));
850:   PetscCall(PCGAMGTruncateProlongator_Private(pc, P));
851:   PetscFunctionReturn(PETSC_SUCCESS);
852: }

854: static PetscErrorCode PCGAMGProlongator_Classical(PC pc, Mat A, PetscCoarsenData *agg_lists, Mat *P)
855: {
856:   PetscErrorCode (*f)(PC, Mat, PetscCoarsenData *, Mat *);
857:   PC_MG             *mg      = (PC_MG *)pc->data;
858:   PC_GAMG           *pc_gamg = (PC_GAMG *)mg->innerctx;
859:   PC_GAMG_Classical *cls     = (PC_GAMG_Classical *)pc_gamg->subctx;

861:   PetscFunctionBegin;
862:   PetscCall(PetscFunctionListFind(PCGAMGClassicalProlongatorList, cls->prolongtype, &f));
863:   PetscCheck(f, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Cannot find PCGAMG Classical prolongator type");
864:   PetscCall((*f)(pc, A, agg_lists, P));
865:   PetscFunctionReturn(PETSC_SUCCESS);
866: }

868: static PetscErrorCode PCGAMGDestroy_Classical(PC pc)
869: {
870:   PC_MG   *mg      = (PC_MG *)pc->data;
871:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;

873:   PetscFunctionBegin;
874:   PetscCall(PetscFree(pc_gamg->subctx));
875:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGClassicalSetType_C", NULL));
876:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGClassicalGetType_C", NULL));
877:   PetscFunctionReturn(PETSC_SUCCESS);
878: }

880: static PetscErrorCode PCGAMGSetFromOptions_Classical(PC pc, PetscOptionItems *PetscOptionsObject)
881: {
882:   PC_MG             *mg      = (PC_MG *)pc->data;
883:   PC_GAMG           *pc_gamg = (PC_GAMG *)mg->innerctx;
884:   PC_GAMG_Classical *cls     = (PC_GAMG_Classical *)pc_gamg->subctx;
885:   char               tname[256];
886:   PetscBool          flg;

888:   PetscFunctionBegin;
889:   PetscOptionsHeadBegin(PetscOptionsObject, "GAMG-Classical options");
890:   PetscCall(PetscOptionsFList("-pc_gamg_classical_type", "Type of Classical AMG prolongation", "PCGAMGClassicalSetType", PCGAMGClassicalProlongatorList, cls->prolongtype, tname, sizeof(tname), &flg));
891:   if (flg) PetscCall(PCGAMGClassicalSetType(pc, tname));
892:   PetscCall(PetscOptionsReal("-pc_gamg_classical_interp_threshold", "Threshold for classical interpolator entries", "", cls->interp_threshold, &cls->interp_threshold, NULL));
893:   PetscCall(PetscOptionsInt("-pc_gamg_classical_nsmooths", "Threshold for classical interpolator entries", "", cls->nsmooths, &cls->nsmooths, NULL));
894:   PetscOptionsHeadEnd();
895:   PetscFunctionReturn(PETSC_SUCCESS);
896: }

898: static PetscErrorCode PCGAMGSetData_Classical(PC pc, Mat A)
899: {
900:   PC_MG   *mg      = (PC_MG *)pc->data;
901:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;

903:   PetscFunctionBegin;
904:   /* no data for classical AMG */
905:   pc_gamg->data           = NULL;
906:   pc_gamg->data_cell_cols = 0;
907:   pc_gamg->data_cell_rows = 0;
908:   pc_gamg->data_sz        = 0;
909:   PetscFunctionReturn(PETSC_SUCCESS);
910: }

912: static PetscErrorCode PCGAMGClassicalFinalizePackage(void)
913: {
914:   PetscFunctionBegin;
915:   PCGAMGClassicalPackageInitialized = PETSC_FALSE;
916:   PetscCall(PetscFunctionListDestroy(&PCGAMGClassicalProlongatorList));
917:   PetscFunctionReturn(PETSC_SUCCESS);
918: }

920: static PetscErrorCode PCGAMGClassicalInitializePackage(void)
921: {
922:   PetscFunctionBegin;
923:   if (PCGAMGClassicalPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
924:   PetscCall(PetscFunctionListAdd(&PCGAMGClassicalProlongatorList, PCGAMGCLASSICALDIRECT, PCGAMGProlongator_Classical_Direct));
925:   PetscCall(PetscFunctionListAdd(&PCGAMGClassicalProlongatorList, PCGAMGCLASSICALSTANDARD, PCGAMGProlongator_Classical_Standard));
926:   PetscCall(PetscRegisterFinalize(PCGAMGClassicalFinalizePackage));
927:   PetscFunctionReturn(PETSC_SUCCESS);
928: }

930: PetscErrorCode PCCreateGAMG_Classical(PC pc)
931: {
932:   PC_MG             *mg      = (PC_MG *)pc->data;
933:   PC_GAMG           *pc_gamg = (PC_GAMG *)mg->innerctx;
934:   PC_GAMG_Classical *pc_gamg_classical;

936:   PetscFunctionBegin;
937:   PetscCall(PCGAMGClassicalInitializePackage());
938:   if (pc_gamg->subctx) {
939:     /* call base class */
940:     PetscCall(PCDestroy_GAMG(pc));
941:   }

943:   /* create sub context for SA */
944:   PetscCall(PetscNew(&pc_gamg_classical));
945:   pc_gamg->subctx         = pc_gamg_classical;
946:   pc->ops->setfromoptions = PCGAMGSetFromOptions_Classical;
947:   /* reset does not do anything; setup not virtual */

949:   /* set internal function pointers */
950:   pc_gamg->ops->destroy        = PCGAMGDestroy_Classical;
951:   pc_gamg->ops->creategraph    = PCGAMGCreateGraph_Classical;
952:   pc_gamg->ops->coarsen        = PCGAMGCoarsen_Classical;
953:   pc_gamg->ops->prolongator    = PCGAMGProlongator_Classical;
954:   pc_gamg->ops->optprolongator = PCGAMGOptProlongator_Classical_Jacobi;
955:   pc_gamg->ops->setfromoptions = PCGAMGSetFromOptions_Classical;

957:   pc_gamg->ops->createdefaultdata     = PCGAMGSetData_Classical;
958:   pc_gamg_classical->interp_threshold = 0.2;
959:   pc_gamg_classical->nsmooths         = 0;
960:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGClassicalSetType_C", PCGAMGClassicalSetType_GAMG));
961:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGClassicalGetType_C", PCGAMGClassicalGetType_GAMG));
962:   PetscCall(PCGAMGClassicalSetType(pc, PCGAMGCLASSICALSTANDARD));
963:   PetscFunctionReturn(PETSC_SUCCESS);
964: }