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: }