Actual source code: ml.c

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

     12: EXTERN_C_BEGIN
     13: /* HAVE_CONFIG_H flag is required by ML include files */
     14: #if !defined(HAVE_CONFIG_H)
     15:   #define HAVE_CONFIG_H
     16: #endif
     17: #include <ml_include.h>
     18: #include <ml_viz_stats.h>
     19: EXTERN_C_END

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

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

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

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

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

     69: 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[])
     70: {
     71:   PetscInt     m, i, j, k = 0, row, *aj;
     72:   PetscScalar *aa;
     73:   FineGridCtx *ml = (FineGridCtx *)ML_Get_MyGetrowData(ML_data);
     74:   Mat_SeqAIJ  *a  = (Mat_SeqAIJ *)ml->Aloc->data;

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

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

    102:   PetscFunctionBegin;
    103:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
    104:   if (size == 1) PetscFunctionReturn(PETSC_SUCCESS);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    324:   MLcomm = mlmat->comm;

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

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

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

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

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

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

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

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

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

    397: /*
    398:    PCSetCoordinates_ML

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

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

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

    419:   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);

    421:   oldarrsz    = pc_ml->dim * pc_ml->nloc;
    422:   pc_ml->dim  = ndm;
    423:   pc_ml->nloc = nloc;
    424:   arrsz       = ndm * nloc;

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

    445: static PetscErrorCode PCReset_ML(PC pc)
    446: {
    447:   PC_MG   *mg    = (PC_MG *)pc->data;
    448:   PC_ML   *pc_ml = (PC_ML *)mg->innerctx;
    449:   PetscInt level, fine_level = pc_ml->Nlevels - 1, dim = pc_ml->dim;

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

    465:   if (pc_ml->PetscMLdata) {
    466:     PetscCall(PetscFree(pc_ml->PetscMLdata->pwork));
    467:     PetscCall(MatDestroy(&pc_ml->PetscMLdata->Aloc));
    468:     PetscCall(VecDestroy(&pc_ml->PetscMLdata->x));
    469:     PetscCall(VecDestroy(&pc_ml->PetscMLdata->y));
    470:   }
    471:   PetscCall(PetscFree(pc_ml->PetscMLdata));

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

    486:   pc_ml->dim  = 0;
    487:   pc_ml->nloc = 0;
    488:   PetscCall(PCReset_MG(pc));
    489:   PetscFunctionReturn(PETSC_SUCCESS);
    490: }

    492: /*
    493:    PCSetUp_ML - Prepares for the use of the ML preconditioner
    494:                     by setting data structures and options.

    496:    Input Parameter:
    497: .  pc - the preconditioner context

    499:    Application Interface Routine: PCSetUp()

    501:    Note:
    502:    The interface routine PCSetUp() is not usually called directly by
    503:    the user, but instead is called by PCApply() if necessary.
    504: */
    505: static PetscErrorCode PCSetUp_ML(PC pc)
    506: {
    507:   PetscMPIInt      size;
    508:   FineGridCtx     *PetscMLdata;
    509:   ML              *ml_object;
    510:   ML_Aggregate    *agg_object;
    511:   ML_Operator     *mlmat;
    512:   PetscInt         nlocal_allcols, Nlevels, mllevel, level, level1, m, fine_level, bs;
    513:   Mat              A, Aloc;
    514:   GridCtx         *gridctx;
    515:   PC_MG           *mg    = (PC_MG *)pc->data;
    516:   PC_ML           *pc_ml = (PC_ML *)mg->innerctx;
    517:   PetscBool        isSeq, isMPI;
    518:   KSP              smoother;
    519:   PC               subpc;
    520:   PetscInt         mesh_level, old_mesh_level;
    521:   MatInfo          info;
    522:   static PetscBool cite = PETSC_FALSE;

    524:   PetscFunctionBegin;
    525:   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 = "
    526:                                    "{SAND2004-4821},\n  year = 2004\n}\n",
    527:                                    &cite));
    528:   A = pc->pmat;
    529:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));

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

    543:       PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &isSeq));
    544:       PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &isMPI));
    545:       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);
    546:       if (isMPI) {
    547:         PetscCall(MatConvert_MPIAIJ_ML(A, NULL, MAT_INITIAL_MATRIX, &Aloc));
    548:       } else {
    549:         Aloc = A;
    550:         PetscCall(PetscObjectReference((PetscObject)Aloc));
    551:       }

    553:       PetscCall(MatGetSize(Aloc, &m, &nlocal_allcols));
    554:       PetscMLdata = pc_ml->PetscMLdata;
    555:       PetscCall(MatDestroy(&PetscMLdata->Aloc));
    556:       PetscMLdata->A    = A;
    557:       PetscMLdata->Aloc = Aloc;
    558:       PetscCallExternalVoid("ML_Aggregate_Destroy", ML_Init_Amatrix(ml_object, 0, m, m, PetscMLdata));
    559:       PetscCallExternalVoid("ML_Set_Amatrix_Matvec", ML_Set_Amatrix_Matvec(ml_object, 0, PetscML_matvec));

    561:       mesh_level = ml_object->ML_finest_level;
    562:       while (ml_object->SingleLevel[mesh_level].Rmat->to) {
    563:         old_mesh_level = mesh_level;
    564:         mesh_level     = ml_object->SingleLevel[mesh_level].Rmat->to->levelnum;

    566:         /* clean and regenerate A */
    567:         mlmat = &ml_object->Amat[mesh_level];
    568:         PetscCallExternalVoid("ML_Operator_Clean", ML_Operator_Clean(mlmat));
    569:         PetscCallExternalVoid("ML_Operator_Init", ML_Operator_Init(mlmat, ml_object->comm));
    570:         PetscCallExternalVoid("ML_Gen_AmatrixRAP", ML_Gen_AmatrixRAP(ml_object, old_mesh_level, mesh_level));
    571:       }

    573:       level = fine_level - 1;
    574:       if (size == 1) { /* convert ML P, R and A into seqaij format */
    575:         for (mllevel = 1; mllevel < Nlevels; mllevel++) {
    576:           mlmat = &ml_object->Amat[mllevel];
    577:           PetscCall(MatWrapML_SeqAIJ(mlmat, MAT_REUSE_MATRIX, &gridctx[level].A));
    578:           level--;
    579:         }
    580:       } else { /* convert ML P and R into shell format, ML A into mpiaij format */
    581:         for (mllevel = 1; mllevel < Nlevels; mllevel++) {
    582:           mlmat = &ml_object->Amat[mllevel];
    583:           PetscCall(MatWrapML_MPIAIJ(mlmat, MAT_REUSE_MATRIX, &gridctx[level].A));
    584:           level--;
    585:         }
    586:       }

    588:       for (level = 0; level < fine_level; level++) {
    589:         if (level > 0) PetscCall(PCMGSetResidual(pc, level, PCMGResidualDefault, gridctx[level].A));
    590:         PetscCall(KSPSetOperators(gridctx[level].ksp, gridctx[level].A, gridctx[level].A));
    591:       }
    592:       PetscCall(PCMGSetResidual(pc, fine_level, PCMGResidualDefault, gridctx[fine_level].A));
    593:       PetscCall(KSPSetOperators(gridctx[fine_level].ksp, gridctx[level].A, gridctx[fine_level].A));

    595:       PetscCall(PCSetUp_MG(pc));
    596:       PetscFunctionReturn(PETSC_SUCCESS);
    597:     } else {
    598:       /* since ML can change the size of vectors/matrices at any level we must destroy everything */
    599:       PetscCall(PCReset_ML(pc));
    600:     }
    601:   }

    603:   /* setup special features of PCML */
    604:   /* convert A to Aloc to be used by ML at fine grid */
    605:   pc_ml->size = size;
    606:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &isSeq));
    607:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &isMPI));
    608:   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);
    609:   if (isMPI) {
    610:     PetscCall(MatConvert_MPIAIJ_ML(A, NULL, MAT_INITIAL_MATRIX, &Aloc));
    611:   } else {
    612:     Aloc = A;
    613:     PetscCall(PetscObjectReference((PetscObject)Aloc));
    614:   }

    616:   /* create and initialize struct 'PetscMLdata' */
    617:   PetscCall(PetscNew(&PetscMLdata));
    618:   pc_ml->PetscMLdata = PetscMLdata;
    619:   PetscCall(PetscMalloc1(Aloc->cmap->n + 1, &PetscMLdata->pwork));

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

    623:   PetscMLdata->A    = A;
    624:   PetscMLdata->Aloc = Aloc;
    625:   if (pc_ml->dim) { /* create vecs around the coordinate data given */
    626:     PetscInt   i, j, dim = pc_ml->dim;
    627:     PetscInt   nloc = pc_ml->nloc, nlocghost;
    628:     PetscReal *ghostedcoords;

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

    647:   /* create ML discretization matrix at fine grid */
    648:   /* ML requires input of fine-grid matrix. It determines nlevels. */
    649:   PetscCall(MatGetSize(Aloc, &m, &nlocal_allcols));
    650:   PetscCall(MatGetBlockSize(A, &bs));
    651:   PetscCallExternalVoid("ML_Create", ML_Create(&ml_object, pc_ml->MaxNlevels));
    652:   PetscCallExternalVoid("ML_Comm_Set_UsrComm", ML_Comm_Set_UsrComm(ml_object->comm, PetscObjectComm((PetscObject)A)));
    653:   pc_ml->ml_object = ml_object;
    654:   PetscCallExternalVoid("ML_Init_Amatrix", ML_Init_Amatrix(ml_object, 0, m, m, PetscMLdata));
    655:   PetscCallExternalVoid("ML_Set_Amatrix_Getrow", ML_Set_Amatrix_Getrow(ml_object, 0, PetscML_getrow, ML_PetscML_comm, nlocal_allcols));
    656:   PetscCallExternalVoid("ML_Set_Amatrix_Matvec", ML_Set_Amatrix_Matvec(ml_object, 0, PetscML_matvec));

    658:   PetscCallExternalVoid("ML_Set_Symmetrize", ML_Set_Symmetrize(ml_object, pc_ml->Symmetrize ? ML_YES : ML_NO));

    660:   /* aggregation */
    661:   PetscCallExternalVoid("ML_Aggregate_Create", ML_Aggregate_Create(&agg_object));
    662:   pc_ml->agg_object = agg_object;

    664:   {
    665:     MatNullSpace mnull;
    666:     PetscCall(MatGetNearNullSpace(A, &mnull));
    667:     if (pc_ml->nulltype == PCML_NULLSPACE_AUTO) {
    668:       if (mnull) pc_ml->nulltype = PCML_NULLSPACE_USER;
    669:       else if (bs > 1) pc_ml->nulltype = PCML_NULLSPACE_BLOCK;
    670:       else pc_ml->nulltype = PCML_NULLSPACE_SCALAR;
    671:     }
    672:     switch (pc_ml->nulltype) {
    673:     case PCML_NULLSPACE_USER: {
    674:       PetscScalar       *nullvec;
    675:       const PetscScalar *v;
    676:       PetscBool          has_const;
    677:       PetscInt           i, j, mlocal, nvec, M;
    678:       const Vec         *vecs;

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

    727:   if (pc_ml->Aux) {
    728:     PetscCheck(pc_ml->dim, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "Auxiliary matrix requires coordinates");
    729:     ml_object->Amat[0].aux_data->threshold = pc_ml->AuxThreshold;
    730:     ml_object->Amat[0].aux_data->enable    = 1;
    731:     ml_object->Amat[0].aux_data->max_level = 10;
    732:     ml_object->Amat[0].num_PDEs            = bs;
    733:   }

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

    738:   if (pc_ml->dim) {
    739:     PetscInt                i, dim = pc_ml->dim;
    740:     ML_Aggregate_Viz_Stats *grid_info;
    741:     PetscInt                nlocghost;

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

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

    769:   /* repartitioning */
    770:   if (pc_ml->Repartition) {
    771:     PetscCallExternalVoid("ML_Repartition_Activate", ML_Repartition_Activate(ml_object));
    772:     PetscCallExternalVoid("ML_Repartition_Set_LargestMinMaxRatio", ML_Repartition_Set_LargestMinMaxRatio(ml_object, pc_ml->MaxMinRatio));
    773:     PetscCallExternalVoid("ML_Repartition_Set_MinPerProc", ML_Repartition_Set_MinPerProc(ml_object, pc_ml->MinPerProc));
    774:     PetscCallExternalVoid("ML_Repartition_Set_PutOnSingleProc", ML_Repartition_Set_PutOnSingleProc(ml_object, pc_ml->PutOnSingleProc));
    775: #if 0 /* Function not yet defined in ml-6.2 */
    776:     /* I'm not sure what compatibility issues might crop up if we partitioned
    777:      * on the finest level, so to be safe repartition starts on the next
    778:      * finest level (reflection default behavior in
    779:      * ml_MultiLevelPreconditioner) */
    780:     PetscCallExternalVoid("ML_Repartition_Set_StartLevel",ML_Repartition_Set_StartLevel(ml_object,1));
    781: #endif

    783:     if (!pc_ml->RepartitionType) {
    784:       PetscInt i;

    786:       PetscCheck(pc_ml->dim, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "ML Zoltan repartitioning requires coordinates");
    787:       PetscCallExternalVoid("ML_Repartition_Set_Partitioner", ML_Repartition_Set_Partitioner(ml_object, ML_USEZOLTAN));
    788:       PetscCallExternalVoid("ML_Aggregate_Set_Dimensions", ML_Aggregate_Set_Dimensions(agg_object, pc_ml->dim));

    790:       for (i = 0; i < ml_object->ML_num_levels; i++) {
    791:         ML_Aggregate_Viz_Stats *grid_info = (ML_Aggregate_Viz_Stats *)ml_object->Grid[i].Grid;
    792:         grid_info->zoltan_type            = pc_ml->ZoltanScheme + 1; /* ml numbers options 1, 2, 3 */
    793:         /* defaults from ml_agg_info.c */
    794:         grid_info->zoltan_estimated_its = 40; /* only relevant to hypergraph / fast hypergraph */
    795:         grid_info->zoltan_timers        = 0;
    796:         grid_info->smoothing_steps      = 4; /* only relevant to hypergraph / fast hypergraph */
    797:       }
    798:     } else {
    799:       PetscCallExternalVoid("ML_Repartition_Set_Partitioner", ML_Repartition_Set_Partitioner(ml_object, ML_USEPARMETIS));
    800:     }
    801:   }

    803:   if (pc_ml->OldHierarchy) {
    804:     PetscCallExternalVoid("ML_Gen_MGHierarchy_UsingAggregation", Nlevels = ML_Gen_MGHierarchy_UsingAggregation(ml_object, 0, ML_INCREASING, agg_object));
    805:   } else {
    806:     PetscCallExternalVoid("ML_Gen_MultiLevelHierarchy_UsingAggregation", Nlevels = ML_Gen_MultiLevelHierarchy_UsingAggregation(ml_object, 0, ML_INCREASING, agg_object));
    807:   }
    808:   PetscCheck(Nlevels > 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Nlevels %d must > 0", Nlevels);
    809:   pc_ml->Nlevels = Nlevels;
    810:   fine_level     = Nlevels - 1;

    812:   PetscCall(PCMGSetLevels(pc, Nlevels, NULL));
    813:   /* set default smoothers */
    814:   for (level = 1; level <= fine_level; level++) {
    815:     PetscCall(PCMGGetSmoother(pc, level, &smoother));
    816:     PetscCall(KSPSetType(smoother, KSPRICHARDSON));
    817:     PetscCall(KSPGetPC(smoother, &subpc));
    818:     PetscCall(PCSetType(subpc, PCSOR));
    819:   }
    820:   PetscObjectOptionsBegin((PetscObject)pc);
    821:   PetscCall(PCSetFromOptions_MG(pc, PetscOptionsObject)); /* should be called in PCSetFromOptions_ML(), but cannot be called prior to PCMGSetLevels() */
    822:   PetscOptionsEnd();

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

    826:   pc_ml->gridctx = gridctx;

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

    832:   level = fine_level - 1;
    833:   /* TODO: support for GPUs */
    834:   if (size == 1) { /* convert ML P, R and A into seqaij format */
    835:     for (mllevel = 1; mllevel < Nlevels; mllevel++) {
    836:       mlmat = &ml_object->Pmat[mllevel];
    837:       PetscCall(MatWrapML_SeqAIJ(mlmat, MAT_INITIAL_MATRIX, &gridctx[level].P));
    838:       mlmat = &ml_object->Rmat[mllevel - 1];
    839:       PetscCall(MatWrapML_SeqAIJ(mlmat, MAT_INITIAL_MATRIX, &gridctx[level].R));

    841:       mlmat = &ml_object->Amat[mllevel];
    842:       PetscCall(MatWrapML_SeqAIJ(mlmat, MAT_INITIAL_MATRIX, &gridctx[level].A));
    843:       level--;
    844:     }
    845:   } else { /* convert ML P and R into shell format, ML A into mpiaij format */
    846:     for (mllevel = 1; mllevel < Nlevels; mllevel++) {
    847:       mlmat = &ml_object->Pmat[mllevel];
    848:       PetscCall(MatWrapML_SHELL(mlmat, MAT_INITIAL_MATRIX, &gridctx[level].P));
    849:       mlmat = &ml_object->Rmat[mllevel - 1];
    850:       PetscCall(MatWrapML_SHELL(mlmat, MAT_INITIAL_MATRIX, &gridctx[level].R));

    852:       mlmat = &ml_object->Amat[mllevel];
    853:       PetscCall(MatWrapML_MPIAIJ(mlmat, MAT_INITIAL_MATRIX, &gridctx[level].A));
    854:       level--;
    855:     }
    856:   }

    858:   /* create vectors and ksp at all levels */
    859:   for (level = 0; level < fine_level; level++) {
    860:     level1 = level + 1;

    862:     PetscCall(MatCreateVecs(gridctx[level].A, &gridctx[level].x, &gridctx[level].b));
    863:     PetscCall(MatCreateVecs(gridctx[level1].A, NULL, &gridctx[level1].r));
    864:     PetscCall(PCMGSetX(pc, level, gridctx[level].x));
    865:     PetscCall(PCMGSetRhs(pc, level, gridctx[level].b));
    866:     PetscCall(PCMGSetR(pc, level1, gridctx[level1].r));

    868:     if (level == 0) {
    869:       PetscCall(PCMGGetCoarseSolve(pc, &gridctx[level].ksp));
    870:     } else {
    871:       PetscCall(PCMGGetSmoother(pc, level, &gridctx[level].ksp));
    872:     }
    873:   }
    874:   PetscCall(PCMGGetSmoother(pc, fine_level, &gridctx[fine_level].ksp));

    876:   /* create coarse level and the interpolation between the levels */
    877:   for (level = 0; level < fine_level; level++) {
    878:     level1 = level + 1;

    880:     PetscCall(PCMGSetInterpolation(pc, level1, gridctx[level].P));
    881:     PetscCall(PCMGSetRestriction(pc, level1, gridctx[level].R));
    882:     if (level > 0) PetscCall(PCMGSetResidual(pc, level, PCMGResidualDefault, gridctx[level].A));
    883:     PetscCall(KSPSetOperators(gridctx[level].ksp, gridctx[level].A, gridctx[level].A));
    884:   }
    885:   PetscCall(PCMGSetResidual(pc, fine_level, PCMGResidualDefault, gridctx[fine_level].A));
    886:   PetscCall(KSPSetOperators(gridctx[fine_level].ksp, gridctx[level].A, gridctx[fine_level].A));

    888:   /* put coordinate info in levels */
    889:   if (pc_ml->dim) {
    890:     PetscInt   i, j, dim = pc_ml->dim;
    891:     PetscInt   bs, nloc;
    892:     PC         subpc;
    893:     PetscReal *array;

    895:     level = fine_level;
    896:     for (mllevel = 0; mllevel < Nlevels; mllevel++) {
    897:       ML_Aggregate_Viz_Stats *grid_info = (ML_Aggregate_Viz_Stats *)ml_object->Amat[mllevel].to->Grid->Grid;
    898:       MPI_Comm                comm      = ((PetscObject)gridctx[level].A)->comm;

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

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

    926:       /* passing coordinates to smoothers/coarse solver, should they need them */
    927:       PetscCall(KSPGetPC(gridctx[level].ksp, &subpc));
    928:       PetscCall(PCSetCoordinates(subpc, dim, nloc, array));
    929:       PetscCall(VecRestoreArray(gridctx[level].coords, &array));
    930:       level--;
    931:     }
    932:   }

    934:   /* setupcalled is set to 0 so that MG is setup from scratch */
    935:   pc->setupcalled = PETSC_FALSE;
    936:   PetscCall(PCSetUp_MG(pc));
    937:   PetscFunctionReturn(PETSC_SUCCESS);
    938: }

    940: /*
    941:    PCDestroy_ML - Destroys the private context for the ML preconditioner
    942:    that was created with PCCreate_ML().

    944:    Input Parameter:
    945: .  pc - the preconditioner context

    947:    Application Interface Routine: PCDestroy()
    948: */
    949: static PetscErrorCode PCDestroy_ML(PC pc)
    950: {
    951:   PC_MG *mg    = (PC_MG *)pc->data;
    952:   PC_ML *pc_ml = (PC_ML *)mg->innerctx;

    954:   PetscFunctionBegin;
    955:   PetscCall(PCReset_ML(pc));
    956:   PetscCall(PetscFree(pc_ml));
    957:   PetscCall(PCDestroy_MG(pc));
    958:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
    959:   PetscFunctionReturn(PETSC_SUCCESS);
    960: }

    962: static PetscErrorCode PCSetFromOptions_ML(PC pc, PetscOptionItems PetscOptionsObject)
    963: {
    964:   PetscInt    indx, PrintLevel, partindx;
    965:   const char *scheme[] = {"Uncoupled", "Coupled", "MIS", "METIS"};
    966:   const char *part[]   = {"Zoltan", "ParMETIS"};
    967: #if defined(HAVE_ML_ZOLTAN)
    968:   const char *zscheme[] = {"RCB", "hypergraph", "fast_hypergraph"};
    969: #endif
    970:   PC_MG      *mg    = (PC_MG *)pc->data;
    971:   PC_ML      *pc_ml = (PC_ML *)mg->innerctx;
    972:   PetscMPIInt size;
    973:   MPI_Comm    comm;

    975:   PetscFunctionBegin;
    976:   PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
    977:   PetscCallMPI(MPI_Comm_size(comm, &size));
    978:   PetscOptionsHeadBegin(PetscOptionsObject, "ML options");

    980:   PrintLevel = 0;
    981:   indx       = 0;
    982:   partindx   = 0;

    984:   PetscCall(PetscOptionsInt("-pc_ml_PrintLevel", "Print level", "ML_Set_PrintLevel", PrintLevel, &PrintLevel, NULL));
    985:   PetscCallExternalVoid("ML_Set_PrintLevel", ML_Set_PrintLevel(PrintLevel));
    986:   PetscCall(PetscOptionsInt("-pc_ml_maxNlevels", "Maximum number of levels", "None", pc_ml->MaxNlevels, &pc_ml->MaxNlevels, NULL));
    987:   PetscCall(PetscOptionsInt("-pc_ml_maxCoarseSize", "Maximum coarsest mesh size", "ML_Aggregate_Set_MaxCoarseSize", pc_ml->MaxCoarseSize, &pc_ml->MaxCoarseSize, NULL));
    988:   PetscCall(PetscOptionsEList("-pc_ml_CoarsenScheme", "Aggregate Coarsen Scheme", "ML_Aggregate_Set_CoarsenScheme_*", scheme, 4, scheme[0], &indx, NULL));

    990:   pc_ml->CoarsenScheme = indx;

    992:   PetscCall(PetscOptionsReal("-pc_ml_DampingFactor", "P damping factor", "ML_Aggregate_Set_DampingFactor", pc_ml->DampingFactor, &pc_ml->DampingFactor, NULL));
    993:   PetscCall(PetscOptionsReal("-pc_ml_Threshold", "Smoother drop tol", "ML_Aggregate_Set_Threshold", pc_ml->Threshold, &pc_ml->Threshold, NULL));
    994:   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));
    995:   PetscCall(PetscOptionsBool("-pc_ml_Symmetrize", "Symmetrize aggregation", "ML_Set_Symmetrize", pc_ml->Symmetrize, &pc_ml->Symmetrize, NULL));
    996:   PetscCall(PetscOptionsBool("-pc_ml_BlockScaling", "Scale all dofs at each node together", "None", pc_ml->BlockScaling, &pc_ml->BlockScaling, NULL));
    997:   PetscCall(PetscOptionsEnum("-pc_ml_nullspace", "Which type of null space information to use", "None", PCMLNullSpaceTypes, (PetscEnum)pc_ml->nulltype, (PetscEnum *)&pc_ml->nulltype, NULL));
    998:   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));
    999:   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));
    1000:   /*
    1001:     The following checks a number of conditions.  If we let this stuff slip by, then ML's error handling will take over.
    1002:     This is suboptimal because it amounts to calling exit(1) so we check for the most common conditions.

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

    1038:     pc_ml->RepartitionType = partindx;
    1039:     if (!partindx) {
    1040:       PetscInt zindx = 0;

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

    1044:       pc_ml->ZoltanScheme = zindx;
    1045:     }
    1046: #else
    1047:     partindx = 1;
    1048:     PetscCall(PetscOptionsEList("-pc_ml_repartitionType", "Repartitioning library to use", "ML_Repartition_Set_Partitioner", part, 2, part[1], &partindx, NULL));
    1049:     pc_ml->RepartitionType = partindx;
    1050:     PetscCheck(partindx, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP_SYS, "ML not compiled with Zoltan");
    1051: #endif
    1052:     PetscCall(PetscOptionsBool("-pc_ml_Aux", "Aggregate using auxiliary coordinate-based laplacian", "None", pc_ml->Aux, &pc_ml->Aux, NULL));
    1053:     PetscCall(PetscOptionsReal("-pc_ml_AuxThreshold", "Auxiliary smoother drop tol", "None", pc_ml->AuxThreshold, &pc_ml->AuxThreshold, NULL));
    1054:   }
    1055:   PetscOptionsHeadEnd();
    1056:   PetscFunctionReturn(PETSC_SUCCESS);
    1057: }

    1059: /*
    1060:    PCCreate_ML - Creates a ML preconditioner context, PC_ML,
    1061:    and sets this as the private data within the generic preconditioning
    1062:    context, PC, that was created within PCCreate().

    1064:    Input Parameter:
    1065: .  pc - the preconditioner context

    1067:    Application Interface Routine: PCCreate()
    1068: */

    1070: /*MC
    1071:      PCML - Use the SNL ML algebraic multigrid preconditioner.

    1073:    Options Database Keys:
    1074:    Multigrid options(inherited):
    1075: +  -pc_mg_cycle_type (v|w)                            - v for V cycle, w for W-cycle (`PCMGSetCycleType()`)
    1076: .  -pc_mg_distinct_smoothup                           - Should one configure the up and down smoothers separately (`PCMGSetDistinctSmoothUp()`)
    1077: -  -pc_mg_type (additive|multiplicative|full|kascade) - type of multgrid

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

    1096:    Level: intermediate

    1098:    Developer Note:
    1099:    The coarser grid matrices and restriction/interpolation
    1100:    operators are computed by ML, with the matrices converted to PETSc matrices in `MATAIJ` format
    1101:    and the restriction/interpolation operators wrapped as PETSc shell matrices.

    1103: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCMGType`, `PCMG`, `PCHYPRE`, `PCGAMG`,
    1104:           `PCMGSetLevels()`, `PCMGGetLevels()`, `PCMGSetType()`, `MPSetCycles()`, `PCMGSetDistinctSmoothUp()`,
    1105:           `PCMGGetCoarseSolve()`, `PCMGSetResidual()`, `PCMGSetInterpolation()`,
    1106:           `PCMGSetRestriction()`, `PCMGGetSmoother()`, `PCMGGetSmootherUp()`, `PCMGGetSmootherDown()`,
    1107:           `PCMGSetCycleTypeOnLevel()`, `PCMGSetRhs()`, `PCMGSetX()`, `PCMGSetR()`
    1108: M*/

    1110: PETSC_EXTERN PetscErrorCode PCCreate_ML(PC pc)
    1111: {
    1112:   PC_ML *pc_ml;
    1113:   PC_MG *mg;

    1115:   PetscFunctionBegin;
    1116:   /* PCML is an inherited class of PCMG. Initialize pc as PCMG */
    1117:   PetscCall(PCSetType(pc, PCMG)); /* calls PCCreate_MG() and MGCreate_Private() */
    1118:   PetscCall(PetscObjectChangeTypeName((PetscObject)pc, PCML));
    1119:   /* Since PCMG tries to use DM associated with PC must delete it */
    1120:   PetscCall(DMDestroy(&pc->dm));
    1121:   PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_EXTERNAL));
    1122:   mg = (PC_MG *)pc->data;

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

    1128:   pc_ml->ml_object                = 0;
    1129:   pc_ml->agg_object               = 0;
    1130:   pc_ml->gridctx                  = 0;
    1131:   pc_ml->PetscMLdata              = 0;
    1132:   pc_ml->Nlevels                  = -1;
    1133:   pc_ml->MaxNlevels               = 10;
    1134:   pc_ml->MaxCoarseSize            = 1;
    1135:   pc_ml->CoarsenScheme            = 1;
    1136:   pc_ml->Threshold                = 0.0;
    1137:   pc_ml->DampingFactor            = 4.0 / 3.0;
    1138:   pc_ml->SpectralNormScheme_Anorm = PETSC_FALSE;
    1139:   pc_ml->size                     = 0;
    1140:   pc_ml->dim                      = 0;
    1141:   pc_ml->nloc                     = 0;
    1142:   pc_ml->coords                   = 0;
    1143:   pc_ml->Repartition              = PETSC_FALSE;
    1144:   pc_ml->MaxMinRatio              = 1.3;
    1145:   pc_ml->MinPerProc               = 512;
    1146:   pc_ml->PutOnSingleProc          = 5000;
    1147:   pc_ml->RepartitionType          = 0;
    1148:   pc_ml->ZoltanScheme             = 0;
    1149:   pc_ml->Aux                      = PETSC_FALSE;
    1150:   pc_ml->AuxThreshold             = 0.0;

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

    1155:   /* overwrite the pointers of PCMG by the functions of PCML */
    1156:   pc->ops->setfromoptions = PCSetFromOptions_ML;
    1157:   pc->ops->setup          = PCSetUp_ML;
    1158:   pc->ops->reset          = PCReset_ML;
    1159:   pc->ops->destroy        = PCDestroy_ML;
    1160:   PetscFunctionReturn(PETSC_SUCCESS);
    1161: }