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