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