1 2 /* 3 Provides an interface to the LLNL package hypre 4 */ 5 6 /* Must use hypre 2.0.0 or more recent. */ 7 8 #include <petsc-private/pcimpl.h> /*I "petscpc.h" I*/ 9 #include <../src/dm/impls/da/hypre/mhyp.h> 10 #include <_hypre_parcsr_ls.h> 11 12 static PetscBool cite = PETSC_FALSE; 13 static const char hypreCitation[] = "@manual{hypre-web-page,\n title = {{\\sl hypre}: High Performance Preconditioners},\n organization = {Lawrence Livermore National Laboratory},\n note = {\\url{http://www.llnl.gov/CASC/hypre/}}\n}\n"; 14 15 /* 16 Private context (data structure) for the preconditioner. 17 */ 18 typedef struct { 19 HYPRE_Solver hsolver; 20 HYPRE_IJMatrix ij; 21 HYPRE_IJVector b,x; 22 23 HYPRE_Int (*destroy)(HYPRE_Solver); 24 HYPRE_Int (*solve)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector); 25 HYPRE_Int (*setup)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector); 26 27 MPI_Comm comm_hypre; 28 char *hypre_type; 29 30 /* options for Pilut and BoomerAMG*/ 31 PetscInt maxiter; 32 double tol; 33 34 /* options for Pilut */ 35 PetscInt factorrowsize; 36 37 /* options for ParaSails */ 38 PetscInt nlevels; 39 double threshhold; 40 double filter; 41 PetscInt sym; 42 double loadbal; 43 PetscInt logging; 44 PetscInt ruse; 45 PetscInt symt; 46 47 /* options for BoomerAMG */ 48 PetscBool printstatistics; 49 50 /* options for BoomerAMG */ 51 PetscInt cycletype; 52 PetscInt maxlevels; 53 double strongthreshold; 54 double maxrowsum; 55 PetscInt gridsweeps[3]; 56 PetscInt coarsentype; 57 PetscInt measuretype; 58 PetscInt relaxtype[3]; 59 double relaxweight; 60 double outerrelaxweight; 61 PetscInt relaxorder; 62 double truncfactor; 63 PetscBool applyrichardson; 64 PetscInt pmax; 65 PetscInt interptype; 66 PetscInt agg_nl; 67 PetscInt agg_num_paths; 68 PetscInt nodal_coarsen; 69 PetscBool nodal_relax; 70 PetscInt nodal_relax_levels; 71 72 /* options for AMS */ 73 PetscInt ams_print; 74 PetscInt ams_max_iter; 75 PetscInt ams_cycle_type; 76 PetscReal ams_tol; 77 PetscInt ams_relax_type; 78 PetscInt ams_relax_times; 79 PetscReal ams_relax_weight; 80 PetscReal ams_omega; 81 PetscInt ams_amg_alpha_opts[5]; /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for vector Poisson */ 82 PetscReal ams_amg_alpha_theta; /* AMG strength for vector Poisson */ 83 PetscInt ams_amg_beta_opts[5]; /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for scalar Poisson */ 84 PetscReal ams_amg_beta_theta; /* AMG strength for scalar Poisson */ 85 86 /* additional data */ 87 HYPRE_IJVector coords[3]; 88 HYPRE_IJVector constants[3]; 89 HYPRE_IJMatrix G; 90 HYPRE_IJMatrix alpha_Poisson; 91 HYPRE_IJMatrix beta_Poisson; 92 PetscBool ams_beta_is_zero; 93 } PC_HYPRE; 94 95 #undef __FUNCT__ 96 #define __FUNCT__ "PCHYPREGetSolver" 97 PetscErrorCode PCHYPREGetSolver(PC pc,HYPRE_Solver *hsolver) 98 { 99 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 100 101 PetscFunctionBegin; 102 *hsolver = jac->hsolver; 103 PetscFunctionReturn(0); 104 } 105 106 #undef __FUNCT__ 107 #define __FUNCT__ "PCSetUp_HYPRE" 108 static PetscErrorCode PCSetUp_HYPRE(PC pc) 109 { 110 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 111 PetscErrorCode ierr; 112 HYPRE_ParCSRMatrix hmat; 113 HYPRE_ParVector bv,xv; 114 PetscInt bs; 115 116 PetscFunctionBegin; 117 if (!jac->hypre_type) { 118 ierr = PCHYPRESetType(pc,"boomeramg");CHKERRQ(ierr); 119 } 120 121 if (pc->setupcalled) { 122 /* always destroy the old matrix and create a new memory; 123 hope this does not churn the memory too much. The problem 124 is I do not know if it is possible to put the matrix back to 125 its initial state so that we can directly copy the values 126 the second time through. */ 127 PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->ij)); 128 jac->ij = 0; 129 } 130 131 if (!jac->ij) { /* create the matrix the first time through */ 132 ierr = MatHYPRE_IJMatrixCreate(pc->pmat,&jac->ij);CHKERRQ(ierr); 133 } 134 if (!jac->b) { /* create the vectors the first time through */ 135 Vec x,b; 136 ierr = MatCreateVecs(pc->pmat,&x,&b);CHKERRQ(ierr); 137 ierr = VecHYPRE_IJVectorCreate(x,&jac->x);CHKERRQ(ierr); 138 ierr = VecHYPRE_IJVectorCreate(b,&jac->b);CHKERRQ(ierr); 139 ierr = VecDestroy(&x);CHKERRQ(ierr); 140 ierr = VecDestroy(&b);CHKERRQ(ierr); 141 } 142 143 /* special case for BoomerAMG */ 144 if (jac->setup == HYPRE_BoomerAMGSetup) { 145 ierr = MatGetBlockSize(pc->pmat,&bs);CHKERRQ(ierr); 146 if (bs > 1) PetscStackCallStandard(HYPRE_BoomerAMGSetNumFunctions,(jac->hsolver,bs)); 147 } 148 /* special case for AMS */ 149 if (jac->setup == HYPRE_AMSSetup) { 150 if (!jac->coords[0] && !jac->constants[0]) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE AMS preconditioner needs either coordinate vectors via PCSetCoordinates() or edge constant vectors via PCHYPRESetEdgeConstantVectors()"); 151 if (!jac->G) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE AMS preconditioner needs discrete gradient operator via PCHYPRESetDiscreteGradient"); 152 } 153 ierr = MatHYPRE_IJMatrixCopy(pc->pmat,jac->ij);CHKERRQ(ierr); 154 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat)); 155 PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&bv)); 156 PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&xv)); 157 PetscStackCall("HYPRE_SetupXXX",ierr = (*jac->setup)(jac->hsolver,hmat,bv,xv);CHKERRQ(ierr);); 158 PetscFunctionReturn(0); 159 } 160 161 /* 162 Replaces the address where the HYPRE vector points to its data with the address of 163 PETSc's data. Saves the old address so it can be reset when we are finished with it. 164 Allows use to get the data into a HYPRE vector without the cost of memcopies 165 */ 166 #define HYPREReplacePointer(b,newvalue,savedvalue) { \ 167 hypre_ParVector *par_vector = (hypre_ParVector*)hypre_IJVectorObject(((hypre_IJVector*)b)); \ 168 hypre_Vector *local_vector = hypre_ParVectorLocalVector(par_vector); \ 169 savedvalue = local_vector->data; \ 170 local_vector->data = newvalue; \ 171 } 172 173 #undef __FUNCT__ 174 #define __FUNCT__ "PCApply_HYPRE" 175 static PetscErrorCode PCApply_HYPRE(PC pc,Vec b,Vec x) 176 { 177 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 178 PetscErrorCode ierr; 179 HYPRE_ParCSRMatrix hmat; 180 PetscScalar *xv; 181 const PetscScalar *bv,*sbv; 182 HYPRE_ParVector jbv,jxv; 183 PetscScalar *sxv; 184 PetscInt hierr; 185 186 PetscFunctionBegin; 187 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 188 if (!jac->applyrichardson) {ierr = VecSet(x,0.0);CHKERRQ(ierr);} 189 ierr = VecGetArrayRead(b,&bv);CHKERRQ(ierr); 190 ierr = VecGetArray(x,&xv);CHKERRQ(ierr); 191 HYPREReplacePointer(jac->b,(PetscScalar*)bv,sbv); 192 HYPREReplacePointer(jac->x,xv,sxv); 193 194 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat)); 195 PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&jbv)); 196 PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&jxv)); 197 PetscStackCall("Hypre solve",hierr = (*jac->solve)(jac->hsolver,hmat,jbv,jxv); 198 if (hierr && hierr != HYPRE_ERROR_CONV) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr); 199 if (hierr) hypre__global_error = 0;); 200 201 HYPREReplacePointer(jac->b,(PetscScalar*)sbv,bv); 202 HYPREReplacePointer(jac->x,sxv,xv); 203 ierr = VecRestoreArray(x,&xv);CHKERRQ(ierr); 204 ierr = VecRestoreArrayRead(b,&bv);CHKERRQ(ierr); 205 PetscFunctionReturn(0); 206 } 207 208 #undef __FUNCT__ 209 #define __FUNCT__ "PCDestroy_HYPRE" 210 static PetscErrorCode PCDestroy_HYPRE(PC pc) 211 { 212 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 213 PetscErrorCode ierr; 214 215 PetscFunctionBegin; 216 if (jac->ij) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->ij)); 217 if (jac->b) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->b)); 218 if (jac->x) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->x)); 219 if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0])); 220 if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1])); 221 if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2])); 222 if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0])); 223 if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1])); 224 if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2])); 225 if (jac->G) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->G)); 226 if (jac->alpha_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->alpha_Poisson)); 227 if (jac->beta_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->beta_Poisson)); 228 if (jac->destroy) PetscStackCall("HYPRE_DestroyXXX",ierr = (*jac->destroy)(jac->hsolver);CHKERRQ(ierr);); 229 ierr = PetscFree(jac->hypre_type);CHKERRQ(ierr); 230 if (jac->comm_hypre != MPI_COMM_NULL) { ierr = MPI_Comm_free(&(jac->comm_hypre));CHKERRQ(ierr);} 231 ierr = PetscFree(pc->data);CHKERRQ(ierr); 232 233 ierr = PetscObjectChangeTypeName((PetscObject)pc,0);CHKERRQ(ierr); 234 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",NULL);CHKERRQ(ierr); 235 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",NULL);CHKERRQ(ierr); 236 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetCoordinates_C",NULL);CHKERRQ(ierr); 237 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",NULL);CHKERRQ(ierr); 238 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetConstantEdgeVectors_C",NULL);CHKERRQ(ierr); 239 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetAlphaPoissonMatrix_C",NULL);CHKERRQ(ierr); 240 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetBetaPoissonMatrix_C",NULL);CHKERRQ(ierr); 241 PetscFunctionReturn(0); 242 } 243 244 /* --------------------------------------------------------------------------------------------*/ 245 #undef __FUNCT__ 246 #define __FUNCT__ "PCSetFromOptions_HYPRE_Pilut" 247 static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PetscOptions *PetscOptionsObject,PC pc) 248 { 249 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 250 PetscErrorCode ierr; 251 PetscBool flag; 252 253 PetscFunctionBegin; 254 ierr = PetscOptionsHead(PetscOptionsObject,"HYPRE Pilut Options");CHKERRQ(ierr); 255 ierr = PetscOptionsInt("-pc_hypre_pilut_maxiter","Number of iterations","None",jac->maxiter,&jac->maxiter,&flag);CHKERRQ(ierr); 256 if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetMaxIter,(jac->hsolver,jac->maxiter)); 257 ierr = PetscOptionsReal("-pc_hypre_pilut_tol","Drop tolerance","None",jac->tol,&jac->tol,&flag);CHKERRQ(ierr); 258 if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetDropTolerance,(jac->hsolver,jac->tol)); 259 ierr = PetscOptionsInt("-pc_hypre_pilut_factorrowsize","FactorRowSize","None",jac->factorrowsize,&jac->factorrowsize,&flag);CHKERRQ(ierr); 260 if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetFactorRowSize,(jac->hsolver,jac->factorrowsize)); 261 ierr = PetscOptionsTail();CHKERRQ(ierr); 262 PetscFunctionReturn(0); 263 } 264 265 #undef __FUNCT__ 266 #define __FUNCT__ "PCView_HYPRE_Pilut" 267 static PetscErrorCode PCView_HYPRE_Pilut(PC pc,PetscViewer viewer) 268 { 269 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 270 PetscErrorCode ierr; 271 PetscBool iascii; 272 273 PetscFunctionBegin; 274 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 275 if (iascii) { 276 ierr = PetscViewerASCIIPrintf(viewer," HYPRE Pilut preconditioning\n");CHKERRQ(ierr); 277 if (jac->maxiter != PETSC_DEFAULT) { 278 ierr = PetscViewerASCIIPrintf(viewer," HYPRE Pilut: maximum number of iterations %d\n",jac->maxiter);CHKERRQ(ierr); 279 } else { 280 ierr = PetscViewerASCIIPrintf(viewer," HYPRE Pilut: default maximum number of iterations \n");CHKERRQ(ierr); 281 } 282 if (jac->tol != PETSC_DEFAULT) { 283 ierr = PetscViewerASCIIPrintf(viewer," HYPRE Pilut: drop tolerance %g\n",(double)jac->tol);CHKERRQ(ierr); 284 } else { 285 ierr = PetscViewerASCIIPrintf(viewer," HYPRE Pilut: default drop tolerance \n");CHKERRQ(ierr); 286 } 287 if (jac->factorrowsize != PETSC_DEFAULT) { 288 ierr = PetscViewerASCIIPrintf(viewer," HYPRE Pilut: factor row size %d\n",jac->factorrowsize);CHKERRQ(ierr); 289 } else { 290 ierr = PetscViewerASCIIPrintf(viewer," HYPRE Pilut: default factor row size \n");CHKERRQ(ierr); 291 } 292 } 293 PetscFunctionReturn(0); 294 } 295 296 /* --------------------------------------------------------------------------------------------*/ 297 298 #undef __FUNCT__ 299 #define __FUNCT__ "PCApplyTranspose_HYPRE_BoomerAMG" 300 static PetscErrorCode PCApplyTranspose_HYPRE_BoomerAMG(PC pc,Vec b,Vec x) 301 { 302 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 303 PetscErrorCode ierr; 304 HYPRE_ParCSRMatrix hmat; 305 PetscScalar *xv; 306 const PetscScalar *bv; 307 HYPRE_ParVector jbv,jxv; 308 PetscScalar *sbv,*sxv; 309 PetscInt hierr; 310 311 PetscFunctionBegin; 312 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 313 ierr = VecSet(x,0.0);CHKERRQ(ierr); 314 ierr = VecGetArrayRead(b,&bv);CHKERRQ(ierr); 315 ierr = VecGetArray(x,&xv);CHKERRQ(ierr); 316 HYPREReplacePointer(jac->b,(PetscScalar*)bv,sbv); 317 HYPREReplacePointer(jac->x,xv,sxv); 318 319 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat)); 320 PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&jbv)); 321 PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&jxv)); 322 323 hierr = HYPRE_BoomerAMGSolveT(jac->hsolver,hmat,jbv,jxv); 324 /* error code of 1 in BoomerAMG merely means convergence not achieved */ 325 if (hierr && (hierr != 1)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr); 326 if (hierr) hypre__global_error = 0; 327 328 HYPREReplacePointer(jac->b,sbv,bv); 329 HYPREReplacePointer(jac->x,sxv,xv); 330 ierr = VecRestoreArray(x,&xv);CHKERRQ(ierr); 331 ierr = VecRestoreArrayRead(b,&bv);CHKERRQ(ierr); 332 PetscFunctionReturn(0); 333 } 334 335 /* static array length */ 336 #define ALEN(a) (sizeof(a)/sizeof((a)[0])) 337 338 static const char *HYPREBoomerAMGCycleType[] = {"","V","W"}; 339 static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP","Ruge-Stueben","","modifiedRuge-Stueben","","","Falgout", "", "PMIS", "", "HMIS"}; 340 static const char *HYPREBoomerAMGMeasureType[] = {"local","global"}; 341 /* The following corresponds to HYPRE_BoomerAMGSetRelaxType which has many missing numbers in the enum */ 342 static const char *HYPREBoomerAMGRelaxType[] = {"Jacobi","sequential-Gauss-Seidel","seqboundary-Gauss-Seidel","SOR/Jacobi","backward-SOR/Jacobi", 343 "" /* [5] hybrid chaotic Gauss-Seidel (works only with OpenMP) */,"symmetric-SOR/Jacobi", 344 "" /* 7 */,"l1scaled-SOR/Jacobi","Gaussian-elimination", 345 "" /* 10 */, "" /* 11 */, "" /* 12 */, "" /* 13 */, "" /* 14 */, 346 "CG" /* non-stationary */,"Chebyshev","FCF-Jacobi","l1scaled-Jacobi"}; 347 static const char *HYPREBoomerAMGInterpType[] = {"classical", "", "", "direct", "multipass", "multipass-wts", "ext+i", 348 "ext+i-cc", "standard", "standard-wts", "", "", "FF", "FF1"}; 349 #undef __FUNCT__ 350 #define __FUNCT__ "PCSetFromOptions_HYPRE_BoomerAMG" 351 static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PetscOptions *PetscOptionsObject,PC pc) 352 { 353 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 354 PetscErrorCode ierr; 355 PetscInt n,indx,level; 356 PetscBool flg, tmp_truth; 357 double tmpdbl, twodbl[2]; 358 359 PetscFunctionBegin; 360 ierr = PetscOptionsHead(PetscOptionsObject,"HYPRE BoomerAMG Options");CHKERRQ(ierr); 361 ierr = PetscOptionsEList("-pc_hypre_boomeramg_cycle_type","Cycle type","None",HYPREBoomerAMGCycleType+1,2,HYPREBoomerAMGCycleType[jac->cycletype],&indx,&flg);CHKERRQ(ierr); 362 if (flg) { 363 jac->cycletype = indx+1; 364 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype)); 365 } 366 ierr = PetscOptionsInt("-pc_hypre_boomeramg_max_levels","Number of levels (of grids) allowed","None",jac->maxlevels,&jac->maxlevels,&flg);CHKERRQ(ierr); 367 if (flg) { 368 if (jac->maxlevels < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be at least two",jac->maxlevels); 369 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels)); 370 } 371 ierr = PetscOptionsInt("-pc_hypre_boomeramg_max_iter","Maximum iterations used PER hypre call","None",jac->maxiter,&jac->maxiter,&flg);CHKERRQ(ierr); 372 if (flg) { 373 if (jac->maxiter < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of iterations %d must be at least one",jac->maxiter); 374 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter)); 375 } 376 ierr = PetscOptionsScalar("-pc_hypre_boomeramg_tol","Convergence tolerance PER hypre call (0.0 = use a fixed number of iterations)","None",jac->tol,&jac->tol,&flg);CHKERRQ(ierr); 377 if (flg) { 378 if (jac->tol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Tolerance %g must be greater than or equal to zero",(double)jac->tol); 379 PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol)); 380 } 381 382 ierr = PetscOptionsScalar("-pc_hypre_boomeramg_truncfactor","Truncation factor for interpolation (0=no truncation)","None",jac->truncfactor,&jac->truncfactor,&flg);CHKERRQ(ierr); 383 if (flg) { 384 if (jac->truncfactor < 0.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Truncation factor %g must be great than or equal zero",(double)jac->truncfactor); 385 PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor)); 386 } 387 388 ierr = PetscOptionsInt("-pc_hypre_boomeramg_P_max","Max elements per row for interpolation operator (0=unlimited)","None",jac->pmax,&jac->pmax,&flg);CHKERRQ(ierr); 389 if (flg) { 390 if (jac->pmax < 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"P_max %g must be greater than or equal to zero",(double)jac->pmax); 391 PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax)); 392 } 393 394 ierr = PetscOptionsInt("-pc_hypre_boomeramg_agg_nl","Number of levels of aggressive coarsening","None",jac->agg_nl,&jac->agg_nl,&flg);CHKERRQ(ierr); 395 if (flg) { 396 if (jac->agg_nl < 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %g must be greater than or equal to zero",(double)jac->agg_nl); 397 398 PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl)); 399 } 400 401 402 ierr = PetscOptionsInt("-pc_hypre_boomeramg_agg_num_paths","Number of paths for aggressive coarsening","None",jac->agg_num_paths,&jac->agg_num_paths,&flg);CHKERRQ(ierr); 403 if (flg) { 404 if (jac->agg_num_paths < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of paths %g must be greater than or equal to 1",(double)jac->agg_num_paths); 405 406 PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths)); 407 } 408 409 410 ierr = PetscOptionsScalar("-pc_hypre_boomeramg_strong_threshold","Threshold for being strongly connected","None",jac->strongthreshold,&jac->strongthreshold,&flg);CHKERRQ(ierr); 411 if (flg) { 412 if (jac->strongthreshold < 0.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Strong threshold %g must be great than or equal zero",(double)jac->strongthreshold); 413 PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold)); 414 } 415 ierr = PetscOptionsScalar("-pc_hypre_boomeramg_max_row_sum","Maximum row sum","None",jac->maxrowsum,&jac->maxrowsum,&flg);CHKERRQ(ierr); 416 if (flg) { 417 if (jac->maxrowsum < 0.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Maximum row sum %g must be greater than zero",(double)jac->maxrowsum); 418 if (jac->maxrowsum > 1.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Maximum row sum %g must be less than or equal one",(double)jac->maxrowsum); 419 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum)); 420 } 421 422 /* Grid sweeps */ 423 ierr = PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all","Number of sweeps for the up and down grid levels","None",jac->gridsweeps[0],&indx,&flg);CHKERRQ(ierr); 424 if (flg) { 425 PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver,indx)); 426 /* modify the jac structure so we can view the updated options with PC_View */ 427 jac->gridsweeps[0] = indx; 428 jac->gridsweeps[1] = indx; 429 /*defaults coarse to 1 */ 430 jac->gridsweeps[2] = 1; 431 } 432 433 ierr = PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down","Number of sweeps for the down cycles","None",jac->gridsweeps[0], &indx,&flg);CHKERRQ(ierr); 434 if (flg) { 435 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 1)); 436 jac->gridsweeps[0] = indx; 437 } 438 ierr = PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up","Number of sweeps for the up cycles","None",jac->gridsweeps[1],&indx,&flg);CHKERRQ(ierr); 439 if (flg) { 440 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 2)); 441 jac->gridsweeps[1] = indx; 442 } 443 ierr = PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse","Number of sweeps for the coarse level","None",jac->gridsweeps[2],&indx,&flg);CHKERRQ(ierr); 444 if (flg) { 445 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 3)); 446 jac->gridsweeps[2] = indx; 447 } 448 449 /* Relax type */ 450 ierr = PetscOptionsEList("-pc_hypre_boomeramg_relax_type_all","Relax type for the up and down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);CHKERRQ(ierr); 451 if (flg) { 452 jac->relaxtype[0] = jac->relaxtype[1] = indx; 453 PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, indx)); 454 /* by default, coarse type set to 9 */ 455 jac->relaxtype[2] = 9; 456 457 } 458 ierr = PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down","Relax type for the down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);CHKERRQ(ierr); 459 if (flg) { 460 jac->relaxtype[0] = indx; 461 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 1)); 462 } 463 ierr = PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up","Relax type for the up cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);CHKERRQ(ierr); 464 if (flg) { 465 jac->relaxtype[1] = indx; 466 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 2)); 467 } 468 ierr = PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse","Relax type on coarse grid","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[9],&indx,&flg);CHKERRQ(ierr); 469 if (flg) { 470 jac->relaxtype[2] = indx; 471 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 3)); 472 } 473 474 /* Relaxation Weight */ 475 ierr = PetscOptionsReal("-pc_hypre_boomeramg_relax_weight_all","Relaxation weight for all levels (0 = hypre estimates, -k = determined with k CG steps)","None",jac->relaxweight, &tmpdbl,&flg);CHKERRQ(ierr); 476 if (flg) { 477 PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxWt,(jac->hsolver,tmpdbl)); 478 jac->relaxweight = tmpdbl; 479 } 480 481 n = 2; 482 twodbl[0] = twodbl[1] = 1.0; 483 ierr = PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level","Set the relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);CHKERRQ(ierr); 484 if (flg) { 485 if (n == 2) { 486 indx = (int)PetscAbsReal(twodbl[1]); 487 PetscStackCallStandard(HYPRE_BoomerAMGSetLevelRelaxWt,(jac->hsolver,twodbl[0],indx)); 488 } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Relax weight level: you must provide 2 values separated by a comma (and no space), you provided %d",n); 489 } 490 491 /* Outer relaxation Weight */ 492 ierr = PetscOptionsReal("-pc_hypre_boomeramg_outer_relax_weight_all","Outer relaxation weight for all levels (-k = determined with k CG steps)","None",jac->outerrelaxweight, &tmpdbl,&flg);CHKERRQ(ierr); 493 if (flg) { 494 PetscStackCallStandard(HYPRE_BoomerAMGSetOuterWt,(jac->hsolver, tmpdbl)); 495 jac->outerrelaxweight = tmpdbl; 496 } 497 498 n = 2; 499 twodbl[0] = twodbl[1] = 1.0; 500 ierr = PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level","Set the outer relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);CHKERRQ(ierr); 501 if (flg) { 502 if (n == 2) { 503 indx = (int)PetscAbsReal(twodbl[1]); 504 PetscStackCallStandard(HYPRE_BoomerAMGSetLevelOuterWt,(jac->hsolver, twodbl[0], indx)); 505 } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Relax weight outer level: You must provide 2 values separated by a comma (and no space), you provided %d",n); 506 } 507 508 /* the Relax Order */ 509 ierr = PetscOptionsBool("-pc_hypre_boomeramg_no_CF", "Do not use CF-relaxation", "None", PETSC_FALSE, &tmp_truth, &flg);CHKERRQ(ierr); 510 511 if (flg && tmp_truth) { 512 jac->relaxorder = 0; 513 PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder)); 514 } 515 ierr = PetscOptionsEList("-pc_hypre_boomeramg_measure_type","Measure type","None",HYPREBoomerAMGMeasureType,ALEN(HYPREBoomerAMGMeasureType),HYPREBoomerAMGMeasureType[0],&indx,&flg);CHKERRQ(ierr); 516 if (flg) { 517 jac->measuretype = indx; 518 PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype)); 519 } 520 /* update list length 3/07 */ 521 ierr = PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type","Coarsen type","None",HYPREBoomerAMGCoarsenType,ALEN(HYPREBoomerAMGCoarsenType),HYPREBoomerAMGCoarsenType[6],&indx,&flg);CHKERRQ(ierr); 522 if (flg) { 523 jac->coarsentype = indx; 524 PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype)); 525 } 526 527 /* new 3/07 */ 528 ierr = PetscOptionsEList("-pc_hypre_boomeramg_interp_type","Interpolation type","None",HYPREBoomerAMGInterpType,ALEN(HYPREBoomerAMGInterpType),HYPREBoomerAMGInterpType[0],&indx,&flg);CHKERRQ(ierr); 529 if (flg) { 530 jac->interptype = indx; 531 PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype)); 532 } 533 534 ierr = PetscOptionsName("-pc_hypre_boomeramg_print_statistics","Print statistics","None",&flg);CHKERRQ(ierr); 535 if (flg) { 536 level = 3; 537 ierr = PetscOptionsInt("-pc_hypre_boomeramg_print_statistics","Print statistics","None",level,&level,NULL);CHKERRQ(ierr); 538 539 jac->printstatistics = PETSC_TRUE; 540 PetscStackCallStandard(HYPRE_BoomerAMGSetPrintLevel,(jac->hsolver,level)); 541 } 542 543 ierr = PetscOptionsName("-pc_hypre_boomeramg_print_debug","Print debug information","None",&flg);CHKERRQ(ierr); 544 if (flg) { 545 level = 3; 546 ierr = PetscOptionsInt("-pc_hypre_boomeramg_print_debug","Print debug information","None",level,&level,NULL);CHKERRQ(ierr); 547 548 jac->printstatistics = PETSC_TRUE; 549 PetscStackCallStandard(HYPRE_BoomerAMGSetDebugFlag,(jac->hsolver,level)); 550 } 551 552 ierr = PetscOptionsBool("-pc_hypre_boomeramg_nodal_coarsen", "HYPRE_BoomerAMGSetNodal()", "None", jac->nodal_coarsen ? PETSC_TRUE : PETSC_FALSE, &tmp_truth, &flg);CHKERRQ(ierr); 553 if (flg && tmp_truth) { 554 jac->nodal_coarsen = 1; 555 PetscStackCallStandard(HYPRE_BoomerAMGSetNodal,(jac->hsolver,1)); 556 } 557 558 ierr = PetscOptionsBool("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg);CHKERRQ(ierr); 559 if (flg && tmp_truth) { 560 PetscInt tmp_int; 561 ierr = PetscOptionsInt("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None",jac->nodal_relax_levels,&tmp_int,&flg);CHKERRQ(ierr); 562 if (flg) jac->nodal_relax_levels = tmp_int; 563 PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,6)); 564 PetscStackCallStandard(HYPRE_BoomerAMGSetDomainType,(jac->hsolver,1)); 565 PetscStackCallStandard(HYPRE_BoomerAMGSetOverlap,(jac->hsolver,0)); 566 PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,jac->nodal_relax_levels)); 567 } 568 569 ierr = PetscOptionsTail();CHKERRQ(ierr); 570 PetscFunctionReturn(0); 571 } 572 573 #undef __FUNCT__ 574 #define __FUNCT__ "PCApplyRichardson_HYPRE_BoomerAMG" 575 static PetscErrorCode PCApplyRichardson_HYPRE_BoomerAMG(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason) 576 { 577 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 578 PetscErrorCode ierr; 579 PetscInt oits; 580 581 PetscFunctionBegin; 582 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 583 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,its*jac->maxiter)); 584 PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,rtol)); 585 jac->applyrichardson = PETSC_TRUE; 586 ierr = PCApply_HYPRE(pc,b,y);CHKERRQ(ierr); 587 jac->applyrichardson = PETSC_FALSE; 588 PetscStackCallStandard(HYPRE_BoomerAMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits)); 589 *outits = oits; 590 if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS; 591 else *reason = PCRICHARDSON_CONVERGED_RTOL; 592 PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol)); 593 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter)); 594 PetscFunctionReturn(0); 595 } 596 597 598 #undef __FUNCT__ 599 #define __FUNCT__ "PCView_HYPRE_BoomerAMG" 600 static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc,PetscViewer viewer) 601 { 602 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 603 PetscErrorCode ierr; 604 PetscBool iascii; 605 606 PetscFunctionBegin; 607 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 608 if (iascii) { 609 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG preconditioning\n");CHKERRQ(ierr); 610 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Cycle type %s\n",HYPREBoomerAMGCycleType[jac->cycletype]);CHKERRQ(ierr); 611 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Maximum number of levels %d\n",jac->maxlevels);CHKERRQ(ierr); 612 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Maximum number of iterations PER hypre call %d\n",jac->maxiter);CHKERRQ(ierr); 613 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Convergence tolerance PER hypre call %g\n",(double)jac->tol);CHKERRQ(ierr); 614 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Threshold for strong coupling %g\n",(double)jac->strongthreshold);CHKERRQ(ierr); 615 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Interpolation truncation factor %g\n",(double)jac->truncfactor);CHKERRQ(ierr); 616 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Interpolation: max elements per row %d\n",jac->pmax);CHKERRQ(ierr); 617 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Number of levels of aggressive coarsening %d\n",jac->agg_nl);CHKERRQ(ierr); 618 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Number of paths for aggressive coarsening %d\n",jac->agg_num_paths);CHKERRQ(ierr); 619 620 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Maximum row sums %g\n",(double)jac->maxrowsum);CHKERRQ(ierr); 621 622 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Sweeps down %d\n",jac->gridsweeps[0]);CHKERRQ(ierr); 623 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Sweeps up %d\n",jac->gridsweeps[1]);CHKERRQ(ierr); 624 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Sweeps on coarse %d\n",jac->gridsweeps[2]);CHKERRQ(ierr); 625 626 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax down %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[0]]);CHKERRQ(ierr); 627 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax up %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[1]]);CHKERRQ(ierr); 628 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax on coarse %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[2]]);CHKERRQ(ierr); 629 630 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax weight (all) %g\n",(double)jac->relaxweight);CHKERRQ(ierr); 631 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Outer relax weight (all) %g\n",(double)jac->outerrelaxweight);CHKERRQ(ierr); 632 633 if (jac->relaxorder) { 634 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Using CF-relaxation\n");CHKERRQ(ierr); 635 } else { 636 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Not using CF-relaxation\n");CHKERRQ(ierr); 637 } 638 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Measure type %s\n",HYPREBoomerAMGMeasureType[jac->measuretype]);CHKERRQ(ierr); 639 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Coarsen type %s\n",HYPREBoomerAMGCoarsenType[jac->coarsentype]);CHKERRQ(ierr); 640 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Interpolation type %s\n",HYPREBoomerAMGInterpType[jac->interptype]);CHKERRQ(ierr); 641 if (jac->nodal_coarsen) { 642 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Using nodal coarsening (with HYPRE_BOOMERAMGSetNodal())\n");CHKERRQ(ierr); 643 } 644 if (jac->nodal_relax) { 645 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Using nodal relaxation via Schwarz smoothing on levels %d\n",jac->nodal_relax_levels);CHKERRQ(ierr); 646 } 647 } 648 PetscFunctionReturn(0); 649 } 650 651 /* --------------------------------------------------------------------------------------------*/ 652 #undef __FUNCT__ 653 #define __FUNCT__ "PCSetFromOptions_HYPRE_ParaSails" 654 static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PetscOptions *PetscOptionsObject,PC pc) 655 { 656 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 657 PetscErrorCode ierr; 658 PetscInt indx; 659 PetscBool flag; 660 const char *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"}; 661 662 PetscFunctionBegin; 663 ierr = PetscOptionsHead(PetscOptionsObject,"HYPRE ParaSails Options");CHKERRQ(ierr); 664 ierr = PetscOptionsInt("-pc_hypre_parasails_nlevels","Number of number of levels","None",jac->nlevels,&jac->nlevels,0);CHKERRQ(ierr); 665 ierr = PetscOptionsReal("-pc_hypre_parasails_thresh","Threshold","None",jac->threshhold,&jac->threshhold,&flag);CHKERRQ(ierr); 666 if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels)); 667 668 ierr = PetscOptionsReal("-pc_hypre_parasails_filter","filter","None",jac->filter,&jac->filter,&flag);CHKERRQ(ierr); 669 if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter)); 670 671 ierr = PetscOptionsReal("-pc_hypre_parasails_loadbal","Load balance","None",jac->loadbal,&jac->loadbal,&flag);CHKERRQ(ierr); 672 if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal)); 673 674 ierr = PetscOptionsBool("-pc_hypre_parasails_logging","Print info to screen","None",(PetscBool)jac->logging,(PetscBool*)&jac->logging,&flag);CHKERRQ(ierr); 675 if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging)); 676 677 ierr = PetscOptionsBool("-pc_hypre_parasails_reuse","Reuse nonzero pattern in preconditioner","None",(PetscBool)jac->ruse,(PetscBool*)&jac->ruse,&flag);CHKERRQ(ierr); 678 if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse)); 679 680 ierr = PetscOptionsEList("-pc_hypre_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,ALEN(symtlist),symtlist[0],&indx,&flag);CHKERRQ(ierr); 681 if (flag) { 682 jac->symt = indx; 683 PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt)); 684 } 685 686 ierr = PetscOptionsTail();CHKERRQ(ierr); 687 PetscFunctionReturn(0); 688 } 689 690 #undef __FUNCT__ 691 #define __FUNCT__ "PCView_HYPRE_ParaSails" 692 static PetscErrorCode PCView_HYPRE_ParaSails(PC pc,PetscViewer viewer) 693 { 694 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 695 PetscErrorCode ierr; 696 PetscBool iascii; 697 const char *symt = 0;; 698 699 PetscFunctionBegin; 700 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 701 if (iascii) { 702 ierr = PetscViewerASCIIPrintf(viewer," HYPRE ParaSails preconditioning\n");CHKERRQ(ierr); 703 ierr = PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: nlevels %d\n",jac->nlevels);CHKERRQ(ierr); 704 ierr = PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: threshold %g\n",(double)jac->threshhold);CHKERRQ(ierr); 705 ierr = PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: filter %g\n",(double)jac->filter);CHKERRQ(ierr); 706 ierr = PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: load balance %g\n",(double)jac->loadbal);CHKERRQ(ierr); 707 ierr = PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: reuse nonzero structure %s\n",PetscBools[jac->ruse]);CHKERRQ(ierr); 708 ierr = PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: print info to screen %s\n",PetscBools[jac->logging]);CHKERRQ(ierr); 709 if (!jac->symt) symt = "nonsymmetric matrix and preconditioner"; 710 else if (jac->symt == 1) symt = "SPD matrix and preconditioner"; 711 else if (jac->symt == 2) symt = "nonsymmetric matrix but SPD preconditioner"; 712 else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Unknown HYPRE ParaSails symmetric option %d",jac->symt); 713 ierr = PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: %s\n",symt);CHKERRQ(ierr); 714 } 715 PetscFunctionReturn(0); 716 } 717 /* --------------------------------------------------------------------------------------------*/ 718 #undef __FUNCT__ 719 #define __FUNCT__ "PCSetFromOptions_HYPRE_AMS" 720 static PetscErrorCode PCSetFromOptions_HYPRE_AMS(PetscOptions *PetscOptionsObject,PC pc) 721 { 722 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 723 PetscErrorCode ierr; 724 PetscInt n; 725 PetscBool flag,flag2,flag3,flag4; 726 727 PetscFunctionBegin; 728 ierr = PetscOptionsHead(PetscOptionsObject,"HYPRE AMS Options");CHKERRQ(ierr); 729 ierr = PetscOptionsInt("-pc_hypre_ams_print_level","Debugging output level for AMS","None",jac->ams_print,&jac->ams_print,&flag);CHKERRQ(ierr); 730 if (flag) PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->ams_print)); 731 ierr = PetscOptionsInt("-pc_hypre_ams_max_iter","Maximum number of AMS multigrid iterations within PCApply","None",jac->ams_max_iter,&jac->ams_max_iter,&flag);CHKERRQ(ierr); 732 if (flag) PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->ams_max_iter)); 733 ierr = PetscOptionsInt("-pc_hypre_ams_cycle_type","Cycle type for AMS multigrid","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag);CHKERRQ(ierr); 734 if (flag) PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type)); 735 ierr = PetscOptionsReal("-pc_hypre_ams_tol","Error tolerance for AMS multigrid","None",jac->ams_tol,&jac->ams_tol,&flag);CHKERRQ(ierr); 736 if (flag) PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->ams_tol)); 737 ierr = PetscOptionsInt("-pc_hypre_ams_relax_type","Relaxation type for AMS smoother","None",jac->ams_relax_type,&jac->ams_relax_type,&flag);CHKERRQ(ierr); 738 ierr = PetscOptionsInt("-pc_hypre_ams_relax_times","Number of relaxation steps for AMS smoother","None",jac->ams_relax_times,&jac->ams_relax_times,&flag2);CHKERRQ(ierr); 739 ierr = PetscOptionsReal("-pc_hypre_ams_relax_weight","Relaxation weight for AMS smoother","None",jac->ams_relax_weight,&jac->ams_relax_weight,&flag3);CHKERRQ(ierr); 740 ierr = PetscOptionsReal("-pc_hypre_ams_omega","SSOR coefficient for AMS smoother","None",jac->ams_omega,&jac->ams_omega,&flag4);CHKERRQ(ierr); 741 if (flag || flag2 || flag3 || flag4) { 742 PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->ams_relax_type, 743 jac->ams_relax_times, 744 jac->ams_relax_weight, 745 jac->ams_omega)); 746 } 747 ierr = PetscOptionsReal("-pc_hypre_ams_amg_alpha_theta","Threshold for strong coupling of vector Poisson AMG solver","None",jac->ams_amg_alpha_theta,&jac->ams_amg_alpha_theta,&flag);CHKERRQ(ierr); 748 n = 5; 749 ierr = PetscOptionsIntArray("-pc_hypre_ams_amg_alpha_options","AMG options for vector Poisson","None",jac->ams_amg_alpha_opts,&n,&flag2);CHKERRQ(ierr); 750 if (flag || flag2) { 751 PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->ams_amg_alpha_opts[0], /* AMG coarsen type */ 752 jac->ams_amg_alpha_opts[1], /* AMG agg_levels */ 753 jac->ams_amg_alpha_opts[2], /* AMG relax_type */ 754 jac->ams_amg_alpha_theta, 755 jac->ams_amg_alpha_opts[3], /* AMG interp_type */ 756 jac->ams_amg_alpha_opts[4])); /* AMG Pmax */ 757 } 758 ierr = PetscOptionsReal("-pc_hypre_ams_amg_beta_theta","Threshold for strong coupling of scalar Poisson AMG solver","None",jac->ams_amg_beta_theta,&jac->ams_amg_beta_theta,&flag);CHKERRQ(ierr); 759 n = 5; 760 ierr = PetscOptionsIntArray("-pc_hypre_ams_amg_beta_options","AMG options for scalar Poisson solver","None",jac->ams_amg_beta_opts,&n,&flag2);CHKERRQ(ierr); 761 if (flag || flag2) { 762 PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->ams_amg_beta_opts[0], /* AMG coarsen type */ 763 jac->ams_amg_beta_opts[1], /* AMG agg_levels */ 764 jac->ams_amg_beta_opts[2], /* AMG relax_type */ 765 jac->ams_amg_beta_theta, 766 jac->ams_amg_beta_opts[3], /* AMG interp_type */ 767 jac->ams_amg_beta_opts[4])); /* AMG Pmax */ 768 } 769 ierr = PetscOptionsTail();CHKERRQ(ierr); 770 PetscFunctionReturn(0); 771 } 772 773 #undef __FUNCT__ 774 #define __FUNCT__ "PCView_HYPRE_AMS" 775 static PetscErrorCode PCView_HYPRE_AMS(PC pc,PetscViewer viewer) 776 { 777 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 778 PetscErrorCode ierr; 779 PetscBool iascii; 780 781 PetscFunctionBegin; 782 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 783 if (iascii) { 784 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS preconditioning\n");CHKERRQ(ierr); 785 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS: subspace iterations per application %d\n",jac->ams_max_iter);CHKERRQ(ierr); 786 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS: subspace cycle type %d\n",jac->ams_cycle_type);CHKERRQ(ierr); 787 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS: subspace iteration tolerance %g\n",jac->ams_tol);CHKERRQ(ierr); 788 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS: smoother type %d\n",jac->ams_relax_type);CHKERRQ(ierr); 789 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS: number of smoothing steps %d\n",jac->ams_relax_times);CHKERRQ(ierr); 790 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS: smoother weight %g\n",jac->ams_relax_weight);CHKERRQ(ierr); 791 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS: smoother omega %g\n",jac->ams_omega);CHKERRQ(ierr); 792 if (jac->alpha_Poisson) { 793 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS: vector Poisson solver (passed in by user)\n");CHKERRQ(ierr); 794 } else { 795 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS: vector Poisson solver (computed) \n");CHKERRQ(ierr); 796 } 797 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS: boomerAMG coarsening type %d\n",jac->ams_amg_alpha_opts[0]);CHKERRQ(ierr); 798 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS: boomerAMG levels of aggressive coarsening %d\n",jac->ams_amg_alpha_opts[1]);CHKERRQ(ierr); 799 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS: boomerAMG relaxation type %d\n",jac->ams_amg_alpha_opts[2]);CHKERRQ(ierr); 800 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS: boomerAMG interpolation type %d\n",jac->ams_amg_alpha_opts[3]);CHKERRQ(ierr); 801 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS: boomerAMG max nonzero elements in interpolation rows %d\n",jac->ams_amg_alpha_opts[4]);CHKERRQ(ierr); 802 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS: boomerAMG strength threshold %g\n",jac->ams_amg_alpha_theta);CHKERRQ(ierr); 803 if (!jac->ams_beta_is_zero) { 804 if (jac->beta_Poisson) { 805 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS: scalar Poisson solver (passed in by user)\n");CHKERRQ(ierr); 806 } else { 807 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS: scalar Poisson solver (computed) \n");CHKERRQ(ierr); 808 } 809 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS: boomerAMG coarsening type %d\n",jac->ams_amg_beta_opts[0]);CHKERRQ(ierr); 810 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS: boomerAMG levels of aggressive coarsening %d\n",jac->ams_amg_beta_opts[1]);CHKERRQ(ierr); 811 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS: boomerAMG relaxation type %d\n",jac->ams_amg_beta_opts[2]);CHKERRQ(ierr); 812 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS: boomerAMG interpolation type %d\n",jac->ams_amg_beta_opts[3]);CHKERRQ(ierr); 813 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS: boomerAMG max nonzero elements in interpolation rows %d\n",jac->ams_amg_beta_opts[4]);CHKERRQ(ierr); 814 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS: boomerAMG strength threshold %g\n",jac->ams_amg_beta_theta);CHKERRQ(ierr); 815 } 816 } 817 PetscFunctionReturn(0); 818 } 819 820 #undef __FUNCT__ 821 #define __FUNCT__ "PCHYPRESetDiscreteGradient_HYPRE_AMS" 822 static PetscErrorCode PCHYPRESetDiscreteGradient_HYPRE_AMS(PC pc, Mat G) 823 { 824 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 825 HYPRE_ParCSRMatrix parcsr_G; 826 PetscErrorCode ierr; 827 828 PetscFunctionBegin; 829 /* throw away any discrete gradient if already set */ 830 if (jac->G) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->G)); 831 ierr = MatHYPRE_IJMatrixCreate(G,&jac->G);CHKERRQ(ierr); 832 ierr = MatHYPRE_IJMatrixCopy(G,jac->G);CHKERRQ(ierr); 833 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->G,(void**)(&parcsr_G))); 834 PetscStackCallStandard(HYPRE_AMSSetDiscreteGradient,(jac->hsolver,parcsr_G)); 835 PetscFunctionReturn(0); 836 } 837 838 #undef __FUNCT__ 839 #define __FUNCT__ "PCHYPRESetDiscreteGradient" 840 /*@ 841 PCHYPRESetDiscreteGradient - Set discrete gradient matrix 842 843 Collective on PC 844 845 Input Parameters: 846 + pc - the preconditioning context 847 - G - the discrete gradient 848 849 Level: intermediate 850 851 Notes: G should have as many rows as the number of edges and as many columns as the number of vertices in the mesh 852 Each row of G has 2 nonzeros, with column indexes being the global indexes of edge's endpoints: matrix values are +1 and -1 depending on edge orientation 853 854 .seealso: 855 @*/ 856 PetscErrorCode PCHYPRESetDiscreteGradient(PC pc, Mat G) 857 { 858 PetscErrorCode ierr; 859 860 PetscFunctionBegin; 861 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 862 PetscValidHeaderSpecific(G,MAT_CLASSID,2); 863 PetscCheckSameComm(pc,1,G,2); 864 ierr = PetscTryMethod(pc,"PCHYPRESetDiscreteGradient_C",(PC,Mat),(pc,G));CHKERRQ(ierr); 865 PetscFunctionReturn(0); 866 } 867 868 #undef __FUNCT__ 869 #define __FUNCT__ "PCHYPRESetAlphaPoissonMatrix_HYPRE_AMS" 870 static PetscErrorCode PCHYPRESetAlphaPoissonMatrix_HYPRE_AMS(PC pc, Mat A) 871 { 872 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 873 HYPRE_ParCSRMatrix parcsr_alpha_Poisson; 874 PetscErrorCode ierr; 875 876 PetscFunctionBegin; 877 /* throw away any matrix if already set */ 878 if (jac->alpha_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->alpha_Poisson)); 879 ierr = MatHYPRE_IJMatrixCreate(A,&jac->alpha_Poisson);CHKERRQ(ierr); 880 ierr = MatHYPRE_IJMatrixCopy(A,jac->alpha_Poisson);CHKERRQ(ierr); 881 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->alpha_Poisson,(void**)(&parcsr_alpha_Poisson))); 882 PetscStackCallStandard(HYPRE_AMSSetAlphaPoissonMatrix,(jac->hsolver,parcsr_alpha_Poisson)); 883 PetscFunctionReturn(0); 884 } 885 886 #undef __FUNCT__ 887 #define __FUNCT__ "PCHYPRESetAlphaPoissonMatrix" 888 /*@ 889 PCHYPRESetAlphaPoissonMatrix - Set vector Poisson matrix 890 891 Collective on PC 892 893 Input Parameters: 894 + pc - the preconditioning context 895 - A - the matrix 896 897 Level: intermediate 898 899 Notes: A should be obtained by discretizing the vector valued Poisson problem with linear finite elements 900 901 .seealso: 902 @*/ 903 PetscErrorCode PCHYPRESetAlphaPoissonMatrix(PC pc, Mat A) 904 { 905 PetscErrorCode ierr; 906 907 PetscFunctionBegin; 908 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 909 PetscValidHeaderSpecific(A,MAT_CLASSID,2); 910 PetscCheckSameComm(pc,1,A,2); 911 ierr = PetscTryMethod(pc,"PCHYPRESetAlphaPoissonMatrix_C",(PC,Mat),(pc,A));CHKERRQ(ierr); 912 PetscFunctionReturn(0); 913 } 914 915 #undef __FUNCT__ 916 #define __FUNCT__ "PCHYPRESetBetaPoissonMatrix_HYPRE_AMS" 917 static PetscErrorCode PCHYPRESetBetaPoissonMatrix_HYPRE_AMS(PC pc, Mat A) 918 { 919 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 920 HYPRE_ParCSRMatrix parcsr_beta_Poisson; 921 PetscErrorCode ierr; 922 923 PetscFunctionBegin; 924 if (!A) { 925 jac->ams_beta_is_zero = PETSC_TRUE; 926 PetscFunctionReturn(0); 927 } 928 jac->ams_beta_is_zero = PETSC_FALSE; 929 /* throw away any matrix if already set */ 930 if (jac->beta_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->beta_Poisson)); 931 ierr = MatHYPRE_IJMatrixCreate(A,&jac->beta_Poisson);CHKERRQ(ierr); 932 ierr = MatHYPRE_IJMatrixCopy(A,jac->beta_Poisson);CHKERRQ(ierr); 933 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->beta_Poisson,(void**)(&parcsr_beta_Poisson))); 934 PetscStackCallStandard(HYPRE_AMSSetBetaPoissonMatrix,(jac->hsolver,parcsr_beta_Poisson)); 935 PetscFunctionReturn(0); 936 } 937 938 #undef __FUNCT__ 939 #define __FUNCT__ "PCHYPRESetBetaPoissonMatrix" 940 /*@ 941 PCHYPRESetBetaPoissonMatrix - Set Poisson matrix 942 943 Collective on PC 944 945 Input Parameters: 946 + pc - the preconditioning context 947 - A - the matrix 948 949 Level: intermediate 950 951 Notes: A should be obtained by discretizing the Poisson problem with linear finite elements. 952 Following HYPRE convention, the scalar Poisson solver of AMS can be turned off by passing NULL. 953 954 .seealso: 955 @*/ 956 PetscErrorCode PCHYPRESetBetaPoissonMatrix(PC pc, Mat A) 957 { 958 PetscErrorCode ierr; 959 960 PetscFunctionBegin; 961 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 962 if (A) { 963 PetscValidHeaderSpecific(A,MAT_CLASSID,2); 964 PetscCheckSameComm(pc,1,A,2); 965 } 966 ierr = PetscTryMethod(pc,"PCHYPRESetBetaPoissonMatrix_C",(PC,Mat),(pc,A));CHKERRQ(ierr); 967 PetscFunctionReturn(0); 968 } 969 970 #undef __FUNCT__ 971 #define __FUNCT__ "PCHYPRESetEdgeConstantVectors_HYPRE_AMS" 972 static PetscErrorCode PCHYPRESetEdgeConstantVectors_HYPRE_AMS(PC pc,Vec ozz, Vec zoz, Vec zzo) 973 { 974 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 975 HYPRE_ParVector par_ozz,par_zoz,par_zzo; 976 PetscInt dim; 977 PetscErrorCode ierr; 978 979 PetscFunctionBegin; 980 /* throw away any vector if already set */ 981 if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0])); 982 if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1])); 983 if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2])); 984 jac->constants[0] = NULL; 985 jac->constants[1] = NULL; 986 jac->constants[2] = NULL; 987 ierr = VecHYPRE_IJVectorCreate(ozz,&jac->constants[0]);CHKERRQ(ierr); 988 ierr = VecHYPRE_IJVectorCopy(ozz,jac->constants[0]);CHKERRQ(ierr); 989 PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[0],(void**)(&par_ozz))); 990 ierr = VecHYPRE_IJVectorCreate(zoz,&jac->constants[1]);CHKERRQ(ierr); 991 ierr = VecHYPRE_IJVectorCopy(zoz,jac->constants[1]);CHKERRQ(ierr); 992 PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[1],(void**)(&par_zoz))); 993 dim = 2; 994 if (zzo) { 995 ierr = VecHYPRE_IJVectorCreate(zzo,&jac->constants[2]);CHKERRQ(ierr); 996 ierr = VecHYPRE_IJVectorCopy(zzo,jac->constants[2]);CHKERRQ(ierr); 997 PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[2],(void**)(&par_zzo))); 998 dim++; 999 } 1000 PetscStackCallStandard(HYPRE_AMSSetEdgeConstantVectors,(jac->hsolver,par_ozz,par_zoz,par_zzo)); 1001 PetscStackCallStandard(HYPRE_AMSSetDimension,(jac->hsolver,dim)); 1002 PetscFunctionReturn(0); 1003 } 1004 1005 #undef __FUNCT__ 1006 #define __FUNCT__ "PCHYPRESetEdgeConstantVectors" 1007 /*@ 1008 PCHYPRESetEdgeConstantVectors - Set the representation of the constant vector fields in edge element basis 1009 1010 Collective on PC 1011 1012 Input Parameters: 1013 + pc - the preconditioning context 1014 - ozz - vector representing (1,0,0) (or (1,0) in 2D) 1015 - zoz - vector representing (0,1,0) (or (0,1) in 2D) 1016 - zzo - vector representing (0,0,1) (use NULL in 2D) 1017 1018 Level: intermediate 1019 1020 Notes: 1021 1022 .seealso: 1023 @*/ 1024 PetscErrorCode PCHYPRESetEdgeConstantVectors(PC pc, Vec ozz, Vec zoz, Vec zzo) 1025 { 1026 PetscErrorCode ierr; 1027 1028 PetscFunctionBegin; 1029 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1030 PetscValidHeaderSpecific(ozz,VEC_CLASSID,2); 1031 PetscValidHeaderSpecific(zoz,VEC_CLASSID,3); 1032 if (zzo) PetscValidHeaderSpecific(zzo,VEC_CLASSID,4); 1033 PetscCheckSameComm(pc,1,ozz,2); 1034 PetscCheckSameComm(pc,1,zoz,3); 1035 if (zzo) PetscCheckSameComm(pc,1,zzo,4); 1036 ierr = PetscTryMethod(pc,"PCHYPRESetEdgeConstantVectors_C",(PC,Vec,Vec,Vec),(pc,ozz,zoz,zzo));CHKERRQ(ierr); 1037 PetscFunctionReturn(0); 1038 } 1039 1040 #undef __FUNCT__ 1041 #define __FUNCT__ "PCSetCoordinates_HYPRE_AMS" 1042 static PetscErrorCode PCSetCoordinates_HYPRE_AMS(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords) 1043 { 1044 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1045 Vec tv; 1046 HYPRE_ParVector par_coords[3]; 1047 PetscInt i; 1048 PetscErrorCode ierr; 1049 1050 PetscFunctionBegin; 1051 /* throw away any coordinate vector if already set */ 1052 if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0])); 1053 if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1])); 1054 if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2])); 1055 /* set problem's dimension */ 1056 PetscStackCallStandard(HYPRE_AMSSetDimension,(jac->hsolver,dim)); 1057 /* compute IJ vector for coordinates */ 1058 ierr = VecCreate(PetscObjectComm((PetscObject)pc),&tv);CHKERRQ(ierr); 1059 ierr = VecSetType(tv,VECSTANDARD);CHKERRQ(ierr); 1060 ierr = VecSetSizes(tv,nloc,PETSC_DECIDE);CHKERRQ(ierr); 1061 for (i=0;i<dim;i++) { 1062 PetscScalar *array; 1063 PetscInt j; 1064 1065 ierr = VecHYPRE_IJVectorCreate(tv,&jac->coords[i]);CHKERRQ(ierr); 1066 ierr = VecGetArray(tv,&array);CHKERRQ(ierr); 1067 for (j=0;j<nloc;j++) { 1068 array[j] = coords[j*dim+i]; 1069 } 1070 PetscStackCallStandard(HYPRE_IJVectorSetValues,(jac->coords[i],nloc,NULL,array)); 1071 PetscStackCallStandard(HYPRE_IJVectorAssemble,(jac->coords[i])); 1072 ierr = VecRestoreArray(tv,&array);CHKERRQ(ierr); 1073 } 1074 ierr = VecDestroy(&tv);CHKERRQ(ierr); 1075 /* pass parCSR vectors to AMS solver */ 1076 par_coords[0] = NULL; 1077 par_coords[1] = NULL; 1078 par_coords[2] = NULL; 1079 if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[0],(void**)(&par_coords[0]))); 1080 if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[1],(void**)(&par_coords[1]))); 1081 if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[2],(void**)(&par_coords[2]))); 1082 PetscStackCallStandard(HYPRE_AMSSetCoordinateVectors,(jac->hsolver,par_coords[0],par_coords[1],par_coords[2])); 1083 PetscFunctionReturn(0); 1084 } 1085 1086 /* ---------------------------------------------------------------------------------*/ 1087 1088 #undef __FUNCT__ 1089 #define __FUNCT__ "PCHYPREGetType_HYPRE" 1090 static PetscErrorCode PCHYPREGetType_HYPRE(PC pc,const char *name[]) 1091 { 1092 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1093 1094 PetscFunctionBegin; 1095 *name = jac->hypre_type; 1096 PetscFunctionReturn(0); 1097 } 1098 1099 #undef __FUNCT__ 1100 #define __FUNCT__ "PCHYPRESetType_HYPRE" 1101 static PetscErrorCode PCHYPRESetType_HYPRE(PC pc,const char name[]) 1102 { 1103 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1104 PetscErrorCode ierr; 1105 PetscBool flag; 1106 1107 PetscFunctionBegin; 1108 if (jac->hypre_type) { 1109 ierr = PetscStrcmp(jac->hypre_type,name,&flag);CHKERRQ(ierr); 1110 if (!flag) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Cannot reset the HYPRE preconditioner type once it has been set"); 1111 PetscFunctionReturn(0); 1112 } else { 1113 ierr = PetscStrallocpy(name, &jac->hypre_type);CHKERRQ(ierr); 1114 } 1115 1116 jac->maxiter = PETSC_DEFAULT; 1117 jac->tol = PETSC_DEFAULT; 1118 jac->printstatistics = PetscLogPrintInfo; 1119 1120 ierr = PetscStrcmp("pilut",jac->hypre_type,&flag);CHKERRQ(ierr); 1121 if (flag) { 1122 PetscStackCallStandard(HYPRE_ParCSRPilutCreate,(jac->comm_hypre,&jac->hsolver)); 1123 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut; 1124 pc->ops->view = PCView_HYPRE_Pilut; 1125 jac->destroy = HYPRE_ParCSRPilutDestroy; 1126 jac->setup = HYPRE_ParCSRPilutSetup; 1127 jac->solve = HYPRE_ParCSRPilutSolve; 1128 jac->factorrowsize = PETSC_DEFAULT; 1129 PetscFunctionReturn(0); 1130 } 1131 ierr = PetscStrcmp("parasails",jac->hypre_type,&flag);CHKERRQ(ierr); 1132 if (flag) { 1133 PetscStackCallStandard(HYPRE_ParaSailsCreate,(jac->comm_hypre,&jac->hsolver)); 1134 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails; 1135 pc->ops->view = PCView_HYPRE_ParaSails; 1136 jac->destroy = HYPRE_ParaSailsDestroy; 1137 jac->setup = HYPRE_ParaSailsSetup; 1138 jac->solve = HYPRE_ParaSailsSolve; 1139 /* initialize */ 1140 jac->nlevels = 1; 1141 jac->threshhold = .1; 1142 jac->filter = .1; 1143 jac->loadbal = 0; 1144 if (PetscLogPrintInfo) jac->logging = (int) PETSC_TRUE; 1145 else jac->logging = (int) PETSC_FALSE; 1146 1147 jac->ruse = (int) PETSC_FALSE; 1148 jac->symt = 0; 1149 PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels)); 1150 PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter)); 1151 PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal)); 1152 PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging)); 1153 PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse)); 1154 PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt)); 1155 PetscFunctionReturn(0); 1156 } 1157 ierr = PetscStrcmp("boomeramg",jac->hypre_type,&flag);CHKERRQ(ierr); 1158 if (flag) { 1159 ierr = HYPRE_BoomerAMGCreate(&jac->hsolver); 1160 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_BoomerAMG; 1161 pc->ops->view = PCView_HYPRE_BoomerAMG; 1162 pc->ops->applytranspose = PCApplyTranspose_HYPRE_BoomerAMG; 1163 pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG; 1164 jac->destroy = HYPRE_BoomerAMGDestroy; 1165 jac->setup = HYPRE_BoomerAMGSetup; 1166 jac->solve = HYPRE_BoomerAMGSolve; 1167 jac->applyrichardson = PETSC_FALSE; 1168 /* these defaults match the hypre defaults */ 1169 jac->cycletype = 1; 1170 jac->maxlevels = 25; 1171 jac->maxiter = 1; 1172 jac->tol = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */ 1173 jac->truncfactor = 0.0; 1174 jac->strongthreshold = .25; 1175 jac->maxrowsum = .9; 1176 jac->coarsentype = 6; 1177 jac->measuretype = 0; 1178 jac->gridsweeps[0] = jac->gridsweeps[1] = jac->gridsweeps[2] = 1; 1179 jac->relaxtype[0] = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a a PC - most likely with CG */ 1180 jac->relaxtype[2] = 9; /*G.E. */ 1181 jac->relaxweight = 1.0; 1182 jac->outerrelaxweight = 1.0; 1183 jac->relaxorder = 1; 1184 jac->interptype = 0; 1185 jac->agg_nl = 0; 1186 jac->pmax = 0; 1187 jac->truncfactor = 0.0; 1188 jac->agg_num_paths = 1; 1189 1190 jac->nodal_coarsen = 0; 1191 jac->nodal_relax = PETSC_FALSE; 1192 jac->nodal_relax_levels = 1; 1193 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype)); 1194 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels)); 1195 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter)); 1196 PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol)); 1197 PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor)); 1198 PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold)); 1199 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum)); 1200 PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype)); 1201 PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype)); 1202 PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder)); 1203 PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype)); 1204 PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl)); 1205 PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax)); 1206 PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths)); 1207 PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, jac->relaxtype[0])); /*defaults coarse to 9*/ 1208 PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver, jac->gridsweeps[0])); /*defaults coarse to 1 */ 1209 PetscFunctionReturn(0); 1210 } 1211 ierr = PetscStrcmp("ams",jac->hypre_type,&flag);CHKERRQ(ierr); 1212 if (flag) { 1213 ierr = HYPRE_AMSCreate(&jac->hsolver); 1214 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_AMS; 1215 pc->ops->view = PCView_HYPRE_AMS; 1216 jac->destroy = HYPRE_AMSDestroy; 1217 jac->setup = HYPRE_AMSSetup; 1218 jac->solve = HYPRE_AMSSolve; 1219 jac->coords[0] = NULL; 1220 jac->coords[1] = NULL; 1221 jac->coords[2] = NULL; 1222 jac->G = NULL; 1223 /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE AMS */ 1224 jac->ams_print = 0; 1225 jac->ams_max_iter = 1; /* used as a preconditioner */ 1226 jac->ams_cycle_type = 13; 1227 jac->ams_tol = 0.; /* used as a preconditioner */ 1228 /* Smoothing options */ 1229 jac->ams_relax_type = 2; 1230 jac->ams_relax_times = 1; 1231 jac->ams_relax_weight = 1.0; 1232 jac->ams_omega = 1.0; 1233 /* Vector valued Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */ 1234 jac->ams_amg_alpha_opts[0] = 10; 1235 jac->ams_amg_alpha_opts[1] = 1; 1236 jac->ams_amg_alpha_opts[2] = 8; 1237 jac->ams_amg_alpha_opts[3] = 6; 1238 jac->ams_amg_alpha_opts[4] = 4; 1239 jac->ams_amg_alpha_theta = 0.25; 1240 /* Scalar Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */ 1241 jac->ams_beta_is_zero = PETSC_FALSE; 1242 jac->ams_amg_beta_opts[0] = 10; 1243 jac->ams_amg_beta_opts[1] = 1; 1244 jac->ams_amg_beta_opts[2] = 8; 1245 jac->ams_amg_beta_opts[3] = 6; 1246 jac->ams_amg_beta_opts[4] = 4; 1247 jac->ams_amg_beta_theta = 0.25; 1248 PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->ams_print)); 1249 PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->ams_max_iter)); 1250 PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type)); 1251 PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->ams_tol)); 1252 PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->ams_relax_type, 1253 jac->ams_relax_times, 1254 jac->ams_relax_weight, 1255 jac->ams_omega)); 1256 PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->ams_amg_alpha_opts[0], /* AMG coarsen type */ 1257 jac->ams_amg_alpha_opts[1], /* AMG agg_levels */ 1258 jac->ams_amg_alpha_opts[2], /* AMG relax_type */ 1259 jac->ams_amg_alpha_theta, 1260 jac->ams_amg_alpha_opts[3], /* AMG interp_type */ 1261 jac->ams_amg_alpha_opts[4])); /* AMG Pmax */ 1262 PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->ams_amg_beta_opts[0], /* AMG coarsen type */ 1263 jac->ams_amg_beta_opts[1], /* AMG agg_levels */ 1264 jac->ams_amg_beta_opts[2], /* AMG relax_type */ 1265 jac->ams_amg_beta_theta, 1266 jac->ams_amg_beta_opts[3], /* AMG interp_type */ 1267 jac->ams_amg_beta_opts[4])); /* AMG Pmax */ 1268 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HYPRE_AMS);CHKERRQ(ierr); 1269 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",PCHYPRESetDiscreteGradient_HYPRE_AMS);CHKERRQ(ierr); 1270 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetEdgeConstantVectors_C",PCHYPRESetEdgeConstantVectors_HYPRE_AMS);CHKERRQ(ierr); 1271 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetAlphaPoissonMatrix_C",PCHYPRESetAlphaPoissonMatrix_HYPRE_AMS);CHKERRQ(ierr); 1272 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetBetaPoissonMatrix_C",PCHYPRESetBetaPoissonMatrix_HYPRE_AMS);CHKERRQ(ierr); 1273 PetscFunctionReturn(0); 1274 } 1275 ierr = PetscFree(jac->hypre_type);CHKERRQ(ierr); 1276 1277 jac->hypre_type = NULL; 1278 SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown HYPRE preconditioner %s; Choices are pilut, parasails, boomeramg, ams",name); 1279 PetscFunctionReturn(0); 1280 } 1281 1282 /* 1283 It only gets here if the HYPRE type has not been set before the call to 1284 ...SetFromOptions() which actually is most of the time 1285 */ 1286 #undef __FUNCT__ 1287 #define __FUNCT__ "PCSetFromOptions_HYPRE" 1288 static PetscErrorCode PCSetFromOptions_HYPRE(PetscOptions *PetscOptionsObject,PC pc) 1289 { 1290 PetscErrorCode ierr; 1291 PetscInt indx; 1292 const char *type[] = {"pilut","parasails","boomeramg","ams"}; 1293 PetscBool flg; 1294 1295 PetscFunctionBegin; 1296 ierr = PetscOptionsHead(PetscOptionsObject,"HYPRE preconditioner options");CHKERRQ(ierr); 1297 ierr = PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,4,"boomeramg",&indx,&flg);CHKERRQ(ierr); 1298 if (flg) { 1299 ierr = PCHYPRESetType_HYPRE(pc,type[indx]);CHKERRQ(ierr); 1300 } else { 1301 ierr = PCHYPRESetType_HYPRE(pc,"boomeramg");CHKERRQ(ierr); 1302 } 1303 if (pc->ops->setfromoptions) { 1304 ierr = pc->ops->setfromoptions(PetscOptionsObject,pc);CHKERRQ(ierr); 1305 } 1306 ierr = PetscOptionsTail();CHKERRQ(ierr); 1307 PetscFunctionReturn(0); 1308 } 1309 1310 #undef __FUNCT__ 1311 #define __FUNCT__ "PCHYPRESetType" 1312 /*@C 1313 PCHYPRESetType - Sets which hypre preconditioner you wish to use 1314 1315 Input Parameters: 1316 + pc - the preconditioner context 1317 - name - either pilut, parasails, boomeramg, ams 1318 1319 Options Database Keys: 1320 -pc_hypre_type - One of pilut, parasails, boomeramg, ams 1321 1322 Level: intermediate 1323 1324 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 1325 PCHYPRE 1326 1327 @*/ 1328 PetscErrorCode PCHYPRESetType(PC pc,const char name[]) 1329 { 1330 PetscErrorCode ierr; 1331 1332 PetscFunctionBegin; 1333 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1334 PetscValidCharPointer(name,2); 1335 ierr = PetscTryMethod(pc,"PCHYPRESetType_C",(PC,const char[]),(pc,name));CHKERRQ(ierr); 1336 PetscFunctionReturn(0); 1337 } 1338 1339 #undef __FUNCT__ 1340 #define __FUNCT__ "PCHYPREGetType" 1341 /*@C 1342 PCHYPREGetType - Gets which hypre preconditioner you are using 1343 1344 Input Parameter: 1345 . pc - the preconditioner context 1346 1347 Output Parameter: 1348 . name - either pilut, parasails, boomeramg, ams 1349 1350 Level: intermediate 1351 1352 .seealso: PCCreate(), PCHYPRESetType(), PCType (for list of available types), PC, 1353 PCHYPRE 1354 1355 @*/ 1356 PetscErrorCode PCHYPREGetType(PC pc,const char *name[]) 1357 { 1358 PetscErrorCode ierr; 1359 1360 PetscFunctionBegin; 1361 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1362 PetscValidPointer(name,2); 1363 ierr = PetscTryMethod(pc,"PCHYPREGetType_C",(PC,const char*[]),(pc,name));CHKERRQ(ierr); 1364 PetscFunctionReturn(0); 1365 } 1366 1367 /*MC 1368 PCHYPRE - Allows you to use the matrix element based preconditioners in the LLNL package hypre 1369 1370 Options Database Keys: 1371 + -pc_hypre_type - One of pilut, parasails, boomeramg, ams 1372 - Too many others to list, run with -pc_type hypre -pc_hypre_type XXX -help to see options for the XXX 1373 preconditioner 1374 1375 Level: intermediate 1376 1377 Notes: Apart from pc_hypre_type (for which there is PCHYPRESetType()), 1378 the many hypre options can ONLY be set via the options database (e.g. the command line 1379 or with PetscOptionsSetValue(), there are no functions to set them) 1380 1381 The options -pc_hypre_boomeramg_max_iter and -pc_hypre_boomeramg_rtol refer to the number of iterations 1382 (V-cycles) and tolerance that boomeramg does EACH time it is called. So for example, if 1383 -pc_hypre_boomeramg_max_iter is set to 2 then 2-V-cycles are being used to define the preconditioner 1384 (-pc_hypre_boomeramg_rtol should be set to 0.0 - the default - to strictly use a fixed number of 1385 iterations per hypre call). -ksp_max_it and -ksp_rtol STILL determine the total number of iterations 1386 and tolerance for the Krylov solver. For example, if -pc_hypre_boomeramg_max_iter is 2 and -ksp_max_it is 10 1387 then AT MOST twenty V-cycles of boomeramg will be called. 1388 1389 Note that the option -pc_hypre_boomeramg_relax_type_all defaults to symmetric relaxation 1390 (symmetric-SOR/Jacobi), which is required for Krylov solvers like CG that expect symmetry. 1391 Otherwise, you may want to use -pc_hypre_boomeramg_relax_type_all SOR/Jacobi. 1392 If you wish to use BoomerAMG WITHOUT a Krylov method use -ksp_type richardson NOT -ksp_type preonly 1393 and use -ksp_max_it to control the number of V-cycles. 1394 (see the PETSc FAQ.html at the PETSc website under the Documentation tab). 1395 1396 2007-02-03 Using HYPRE-1.11.1b, the routine HYPRE_BoomerAMGSolveT and the option 1397 -pc_hypre_parasails_reuse were failing with SIGSEGV. Dalcin L. 1398 1399 See PCPFMG for access to the hypre Struct PFMG solver 1400 1401 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 1402 PCHYPRESetType(), PCPFMG 1403 1404 M*/ 1405 1406 #undef __FUNCT__ 1407 #define __FUNCT__ "PCCreate_HYPRE" 1408 PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc) 1409 { 1410 PC_HYPRE *jac; 1411 PetscErrorCode ierr; 1412 1413 PetscFunctionBegin; 1414 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 1415 1416 pc->data = jac; 1417 pc->ops->destroy = PCDestroy_HYPRE; 1418 pc->ops->setfromoptions = PCSetFromOptions_HYPRE; 1419 pc->ops->setup = PCSetUp_HYPRE; 1420 pc->ops->apply = PCApply_HYPRE; 1421 jac->comm_hypre = MPI_COMM_NULL; 1422 jac->hypre_type = NULL; 1423 jac->coords[0] = NULL; 1424 jac->coords[1] = NULL; 1425 jac->coords[2] = NULL; 1426 jac->constants[0] = NULL; 1427 jac->constants[1] = NULL; 1428 jac->constants[2] = NULL; 1429 /* duplicate communicator for hypre */ 1430 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));CHKERRQ(ierr); 1431 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",PCHYPRESetType_HYPRE);CHKERRQ(ierr); 1432 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",PCHYPREGetType_HYPRE);CHKERRQ(ierr); 1433 PetscFunctionReturn(0); 1434 } 1435 1436 /* ---------------------------------------------------------------------------------------------------------------------------------*/ 1437 1438 /* this include is needed ONLY to allow access to the private data inside the Mat object specific to hypre */ 1439 #include <petsc-private/matimpl.h> 1440 1441 typedef struct { 1442 MPI_Comm hcomm; /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */ 1443 HYPRE_StructSolver hsolver; 1444 1445 /* keep copy of PFMG options used so may view them */ 1446 PetscInt its; 1447 double tol; 1448 PetscInt relax_type; 1449 PetscInt rap_type; 1450 PetscInt num_pre_relax,num_post_relax; 1451 PetscInt max_levels; 1452 } PC_PFMG; 1453 1454 #undef __FUNCT__ 1455 #define __FUNCT__ "PCDestroy_PFMG" 1456 PetscErrorCode PCDestroy_PFMG(PC pc) 1457 { 1458 PetscErrorCode ierr; 1459 PC_PFMG *ex = (PC_PFMG*) pc->data; 1460 1461 PetscFunctionBegin; 1462 if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver)); 1463 ierr = MPI_Comm_free(&ex->hcomm);CHKERRQ(ierr); 1464 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1465 PetscFunctionReturn(0); 1466 } 1467 1468 static const char *PFMGRelaxType[] = {"Jacobi","Weighted-Jacobi","symmetric-Red/Black-Gauss-Seidel","Red/Black-Gauss-Seidel"}; 1469 static const char *PFMGRAPType[] = {"Galerkin","non-Galerkin"}; 1470 1471 #undef __FUNCT__ 1472 #define __FUNCT__ "PCView_PFMG" 1473 PetscErrorCode PCView_PFMG(PC pc,PetscViewer viewer) 1474 { 1475 PetscErrorCode ierr; 1476 PetscBool iascii; 1477 PC_PFMG *ex = (PC_PFMG*) pc->data; 1478 1479 PetscFunctionBegin; 1480 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1481 if (iascii) { 1482 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG preconditioning\n");CHKERRQ(ierr); 1483 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG: max iterations %d\n",ex->its);CHKERRQ(ierr); 1484 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG: tolerance %g\n",ex->tol);CHKERRQ(ierr); 1485 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);CHKERRQ(ierr); 1486 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG: RAP type %s\n",PFMGRAPType[ex->rap_type]);CHKERRQ(ierr); 1487 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);CHKERRQ(ierr); 1488 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG: max levels %d\n",ex->max_levels);CHKERRQ(ierr); 1489 } 1490 PetscFunctionReturn(0); 1491 } 1492 1493 1494 #undef __FUNCT__ 1495 #define __FUNCT__ "PCSetFromOptions_PFMG" 1496 PetscErrorCode PCSetFromOptions_PFMG(PetscOptions *PetscOptionsObject,PC pc) 1497 { 1498 PetscErrorCode ierr; 1499 PC_PFMG *ex = (PC_PFMG*) pc->data; 1500 PetscBool flg = PETSC_FALSE; 1501 1502 PetscFunctionBegin; 1503 ierr = PetscOptionsHead(PetscOptionsObject,"PFMG options");CHKERRQ(ierr); 1504 ierr = PetscOptionsBool("-pc_pfmg_print_statistics","Print statistics","HYPRE_StructPFMGSetPrintLevel",flg,&flg,NULL);CHKERRQ(ierr); 1505 if (flg) { 1506 int level=3; 1507 PetscStackCallStandard(HYPRE_StructPFMGSetPrintLevel,(ex->hsolver,level)); 1508 } 1509 ierr = PetscOptionsInt("-pc_pfmg_its","Number of iterations of PFMG to use as preconditioner","HYPRE_StructPFMGSetMaxIter",ex->its,&ex->its,NULL);CHKERRQ(ierr); 1510 PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(ex->hsolver,ex->its)); 1511 ierr = PetscOptionsInt("-pc_pfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_StructPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);CHKERRQ(ierr); 1512 PetscStackCallStandard(HYPRE_StructPFMGSetNumPreRelax,(ex->hsolver,ex->num_pre_relax)); 1513 ierr = PetscOptionsInt("-pc_pfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_StructPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);CHKERRQ(ierr); 1514 PetscStackCallStandard(HYPRE_StructPFMGSetNumPostRelax,(ex->hsolver,ex->num_post_relax)); 1515 1516 ierr = PetscOptionsInt("-pc_pfmg_max_levels","Max Levels for MG hierarchy","HYPRE_StructPFMGSetMaxLevels",ex->max_levels,&ex->max_levels,NULL);CHKERRQ(ierr); 1517 PetscStackCallStandard(HYPRE_StructPFMGSetMaxLevels,(ex->hsolver,ex->max_levels)); 1518 1519 ierr = PetscOptionsReal("-pc_pfmg_tol","Tolerance of PFMG","HYPRE_StructPFMGSetTol",ex->tol,&ex->tol,NULL);CHKERRQ(ierr); 1520 PetscStackCallStandard(HYPRE_StructPFMGSetTol,(ex->hsolver,ex->tol)); 1521 ierr = PetscOptionsEList("-pc_pfmg_relax_type","Relax type for the up and down cycles","HYPRE_StructPFMGSetRelaxType",PFMGRelaxType,ALEN(PFMGRelaxType),PFMGRelaxType[ex->relax_type],&ex->relax_type,NULL);CHKERRQ(ierr); 1522 PetscStackCallStandard(HYPRE_StructPFMGSetRelaxType,(ex->hsolver, ex->relax_type)); 1523 ierr = PetscOptionsEList("-pc_pfmg_rap_type","RAP type","HYPRE_StructPFMGSetRAPType",PFMGRAPType,ALEN(PFMGRAPType),PFMGRAPType[ex->rap_type],&ex->rap_type,NULL);CHKERRQ(ierr); 1524 PetscStackCallStandard(HYPRE_StructPFMGSetRAPType,(ex->hsolver, ex->rap_type)); 1525 ierr = PetscOptionsTail();CHKERRQ(ierr); 1526 PetscFunctionReturn(0); 1527 } 1528 1529 #undef __FUNCT__ 1530 #define __FUNCT__ "PCApply_PFMG" 1531 PetscErrorCode PCApply_PFMG(PC pc,Vec x,Vec y) 1532 { 1533 PetscErrorCode ierr; 1534 PC_PFMG *ex = (PC_PFMG*) pc->data; 1535 PetscScalar *yy; 1536 const PetscScalar *xx; 1537 PetscInt ilower[3],iupper[3]; 1538 Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data); 1539 1540 PetscFunctionBegin; 1541 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 1542 ierr = DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr); 1543 iupper[0] += ilower[0] - 1; 1544 iupper[1] += ilower[1] - 1; 1545 iupper[2] += ilower[2] - 1; 1546 1547 /* copy x values over to hypre */ 1548 PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0)); 1549 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1550 PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,(PetscScalar*)xx)); 1551 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1552 PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb)); 1553 PetscStackCallStandard(HYPRE_StructPFMGSolve,(ex->hsolver,mx->hmat,mx->hb,mx->hx)); 1554 1555 /* copy solution values back to PETSc */ 1556 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 1557 PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,yy)); 1558 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 1559 PetscFunctionReturn(0); 1560 } 1561 1562 #undef __FUNCT__ 1563 #define __FUNCT__ "PCApplyRichardson_PFMG" 1564 static PetscErrorCode PCApplyRichardson_PFMG(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason) 1565 { 1566 PC_PFMG *jac = (PC_PFMG*)pc->data; 1567 PetscErrorCode ierr; 1568 PetscInt oits; 1569 1570 PetscFunctionBegin; 1571 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 1572 PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,its*jac->its)); 1573 PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,rtol)); 1574 1575 ierr = PCApply_PFMG(pc,b,y);CHKERRQ(ierr); 1576 PetscStackCallStandard(HYPRE_StructPFMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits)); 1577 *outits = oits; 1578 if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS; 1579 else *reason = PCRICHARDSON_CONVERGED_RTOL; 1580 PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,jac->tol)); 1581 PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,jac->its)); 1582 PetscFunctionReturn(0); 1583 } 1584 1585 1586 #undef __FUNCT__ 1587 #define __FUNCT__ "PCSetUp_PFMG" 1588 PetscErrorCode PCSetUp_PFMG(PC pc) 1589 { 1590 PetscErrorCode ierr; 1591 PC_PFMG *ex = (PC_PFMG*) pc->data; 1592 Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data); 1593 PetscBool flg; 1594 1595 PetscFunctionBegin; 1596 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESTRUCT,&flg);CHKERRQ(ierr); 1597 if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESTRUCT with this preconditioner"); 1598 1599 /* create the hypre solver object and set its information */ 1600 if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver)); 1601 PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver)); 1602 PetscStackCallStandard(HYPRE_StructPFMGSetup,(ex->hsolver,mx->hmat,mx->hb,mx->hx)); 1603 PetscStackCallStandard(HYPRE_StructPFMGSetZeroGuess,(ex->hsolver)); 1604 PetscFunctionReturn(0); 1605 } 1606 1607 1608 /*MC 1609 PCPFMG - the hypre PFMG multigrid solver 1610 1611 Level: advanced 1612 1613 Options Database: 1614 + -pc_pfmg_its <its> number of iterations of PFMG to use as preconditioner 1615 . -pc_pfmg_num_pre_relax <steps> number of smoothing steps before coarse grid 1616 . -pc_pfmg_num_post_relax <steps> number of smoothing steps after coarse grid 1617 . -pc_pfmg_tol <tol> tolerance of PFMG 1618 . -pc_pfmg_relax_type -relaxation type for the up and down cycles, one of Jacobi,Weighted-Jacobi,symmetric-Red/Black-Gauss-Seidel,Red/Black-Gauss-Seidel 1619 - -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin 1620 1621 Notes: This is for CELL-centered descretizations 1622 1623 This must be used with the MATHYPRESTRUCT matrix type. 1624 This is less general than in hypre, it supports only one block per process defined by a PETSc DMDA. 1625 1626 .seealso: PCMG, MATHYPRESTRUCT 1627 M*/ 1628 1629 #undef __FUNCT__ 1630 #define __FUNCT__ "PCCreate_PFMG" 1631 PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc) 1632 { 1633 PetscErrorCode ierr; 1634 PC_PFMG *ex; 1635 1636 PetscFunctionBegin; 1637 ierr = PetscNew(&ex);CHKERRQ(ierr); \ 1638 pc->data = ex; 1639 1640 ex->its = 1; 1641 ex->tol = 1.e-8; 1642 ex->relax_type = 1; 1643 ex->rap_type = 0; 1644 ex->num_pre_relax = 1; 1645 ex->num_post_relax = 1; 1646 ex->max_levels = 0; 1647 1648 pc->ops->setfromoptions = PCSetFromOptions_PFMG; 1649 pc->ops->view = PCView_PFMG; 1650 pc->ops->destroy = PCDestroy_PFMG; 1651 pc->ops->apply = PCApply_PFMG; 1652 pc->ops->applyrichardson = PCApplyRichardson_PFMG; 1653 pc->ops->setup = PCSetUp_PFMG; 1654 1655 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));CHKERRQ(ierr); 1656 PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver)); 1657 PetscFunctionReturn(0); 1658 } 1659 1660 /* ---------------------------------------------------------------------------------------------------------------------------------------------------*/ 1661 1662 /* we know we are working with a HYPRE_SStructMatrix */ 1663 typedef struct { 1664 MPI_Comm hcomm; /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */ 1665 HYPRE_SStructSolver ss_solver; 1666 1667 /* keep copy of SYSPFMG options used so may view them */ 1668 PetscInt its; 1669 double tol; 1670 PetscInt relax_type; 1671 PetscInt num_pre_relax,num_post_relax; 1672 } PC_SysPFMG; 1673 1674 #undef __FUNCT__ 1675 #define __FUNCT__ "PCDestroy_SysPFMG" 1676 PetscErrorCode PCDestroy_SysPFMG(PC pc) 1677 { 1678 PetscErrorCode ierr; 1679 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 1680 1681 PetscFunctionBegin; 1682 if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver)); 1683 ierr = MPI_Comm_free(&ex->hcomm);CHKERRQ(ierr); 1684 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1685 PetscFunctionReturn(0); 1686 } 1687 1688 static const char *SysPFMGRelaxType[] = {"Weighted-Jacobi","Red/Black-Gauss-Seidel"}; 1689 1690 #undef __FUNCT__ 1691 #define __FUNCT__ "PCView_SysPFMG" 1692 PetscErrorCode PCView_SysPFMG(PC pc,PetscViewer viewer) 1693 { 1694 PetscErrorCode ierr; 1695 PetscBool iascii; 1696 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 1697 1698 PetscFunctionBegin; 1699 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1700 if (iascii) { 1701 ierr = PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG preconditioning\n");CHKERRQ(ierr); 1702 ierr = PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: max iterations %d\n",ex->its);CHKERRQ(ierr); 1703 ierr = PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: tolerance %g\n",ex->tol);CHKERRQ(ierr); 1704 ierr = PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);CHKERRQ(ierr); 1705 ierr = PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);CHKERRQ(ierr); 1706 } 1707 PetscFunctionReturn(0); 1708 } 1709 1710 1711 #undef __FUNCT__ 1712 #define __FUNCT__ "PCSetFromOptions_SysPFMG" 1713 PetscErrorCode PCSetFromOptions_SysPFMG(PetscOptions *PetscOptionsObject,PC pc) 1714 { 1715 PetscErrorCode ierr; 1716 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 1717 PetscBool flg = PETSC_FALSE; 1718 1719 PetscFunctionBegin; 1720 ierr = PetscOptionsHead(PetscOptionsObject,"SysPFMG options");CHKERRQ(ierr); 1721 ierr = PetscOptionsBool("-pc_syspfmg_print_statistics","Print statistics","HYPRE_SStructSysPFMGSetPrintLevel",flg,&flg,NULL);CHKERRQ(ierr); 1722 if (flg) { 1723 int level=3; 1724 PetscStackCallStandard(HYPRE_SStructSysPFMGSetPrintLevel,(ex->ss_solver,level)); 1725 } 1726 ierr = PetscOptionsInt("-pc_syspfmg_its","Number of iterations of SysPFMG to use as preconditioner","HYPRE_SStructSysPFMGSetMaxIter",ex->its,&ex->its,NULL);CHKERRQ(ierr); 1727 PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(ex->ss_solver,ex->its)); 1728 ierr = PetscOptionsInt("-pc_syspfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_SStructSysPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);CHKERRQ(ierr); 1729 PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPreRelax,(ex->ss_solver,ex->num_pre_relax)); 1730 ierr = PetscOptionsInt("-pc_syspfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_SStructSysPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);CHKERRQ(ierr); 1731 PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPostRelax,(ex->ss_solver,ex->num_post_relax)); 1732 1733 ierr = PetscOptionsReal("-pc_syspfmg_tol","Tolerance of SysPFMG","HYPRE_SStructSysPFMGSetTol",ex->tol,&ex->tol,NULL);CHKERRQ(ierr); 1734 PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(ex->ss_solver,ex->tol)); 1735 ierr = PetscOptionsEList("-pc_syspfmg_relax_type","Relax type for the up and down cycles","HYPRE_SStructSysPFMGSetRelaxType",SysPFMGRelaxType,4,SysPFMGRelaxType[ex->relax_type],&ex->relax_type,NULL);CHKERRQ(ierr); 1736 PetscStackCallStandard(HYPRE_SStructSysPFMGSetRelaxType,(ex->ss_solver, ex->relax_type)); 1737 ierr = PetscOptionsTail();CHKERRQ(ierr); 1738 PetscFunctionReturn(0); 1739 } 1740 1741 #undef __FUNCT__ 1742 #define __FUNCT__ "PCApply_SysPFMG" 1743 PetscErrorCode PCApply_SysPFMG(PC pc,Vec x,Vec y) 1744 { 1745 PetscErrorCode ierr; 1746 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 1747 PetscScalar *yy; 1748 const PetscScalar *xx; 1749 PetscInt ilower[3],iupper[3]; 1750 Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data); 1751 PetscInt ordering= mx->dofs_order; 1752 PetscInt nvars = mx->nvars; 1753 PetscInt part = 0; 1754 PetscInt size; 1755 PetscInt i; 1756 1757 PetscFunctionBegin; 1758 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 1759 ierr = DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr); 1760 iupper[0] += ilower[0] - 1; 1761 iupper[1] += ilower[1] - 1; 1762 iupper[2] += ilower[2] - 1; 1763 1764 size = 1; 1765 for (i= 0; i< 3; i++) size *= (iupper[i]-ilower[i]+1); 1766 1767 /* copy x values over to hypre for variable ordering */ 1768 if (ordering) { 1769 PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0)); 1770 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1771 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,(PetscScalar*)xx+(size*i))); 1772 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1773 PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b)); 1774 PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x)); 1775 PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x)); 1776 1777 /* copy solution values back to PETSc */ 1778 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 1779 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,yy+(size*i))); 1780 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 1781 } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */ 1782 PetscScalar *z; 1783 PetscInt j, k; 1784 1785 ierr = PetscMalloc1(nvars*size,&z);CHKERRQ(ierr); 1786 PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0)); 1787 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1788 1789 /* transform nodal to hypre's variable ordering for sys_pfmg */ 1790 for (i= 0; i< size; i++) { 1791 k= i*nvars; 1792 for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j]; 1793 } 1794 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i))); 1795 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1796 PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b)); 1797 PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x)); 1798 1799 /* copy solution values back to PETSc */ 1800 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 1801 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i))); 1802 /* transform hypre's variable ordering for sys_pfmg to nodal ordering */ 1803 for (i= 0; i< size; i++) { 1804 k= i*nvars; 1805 for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i]; 1806 } 1807 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 1808 ierr = PetscFree(z);CHKERRQ(ierr); 1809 } 1810 PetscFunctionReturn(0); 1811 } 1812 1813 #undef __FUNCT__ 1814 #define __FUNCT__ "PCApplyRichardson_SysPFMG" 1815 static PetscErrorCode PCApplyRichardson_SysPFMG(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason) 1816 { 1817 PC_SysPFMG *jac = (PC_SysPFMG*)pc->data; 1818 PetscErrorCode ierr; 1819 PetscInt oits; 1820 1821 PetscFunctionBegin; 1822 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 1823 PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,its*jac->its)); 1824 PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,rtol)); 1825 ierr = PCApply_SysPFMG(pc,b,y);CHKERRQ(ierr); 1826 PetscStackCallStandard(HYPRE_SStructSysPFMGGetNumIterations,(jac->ss_solver,(HYPRE_Int *)&oits)); 1827 *outits = oits; 1828 if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS; 1829 else *reason = PCRICHARDSON_CONVERGED_RTOL; 1830 PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,jac->tol)); 1831 PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,jac->its)); 1832 PetscFunctionReturn(0); 1833 } 1834 1835 1836 #undef __FUNCT__ 1837 #define __FUNCT__ "PCSetUp_SysPFMG" 1838 PetscErrorCode PCSetUp_SysPFMG(PC pc) 1839 { 1840 PetscErrorCode ierr; 1841 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 1842 Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data); 1843 PetscBool flg; 1844 1845 PetscFunctionBegin; 1846 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESSTRUCT,&flg);CHKERRQ(ierr); 1847 if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESSTRUCT with this preconditioner"); 1848 1849 /* create the hypre sstruct solver object and set its information */ 1850 if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver)); 1851 PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver)); 1852 PetscStackCallStandard(HYPRE_SStructSysPFMGSetZeroGuess,(ex->ss_solver)); 1853 PetscStackCallStandard(HYPRE_SStructSysPFMGSetup,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x)); 1854 PetscFunctionReturn(0); 1855 } 1856 1857 1858 /*MC 1859 PCSysPFMG - the hypre SysPFMG multigrid solver 1860 1861 Level: advanced 1862 1863 Options Database: 1864 + -pc_syspfmg_its <its> number of iterations of SysPFMG to use as preconditioner 1865 . -pc_syspfmg_num_pre_relax <steps> number of smoothing steps before coarse grid 1866 . -pc_syspfmg_num_post_relax <steps> number of smoothing steps after coarse grid 1867 . -pc_syspfmg_tol <tol> tolerance of SysPFMG 1868 . -pc_syspfmg_relax_type -relaxation type for the up and down cycles, one of Weighted-Jacobi,Red/Black-Gauss-Seidel 1869 1870 Notes: This is for CELL-centered descretizations 1871 1872 This must be used with the MATHYPRESSTRUCT matrix type. 1873 This is less general than in hypre, it supports only one part, and one block per process defined by a PETSc DMDA. 1874 Also, only cell-centered variables. 1875 1876 .seealso: PCMG, MATHYPRESSTRUCT 1877 M*/ 1878 1879 #undef __FUNCT__ 1880 #define __FUNCT__ "PCCreate_SysPFMG" 1881 PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc) 1882 { 1883 PetscErrorCode ierr; 1884 PC_SysPFMG *ex; 1885 1886 PetscFunctionBegin; 1887 ierr = PetscNew(&ex);CHKERRQ(ierr); \ 1888 pc->data = ex; 1889 1890 ex->its = 1; 1891 ex->tol = 1.e-8; 1892 ex->relax_type = 1; 1893 ex->num_pre_relax = 1; 1894 ex->num_post_relax = 1; 1895 1896 pc->ops->setfromoptions = PCSetFromOptions_SysPFMG; 1897 pc->ops->view = PCView_SysPFMG; 1898 pc->ops->destroy = PCDestroy_SysPFMG; 1899 pc->ops->apply = PCApply_SysPFMG; 1900 pc->ops->applyrichardson = PCApplyRichardson_SysPFMG; 1901 pc->ops->setup = PCSetUp_SysPFMG; 1902 1903 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));CHKERRQ(ierr); 1904 PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver)); 1905 PetscFunctionReturn(0); 1906 } 1907