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:   PetscFunctionBegin;
104:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
105:   if (size == 1) PetscFunctionReturn(PETSC_SUCCESS);

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

117: /*
118:   Needed since ML expects an int (*)(double *, void *) but PetscErrorCode may be an
119:   enum. Instead of modifying PetscML_comm() it is easier to just wrap it
120: */
121: static int ML_PetscML_comm(double p[], void *ML_data)
122: {
123:   return (int)PetscML_comm(p, ML_data);
124: }

126: static int PetscML_matvec(ML_Operator *ML_data, int in_length, double p[], int out_length, double ap[])
127: {
128:   FineGridCtx *ml = (FineGridCtx *)ML_Get_MyMatvecData(ML_data);
129:   Mat          A = ml->A, Aloc = ml->Aloc;
130:   PetscMPIInt  size;
131:   PetscScalar *pwork = ml->pwork;
132:   PetscInt     i;

134:   PetscFunctionBegin;
135:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
136:   if (size == 1) {
137:     PetscCall(VecPlaceArray(ml->x, p));
138:   } else {
139:     for (i = 0; i < in_length; i++) pwork[i] = p[i];
140:     PetscCall(PetscML_comm(pwork, ml));
141:     PetscCall(VecPlaceArray(ml->x, pwork));
142:   }
143:   PetscCall(VecPlaceArray(ml->y, ap));
144:   PetscCall(MatMult(Aloc, ml->x, ml->y));
145:   PetscCall(VecResetArray(ml->x));
146:   PetscCall(VecResetArray(ml->y));
147:   PetscFunctionReturn(PETSC_SUCCESS);
148: }

150: static PetscErrorCode MatMult_ML(Mat A, Vec x, Vec y)
151: {
152:   Mat_MLShell       *shell;
153:   PetscScalar       *yarray;
154:   const PetscScalar *xarray;
155:   PetscInt           x_length, y_length;

157:   PetscFunctionBegin;
158:   PetscCall(MatShellGetContext(A, &shell));
159:   PetscCall(VecGetArrayRead(x, &xarray));
160:   PetscCall(VecGetArray(y, &yarray));
161:   x_length = shell->mlmat->invec_leng;
162:   y_length = shell->mlmat->outvec_leng;
163:   PetscStackCallExternalVoid("ML_Operator_Apply", ML_Operator_Apply(shell->mlmat, x_length, (PetscScalar *)xarray, y_length, yarray));
164:   PetscCall(VecRestoreArrayRead(x, &xarray));
165:   PetscCall(VecRestoreArray(y, &yarray));
166:   PetscFunctionReturn(PETSC_SUCCESS);
167: }

169: /* newtype is ignored since only handles one case */
170: static PetscErrorCode MatConvert_MPIAIJ_ML(Mat A, MatType newtype, MatReuse scall, Mat *Aloc)
171: {
172:   Mat_MPIAIJ  *mpimat = (Mat_MPIAIJ *)A->data;
173:   Mat_SeqAIJ  *mat, *a = (Mat_SeqAIJ *)(mpimat->A)->data, *b = (Mat_SeqAIJ *)(mpimat->B)->data;
174:   PetscInt    *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j;
175:   PetscScalar *aa, *ba, *ca;
176:   PetscInt     am = A->rmap->n, an = A->cmap->n, i, j, k;
177:   PetscInt    *ci, *cj, ncols;

179:   PetscFunctionBegin;
180:   PetscCheck(am == an, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "A must have a square diagonal portion, am: %d != an: %d", am, an);
181:   PetscCall(MatSeqAIJGetArrayRead(mpimat->A, (const PetscScalar **)&aa));
182:   PetscCall(MatSeqAIJGetArrayRead(mpimat->B, (const PetscScalar **)&ba));
183:   if (scall == MAT_INITIAL_MATRIX) {
184:     PetscCall(PetscMalloc1(1 + am, &ci));
185:     ci[0] = 0;
186:     for (i = 0; i < am; i++) ci[i + 1] = ci[i] + (ai[i + 1] - ai[i]) + (bi[i + 1] - bi[i]);
187:     PetscCall(PetscMalloc1(1 + ci[am], &cj));
188:     PetscCall(PetscMalloc1(1 + ci[am], &ca));

190:     k = 0;
191:     for (i = 0; i < am; i++) {
192:       /* diagonal portion of A */
193:       ncols = ai[i + 1] - ai[i];
194:       for (j = 0; j < ncols; j++) {
195:         cj[k]   = *aj++;
196:         ca[k++] = *aa++;
197:       }
198:       /* off-diagonal portion of A */
199:       ncols = bi[i + 1] - bi[i];
200:       for (j = 0; j < ncols; j++) {
201:         cj[k] = an + (*bj);
202:         bj++;
203:         ca[k++] = *ba++;
204:       }
205:     }
206:     PetscCheck(k == ci[am], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "k: %d != ci[am]: %d", k, ci[am]);

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

212:     /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
213:     /* Since these are PETSc arrays, change flags to free them as necessary. */
214:     mat          = (Mat_SeqAIJ *)(*Aloc)->data;
215:     mat->free_a  = PETSC_TRUE;
216:     mat->free_ij = PETSC_TRUE;

218:     mat->nonew = 0;
219:   } else if (scall == MAT_REUSE_MATRIX) {
220:     mat = (Mat_SeqAIJ *)(*Aloc)->data;
221:     ci  = mat->i;
222:     cj  = mat->j;
223:     ca  = mat->a;
224:     for (i = 0; i < am; i++) {
225:       /* diagonal portion of A */
226:       ncols = ai[i + 1] - ai[i];
227:       for (j = 0; j < ncols; j++) *ca++ = *aa++;
228:       /* off-diagonal portion of A */
229:       ncols = bi[i + 1] - bi[i];
230:       for (j = 0; j < ncols; j++) *ca++ = *ba++;
231:     }
232:   } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid MatReuse %d", (int)scall);
233:   PetscCall(MatSeqAIJRestoreArrayRead(mpimat->A, (const PetscScalar **)&aa));
234:   PetscCall(MatSeqAIJRestoreArrayRead(mpimat->B, (const PetscScalar **)&ba));
235:   PetscFunctionReturn(PETSC_SUCCESS);
236: }

238: static PetscErrorCode MatDestroy_ML(Mat A)
239: {
240:   Mat_MLShell *shell;

242:   PetscFunctionBegin;
243:   PetscCall(MatShellGetContext(A, &shell));
244:   PetscCall(PetscFree(shell));
245:   PetscFunctionReturn(PETSC_SUCCESS);
246: }

248: static PetscErrorCode MatWrapML_SeqAIJ(ML_Operator *mlmat, MatReuse reuse, Mat *newmat)
249: {
250:   struct ML_CSR_MSRdata *matdata = (struct ML_CSR_MSRdata *)mlmat->data;
251:   PetscInt               m = mlmat->outvec_leng, n = mlmat->invec_leng, *nnz = NULL, nz_max;
252:   PetscInt              *ml_cols = matdata->columns, *ml_rowptr = matdata->rowptr, *aj, i;
253:   PetscScalar           *ml_vals = matdata->values, *aa;

255:   PetscFunctionBegin;
256:   PetscCheck(mlmat->getrow, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "mlmat->getrow = NULL");
257:   if (m != n) { /* ML Pmat and Rmat are in CSR format. Pass array pointers into SeqAIJ matrix */
258:     if (reuse) {
259:       Mat_SeqAIJ *aij = (Mat_SeqAIJ *)(*newmat)->data;
260:       aij->i          = ml_rowptr;
261:       aij->j          = ml_cols;
262:       aij->a          = ml_vals;
263:     } else {
264:       /* sort ml_cols and ml_vals */
265:       PetscCall(PetscMalloc1(m + 1, &nnz));
266:       for (i = 0; i < m; i++) nnz[i] = ml_rowptr[i + 1] - ml_rowptr[i];
267:       aj = ml_cols;
268:       aa = ml_vals;
269:       for (i = 0; i < m; i++) {
270:         PetscCall(PetscSortIntWithScalarArray(nnz[i], aj, aa));
271:         aj += nnz[i];
272:         aa += nnz[i];
273:       }
274:       PetscCall(MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, m, n, ml_rowptr, ml_cols, ml_vals, newmat));
275:       PetscCall(PetscFree(nnz));
276:     }
277:     PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY));
278:     PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY));
279:     PetscFunctionReturn(PETSC_SUCCESS);
280:   }

282:   nz_max = PetscMax(1, mlmat->max_nz_per_row);
283:   PetscCall(PetscMalloc2(nz_max, &aa, nz_max, &aj));
284:   if (!reuse) {
285:     PetscCall(MatCreate(PETSC_COMM_SELF, newmat));
286:     PetscCall(MatSetSizes(*newmat, m, n, PETSC_DECIDE, PETSC_DECIDE));
287:     PetscCall(MatSetType(*newmat, MATSEQAIJ));
288:     /* keep track of block size for A matrices */
289:     PetscCall(MatSetBlockSize(*newmat, mlmat->num_PDEs));

291:     PetscCall(PetscMalloc1(m, &nnz));
292:     for (i = 0; i < m; i++) PetscStackCallExternalVoid("ML_Operator_Getrow", ML_Operator_Getrow(mlmat, 1, &i, nz_max, aj, aa, &nnz[i]));
293:     PetscCall(MatSeqAIJSetPreallocation(*newmat, 0, nnz));
294:   }
295:   for (i = 0; i < m; i++) {
296:     PetscInt ncols;

298:     PetscStackCallExternalVoid("ML_Operator_Getrow", ML_Operator_Getrow(mlmat, 1, &i, nz_max, aj, aa, &ncols));
299:     PetscCall(MatSetValues(*newmat, 1, &i, ncols, aj, aa, INSERT_VALUES));
300:   }
301:   PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY));
302:   PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY));

304:   PetscCall(PetscFree2(aa, aj));
305:   PetscCall(PetscFree(nnz));
306:   PetscFunctionReturn(PETSC_SUCCESS);
307: }

309: static PetscErrorCode MatWrapML_SHELL(ML_Operator *mlmat, MatReuse reuse, Mat *newmat)
310: {
311:   PetscInt     m, n;
312:   ML_Comm     *MLcomm;
313:   Mat_MLShell *shellctx;

315:   PetscFunctionBegin;
316:   m = mlmat->outvec_leng;
317:   n = mlmat->invec_leng;

319:   if (reuse) {
320:     PetscCall(MatShellGetContext(*newmat, &shellctx));
321:     shellctx->mlmat = mlmat;
322:     PetscFunctionReturn(PETSC_SUCCESS);
323:   }

325:   MLcomm = mlmat->comm;

327:   PetscCall(PetscNew(&shellctx));
328:   PetscCall(MatCreateShell(MLcomm->USR_comm, m, n, PETSC_DETERMINE, PETSC_DETERMINE, shellctx, newmat));
329:   PetscCall(MatShellSetOperation(*newmat, MATOP_MULT, (void (*)(void))MatMult_ML));
330:   PetscCall(MatShellSetOperation(*newmat, MATOP_DESTROY, (void (*)(void))MatDestroy_ML));

332:   shellctx->A     = *newmat;
333:   shellctx->mlmat = mlmat;
334:   PetscFunctionReturn(PETSC_SUCCESS);
335: }

337: static PetscErrorCode MatWrapML_MPIAIJ(ML_Operator *mlmat, MatReuse reuse, Mat *newmat)
338: {
339:   PetscInt    *aj;
340:   PetscScalar *aa;
341:   PetscInt     i, j, *gordering;
342:   PetscInt     m = mlmat->outvec_leng, n, nz_max, row;
343:   Mat          A;

345:   PetscFunctionBegin;
346:   PetscCheck(mlmat->getrow, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "mlmat->getrow = NULL");
347:   n = mlmat->invec_leng;
348:   PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "m %d must equal to n %d", m, n);

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

353:   nz_max = PetscMax(1, mlmat->max_nz_per_row) + 1;
354:   PetscCall(PetscMalloc2(nz_max, &aa, nz_max, &aj));
355:   if (reuse) {
356:     A = *newmat;
357:   } else {
358:     PetscInt *nnzA, *nnzB, *nnz;
359:     PetscInt  rstart;
360:     PetscCall(MatCreate(mlmat->comm->USR_comm, &A));
361:     PetscCall(MatSetSizes(A, m, n, PETSC_DECIDE, PETSC_DECIDE));
362:     PetscCall(MatSetType(A, MATMPIAIJ));
363:     /* keep track of block size for A matrices */
364:     PetscCall(MatSetBlockSize(A, mlmat->num_PDEs));
365:     PetscCall(PetscMalloc3(m, &nnzA, m, &nnzB, m, &nnz));
366:     PetscCallMPI(MPI_Scan(&m, &rstart, 1, MPIU_INT, MPI_SUM, mlmat->comm->USR_comm));
367:     rstart -= m;

369:     for (i = 0; i < m; i++) {
370:       row = gordering[i] - rstart;
371:       PetscStackCallExternalVoid("ML_Operator_Getrow", ML_Operator_Getrow(mlmat, 1, &i, nz_max, aj, aa, &nnz[i]));
372:       nnzA[row] = 0;
373:       for (j = 0; j < nnz[i]; j++) {
374:         if (aj[j] < m) nnzA[row]++;
375:       }
376:       nnzB[row] = nnz[i] - nnzA[row];
377:     }
378:     PetscCall(MatMPIAIJSetPreallocation(A, 0, nnzA, 0, nnzB));
379:     PetscCall(PetscFree3(nnzA, nnzB, nnz));
380:   }
381:   for (i = 0; i < m; i++) {
382:     PetscInt ncols;
383:     row = gordering[i];

385:     PetscStackCallExternalVoid(",ML_Operator_Getrow", ML_Operator_Getrow(mlmat, 1, &i, nz_max, aj, aa, &ncols));
386:     for (j = 0; j < ncols; j++) aj[j] = gordering[aj[j]];
387:     PetscCall(MatSetValues(A, 1, &row, ncols, aj, aa, INSERT_VALUES));
388:   }
389:   PetscStackCallExternalVoid("ML_free", ML_free(gordering));
390:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
391:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
392:   *newmat = A;

394:   PetscCall(PetscFree2(aa, aj));
395:   PetscFunctionReturn(PETSC_SUCCESS);
396: }

398: /* -------------------------------------------------------------------------- */
399: /*
400:    PCSetCoordinates_ML

402:    Input Parameter:
403:    .  pc - the preconditioner context
404: */
405: static PetscErrorCode PCSetCoordinates_ML(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords)
406: {
407:   PC_MG   *mg    = (PC_MG *)pc->data;
408:   PC_ML   *pc_ml = (PC_ML *)mg->innerctx;
409:   PetscInt arrsz, oldarrsz, bs, my0, kk, ii, nloc, Iend, aloc;
410:   Mat      Amat = pc->pmat;

412:   /* this function copied and modified from PCSetCoordinates_GEO -TGI */
413:   PetscFunctionBegin;
415:   PetscCall(MatGetBlockSize(Amat, &bs));

417:   PetscCall(MatGetOwnershipRange(Amat, &my0, &Iend));
418:   aloc = (Iend - my0);
419:   nloc = (Iend - my0) / bs;

421:   PetscCheck((nloc == a_nloc) || (aloc == a_nloc), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of local blocks %" PetscInt_FMT " must be %" PetscInt_FMT " or %" PetscInt_FMT ".", a_nloc, nloc, aloc);

423:   oldarrsz    = pc_ml->dim * pc_ml->nloc;
424:   pc_ml->dim  = ndm;
425:   pc_ml->nloc = nloc;
426:   arrsz       = ndm * nloc;

428:   /* create data - syntactic sugar that should be refactored at some point */
429:   if (pc_ml->coords == 0 || (oldarrsz != arrsz)) {
430:     PetscCall(PetscFree(pc_ml->coords));
431:     PetscCall(PetscMalloc1(arrsz, &pc_ml->coords));
432:   }
433:   for (kk = 0; kk < arrsz; kk++) pc_ml->coords[kk] = -999.;
434:   /* copy data in - column oriented */
435:   if (nloc == a_nloc) {
436:     for (kk = 0; kk < nloc; kk++) {
437:       for (ii = 0; ii < ndm; ii++) pc_ml->coords[ii * nloc + kk] = coords[kk * ndm + ii];
438:     }
439:   } else { /* assumes the coordinates are blocked */
440:     for (kk = 0; kk < nloc; kk++) {
441:       for (ii = 0; ii < ndm; ii++) pc_ml->coords[ii * nloc + kk] = coords[bs * kk * ndm + ii];
442:     }
443:   }
444:   PetscFunctionReturn(PETSC_SUCCESS);
445: }

447: /* -----------------------------------------------------------------------------*/
448: extern PetscErrorCode PCReset_MG(PC);
449: PetscErrorCode        PCReset_ML(PC pc)
450: {
451:   PC_MG   *mg    = (PC_MG *)pc->data;
452:   PC_ML   *pc_ml = (PC_ML *)mg->innerctx;
453:   PetscInt level, fine_level = pc_ml->Nlevels - 1, dim = pc_ml->dim;

455:   PetscFunctionBegin;
456:   if (dim) {
457:     for (level = 0; level <= fine_level; level++) PetscCall(VecDestroy(&pc_ml->gridctx[level].coords));
458:     if (pc_ml->ml_object && pc_ml->ml_object->Grid) {
459:       ML_Aggregate_Viz_Stats *grid_info = (ML_Aggregate_Viz_Stats *)pc_ml->ml_object->Grid[0].Grid;
460:       grid_info->x                      = 0; /* do this so ML doesn't try to free coordinates */
461:       grid_info->y                      = 0;
462:       grid_info->z                      = 0;
463:       PetscStackCallExternalVoid("ML_Operator_Getrow", ML_Aggregate_VizAndStats_Clean(pc_ml->ml_object));
464:     }
465:   }
466:   PetscStackCallExternalVoid("ML_Aggregate_Destroy", ML_Aggregate_Destroy(&pc_ml->agg_object));
467:   PetscStackCallExternalVoid("ML_Aggregate_Destroy", ML_Destroy(&pc_ml->ml_object));

469:   if (pc_ml->PetscMLdata) {
470:     PetscCall(PetscFree(pc_ml->PetscMLdata->pwork));
471:     PetscCall(MatDestroy(&pc_ml->PetscMLdata->Aloc));
472:     PetscCall(VecDestroy(&pc_ml->PetscMLdata->x));
473:     PetscCall(VecDestroy(&pc_ml->PetscMLdata->y));
474:   }
475:   PetscCall(PetscFree(pc_ml->PetscMLdata));

477:   if (pc_ml->gridctx) {
478:     for (level = 0; level < fine_level; level++) {
479:       if (pc_ml->gridctx[level].A) PetscCall(MatDestroy(&pc_ml->gridctx[level].A));
480:       if (pc_ml->gridctx[level].P) PetscCall(MatDestroy(&pc_ml->gridctx[level].P));
481:       if (pc_ml->gridctx[level].R) PetscCall(MatDestroy(&pc_ml->gridctx[level].R));
482:       if (pc_ml->gridctx[level].x) PetscCall(VecDestroy(&pc_ml->gridctx[level].x));
483:       if (pc_ml->gridctx[level].b) PetscCall(VecDestroy(&pc_ml->gridctx[level].b));
484:       if (pc_ml->gridctx[level + 1].r) PetscCall(VecDestroy(&pc_ml->gridctx[level + 1].r));
485:     }
486:   }
487:   PetscCall(PetscFree(pc_ml->gridctx));
488:   PetscCall(PetscFree(pc_ml->coords));

490:   pc_ml->dim  = 0;
491:   pc_ml->nloc = 0;
492:   PetscCall(PCReset_MG(pc));
493:   PetscFunctionReturn(PETSC_SUCCESS);
494: }
495: /* -------------------------------------------------------------------------- */
496: /*
497:    PCSetUp_ML - Prepares for the use of the ML preconditioner
498:                     by setting data structures and options.

500:    Input Parameter:
501: .  pc - the preconditioner context

503:    Application Interface Routine: PCSetUp()

505:    Note:
506:    The interface routine PCSetUp() is not usually called directly by
507:    the user, but instead is called by PCApply() if necessary.
508: */
509: extern PetscErrorCode PCSetFromOptions_MG(PC, PetscOptionItems *PetscOptionsObject);
510: extern PetscErrorCode PCReset_MG(PC);

512: PetscErrorCode PCSetUp_ML(PC pc)
513: {
514:   PetscMPIInt      size;
515:   FineGridCtx     *PetscMLdata;
516:   ML              *ml_object;
517:   ML_Aggregate    *agg_object;
518:   ML_Operator     *mlmat;
519:   PetscInt         nlocal_allcols, Nlevels, mllevel, level, level1, m, fine_level, bs;
520:   Mat              A, Aloc;
521:   GridCtx         *gridctx;
522:   PC_MG           *mg    = (PC_MG *)pc->data;
523:   PC_ML           *pc_ml = (PC_ML *)mg->innerctx;
524:   PetscBool        isSeq, isMPI;
525:   KSP              smoother;
526:   PC               subpc;
527:   PetscInt         mesh_level, old_mesh_level;
528:   MatInfo          info;
529:   static PetscBool cite = PETSC_FALSE;

531:   PetscFunctionBegin;
532:   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 = "
533:                                    "{SAND2004-4821},\n  year = 2004\n}\n",
534:                                    &cite));
535:   A = pc->pmat;
536:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));

538:   if (pc->setupcalled) {
539:     if (pc->flag == SAME_NONZERO_PATTERN && pc_ml->reuse_interpolation) {
540:       /*
541:        Reuse interpolaton instead of recomputing aggregates and updating the whole hierarchy. This is less expensive for
542:        multiple solves in which the matrix is not changing too quickly.
543:        */
544:       ml_object             = pc_ml->ml_object;
545:       gridctx               = pc_ml->gridctx;
546:       Nlevels               = pc_ml->Nlevels;
547:       fine_level            = Nlevels - 1;
548:       gridctx[fine_level].A = A;

550:       PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &isSeq));
551:       PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &isMPI));
552:       PetscCheck(isMPI || isSeq, 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);
553:       if (isMPI) {
554:         PetscCall(MatConvert_MPIAIJ_ML(A, NULL, MAT_INITIAL_MATRIX, &Aloc));
555:       } else {
556:         Aloc = A;
557:         PetscCall(PetscObjectReference((PetscObject)Aloc));
558:       }

560:       PetscCall(MatGetSize(Aloc, &m, &nlocal_allcols));
561:       PetscMLdata = pc_ml->PetscMLdata;
562:       PetscCall(MatDestroy(&PetscMLdata->Aloc));
563:       PetscMLdata->A    = A;
564:       PetscMLdata->Aloc = Aloc;
565:       PetscStackCallExternalVoid("ML_Aggregate_Destroy", ML_Init_Amatrix(ml_object, 0, m, m, PetscMLdata));
566:       PetscStackCallExternalVoid("ML_Set_Amatrix_Matvec", ML_Set_Amatrix_Matvec(ml_object, 0, PetscML_matvec));

568:       mesh_level = ml_object->ML_finest_level;
569:       while (ml_object->SingleLevel[mesh_level].Rmat->to) {
570:         old_mesh_level = mesh_level;
571:         mesh_level     = ml_object->SingleLevel[mesh_level].Rmat->to->levelnum;

573:         /* clean and regenerate A */
574:         mlmat = &(ml_object->Amat[mesh_level]);
575:         PetscStackCallExternalVoid("ML_Operator_Clean", ML_Operator_Clean(mlmat));
576:         PetscStackCallExternalVoid("ML_Operator_Init", ML_Operator_Init(mlmat, ml_object->comm));
577:         PetscStackCallExternalVoid("ML_Gen_AmatrixRAP", ML_Gen_AmatrixRAP(ml_object, old_mesh_level, mesh_level));
578:       }

580:       level = fine_level - 1;
581:       if (size == 1) { /* convert ML P, R and A into seqaij format */
582:         for (mllevel = 1; mllevel < Nlevels; mllevel++) {
583:           mlmat = &(ml_object->Amat[mllevel]);
584:           PetscCall(MatWrapML_SeqAIJ(mlmat, MAT_REUSE_MATRIX, &gridctx[level].A));
585:           level--;
586:         }
587:       } else { /* convert ML P and R into shell format, ML A into mpiaij format */
588:         for (mllevel = 1; mllevel < Nlevels; mllevel++) {
589:           mlmat = &(ml_object->Amat[mllevel]);
590:           PetscCall(MatWrapML_MPIAIJ(mlmat, MAT_REUSE_MATRIX, &gridctx[level].A));
591:           level--;
592:         }
593:       }

595:       for (level = 0; level < fine_level; level++) {
596:         if (level > 0) PetscCall(PCMGSetResidual(pc, level, PCMGResidualDefault, gridctx[level].A));
597:         PetscCall(KSPSetOperators(gridctx[level].ksp, gridctx[level].A, gridctx[level].A));
598:       }
599:       PetscCall(PCMGSetResidual(pc, fine_level, PCMGResidualDefault, gridctx[fine_level].A));
600:       PetscCall(KSPSetOperators(gridctx[fine_level].ksp, gridctx[level].A, gridctx[fine_level].A));

602:       PetscCall(PCSetUp_MG(pc));
603:       PetscFunctionReturn(PETSC_SUCCESS);
604:     } else {
605:       /* since ML can change the size of vectors/matrices at any level we must destroy everything */
606:       PetscCall(PCReset_ML(pc));
607:     }
608:   }

610:   /* setup special features of PCML */
611:   /*--------------------------------*/
612:   /* convert A to Aloc to be used by ML at fine grid */
613:   pc_ml->size = size;
614:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &isSeq));
615:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &isMPI));
616:   PetscCheck(isMPI || isSeq, 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);
617:   if (isMPI) {
618:     PetscCall(MatConvert_MPIAIJ_ML(A, NULL, MAT_INITIAL_MATRIX, &Aloc));
619:   } else {
620:     Aloc = A;
621:     PetscCall(PetscObjectReference((PetscObject)Aloc));
622:   }

624:   /* create and initialize struct 'PetscMLdata' */
625:   PetscCall(PetscNew(&PetscMLdata));
626:   pc_ml->PetscMLdata = PetscMLdata;
627:   PetscCall(PetscMalloc1(Aloc->cmap->n + 1, &PetscMLdata->pwork));

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

631:   PetscMLdata->A    = A;
632:   PetscMLdata->Aloc = Aloc;
633:   if (pc_ml->dim) { /* create vecs around the coordinate data given */
634:     PetscInt   i, j, dim = pc_ml->dim;
635:     PetscInt   nloc = pc_ml->nloc, nlocghost;
636:     PetscReal *ghostedcoords;

638:     PetscCall(MatGetBlockSize(A, &bs));
639:     nlocghost = Aloc->cmap->n / bs;
640:     PetscCall(PetscMalloc1(dim * nlocghost, &ghostedcoords));
641:     for (i = 0; i < dim; i++) {
642:       /* copy coordinate values into first component of pwork */
643:       for (j = 0; j < nloc; j++) PetscMLdata->pwork[bs * j] = pc_ml->coords[nloc * i + j];
644:       /* get the ghost values */
645:       PetscCall(PetscML_comm(PetscMLdata->pwork, PetscMLdata));
646:       /* write into the vector */
647:       for (j = 0; j < nlocghost; j++) ghostedcoords[i * nlocghost + j] = PetscMLdata->pwork[bs * j];
648:     }
649:     /* replace the original coords with the ghosted coords, because these are
650:      * what ML needs */
651:     PetscCall(PetscFree(pc_ml->coords));
652:     pc_ml->coords = ghostedcoords;
653:   }

655:   /* create ML discretization matrix at fine grid */
656:   /* ML requires input of fine-grid matrix. It determines nlevels. */
657:   PetscCall(MatGetSize(Aloc, &m, &nlocal_allcols));
658:   PetscCall(MatGetBlockSize(A, &bs));
659:   PetscStackCallExternalVoid("ML_Create", ML_Create(&ml_object, pc_ml->MaxNlevels));
660:   PetscStackCallExternalVoid("ML_Comm_Set_UsrComm", ML_Comm_Set_UsrComm(ml_object->comm, PetscObjectComm((PetscObject)A)));
661:   pc_ml->ml_object = ml_object;
662:   PetscStackCallExternalVoid("ML_Init_Amatrix", ML_Init_Amatrix(ml_object, 0, m, m, PetscMLdata));
663:   PetscStackCallExternalVoid("ML_Set_Amatrix_Getrow", ML_Set_Amatrix_Getrow(ml_object, 0, PetscML_getrow, ML_PetscML_comm, nlocal_allcols));
664:   PetscStackCallExternalVoid("ML_Set_Amatrix_Matvec", ML_Set_Amatrix_Matvec(ml_object, 0, PetscML_matvec));

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

668:   /* aggregation */
669:   PetscStackCallExternalVoid("ML_Aggregate_Create", ML_Aggregate_Create(&agg_object));
670:   pc_ml->agg_object = agg_object;

672:   {
673:     MatNullSpace mnull;
674:     PetscCall(MatGetNearNullSpace(A, &mnull));
675:     if (pc_ml->nulltype == PCML_NULLSPACE_AUTO) {
676:       if (mnull) pc_ml->nulltype = PCML_NULLSPACE_USER;
677:       else if (bs > 1) pc_ml->nulltype = PCML_NULLSPACE_BLOCK;
678:       else pc_ml->nulltype = PCML_NULLSPACE_SCALAR;
679:     }
680:     switch (pc_ml->nulltype) {
681:     case PCML_NULLSPACE_USER: {
682:       PetscScalar       *nullvec;
683:       const PetscScalar *v;
684:       PetscBool          has_const;
685:       PetscInt           i, j, mlocal, nvec, M;
686:       const Vec         *vecs;

688:       PetscCheck(mnull, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "Must provide explicit null space using MatSetNearNullSpace() to use user-specified null space");
689:       PetscCall(MatGetSize(A, &M, NULL));
690:       PetscCall(MatGetLocalSize(Aloc, &mlocal, NULL));
691:       PetscCall(MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs));
692:       PetscCall(PetscMalloc1((nvec + !!has_const) * mlocal, &nullvec));
693:       if (has_const)
694:         for (i = 0; i < mlocal; i++) nullvec[i] = 1.0 / M;
695:       for (i = 0; i < nvec; i++) {
696:         PetscCall(VecGetArrayRead(vecs[i], &v));
697:         for (j = 0; j < mlocal; j++) nullvec[(i + !!has_const) * mlocal + j] = v[j];
698:         PetscCall(VecRestoreArrayRead(vecs[i], &v));
699:       }
700:       PetscStackCallExternalVoid("ML_Aggregate_Create", PetscCall(ML_Aggregate_Set_NullSpace(agg_object, bs, nvec + !!has_const, nullvec, mlocal)));
701:       PetscCall(PetscFree(nullvec));
702:     } break;
703:     case PCML_NULLSPACE_BLOCK:
704:       PetscStackCallExternalVoid("ML_Aggregate_Set_NullSpace", PetscCall(ML_Aggregate_Set_NullSpace(agg_object, bs, bs, 0, 0)));
705:       break;
706:     case PCML_NULLSPACE_SCALAR:
707:       break;
708:     default:
709:       SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Unknown null space type");
710:     }
711:   }
712:   PetscStackCallExternalVoid("ML_Aggregate_Set_MaxCoarseSize", ML_Aggregate_Set_MaxCoarseSize(agg_object, pc_ml->MaxCoarseSize));
713:   /* set options */
714:   switch (pc_ml->CoarsenScheme) {
715:   case 1:
716:     PetscStackCallExternalVoid("ML_Aggregate_Set_CoarsenScheme_Coupled", ML_Aggregate_Set_CoarsenScheme_Coupled(agg_object));
717:     break;
718:   case 2:
719:     PetscStackCallExternalVoid("ML_Aggregate_Set_CoarsenScheme_MIS", ML_Aggregate_Set_CoarsenScheme_MIS(agg_object));
720:     break;
721:   case 3:
722:     PetscStackCallExternalVoid("ML_Aggregate_Set_CoarsenScheme_METIS", ML_Aggregate_Set_CoarsenScheme_METIS(agg_object));
723:     break;
724:   }
725:   PetscStackCallExternalVoid("ML_Aggregate_Set_Threshold", ML_Aggregate_Set_Threshold(agg_object, pc_ml->Threshold));
726:   PetscStackCallExternalVoid("ML_Aggregate_Set_DampingFactor", ML_Aggregate_Set_DampingFactor(agg_object, pc_ml->DampingFactor));
727:   if (pc_ml->SpectralNormScheme_Anorm) PetscStackCallExternalVoid("ML_Set_SpectralNormScheme_Anorm", ML_Set_SpectralNormScheme_Anorm(ml_object));
728:   agg_object->keep_agg_information      = (int)pc_ml->KeepAggInfo;
729:   agg_object->keep_P_tentative          = (int)pc_ml->Reusable;
730:   agg_object->block_scaled_SA           = (int)pc_ml->BlockScaling;
731:   agg_object->minimizing_energy         = (int)pc_ml->EnergyMinimization;
732:   agg_object->minimizing_energy_droptol = (double)pc_ml->EnergyMinimizationDropTol;
733:   agg_object->cheap_minimizing_energy   = (int)pc_ml->EnergyMinimizationCheap;

735:   if (pc_ml->Aux) {
736:     PetscCheck(pc_ml->dim, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "Auxiliary matrix requires coordinates");
737:     ml_object->Amat[0].aux_data->threshold = pc_ml->AuxThreshold;
738:     ml_object->Amat[0].aux_data->enable    = 1;
739:     ml_object->Amat[0].aux_data->max_level = 10;
740:     ml_object->Amat[0].num_PDEs            = bs;
741:   }

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

746:   if (pc_ml->dim) {
747:     PetscInt                i, dim = pc_ml->dim;
748:     ML_Aggregate_Viz_Stats *grid_info;
749:     PetscInt                nlocghost;

751:     PetscCall(MatGetBlockSize(A, &bs));
752:     nlocghost = Aloc->cmap->n / bs;

754:     PetscStackCallExternalVoid("ML_Aggregate_VizAndStats_Setup(", ML_Aggregate_VizAndStats_Setup(ml_object)); /* create ml info for coords */
755:     grid_info = (ML_Aggregate_Viz_Stats *)ml_object->Grid[0].Grid;
756:     for (i = 0; i < dim; i++) {
757:       /* set the finest level coordinates to point to the column-order array
758:        * in pc_ml */
759:       /* NOTE: must point away before VizAndStats_Clean so ML doesn't free */
760:       switch (i) {
761:       case 0:
762:         grid_info->x = pc_ml->coords + nlocghost * i;
763:         break;
764:       case 1:
765:         grid_info->y = pc_ml->coords + nlocghost * i;
766:         break;
767:       case 2:
768:         grid_info->z = pc_ml->coords + nlocghost * i;
769:         break;
770:       default:
771:         SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_SIZ, "PCML coordinate dimension must be <= 3");
772:       }
773:     }
774:     grid_info->Ndim = dim;
775:   }

777:   /* repartitioning */
778:   if (pc_ml->Repartition) {
779:     PetscStackCallExternalVoid("ML_Repartition_Activate", ML_Repartition_Activate(ml_object));
780:     PetscStackCallExternalVoid("ML_Repartition_Set_LargestMinMaxRatio", ML_Repartition_Set_LargestMinMaxRatio(ml_object, pc_ml->MaxMinRatio));
781:     PetscStackCallExternalVoid("ML_Repartition_Set_MinPerProc", ML_Repartition_Set_MinPerProc(ml_object, pc_ml->MinPerProc));
782:     PetscStackCallExternalVoid("ML_Repartition_Set_PutOnSingleProc", ML_Repartition_Set_PutOnSingleProc(ml_object, pc_ml->PutOnSingleProc));
783: #if 0 /* Function not yet defined in ml-6.2 */
784:     /* I'm not sure what compatibility issues might crop up if we partitioned
785:      * on the finest level, so to be safe repartition starts on the next
786:      * finest level (reflection default behavior in
787:      * ml_MultiLevelPreconditioner) */
788:     PetscStackCallExternalVoid("ML_Repartition_Set_StartLevel",ML_Repartition_Set_StartLevel(ml_object,1));
789: #endif

791:     if (!pc_ml->RepartitionType) {
792:       PetscInt i;

794:       PetscCheck(pc_ml->dim, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "ML Zoltan repartitioning requires coordinates");
795:       PetscStackCallExternalVoid("ML_Repartition_Set_Partitioner", ML_Repartition_Set_Partitioner(ml_object, ML_USEZOLTAN));
796:       PetscStackCallExternalVoid("ML_Aggregate_Set_Dimensions", ML_Aggregate_Set_Dimensions(agg_object, pc_ml->dim));

798:       for (i = 0; i < ml_object->ML_num_levels; i++) {
799:         ML_Aggregate_Viz_Stats *grid_info = (ML_Aggregate_Viz_Stats *)ml_object->Grid[i].Grid;
800:         grid_info->zoltan_type            = pc_ml->ZoltanScheme + 1; /* ml numbers options 1, 2, 3 */
801:         /* defaults from ml_agg_info.c */
802:         grid_info->zoltan_estimated_its = 40; /* only relevant to hypergraph / fast hypergraph */
803:         grid_info->zoltan_timers        = 0;
804:         grid_info->smoothing_steps      = 4; /* only relevant to hypergraph / fast hypergraph */
805:       }
806:     } else {
807:       PetscStackCallExternalVoid("ML_Repartition_Set_Partitioner", ML_Repartition_Set_Partitioner(ml_object, ML_USEPARMETIS));
808:     }
809:   }

811:   if (pc_ml->OldHierarchy) {
812:     PetscStackCallExternalVoid("ML_Gen_MGHierarchy_UsingAggregation", Nlevels = ML_Gen_MGHierarchy_UsingAggregation(ml_object, 0, ML_INCREASING, agg_object));
813:   } else {
814:     PetscStackCallExternalVoid("ML_Gen_MultiLevelHierarchy_UsingAggregation", Nlevels = ML_Gen_MultiLevelHierarchy_UsingAggregation(ml_object, 0, ML_INCREASING, agg_object));
815:   }
816:   PetscCheck(Nlevels > 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Nlevels %d must > 0", Nlevels);
817:   pc_ml->Nlevels = Nlevels;
818:   fine_level     = Nlevels - 1;

820:   PetscCall(PCMGSetLevels(pc, Nlevels, NULL));
821:   /* set default smoothers */
822:   for (level = 1; level <= fine_level; level++) {
823:     PetscCall(PCMGGetSmoother(pc, level, &smoother));
824:     PetscCall(KSPSetType(smoother, KSPRICHARDSON));
825:     PetscCall(KSPGetPC(smoother, &subpc));
826:     PetscCall(PCSetType(subpc, PCSOR));
827:   }
828:   PetscObjectOptionsBegin((PetscObject)pc);
829:   PetscCall(PCSetFromOptions_MG(pc, PetscOptionsObject)); /* should be called in PCSetFromOptions_ML(), but cannot be called prior to PCMGSetLevels() */
830:   PetscOptionsEnd();

832:   PetscCall(PetscMalloc1(Nlevels, &gridctx));

834:   pc_ml->gridctx = gridctx;

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

840:   level = fine_level - 1;
841:   /* TODO: support for GPUs */
842:   if (size == 1) { /* convert ML P, R and A into seqaij format */
843:     for (mllevel = 1; mllevel < Nlevels; mllevel++) {
844:       mlmat = &(ml_object->Pmat[mllevel]);
845:       PetscCall(MatWrapML_SeqAIJ(mlmat, MAT_INITIAL_MATRIX, &gridctx[level].P));
846:       mlmat = &(ml_object->Rmat[mllevel - 1]);
847:       PetscCall(MatWrapML_SeqAIJ(mlmat, MAT_INITIAL_MATRIX, &gridctx[level].R));

849:       mlmat = &(ml_object->Amat[mllevel]);
850:       PetscCall(MatWrapML_SeqAIJ(mlmat, MAT_INITIAL_MATRIX, &gridctx[level].A));
851:       level--;
852:     }
853:   } else { /* convert ML P and R into shell format, ML A into mpiaij format */
854:     for (mllevel = 1; mllevel < Nlevels; mllevel++) {
855:       mlmat = &(ml_object->Pmat[mllevel]);
856:       PetscCall(MatWrapML_SHELL(mlmat, MAT_INITIAL_MATRIX, &gridctx[level].P));
857:       mlmat = &(ml_object->Rmat[mllevel - 1]);
858:       PetscCall(MatWrapML_SHELL(mlmat, MAT_INITIAL_MATRIX, &gridctx[level].R));

860:       mlmat = &(ml_object->Amat[mllevel]);
861:       PetscCall(MatWrapML_MPIAIJ(mlmat, MAT_INITIAL_MATRIX, &gridctx[level].A));
862:       level--;
863:     }
864:   }

866:   /* create vectors and ksp at all levels */
867:   for (level = 0; level < fine_level; level++) {
868:     level1 = level + 1;

870:     PetscCall(MatCreateVecs(gridctx[level].A, &gridctx[level].x, &gridctx[level].b));
871:     PetscCall(MatCreateVecs(gridctx[level1].A, NULL, &gridctx[level1].r));
872:     PetscCall(PCMGSetX(pc, level, gridctx[level].x));
873:     PetscCall(PCMGSetRhs(pc, level, gridctx[level].b));
874:     PetscCall(PCMGSetR(pc, level1, gridctx[level1].r));

876:     if (level == 0) {
877:       PetscCall(PCMGGetCoarseSolve(pc, &gridctx[level].ksp));
878:     } else {
879:       PetscCall(PCMGGetSmoother(pc, level, &gridctx[level].ksp));
880:     }
881:   }
882:   PetscCall(PCMGGetSmoother(pc, fine_level, &gridctx[fine_level].ksp));

884:   /* create coarse level and the interpolation between the levels */
885:   for (level = 0; level < fine_level; level++) {
886:     level1 = level + 1;

888:     PetscCall(PCMGSetInterpolation(pc, level1, gridctx[level].P));
889:     PetscCall(PCMGSetRestriction(pc, level1, gridctx[level].R));
890:     if (level > 0) PetscCall(PCMGSetResidual(pc, level, PCMGResidualDefault, gridctx[level].A));
891:     PetscCall(KSPSetOperators(gridctx[level].ksp, gridctx[level].A, gridctx[level].A));
892:   }
893:   PetscCall(PCMGSetResidual(pc, fine_level, PCMGResidualDefault, gridctx[fine_level].A));
894:   PetscCall(KSPSetOperators(gridctx[fine_level].ksp, gridctx[level].A, gridctx[fine_level].A));

896:   /* put coordinate info in levels */
897:   if (pc_ml->dim) {
898:     PetscInt   i, j, dim = pc_ml->dim;
899:     PetscInt   bs, nloc;
900:     PC         subpc;
901:     PetscReal *array;

903:     level = fine_level;
904:     for (mllevel = 0; mllevel < Nlevels; mllevel++) {
905:       ML_Aggregate_Viz_Stats *grid_info = (ML_Aggregate_Viz_Stats *)ml_object->Amat[mllevel].to->Grid->Grid;
906:       MPI_Comm                comm      = ((PetscObject)gridctx[level].A)->comm;

908:       PetscCall(MatGetBlockSize(gridctx[level].A, &bs));
909:       PetscCall(MatGetLocalSize(gridctx[level].A, NULL, &nloc));
910:       nloc /= bs; /* number of local nodes */

912:       PetscCall(VecCreate(comm, &gridctx[level].coords));
913:       PetscCall(VecSetSizes(gridctx[level].coords, dim * nloc, PETSC_DECIDE));
914:       PetscCall(VecSetType(gridctx[level].coords, VECMPI));
915:       PetscCall(VecGetArray(gridctx[level].coords, &array));
916:       for (j = 0; j < nloc; j++) {
917:         for (i = 0; i < dim; i++) {
918:           switch (i) {
919:           case 0:
920:             array[dim * j + i] = grid_info->x[j];
921:             break;
922:           case 1:
923:             array[dim * j + i] = grid_info->y[j];
924:             break;
925:           case 2:
926:             array[dim * j + i] = grid_info->z[j];
927:             break;
928:           default:
929:             SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_SIZ, "PCML coordinate dimension must be <= 3");
930:           }
931:         }
932:       }

934:       /* passing coordinates to smoothers/coarse solver, should they need them */
935:       PetscCall(KSPGetPC(gridctx[level].ksp, &subpc));
936:       PetscCall(PCSetCoordinates(subpc, dim, nloc, array));
937:       PetscCall(VecRestoreArray(gridctx[level].coords, &array));
938:       level--;
939:     }
940:   }

942:   /* setupcalled is set to 0 so that MG is setup from scratch */
943:   pc->setupcalled = 0;
944:   PetscCall(PCSetUp_MG(pc));
945:   PetscFunctionReturn(PETSC_SUCCESS);
946: }

948: /* -------------------------------------------------------------------------- */
949: /*
950:    PCDestroy_ML - Destroys the private context for the ML preconditioner
951:    that was created with PCCreate_ML().

953:    Input Parameter:
954: .  pc - the preconditioner context

956:    Application Interface Routine: PCDestroy()
957: */
958: PetscErrorCode PCDestroy_ML(PC pc)
959: {
960:   PC_MG *mg    = (PC_MG *)pc->data;
961:   PC_ML *pc_ml = (PC_ML *)mg->innerctx;

963:   PetscFunctionBegin;
964:   PetscCall(PCReset_ML(pc));
965:   PetscCall(PetscFree(pc_ml));
966:   PetscCall(PCDestroy_MG(pc));
967:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
968:   PetscFunctionReturn(PETSC_SUCCESS);
969: }

971: PetscErrorCode PCSetFromOptions_ML(PC pc, PetscOptionItems *PetscOptionsObject)
972: {
973:   PetscInt    indx, PrintLevel, partindx;
974:   const char *scheme[] = {"Uncoupled", "Coupled", "MIS", "METIS"};
975:   const char *part[]   = {"Zoltan", "ParMETIS"};
976: #if defined(HAVE_ML_ZOLTAN)
977:   const char *zscheme[] = {"RCB", "hypergraph", "fast_hypergraph"};
978: #endif
979:   PC_MG      *mg    = (PC_MG *)pc->data;
980:   PC_ML      *pc_ml = (PC_ML *)mg->innerctx;
981:   PetscMPIInt size;
982:   MPI_Comm    comm;

984:   PetscFunctionBegin;
985:   PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
986:   PetscCallMPI(MPI_Comm_size(comm, &size));
987:   PetscOptionsHeadBegin(PetscOptionsObject, "ML options");

989:   PrintLevel = 0;
990:   indx       = 0;
991:   partindx   = 0;

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

999:   pc_ml->CoarsenScheme = indx;

1001:   PetscCall(PetscOptionsReal("-pc_ml_DampingFactor", "P damping factor", "ML_Aggregate_Set_DampingFactor", pc_ml->DampingFactor, &pc_ml->DampingFactor, NULL));
1002:   PetscCall(PetscOptionsReal("-pc_ml_Threshold", "Smoother drop tol", "ML_Aggregate_Set_Threshold", pc_ml->Threshold, &pc_ml->Threshold, NULL));
1003:   PetscCall(PetscOptionsBool("-pc_ml_SpectralNormScheme_Anorm", "Method used for estimating spectral radius", "ML_Set_SpectralNormScheme_Anorm", pc_ml->SpectralNormScheme_Anorm, &pc_ml->SpectralNormScheme_Anorm, NULL));
1004:   PetscCall(PetscOptionsBool("-pc_ml_Symmetrize", "Symmetrize aggregation", "ML_Set_Symmetrize", pc_ml->Symmetrize, &pc_ml->Symmetrize, NULL));
1005:   PetscCall(PetscOptionsBool("-pc_ml_BlockScaling", "Scale all dofs at each node together", "None", pc_ml->BlockScaling, &pc_ml->BlockScaling, NULL));
1006:   PetscCall(PetscOptionsEnum("-pc_ml_nullspace", "Which type of null space information to use", "None", PCMLNullSpaceTypes, (PetscEnum)pc_ml->nulltype, (PetscEnum *)&pc_ml->nulltype, NULL));
1007:   PetscCall(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));
1008:   PetscCall(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));
1009:   /*
1010:     The following checks a number of conditions.  If we let this stuff slip by, then ML's error handling will take over.
1011:     This is suboptimal because it amounts to calling exit(1) so we check for the most common conditions.

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

1047:     pc_ml->RepartitionType = partindx;
1048:     if (!partindx) {
1049:       PetscInt zindx = 0;

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

1053:       pc_ml->ZoltanScheme = zindx;
1054:     }
1055: #else
1056:     partindx = 1;
1057:     PetscCall(PetscOptionsEList("-pc_ml_repartitionType", "Repartitioning library to use", "ML_Repartition_Set_Partitioner", part, 2, part[1], &partindx, NULL));
1058:     pc_ml->RepartitionType = partindx;
1059:     PetscCheck(partindx, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP_SYS, "ML not compiled with Zoltan");
1060: #endif
1061:     PetscCall(PetscOptionsBool("-pc_ml_Aux", "Aggregate using auxiliary coordinate-based laplacian", "None", pc_ml->Aux, &pc_ml->Aux, NULL));
1062:     PetscCall(PetscOptionsReal("-pc_ml_AuxThreshold", "Auxiliary smoother drop tol", "None", pc_ml->AuxThreshold, &pc_ml->AuxThreshold, NULL));
1063:   }
1064:   PetscOptionsHeadEnd();
1065:   PetscFunctionReturn(PETSC_SUCCESS);
1066: }

1068: /*
1069:    PCCreate_ML - Creates a ML preconditioner context, PC_ML,
1070:    and sets this as the private data within the generic preconditioning
1071:    context, PC, that was created within PCCreate().

1073:    Input Parameter:
1074: .  pc - the preconditioner context

1076:    Application Interface Routine: PCCreate()
1077: */

1079: /*MC
1080:      PCML - Use the SNL ML algebraic multigrid preconditioner.

1082:    Options Database Keys:
1083:    Multigrid options(inherited):
1084: +  -pc_mg_cycle_type <v> - v for V cycle, w for W-cycle (`PCMGSetCycleType()`)
1085: .  -pc_mg_distinct_smoothup - Should one configure the up and down smoothers separately (`PCMGSetDistinctSmoothUp()`)
1086: -  -pc_mg_type <multiplicative> - (one of) additive multiplicative full kascade

1088:    ML Options Database Key:
1089: +  -pc_ml_PrintLevel <0> - Print level (`ML_Set_PrintLevel()`)
1090: .  -pc_ml_maxNlevels <10> - Maximum number of levels (None)
1091: .  -pc_ml_maxCoarseSize <1> - Maximum coarsest mesh size (`ML_Aggregate_Set_MaxCoarseSize()`)
1092: .  -pc_ml_CoarsenScheme <Uncoupled> - (one of) Uncoupled Coupled MIS METIS
1093: .  -pc_ml_DampingFactor <1.33333> - P damping factor (`ML_Aggregate_Set_DampingFactor()`)
1094: .  -pc_ml_Threshold <0> - Smoother drop tol (`ML_Aggregate_Set_Threshold()`)
1095: .  -pc_ml_SpectralNormScheme_Anorm <false> - Method used for estimating spectral radius (`ML_Set_SpectralNormScheme_Anorm()`)
1096: .  -pc_ml_repartition <false> - Allow ML to repartition levels of the hierarchy (`ML_Repartition_Activate()`)
1097: .  -pc_ml_repartitionMaxMinRatio <1.3> - Acceptable ratio of repartitioned sizes (`ML_Repartition_Set_LargestMinMaxRatio()`)
1098: .  -pc_ml_repartitionMinPerProc <512> - Smallest repartitioned size (`ML_Repartition_Set_MinPerProc()`)
1099: .  -pc_ml_repartitionPutOnSingleProc <5000> - Problem size automatically repartitioned to one processor (`ML_Repartition_Set_PutOnSingleProc()`)
1100: .  -pc_ml_repartitionType <Zoltan> - Repartitioning library to use (`ML_Repartition_Set_Partitioner()`)
1101: .  -pc_ml_repartitionZoltanScheme <RCB> - Repartitioning scheme to use (None)
1102: .  -pc_ml_Aux <false> - Aggregate using auxiliary coordinate-based Laplacian (None)
1103: -  -pc_ml_AuxThreshold <0.0> - Auxiliary smoother drop tol (None)

1105:    Level: intermediate

1107:    Developer Note:
1108:    The coarser grid matrices and restriction/interpolation
1109:    operators are computed by ML, with the matrices converted to PETSc matrices in `MATAIJ` format
1110:    and the restriction/interpolation operators wrapped as PETSc shell matrices.

1112: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCMGType`, `PCMG`, `PCHYPRE`, `PCGAMG`,
1113:           `PCMGSetLevels()`, `PCMGGetLevels()`, `PCMGSetType()`, `MPSetCycles()`, `PCMGSetDistinctSmoothUp()`,
1114:           `PCMGGetCoarseSolve()`, `PCMGSetResidual()`, `PCMGSetInterpolation()`,
1115:           `PCMGSetRestriction()`, `PCMGGetSmoother()`, `PCMGGetSmootherUp()`, `PCMGGetSmootherDown()`,
1116:           `PCMGSetCycleTypeOnLevel()`, `PCMGSetRhs()`, `PCMGSetX()`, `PCMGSetR()`
1117: M*/

1119: PETSC_EXTERN PetscErrorCode PCCreate_ML(PC pc)
1120: {
1121:   PC_ML *pc_ml;
1122:   PC_MG *mg;

1124:   PetscFunctionBegin;
1125:   /* PCML is an inherited class of PCMG. Initialize pc as PCMG */
1126:   PetscCall(PCSetType(pc, PCMG)); /* calls PCCreate_MG() and MGCreate_Private() */
1127:   PetscCall(PetscObjectChangeTypeName((PetscObject)pc, PCML));
1128:   /* Since PCMG tries to use DM associated with PC must delete it */
1129:   PetscCall(DMDestroy(&pc->dm));
1130:   PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_EXTERNAL));
1131:   mg = (PC_MG *)pc->data;

1133:   /* create a supporting struct and attach it to pc */
1134:   PetscCall(PetscNew(&pc_ml));
1135:   mg->innerctx = pc_ml;

1137:   pc_ml->ml_object                = 0;
1138:   pc_ml->agg_object               = 0;
1139:   pc_ml->gridctx                  = 0;
1140:   pc_ml->PetscMLdata              = 0;
1141:   pc_ml->Nlevels                  = -1;
1142:   pc_ml->MaxNlevels               = 10;
1143:   pc_ml->MaxCoarseSize            = 1;
1144:   pc_ml->CoarsenScheme            = 1;
1145:   pc_ml->Threshold                = 0.0;
1146:   pc_ml->DampingFactor            = 4.0 / 3.0;
1147:   pc_ml->SpectralNormScheme_Anorm = PETSC_FALSE;
1148:   pc_ml->size                     = 0;
1149:   pc_ml->dim                      = 0;
1150:   pc_ml->nloc                     = 0;
1151:   pc_ml->coords                   = 0;
1152:   pc_ml->Repartition              = PETSC_FALSE;
1153:   pc_ml->MaxMinRatio              = 1.3;
1154:   pc_ml->MinPerProc               = 512;
1155:   pc_ml->PutOnSingleProc          = 5000;
1156:   pc_ml->RepartitionType          = 0;
1157:   pc_ml->ZoltanScheme             = 0;
1158:   pc_ml->Aux                      = PETSC_FALSE;
1159:   pc_ml->AuxThreshold             = 0.0;

1161:   /* allow for coordinates to be passed */
1162:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_ML));

1164:   /* overwrite the pointers of PCMG by the functions of PCML */
1165:   pc->ops->setfromoptions = PCSetFromOptions_ML;
1166:   pc->ops->setup          = PCSetUp_ML;
1167:   pc->ops->reset          = PCReset_ML;
1168:   pc->ops->destroy        = PCDestroy_ML;
1169:   PetscFunctionReturn(PETSC_SUCCESS);
1170: }