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