Actual source code: ml.c
2: /*
3: Provides an interface to the ML smoothed Aggregation
4: Note: Something non-obvious breaks -pc_mg_type ADDITIVE for parallel runs
5: Jed Brown, see [PETSC #18321, #18449].
6: */
7: #include <petsc/private/pcimpl.h>
8: #include <petsc/private/pcmgimpl.h>
9: #include <../src/mat/impls/aij/seq/aij.h>
10: #include <../src/mat/impls/aij/mpi/mpiaij.h>
11: #include <petscdm.h>
13: EXTERN_C_BEGIN
14: /* HAVE_CONFIG_H flag is required by ML include files */
15: #ifndef HAVE_CONFIG_H
16: #define HAVE_CONFIG_H
17: #endif
18: #include <ml_include.h>
19: #include <ml_viz_stats.h>
20: EXTERN_C_END
22: typedef enum {
23: PCML_NULLSPACE_AUTO,
24: PCML_NULLSPACE_USER,
25: PCML_NULLSPACE_BLOCK,
26: PCML_NULLSPACE_SCALAR
27: } PCMLNullSpaceType;
28: static const char *const PCMLNullSpaceTypes[] = {"AUTO", "USER", "BLOCK", "SCALAR", "PCMLNullSpaceType", "PCML_NULLSPACE_", 0};
30: /* The context (data structure) at each grid level */
31: typedef struct {
32: Vec x, b, r; /* global vectors */
33: Mat A, P, R;
34: KSP ksp;
35: Vec coords; /* projected by ML, if PCSetCoordinates is called; values packed by node */
36: } GridCtx;
38: /* The context used to input PETSc matrix into ML at fine grid */
39: typedef struct {
40: Mat A; /* Petsc matrix in aij format */
41: Mat Aloc; /* local portion of A to be used by ML */
42: Vec x, y;
43: ML_Operator *mlmat;
44: PetscScalar *pwork; /* tmp array used by PetscML_comm() */
45: } FineGridCtx;
47: /* The context associates a ML matrix with a PETSc shell matrix */
48: typedef struct {
49: Mat A; /* PETSc shell matrix associated with mlmat */
50: ML_Operator *mlmat; /* ML matrix assorciated with A */
51: } Mat_MLShell;
53: /* Private context for the ML preconditioner */
54: typedef struct {
55: ML *ml_object;
56: ML_Aggregate *agg_object;
57: GridCtx *gridctx;
58: FineGridCtx *PetscMLdata;
59: PetscInt Nlevels, MaxNlevels, MaxCoarseSize, CoarsenScheme, EnergyMinimization, MinPerProc, PutOnSingleProc, RepartitionType, ZoltanScheme;
60: PetscReal Threshold, DampingFactor, EnergyMinimizationDropTol, MaxMinRatio, AuxThreshold;
61: PetscBool SpectralNormScheme_Anorm, BlockScaling, EnergyMinimizationCheap, Symmetrize, OldHierarchy, KeepAggInfo, Reusable, Repartition, Aux;
62: PetscBool reuse_interpolation;
63: PCMLNullSpaceType nulltype;
64: PetscMPIInt size; /* size of communicator for pc->pmat */
65: PetscInt dim; /* data from PCSetCoordinates(_ML) */
66: PetscInt nloc;
67: PetscReal *coords; /* ML has a grid object for each level: the finest grid will point into coords */
68: } PC_ML;
70: static int PetscML_getrow(ML_Operator *ML_data, int N_requested_rows, int requested_rows[], int allocated_space, int columns[], double values[], int row_lengths[])
71: {
72: PetscInt m, i, j, k = 0, row, *aj;
73: PetscScalar *aa;
74: FineGridCtx *ml = (FineGridCtx *)ML_Get_MyGetrowData(ML_data);
75: Mat_SeqAIJ *a = (Mat_SeqAIJ *)ml->Aloc->data;
77: if (MatGetSize(ml->Aloc, &m, NULL)) return (0);
78: for (i = 0; i < N_requested_rows; i++) {
79: row = requested_rows[i];
80: row_lengths[i] = a->ilen[row];
81: if (allocated_space < k + row_lengths[i]) return (0);
82: if ((row >= 0) || (row <= (m - 1))) {
83: aj = a->j + a->i[row];
84: aa = a->a + a->i[row];
85: for (j = 0; j < row_lengths[i]; j++) {
86: columns[k] = aj[j];
87: values[k++] = aa[j];
88: }
89: }
90: }
91: return (1);
92: }
94: static PetscErrorCode PetscML_comm(double p[], void *ML_data)
95: {
96: FineGridCtx *ml = (FineGridCtx *)ML_data;
97: Mat A = ml->A;
98: Mat_MPIAIJ *a = (Mat_MPIAIJ *)A->data;
99: PetscMPIInt size;
100: PetscInt i, in_length = A->rmap->n, out_length = ml->Aloc->cmap->n;
101: const PetscScalar *array;
103: MPI_Comm_size(PetscObjectComm((PetscObject)A), &size);
104: if (size == 1) return 0;
106: VecPlaceArray(ml->y, p);
107: VecScatterBegin(a->Mvctx, ml->y, a->lvec, INSERT_VALUES, SCATTER_FORWARD);
108: VecScatterEnd(a->Mvctx, ml->y, a->lvec, INSERT_VALUES, SCATTER_FORWARD);
109: VecResetArray(ml->y);
110: VecGetArrayRead(a->lvec, &array);
111: for (i = in_length; i < out_length; i++) p[i] = array[i - in_length];
112: VecRestoreArrayRead(a->lvec, &array);
113: return 0;
114: }
116: static int PetscML_matvec(ML_Operator *ML_data, int in_length, double p[], int out_length, double ap[])
117: {
118: FineGridCtx *ml = (FineGridCtx *)ML_Get_MyMatvecData(ML_data);
119: Mat A = ml->A, Aloc = ml->Aloc;
120: PetscMPIInt size;
121: PetscScalar *pwork = ml->pwork;
122: PetscInt i;
124: MPI_Comm_size(PetscObjectComm((PetscObject)A), &size);
125: if (size == 1) {
126: VecPlaceArray(ml->x, p);
127: } else {
128: for (i = 0; i < in_length; i++) pwork[i] = p[i];
129: PetscML_comm(pwork, ml);
130: VecPlaceArray(ml->x, pwork);
131: }
132: VecPlaceArray(ml->y, ap);
133: MatMult(Aloc, ml->x, ml->y);
134: VecResetArray(ml->x);
135: VecResetArray(ml->y);
136: return 0;
137: }
139: static PetscErrorCode MatMult_ML(Mat A, Vec x, Vec y)
140: {
141: Mat_MLShell *shell;
142: PetscScalar *yarray;
143: const PetscScalar *xarray;
144: PetscInt x_length, y_length;
146: MatShellGetContext(A, &shell);
147: VecGetArrayRead(x, &xarray);
148: VecGetArray(y, &yarray);
149: x_length = shell->mlmat->invec_leng;
150: y_length = shell->mlmat->outvec_leng;
151: PetscStackCallExternalVoid("ML_Operator_Apply", ML_Operator_Apply(shell->mlmat, x_length, (PetscScalar *)xarray, y_length, yarray));
152: VecRestoreArrayRead(x, &xarray);
153: VecRestoreArray(y, &yarray);
154: return 0;
155: }
157: /* newtype is ignored since only handles one case */
158: static PetscErrorCode MatConvert_MPIAIJ_ML(Mat A, MatType newtype, MatReuse scall, Mat *Aloc)
159: {
160: Mat_MPIAIJ *mpimat = (Mat_MPIAIJ *)A->data;
161: Mat_SeqAIJ *mat, *a = (Mat_SeqAIJ *)(mpimat->A)->data, *b = (Mat_SeqAIJ *)(mpimat->B)->data;
162: PetscInt *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j;
163: PetscScalar *aa, *ba, *ca;
164: PetscInt am = A->rmap->n, an = A->cmap->n, i, j, k;
165: PetscInt *ci, *cj, ncols;
168: MatSeqAIJGetArrayRead(mpimat->A, (const PetscScalar **)&aa);
169: MatSeqAIJGetArrayRead(mpimat->B, (const PetscScalar **)&ba);
170: if (scall == MAT_INITIAL_MATRIX) {
171: PetscMalloc1(1 + am, &ci);
172: ci[0] = 0;
173: for (i = 0; i < am; i++) ci[i + 1] = ci[i] + (ai[i + 1] - ai[i]) + (bi[i + 1] - bi[i]);
174: PetscMalloc1(1 + ci[am], &cj);
175: PetscMalloc1(1 + ci[am], &ca);
177: k = 0;
178: for (i = 0; i < am; i++) {
179: /* diagonal portion of A */
180: ncols = ai[i + 1] - ai[i];
181: for (j = 0; j < ncols; j++) {
182: cj[k] = *aj++;
183: ca[k++] = *aa++;
184: }
185: /* off-diagonal portion of A */
186: ncols = bi[i + 1] - bi[i];
187: for (j = 0; j < ncols; j++) {
188: cj[k] = an + (*bj);
189: bj++;
190: ca[k++] = *ba++;
191: }
192: }
195: /* put together the new matrix */
196: an = mpimat->A->cmap->n + mpimat->B->cmap->n;
197: MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, am, an, ci, cj, ca, Aloc);
199: /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
200: /* Since these are PETSc arrays, change flags to free them as necessary. */
201: mat = (Mat_SeqAIJ *)(*Aloc)->data;
202: mat->free_a = PETSC_TRUE;
203: mat->free_ij = PETSC_TRUE;
205: mat->nonew = 0;
206: } else if (scall == MAT_REUSE_MATRIX) {
207: mat = (Mat_SeqAIJ *)(*Aloc)->data;
208: ci = mat->i;
209: cj = mat->j;
210: ca = mat->a;
211: for (i = 0; i < am; i++) {
212: /* diagonal portion of A */
213: ncols = ai[i + 1] - ai[i];
214: for (j = 0; j < ncols; j++) *ca++ = *aa++;
215: /* off-diagonal portion of A */
216: ncols = bi[i + 1] - bi[i];
217: for (j = 0; j < ncols; j++) *ca++ = *ba++;
218: }
219: } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid MatReuse %d", (int)scall);
220: MatSeqAIJRestoreArrayRead(mpimat->A, (const PetscScalar **)&aa);
221: MatSeqAIJRestoreArrayRead(mpimat->B, (const PetscScalar **)&ba);
222: return 0;
223: }
225: static PetscErrorCode MatDestroy_ML(Mat A)
226: {
227: Mat_MLShell *shell;
229: MatShellGetContext(A, &shell);
230: PetscFree(shell);
231: return 0;
232: }
234: static PetscErrorCode MatWrapML_SeqAIJ(ML_Operator *mlmat, MatReuse reuse, Mat *newmat)
235: {
236: struct ML_CSR_MSRdata *matdata = (struct ML_CSR_MSRdata *)mlmat->data;
237: PetscInt m = mlmat->outvec_leng, n = mlmat->invec_leng, *nnz = NULL, nz_max;
238: PetscInt *ml_cols = matdata->columns, *ml_rowptr = matdata->rowptr, *aj, i;
239: PetscScalar *ml_vals = matdata->values, *aa;
242: if (m != n) { /* ML Pmat and Rmat are in CSR format. Pass array pointers into SeqAIJ matrix */
243: if (reuse) {
244: Mat_SeqAIJ *aij = (Mat_SeqAIJ *)(*newmat)->data;
245: aij->i = ml_rowptr;
246: aij->j = ml_cols;
247: aij->a = ml_vals;
248: } else {
249: /* sort ml_cols and ml_vals */
250: PetscMalloc1(m + 1, &nnz);
251: for (i = 0; i < m; i++) nnz[i] = ml_rowptr[i + 1] - ml_rowptr[i];
252: aj = ml_cols;
253: aa = ml_vals;
254: for (i = 0; i < m; i++) {
255: PetscSortIntWithScalarArray(nnz[i], aj, aa);
256: aj += nnz[i];
257: aa += nnz[i];
258: }
259: MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, m, n, ml_rowptr, ml_cols, ml_vals, newmat);
260: PetscFree(nnz);
261: }
262: MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY);
263: MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY);
264: return 0;
265: }
267: nz_max = PetscMax(1, mlmat->max_nz_per_row);
268: PetscMalloc2(nz_max, &aa, nz_max, &aj);
269: if (!reuse) {
270: MatCreate(PETSC_COMM_SELF, newmat);
271: MatSetSizes(*newmat, m, n, PETSC_DECIDE, PETSC_DECIDE);
272: MatSetType(*newmat, MATSEQAIJ);
273: /* keep track of block size for A matrices */
274: MatSetBlockSize(*newmat, mlmat->num_PDEs);
276: PetscMalloc1(m, &nnz);
277: for (i = 0; i < m; i++) PetscStackCallExternalVoid("ML_Operator_Getrow", ML_Operator_Getrow(mlmat, 1, &i, nz_max, aj, aa, &nnz[i]));
278: MatSeqAIJSetPreallocation(*newmat, 0, nnz);
279: }
280: for (i = 0; i < m; i++) {
281: PetscInt ncols;
283: PetscStackCallExternalVoid("ML_Operator_Getrow", ML_Operator_Getrow(mlmat, 1, &i, nz_max, aj, aa, &ncols));
284: MatSetValues(*newmat, 1, &i, ncols, aj, aa, INSERT_VALUES);
285: }
286: MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY);
287: MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY);
289: PetscFree2(aa, aj);
290: PetscFree(nnz);
291: return 0;
292: }
294: static PetscErrorCode MatWrapML_SHELL(ML_Operator *mlmat, MatReuse reuse, Mat *newmat)
295: {
296: PetscInt m, n;
297: ML_Comm *MLcomm;
298: Mat_MLShell *shellctx;
300: m = mlmat->outvec_leng;
301: n = mlmat->invec_leng;
303: if (reuse) {
304: MatShellGetContext(*newmat, &shellctx);
305: shellctx->mlmat = mlmat;
306: return 0;
307: }
309: MLcomm = mlmat->comm;
311: PetscNew(&shellctx);
312: MatCreateShell(MLcomm->USR_comm, m, n, PETSC_DETERMINE, PETSC_DETERMINE, shellctx, newmat);
313: MatShellSetOperation(*newmat, MATOP_MULT, (void (*)(void))MatMult_ML);
314: MatShellSetOperation(*newmat, MATOP_DESTROY, (void (*)(void))MatDestroy_ML);
316: shellctx->A = *newmat;
317: shellctx->mlmat = mlmat;
318: return 0;
319: }
321: static PetscErrorCode MatWrapML_MPIAIJ(ML_Operator *mlmat, MatReuse reuse, Mat *newmat)
322: {
323: PetscInt *aj;
324: PetscScalar *aa;
325: PetscInt i, j, *gordering;
326: PetscInt m = mlmat->outvec_leng, n, nz_max, row;
327: Mat A;
330: n = mlmat->invec_leng;
333: /* create global row numbering for a ML_Operator */
334: PetscStackCallExternalVoid("ML_build_global_numbering", ML_build_global_numbering(mlmat, &gordering, "rows"));
336: nz_max = PetscMax(1, mlmat->max_nz_per_row) + 1;
337: PetscMalloc2(nz_max, &aa, nz_max, &aj);
338: if (reuse) {
339: A = *newmat;
340: } else {
341: PetscInt *nnzA, *nnzB, *nnz;
342: PetscInt rstart;
343: MatCreate(mlmat->comm->USR_comm, &A);
344: MatSetSizes(A, m, n, PETSC_DECIDE, PETSC_DECIDE);
345: MatSetType(A, MATMPIAIJ);
346: /* keep track of block size for A matrices */
347: MatSetBlockSize(A, mlmat->num_PDEs);
348: PetscMalloc3(m, &nnzA, m, &nnzB, m, &nnz);
349: MPI_Scan(&m, &rstart, 1, MPIU_INT, MPI_SUM, mlmat->comm->USR_comm);
350: rstart -= m;
352: for (i = 0; i < m; i++) {
353: row = gordering[i] - rstart;
354: PetscStackCallExternalVoid("ML_Operator_Getrow", ML_Operator_Getrow(mlmat, 1, &i, nz_max, aj, aa, &nnz[i]));
355: nnzA[row] = 0;
356: for (j = 0; j < nnz[i]; j++) {
357: if (aj[j] < m) nnzA[row]++;
358: }
359: nnzB[row] = nnz[i] - nnzA[row];
360: }
361: MatMPIAIJSetPreallocation(A, 0, nnzA, 0, nnzB);
362: PetscFree3(nnzA, nnzB, nnz);
363: }
364: for (i = 0; i < m; i++) {
365: PetscInt ncols;
366: row = gordering[i];
368: PetscStackCallExternalVoid(",ML_Operator_Getrow", ML_Operator_Getrow(mlmat, 1, &i, nz_max, aj, aa, &ncols));
369: for (j = 0; j < ncols; j++) aj[j] = gordering[aj[j]];
370: MatSetValues(A, 1, &row, ncols, aj, aa, INSERT_VALUES);
371: }
372: PetscStackCallExternalVoid("ML_free", ML_free(gordering));
373: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
374: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
375: *newmat = A;
377: PetscFree2(aa, aj);
378: return 0;
379: }
381: /* -------------------------------------------------------------------------- */
382: /*
383: PCSetCoordinates_ML
385: Input Parameter:
386: . pc - the preconditioner context
387: */
388: static PetscErrorCode PCSetCoordinates_ML(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords)
389: {
390: PC_MG *mg = (PC_MG *)pc->data;
391: PC_ML *pc_ml = (PC_ML *)mg->innerctx;
392: PetscInt arrsz, oldarrsz, bs, my0, kk, ii, nloc, Iend, aloc;
393: Mat Amat = pc->pmat;
395: /* this function copied and modified from PCSetCoordinates_GEO -TGI */
397: MatGetBlockSize(Amat, &bs);
399: MatGetOwnershipRange(Amat, &my0, &Iend);
400: aloc = (Iend - my0);
401: nloc = (Iend - my0) / bs;
405: oldarrsz = pc_ml->dim * pc_ml->nloc;
406: pc_ml->dim = ndm;
407: pc_ml->nloc = nloc;
408: arrsz = ndm * nloc;
410: /* create data - syntactic sugar that should be refactored at some point */
411: if (pc_ml->coords == 0 || (oldarrsz != arrsz)) {
412: PetscFree(pc_ml->coords);
413: PetscMalloc1(arrsz, &pc_ml->coords);
414: }
415: for (kk = 0; kk < arrsz; kk++) pc_ml->coords[kk] = -999.;
416: /* copy data in - column oriented */
417: if (nloc == a_nloc) {
418: for (kk = 0; kk < nloc; kk++) {
419: for (ii = 0; ii < ndm; ii++) pc_ml->coords[ii * nloc + kk] = coords[kk * ndm + ii];
420: }
421: } else { /* assumes the coordinates are blocked */
422: for (kk = 0; kk < nloc; kk++) {
423: for (ii = 0; ii < ndm; ii++) pc_ml->coords[ii * nloc + kk] = coords[bs * kk * ndm + ii];
424: }
425: }
426: return 0;
427: }
429: /* -----------------------------------------------------------------------------*/
430: extern PetscErrorCode PCReset_MG(PC);
431: PetscErrorCode PCReset_ML(PC pc)
432: {
433: PC_MG *mg = (PC_MG *)pc->data;
434: PC_ML *pc_ml = (PC_ML *)mg->innerctx;
435: PetscInt level, fine_level = pc_ml->Nlevels - 1, dim = pc_ml->dim;
437: if (dim) {
438: for (level = 0; level <= fine_level; level++) VecDestroy(&pc_ml->gridctx[level].coords);
439: if (pc_ml->ml_object && pc_ml->ml_object->Grid) {
440: ML_Aggregate_Viz_Stats *grid_info = (ML_Aggregate_Viz_Stats *)pc_ml->ml_object->Grid[0].Grid;
441: grid_info->x = 0; /* do this so ML doesn't try to free coordinates */
442: grid_info->y = 0;
443: grid_info->z = 0;
444: PetscStackCallExternalVoid("ML_Operator_Getrow", ML_Aggregate_VizAndStats_Clean(pc_ml->ml_object));
445: }
446: }
447: PetscStackCallExternalVoid("ML_Aggregate_Destroy", ML_Aggregate_Destroy(&pc_ml->agg_object));
448: PetscStackCallExternalVoid("ML_Aggregate_Destroy", ML_Destroy(&pc_ml->ml_object));
450: if (pc_ml->PetscMLdata) {
451: PetscFree(pc_ml->PetscMLdata->pwork);
452: MatDestroy(&pc_ml->PetscMLdata->Aloc);
453: VecDestroy(&pc_ml->PetscMLdata->x);
454: VecDestroy(&pc_ml->PetscMLdata->y);
455: }
456: PetscFree(pc_ml->PetscMLdata);
458: if (pc_ml->gridctx) {
459: for (level = 0; level < fine_level; level++) {
460: if (pc_ml->gridctx[level].A) MatDestroy(&pc_ml->gridctx[level].A);
461: if (pc_ml->gridctx[level].P) MatDestroy(&pc_ml->gridctx[level].P);
462: if (pc_ml->gridctx[level].R) MatDestroy(&pc_ml->gridctx[level].R);
463: if (pc_ml->gridctx[level].x) VecDestroy(&pc_ml->gridctx[level].x);
464: if (pc_ml->gridctx[level].b) VecDestroy(&pc_ml->gridctx[level].b);
465: if (pc_ml->gridctx[level + 1].r) VecDestroy(&pc_ml->gridctx[level + 1].r);
466: }
467: }
468: PetscFree(pc_ml->gridctx);
469: PetscFree(pc_ml->coords);
471: pc_ml->dim = 0;
472: pc_ml->nloc = 0;
473: PCReset_MG(pc);
474: return 0;
475: }
476: /* -------------------------------------------------------------------------- */
477: /*
478: PCSetUp_ML - Prepares for the use of the ML preconditioner
479: by setting data structures and options.
481: Input Parameter:
482: . pc - the preconditioner context
484: Application Interface Routine: PCSetUp()
486: Note:
487: The interface routine PCSetUp() is not usually called directly by
488: the user, but instead is called by PCApply() if necessary.
489: */
490: extern PetscErrorCode PCSetFromOptions_MG(PC, PetscOptionItems *PetscOptionsObject);
491: extern PetscErrorCode PCReset_MG(PC);
493: PetscErrorCode PCSetUp_ML(PC pc)
494: {
495: PetscMPIInt size;
496: FineGridCtx *PetscMLdata;
497: ML *ml_object;
498: ML_Aggregate *agg_object;
499: ML_Operator *mlmat;
500: PetscInt nlocal_allcols, Nlevels, mllevel, level, level1, m, fine_level, bs;
501: Mat A, Aloc;
502: GridCtx *gridctx;
503: PC_MG *mg = (PC_MG *)pc->data;
504: PC_ML *pc_ml = (PC_ML *)mg->innerctx;
505: PetscBool isSeq, isMPI;
506: KSP smoother;
507: PC subpc;
508: PetscInt mesh_level, old_mesh_level;
509: MatInfo info;
510: static PetscBool cite = PETSC_FALSE;
512: PetscCall(PetscCitationsRegister("@TechReport{ml_users_guide,\n author = {M. Sala and J.J. Hu and R.S. Tuminaro},\n title = {{ML}3.1 {S}moothed {A}ggregation {U}ser's {G}uide},\n institution = {Sandia National Laboratories},\n number = "
513: "{SAND2004-4821},\n year = 2004\n}\n",
514: &cite));
515: A = pc->pmat;
516: MPI_Comm_size(PetscObjectComm((PetscObject)A), &size);
518: if (pc->setupcalled) {
519: if (pc->flag == SAME_NONZERO_PATTERN && pc_ml->reuse_interpolation) {
520: /*
521: Reuse interpolaton instead of recomputing aggregates and updating the whole hierarchy. This is less expensive for
522: multiple solves in which the matrix is not changing too quickly.
523: */
524: ml_object = pc_ml->ml_object;
525: gridctx = pc_ml->gridctx;
526: Nlevels = pc_ml->Nlevels;
527: fine_level = Nlevels - 1;
528: gridctx[fine_level].A = A;
530: PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &isSeq);
531: PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &isMPI);
532: if (isMPI) {
533: MatConvert_MPIAIJ_ML(A, NULL, MAT_INITIAL_MATRIX, &Aloc);
534: } else if (isSeq) {
535: Aloc = A;
536: PetscObjectReference((PetscObject)Aloc);
537: } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Matrix type '%s' cannot be used with ML. ML can only handle AIJ matrices.", ((PetscObject)A)->type_name);
539: MatGetSize(Aloc, &m, &nlocal_allcols);
540: PetscMLdata = pc_ml->PetscMLdata;
541: MatDestroy(&PetscMLdata->Aloc);
542: PetscMLdata->A = A;
543: PetscMLdata->Aloc = Aloc;
544: PetscStackCallExternalVoid("ML_Aggregate_Destroy", ML_Init_Amatrix(ml_object, 0, m, m, PetscMLdata));
545: PetscStackCallExternalVoid("ML_Set_Amatrix_Matvec", ML_Set_Amatrix_Matvec(ml_object, 0, PetscML_matvec));
547: mesh_level = ml_object->ML_finest_level;
548: while (ml_object->SingleLevel[mesh_level].Rmat->to) {
549: old_mesh_level = mesh_level;
550: mesh_level = ml_object->SingleLevel[mesh_level].Rmat->to->levelnum;
552: /* clean and regenerate A */
553: mlmat = &(ml_object->Amat[mesh_level]);
554: PetscStackCallExternalVoid("ML_Operator_Clean", ML_Operator_Clean(mlmat));
555: PetscStackCallExternalVoid("ML_Operator_Init", ML_Operator_Init(mlmat, ml_object->comm));
556: PetscStackCallExternalVoid("ML_Gen_AmatrixRAP", ML_Gen_AmatrixRAP(ml_object, old_mesh_level, mesh_level));
557: }
559: level = fine_level - 1;
560: if (size == 1) { /* convert ML P, R and A into seqaij format */
561: for (mllevel = 1; mllevel < Nlevels; mllevel++) {
562: mlmat = &(ml_object->Amat[mllevel]);
563: MatWrapML_SeqAIJ(mlmat, MAT_REUSE_MATRIX, &gridctx[level].A);
564: level--;
565: }
566: } else { /* convert ML P and R into shell format, ML A into mpiaij format */
567: for (mllevel = 1; mllevel < Nlevels; mllevel++) {
568: mlmat = &(ml_object->Amat[mllevel]);
569: MatWrapML_MPIAIJ(mlmat, MAT_REUSE_MATRIX, &gridctx[level].A);
570: level--;
571: }
572: }
574: for (level = 0; level < fine_level; level++) {
575: if (level > 0) PCMGSetResidual(pc, level, PCMGResidualDefault, gridctx[level].A);
576: KSPSetOperators(gridctx[level].ksp, gridctx[level].A, gridctx[level].A);
577: }
578: PCMGSetResidual(pc, fine_level, PCMGResidualDefault, gridctx[fine_level].A);
579: KSPSetOperators(gridctx[fine_level].ksp, gridctx[level].A, gridctx[fine_level].A);
581: PCSetUp_MG(pc);
582: return 0;
583: } else {
584: /* since ML can change the size of vectors/matrices at any level we must destroy everything */
585: PCReset_ML(pc);
586: }
587: }
589: /* setup special features of PCML */
590: /*--------------------------------*/
591: /* convert A to Aloc to be used by ML at fine grid */
592: pc_ml->size = size;
593: PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &isSeq);
594: PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &isMPI);
595: if (isMPI) {
596: MatConvert_MPIAIJ_ML(A, NULL, MAT_INITIAL_MATRIX, &Aloc);
597: } else if (isSeq) {
598: Aloc = A;
599: PetscObjectReference((PetscObject)Aloc);
600: } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Matrix type '%s' cannot be used with ML. ML can only handle AIJ matrices.", ((PetscObject)A)->type_name);
602: /* create and initialize struct 'PetscMLdata' */
603: PetscNew(&PetscMLdata);
604: pc_ml->PetscMLdata = PetscMLdata;
605: PetscMalloc1(Aloc->cmap->n + 1, &PetscMLdata->pwork);
607: MatCreateVecs(Aloc, &PetscMLdata->x, &PetscMLdata->y);
609: PetscMLdata->A = A;
610: PetscMLdata->Aloc = Aloc;
611: if (pc_ml->dim) { /* create vecs around the coordinate data given */
612: PetscInt i, j, dim = pc_ml->dim;
613: PetscInt nloc = pc_ml->nloc, nlocghost;
614: PetscReal *ghostedcoords;
616: MatGetBlockSize(A, &bs);
617: nlocghost = Aloc->cmap->n / bs;
618: PetscMalloc1(dim * nlocghost, &ghostedcoords);
619: for (i = 0; i < dim; i++) {
620: /* copy coordinate values into first component of pwork */
621: for (j = 0; j < nloc; j++) PetscMLdata->pwork[bs * j] = pc_ml->coords[nloc * i + j];
622: /* get the ghost values */
623: PetscML_comm(PetscMLdata->pwork, PetscMLdata);
624: /* write into the vector */
625: for (j = 0; j < nlocghost; j++) ghostedcoords[i * nlocghost + j] = PetscMLdata->pwork[bs * j];
626: }
627: /* replace the original coords with the ghosted coords, because these are
628: * what ML needs */
629: PetscFree(pc_ml->coords);
630: pc_ml->coords = ghostedcoords;
631: }
633: /* create ML discretization matrix at fine grid */
634: /* ML requires input of fine-grid matrix. It determines nlevels. */
635: MatGetSize(Aloc, &m, &nlocal_allcols);
636: MatGetBlockSize(A, &bs);
637: PetscStackCallExternalVoid("ML_Create", ML_Create(&ml_object, pc_ml->MaxNlevels));
638: PetscStackCallExternalVoid("ML_Comm_Set_UsrComm", ML_Comm_Set_UsrComm(ml_object->comm, PetscObjectComm((PetscObject)A)));
639: pc_ml->ml_object = ml_object;
640: PetscStackCallExternalVoid("ML_Init_Amatrix", ML_Init_Amatrix(ml_object, 0, m, m, PetscMLdata));
641: PetscStackCallExternalVoid("ML_Set_Amatrix_Getrow", ML_Set_Amatrix_Getrow(ml_object, 0, PetscML_getrow, PetscML_comm, nlocal_allcols));
642: PetscStackCallExternalVoid("ML_Set_Amatrix_Matvec", ML_Set_Amatrix_Matvec(ml_object, 0, PetscML_matvec));
644: PetscStackCallExternalVoid("ML_Set_Symmetrize", ML_Set_Symmetrize(ml_object, pc_ml->Symmetrize ? ML_YES : ML_NO));
646: /* aggregation */
647: PetscStackCallExternalVoid("ML_Aggregate_Create", ML_Aggregate_Create(&agg_object));
648: pc_ml->agg_object = agg_object;
650: {
651: MatNullSpace mnull;
652: MatGetNearNullSpace(A, &mnull);
653: if (pc_ml->nulltype == PCML_NULLSPACE_AUTO) {
654: if (mnull) pc_ml->nulltype = PCML_NULLSPACE_USER;
655: else if (bs > 1) pc_ml->nulltype = PCML_NULLSPACE_BLOCK;
656: else pc_ml->nulltype = PCML_NULLSPACE_SCALAR;
657: }
658: switch (pc_ml->nulltype) {
659: case PCML_NULLSPACE_USER: {
660: PetscScalar *nullvec;
661: const PetscScalar *v;
662: PetscBool has_const;
663: PetscInt i, j, mlocal, nvec, M;
664: const Vec *vecs;
667: MatGetSize(A, &M, NULL);
668: MatGetLocalSize(Aloc, &mlocal, NULL);
669: MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs);
670: PetscMalloc1((nvec + !!has_const) * mlocal, &nullvec);
671: if (has_const)
672: for (i = 0; i < mlocal; i++) nullvec[i] = 1.0 / M;
673: for (i = 0; i < nvec; i++) {
674: VecGetArrayRead(vecs[i], &v);
675: for (j = 0; j < mlocal; j++) nullvec[(i + !!has_const) * mlocal + j] = v[j];
676: VecRestoreArrayRead(vecs[i], &v);
677: }
678: PetscStackCallExternalVoid("ML_Aggregate_Create", ML_Aggregate_Set_NullSpace(agg_object, bs, nvec + !!has_const, nullvec, mlocal));
679: PetscFree(nullvec);
680: } break;
681: case PCML_NULLSPACE_BLOCK:
682: PetscStackCallExternalVoid("ML_Aggregate_Set_NullSpace", ML_Aggregate_Set_NullSpace(agg_object, bs, bs, 0, 0));
683: break;
684: case PCML_NULLSPACE_SCALAR:
685: break;
686: default:
687: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Unknown null space type");
688: }
689: }
690: PetscStackCallExternalVoid("ML_Aggregate_Set_MaxCoarseSize", ML_Aggregate_Set_MaxCoarseSize(agg_object, pc_ml->MaxCoarseSize));
691: /* set options */
692: switch (pc_ml->CoarsenScheme) {
693: case 1:
694: PetscStackCallExternalVoid("ML_Aggregate_Set_CoarsenScheme_Coupled", ML_Aggregate_Set_CoarsenScheme_Coupled(agg_object));
695: break;
696: case 2:
697: PetscStackCallExternalVoid("ML_Aggregate_Set_CoarsenScheme_MIS", ML_Aggregate_Set_CoarsenScheme_MIS(agg_object));
698: break;
699: case 3:
700: PetscStackCallExternalVoid("ML_Aggregate_Set_CoarsenScheme_METIS", ML_Aggregate_Set_CoarsenScheme_METIS(agg_object));
701: break;
702: }
703: PetscStackCallExternalVoid("ML_Aggregate_Set_Threshold", ML_Aggregate_Set_Threshold(agg_object, pc_ml->Threshold));
704: PetscStackCallExternalVoid("ML_Aggregate_Set_DampingFactor", ML_Aggregate_Set_DampingFactor(agg_object, pc_ml->DampingFactor));
705: if (pc_ml->SpectralNormScheme_Anorm) PetscStackCallExternalVoid("ML_Set_SpectralNormScheme_Anorm", ML_Set_SpectralNormScheme_Anorm(ml_object));
706: agg_object->keep_agg_information = (int)pc_ml->KeepAggInfo;
707: agg_object->keep_P_tentative = (int)pc_ml->Reusable;
708: agg_object->block_scaled_SA = (int)pc_ml->BlockScaling;
709: agg_object->minimizing_energy = (int)pc_ml->EnergyMinimization;
710: agg_object->minimizing_energy_droptol = (double)pc_ml->EnergyMinimizationDropTol;
711: agg_object->cheap_minimizing_energy = (int)pc_ml->EnergyMinimizationCheap;
713: if (pc_ml->Aux) {
715: ml_object->Amat[0].aux_data->threshold = pc_ml->AuxThreshold;
716: ml_object->Amat[0].aux_data->enable = 1;
717: ml_object->Amat[0].aux_data->max_level = 10;
718: ml_object->Amat[0].num_PDEs = bs;
719: }
721: MatGetInfo(A, MAT_LOCAL, &info);
722: ml_object->Amat[0].N_nonzeros = (int)info.nz_used;
724: if (pc_ml->dim) {
725: PetscInt i, dim = pc_ml->dim;
726: ML_Aggregate_Viz_Stats *grid_info;
727: PetscInt nlocghost;
729: MatGetBlockSize(A, &bs);
730: nlocghost = Aloc->cmap->n / bs;
732: PetscStackCallExternalVoid("ML_Aggregate_VizAndStats_Setup(", ML_Aggregate_VizAndStats_Setup(ml_object)); /* create ml info for coords */
733: grid_info = (ML_Aggregate_Viz_Stats *)ml_object->Grid[0].Grid;
734: for (i = 0; i < dim; i++) {
735: /* set the finest level coordinates to point to the column-order array
736: * in pc_ml */
737: /* NOTE: must point away before VizAndStats_Clean so ML doesn't free */
738: switch (i) {
739: case 0:
740: grid_info->x = pc_ml->coords + nlocghost * i;
741: break;
742: case 1:
743: grid_info->y = pc_ml->coords + nlocghost * i;
744: break;
745: case 2:
746: grid_info->z = pc_ml->coords + nlocghost * i;
747: break;
748: default:
749: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_SIZ, "PCML coordinate dimension must be <= 3");
750: }
751: }
752: grid_info->Ndim = dim;
753: }
755: /* repartitioning */
756: if (pc_ml->Repartition) {
757: PetscStackCallExternalVoid("ML_Repartition_Activate", ML_Repartition_Activate(ml_object));
758: PetscStackCallExternalVoid("ML_Repartition_Set_LargestMinMaxRatio", ML_Repartition_Set_LargestMinMaxRatio(ml_object, pc_ml->MaxMinRatio));
759: PetscStackCallExternalVoid("ML_Repartition_Set_MinPerProc", ML_Repartition_Set_MinPerProc(ml_object, pc_ml->MinPerProc));
760: PetscStackCallExternalVoid("ML_Repartition_Set_PutOnSingleProc", ML_Repartition_Set_PutOnSingleProc(ml_object, pc_ml->PutOnSingleProc));
761: #if 0 /* Function not yet defined in ml-6.2 */
762: /* I'm not sure what compatibility issues might crop up if we partitioned
763: * on the finest level, so to be safe repartition starts on the next
764: * finest level (reflection default behavior in
765: * ml_MultiLevelPreconditioner) */
766: PetscStackCallExternalVoid("ML_Repartition_Set_StartLevel",ML_Repartition_Set_StartLevel(ml_object,1));
767: #endif
769: if (!pc_ml->RepartitionType) {
770: PetscInt i;
773: PetscStackCallExternalVoid("ML_Repartition_Set_Partitioner", ML_Repartition_Set_Partitioner(ml_object, ML_USEZOLTAN));
774: PetscStackCallExternalVoid("ML_Aggregate_Set_Dimensions", ML_Aggregate_Set_Dimensions(agg_object, pc_ml->dim));
776: for (i = 0; i < ml_object->ML_num_levels; i++) {
777: ML_Aggregate_Viz_Stats *grid_info = (ML_Aggregate_Viz_Stats *)ml_object->Grid[i].Grid;
778: grid_info->zoltan_type = pc_ml->ZoltanScheme + 1; /* ml numbers options 1, 2, 3 */
779: /* defaults from ml_agg_info.c */
780: grid_info->zoltan_estimated_its = 40; /* only relevant to hypergraph / fast hypergraph */
781: grid_info->zoltan_timers = 0;
782: grid_info->smoothing_steps = 4; /* only relevant to hypergraph / fast hypergraph */
783: }
784: } else {
785: PetscStackCallExternalVoid("ML_Repartition_Set_Partitioner", ML_Repartition_Set_Partitioner(ml_object, ML_USEPARMETIS));
786: }
787: }
789: if (pc_ml->OldHierarchy) {
790: PetscStackCallExternalVoid("ML_Gen_MGHierarchy_UsingAggregation", Nlevels = ML_Gen_MGHierarchy_UsingAggregation(ml_object, 0, ML_INCREASING, agg_object));
791: } else {
792: PetscStackCallExternalVoid("ML_Gen_MultiLevelHierarchy_UsingAggregation", Nlevels = ML_Gen_MultiLevelHierarchy_UsingAggregation(ml_object, 0, ML_INCREASING, agg_object));
793: }
795: pc_ml->Nlevels = Nlevels;
796: fine_level = Nlevels - 1;
798: PCMGSetLevels(pc, Nlevels, NULL);
799: /* set default smoothers */
800: for (level = 1; level <= fine_level; level++) {
801: PCMGGetSmoother(pc, level, &smoother);
802: KSPSetType(smoother, KSPRICHARDSON);
803: KSPGetPC(smoother, &subpc);
804: PCSetType(subpc, PCSOR);
805: }
806: PetscObjectOptionsBegin((PetscObject)pc);
807: PCSetFromOptions_MG(pc, PetscOptionsObject); /* should be called in PCSetFromOptions_ML(), but cannot be called prior to PCMGSetLevels() */
808: PetscOptionsEnd();
810: PetscMalloc1(Nlevels, &gridctx);
812: pc_ml->gridctx = gridctx;
814: /* wrap ML matrices by PETSc shell matrices at coarsened grids.
815: Level 0 is the finest grid for ML, but coarsest for PETSc! */
816: gridctx[fine_level].A = A;
818: level = fine_level - 1;
819: /* TODO: support for GPUs */
820: if (size == 1) { /* convert ML P, R and A into seqaij format */
821: for (mllevel = 1; mllevel < Nlevels; mllevel++) {
822: mlmat = &(ml_object->Pmat[mllevel]);
823: MatWrapML_SeqAIJ(mlmat, MAT_INITIAL_MATRIX, &gridctx[level].P);
824: mlmat = &(ml_object->Rmat[mllevel - 1]);
825: MatWrapML_SeqAIJ(mlmat, MAT_INITIAL_MATRIX, &gridctx[level].R);
827: mlmat = &(ml_object->Amat[mllevel]);
828: MatWrapML_SeqAIJ(mlmat, MAT_INITIAL_MATRIX, &gridctx[level].A);
829: level--;
830: }
831: } else { /* convert ML P and R into shell format, ML A into mpiaij format */
832: for (mllevel = 1; mllevel < Nlevels; mllevel++) {
833: mlmat = &(ml_object->Pmat[mllevel]);
834: MatWrapML_SHELL(mlmat, MAT_INITIAL_MATRIX, &gridctx[level].P);
835: mlmat = &(ml_object->Rmat[mllevel - 1]);
836: MatWrapML_SHELL(mlmat, MAT_INITIAL_MATRIX, &gridctx[level].R);
838: mlmat = &(ml_object->Amat[mllevel]);
839: MatWrapML_MPIAIJ(mlmat, MAT_INITIAL_MATRIX, &gridctx[level].A);
840: level--;
841: }
842: }
844: /* create vectors and ksp at all levels */
845: for (level = 0; level < fine_level; level++) {
846: level1 = level + 1;
848: MatCreateVecs(gridctx[level].A, &gridctx[level].x, &gridctx[level].b);
849: MatCreateVecs(gridctx[level1].A, NULL, &gridctx[level1].r);
850: PCMGSetX(pc, level, gridctx[level].x);
851: PCMGSetRhs(pc, level, gridctx[level].b);
852: PCMGSetR(pc, level1, gridctx[level1].r);
854: if (level == 0) {
855: PCMGGetCoarseSolve(pc, &gridctx[level].ksp);
856: } else {
857: PCMGGetSmoother(pc, level, &gridctx[level].ksp);
858: }
859: }
860: PCMGGetSmoother(pc, fine_level, &gridctx[fine_level].ksp);
862: /* create coarse level and the interpolation between the levels */
863: for (level = 0; level < fine_level; level++) {
864: level1 = level + 1;
866: PCMGSetInterpolation(pc, level1, gridctx[level].P);
867: PCMGSetRestriction(pc, level1, gridctx[level].R);
868: if (level > 0) PCMGSetResidual(pc, level, PCMGResidualDefault, gridctx[level].A);
869: KSPSetOperators(gridctx[level].ksp, gridctx[level].A, gridctx[level].A);
870: }
871: PCMGSetResidual(pc, fine_level, PCMGResidualDefault, gridctx[fine_level].A);
872: KSPSetOperators(gridctx[fine_level].ksp, gridctx[level].A, gridctx[fine_level].A);
874: /* put coordinate info in levels */
875: if (pc_ml->dim) {
876: PetscInt i, j, dim = pc_ml->dim;
877: PetscInt bs, nloc;
878: PC subpc;
879: PetscReal *array;
881: level = fine_level;
882: for (mllevel = 0; mllevel < Nlevels; mllevel++) {
883: ML_Aggregate_Viz_Stats *grid_info = (ML_Aggregate_Viz_Stats *)ml_object->Amat[mllevel].to->Grid->Grid;
884: MPI_Comm comm = ((PetscObject)gridctx[level].A)->comm;
886: MatGetBlockSize(gridctx[level].A, &bs);
887: MatGetLocalSize(gridctx[level].A, NULL, &nloc);
888: nloc /= bs; /* number of local nodes */
890: VecCreate(comm, &gridctx[level].coords);
891: VecSetSizes(gridctx[level].coords, dim * nloc, PETSC_DECIDE);
892: VecSetType(gridctx[level].coords, VECMPI);
893: VecGetArray(gridctx[level].coords, &array);
894: for (j = 0; j < nloc; j++) {
895: for (i = 0; i < dim; i++) {
896: switch (i) {
897: case 0:
898: array[dim * j + i] = grid_info->x[j];
899: break;
900: case 1:
901: array[dim * j + i] = grid_info->y[j];
902: break;
903: case 2:
904: array[dim * j + i] = grid_info->z[j];
905: break;
906: default:
907: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_SIZ, "PCML coordinate dimension must be <= 3");
908: }
909: }
910: }
912: /* passing coordinates to smoothers/coarse solver, should they need them */
913: KSPGetPC(gridctx[level].ksp, &subpc);
914: PCSetCoordinates(subpc, dim, nloc, array);
915: VecRestoreArray(gridctx[level].coords, &array);
916: level--;
917: }
918: }
920: /* setupcalled is set to 0 so that MG is setup from scratch */
921: pc->setupcalled = 0;
922: PCSetUp_MG(pc);
923: return 0;
924: }
926: /* -------------------------------------------------------------------------- */
927: /*
928: PCDestroy_ML - Destroys the private context for the ML preconditioner
929: that was created with PCCreate_ML().
931: Input Parameter:
932: . pc - the preconditioner context
934: Application Interface Routine: PCDestroy()
935: */
936: PetscErrorCode PCDestroy_ML(PC pc)
937: {
938: PC_MG *mg = (PC_MG *)pc->data;
939: PC_ML *pc_ml = (PC_ML *)mg->innerctx;
941: PCReset_ML(pc);
942: PetscFree(pc_ml);
943: PCDestroy_MG(pc);
944: PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL);
945: return 0;
946: }
948: PetscErrorCode PCSetFromOptions_ML(PC pc, PetscOptionItems *PetscOptionsObject)
949: {
950: PetscInt indx, PrintLevel, partindx;
951: const char *scheme[] = {"Uncoupled", "Coupled", "MIS", "METIS"};
952: const char *part[] = {"Zoltan", "ParMETIS"};
953: #if defined(HAVE_ML_ZOLTAN)
954: const char *zscheme[] = {"RCB", "hypergraph", "fast_hypergraph"};
955: #endif
956: PC_MG *mg = (PC_MG *)pc->data;
957: PC_ML *pc_ml = (PC_ML *)mg->innerctx;
958: PetscMPIInt size;
959: MPI_Comm comm;
961: PetscObjectGetComm((PetscObject)pc, &comm);
962: MPI_Comm_size(comm, &size);
963: PetscOptionsHeadBegin(PetscOptionsObject, "ML options");
965: PrintLevel = 0;
966: indx = 0;
967: partindx = 0;
969: PetscOptionsInt("-pc_ml_PrintLevel", "Print level", "ML_Set_PrintLevel", PrintLevel, &PrintLevel, NULL);
970: PetscStackCallExternalVoid("ML_Set_PrintLevel", ML_Set_PrintLevel(PrintLevel));
971: PetscOptionsInt("-pc_ml_maxNlevels", "Maximum number of levels", "None", pc_ml->MaxNlevels, &pc_ml->MaxNlevels, NULL);
972: PetscOptionsInt("-pc_ml_maxCoarseSize", "Maximum coarsest mesh size", "ML_Aggregate_Set_MaxCoarseSize", pc_ml->MaxCoarseSize, &pc_ml->MaxCoarseSize, NULL);
973: PetscOptionsEList("-pc_ml_CoarsenScheme", "Aggregate Coarsen Scheme", "ML_Aggregate_Set_CoarsenScheme_*", scheme, 4, scheme[0], &indx, NULL);
975: pc_ml->CoarsenScheme = indx;
977: PetscOptionsReal("-pc_ml_DampingFactor", "P damping factor", "ML_Aggregate_Set_DampingFactor", pc_ml->DampingFactor, &pc_ml->DampingFactor, NULL);
978: PetscOptionsReal("-pc_ml_Threshold", "Smoother drop tol", "ML_Aggregate_Set_Threshold", pc_ml->Threshold, &pc_ml->Threshold, NULL);
979: PetscOptionsBool("-pc_ml_SpectralNormScheme_Anorm", "Method used for estimating spectral radius", "ML_Set_SpectralNormScheme_Anorm", pc_ml->SpectralNormScheme_Anorm, &pc_ml->SpectralNormScheme_Anorm, NULL);
980: PetscOptionsBool("-pc_ml_Symmetrize", "Symmetrize aggregation", "ML_Set_Symmetrize", pc_ml->Symmetrize, &pc_ml->Symmetrize, NULL);
981: PetscOptionsBool("-pc_ml_BlockScaling", "Scale all dofs at each node together", "None", pc_ml->BlockScaling, &pc_ml->BlockScaling, NULL);
982: PetscOptionsEnum("-pc_ml_nullspace", "Which type of null space information to use", "None", PCMLNullSpaceTypes, (PetscEnum)pc_ml->nulltype, (PetscEnum *)&pc_ml->nulltype, NULL);
983: PetscOptionsInt("-pc_ml_EnergyMinimization", "Energy minimization norm type (0=no minimization; see ML manual for 1,2,3; -1 and 4 undocumented)", "None", pc_ml->EnergyMinimization, &pc_ml->EnergyMinimization, NULL);
984: PetscOptionsBool("-pc_ml_reuse_interpolation", "Reuse the interpolation operators when possible (cheaper, weaker when matrix entries change a lot)", "None", pc_ml->reuse_interpolation, &pc_ml->reuse_interpolation, NULL);
985: /*
986: The following checks a number of conditions. If we let this stuff slip by, then ML's error handling will take over.
987: This is suboptimal because it amounts to calling exit(1) so we check for the most common conditions.
989: We also try to set some sane defaults when energy minimization is activated, otherwise it's hard to find a working
990: combination of options and ML's exit(1) explanations don't help matters.
991: */
994: if (pc_ml->EnergyMinimization == 4) PetscInfo(pc, "Mandel's energy minimization scheme is experimental and broken in ML-6.2\n");
995: if (pc_ml->EnergyMinimization) PetscOptionsReal("-pc_ml_EnergyMinimizationDropTol", "Energy minimization drop tolerance", "None", pc_ml->EnergyMinimizationDropTol, &pc_ml->EnergyMinimizationDropTol, NULL);
996: if (pc_ml->EnergyMinimization == 2) {
997: /* According to ml_MultiLevelPreconditioner.cpp, this option is only meaningful for norm type (2) */
998: PetscOptionsBool("-pc_ml_EnergyMinimizationCheap", "Use cheaper variant of norm type 2", "None", pc_ml->EnergyMinimizationCheap, &pc_ml->EnergyMinimizationCheap, NULL);
999: }
1000: /* energy minimization sometimes breaks if this is turned off, the more classical stuff should be okay without it */
1001: if (pc_ml->EnergyMinimization) pc_ml->KeepAggInfo = PETSC_TRUE;
1002: PetscOptionsBool("-pc_ml_KeepAggInfo", "Allows the preconditioner to be reused, or auxiliary matrices to be generated", "None", pc_ml->KeepAggInfo, &pc_ml->KeepAggInfo, NULL);
1003: /* Option (-1) doesn't work at all (calls exit(1)) if the tentative restriction operator isn't stored. */
1004: if (pc_ml->EnergyMinimization == -1) pc_ml->Reusable = PETSC_TRUE;
1005: PetscOptionsBool("-pc_ml_Reusable", "Store intermedaiate data structures so that the multilevel hierarchy is reusable", "None", pc_ml->Reusable, &pc_ml->Reusable, NULL);
1006: /*
1007: ML's C API is severely underdocumented and lacks significant functionality. The C++ API calls
1008: ML_Gen_MultiLevelHierarchy_UsingAggregation() which is a modified copy (!?) of the documented function
1009: ML_Gen_MGHierarchy_UsingAggregation(). This modification, however, does not provide a strict superset of the
1010: functionality in the old function, so some users may still want to use it. Note that many options are ignored in
1011: this context, but ML doesn't provide a way to find out which ones.
1012: */
1013: PetscOptionsBool("-pc_ml_OldHierarchy", "Use old routine to generate hierarchy", "None", pc_ml->OldHierarchy, &pc_ml->OldHierarchy, NULL);
1014: PetscOptionsBool("-pc_ml_repartition", "Allow ML to repartition levels of the hierarchy", "ML_Repartition_Activate", pc_ml->Repartition, &pc_ml->Repartition, NULL);
1015: if (pc_ml->Repartition) {
1016: PetscOptionsReal("-pc_ml_repartitionMaxMinRatio", "Acceptable ratio of repartitioned sizes", "ML_Repartition_Set_LargestMinMaxRatio", pc_ml->MaxMinRatio, &pc_ml->MaxMinRatio, NULL);
1017: PetscOptionsInt("-pc_ml_repartitionMinPerProc", "Smallest repartitioned size", "ML_Repartition_Set_MinPerProc", pc_ml->MinPerProc, &pc_ml->MinPerProc, NULL);
1018: PetscOptionsInt("-pc_ml_repartitionPutOnSingleProc", "Problem size automatically repartitioned to one processor", "ML_Repartition_Set_PutOnSingleProc", pc_ml->PutOnSingleProc, &pc_ml->PutOnSingleProc, NULL);
1019: #if defined(HAVE_ML_ZOLTAN)
1020: partindx = 0;
1021: PetscOptionsEList("-pc_ml_repartitionType", "Repartitioning library to use", "ML_Repartition_Set_Partitioner", part, 2, part[0], &partindx, NULL);
1023: pc_ml->RepartitionType = partindx;
1024: if (!partindx) {
1025: PetscInt zindx = 0;
1027: PetscOptionsEList("-pc_ml_repartitionZoltanScheme", "Repartitioning scheme to use", "None", zscheme, 3, zscheme[0], &zindx, NULL);
1029: pc_ml->ZoltanScheme = zindx;
1030: }
1031: #else
1032: partindx = 1;
1033: PetscOptionsEList("-pc_ml_repartitionType", "Repartitioning library to use", "ML_Repartition_Set_Partitioner", part, 2, part[1], &partindx, NULL);
1034: pc_ml->RepartitionType = partindx;
1036: #endif
1037: PetscOptionsBool("-pc_ml_Aux", "Aggregate using auxiliary coordinate-based laplacian", "None", pc_ml->Aux, &pc_ml->Aux, NULL);
1038: PetscOptionsReal("-pc_ml_AuxThreshold", "Auxiliary smoother drop tol", "None", pc_ml->AuxThreshold, &pc_ml->AuxThreshold, NULL);
1039: }
1040: PetscOptionsHeadEnd();
1041: return 0;
1042: }
1044: /*
1045: PCCreate_ML - Creates a ML preconditioner context, PC_ML,
1046: and sets this as the private data within the generic preconditioning
1047: context, PC, that was created within PCCreate().
1049: Input Parameter:
1050: . pc - the preconditioner context
1052: Application Interface Routine: PCCreate()
1053: */
1055: /*MC
1056: PCML - Use the SNL ML algebraic multigrid preconditioner.
1058: Options Database Keys:
1059: Multigrid options(inherited):
1060: + -pc_mg_cycles <1> - 1 for V cycle, 2 for W-cycle (`MGSetCycles()`)
1061: . -pc_mg_distinct_smoothup - Should one configure the up and down smoothers separately (`PCMGSetDistinctSmoothUp()`)
1062: - -pc_mg_type <multiplicative> - (one of) additive multiplicative full kascade
1064: ML Options Database Key:
1065: + -pc_ml_PrintLevel <0> - Print level (`ML_Set_PrintLevel()`)
1066: . -pc_ml_maxNlevels <10> - Maximum number of levels (None)
1067: . -pc_ml_maxCoarseSize <1> - Maximum coarsest mesh size (`ML_Aggregate_Set_MaxCoarseSize()`)
1068: . -pc_ml_CoarsenScheme <Uncoupled> - (one of) Uncoupled Coupled MIS METIS
1069: . -pc_ml_DampingFactor <1.33333> - P damping factor (`ML_Aggregate_Set_DampingFactor()`)
1070: . -pc_ml_Threshold <0> - Smoother drop tol (`ML_Aggregate_Set_Threshold()`)
1071: . -pc_ml_SpectralNormScheme_Anorm <false> - Method used for estimating spectral radius (`ML_Set_SpectralNormScheme_Anorm()`)
1072: . -pc_ml_repartition <false> - Allow ML to repartition levels of the hierarchy (`ML_Repartition_Activate()`)
1073: . -pc_ml_repartitionMaxMinRatio <1.3> - Acceptable ratio of repartitioned sizes (`ML_Repartition_Set_LargestMinMaxRatio()`)
1074: . -pc_ml_repartitionMinPerProc <512> - Smallest repartitioned size (`ML_Repartition_Set_MinPerProc()`)
1075: . -pc_ml_repartitionPutOnSingleProc <5000> - Problem size automatically repartitioned to one processor (`ML_Repartition_Set_PutOnSingleProc()`)
1076: . -pc_ml_repartitionType <Zoltan> - Repartitioning library to use (`ML_Repartition_Set_Partitioner()`)
1077: . -pc_ml_repartitionZoltanScheme <RCB> - Repartitioning scheme to use (None)
1078: . -pc_ml_Aux <false> - Aggregate using auxiliary coordinate-based Laplacian (None)
1079: - -pc_ml_AuxThreshold <0.0> - Auxiliary smoother drop tol (None)
1081: Level: intermediate
1083: Developer Note:
1084: The coarser grid matrices and restriction/interpolation
1085: operators are computed by ML, with the matrices converted to PETSc matrices in `MATAIJ` format
1086: and the restriction/interpolation operators wrapped as PETSc shell matrices.
1088: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCMGType`, `PCMG`, `PCHYPRE`, `PCGAMG`,
1089: `PCMGSetLevels()`, `PCMGGetLevels()`, `PCMGSetType()`, `MPSetCycles()`, `PCMGSetDistinctSmoothUp()`,
1090: `PCMGGetCoarseSolve()`, `PCMGSetResidual()`, `PCMGSetInterpolation()`,
1091: `PCMGSetRestriction()`, `PCMGGetSmoother()`, `PCMGGetSmootherUp()`, `PCMGGetSmootherDown()`,
1092: `PCMGSetCycleTypeOnLevel()`, `PCMGSetRhs()`, `PCMGSetX()`, `PCMGSetR()`
1093: M*/
1095: PETSC_EXTERN PetscErrorCode PCCreate_ML(PC pc)
1096: {
1097: PC_ML *pc_ml;
1098: PC_MG *mg;
1100: /* PCML is an inherited class of PCMG. Initialize pc as PCMG */
1101: PCSetType(pc, PCMG); /* calls PCCreate_MG() and MGCreate_Private() */
1102: PetscObjectChangeTypeName((PetscObject)pc, PCML);
1103: /* Since PCMG tries to use DM associated with PC must delete it */
1104: DMDestroy(&pc->dm);
1105: PCMGSetGalerkin(pc, PC_MG_GALERKIN_EXTERNAL);
1106: mg = (PC_MG *)pc->data;
1108: /* create a supporting struct and attach it to pc */
1109: PetscNew(&pc_ml);
1110: mg->innerctx = pc_ml;
1112: pc_ml->ml_object = 0;
1113: pc_ml->agg_object = 0;
1114: pc_ml->gridctx = 0;
1115: pc_ml->PetscMLdata = 0;
1116: pc_ml->Nlevels = -1;
1117: pc_ml->MaxNlevels = 10;
1118: pc_ml->MaxCoarseSize = 1;
1119: pc_ml->CoarsenScheme = 1;
1120: pc_ml->Threshold = 0.0;
1121: pc_ml->DampingFactor = 4.0 / 3.0;
1122: pc_ml->SpectralNormScheme_Anorm = PETSC_FALSE;
1123: pc_ml->size = 0;
1124: pc_ml->dim = 0;
1125: pc_ml->nloc = 0;
1126: pc_ml->coords = 0;
1127: pc_ml->Repartition = PETSC_FALSE;
1128: pc_ml->MaxMinRatio = 1.3;
1129: pc_ml->MinPerProc = 512;
1130: pc_ml->PutOnSingleProc = 5000;
1131: pc_ml->RepartitionType = 0;
1132: pc_ml->ZoltanScheme = 0;
1133: pc_ml->Aux = PETSC_FALSE;
1134: pc_ml->AuxThreshold = 0.0;
1136: /* allow for coordinates to be passed */
1137: PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_ML);
1139: /* overwrite the pointers of PCMG by the functions of PCML */
1140: pc->ops->setfromoptions = PCSetFromOptions_ML;
1141: pc->ops->setup = PCSetUp_ML;
1142: pc->ops->reset = PCReset_ML;
1143: pc->ops->destroy = PCDestroy_ML;
1144: return 0;
1145: }