Actual source code: ml.c


  2: /*
  3:    Provides an interface to the ML smoothed Aggregation
  4:    Note: Something non-obvious breaks -pc_mg_type ADDITIVE for parallel runs
  5:                                     Jed Brown, see [PETSC #18321, #18449].
  6: */
  7: #include <petsc/private/pcimpl.h>
  8: #include <petsc/private/pcmgimpl.h>
  9: #include <../src/mat/impls/aij/seq/aij.h>
 10: #include <../src/mat/impls/aij/mpi/mpiaij.h>
 11: #include <petscdm.h>

 13: EXTERN_C_BEGIN
 14: /* HAVE_CONFIG_H flag is required by ML include files */
 15: #ifndef HAVE_CONFIG_H
 16:   #define HAVE_CONFIG_H
 17: #endif
 18: #include <ml_include.h>
 19: #include <ml_viz_stats.h>
 20: EXTERN_C_END

 22: typedef enum {
 23:   PCML_NULLSPACE_AUTO,
 24:   PCML_NULLSPACE_USER,
 25:   PCML_NULLSPACE_BLOCK,
 26:   PCML_NULLSPACE_SCALAR
 27: } PCMLNullSpaceType;
 28: static const char *const PCMLNullSpaceTypes[] = {"AUTO", "USER", "BLOCK", "SCALAR", "PCMLNullSpaceType", "PCML_NULLSPACE_", 0};

 30: /* The context (data structure) at each grid level */
 31: typedef struct {
 32:   Vec x, b, r; /* global vectors */
 33:   Mat A, P, R;
 34:   KSP ksp;
 35:   Vec coords; /* projected by ML, if PCSetCoordinates is called; values packed by node */
 36: } GridCtx;

 38: /* The context used to input PETSc matrix into ML at fine grid */
 39: typedef struct {
 40:   Mat          A;    /* Petsc matrix in aij format */
 41:   Mat          Aloc; /* local portion of A to be used by ML */
 42:   Vec          x, y;
 43:   ML_Operator *mlmat;
 44:   PetscScalar *pwork; /* tmp array used by PetscML_comm() */
 45: } FineGridCtx;

 47: /* The context associates a ML matrix with a PETSc shell matrix */
 48: typedef struct {
 49:   Mat          A;     /* PETSc shell matrix associated with mlmat */
 50:   ML_Operator *mlmat; /* ML matrix assorciated with A */
 51: } Mat_MLShell;

 53: /* Private context for the ML preconditioner */
 54: typedef struct {
 55:   ML               *ml_object;
 56:   ML_Aggregate     *agg_object;
 57:   GridCtx          *gridctx;
 58:   FineGridCtx      *PetscMLdata;
 59:   PetscInt          Nlevels, MaxNlevels, MaxCoarseSize, CoarsenScheme, EnergyMinimization, MinPerProc, PutOnSingleProc, RepartitionType, ZoltanScheme;
 60:   PetscReal         Threshold, DampingFactor, EnergyMinimizationDropTol, MaxMinRatio, AuxThreshold;
 61:   PetscBool         SpectralNormScheme_Anorm, BlockScaling, EnergyMinimizationCheap, Symmetrize, OldHierarchy, KeepAggInfo, Reusable, Repartition, Aux;
 62:   PetscBool         reuse_interpolation;
 63:   PCMLNullSpaceType nulltype;
 64:   PetscMPIInt       size; /* size of communicator for pc->pmat */
 65:   PetscInt          dim;  /* data from PCSetCoordinates(_ML) */
 66:   PetscInt          nloc;
 67:   PetscReal        *coords; /* ML has a grid object for each level: the finest grid will point into coords */
 68: } PC_ML;

 70: static int PetscML_getrow(ML_Operator *ML_data, int N_requested_rows, int requested_rows[], int allocated_space, int columns[], double values[], int row_lengths[])
 71: {
 72:   PetscInt     m, i, j, k = 0, row, *aj;
 73:   PetscScalar *aa;
 74:   FineGridCtx *ml = (FineGridCtx *)ML_Get_MyGetrowData(ML_data);
 75:   Mat_SeqAIJ  *a  = (Mat_SeqAIJ *)ml->Aloc->data;

 77:   if (MatGetSize(ml->Aloc, &m, NULL)) return (0);
 78:   for (i = 0; i < N_requested_rows; i++) {
 79:     row            = requested_rows[i];
 80:     row_lengths[i] = a->ilen[row];
 81:     if (allocated_space < k + row_lengths[i]) return (0);
 82:     if ((row >= 0) || (row <= (m - 1))) {
 83:       aj = a->j + a->i[row];
 84:       aa = a->a + a->i[row];
 85:       for (j = 0; j < row_lengths[i]; j++) {
 86:         columns[k]  = aj[j];
 87:         values[k++] = aa[j];
 88:       }
 89:     }
 90:   }
 91:   return (1);
 92: }

 94: static PetscErrorCode PetscML_comm(double p[], void *ML_data)
 95: {
 96:   FineGridCtx       *ml = (FineGridCtx *)ML_data;
 97:   Mat                A  = ml->A;
 98:   Mat_MPIAIJ        *a  = (Mat_MPIAIJ *)A->data;
 99:   PetscMPIInt        size;
100:   PetscInt           i, in_length = A->rmap->n, out_length = ml->Aloc->cmap->n;
101:   const PetscScalar *array;

103:   MPI_Comm_size(PetscObjectComm((PetscObject)A), &size);
104:   if (size == 1) return 0;

106:   VecPlaceArray(ml->y, p);
107:   VecScatterBegin(a->Mvctx, ml->y, a->lvec, INSERT_VALUES, SCATTER_FORWARD);
108:   VecScatterEnd(a->Mvctx, ml->y, a->lvec, INSERT_VALUES, SCATTER_FORWARD);
109:   VecResetArray(ml->y);
110:   VecGetArrayRead(a->lvec, &array);
111:   for (i = in_length; i < out_length; i++) p[i] = array[i - in_length];
112:   VecRestoreArrayRead(a->lvec, &array);
113:   return 0;
114: }

116: static int PetscML_matvec(ML_Operator *ML_data, int in_length, double p[], int out_length, double ap[])
117: {
118:   FineGridCtx *ml = (FineGridCtx *)ML_Get_MyMatvecData(ML_data);
119:   Mat          A = ml->A, Aloc = ml->Aloc;
120:   PetscMPIInt  size;
121:   PetscScalar *pwork = ml->pwork;
122:   PetscInt     i;

124:   MPI_Comm_size(PetscObjectComm((PetscObject)A), &size);
125:   if (size == 1) {
126:     VecPlaceArray(ml->x, p);
127:   } else {
128:     for (i = 0; i < in_length; i++) pwork[i] = p[i];
129:     PetscML_comm(pwork, ml);
130:     VecPlaceArray(ml->x, pwork);
131:   }
132:   VecPlaceArray(ml->y, ap);
133:   MatMult(Aloc, ml->x, ml->y);
134:   VecResetArray(ml->x);
135:   VecResetArray(ml->y);
136:   return 0;
137: }

139: static PetscErrorCode MatMult_ML(Mat A, Vec x, Vec y)
140: {
141:   Mat_MLShell       *shell;
142:   PetscScalar       *yarray;
143:   const PetscScalar *xarray;
144:   PetscInt           x_length, y_length;

146:   MatShellGetContext(A, &shell);
147:   VecGetArrayRead(x, &xarray);
148:   VecGetArray(y, &yarray);
149:   x_length = shell->mlmat->invec_leng;
150:   y_length = shell->mlmat->outvec_leng;
151:   PetscStackCallExternalVoid("ML_Operator_Apply", ML_Operator_Apply(shell->mlmat, x_length, (PetscScalar *)xarray, y_length, yarray));
152:   VecRestoreArrayRead(x, &xarray);
153:   VecRestoreArray(y, &yarray);
154:   return 0;
155: }

157: /* newtype is ignored since only handles one case */
158: static PetscErrorCode MatConvert_MPIAIJ_ML(Mat A, MatType newtype, MatReuse scall, Mat *Aloc)
159: {
160:   Mat_MPIAIJ  *mpimat = (Mat_MPIAIJ *)A->data;
161:   Mat_SeqAIJ  *mat, *a = (Mat_SeqAIJ *)(mpimat->A)->data, *b = (Mat_SeqAIJ *)(mpimat->B)->data;
162:   PetscInt    *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j;
163:   PetscScalar *aa, *ba, *ca;
164:   PetscInt     am = A->rmap->n, an = A->cmap->n, i, j, k;
165:   PetscInt    *ci, *cj, ncols;

168:   MatSeqAIJGetArrayRead(mpimat->A, (const PetscScalar **)&aa);
169:   MatSeqAIJGetArrayRead(mpimat->B, (const PetscScalar **)&ba);
170:   if (scall == MAT_INITIAL_MATRIX) {
171:     PetscMalloc1(1 + am, &ci);
172:     ci[0] = 0;
173:     for (i = 0; i < am; i++) ci[i + 1] = ci[i] + (ai[i + 1] - ai[i]) + (bi[i + 1] - bi[i]);
174:     PetscMalloc1(1 + ci[am], &cj);
175:     PetscMalloc1(1 + ci[am], &ca);

177:     k = 0;
178:     for (i = 0; i < am; i++) {
179:       /* diagonal portion of A */
180:       ncols = ai[i + 1] - ai[i];
181:       for (j = 0; j < ncols; j++) {
182:         cj[k]   = *aj++;
183:         ca[k++] = *aa++;
184:       }
185:       /* off-diagonal portion of A */
186:       ncols = bi[i + 1] - bi[i];
187:       for (j = 0; j < ncols; j++) {
188:         cj[k] = an + (*bj);
189:         bj++;
190:         ca[k++] = *ba++;
191:       }
192:     }

195:     /* put together the new matrix */
196:     an = mpimat->A->cmap->n + mpimat->B->cmap->n;
197:     MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, am, an, ci, cj, ca, Aloc);

199:     /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
200:     /* Since these are PETSc arrays, change flags to free them as necessary. */
201:     mat          = (Mat_SeqAIJ *)(*Aloc)->data;
202:     mat->free_a  = PETSC_TRUE;
203:     mat->free_ij = PETSC_TRUE;

205:     mat->nonew = 0;
206:   } else if (scall == MAT_REUSE_MATRIX) {
207:     mat = (Mat_SeqAIJ *)(*Aloc)->data;
208:     ci  = mat->i;
209:     cj  = mat->j;
210:     ca  = mat->a;
211:     for (i = 0; i < am; i++) {
212:       /* diagonal portion of A */
213:       ncols = ai[i + 1] - ai[i];
214:       for (j = 0; j < ncols; j++) *ca++ = *aa++;
215:       /* off-diagonal portion of A */
216:       ncols = bi[i + 1] - bi[i];
217:       for (j = 0; j < ncols; j++) *ca++ = *ba++;
218:     }
219:   } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid MatReuse %d", (int)scall);
220:   MatSeqAIJRestoreArrayRead(mpimat->A, (const PetscScalar **)&aa);
221:   MatSeqAIJRestoreArrayRead(mpimat->B, (const PetscScalar **)&ba);
222:   return 0;
223: }

225: static PetscErrorCode MatDestroy_ML(Mat A)
226: {
227:   Mat_MLShell *shell;

229:   MatShellGetContext(A, &shell);
230:   PetscFree(shell);
231:   return 0;
232: }

234: static PetscErrorCode MatWrapML_SeqAIJ(ML_Operator *mlmat, MatReuse reuse, Mat *newmat)
235: {
236:   struct ML_CSR_MSRdata *matdata = (struct ML_CSR_MSRdata *)mlmat->data;
237:   PetscInt               m = mlmat->outvec_leng, n = mlmat->invec_leng, *nnz = NULL, nz_max;
238:   PetscInt              *ml_cols = matdata->columns, *ml_rowptr = matdata->rowptr, *aj, i;
239:   PetscScalar           *ml_vals = matdata->values, *aa;

242:   if (m != n) { /* ML Pmat and Rmat are in CSR format. Pass array pointers into SeqAIJ matrix */
243:     if (reuse) {
244:       Mat_SeqAIJ *aij = (Mat_SeqAIJ *)(*newmat)->data;
245:       aij->i          = ml_rowptr;
246:       aij->j          = ml_cols;
247:       aij->a          = ml_vals;
248:     } else {
249:       /* sort ml_cols and ml_vals */
250:       PetscMalloc1(m + 1, &nnz);
251:       for (i = 0; i < m; i++) nnz[i] = ml_rowptr[i + 1] - ml_rowptr[i];
252:       aj = ml_cols;
253:       aa = ml_vals;
254:       for (i = 0; i < m; i++) {
255:         PetscSortIntWithScalarArray(nnz[i], aj, aa);
256:         aj += nnz[i];
257:         aa += nnz[i];
258:       }
259:       MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, m, n, ml_rowptr, ml_cols, ml_vals, newmat);
260:       PetscFree(nnz);
261:     }
262:     MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY);
263:     MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY);
264:     return 0;
265:   }

267:   nz_max = PetscMax(1, mlmat->max_nz_per_row);
268:   PetscMalloc2(nz_max, &aa, nz_max, &aj);
269:   if (!reuse) {
270:     MatCreate(PETSC_COMM_SELF, newmat);
271:     MatSetSizes(*newmat, m, n, PETSC_DECIDE, PETSC_DECIDE);
272:     MatSetType(*newmat, MATSEQAIJ);
273:     /* keep track of block size for A matrices */
274:     MatSetBlockSize(*newmat, mlmat->num_PDEs);

276:     PetscMalloc1(m, &nnz);
277:     for (i = 0; i < m; i++) PetscStackCallExternalVoid("ML_Operator_Getrow", ML_Operator_Getrow(mlmat, 1, &i, nz_max, aj, aa, &nnz[i]));
278:     MatSeqAIJSetPreallocation(*newmat, 0, nnz);
279:   }
280:   for (i = 0; i < m; i++) {
281:     PetscInt ncols;

283:     PetscStackCallExternalVoid("ML_Operator_Getrow", ML_Operator_Getrow(mlmat, 1, &i, nz_max, aj, aa, &ncols));
284:     MatSetValues(*newmat, 1, &i, ncols, aj, aa, INSERT_VALUES);
285:   }
286:   MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY);
287:   MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY);

289:   PetscFree2(aa, aj);
290:   PetscFree(nnz);
291:   return 0;
292: }

294: static PetscErrorCode MatWrapML_SHELL(ML_Operator *mlmat, MatReuse reuse, Mat *newmat)
295: {
296:   PetscInt     m, n;
297:   ML_Comm     *MLcomm;
298:   Mat_MLShell *shellctx;

300:   m = mlmat->outvec_leng;
301:   n = mlmat->invec_leng;

303:   if (reuse) {
304:     MatShellGetContext(*newmat, &shellctx);
305:     shellctx->mlmat = mlmat;
306:     return 0;
307:   }

309:   MLcomm = mlmat->comm;

311:   PetscNew(&shellctx);
312:   MatCreateShell(MLcomm->USR_comm, m, n, PETSC_DETERMINE, PETSC_DETERMINE, shellctx, newmat);
313:   MatShellSetOperation(*newmat, MATOP_MULT, (void (*)(void))MatMult_ML);
314:   MatShellSetOperation(*newmat, MATOP_DESTROY, (void (*)(void))MatDestroy_ML);

316:   shellctx->A     = *newmat;
317:   shellctx->mlmat = mlmat;
318:   return 0;
319: }

321: static PetscErrorCode MatWrapML_MPIAIJ(ML_Operator *mlmat, MatReuse reuse, Mat *newmat)
322: {
323:   PetscInt    *aj;
324:   PetscScalar *aa;
325:   PetscInt     i, j, *gordering;
326:   PetscInt     m = mlmat->outvec_leng, n, nz_max, row;
327:   Mat          A;

330:   n = mlmat->invec_leng;

333:   /* create global row numbering for a ML_Operator */
334:   PetscStackCallExternalVoid("ML_build_global_numbering", ML_build_global_numbering(mlmat, &gordering, "rows"));

336:   nz_max = PetscMax(1, mlmat->max_nz_per_row) + 1;
337:   PetscMalloc2(nz_max, &aa, nz_max, &aj);
338:   if (reuse) {
339:     A = *newmat;
340:   } else {
341:     PetscInt *nnzA, *nnzB, *nnz;
342:     PetscInt  rstart;
343:     MatCreate(mlmat->comm->USR_comm, &A);
344:     MatSetSizes(A, m, n, PETSC_DECIDE, PETSC_DECIDE);
345:     MatSetType(A, MATMPIAIJ);
346:     /* keep track of block size for A matrices */
347:     MatSetBlockSize(A, mlmat->num_PDEs);
348:     PetscMalloc3(m, &nnzA, m, &nnzB, m, &nnz);
349:     MPI_Scan(&m, &rstart, 1, MPIU_INT, MPI_SUM, mlmat->comm->USR_comm);
350:     rstart -= m;

352:     for (i = 0; i < m; i++) {
353:       row = gordering[i] - rstart;
354:       PetscStackCallExternalVoid("ML_Operator_Getrow", ML_Operator_Getrow(mlmat, 1, &i, nz_max, aj, aa, &nnz[i]));
355:       nnzA[row] = 0;
356:       for (j = 0; j < nnz[i]; j++) {
357:         if (aj[j] < m) nnzA[row]++;
358:       }
359:       nnzB[row] = nnz[i] - nnzA[row];
360:     }
361:     MatMPIAIJSetPreallocation(A, 0, nnzA, 0, nnzB);
362:     PetscFree3(nnzA, nnzB, nnz);
363:   }
364:   for (i = 0; i < m; i++) {
365:     PetscInt ncols;
366:     row = gordering[i];

368:     PetscStackCallExternalVoid(",ML_Operator_Getrow", ML_Operator_Getrow(mlmat, 1, &i, nz_max, aj, aa, &ncols));
369:     for (j = 0; j < ncols; j++) aj[j] = gordering[aj[j]];
370:     MatSetValues(A, 1, &row, ncols, aj, aa, INSERT_VALUES);
371:   }
372:   PetscStackCallExternalVoid("ML_free", ML_free(gordering));
373:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
374:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
375:   *newmat = A;

377:   PetscFree2(aa, aj);
378:   return 0;
379: }

381: /* -------------------------------------------------------------------------- */
382: /*
383:    PCSetCoordinates_ML

385:    Input Parameter:
386:    .  pc - the preconditioner context
387: */
388: static PetscErrorCode PCSetCoordinates_ML(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords)
389: {
390:   PC_MG   *mg    = (PC_MG *)pc->data;
391:   PC_ML   *pc_ml = (PC_ML *)mg->innerctx;
392:   PetscInt arrsz, oldarrsz, bs, my0, kk, ii, nloc, Iend, aloc;
393:   Mat      Amat = pc->pmat;

395:   /* this function copied and modified from PCSetCoordinates_GEO -TGI */
397:   MatGetBlockSize(Amat, &bs);

399:   MatGetOwnershipRange(Amat, &my0, &Iend);
400:   aloc = (Iend - my0);
401:   nloc = (Iend - my0) / bs;


405:   oldarrsz    = pc_ml->dim * pc_ml->nloc;
406:   pc_ml->dim  = ndm;
407:   pc_ml->nloc = nloc;
408:   arrsz       = ndm * nloc;

410:   /* create data - syntactic sugar that should be refactored at some point */
411:   if (pc_ml->coords == 0 || (oldarrsz != arrsz)) {
412:     PetscFree(pc_ml->coords);
413:     PetscMalloc1(arrsz, &pc_ml->coords);
414:   }
415:   for (kk = 0; kk < arrsz; kk++) pc_ml->coords[kk] = -999.;
416:   /* copy data in - column oriented */
417:   if (nloc == a_nloc) {
418:     for (kk = 0; kk < nloc; kk++) {
419:       for (ii = 0; ii < ndm; ii++) pc_ml->coords[ii * nloc + kk] = coords[kk * ndm + ii];
420:     }
421:   } else { /* assumes the coordinates are blocked */
422:     for (kk = 0; kk < nloc; kk++) {
423:       for (ii = 0; ii < ndm; ii++) pc_ml->coords[ii * nloc + kk] = coords[bs * kk * ndm + ii];
424:     }
425:   }
426:   return 0;
427: }

429: /* -----------------------------------------------------------------------------*/
430: extern PetscErrorCode PCReset_MG(PC);
431: PetscErrorCode        PCReset_ML(PC pc)
432: {
433:   PC_MG   *mg    = (PC_MG *)pc->data;
434:   PC_ML   *pc_ml = (PC_ML *)mg->innerctx;
435:   PetscInt level, fine_level = pc_ml->Nlevels - 1, dim = pc_ml->dim;

437:   if (dim) {
438:     for (level = 0; level <= fine_level; level++) VecDestroy(&pc_ml->gridctx[level].coords);
439:     if (pc_ml->ml_object && pc_ml->ml_object->Grid) {
440:       ML_Aggregate_Viz_Stats *grid_info = (ML_Aggregate_Viz_Stats *)pc_ml->ml_object->Grid[0].Grid;
441:       grid_info->x                      = 0; /* do this so ML doesn't try to free coordinates */
442:       grid_info->y                      = 0;
443:       grid_info->z                      = 0;
444:       PetscStackCallExternalVoid("ML_Operator_Getrow", ML_Aggregate_VizAndStats_Clean(pc_ml->ml_object));
445:     }
446:   }
447:   PetscStackCallExternalVoid("ML_Aggregate_Destroy", ML_Aggregate_Destroy(&pc_ml->agg_object));
448:   PetscStackCallExternalVoid("ML_Aggregate_Destroy", ML_Destroy(&pc_ml->ml_object));

450:   if (pc_ml->PetscMLdata) {
451:     PetscFree(pc_ml->PetscMLdata->pwork);
452:     MatDestroy(&pc_ml->PetscMLdata->Aloc);
453:     VecDestroy(&pc_ml->PetscMLdata->x);
454:     VecDestroy(&pc_ml->PetscMLdata->y);
455:   }
456:   PetscFree(pc_ml->PetscMLdata);

458:   if (pc_ml->gridctx) {
459:     for (level = 0; level < fine_level; level++) {
460:       if (pc_ml->gridctx[level].A) MatDestroy(&pc_ml->gridctx[level].A);
461:       if (pc_ml->gridctx[level].P) MatDestroy(&pc_ml->gridctx[level].P);
462:       if (pc_ml->gridctx[level].R) MatDestroy(&pc_ml->gridctx[level].R);
463:       if (pc_ml->gridctx[level].x) VecDestroy(&pc_ml->gridctx[level].x);
464:       if (pc_ml->gridctx[level].b) VecDestroy(&pc_ml->gridctx[level].b);
465:       if (pc_ml->gridctx[level + 1].r) VecDestroy(&pc_ml->gridctx[level + 1].r);
466:     }
467:   }
468:   PetscFree(pc_ml->gridctx);
469:   PetscFree(pc_ml->coords);

471:   pc_ml->dim  = 0;
472:   pc_ml->nloc = 0;
473:   PCReset_MG(pc);
474:   return 0;
475: }
476: /* -------------------------------------------------------------------------- */
477: /*
478:    PCSetUp_ML - Prepares for the use of the ML preconditioner
479:                     by setting data structures and options.

481:    Input Parameter:
482: .  pc - the preconditioner context

484:    Application Interface Routine: PCSetUp()

486:    Note:
487:    The interface routine PCSetUp() is not usually called directly by
488:    the user, but instead is called by PCApply() if necessary.
489: */
490: extern PetscErrorCode PCSetFromOptions_MG(PC, PetscOptionItems *PetscOptionsObject);
491: extern PetscErrorCode PCReset_MG(PC);

493: PetscErrorCode PCSetUp_ML(PC pc)
494: {
495:   PetscMPIInt      size;
496:   FineGridCtx     *PetscMLdata;
497:   ML              *ml_object;
498:   ML_Aggregate    *agg_object;
499:   ML_Operator     *mlmat;
500:   PetscInt         nlocal_allcols, Nlevels, mllevel, level, level1, m, fine_level, bs;
501:   Mat              A, Aloc;
502:   GridCtx         *gridctx;
503:   PC_MG           *mg    = (PC_MG *)pc->data;
504:   PC_ML           *pc_ml = (PC_ML *)mg->innerctx;
505:   PetscBool        isSeq, isMPI;
506:   KSP              smoother;
507:   PC               subpc;
508:   PetscInt         mesh_level, old_mesh_level;
509:   MatInfo          info;
510:   static PetscBool cite = PETSC_FALSE;

512:   PetscCall(PetscCitationsRegister("@TechReport{ml_users_guide,\n  author = {M. Sala and J.J. Hu and R.S. Tuminaro},\n  title = {{ML}3.1 {S}moothed {A}ggregation {U}ser's {G}uide},\n  institution =  {Sandia National Laboratories},\n  number = "
513:                                    "{SAND2004-4821},\n  year = 2004\n}\n",
514:                                    &cite));
515:   A = pc->pmat;
516:   MPI_Comm_size(PetscObjectComm((PetscObject)A), &size);

518:   if (pc->setupcalled) {
519:     if (pc->flag == SAME_NONZERO_PATTERN && pc_ml->reuse_interpolation) {
520:       /*
521:        Reuse interpolaton instead of recomputing aggregates and updating the whole hierarchy. This is less expensive for
522:        multiple solves in which the matrix is not changing too quickly.
523:        */
524:       ml_object             = pc_ml->ml_object;
525:       gridctx               = pc_ml->gridctx;
526:       Nlevels               = pc_ml->Nlevels;
527:       fine_level            = Nlevels - 1;
528:       gridctx[fine_level].A = A;

530:       PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &isSeq);
531:       PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &isMPI);
532:       if (isMPI) {
533:         MatConvert_MPIAIJ_ML(A, NULL, MAT_INITIAL_MATRIX, &Aloc);
534:       } else if (isSeq) {
535:         Aloc = A;
536:         PetscObjectReference((PetscObject)Aloc);
537:       } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Matrix type '%s' cannot be used with ML. ML can only handle AIJ matrices.", ((PetscObject)A)->type_name);

539:       MatGetSize(Aloc, &m, &nlocal_allcols);
540:       PetscMLdata = pc_ml->PetscMLdata;
541:       MatDestroy(&PetscMLdata->Aloc);
542:       PetscMLdata->A    = A;
543:       PetscMLdata->Aloc = Aloc;
544:       PetscStackCallExternalVoid("ML_Aggregate_Destroy", ML_Init_Amatrix(ml_object, 0, m, m, PetscMLdata));
545:       PetscStackCallExternalVoid("ML_Set_Amatrix_Matvec", ML_Set_Amatrix_Matvec(ml_object, 0, PetscML_matvec));

547:       mesh_level = ml_object->ML_finest_level;
548:       while (ml_object->SingleLevel[mesh_level].Rmat->to) {
549:         old_mesh_level = mesh_level;
550:         mesh_level     = ml_object->SingleLevel[mesh_level].Rmat->to->levelnum;

552:         /* clean and regenerate A */
553:         mlmat = &(ml_object->Amat[mesh_level]);
554:         PetscStackCallExternalVoid("ML_Operator_Clean", ML_Operator_Clean(mlmat));
555:         PetscStackCallExternalVoid("ML_Operator_Init", ML_Operator_Init(mlmat, ml_object->comm));
556:         PetscStackCallExternalVoid("ML_Gen_AmatrixRAP", ML_Gen_AmatrixRAP(ml_object, old_mesh_level, mesh_level));
557:       }

559:       level = fine_level - 1;
560:       if (size == 1) { /* convert ML P, R and A into seqaij format */
561:         for (mllevel = 1; mllevel < Nlevels; mllevel++) {
562:           mlmat = &(ml_object->Amat[mllevel]);
563:           MatWrapML_SeqAIJ(mlmat, MAT_REUSE_MATRIX, &gridctx[level].A);
564:           level--;
565:         }
566:       } else { /* convert ML P and R into shell format, ML A into mpiaij format */
567:         for (mllevel = 1; mllevel < Nlevels; mllevel++) {
568:           mlmat = &(ml_object->Amat[mllevel]);
569:           MatWrapML_MPIAIJ(mlmat, MAT_REUSE_MATRIX, &gridctx[level].A);
570:           level--;
571:         }
572:       }

574:       for (level = 0; level < fine_level; level++) {
575:         if (level > 0) PCMGSetResidual(pc, level, PCMGResidualDefault, gridctx[level].A);
576:         KSPSetOperators(gridctx[level].ksp, gridctx[level].A, gridctx[level].A);
577:       }
578:       PCMGSetResidual(pc, fine_level, PCMGResidualDefault, gridctx[fine_level].A);
579:       KSPSetOperators(gridctx[fine_level].ksp, gridctx[level].A, gridctx[fine_level].A);

581:       PCSetUp_MG(pc);
582:       return 0;
583:     } else {
584:       /* since ML can change the size of vectors/matrices at any level we must destroy everything */
585:       PCReset_ML(pc);
586:     }
587:   }

589:   /* setup special features of PCML */
590:   /*--------------------------------*/
591:   /* convert A to Aloc to be used by ML at fine grid */
592:   pc_ml->size = size;
593:   PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &isSeq);
594:   PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &isMPI);
595:   if (isMPI) {
596:     MatConvert_MPIAIJ_ML(A, NULL, MAT_INITIAL_MATRIX, &Aloc);
597:   } else if (isSeq) {
598:     Aloc = A;
599:     PetscObjectReference((PetscObject)Aloc);
600:   } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Matrix type '%s' cannot be used with ML. ML can only handle AIJ matrices.", ((PetscObject)A)->type_name);

602:   /* create and initialize struct 'PetscMLdata' */
603:   PetscNew(&PetscMLdata);
604:   pc_ml->PetscMLdata = PetscMLdata;
605:   PetscMalloc1(Aloc->cmap->n + 1, &PetscMLdata->pwork);

607:   MatCreateVecs(Aloc, &PetscMLdata->x, &PetscMLdata->y);

609:   PetscMLdata->A    = A;
610:   PetscMLdata->Aloc = Aloc;
611:   if (pc_ml->dim) { /* create vecs around the coordinate data given */
612:     PetscInt   i, j, dim = pc_ml->dim;
613:     PetscInt   nloc = pc_ml->nloc, nlocghost;
614:     PetscReal *ghostedcoords;

616:     MatGetBlockSize(A, &bs);
617:     nlocghost = Aloc->cmap->n / bs;
618:     PetscMalloc1(dim * nlocghost, &ghostedcoords);
619:     for (i = 0; i < dim; i++) {
620:       /* copy coordinate values into first component of pwork */
621:       for (j = 0; j < nloc; j++) PetscMLdata->pwork[bs * j] = pc_ml->coords[nloc * i + j];
622:       /* get the ghost values */
623:       PetscML_comm(PetscMLdata->pwork, PetscMLdata);
624:       /* write into the vector */
625:       for (j = 0; j < nlocghost; j++) ghostedcoords[i * nlocghost + j] = PetscMLdata->pwork[bs * j];
626:     }
627:     /* replace the original coords with the ghosted coords, because these are
628:      * what ML needs */
629:     PetscFree(pc_ml->coords);
630:     pc_ml->coords = ghostedcoords;
631:   }

633:   /* create ML discretization matrix at fine grid */
634:   /* ML requires input of fine-grid matrix. It determines nlevels. */
635:   MatGetSize(Aloc, &m, &nlocal_allcols);
636:   MatGetBlockSize(A, &bs);
637:   PetscStackCallExternalVoid("ML_Create", ML_Create(&ml_object, pc_ml->MaxNlevels));
638:   PetscStackCallExternalVoid("ML_Comm_Set_UsrComm", ML_Comm_Set_UsrComm(ml_object->comm, PetscObjectComm((PetscObject)A)));
639:   pc_ml->ml_object = ml_object;
640:   PetscStackCallExternalVoid("ML_Init_Amatrix", ML_Init_Amatrix(ml_object, 0, m, m, PetscMLdata));
641:   PetscStackCallExternalVoid("ML_Set_Amatrix_Getrow", ML_Set_Amatrix_Getrow(ml_object, 0, PetscML_getrow, PetscML_comm, nlocal_allcols));
642:   PetscStackCallExternalVoid("ML_Set_Amatrix_Matvec", ML_Set_Amatrix_Matvec(ml_object, 0, PetscML_matvec));

644:   PetscStackCallExternalVoid("ML_Set_Symmetrize", ML_Set_Symmetrize(ml_object, pc_ml->Symmetrize ? ML_YES : ML_NO));

646:   /* aggregation */
647:   PetscStackCallExternalVoid("ML_Aggregate_Create", ML_Aggregate_Create(&agg_object));
648:   pc_ml->agg_object = agg_object;

650:   {
651:     MatNullSpace mnull;
652:     MatGetNearNullSpace(A, &mnull);
653:     if (pc_ml->nulltype == PCML_NULLSPACE_AUTO) {
654:       if (mnull) pc_ml->nulltype = PCML_NULLSPACE_USER;
655:       else if (bs > 1) pc_ml->nulltype = PCML_NULLSPACE_BLOCK;
656:       else pc_ml->nulltype = PCML_NULLSPACE_SCALAR;
657:     }
658:     switch (pc_ml->nulltype) {
659:     case PCML_NULLSPACE_USER: {
660:       PetscScalar       *nullvec;
661:       const PetscScalar *v;
662:       PetscBool          has_const;
663:       PetscInt           i, j, mlocal, nvec, M;
664:       const Vec         *vecs;

667:       MatGetSize(A, &M, NULL);
668:       MatGetLocalSize(Aloc, &mlocal, NULL);
669:       MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs);
670:       PetscMalloc1((nvec + !!has_const) * mlocal, &nullvec);
671:       if (has_const)
672:         for (i = 0; i < mlocal; i++) nullvec[i] = 1.0 / M;
673:       for (i = 0; i < nvec; i++) {
674:         VecGetArrayRead(vecs[i], &v);
675:         for (j = 0; j < mlocal; j++) nullvec[(i + !!has_const) * mlocal + j] = v[j];
676:         VecRestoreArrayRead(vecs[i], &v);
677:       }
678:       PetscStackCallExternalVoid("ML_Aggregate_Create", ML_Aggregate_Set_NullSpace(agg_object, bs, nvec + !!has_const, nullvec, mlocal));
679:       PetscFree(nullvec);
680:     } break;
681:     case PCML_NULLSPACE_BLOCK:
682:       PetscStackCallExternalVoid("ML_Aggregate_Set_NullSpace", ML_Aggregate_Set_NullSpace(agg_object, bs, bs, 0, 0));
683:       break;
684:     case PCML_NULLSPACE_SCALAR:
685:       break;
686:     default:
687:       SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Unknown null space type");
688:     }
689:   }
690:   PetscStackCallExternalVoid("ML_Aggregate_Set_MaxCoarseSize", ML_Aggregate_Set_MaxCoarseSize(agg_object, pc_ml->MaxCoarseSize));
691:   /* set options */
692:   switch (pc_ml->CoarsenScheme) {
693:   case 1:
694:     PetscStackCallExternalVoid("ML_Aggregate_Set_CoarsenScheme_Coupled", ML_Aggregate_Set_CoarsenScheme_Coupled(agg_object));
695:     break;
696:   case 2:
697:     PetscStackCallExternalVoid("ML_Aggregate_Set_CoarsenScheme_MIS", ML_Aggregate_Set_CoarsenScheme_MIS(agg_object));
698:     break;
699:   case 3:
700:     PetscStackCallExternalVoid("ML_Aggregate_Set_CoarsenScheme_METIS", ML_Aggregate_Set_CoarsenScheme_METIS(agg_object));
701:     break;
702:   }
703:   PetscStackCallExternalVoid("ML_Aggregate_Set_Threshold", ML_Aggregate_Set_Threshold(agg_object, pc_ml->Threshold));
704:   PetscStackCallExternalVoid("ML_Aggregate_Set_DampingFactor", ML_Aggregate_Set_DampingFactor(agg_object, pc_ml->DampingFactor));
705:   if (pc_ml->SpectralNormScheme_Anorm) PetscStackCallExternalVoid("ML_Set_SpectralNormScheme_Anorm", ML_Set_SpectralNormScheme_Anorm(ml_object));
706:   agg_object->keep_agg_information      = (int)pc_ml->KeepAggInfo;
707:   agg_object->keep_P_tentative          = (int)pc_ml->Reusable;
708:   agg_object->block_scaled_SA           = (int)pc_ml->BlockScaling;
709:   agg_object->minimizing_energy         = (int)pc_ml->EnergyMinimization;
710:   agg_object->minimizing_energy_droptol = (double)pc_ml->EnergyMinimizationDropTol;
711:   agg_object->cheap_minimizing_energy   = (int)pc_ml->EnergyMinimizationCheap;

713:   if (pc_ml->Aux) {
715:     ml_object->Amat[0].aux_data->threshold = pc_ml->AuxThreshold;
716:     ml_object->Amat[0].aux_data->enable    = 1;
717:     ml_object->Amat[0].aux_data->max_level = 10;
718:     ml_object->Amat[0].num_PDEs            = bs;
719:   }

721:   MatGetInfo(A, MAT_LOCAL, &info);
722:   ml_object->Amat[0].N_nonzeros = (int)info.nz_used;

724:   if (pc_ml->dim) {
725:     PetscInt                i, dim = pc_ml->dim;
726:     ML_Aggregate_Viz_Stats *grid_info;
727:     PetscInt                nlocghost;

729:     MatGetBlockSize(A, &bs);
730:     nlocghost = Aloc->cmap->n / bs;

732:     PetscStackCallExternalVoid("ML_Aggregate_VizAndStats_Setup(", ML_Aggregate_VizAndStats_Setup(ml_object)); /* create ml info for coords */
733:     grid_info = (ML_Aggregate_Viz_Stats *)ml_object->Grid[0].Grid;
734:     for (i = 0; i < dim; i++) {
735:       /* set the finest level coordinates to point to the column-order array
736:        * in pc_ml */
737:       /* NOTE: must point away before VizAndStats_Clean so ML doesn't free */
738:       switch (i) {
739:       case 0:
740:         grid_info->x = pc_ml->coords + nlocghost * i;
741:         break;
742:       case 1:
743:         grid_info->y = pc_ml->coords + nlocghost * i;
744:         break;
745:       case 2:
746:         grid_info->z = pc_ml->coords + nlocghost * i;
747:         break;
748:       default:
749:         SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_SIZ, "PCML coordinate dimension must be <= 3");
750:       }
751:     }
752:     grid_info->Ndim = dim;
753:   }

755:   /* repartitioning */
756:   if (pc_ml->Repartition) {
757:     PetscStackCallExternalVoid("ML_Repartition_Activate", ML_Repartition_Activate(ml_object));
758:     PetscStackCallExternalVoid("ML_Repartition_Set_LargestMinMaxRatio", ML_Repartition_Set_LargestMinMaxRatio(ml_object, pc_ml->MaxMinRatio));
759:     PetscStackCallExternalVoid("ML_Repartition_Set_MinPerProc", ML_Repartition_Set_MinPerProc(ml_object, pc_ml->MinPerProc));
760:     PetscStackCallExternalVoid("ML_Repartition_Set_PutOnSingleProc", ML_Repartition_Set_PutOnSingleProc(ml_object, pc_ml->PutOnSingleProc));
761: #if 0 /* Function not yet defined in ml-6.2 */
762:     /* I'm not sure what compatibility issues might crop up if we partitioned
763:      * on the finest level, so to be safe repartition starts on the next
764:      * finest level (reflection default behavior in
765:      * ml_MultiLevelPreconditioner) */
766:     PetscStackCallExternalVoid("ML_Repartition_Set_StartLevel",ML_Repartition_Set_StartLevel(ml_object,1));
767: #endif

769:     if (!pc_ml->RepartitionType) {
770:       PetscInt i;

773:       PetscStackCallExternalVoid("ML_Repartition_Set_Partitioner", ML_Repartition_Set_Partitioner(ml_object, ML_USEZOLTAN));
774:       PetscStackCallExternalVoid("ML_Aggregate_Set_Dimensions", ML_Aggregate_Set_Dimensions(agg_object, pc_ml->dim));

776:       for (i = 0; i < ml_object->ML_num_levels; i++) {
777:         ML_Aggregate_Viz_Stats *grid_info = (ML_Aggregate_Viz_Stats *)ml_object->Grid[i].Grid;
778:         grid_info->zoltan_type            = pc_ml->ZoltanScheme + 1; /* ml numbers options 1, 2, 3 */
779:         /* defaults from ml_agg_info.c */
780:         grid_info->zoltan_estimated_its = 40; /* only relevant to hypergraph / fast hypergraph */
781:         grid_info->zoltan_timers        = 0;
782:         grid_info->smoothing_steps      = 4; /* only relevant to hypergraph / fast hypergraph */
783:       }
784:     } else {
785:       PetscStackCallExternalVoid("ML_Repartition_Set_Partitioner", ML_Repartition_Set_Partitioner(ml_object, ML_USEPARMETIS));
786:     }
787:   }

789:   if (pc_ml->OldHierarchy) {
790:     PetscStackCallExternalVoid("ML_Gen_MGHierarchy_UsingAggregation", Nlevels = ML_Gen_MGHierarchy_UsingAggregation(ml_object, 0, ML_INCREASING, agg_object));
791:   } else {
792:     PetscStackCallExternalVoid("ML_Gen_MultiLevelHierarchy_UsingAggregation", Nlevels = ML_Gen_MultiLevelHierarchy_UsingAggregation(ml_object, 0, ML_INCREASING, agg_object));
793:   }
795:   pc_ml->Nlevels = Nlevels;
796:   fine_level     = Nlevels - 1;

798:   PCMGSetLevels(pc, Nlevels, NULL);
799:   /* set default smoothers */
800:   for (level = 1; level <= fine_level; level++) {
801:     PCMGGetSmoother(pc, level, &smoother);
802:     KSPSetType(smoother, KSPRICHARDSON);
803:     KSPGetPC(smoother, &subpc);
804:     PCSetType(subpc, PCSOR);
805:   }
806:   PetscObjectOptionsBegin((PetscObject)pc);
807:   PCSetFromOptions_MG(pc, PetscOptionsObject); /* should be called in PCSetFromOptions_ML(), but cannot be called prior to PCMGSetLevels() */
808:   PetscOptionsEnd();

810:   PetscMalloc1(Nlevels, &gridctx);

812:   pc_ml->gridctx = gridctx;

814:   /* wrap ML matrices by PETSc shell matrices at coarsened grids.
815:      Level 0 is the finest grid for ML, but coarsest for PETSc! */
816:   gridctx[fine_level].A = A;

818:   level = fine_level - 1;
819:   /* TODO: support for GPUs */
820:   if (size == 1) { /* convert ML P, R and A into seqaij format */
821:     for (mllevel = 1; mllevel < Nlevels; mllevel++) {
822:       mlmat = &(ml_object->Pmat[mllevel]);
823:       MatWrapML_SeqAIJ(mlmat, MAT_INITIAL_MATRIX, &gridctx[level].P);
824:       mlmat = &(ml_object->Rmat[mllevel - 1]);
825:       MatWrapML_SeqAIJ(mlmat, MAT_INITIAL_MATRIX, &gridctx[level].R);

827:       mlmat = &(ml_object->Amat[mllevel]);
828:       MatWrapML_SeqAIJ(mlmat, MAT_INITIAL_MATRIX, &gridctx[level].A);
829:       level--;
830:     }
831:   } else { /* convert ML P and R into shell format, ML A into mpiaij format */
832:     for (mllevel = 1; mllevel < Nlevels; mllevel++) {
833:       mlmat = &(ml_object->Pmat[mllevel]);
834:       MatWrapML_SHELL(mlmat, MAT_INITIAL_MATRIX, &gridctx[level].P);
835:       mlmat = &(ml_object->Rmat[mllevel - 1]);
836:       MatWrapML_SHELL(mlmat, MAT_INITIAL_MATRIX, &gridctx[level].R);

838:       mlmat = &(ml_object->Amat[mllevel]);
839:       MatWrapML_MPIAIJ(mlmat, MAT_INITIAL_MATRIX, &gridctx[level].A);
840:       level--;
841:     }
842:   }

844:   /* create vectors and ksp at all levels */
845:   for (level = 0; level < fine_level; level++) {
846:     level1 = level + 1;

848:     MatCreateVecs(gridctx[level].A, &gridctx[level].x, &gridctx[level].b);
849:     MatCreateVecs(gridctx[level1].A, NULL, &gridctx[level1].r);
850:     PCMGSetX(pc, level, gridctx[level].x);
851:     PCMGSetRhs(pc, level, gridctx[level].b);
852:     PCMGSetR(pc, level1, gridctx[level1].r);

854:     if (level == 0) {
855:       PCMGGetCoarseSolve(pc, &gridctx[level].ksp);
856:     } else {
857:       PCMGGetSmoother(pc, level, &gridctx[level].ksp);
858:     }
859:   }
860:   PCMGGetSmoother(pc, fine_level, &gridctx[fine_level].ksp);

862:   /* create coarse level and the interpolation between the levels */
863:   for (level = 0; level < fine_level; level++) {
864:     level1 = level + 1;

866:     PCMGSetInterpolation(pc, level1, gridctx[level].P);
867:     PCMGSetRestriction(pc, level1, gridctx[level].R);
868:     if (level > 0) PCMGSetResidual(pc, level, PCMGResidualDefault, gridctx[level].A);
869:     KSPSetOperators(gridctx[level].ksp, gridctx[level].A, gridctx[level].A);
870:   }
871:   PCMGSetResidual(pc, fine_level, PCMGResidualDefault, gridctx[fine_level].A);
872:   KSPSetOperators(gridctx[fine_level].ksp, gridctx[level].A, gridctx[fine_level].A);

874:   /* put coordinate info in levels */
875:   if (pc_ml->dim) {
876:     PetscInt   i, j, dim = pc_ml->dim;
877:     PetscInt   bs, nloc;
878:     PC         subpc;
879:     PetscReal *array;

881:     level = fine_level;
882:     for (mllevel = 0; mllevel < Nlevels; mllevel++) {
883:       ML_Aggregate_Viz_Stats *grid_info = (ML_Aggregate_Viz_Stats *)ml_object->Amat[mllevel].to->Grid->Grid;
884:       MPI_Comm                comm      = ((PetscObject)gridctx[level].A)->comm;

886:       MatGetBlockSize(gridctx[level].A, &bs);
887:       MatGetLocalSize(gridctx[level].A, NULL, &nloc);
888:       nloc /= bs; /* number of local nodes */

890:       VecCreate(comm, &gridctx[level].coords);
891:       VecSetSizes(gridctx[level].coords, dim * nloc, PETSC_DECIDE);
892:       VecSetType(gridctx[level].coords, VECMPI);
893:       VecGetArray(gridctx[level].coords, &array);
894:       for (j = 0; j < nloc; j++) {
895:         for (i = 0; i < dim; i++) {
896:           switch (i) {
897:           case 0:
898:             array[dim * j + i] = grid_info->x[j];
899:             break;
900:           case 1:
901:             array[dim * j + i] = grid_info->y[j];
902:             break;
903:           case 2:
904:             array[dim * j + i] = grid_info->z[j];
905:             break;
906:           default:
907:             SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_SIZ, "PCML coordinate dimension must be <= 3");
908:           }
909:         }
910:       }

912:       /* passing coordinates to smoothers/coarse solver, should they need them */
913:       KSPGetPC(gridctx[level].ksp, &subpc);
914:       PCSetCoordinates(subpc, dim, nloc, array);
915:       VecRestoreArray(gridctx[level].coords, &array);
916:       level--;
917:     }
918:   }

920:   /* setupcalled is set to 0 so that MG is setup from scratch */
921:   pc->setupcalled = 0;
922:   PCSetUp_MG(pc);
923:   return 0;
924: }

926: /* -------------------------------------------------------------------------- */
927: /*
928:    PCDestroy_ML - Destroys the private context for the ML preconditioner
929:    that was created with PCCreate_ML().

931:    Input Parameter:
932: .  pc - the preconditioner context

934:    Application Interface Routine: PCDestroy()
935: */
936: PetscErrorCode PCDestroy_ML(PC pc)
937: {
938:   PC_MG *mg    = (PC_MG *)pc->data;
939:   PC_ML *pc_ml = (PC_ML *)mg->innerctx;

941:   PCReset_ML(pc);
942:   PetscFree(pc_ml);
943:   PCDestroy_MG(pc);
944:   PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL);
945:   return 0;
946: }

948: PetscErrorCode PCSetFromOptions_ML(PC pc, PetscOptionItems *PetscOptionsObject)
949: {
950:   PetscInt    indx, PrintLevel, partindx;
951:   const char *scheme[] = {"Uncoupled", "Coupled", "MIS", "METIS"};
952:   const char *part[]   = {"Zoltan", "ParMETIS"};
953: #if defined(HAVE_ML_ZOLTAN)
954:   const char *zscheme[] = {"RCB", "hypergraph", "fast_hypergraph"};
955: #endif
956:   PC_MG      *mg    = (PC_MG *)pc->data;
957:   PC_ML      *pc_ml = (PC_ML *)mg->innerctx;
958:   PetscMPIInt size;
959:   MPI_Comm    comm;

961:   PetscObjectGetComm((PetscObject)pc, &comm);
962:   MPI_Comm_size(comm, &size);
963:   PetscOptionsHeadBegin(PetscOptionsObject, "ML options");

965:   PrintLevel = 0;
966:   indx       = 0;
967:   partindx   = 0;

969:   PetscOptionsInt("-pc_ml_PrintLevel", "Print level", "ML_Set_PrintLevel", PrintLevel, &PrintLevel, NULL);
970:   PetscStackCallExternalVoid("ML_Set_PrintLevel", ML_Set_PrintLevel(PrintLevel));
971:   PetscOptionsInt("-pc_ml_maxNlevels", "Maximum number of levels", "None", pc_ml->MaxNlevels, &pc_ml->MaxNlevels, NULL);
972:   PetscOptionsInt("-pc_ml_maxCoarseSize", "Maximum coarsest mesh size", "ML_Aggregate_Set_MaxCoarseSize", pc_ml->MaxCoarseSize, &pc_ml->MaxCoarseSize, NULL);
973:   PetscOptionsEList("-pc_ml_CoarsenScheme", "Aggregate Coarsen Scheme", "ML_Aggregate_Set_CoarsenScheme_*", scheme, 4, scheme[0], &indx, NULL);

975:   pc_ml->CoarsenScheme = indx;

977:   PetscOptionsReal("-pc_ml_DampingFactor", "P damping factor", "ML_Aggregate_Set_DampingFactor", pc_ml->DampingFactor, &pc_ml->DampingFactor, NULL);
978:   PetscOptionsReal("-pc_ml_Threshold", "Smoother drop tol", "ML_Aggregate_Set_Threshold", pc_ml->Threshold, &pc_ml->Threshold, NULL);
979:   PetscOptionsBool("-pc_ml_SpectralNormScheme_Anorm", "Method used for estimating spectral radius", "ML_Set_SpectralNormScheme_Anorm", pc_ml->SpectralNormScheme_Anorm, &pc_ml->SpectralNormScheme_Anorm, NULL);
980:   PetscOptionsBool("-pc_ml_Symmetrize", "Symmetrize aggregation", "ML_Set_Symmetrize", pc_ml->Symmetrize, &pc_ml->Symmetrize, NULL);
981:   PetscOptionsBool("-pc_ml_BlockScaling", "Scale all dofs at each node together", "None", pc_ml->BlockScaling, &pc_ml->BlockScaling, NULL);
982:   PetscOptionsEnum("-pc_ml_nullspace", "Which type of null space information to use", "None", PCMLNullSpaceTypes, (PetscEnum)pc_ml->nulltype, (PetscEnum *)&pc_ml->nulltype, NULL);
983:   PetscOptionsInt("-pc_ml_EnergyMinimization", "Energy minimization norm type (0=no minimization; see ML manual for 1,2,3; -1 and 4 undocumented)", "None", pc_ml->EnergyMinimization, &pc_ml->EnergyMinimization, NULL);
984:   PetscOptionsBool("-pc_ml_reuse_interpolation", "Reuse the interpolation operators when possible (cheaper, weaker when matrix entries change a lot)", "None", pc_ml->reuse_interpolation, &pc_ml->reuse_interpolation, NULL);
985:   /*
986:     The following checks a number of conditions.  If we let this stuff slip by, then ML's error handling will take over.
987:     This is suboptimal because it amounts to calling exit(1) so we check for the most common conditions.

989:     We also try to set some sane defaults when energy minimization is activated, otherwise it's hard to find a working
990:     combination of options and ML's exit(1) explanations don't help matters.
991:   */
994:   if (pc_ml->EnergyMinimization == 4) PetscInfo(pc, "Mandel's energy minimization scheme is experimental and broken in ML-6.2\n");
995:   if (pc_ml->EnergyMinimization) PetscOptionsReal("-pc_ml_EnergyMinimizationDropTol", "Energy minimization drop tolerance", "None", pc_ml->EnergyMinimizationDropTol, &pc_ml->EnergyMinimizationDropTol, NULL);
996:   if (pc_ml->EnergyMinimization == 2) {
997:     /* According to ml_MultiLevelPreconditioner.cpp, this option is only meaningful for norm type (2) */
998:     PetscOptionsBool("-pc_ml_EnergyMinimizationCheap", "Use cheaper variant of norm type 2", "None", pc_ml->EnergyMinimizationCheap, &pc_ml->EnergyMinimizationCheap, NULL);
999:   }
1000:   /* energy minimization sometimes breaks if this is turned off, the more classical stuff should be okay without it */
1001:   if (pc_ml->EnergyMinimization) pc_ml->KeepAggInfo = PETSC_TRUE;
1002:   PetscOptionsBool("-pc_ml_KeepAggInfo", "Allows the preconditioner to be reused, or auxiliary matrices to be generated", "None", pc_ml->KeepAggInfo, &pc_ml->KeepAggInfo, NULL);
1003:   /* Option (-1) doesn't work at all (calls exit(1)) if the tentative restriction operator isn't stored. */
1004:   if (pc_ml->EnergyMinimization == -1) pc_ml->Reusable = PETSC_TRUE;
1005:   PetscOptionsBool("-pc_ml_Reusable", "Store intermedaiate data structures so that the multilevel hierarchy is reusable", "None", pc_ml->Reusable, &pc_ml->Reusable, NULL);
1006:   /*
1007:     ML's C API is severely underdocumented and lacks significant functionality.  The C++ API calls
1008:     ML_Gen_MultiLevelHierarchy_UsingAggregation() which is a modified copy (!?) of the documented function
1009:     ML_Gen_MGHierarchy_UsingAggregation().  This modification, however, does not provide a strict superset of the
1010:     functionality in the old function, so some users may still want to use it.  Note that many options are ignored in
1011:     this context, but ML doesn't provide a way to find out which ones.
1012:    */
1013:   PetscOptionsBool("-pc_ml_OldHierarchy", "Use old routine to generate hierarchy", "None", pc_ml->OldHierarchy, &pc_ml->OldHierarchy, NULL);
1014:   PetscOptionsBool("-pc_ml_repartition", "Allow ML to repartition levels of the hierarchy", "ML_Repartition_Activate", pc_ml->Repartition, &pc_ml->Repartition, NULL);
1015:   if (pc_ml->Repartition) {
1016:     PetscOptionsReal("-pc_ml_repartitionMaxMinRatio", "Acceptable ratio of repartitioned sizes", "ML_Repartition_Set_LargestMinMaxRatio", pc_ml->MaxMinRatio, &pc_ml->MaxMinRatio, NULL);
1017:     PetscOptionsInt("-pc_ml_repartitionMinPerProc", "Smallest repartitioned size", "ML_Repartition_Set_MinPerProc", pc_ml->MinPerProc, &pc_ml->MinPerProc, NULL);
1018:     PetscOptionsInt("-pc_ml_repartitionPutOnSingleProc", "Problem size automatically repartitioned to one processor", "ML_Repartition_Set_PutOnSingleProc", pc_ml->PutOnSingleProc, &pc_ml->PutOnSingleProc, NULL);
1019: #if defined(HAVE_ML_ZOLTAN)
1020:     partindx = 0;
1021:     PetscOptionsEList("-pc_ml_repartitionType", "Repartitioning library to use", "ML_Repartition_Set_Partitioner", part, 2, part[0], &partindx, NULL);

1023:     pc_ml->RepartitionType = partindx;
1024:     if (!partindx) {
1025:       PetscInt zindx = 0;

1027:       PetscOptionsEList("-pc_ml_repartitionZoltanScheme", "Repartitioning scheme to use", "None", zscheme, 3, zscheme[0], &zindx, NULL);

1029:       pc_ml->ZoltanScheme = zindx;
1030:     }
1031: #else
1032:     partindx = 1;
1033:     PetscOptionsEList("-pc_ml_repartitionType", "Repartitioning library to use", "ML_Repartition_Set_Partitioner", part, 2, part[1], &partindx, NULL);
1034:     pc_ml->RepartitionType = partindx;
1036: #endif
1037:     PetscOptionsBool("-pc_ml_Aux", "Aggregate using auxiliary coordinate-based laplacian", "None", pc_ml->Aux, &pc_ml->Aux, NULL);
1038:     PetscOptionsReal("-pc_ml_AuxThreshold", "Auxiliary smoother drop tol", "None", pc_ml->AuxThreshold, &pc_ml->AuxThreshold, NULL);
1039:   }
1040:   PetscOptionsHeadEnd();
1041:   return 0;
1042: }

1044: /*
1045:    PCCreate_ML - Creates a ML preconditioner context, PC_ML,
1046:    and sets this as the private data within the generic preconditioning
1047:    context, PC, that was created within PCCreate().

1049:    Input Parameter:
1050: .  pc - the preconditioner context

1052:    Application Interface Routine: PCCreate()
1053: */

1055: /*MC
1056:      PCML - Use the SNL ML algebraic multigrid preconditioner.

1058:    Options Database Keys:
1059:    Multigrid options(inherited):
1060: +  -pc_mg_cycles <1> - 1 for V cycle, 2 for W-cycle (`MGSetCycles()`)
1061: .  -pc_mg_distinct_smoothup - Should one configure the up and down smoothers separately (`PCMGSetDistinctSmoothUp()`)
1062: -  -pc_mg_type <multiplicative> - (one of) additive multiplicative full kascade

1064:    ML Options Database Key:
1065: +  -pc_ml_PrintLevel <0> - Print level (`ML_Set_PrintLevel()`)
1066: .  -pc_ml_maxNlevels <10> - Maximum number of levels (None)
1067: .  -pc_ml_maxCoarseSize <1> - Maximum coarsest mesh size (`ML_Aggregate_Set_MaxCoarseSize()`)
1068: .  -pc_ml_CoarsenScheme <Uncoupled> - (one of) Uncoupled Coupled MIS METIS
1069: .  -pc_ml_DampingFactor <1.33333> - P damping factor (`ML_Aggregate_Set_DampingFactor()`)
1070: .  -pc_ml_Threshold <0> - Smoother drop tol (`ML_Aggregate_Set_Threshold()`)
1071: .  -pc_ml_SpectralNormScheme_Anorm <false> - Method used for estimating spectral radius (`ML_Set_SpectralNormScheme_Anorm()`)
1072: .  -pc_ml_repartition <false> - Allow ML to repartition levels of the hierarchy (`ML_Repartition_Activate()`)
1073: .  -pc_ml_repartitionMaxMinRatio <1.3> - Acceptable ratio of repartitioned sizes (`ML_Repartition_Set_LargestMinMaxRatio()`)
1074: .  -pc_ml_repartitionMinPerProc <512> - Smallest repartitioned size (`ML_Repartition_Set_MinPerProc()`)
1075: .  -pc_ml_repartitionPutOnSingleProc <5000> - Problem size automatically repartitioned to one processor (`ML_Repartition_Set_PutOnSingleProc()`)
1076: .  -pc_ml_repartitionType <Zoltan> - Repartitioning library to use (`ML_Repartition_Set_Partitioner()`)
1077: .  -pc_ml_repartitionZoltanScheme <RCB> - Repartitioning scheme to use (None)
1078: .  -pc_ml_Aux <false> - Aggregate using auxiliary coordinate-based Laplacian (None)
1079: -  -pc_ml_AuxThreshold <0.0> - Auxiliary smoother drop tol (None)

1081:    Level: intermediate

1083:    Developer Note:
1084:    The coarser grid matrices and restriction/interpolation
1085:    operators are computed by ML, with the matrices converted to PETSc matrices in `MATAIJ` format
1086:    and the restriction/interpolation operators wrapped as PETSc shell matrices.

1088: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCMGType`, `PCMG`, `PCHYPRE`, `PCGAMG`,
1089:           `PCMGSetLevels()`, `PCMGGetLevels()`, `PCMGSetType()`, `MPSetCycles()`, `PCMGSetDistinctSmoothUp()`,
1090:           `PCMGGetCoarseSolve()`, `PCMGSetResidual()`, `PCMGSetInterpolation()`,
1091:           `PCMGSetRestriction()`, `PCMGGetSmoother()`, `PCMGGetSmootherUp()`, `PCMGGetSmootherDown()`,
1092:           `PCMGSetCycleTypeOnLevel()`, `PCMGSetRhs()`, `PCMGSetX()`, `PCMGSetR()`
1093: M*/

1095: PETSC_EXTERN PetscErrorCode PCCreate_ML(PC pc)
1096: {
1097:   PC_ML *pc_ml;
1098:   PC_MG *mg;

1100:   /* PCML is an inherited class of PCMG. Initialize pc as PCMG */
1101:   PCSetType(pc, PCMG); /* calls PCCreate_MG() and MGCreate_Private() */
1102:   PetscObjectChangeTypeName((PetscObject)pc, PCML);
1103:   /* Since PCMG tries to use DM associated with PC must delete it */
1104:   DMDestroy(&pc->dm);
1105:   PCMGSetGalerkin(pc, PC_MG_GALERKIN_EXTERNAL);
1106:   mg = (PC_MG *)pc->data;

1108:   /* create a supporting struct and attach it to pc */
1109:   PetscNew(&pc_ml);
1110:   mg->innerctx = pc_ml;

1112:   pc_ml->ml_object                = 0;
1113:   pc_ml->agg_object               = 0;
1114:   pc_ml->gridctx                  = 0;
1115:   pc_ml->PetscMLdata              = 0;
1116:   pc_ml->Nlevels                  = -1;
1117:   pc_ml->MaxNlevels               = 10;
1118:   pc_ml->MaxCoarseSize            = 1;
1119:   pc_ml->CoarsenScheme            = 1;
1120:   pc_ml->Threshold                = 0.0;
1121:   pc_ml->DampingFactor            = 4.0 / 3.0;
1122:   pc_ml->SpectralNormScheme_Anorm = PETSC_FALSE;
1123:   pc_ml->size                     = 0;
1124:   pc_ml->dim                      = 0;
1125:   pc_ml->nloc                     = 0;
1126:   pc_ml->coords                   = 0;
1127:   pc_ml->Repartition              = PETSC_FALSE;
1128:   pc_ml->MaxMinRatio              = 1.3;
1129:   pc_ml->MinPerProc               = 512;
1130:   pc_ml->PutOnSingleProc          = 5000;
1131:   pc_ml->RepartitionType          = 0;
1132:   pc_ml->ZoltanScheme             = 0;
1133:   pc_ml->Aux                      = PETSC_FALSE;
1134:   pc_ml->AuxThreshold             = 0.0;

1136:   /* allow for coordinates to be passed */
1137:   PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_ML);

1139:   /* overwrite the pointers of PCMG by the functions of PCML */
1140:   pc->ops->setfromoptions = PCSetFromOptions_ML;
1141:   pc->ops->setup          = PCSetUp_ML;
1142:   pc->ops->reset          = PCReset_ML;
1143:   pc->ops->destroy        = PCDestroy_ML;
1144:   return 0;
1145: }