1 2 /* 3 Defines the multigrid preconditioner interface. 4 */ 5 #include <petsc/private/pcmgimpl.h> /*I "petscksp.h" I*/ 6 #include <petsc/private/kspimpl.h> 7 #include <petscdm.h> 8 PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC,PetscBool*); 9 10 /* 11 Contains the list of registered coarse space construction routines 12 */ 13 PetscFunctionList PCMGCoarseList = NULL; 14 15 PetscErrorCode PCMGMCycle_Private(PC pc,PC_MG_Levels **mglevelsin,PetscBool transpose,PetscBool matapp,PCRichardsonConvergedReason *reason) 16 { 17 PC_MG *mg = (PC_MG*)pc->data; 18 PC_MG_Levels *mgc,*mglevels = *mglevelsin; 19 PetscErrorCode ierr; 20 PetscInt cycles = (mglevels->level == 1) ? 1 : (PetscInt) mglevels->cycles; 21 22 PetscFunctionBegin; 23 if (mglevels->eventsmoothsolve) {ierr = PetscLogEventBegin(mglevels->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);} 24 if (!transpose) { 25 if (matapp) { 26 ierr = KSPMatSolve(mglevels->smoothd,mglevels->B,mglevels->X);CHKERRQ(ierr); /* pre-smooth */ 27 ierr = KSPCheckSolve(mglevels->smoothd,pc,NULL);CHKERRQ(ierr); 28 } else { 29 ierr = KSPSolve(mglevels->smoothd,mglevels->b,mglevels->x);CHKERRQ(ierr); /* pre-smooth */ 30 ierr = KSPCheckSolve(mglevels->smoothd,pc,mglevels->x);CHKERRQ(ierr); 31 } 32 } else { 33 if (matapp) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Not supported"); 34 ierr = KSPSolveTranspose(mglevels->smoothu,mglevels->b,mglevels->x);CHKERRQ(ierr); /* transpose of post-smooth */ 35 ierr = KSPCheckSolve(mglevels->smoothu,pc,mglevels->x);CHKERRQ(ierr); 36 } 37 if (mglevels->eventsmoothsolve) {ierr = PetscLogEventEnd(mglevels->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);} 38 if (mglevels->level) { /* not the coarsest grid */ 39 if (mglevels->eventresidual) {ierr = PetscLogEventBegin(mglevels->eventresidual,0,0,0,0);CHKERRQ(ierr);} 40 if (matapp && !mglevels->R) { 41 ierr = MatDuplicate(mglevels->B,MAT_DO_NOT_COPY_VALUES,&mglevels->R);CHKERRQ(ierr); 42 } 43 if (!transpose) { 44 if (matapp) { ierr = (*mglevels->matresidual)(mglevels->A,mglevels->B,mglevels->X,mglevels->R);CHKERRQ(ierr); } 45 else { ierr = (*mglevels->residual)(mglevels->A,mglevels->b,mglevels->x,mglevels->r);CHKERRQ(ierr); } 46 } else { 47 if (matapp) { ierr = (*mglevels->matresidualtranspose)(mglevels->A,mglevels->B,mglevels->X,mglevels->R);CHKERRQ(ierr); } 48 else { ierr = (*mglevels->residualtranspose)(mglevels->A,mglevels->b,mglevels->x,mglevels->r);CHKERRQ(ierr); } 49 } 50 if (mglevels->eventresidual) {ierr = PetscLogEventEnd(mglevels->eventresidual,0,0,0,0);CHKERRQ(ierr);} 51 52 /* if on finest level and have convergence criteria set */ 53 if (mglevels->level == mglevels->levels-1 && mg->ttol && reason) { 54 PetscReal rnorm; 55 ierr = VecNorm(mglevels->r,NORM_2,&rnorm);CHKERRQ(ierr); 56 if (rnorm <= mg->ttol) { 57 if (rnorm < mg->abstol) { 58 *reason = PCRICHARDSON_CONVERGED_ATOL; 59 ierr = PetscInfo2(pc,"Linear solver has converged. Residual norm %g is less than absolute tolerance %g\n",(double)rnorm,(double)mg->abstol);CHKERRQ(ierr); 60 } else { 61 *reason = PCRICHARDSON_CONVERGED_RTOL; 62 ierr = PetscInfo2(pc,"Linear solver has converged. Residual norm %g is less than relative tolerance times initial residual norm %g\n",(double)rnorm,(double)mg->ttol);CHKERRQ(ierr); 63 } 64 PetscFunctionReturn(0); 65 } 66 } 67 68 mgc = *(mglevelsin - 1); 69 if (mglevels->eventinterprestrict) {ierr = PetscLogEventBegin(mglevels->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);} 70 if (!transpose) { 71 if (matapp) { ierr = MatMatRestrict(mglevels->restrct,mglevels->R,&mgc->B);CHKERRQ(ierr); } 72 else { ierr = MatRestrict(mglevels->restrct,mglevels->r,mgc->b);CHKERRQ(ierr); } 73 } else { 74 if (matapp) { ierr = MatMatRestrict(mglevels->interpolate,mglevels->R,&mgc->B);CHKERRQ(ierr); } 75 else { ierr = MatRestrict(mglevels->interpolate,mglevels->r,mgc->b);CHKERRQ(ierr); } 76 } 77 if (mglevels->eventinterprestrict) {ierr = PetscLogEventEnd(mglevels->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);} 78 if (matapp) { 79 if (!mgc->X) { 80 ierr = MatDuplicate(mgc->B,MAT_DO_NOT_COPY_VALUES,&mgc->X);CHKERRQ(ierr); 81 } else { 82 ierr = MatZeroEntries(mgc->X);CHKERRQ(ierr); 83 } 84 } else { 85 ierr = VecZeroEntries(mgc->x);CHKERRQ(ierr); 86 } 87 while (cycles--) { 88 ierr = PCMGMCycle_Private(pc,mglevelsin-1,transpose,matapp,reason);CHKERRQ(ierr); 89 } 90 if (mglevels->eventinterprestrict) {ierr = PetscLogEventBegin(mglevels->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);} 91 if (!transpose) { 92 if (matapp) { ierr = MatMatInterpolateAdd(mglevels->interpolate,mgc->X,mglevels->X,&mglevels->X);CHKERRQ(ierr); } 93 else { ierr = MatInterpolateAdd(mglevels->interpolate,mgc->x,mglevels->x,mglevels->x);CHKERRQ(ierr); } 94 } else { 95 ierr = MatInterpolateAdd(mglevels->restrct,mgc->x,mglevels->x,mglevels->x);CHKERRQ(ierr); 96 } 97 if (mglevels->eventinterprestrict) {ierr = PetscLogEventEnd(mglevels->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);} 98 if (mglevels->eventsmoothsolve) {ierr = PetscLogEventBegin(mglevels->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);} 99 if (!transpose) { 100 if (matapp) { 101 ierr = KSPMatSolve(mglevels->smoothu,mglevels->B,mglevels->X);CHKERRQ(ierr); /* post smooth */ 102 ierr = KSPCheckSolve(mglevels->smoothu,pc,NULL);CHKERRQ(ierr); 103 } else { 104 ierr = KSPSolve(mglevels->smoothu,mglevels->b,mglevels->x);CHKERRQ(ierr); /* post smooth */ 105 ierr = KSPCheckSolve(mglevels->smoothu,pc,mglevels->x);CHKERRQ(ierr); 106 } 107 } else { 108 if (matapp) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Not supported"); 109 ierr = KSPSolveTranspose(mglevels->smoothd,mglevels->b,mglevels->x);CHKERRQ(ierr); /* post smooth */ 110 ierr = KSPCheckSolve(mglevels->smoothd,pc,mglevels->x);CHKERRQ(ierr); 111 } 112 if (mglevels->cr) { 113 if (matapp) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Not supported"); 114 /* TODO Turn on copy and turn off noisy if we have an exact solution 115 ierr = VecCopy(mglevels->x, mglevels->crx);CHKERRQ(ierr); 116 ierr = VecCopy(mglevels->b, mglevels->crb);CHKERRQ(ierr); */ 117 ierr = KSPSetNoisy_Private(mglevels->crx);CHKERRQ(ierr); 118 ierr = KSPSolve(mglevels->cr,mglevels->crb,mglevels->crx);CHKERRQ(ierr); /* compatible relaxation */ 119 ierr = KSPCheckSolve(mglevels->cr,pc,mglevels->crx);CHKERRQ(ierr); 120 } 121 if (mglevels->eventsmoothsolve) {ierr = PetscLogEventEnd(mglevels->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);} 122 } 123 PetscFunctionReturn(0); 124 } 125 126 static PetscErrorCode PCApplyRichardson_MG(PC pc,Vec b,Vec x,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool zeroguess,PetscInt *outits,PCRichardsonConvergedReason *reason) 127 { 128 PC_MG *mg = (PC_MG*)pc->data; 129 PC_MG_Levels **mglevels = mg->levels; 130 PetscErrorCode ierr; 131 PC tpc; 132 PetscBool changeu,changed; 133 PetscInt levels = mglevels[0]->levels,i; 134 135 PetscFunctionBegin; 136 /* When the DM is supplying the matrix then it will not exist until here */ 137 for (i=0; i<levels; i++) { 138 if (!mglevels[i]->A) { 139 ierr = KSPGetOperators(mglevels[i]->smoothu,&mglevels[i]->A,NULL);CHKERRQ(ierr); 140 ierr = PetscObjectReference((PetscObject)mglevels[i]->A);CHKERRQ(ierr); 141 } 142 } 143 144 ierr = KSPGetPC(mglevels[levels-1]->smoothd,&tpc);CHKERRQ(ierr); 145 ierr = PCPreSolveChangeRHS(tpc,&changed);CHKERRQ(ierr); 146 ierr = KSPGetPC(mglevels[levels-1]->smoothu,&tpc);CHKERRQ(ierr); 147 ierr = PCPreSolveChangeRHS(tpc,&changeu);CHKERRQ(ierr); 148 if (!changed && !changeu) { 149 ierr = VecDestroy(&mglevels[levels-1]->b);CHKERRQ(ierr); 150 mglevels[levels-1]->b = b; 151 } else { /* if the smoother changes the rhs during PreSolve, we cannot use the input vector */ 152 if (!mglevels[levels-1]->b) { 153 Vec *vec; 154 155 ierr = KSPCreateVecs(mglevels[levels-1]->smoothd,1,&vec,0,NULL);CHKERRQ(ierr); 156 mglevels[levels-1]->b = *vec; 157 ierr = PetscFree(vec);CHKERRQ(ierr); 158 } 159 ierr = VecCopy(b,mglevels[levels-1]->b);CHKERRQ(ierr); 160 } 161 mglevels[levels-1]->x = x; 162 163 mg->rtol = rtol; 164 mg->abstol = abstol; 165 mg->dtol = dtol; 166 if (rtol) { 167 /* compute initial residual norm for relative convergence test */ 168 PetscReal rnorm; 169 if (zeroguess) { 170 ierr = VecNorm(b,NORM_2,&rnorm);CHKERRQ(ierr); 171 } else { 172 ierr = (*mglevels[levels-1]->residual)(mglevels[levels-1]->A,b,x,w);CHKERRQ(ierr); 173 ierr = VecNorm(w,NORM_2,&rnorm);CHKERRQ(ierr); 174 } 175 mg->ttol = PetscMax(rtol*rnorm,abstol); 176 } else if (abstol) mg->ttol = abstol; 177 else mg->ttol = 0.0; 178 179 /* since smoother is applied to full system, not just residual we need to make sure that smoothers don't 180 stop prematurely due to small residual */ 181 for (i=1; i<levels; i++) { 182 ierr = KSPSetTolerances(mglevels[i]->smoothu,0,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 183 if (mglevels[i]->smoothu != mglevels[i]->smoothd) { 184 /* For Richardson the initial guess is nonzero since it is solving in each cycle the original system not just applying as a preconditioner */ 185 ierr = KSPSetInitialGuessNonzero(mglevels[i]->smoothd,PETSC_TRUE);CHKERRQ(ierr); 186 ierr = KSPSetTolerances(mglevels[i]->smoothd,0,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 187 } 188 } 189 190 *reason = (PCRichardsonConvergedReason)0; 191 for (i=0; i<its; i++) { 192 ierr = PCMGMCycle_Private(pc,mglevels+levels-1,PETSC_FALSE,PETSC_FALSE,reason);CHKERRQ(ierr); 193 if (*reason) break; 194 } 195 if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS; 196 *outits = i; 197 if (!changed && !changeu) mglevels[levels-1]->b = NULL; 198 PetscFunctionReturn(0); 199 } 200 201 PetscErrorCode PCReset_MG(PC pc) 202 { 203 PC_MG *mg = (PC_MG*)pc->data; 204 PC_MG_Levels **mglevels = mg->levels; 205 PetscErrorCode ierr; 206 PetscInt i,c,n; 207 208 PetscFunctionBegin; 209 if (mglevels) { 210 n = mglevels[0]->levels; 211 for (i=0; i<n-1; i++) { 212 ierr = VecDestroy(&mglevels[i+1]->r);CHKERRQ(ierr); 213 ierr = VecDestroy(&mglevels[i]->b);CHKERRQ(ierr); 214 ierr = VecDestroy(&mglevels[i]->x);CHKERRQ(ierr); 215 ierr = MatDestroy(&mglevels[i+1]->R);CHKERRQ(ierr); 216 ierr = MatDestroy(&mglevels[i]->B);CHKERRQ(ierr); 217 ierr = MatDestroy(&mglevels[i]->X);CHKERRQ(ierr); 218 ierr = VecDestroy(&mglevels[i]->crx);CHKERRQ(ierr); 219 ierr = VecDestroy(&mglevels[i]->crb);CHKERRQ(ierr); 220 ierr = MatDestroy(&mglevels[i+1]->restrct);CHKERRQ(ierr); 221 ierr = MatDestroy(&mglevels[i+1]->interpolate);CHKERRQ(ierr); 222 ierr = MatDestroy(&mglevels[i+1]->inject);CHKERRQ(ierr); 223 ierr = VecDestroy(&mglevels[i+1]->rscale);CHKERRQ(ierr); 224 } 225 ierr = VecDestroy(&mglevels[n-1]->crx);CHKERRQ(ierr); 226 ierr = VecDestroy(&mglevels[n-1]->crb);CHKERRQ(ierr); 227 /* this is not null only if the smoother on the finest level 228 changes the rhs during PreSolve */ 229 ierr = VecDestroy(&mglevels[n-1]->b);CHKERRQ(ierr); 230 ierr = MatDestroy(&mglevels[n-1]->B);CHKERRQ(ierr); 231 232 for (i=0; i<n; i++) { 233 if (mglevels[i]->coarseSpace) for (c = 0; c < mg->Nc; ++c) {ierr = VecDestroy(&mglevels[i]->coarseSpace[c]);CHKERRQ(ierr);} 234 ierr = PetscFree(mglevels[i]->coarseSpace);CHKERRQ(ierr); 235 mglevels[i]->coarseSpace = NULL; 236 ierr = MatDestroy(&mglevels[i]->A);CHKERRQ(ierr); 237 if (mglevels[i]->smoothd != mglevels[i]->smoothu) { 238 ierr = KSPReset(mglevels[i]->smoothd);CHKERRQ(ierr); 239 } 240 ierr = KSPReset(mglevels[i]->smoothu);CHKERRQ(ierr); 241 if (mglevels[i]->cr) {ierr = KSPReset(mglevels[i]->cr);CHKERRQ(ierr);} 242 } 243 mg->Nc = 0; 244 } 245 PetscFunctionReturn(0); 246 } 247 248 /* Implementing CR 249 250 We only want to make corrections that ``do not change'' the coarse solution. What we mean by not changing is that if I prolong my coarse solution to the fine grid and then inject that fine solution back to the coarse grid, I get the same answer. Injection is what Brannick calls R. We want the complementary projector to Inj, which we will call S, after Brannick, so that Inj S = 0. Now the orthogonal projector onto the range of Inj^T is 251 252 Inj^T (Inj Inj^T)^{-1} Inj 253 254 and if Inj is a VecScatter, as it is now in PETSc, we have 255 256 Inj^T Inj 257 258 and 259 260 S = I - Inj^T Inj 261 262 since 263 264 Inj S = Inj - (Inj Inj^T) Inj = 0. 265 266 Brannick suggests 267 268 A \to S^T A S \qquad\mathrm{and}\qquad M \to S^T M S 269 270 but I do not think his :math:`S^T S = I` is correct. Our S is an orthogonal projector, so :math:`S^T S = S^2 = S`. We will use 271 272 M^{-1} A \to S M^{-1} A S 273 274 In fact, since it is somewhat hard in PETSc to do the symmetric application, we will just apply S on the left. 275 276 Check: || Inj P - I ||_F < tol 277 Check: In general, Inj Inj^T = I 278 */ 279 280 typedef struct { 281 PC mg; /* The PCMG object */ 282 PetscInt l; /* The multigrid level for this solver */ 283 Mat Inj; /* The injection matrix */ 284 Mat S; /* I - Inj^T Inj */ 285 } CRContext; 286 287 static PetscErrorCode CRSetup_Private(PC pc) 288 { 289 CRContext *ctx; 290 Mat It; 291 PetscErrorCode ierr; 292 293 PetscFunctionBeginUser; 294 ierr = PCShellGetContext(pc, (void **) &ctx);CHKERRQ(ierr); 295 ierr = PCMGGetInjection(ctx->mg, ctx->l, &It);CHKERRQ(ierr); 296 if (!It) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONGSTATE, "CR requires that injection be defined for this PCMG"); 297 ierr = MatCreateTranspose(It, &ctx->Inj);CHKERRQ(ierr); 298 ierr = MatCreateNormal(ctx->Inj, &ctx->S);CHKERRQ(ierr); 299 ierr = MatScale(ctx->S, -1.0);CHKERRQ(ierr); 300 ierr = MatShift(ctx->S, 1.0);CHKERRQ(ierr); 301 PetscFunctionReturn(0); 302 } 303 304 static PetscErrorCode CRApply_Private(PC pc, Vec x, Vec y) 305 { 306 CRContext *ctx; 307 PetscErrorCode ierr; 308 309 PetscFunctionBeginUser; 310 ierr = PCShellGetContext(pc, (void **) &ctx);CHKERRQ(ierr); 311 ierr = MatMult(ctx->S, x, y);CHKERRQ(ierr); 312 PetscFunctionReturn(0); 313 } 314 315 static PetscErrorCode CRDestroy_Private(PC pc) 316 { 317 CRContext *ctx; 318 PetscErrorCode ierr; 319 320 PetscFunctionBeginUser; 321 ierr = PCShellGetContext(pc, (void **) &ctx);CHKERRQ(ierr); 322 ierr = MatDestroy(&ctx->Inj);CHKERRQ(ierr); 323 ierr = MatDestroy(&ctx->S);CHKERRQ(ierr); 324 ierr = PetscFree(ctx);CHKERRQ(ierr); 325 ierr = PCShellSetContext(pc, NULL);CHKERRQ(ierr); 326 PetscFunctionReturn(0); 327 } 328 329 static PetscErrorCode CreateCR_Private(PC pc, PetscInt l, PC *cr) 330 { 331 CRContext *ctx; 332 PetscErrorCode ierr; 333 334 PetscFunctionBeginUser; 335 ierr = PCCreate(PetscObjectComm((PetscObject) pc), cr);CHKERRQ(ierr); 336 ierr = PetscObjectSetName((PetscObject) *cr, "S (complementary projector to injection)");CHKERRQ(ierr); 337 ierr = PetscCalloc1(1, &ctx);CHKERRQ(ierr); 338 ctx->mg = pc; 339 ctx->l = l; 340 ierr = PCSetType(*cr, PCSHELL);CHKERRQ(ierr); 341 ierr = PCShellSetContext(*cr, ctx);CHKERRQ(ierr); 342 ierr = PCShellSetApply(*cr, CRApply_Private);CHKERRQ(ierr); 343 ierr = PCShellSetSetUp(*cr, CRSetup_Private);CHKERRQ(ierr); 344 ierr = PCShellSetDestroy(*cr, CRDestroy_Private);CHKERRQ(ierr); 345 PetscFunctionReturn(0); 346 } 347 348 PetscErrorCode PCMGSetLevels_MG(PC pc,PetscInt levels,MPI_Comm *comms) 349 { 350 PetscErrorCode ierr; 351 PC_MG *mg = (PC_MG*)pc->data; 352 MPI_Comm comm; 353 PC_MG_Levels **mglevels = mg->levels; 354 PCMGType mgtype = mg->am; 355 PetscInt mgctype = (PetscInt) PC_MG_CYCLE_V; 356 PetscInt i; 357 PetscMPIInt size; 358 const char *prefix; 359 PC ipc; 360 PetscInt n; 361 362 PetscFunctionBegin; 363 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 364 PetscValidLogicalCollectiveInt(pc,levels,2); 365 if (mg->nlevels == levels) PetscFunctionReturn(0); 366 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 367 if (mglevels) { 368 mgctype = mglevels[0]->cycles; 369 /* changing the number of levels so free up the previous stuff */ 370 ierr = PCReset_MG(pc);CHKERRQ(ierr); 371 n = mglevels[0]->levels; 372 for (i=0; i<n; i++) { 373 if (mglevels[i]->smoothd != mglevels[i]->smoothu) { 374 ierr = KSPDestroy(&mglevels[i]->smoothd);CHKERRQ(ierr); 375 } 376 ierr = KSPDestroy(&mglevels[i]->smoothu);CHKERRQ(ierr); 377 ierr = KSPDestroy(&mglevels[i]->cr);CHKERRQ(ierr); 378 ierr = PetscFree(mglevels[i]);CHKERRQ(ierr); 379 } 380 ierr = PetscFree(mg->levels);CHKERRQ(ierr); 381 } 382 383 mg->nlevels = levels; 384 385 ierr = PetscMalloc1(levels,&mglevels);CHKERRQ(ierr); 386 ierr = PetscLogObjectMemory((PetscObject)pc,levels*(sizeof(PC_MG*)));CHKERRQ(ierr); 387 388 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr); 389 390 mg->stageApply = 0; 391 for (i=0; i<levels; i++) { 392 ierr = PetscNewLog(pc,&mglevels[i]);CHKERRQ(ierr); 393 394 mglevels[i]->level = i; 395 mglevels[i]->levels = levels; 396 mglevels[i]->cycles = mgctype; 397 mg->default_smoothu = 2; 398 mg->default_smoothd = 2; 399 mglevels[i]->eventsmoothsetup = 0; 400 mglevels[i]->eventsmoothsolve = 0; 401 mglevels[i]->eventresidual = 0; 402 mglevels[i]->eventinterprestrict = 0; 403 404 if (comms) comm = comms[i]; 405 if (comm != MPI_COMM_NULL) { 406 ierr = KSPCreate(comm,&mglevels[i]->smoothd);CHKERRQ(ierr); 407 ierr = KSPSetErrorIfNotConverged(mglevels[i]->smoothd,pc->erroriffailure);CHKERRQ(ierr); 408 ierr = PetscObjectIncrementTabLevel((PetscObject)mglevels[i]->smoothd,(PetscObject)pc,levels-i);CHKERRQ(ierr); 409 ierr = KSPSetOptionsPrefix(mglevels[i]->smoothd,prefix);CHKERRQ(ierr); 410 ierr = PetscObjectComposedDataSetInt((PetscObject) mglevels[i]->smoothd, PetscMGLevelId, mglevels[i]->level);CHKERRQ(ierr); 411 if (i || levels == 1) { 412 char tprefix[128]; 413 414 ierr = KSPSetType(mglevels[i]->smoothd,KSPCHEBYSHEV);CHKERRQ(ierr); 415 ierr = KSPSetConvergenceTest(mglevels[i]->smoothd,KSPConvergedSkip,NULL,NULL);CHKERRQ(ierr); 416 ierr = KSPSetNormType(mglevels[i]->smoothd,KSP_NORM_NONE);CHKERRQ(ierr); 417 ierr = KSPGetPC(mglevels[i]->smoothd,&ipc);CHKERRQ(ierr); 418 ierr = PCSetType(ipc,PCSOR);CHKERRQ(ierr); 419 ierr = KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT, mg->default_smoothd);CHKERRQ(ierr); 420 421 ierr = PetscSNPrintf(tprefix,128,"mg_levels_%d_",(int)i);CHKERRQ(ierr); 422 ierr = KSPAppendOptionsPrefix(mglevels[i]->smoothd,tprefix);CHKERRQ(ierr); 423 } else { 424 ierr = KSPAppendOptionsPrefix(mglevels[0]->smoothd,"mg_coarse_");CHKERRQ(ierr); 425 426 /* coarse solve is (redundant) LU by default; set shifttype NONZERO to avoid annoying zero-pivot in LU preconditioner */ 427 ierr = KSPSetType(mglevels[0]->smoothd,KSPPREONLY);CHKERRQ(ierr); 428 ierr = KSPGetPC(mglevels[0]->smoothd,&ipc);CHKERRQ(ierr); 429 ierr = MPI_Comm_size(comm,&size);CHKERRMPI(ierr); 430 if (size > 1) { 431 ierr = PCSetType(ipc,PCREDUNDANT);CHKERRQ(ierr); 432 } else { 433 ierr = PCSetType(ipc,PCLU);CHKERRQ(ierr); 434 } 435 ierr = PCFactorSetShiftType(ipc,MAT_SHIFT_INBLOCKS);CHKERRQ(ierr); 436 } 437 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)mglevels[i]->smoothd);CHKERRQ(ierr); 438 } 439 mglevels[i]->smoothu = mglevels[i]->smoothd; 440 mg->rtol = 0.0; 441 mg->abstol = 0.0; 442 mg->dtol = 0.0; 443 mg->ttol = 0.0; 444 mg->cyclesperpcapply = 1; 445 } 446 mg->levels = mglevels; 447 ierr = PCMGSetType(pc,mgtype);CHKERRQ(ierr); 448 PetscFunctionReturn(0); 449 } 450 451 /*@C 452 PCMGSetLevels - Sets the number of levels to use with MG. 453 Must be called before any other MG routine. 454 455 Logically Collective on PC 456 457 Input Parameters: 458 + pc - the preconditioner context 459 . levels - the number of levels 460 - comms - optional communicators for each level; this is to allow solving the coarser problems 461 on smaller sets of processes. For processes that are not included in the computation 462 you must pass MPI_COMM_NULL. Use comms = NULL to specify that all processes 463 should participate in each level of problem. 464 465 Level: intermediate 466 467 Notes: 468 If the number of levels is one then the multigrid uses the -mg_levels prefix 469 for setting the level options rather than the -mg_coarse prefix. 470 471 You can free the information in comms after this routine is called. 472 473 The array of MPI communicators must contain MPI_COMM_NULL for those ranks that at each level 474 are not participating in the coarser solve. For example, with 2 levels and 1 and 2 ranks on 475 the two levels, rank 0 in the original communicator will pass in an array of 2 communicators 476 of size 2 and 1, while rank 1 in the original communicator will pass in array of 2 communicators 477 the first of size 2 and the second of value MPI_COMM_NULL since the rank 1 does not participate 478 in the coarse grid solve. 479 480 Since each coarser level may have a new MPI_Comm with fewer ranks than the previous, one 481 must take special care in providing the restriction and interpolation operation. We recommend 482 providing these as two step operations; first perform a standard restriction or interpolation on 483 the full number of ranks for that level and then use an MPI call to copy the resulting vector 484 array entries (after calls to VecGetArray()) to the smaller or larger number of ranks, note in both 485 cases the MPI calls must be made on the larger of the two communicators. Traditional MPI send and 486 recieves or MPI_AlltoAllv() could be used to do the reshuffling of the vector entries. 487 488 Fortran Notes: 489 Use comms = PETSC_NULL_MPI_COMM as the equivalent of NULL in the C interface. Note PETSC_NULL_MPI_COMM 490 is not MPI_COMM_NULL. It is more like PETSC_NULL_INTEGER, PETSC_NULL_REAL etc. 491 492 .seealso: PCMGSetType(), PCMGGetLevels() 493 @*/ 494 PetscErrorCode PCMGSetLevels(PC pc,PetscInt levels,MPI_Comm *comms) 495 { 496 PetscErrorCode ierr; 497 498 PetscFunctionBegin; 499 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 500 if (comms) PetscValidPointer(comms,3); 501 ierr = PetscTryMethod(pc,"PCMGSetLevels_C",(PC,PetscInt,MPI_Comm*),(pc,levels,comms));CHKERRQ(ierr); 502 PetscFunctionReturn(0); 503 } 504 505 506 PetscErrorCode PCDestroy_MG(PC pc) 507 { 508 PetscErrorCode ierr; 509 PC_MG *mg = (PC_MG*)pc->data; 510 PC_MG_Levels **mglevels = mg->levels; 511 PetscInt i,n; 512 513 PetscFunctionBegin; 514 ierr = PCReset_MG(pc);CHKERRQ(ierr); 515 if (mglevels) { 516 n = mglevels[0]->levels; 517 for (i=0; i<n; i++) { 518 if (mglevels[i]->smoothd != mglevels[i]->smoothu) { 519 ierr = KSPDestroy(&mglevels[i]->smoothd);CHKERRQ(ierr); 520 } 521 ierr = KSPDestroy(&mglevels[i]->smoothu);CHKERRQ(ierr); 522 ierr = KSPDestroy(&mglevels[i]->cr);CHKERRQ(ierr); 523 ierr = PetscFree(mglevels[i]);CHKERRQ(ierr); 524 } 525 ierr = PetscFree(mg->levels);CHKERRQ(ierr); 526 } 527 ierr = PetscFree(pc->data);CHKERRQ(ierr); 528 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGetInterpolations_C",NULL);CHKERRQ(ierr); 529 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGetCoarseOperators_C",NULL);CHKERRQ(ierr); 530 PetscFunctionReturn(0); 531 } 532 533 534 /* 535 PCApply_MG - Runs either an additive, multiplicative, Kaskadic 536 or full cycle of multigrid. 537 538 Note: 539 A simple wrapper which calls PCMGMCycle(),PCMGACycle(), or PCMGFCycle(). 540 */ 541 static PetscErrorCode PCApply_MG_Internal(PC pc,Vec b,Vec x,Mat B,Mat X,PetscBool transpose) 542 { 543 PC_MG *mg = (PC_MG*)pc->data; 544 PC_MG_Levels **mglevels = mg->levels; 545 PetscErrorCode ierr; 546 PC tpc; 547 PetscInt levels = mglevels[0]->levels,i; 548 PetscBool changeu,changed,matapp; 549 550 PetscFunctionBegin; 551 matapp = (PetscBool)(B && X); 552 if (mg->stageApply) {ierr = PetscLogStagePush(mg->stageApply);CHKERRQ(ierr);} 553 /* When the DM is supplying the matrix then it will not exist until here */ 554 for (i=0; i<levels; i++) { 555 if (!mglevels[i]->A) { 556 ierr = KSPGetOperators(mglevels[i]->smoothu,&mglevels[i]->A,NULL);CHKERRQ(ierr); 557 ierr = PetscObjectReference((PetscObject)mglevels[i]->A);CHKERRQ(ierr); 558 } 559 } 560 561 ierr = KSPGetPC(mglevels[levels-1]->smoothd,&tpc);CHKERRQ(ierr); 562 ierr = PCPreSolveChangeRHS(tpc,&changed);CHKERRQ(ierr); 563 ierr = KSPGetPC(mglevels[levels-1]->smoothu,&tpc);CHKERRQ(ierr); 564 ierr = PCPreSolveChangeRHS(tpc,&changeu);CHKERRQ(ierr); 565 if (!changeu && !changed) { 566 if (matapp) { 567 ierr = MatDestroy(&mglevels[levels-1]->B);CHKERRQ(ierr); 568 mglevels[levels-1]->B = B; 569 } else { 570 ierr = VecDestroy(&mglevels[levels-1]->b);CHKERRQ(ierr); 571 mglevels[levels-1]->b = b; 572 } 573 } else { /* if the smoother changes the rhs during PreSolve, we cannot use the input vector */ 574 if (matapp) { 575 if (mglevels[levels-1]->B) { 576 PetscInt N1,N2; 577 PetscBool flg; 578 579 ierr = MatGetSize(mglevels[levels-1]->B,NULL,&N1);CHKERRQ(ierr); 580 ierr = MatGetSize(B,NULL,&N2);CHKERRQ(ierr); 581 ierr = PetscObjectTypeCompare((PetscObject)mglevels[levels-1]->B,((PetscObject)B)->type_name,&flg);CHKERRQ(ierr); 582 if (N1 != N2 || !flg) { 583 ierr = MatDestroy(&mglevels[levels-1]->B);CHKERRQ(ierr); 584 } 585 } 586 if (!mglevels[levels-1]->B) { 587 ierr = MatDuplicate(B,MAT_COPY_VALUES,&mglevels[levels-1]->B);CHKERRQ(ierr); 588 } else { 589 ierr = MatCopy(B,mglevels[levels-1]->B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 590 } 591 } else { 592 if (!mglevels[levels-1]->b) { 593 Vec *vec; 594 595 ierr = KSPCreateVecs(mglevels[levels-1]->smoothd,1,&vec,0,NULL);CHKERRQ(ierr); 596 mglevels[levels-1]->b = *vec; 597 ierr = PetscFree(vec);CHKERRQ(ierr); 598 } 599 ierr = VecCopy(b,mglevels[levels-1]->b);CHKERRQ(ierr); 600 } 601 } 602 if (matapp) { mglevels[levels-1]->X = X; } 603 else { mglevels[levels-1]->x = x; } 604 605 /* If coarser Xs are present, it means we have already block applied the PC at least once 606 Reset operators if sizes/type do no match */ 607 if (matapp && levels > 1 && mglevels[levels-2]->X) { 608 PetscInt Xc,Bc; 609 PetscBool flg; 610 611 ierr = MatGetSize(mglevels[levels-2]->X,NULL,&Xc);CHKERRQ(ierr); 612 ierr = MatGetSize(mglevels[levels-1]->B,NULL,&Bc);CHKERRQ(ierr); 613 ierr = PetscObjectTypeCompare((PetscObject)mglevels[levels-2]->X,((PetscObject)mglevels[levels-1]->X)->type_name,&flg);CHKERRQ(ierr); 614 if (Xc != Bc || !flg) { 615 ierr = MatDestroy(&mglevels[levels-1]->R);CHKERRQ(ierr); 616 for (i=0;i<levels-1;i++) { 617 ierr = MatDestroy(&mglevels[i]->R);CHKERRQ(ierr); 618 ierr = MatDestroy(&mglevels[i]->B);CHKERRQ(ierr); 619 ierr = MatDestroy(&mglevels[i]->X);CHKERRQ(ierr); 620 } 621 } 622 } 623 624 if (mg->am == PC_MG_MULTIPLICATIVE) { 625 if (matapp) { ierr = MatZeroEntries(X);CHKERRQ(ierr); } 626 else { ierr = VecZeroEntries(x);CHKERRQ(ierr); } 627 for (i=0; i<mg->cyclesperpcapply; i++) { 628 ierr = PCMGMCycle_Private(pc,mglevels+levels-1,transpose,matapp,NULL);CHKERRQ(ierr); 629 } 630 } else if (mg->am == PC_MG_ADDITIVE) { 631 ierr = PCMGACycle_Private(pc,mglevels,transpose,matapp);CHKERRQ(ierr); 632 } else if (mg->am == PC_MG_KASKADE) { 633 ierr = PCMGKCycle_Private(pc,mglevels,transpose,matapp);CHKERRQ(ierr); 634 } else { 635 ierr = PCMGFCycle_Private(pc,mglevels,transpose,matapp);CHKERRQ(ierr); 636 } 637 if (mg->stageApply) {ierr = PetscLogStagePop();CHKERRQ(ierr);} 638 if (!changeu && !changed) { 639 if (matapp) { mglevels[levels-1]->B = NULL; } 640 else { mglevels[levels-1]->b = NULL; } 641 } 642 PetscFunctionReturn(0); 643 } 644 645 static PetscErrorCode PCApply_MG(PC pc,Vec b,Vec x) 646 { 647 PetscErrorCode ierr; 648 649 PetscFunctionBegin; 650 ierr = PCApply_MG_Internal(pc,b,x,NULL,NULL,PETSC_FALSE);CHKERRQ(ierr); 651 PetscFunctionReturn(0); 652 } 653 654 static PetscErrorCode PCApplyTranspose_MG(PC pc,Vec b,Vec x) 655 { 656 PetscErrorCode ierr; 657 658 PetscFunctionBegin; 659 ierr = PCApply_MG_Internal(pc,b,x,NULL,NULL,PETSC_TRUE);CHKERRQ(ierr); 660 PetscFunctionReturn(0); 661 } 662 663 static PetscErrorCode PCMatApply_MG(PC pc,Mat b,Mat x) 664 { 665 PetscErrorCode ierr; 666 667 PetscFunctionBegin; 668 ierr = PCApply_MG_Internal(pc,NULL,NULL,b,x,PETSC_FALSE);CHKERRQ(ierr); 669 PetscFunctionReturn(0); 670 } 671 672 PetscErrorCode PCSetFromOptions_MG(PetscOptionItems *PetscOptionsObject,PC pc) 673 { 674 PetscErrorCode ierr; 675 PetscInt levels,cycles; 676 PetscBool flg, flg2; 677 PC_MG *mg = (PC_MG*)pc->data; 678 PC_MG_Levels **mglevels; 679 PCMGType mgtype; 680 PCMGCycleType mgctype; 681 PCMGGalerkinType gtype; 682 683 PetscFunctionBegin; 684 levels = PetscMax(mg->nlevels,1); 685 ierr = PetscOptionsHead(PetscOptionsObject,"Multigrid options");CHKERRQ(ierr); 686 ierr = PetscOptionsInt("-pc_mg_levels","Number of Levels","PCMGSetLevels",levels,&levels,&flg);CHKERRQ(ierr); 687 if (!flg && !mg->levels && pc->dm) { 688 ierr = DMGetRefineLevel(pc->dm,&levels);CHKERRQ(ierr); 689 levels++; 690 mg->usedmfornumberoflevels = PETSC_TRUE; 691 } 692 ierr = PCMGSetLevels(pc,levels,NULL);CHKERRQ(ierr); 693 mglevels = mg->levels; 694 695 mgctype = (PCMGCycleType) mglevels[0]->cycles; 696 ierr = PetscOptionsEnum("-pc_mg_cycle_type","V cycle or for W-cycle","PCMGSetCycleType",PCMGCycleTypes,(PetscEnum)mgctype,(PetscEnum*)&mgctype,&flg);CHKERRQ(ierr); 697 if (flg) { 698 ierr = PCMGSetCycleType(pc,mgctype);CHKERRQ(ierr); 699 } 700 gtype = mg->galerkin; 701 ierr = PetscOptionsEnum("-pc_mg_galerkin","Use Galerkin process to compute coarser operators","PCMGSetGalerkin",PCMGGalerkinTypes,(PetscEnum)gtype,(PetscEnum*)>ype,&flg);CHKERRQ(ierr); 702 if (flg) { 703 ierr = PCMGSetGalerkin(pc,gtype);CHKERRQ(ierr); 704 } 705 flg2 = PETSC_FALSE; 706 ierr = PetscOptionsBool("-pc_mg_adapt_interp","Adapt interpolation using some coarse space","PCMGSetAdaptInterpolation",PETSC_FALSE,&flg2,&flg);CHKERRQ(ierr); 707 if (flg) {ierr = PCMGSetAdaptInterpolation(pc, flg2);CHKERRQ(ierr);} 708 ierr = PetscOptionsInt("-pc_mg_adapt_interp_n","Size of the coarse space for adaptive interpolation","PCMGSetCoarseSpace",mg->Nc,&mg->Nc,&flg);CHKERRQ(ierr); 709 ierr = PetscOptionsEnum("-pc_mg_adapt_interp_coarse_space","Type of coarse space: polynomial, harmonic, eigenvector, generalized_eigenvector","PCMGSetAdaptCoarseSpaceType",PCMGCoarseSpaceTypes,(PetscEnum)mg->coarseSpaceType,(PetscEnum*)&mg->coarseSpaceType,&flg);CHKERRQ(ierr); 710 ierr = PetscOptionsBool("-pc_mg_mesp_monitor","Monitor the multilevel eigensolver","PCMGSetAdaptInterpolation",PETSC_FALSE,&mg->mespMonitor,&flg);CHKERRQ(ierr); 711 flg2 = PETSC_FALSE; 712 ierr = PetscOptionsBool("-pc_mg_adapt_cr","Monitor coarse space quality using Compatible Relaxation (CR)","PCMGSetAdaptCR",PETSC_FALSE,&flg2,&flg);CHKERRQ(ierr); 713 if (flg) {ierr = PCMGSetAdaptCR(pc, flg2);CHKERRQ(ierr);} 714 flg = PETSC_FALSE; 715 ierr = PetscOptionsBool("-pc_mg_distinct_smoothup","Create separate smoothup KSP and append the prefix _up","PCMGSetDistinctSmoothUp",PETSC_FALSE,&flg,NULL);CHKERRQ(ierr); 716 if (flg) { 717 ierr = PCMGSetDistinctSmoothUp(pc);CHKERRQ(ierr); 718 } 719 mgtype = mg->am; 720 ierr = PetscOptionsEnum("-pc_mg_type","Multigrid type","PCMGSetType",PCMGTypes,(PetscEnum)mgtype,(PetscEnum*)&mgtype,&flg);CHKERRQ(ierr); 721 if (flg) { 722 ierr = PCMGSetType(pc,mgtype);CHKERRQ(ierr); 723 } 724 if (mg->am == PC_MG_MULTIPLICATIVE) { 725 ierr = PetscOptionsInt("-pc_mg_multiplicative_cycles","Number of cycles for each preconditioner step","PCMGMultiplicativeSetCycles",mg->cyclesperpcapply,&cycles,&flg);CHKERRQ(ierr); 726 if (flg) { 727 ierr = PCMGMultiplicativeSetCycles(pc,cycles);CHKERRQ(ierr); 728 } 729 } 730 flg = PETSC_FALSE; 731 ierr = PetscOptionsBool("-pc_mg_log","Log times for each multigrid level","None",flg,&flg,NULL);CHKERRQ(ierr); 732 if (flg) { 733 PetscInt i; 734 char eventname[128]; 735 736 levels = mglevels[0]->levels; 737 for (i=0; i<levels; i++) { 738 sprintf(eventname,"MGSetup Level %d",(int)i); 739 ierr = PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventsmoothsetup);CHKERRQ(ierr); 740 sprintf(eventname,"MGSmooth Level %d",(int)i); 741 ierr = PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventsmoothsolve);CHKERRQ(ierr); 742 if (i) { 743 sprintf(eventname,"MGResid Level %d",(int)i); 744 ierr = PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventresidual);CHKERRQ(ierr); 745 sprintf(eventname,"MGInterp Level %d",(int)i); 746 ierr = PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventinterprestrict);CHKERRQ(ierr); 747 } 748 } 749 750 #if defined(PETSC_USE_LOG) 751 { 752 const char *sname = "MG Apply"; 753 PetscStageLog stageLog; 754 PetscInt st; 755 756 ierr = PetscLogGetStageLog(&stageLog);CHKERRQ(ierr); 757 for (st = 0; st < stageLog->numStages; ++st) { 758 PetscBool same; 759 760 ierr = PetscStrcmp(stageLog->stageInfo[st].name, sname, &same);CHKERRQ(ierr); 761 if (same) mg->stageApply = st; 762 } 763 if (!mg->stageApply) { 764 ierr = PetscLogStageRegister(sname, &mg->stageApply);CHKERRQ(ierr); 765 } 766 } 767 #endif 768 } 769 ierr = PetscOptionsTail();CHKERRQ(ierr); 770 /* Check option consistency */ 771 ierr = PCMGGetGalerkin(pc, >ype);CHKERRQ(ierr); 772 ierr = PCMGGetAdaptInterpolation(pc, &flg);CHKERRQ(ierr); 773 if (flg && (gtype >= PC_MG_GALERKIN_NONE)) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_INCOMP, "Must use Galerkin coarse operators when adapting the interpolator"); 774 PetscFunctionReturn(0); 775 } 776 777 const char *const PCMGTypes[] = {"MULTIPLICATIVE","ADDITIVE","FULL","KASKADE","PCMGType","PC_MG",NULL}; 778 const char *const PCMGCycleTypes[] = {"invalid","v","w","PCMGCycleType","PC_MG_CYCLE",NULL}; 779 const char *const PCMGGalerkinTypes[] = {"both","pmat","mat","none","external","PCMGGalerkinType","PC_MG_GALERKIN",NULL}; 780 const char *const PCMGCoarseSpaceTypes[] = {"polynomial","harmonic","eigenvector","generalized_eigenvector","PCMGCoarseSpaceType","PCMG_POLYNOMIAL",NULL}; 781 782 #include <petscdraw.h> 783 PetscErrorCode PCView_MG(PC pc,PetscViewer viewer) 784 { 785 PC_MG *mg = (PC_MG*)pc->data; 786 PC_MG_Levels **mglevels = mg->levels; 787 PetscErrorCode ierr; 788 PetscInt levels = mglevels ? mglevels[0]->levels : 0,i; 789 PetscBool iascii,isbinary,isdraw; 790 791 PetscFunctionBegin; 792 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 793 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 794 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 795 if (iascii) { 796 const char *cyclename = levels ? (mglevels[0]->cycles == PC_MG_CYCLE_V ? "v" : "w") : "unknown"; 797 ierr = PetscViewerASCIIPrintf(viewer," type is %s, levels=%D cycles=%s\n", PCMGTypes[mg->am],levels,cyclename);CHKERRQ(ierr); 798 if (mg->am == PC_MG_MULTIPLICATIVE) { 799 ierr = PetscViewerASCIIPrintf(viewer," Cycles per PCApply=%d\n",mg->cyclesperpcapply);CHKERRQ(ierr); 800 } 801 if (mg->galerkin == PC_MG_GALERKIN_BOTH) { 802 ierr = PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid matrices\n");CHKERRQ(ierr); 803 } else if (mg->galerkin == PC_MG_GALERKIN_PMAT) { 804 ierr = PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid matrices for pmat\n");CHKERRQ(ierr); 805 } else if (mg->galerkin == PC_MG_GALERKIN_MAT) { 806 ierr = PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid matrices for mat\n");CHKERRQ(ierr); 807 } else if (mg->galerkin == PC_MG_GALERKIN_EXTERNAL) { 808 ierr = PetscViewerASCIIPrintf(viewer," Using externally compute Galerkin coarse grid matrices\n");CHKERRQ(ierr); 809 } else { 810 ierr = PetscViewerASCIIPrintf(viewer," Not using Galerkin computed coarse grid matrices\n");CHKERRQ(ierr); 811 } 812 if (mg->view){ 813 ierr = (*mg->view)(pc,viewer);CHKERRQ(ierr); 814 } 815 for (i=0; i<levels; i++) { 816 if (!i) { 817 ierr = PetscViewerASCIIPrintf(viewer,"Coarse grid solver -- level -------------------------------\n",i);CHKERRQ(ierr); 818 } else { 819 ierr = PetscViewerASCIIPrintf(viewer,"Down solver (pre-smoother) on level %D -------------------------------\n",i);CHKERRQ(ierr); 820 } 821 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 822 ierr = KSPView(mglevels[i]->smoothd,viewer);CHKERRQ(ierr); 823 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 824 if (i && mglevels[i]->smoothd == mglevels[i]->smoothu) { 825 ierr = PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) same as down solver (pre-smoother)\n");CHKERRQ(ierr); 826 } else if (i) { 827 ierr = PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) on level %D -------------------------------\n",i);CHKERRQ(ierr); 828 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 829 ierr = KSPView(mglevels[i]->smoothu,viewer);CHKERRQ(ierr); 830 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 831 } 832 if (i && mglevels[i]->cr) { 833 ierr = PetscViewerASCIIPrintf(viewer,"CR solver on level %D -------------------------------\n",i);CHKERRQ(ierr); 834 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 835 ierr = KSPView(mglevels[i]->cr,viewer);CHKERRQ(ierr); 836 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 837 } 838 } 839 } else if (isbinary) { 840 for (i=levels-1; i>=0; i--) { 841 ierr = KSPView(mglevels[i]->smoothd,viewer);CHKERRQ(ierr); 842 if (i && mglevels[i]->smoothd != mglevels[i]->smoothu) { 843 ierr = KSPView(mglevels[i]->smoothu,viewer);CHKERRQ(ierr); 844 } 845 } 846 } else if (isdraw) { 847 PetscDraw draw; 848 PetscReal x,w,y,bottom,th; 849 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 850 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr); 851 ierr = PetscDrawStringGetSize(draw,NULL,&th);CHKERRQ(ierr); 852 bottom = y - th; 853 for (i=levels-1; i>=0; i--) { 854 if (!mglevels[i]->smoothu || (mglevels[i]->smoothu == mglevels[i]->smoothd)) { 855 ierr = PetscDrawPushCurrentPoint(draw,x,bottom);CHKERRQ(ierr); 856 ierr = KSPView(mglevels[i]->smoothd,viewer);CHKERRQ(ierr); 857 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 858 } else { 859 w = 0.5*PetscMin(1.0-x,x); 860 ierr = PetscDrawPushCurrentPoint(draw,x+w,bottom);CHKERRQ(ierr); 861 ierr = KSPView(mglevels[i]->smoothd,viewer);CHKERRQ(ierr); 862 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 863 ierr = PetscDrawPushCurrentPoint(draw,x-w,bottom);CHKERRQ(ierr); 864 ierr = KSPView(mglevels[i]->smoothu,viewer);CHKERRQ(ierr); 865 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 866 } 867 ierr = PetscDrawGetBoundingBox(draw,NULL,&bottom,NULL,NULL);CHKERRQ(ierr); 868 bottom -= th; 869 } 870 } 871 PetscFunctionReturn(0); 872 } 873 874 #include <petsc/private/kspimpl.h> 875 876 /* 877 Calls setup for the KSP on each level 878 */ 879 PetscErrorCode PCSetUp_MG(PC pc) 880 { 881 PC_MG *mg = (PC_MG*)pc->data; 882 PC_MG_Levels **mglevels = mg->levels; 883 PetscErrorCode ierr; 884 PetscInt i,n; 885 PC cpc; 886 PetscBool dump = PETSC_FALSE,opsset,use_amat,missinginterpolate = PETSC_FALSE; 887 Mat dA,dB; 888 Vec tvec; 889 DM *dms; 890 PetscViewer viewer = NULL; 891 PetscBool dAeqdB = PETSC_FALSE, needRestricts = PETSC_FALSE, doCR = PETSC_FALSE; 892 893 PetscFunctionBegin; 894 if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels with PCMGSetLevels() before setting up"); 895 n = mglevels[0]->levels; 896 /* FIX: Move this to PCSetFromOptions_MG? */ 897 if (mg->usedmfornumberoflevels) { 898 PetscInt levels; 899 ierr = DMGetRefineLevel(pc->dm,&levels);CHKERRQ(ierr); 900 levels++; 901 if (levels > n) { /* the problem is now being solved on a finer grid */ 902 ierr = PCMGSetLevels(pc,levels,NULL);CHKERRQ(ierr); 903 n = levels; 904 ierr = PCSetFromOptions(pc);CHKERRQ(ierr); /* it is bad to call this here, but otherwise will never be called for the new hierarchy */ 905 mglevels = mg->levels; 906 } 907 } 908 ierr = KSPGetPC(mglevels[0]->smoothd,&cpc);CHKERRQ(ierr); 909 910 911 /* If user did not provide fine grid operators OR operator was not updated since last global KSPSetOperators() */ 912 /* so use those from global PC */ 913 /* Is this what we always want? What if user wants to keep old one? */ 914 ierr = KSPGetOperatorsSet(mglevels[n-1]->smoothd,NULL,&opsset);CHKERRQ(ierr); 915 if (opsset) { 916 Mat mmat; 917 ierr = KSPGetOperators(mglevels[n-1]->smoothd,NULL,&mmat);CHKERRQ(ierr); 918 if (mmat == pc->pmat) opsset = PETSC_FALSE; 919 } 920 921 /* Create CR solvers */ 922 ierr = PCMGGetAdaptCR(pc, &doCR);CHKERRQ(ierr); 923 if (doCR) { 924 const char *prefix; 925 926 ierr = PCGetOptionsPrefix(pc, &prefix);CHKERRQ(ierr); 927 for (i = 1; i < n; ++i) { 928 PC ipc, cr; 929 char crprefix[128]; 930 931 ierr = KSPCreate(PetscObjectComm((PetscObject) pc), &mglevels[i]->cr);CHKERRQ(ierr); 932 ierr = KSPSetErrorIfNotConverged(mglevels[i]->cr, PETSC_FALSE);CHKERRQ(ierr); 933 ierr = PetscObjectIncrementTabLevel((PetscObject) mglevels[i]->cr, (PetscObject) pc, n-i);CHKERRQ(ierr); 934 ierr = KSPSetOptionsPrefix(mglevels[i]->cr, prefix);CHKERRQ(ierr); 935 ierr = PetscObjectComposedDataSetInt((PetscObject) mglevels[i]->cr, PetscMGLevelId, mglevels[i]->level);CHKERRQ(ierr); 936 ierr = KSPSetType(mglevels[i]->cr, KSPCHEBYSHEV);CHKERRQ(ierr); 937 ierr = KSPSetConvergenceTest(mglevels[i]->cr, KSPConvergedSkip, NULL, NULL);CHKERRQ(ierr); 938 ierr = KSPSetNormType(mglevels[i]->cr, KSP_NORM_PRECONDITIONED);CHKERRQ(ierr); 939 ierr = KSPGetPC(mglevels[i]->cr, &ipc);CHKERRQ(ierr); 940 941 ierr = PCSetType(ipc, PCCOMPOSITE);CHKERRQ(ierr); 942 ierr = PCCompositeSetType(ipc, PC_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr); 943 ierr = PCCompositeAddPCType(ipc, PCSOR);CHKERRQ(ierr); 944 ierr = CreateCR_Private(pc, i, &cr);CHKERRQ(ierr); 945 ierr = PCCompositeAddPC(ipc, cr);CHKERRQ(ierr); 946 ierr = PCDestroy(&cr);CHKERRQ(ierr); 947 948 ierr = KSPSetTolerances(mglevels[i]->cr, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, mg->default_smoothd);CHKERRQ(ierr); 949 ierr = KSPSetInitialGuessNonzero(mglevels[i]->cr, PETSC_TRUE);CHKERRQ(ierr); 950 ierr = PetscSNPrintf(crprefix, 128, "mg_levels_%d_cr_", (int) i);CHKERRQ(ierr); 951 ierr = KSPAppendOptionsPrefix(mglevels[i]->cr, crprefix);CHKERRQ(ierr); 952 ierr = PetscLogObjectParent((PetscObject) pc, (PetscObject) mglevels[i]->cr);CHKERRQ(ierr); 953 } 954 } 955 956 if (!opsset) { 957 ierr = PCGetUseAmat(pc,&use_amat);CHKERRQ(ierr); 958 if (use_amat) { 959 ierr = PetscInfo(pc,"Using outer operators to define finest grid operator \n because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");CHKERRQ(ierr); 960 ierr = KSPSetOperators(mglevels[n-1]->smoothd,pc->mat,pc->pmat);CHKERRQ(ierr); 961 } else { 962 ierr = PetscInfo(pc,"Using matrix (pmat) operators to define finest grid operator \n because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");CHKERRQ(ierr); 963 ierr = KSPSetOperators(mglevels[n-1]->smoothd,pc->pmat,pc->pmat);CHKERRQ(ierr); 964 } 965 } 966 967 for (i=n-1; i>0; i--) { 968 if (!(mglevels[i]->interpolate || mglevels[i]->restrct)) { 969 missinginterpolate = PETSC_TRUE; 970 continue; 971 } 972 } 973 974 ierr = KSPGetOperators(mglevels[n-1]->smoothd,&dA,&dB);CHKERRQ(ierr); 975 if (dA == dB) dAeqdB = PETSC_TRUE; 976 if ((mg->galerkin == PC_MG_GALERKIN_NONE) || (((mg->galerkin == PC_MG_GALERKIN_PMAT) || (mg->galerkin == PC_MG_GALERKIN_MAT)) && !dAeqdB)) { 977 needRestricts = PETSC_TRUE; /* user must compute either mat, pmat, or both so must restrict x to coarser levels */ 978 } 979 980 981 /* 982 Skipping if user has provided all interpolation/restriction needed (since DM might not be able to produce them (when coming from SNES/TS) 983 Skipping for galerkin==2 (externally managed hierarchy such as ML and GAMG). Cleaner logic here would be great. Wrap ML/GAMG as DMs? 984 */ 985 if (missinginterpolate && pc->dm && mg->galerkin != PC_MG_GALERKIN_EXTERNAL && !pc->setupcalled) { 986 /* construct the interpolation from the DMs */ 987 Mat p; 988 Vec rscale; 989 ierr = PetscMalloc1(n,&dms);CHKERRQ(ierr); 990 dms[n-1] = pc->dm; 991 /* Separately create them so we do not get DMKSP interference between levels */ 992 for (i=n-2; i>-1; i--) {ierr = DMCoarsen(dms[i+1],MPI_COMM_NULL,&dms[i]);CHKERRQ(ierr);} 993 /* 994 Force the mat type of coarse level operator to be AIJ because usually we want to use LU for coarse level. 995 Notice that it can be overwritten by -mat_type because KSPSetUp() reads command line options. 996 But it is safe to use -dm_mat_type. 997 998 The mat type should not be hardcoded like this, we need to find a better way. 999 ierr = DMSetMatType(dms[0],MATAIJ);CHKERRQ(ierr); 1000 */ 1001 for (i=n-2; i>-1; i--) { 1002 DMKSP kdm; 1003 PetscBool dmhasrestrict, dmhasinject; 1004 ierr = KSPSetDM(mglevels[i]->smoothd,dms[i]);CHKERRQ(ierr); 1005 if (!needRestricts) {ierr = KSPSetDMActive(mglevels[i]->smoothd,PETSC_FALSE);CHKERRQ(ierr);} 1006 if (mglevels[i]->smoothd != mglevels[i]->smoothu) { 1007 ierr = KSPSetDM(mglevels[i]->smoothu,dms[i]);CHKERRQ(ierr); 1008 if (!needRestricts) {ierr = KSPSetDMActive(mglevels[i]->smoothu,PETSC_FALSE);CHKERRQ(ierr);} 1009 } 1010 if (mglevels[i]->cr) { 1011 ierr = KSPSetDM(mglevels[i]->cr,dms[i]);CHKERRQ(ierr); 1012 if (!needRestricts) {ierr = KSPSetDMActive(mglevels[i]->cr,PETSC_FALSE);CHKERRQ(ierr);} 1013 } 1014 ierr = DMGetDMKSPWrite(dms[i],&kdm);CHKERRQ(ierr); 1015 /* Ugly hack so that the next KSPSetUp() will use the RHS that we set. A better fix is to change dmActive to take 1016 * a bitwise OR of computing the matrix, RHS, and initial iterate. */ 1017 kdm->ops->computerhs = NULL; 1018 kdm->rhsctx = NULL; 1019 if (!mglevels[i+1]->interpolate) { 1020 ierr = DMCreateInterpolation(dms[i],dms[i+1],&p,&rscale);CHKERRQ(ierr); 1021 ierr = PCMGSetInterpolation(pc,i+1,p);CHKERRQ(ierr); 1022 if (rscale) {ierr = PCMGSetRScale(pc,i+1,rscale);CHKERRQ(ierr);} 1023 ierr = VecDestroy(&rscale);CHKERRQ(ierr); 1024 ierr = MatDestroy(&p);CHKERRQ(ierr); 1025 } 1026 ierr = DMHasCreateRestriction(dms[i],&dmhasrestrict);CHKERRQ(ierr); 1027 if (dmhasrestrict && !mglevels[i+1]->restrct){ 1028 ierr = DMCreateRestriction(dms[i],dms[i+1],&p);CHKERRQ(ierr); 1029 ierr = PCMGSetRestriction(pc,i+1,p);CHKERRQ(ierr); 1030 ierr = MatDestroy(&p);CHKERRQ(ierr); 1031 } 1032 ierr = DMHasCreateInjection(dms[i],&dmhasinject);CHKERRQ(ierr); 1033 if (dmhasinject && !mglevels[i+1]->inject){ 1034 ierr = DMCreateInjection(dms[i],dms[i+1],&p);CHKERRQ(ierr); 1035 ierr = PCMGSetInjection(pc,i+1,p);CHKERRQ(ierr); 1036 ierr = MatDestroy(&p);CHKERRQ(ierr); 1037 } 1038 } 1039 1040 for (i=n-2; i>-1; i--) {ierr = DMDestroy(&dms[i]);CHKERRQ(ierr);} 1041 ierr = PetscFree(dms);CHKERRQ(ierr); 1042 } 1043 1044 if (pc->dm && !pc->setupcalled) { 1045 /* finest smoother also gets DM but it is not active, independent of whether galerkin==PC_MG_GALERKIN_EXTERNAL */ 1046 ierr = KSPSetDM(mglevels[n-1]->smoothd,pc->dm);CHKERRQ(ierr); 1047 ierr = KSPSetDMActive(mglevels[n-1]->smoothd,PETSC_FALSE);CHKERRQ(ierr); 1048 if (mglevels[n-1]->smoothd != mglevels[n-1]->smoothu) { 1049 ierr = KSPSetDM(mglevels[n-1]->smoothu,pc->dm);CHKERRQ(ierr); 1050 ierr = KSPSetDMActive(mglevels[n-1]->smoothu,PETSC_FALSE);CHKERRQ(ierr); 1051 } 1052 if (mglevels[n-1]->cr) { 1053 ierr = KSPSetDM(mglevels[n-1]->cr,pc->dm);CHKERRQ(ierr); 1054 ierr = KSPSetDMActive(mglevels[n-1]->cr,PETSC_FALSE);CHKERRQ(ierr); 1055 } 1056 } 1057 1058 if (mg->galerkin < PC_MG_GALERKIN_NONE) { 1059 Mat A,B; 1060 PetscBool doA = PETSC_FALSE,doB = PETSC_FALSE; 1061 MatReuse reuse = MAT_INITIAL_MATRIX; 1062 1063 if ((mg->galerkin == PC_MG_GALERKIN_PMAT) || (mg->galerkin == PC_MG_GALERKIN_BOTH)) doB = PETSC_TRUE; 1064 if ((mg->galerkin == PC_MG_GALERKIN_MAT) || ((mg->galerkin == PC_MG_GALERKIN_BOTH) && (dA != dB))) doA = PETSC_TRUE; 1065 if (pc->setupcalled) reuse = MAT_REUSE_MATRIX; 1066 for (i=n-2; i>-1; i--) { 1067 if (!mglevels[i+1]->restrct && !mglevels[i+1]->interpolate) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must provide interpolation or restriction for each MG level except level 0"); 1068 if (!mglevels[i+1]->interpolate) { 1069 ierr = PCMGSetInterpolation(pc,i+1,mglevels[i+1]->restrct);CHKERRQ(ierr); 1070 } 1071 if (!mglevels[i+1]->restrct) { 1072 ierr = PCMGSetRestriction(pc,i+1,mglevels[i+1]->interpolate);CHKERRQ(ierr); 1073 } 1074 if (reuse == MAT_REUSE_MATRIX) { 1075 ierr = KSPGetOperators(mglevels[i]->smoothd,&A,&B);CHKERRQ(ierr); 1076 } 1077 if (doA) { 1078 ierr = MatGalerkin(mglevels[i+1]->restrct,dA,mglevels[i+1]->interpolate,reuse,1.0,&A);CHKERRQ(ierr); 1079 } 1080 if (doB) { 1081 ierr = MatGalerkin(mglevels[i+1]->restrct,dB,mglevels[i+1]->interpolate,reuse,1.0,&B);CHKERRQ(ierr); 1082 } 1083 /* the management of the PetscObjectReference() and PetscObjecDereference() below is rather delicate */ 1084 if (!doA && dAeqdB) { 1085 if (reuse == MAT_INITIAL_MATRIX) {ierr = PetscObjectReference((PetscObject)B);CHKERRQ(ierr);} 1086 A = B; 1087 } else if (!doA && reuse == MAT_INITIAL_MATRIX) { 1088 ierr = KSPGetOperators(mglevels[i]->smoothd,&A,NULL);CHKERRQ(ierr); 1089 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 1090 } 1091 if (!doB && dAeqdB) { 1092 if (reuse == MAT_INITIAL_MATRIX) {ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);} 1093 B = A; 1094 } else if (!doB && reuse == MAT_INITIAL_MATRIX) { 1095 ierr = KSPGetOperators(mglevels[i]->smoothd,NULL,&B);CHKERRQ(ierr); 1096 ierr = PetscObjectReference((PetscObject)B);CHKERRQ(ierr); 1097 } 1098 if (reuse == MAT_INITIAL_MATRIX) { 1099 ierr = KSPSetOperators(mglevels[i]->smoothd,A,B);CHKERRQ(ierr); 1100 ierr = PetscObjectDereference((PetscObject)A);CHKERRQ(ierr); 1101 ierr = PetscObjectDereference((PetscObject)B);CHKERRQ(ierr); 1102 } 1103 dA = A; 1104 dB = B; 1105 } 1106 } 1107 1108 1109 /* Adapt interpolation matrices */ 1110 if (mg->adaptInterpolation) { 1111 mg->Nc = mg->Nc < 0 ? 6 : mg->Nc; /* Default to 6 modes */ 1112 for (i = 0; i < n; ++i) { 1113 ierr = PCMGComputeCoarseSpace_Internal(pc, i, mg->coarseSpaceType, mg->Nc, !i ? NULL : mglevels[i-1]->coarseSpace, &mglevels[i]->coarseSpace);CHKERRQ(ierr); 1114 if (i) {ierr = PCMGAdaptInterpolator_Internal(pc, i, mglevels[i-1]->smoothu, mglevels[i]->smoothu, mg->Nc, mglevels[i-1]->coarseSpace, mglevels[i]->coarseSpace);CHKERRQ(ierr);} 1115 } 1116 for (i = n-2; i > -1; --i) { 1117 ierr = PCMGRecomputeLevelOperators_Internal(pc, i);CHKERRQ(ierr); 1118 } 1119 } 1120 1121 if (needRestricts && pc->dm) { 1122 for (i=n-2; i>=0; i--) { 1123 DM dmfine,dmcoarse; 1124 Mat Restrict,Inject; 1125 Vec rscale; 1126 ierr = KSPGetDM(mglevels[i+1]->smoothd,&dmfine);CHKERRQ(ierr); 1127 ierr = KSPGetDM(mglevels[i]->smoothd,&dmcoarse);CHKERRQ(ierr); 1128 ierr = PCMGGetRestriction(pc,i+1,&Restrict);CHKERRQ(ierr); 1129 ierr = PCMGGetRScale(pc,i+1,&rscale);CHKERRQ(ierr); 1130 ierr = PCMGGetInjection(pc,i+1,&Inject);CHKERRQ(ierr); 1131 ierr = DMRestrict(dmfine,Restrict,rscale,Inject,dmcoarse);CHKERRQ(ierr); 1132 } 1133 } 1134 1135 if (!pc->setupcalled) { 1136 for (i=0; i<n; i++) { 1137 ierr = KSPSetFromOptions(mglevels[i]->smoothd);CHKERRQ(ierr); 1138 } 1139 for (i=1; i<n; i++) { 1140 if (mglevels[i]->smoothu && (mglevels[i]->smoothu != mglevels[i]->smoothd)) { 1141 ierr = KSPSetFromOptions(mglevels[i]->smoothu);CHKERRQ(ierr); 1142 } 1143 if (mglevels[i]->cr) { 1144 ierr = KSPSetFromOptions(mglevels[i]->cr);CHKERRQ(ierr); 1145 } 1146 } 1147 /* insure that if either interpolation or restriction is set the other other one is set */ 1148 for (i=1; i<n; i++) { 1149 ierr = PCMGGetInterpolation(pc,i,NULL);CHKERRQ(ierr); 1150 ierr = PCMGGetRestriction(pc,i,NULL);CHKERRQ(ierr); 1151 } 1152 for (i=0; i<n-1; i++) { 1153 if (!mglevels[i]->b) { 1154 Vec *vec; 1155 ierr = KSPCreateVecs(mglevels[i]->smoothd,1,&vec,0,NULL);CHKERRQ(ierr); 1156 ierr = PCMGSetRhs(pc,i,*vec);CHKERRQ(ierr); 1157 ierr = VecDestroy(vec);CHKERRQ(ierr); 1158 ierr = PetscFree(vec);CHKERRQ(ierr); 1159 } 1160 if (!mglevels[i]->r && i) { 1161 ierr = VecDuplicate(mglevels[i]->b,&tvec);CHKERRQ(ierr); 1162 ierr = PCMGSetR(pc,i,tvec);CHKERRQ(ierr); 1163 ierr = VecDestroy(&tvec);CHKERRQ(ierr); 1164 } 1165 if (!mglevels[i]->x) { 1166 ierr = VecDuplicate(mglevels[i]->b,&tvec);CHKERRQ(ierr); 1167 ierr = PCMGSetX(pc,i,tvec);CHKERRQ(ierr); 1168 ierr = VecDestroy(&tvec);CHKERRQ(ierr); 1169 } 1170 if (doCR) { 1171 ierr = VecDuplicate(mglevels[i]->b,&mglevels[i]->crx);CHKERRQ(ierr); 1172 ierr = VecDuplicate(mglevels[i]->b,&mglevels[i]->crb);CHKERRQ(ierr); 1173 ierr = VecZeroEntries(mglevels[i]->crb);CHKERRQ(ierr); 1174 } 1175 } 1176 if (n != 1 && !mglevels[n-1]->r) { 1177 /* PCMGSetR() on the finest level if user did not supply it */ 1178 Vec *vec; 1179 ierr = KSPCreateVecs(mglevels[n-1]->smoothd,1,&vec,0,NULL);CHKERRQ(ierr); 1180 ierr = PCMGSetR(pc,n-1,*vec);CHKERRQ(ierr); 1181 ierr = VecDestroy(vec);CHKERRQ(ierr); 1182 ierr = PetscFree(vec);CHKERRQ(ierr); 1183 } 1184 if (doCR) { 1185 ierr = VecDuplicate(mglevels[n-1]->r, &mglevels[n-1]->crx);CHKERRQ(ierr); 1186 ierr = VecDuplicate(mglevels[n-1]->r, &mglevels[n-1]->crb);CHKERRQ(ierr); 1187 ierr = VecZeroEntries(mglevels[n-1]->crb);CHKERRQ(ierr); 1188 } 1189 } 1190 1191 if (pc->dm) { 1192 /* need to tell all the coarser levels to rebuild the matrix using the DM for that level */ 1193 for (i=0; i<n-1; i++) { 1194 if (mglevels[i]->smoothd->setupstage != KSP_SETUP_NEW) mglevels[i]->smoothd->setupstage = KSP_SETUP_NEWMATRIX; 1195 } 1196 } 1197 1198 for (i=1; i<n; i++) { 1199 if (mglevels[i]->smoothu == mglevels[i]->smoothd || mg->am == PC_MG_FULL || mg->am == PC_MG_KASKADE || mg->cyclesperpcapply > 1){ 1200 /* if doing only down then initial guess is zero */ 1201 ierr = KSPSetInitialGuessNonzero(mglevels[i]->smoothd,PETSC_TRUE);CHKERRQ(ierr); 1202 } 1203 if (mglevels[i]->cr) {ierr = KSPSetInitialGuessNonzero(mglevels[i]->cr,PETSC_TRUE);CHKERRQ(ierr);} 1204 if (mglevels[i]->eventsmoothsetup) {ierr = PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);} 1205 ierr = KSPSetUp(mglevels[i]->smoothd);CHKERRQ(ierr); 1206 if (mglevels[i]->smoothd->reason == KSP_DIVERGED_PC_FAILED) { 1207 pc->failedreason = PC_SUBPC_ERROR; 1208 } 1209 if (mglevels[i]->eventsmoothsetup) {ierr = PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);} 1210 if (!mglevels[i]->residual) { 1211 Mat mat; 1212 ierr = KSPGetOperators(mglevels[i]->smoothd,&mat,NULL);CHKERRQ(ierr); 1213 ierr = PCMGSetResidual(pc,i,PCMGResidualDefault,mat);CHKERRQ(ierr); 1214 } 1215 if (!mglevels[i]->residualtranspose) { 1216 Mat mat; 1217 ierr = KSPGetOperators(mglevels[i]->smoothd,&mat,NULL);CHKERRQ(ierr); 1218 ierr = PCMGSetResidualTranspose(pc,i,PCMGResidualTransposeDefault,mat);CHKERRQ(ierr); 1219 } 1220 } 1221 for (i=1; i<n; i++) { 1222 if (mglevels[i]->smoothu && mglevels[i]->smoothu != mglevels[i]->smoothd) { 1223 Mat downmat,downpmat; 1224 1225 /* check if operators have been set for up, if not use down operators to set them */ 1226 ierr = KSPGetOperatorsSet(mglevels[i]->smoothu,&opsset,NULL);CHKERRQ(ierr); 1227 if (!opsset) { 1228 ierr = KSPGetOperators(mglevels[i]->smoothd,&downmat,&downpmat);CHKERRQ(ierr); 1229 ierr = KSPSetOperators(mglevels[i]->smoothu,downmat,downpmat);CHKERRQ(ierr); 1230 } 1231 1232 ierr = KSPSetInitialGuessNonzero(mglevels[i]->smoothu,PETSC_TRUE);CHKERRQ(ierr); 1233 if (mglevels[i]->eventsmoothsetup) {ierr = PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);} 1234 ierr = KSPSetUp(mglevels[i]->smoothu);CHKERRQ(ierr); 1235 if (mglevels[i]->smoothu->reason == KSP_DIVERGED_PC_FAILED) { 1236 pc->failedreason = PC_SUBPC_ERROR; 1237 } 1238 if (mglevels[i]->eventsmoothsetup) {ierr = PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);} 1239 } 1240 if (mglevels[i]->cr) { 1241 Mat downmat,downpmat; 1242 1243 /* check if operators have been set for up, if not use down operators to set them */ 1244 ierr = KSPGetOperatorsSet(mglevels[i]->cr,&opsset,NULL);CHKERRQ(ierr); 1245 if (!opsset) { 1246 ierr = KSPGetOperators(mglevels[i]->smoothd,&downmat,&downpmat);CHKERRQ(ierr); 1247 ierr = KSPSetOperators(mglevels[i]->cr,downmat,downpmat);CHKERRQ(ierr); 1248 } 1249 1250 ierr = KSPSetInitialGuessNonzero(mglevels[i]->cr,PETSC_TRUE);CHKERRQ(ierr); 1251 if (mglevels[i]->eventsmoothsetup) {ierr = PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);} 1252 ierr = KSPSetUp(mglevels[i]->cr);CHKERRQ(ierr); 1253 if (mglevels[i]->cr->reason == KSP_DIVERGED_PC_FAILED) { 1254 pc->failedreason = PC_SUBPC_ERROR; 1255 } 1256 if (mglevels[i]->eventsmoothsetup) {ierr = PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);} 1257 } 1258 } 1259 1260 if (mglevels[0]->eventsmoothsetup) {ierr = PetscLogEventBegin(mglevels[0]->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);} 1261 ierr = KSPSetUp(mglevels[0]->smoothd);CHKERRQ(ierr); 1262 if (mglevels[0]->smoothd->reason == KSP_DIVERGED_PC_FAILED) { 1263 pc->failedreason = PC_SUBPC_ERROR; 1264 } 1265 if (mglevels[0]->eventsmoothsetup) {ierr = PetscLogEventEnd(mglevels[0]->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);} 1266 1267 /* 1268 Dump the interpolation/restriction matrices plus the 1269 Jacobian/stiffness on each level. This allows MATLAB users to 1270 easily check if the Galerkin condition A_c = R A_f R^T is satisfied. 1271 1272 Only support one or the other at the same time. 1273 */ 1274 #if defined(PETSC_USE_SOCKET_VIEWER) 1275 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_mg_dump_matlab",&dump,NULL);CHKERRQ(ierr); 1276 if (dump) viewer = PETSC_VIEWER_SOCKET_(PetscObjectComm((PetscObject)pc)); 1277 dump = PETSC_FALSE; 1278 #endif 1279 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_mg_dump_binary",&dump,NULL);CHKERRQ(ierr); 1280 if (dump) viewer = PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)pc)); 1281 1282 if (viewer) { 1283 for (i=1; i<n; i++) { 1284 ierr = MatView(mglevels[i]->restrct,viewer);CHKERRQ(ierr); 1285 } 1286 for (i=0; i<n; i++) { 1287 ierr = KSPGetPC(mglevels[i]->smoothd,&pc);CHKERRQ(ierr); 1288 ierr = MatView(pc->mat,viewer);CHKERRQ(ierr); 1289 } 1290 } 1291 PetscFunctionReturn(0); 1292 } 1293 1294 /* -------------------------------------------------------------------------------------*/ 1295 1296 PetscErrorCode PCMGGetLevels_MG(PC pc, PetscInt *levels) 1297 { 1298 PC_MG *mg = (PC_MG *) pc->data; 1299 1300 PetscFunctionBegin; 1301 *levels = mg->nlevels; 1302 PetscFunctionReturn(0); 1303 } 1304 1305 /*@ 1306 PCMGGetLevels - Gets the number of levels to use with MG. 1307 1308 Not Collective 1309 1310 Input Parameter: 1311 . pc - the preconditioner context 1312 1313 Output parameter: 1314 . levels - the number of levels 1315 1316 Level: advanced 1317 1318 .seealso: PCMGSetLevels() 1319 @*/ 1320 PetscErrorCode PCMGGetLevels(PC pc,PetscInt *levels) 1321 { 1322 PetscErrorCode ierr; 1323 1324 PetscFunctionBegin; 1325 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1326 PetscValidIntPointer(levels,2); 1327 *levels = 0; 1328 ierr = PetscTryMethod(pc,"PCMGGetLevels_C",(PC,PetscInt*),(pc,levels));CHKERRQ(ierr); 1329 PetscFunctionReturn(0); 1330 } 1331 1332 /*@ 1333 PCMGSetType - Determines the form of multigrid to use: 1334 multiplicative, additive, full, or the Kaskade algorithm. 1335 1336 Logically Collective on PC 1337 1338 Input Parameters: 1339 + pc - the preconditioner context 1340 - form - multigrid form, one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE, 1341 PC_MG_FULL, PC_MG_KASKADE 1342 1343 Options Database Key: 1344 . -pc_mg_type <form> - Sets <form>, one of multiplicative, 1345 additive, full, kaskade 1346 1347 Level: advanced 1348 1349 .seealso: PCMGSetLevels() 1350 @*/ 1351 PetscErrorCode PCMGSetType(PC pc,PCMGType form) 1352 { 1353 PC_MG *mg = (PC_MG*)pc->data; 1354 1355 PetscFunctionBegin; 1356 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1357 PetscValidLogicalCollectiveEnum(pc,form,2); 1358 mg->am = form; 1359 if (form == PC_MG_MULTIPLICATIVE) pc->ops->applyrichardson = PCApplyRichardson_MG; 1360 else pc->ops->applyrichardson = NULL; 1361 PetscFunctionReturn(0); 1362 } 1363 1364 /*@ 1365 PCMGGetType - Determines the form of multigrid to use: 1366 multiplicative, additive, full, or the Kaskade algorithm. 1367 1368 Logically Collective on PC 1369 1370 Input Parameter: 1371 . pc - the preconditioner context 1372 1373 Output Parameter: 1374 . type - one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,PC_MG_FULL, PC_MG_KASKADE 1375 1376 1377 Level: advanced 1378 1379 .seealso: PCMGSetLevels() 1380 @*/ 1381 PetscErrorCode PCMGGetType(PC pc,PCMGType *type) 1382 { 1383 PC_MG *mg = (PC_MG*)pc->data; 1384 1385 PetscFunctionBegin; 1386 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1387 *type = mg->am; 1388 PetscFunctionReturn(0); 1389 } 1390 1391 /*@ 1392 PCMGSetCycleType - Sets the type cycles to use. Use PCMGSetCycleTypeOnLevel() for more 1393 complicated cycling. 1394 1395 Logically Collective on PC 1396 1397 Input Parameters: 1398 + pc - the multigrid context 1399 - n - either PC_MG_CYCLE_V or PC_MG_CYCLE_W 1400 1401 Options Database Key: 1402 . -pc_mg_cycle_type <v,w> - provide the cycle desired 1403 1404 Level: advanced 1405 1406 .seealso: PCMGSetCycleTypeOnLevel() 1407 @*/ 1408 PetscErrorCode PCMGSetCycleType(PC pc,PCMGCycleType n) 1409 { 1410 PC_MG *mg = (PC_MG*)pc->data; 1411 PC_MG_Levels **mglevels = mg->levels; 1412 PetscInt i,levels; 1413 1414 PetscFunctionBegin; 1415 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1416 PetscValidLogicalCollectiveEnum(pc,n,2); 1417 if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Must set MG levels with PCMGSetLevels() before calling"); 1418 levels = mglevels[0]->levels; 1419 for (i=0; i<levels; i++) mglevels[i]->cycles = n; 1420 PetscFunctionReturn(0); 1421 } 1422 1423 /*@ 1424 PCMGMultiplicativeSetCycles - Sets the number of cycles to use for each preconditioner step 1425 of multigrid when PCMGType of PC_MG_MULTIPLICATIVE is used 1426 1427 Logically Collective on PC 1428 1429 Input Parameters: 1430 + pc - the multigrid context 1431 - n - number of cycles (default is 1) 1432 1433 Options Database Key: 1434 . -pc_mg_multiplicative_cycles n 1435 1436 Level: advanced 1437 1438 Notes: 1439 This is not associated with setting a v or w cycle, that is set with PCMGSetCycleType() 1440 1441 .seealso: PCMGSetCycleTypeOnLevel(), PCMGSetCycleType() 1442 @*/ 1443 PetscErrorCode PCMGMultiplicativeSetCycles(PC pc,PetscInt n) 1444 { 1445 PC_MG *mg = (PC_MG*)pc->data; 1446 1447 PetscFunctionBegin; 1448 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1449 PetscValidLogicalCollectiveInt(pc,n,2); 1450 mg->cyclesperpcapply = n; 1451 PetscFunctionReturn(0); 1452 } 1453 1454 PetscErrorCode PCMGSetGalerkin_MG(PC pc,PCMGGalerkinType use) 1455 { 1456 PC_MG *mg = (PC_MG*)pc->data; 1457 1458 PetscFunctionBegin; 1459 mg->galerkin = use; 1460 PetscFunctionReturn(0); 1461 } 1462 1463 /*@ 1464 PCMGSetGalerkin - Causes the coarser grid matrices to be computed from the 1465 finest grid via the Galerkin process: A_i-1 = r_i * A_i * p_i 1466 1467 Logically Collective on PC 1468 1469 Input Parameters: 1470 + pc - the multigrid context 1471 - use - one of PC_MG_GALERKIN_BOTH,PC_MG_GALERKIN_PMAT,PC_MG_GALERKIN_MAT, or PC_MG_GALERKIN_NONE 1472 1473 Options Database Key: 1474 . -pc_mg_galerkin <both,pmat,mat,none> 1475 1476 Level: intermediate 1477 1478 Notes: 1479 Some codes that use PCMG such as PCGAMG use Galerkin internally while constructing the hierarchy and thus do not 1480 use the PCMG construction of the coarser grids. 1481 1482 .seealso: PCMGGetGalerkin(), PCMGGalerkinType 1483 1484 @*/ 1485 PetscErrorCode PCMGSetGalerkin(PC pc,PCMGGalerkinType use) 1486 { 1487 PetscErrorCode ierr; 1488 1489 PetscFunctionBegin; 1490 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1491 ierr = PetscTryMethod(pc,"PCMGSetGalerkin_C",(PC,PCMGGalerkinType),(pc,use));CHKERRQ(ierr); 1492 PetscFunctionReturn(0); 1493 } 1494 1495 /*@ 1496 PCMGGetGalerkin - Checks if Galerkin multigrid is being used, i.e. 1497 A_i-1 = r_i * A_i * p_i 1498 1499 Not Collective 1500 1501 Input Parameter: 1502 . pc - the multigrid context 1503 1504 Output Parameter: 1505 . galerkin - one of PC_MG_GALERKIN_BOTH,PC_MG_GALERKIN_PMAT,PC_MG_GALERKIN_MAT, PC_MG_GALERKIN_NONE, or PC_MG_GALERKIN_EXTERNAL 1506 1507 Level: intermediate 1508 1509 .seealso: PCMGSetGalerkin(), PCMGGalerkinType 1510 1511 @*/ 1512 PetscErrorCode PCMGGetGalerkin(PC pc,PCMGGalerkinType *galerkin) 1513 { 1514 PC_MG *mg = (PC_MG*)pc->data; 1515 1516 PetscFunctionBegin; 1517 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1518 *galerkin = mg->galerkin; 1519 PetscFunctionReturn(0); 1520 } 1521 1522 PetscErrorCode PCMGSetAdaptInterpolation_MG(PC pc, PetscBool adapt) 1523 { 1524 PC_MG *mg = (PC_MG *) pc->data; 1525 1526 PetscFunctionBegin; 1527 mg->adaptInterpolation = adapt; 1528 PetscFunctionReturn(0); 1529 } 1530 1531 PetscErrorCode PCMGGetAdaptInterpolation_MG(PC pc, PetscBool *adapt) 1532 { 1533 PC_MG *mg = (PC_MG *) pc->data; 1534 1535 PetscFunctionBegin; 1536 *adapt = mg->adaptInterpolation; 1537 PetscFunctionReturn(0); 1538 } 1539 1540 PetscErrorCode PCMGSetAdaptCR_MG(PC pc, PetscBool cr) 1541 { 1542 PC_MG *mg = (PC_MG *) pc->data; 1543 1544 PetscFunctionBegin; 1545 mg->compatibleRelaxation = cr; 1546 PetscFunctionReturn(0); 1547 } 1548 1549 PetscErrorCode PCMGGetAdaptCR_MG(PC pc, PetscBool *cr) 1550 { 1551 PC_MG *mg = (PC_MG *) pc->data; 1552 1553 PetscFunctionBegin; 1554 *cr = mg->compatibleRelaxation; 1555 PetscFunctionReturn(0); 1556 } 1557 1558 /*@ 1559 PCMGSetAdaptInterpolation - Adapt the interpolator based upon a vector space which should be accurately captured by the next coarser mesh, and thus accurately interpolated. 1560 1561 Logically Collective on PC 1562 1563 Input Parameters: 1564 + pc - the multigrid context 1565 - adapt - flag for adaptation of the interpolator 1566 1567 Options Database Keys: 1568 + -pc_mg_adapt_interp - Turn on adaptation 1569 . -pc_mg_adapt_interp_n <int> - The number of modes to use, should be divisible by dimension 1570 - -pc_mg_adapt_interp_coarse_space <type> - The type of coarse space: polynomial, harmonic, eigenvector, generalized_eigenvector 1571 1572 Level: intermediate 1573 1574 .keywords: MG, set, Galerkin 1575 .seealso: PCMGGetAdaptInterpolation(), PCMGSetGalerkin() 1576 @*/ 1577 PetscErrorCode PCMGSetAdaptInterpolation(PC pc, PetscBool adapt) 1578 { 1579 PetscErrorCode ierr; 1580 1581 PetscFunctionBegin; 1582 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1583 ierr = PetscTryMethod(pc,"PCMGSetAdaptInterpolation_C",(PC,PetscBool),(pc,adapt));CHKERRQ(ierr); 1584 PetscFunctionReturn(0); 1585 } 1586 1587 /*@ 1588 PCMGGetAdaptInterpolation - Get the flag to adapt the interpolator based upon a vector space which should be accurately captured by the next coarser mesh, and thus accurately interpolated. 1589 1590 Not collective 1591 1592 Input Parameter: 1593 . pc - the multigrid context 1594 1595 Output Parameter: 1596 . adapt - flag for adaptation of the interpolator 1597 1598 Level: intermediate 1599 1600 .keywords: MG, set, Galerkin 1601 .seealso: PCMGSetAdaptInterpolation(), PCMGSetGalerkin() 1602 @*/ 1603 PetscErrorCode PCMGGetAdaptInterpolation(PC pc, PetscBool *adapt) 1604 { 1605 PetscErrorCode ierr; 1606 1607 PetscFunctionBegin; 1608 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1609 PetscValidPointer(adapt, 2); 1610 ierr = PetscUseMethod(pc,"PCMGGetAdaptInterpolation_C",(PC,PetscBool*),(pc,adapt));CHKERRQ(ierr); 1611 PetscFunctionReturn(0); 1612 } 1613 1614 /*@ 1615 PCMGSetAdaptCR - Monitor the coarse space quality using an auxiliary solve with compatible relaxation. 1616 1617 Logically Collective on PC 1618 1619 Input Parameters: 1620 + pc - the multigrid context 1621 - cr - flag for compatible relaxation 1622 1623 Options Database Keys: 1624 . -pc_mg_adapt_cr - Turn on compatible relaxation 1625 1626 Level: intermediate 1627 1628 .keywords: MG, set, Galerkin 1629 .seealso: PCMGGetAdaptCR(), PCMGSetAdaptInterpolation(), PCMGSetGalerkin() 1630 @*/ 1631 PetscErrorCode PCMGSetAdaptCR(PC pc, PetscBool cr) 1632 { 1633 PetscErrorCode ierr; 1634 1635 PetscFunctionBegin; 1636 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1637 ierr = PetscTryMethod(pc,"PCMGSetAdaptCR_C",(PC,PetscBool),(pc,cr));CHKERRQ(ierr); 1638 PetscFunctionReturn(0); 1639 } 1640 1641 /*@ 1642 PCMGGetAdaptCR - Get the flag to monitor coarse space quality using an auxiliary solve with compatible relaxation. 1643 1644 Not collective 1645 1646 Input Parameter: 1647 . pc - the multigrid context 1648 1649 Output Parameter: 1650 . cr - flag for compatible relaxaion 1651 1652 Level: intermediate 1653 1654 .keywords: MG, set, Galerkin 1655 .seealso: PCMGSetAdaptCR(), PCMGGetAdaptInterpolation(), PCMGSetGalerkin() 1656 @*/ 1657 PetscErrorCode PCMGGetAdaptCR(PC pc, PetscBool *cr) 1658 { 1659 PetscErrorCode ierr; 1660 1661 PetscFunctionBegin; 1662 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1663 PetscValidPointer(cr, 2); 1664 ierr = PetscUseMethod(pc,"PCMGGetAdaptCR_C",(PC,PetscBool*),(pc,cr));CHKERRQ(ierr); 1665 PetscFunctionReturn(0); 1666 } 1667 1668 /*@ 1669 PCMGSetNumberSmooth - Sets the number of pre and post-smoothing steps to use 1670 on all levels. Use PCMGDistinctSmoothUp() to create separate up and down smoothers if you want different numbers of 1671 pre- and post-smoothing steps. 1672 1673 Logically Collective on PC 1674 1675 Input Parameters: 1676 + mg - the multigrid context 1677 - n - the number of smoothing steps 1678 1679 Options Database Key: 1680 . -mg_levels_ksp_max_it <n> - Sets number of pre and post-smoothing steps 1681 1682 Level: advanced 1683 1684 Notes: 1685 this does not set a value on the coarsest grid, since we assume that 1686 there is no separate smooth up on the coarsest grid. 1687 1688 .seealso: PCMGSetDistinctSmoothUp() 1689 @*/ 1690 PetscErrorCode PCMGSetNumberSmooth(PC pc,PetscInt n) 1691 { 1692 PC_MG *mg = (PC_MG*)pc->data; 1693 PC_MG_Levels **mglevels = mg->levels; 1694 PetscErrorCode ierr; 1695 PetscInt i,levels; 1696 1697 PetscFunctionBegin; 1698 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1699 PetscValidLogicalCollectiveInt(pc,n,2); 1700 if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Must set MG levels with PCMGSetLevels() before calling"); 1701 levels = mglevels[0]->levels; 1702 1703 for (i=1; i<levels; i++) { 1704 ierr = KSPSetTolerances(mglevels[i]->smoothu,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);CHKERRQ(ierr); 1705 ierr = KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);CHKERRQ(ierr); 1706 mg->default_smoothu = n; 1707 mg->default_smoothd = n; 1708 } 1709 PetscFunctionReturn(0); 1710 } 1711 1712 /*@ 1713 PCMGSetDistinctSmoothUp - sets the up (post) smoother to be a separate KSP from the down (pre) smoother on all levels 1714 and adds the suffix _up to the options name 1715 1716 Logically Collective on PC 1717 1718 Input Parameters: 1719 . pc - the preconditioner context 1720 1721 Options Database Key: 1722 . -pc_mg_distinct_smoothup 1723 1724 Level: advanced 1725 1726 Notes: 1727 this does not set a value on the coarsest grid, since we assume that 1728 there is no separate smooth up on the coarsest grid. 1729 1730 .seealso: PCMGSetNumberSmooth() 1731 @*/ 1732 PetscErrorCode PCMGSetDistinctSmoothUp(PC pc) 1733 { 1734 PC_MG *mg = (PC_MG*)pc->data; 1735 PC_MG_Levels **mglevels = mg->levels; 1736 PetscErrorCode ierr; 1737 PetscInt i,levels; 1738 KSP subksp; 1739 1740 PetscFunctionBegin; 1741 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1742 if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Must set MG levels with PCMGSetLevels() before calling"); 1743 levels = mglevels[0]->levels; 1744 1745 for (i=1; i<levels; i++) { 1746 const char *prefix = NULL; 1747 /* make sure smoother up and down are different */ 1748 ierr = PCMGGetSmootherUp(pc,i,&subksp);CHKERRQ(ierr); 1749 ierr = KSPGetOptionsPrefix(mglevels[i]->smoothd,&prefix);CHKERRQ(ierr); 1750 ierr = KSPSetOptionsPrefix(subksp,prefix);CHKERRQ(ierr); 1751 ierr = KSPAppendOptionsPrefix(subksp,"up_");CHKERRQ(ierr); 1752 } 1753 PetscFunctionReturn(0); 1754 } 1755 1756 /* No new matrices are created, and the coarse operator matrices are the references to the original ones */ 1757 PetscErrorCode PCGetInterpolations_MG(PC pc,PetscInt *num_levels,Mat *interpolations[]) 1758 { 1759 PC_MG *mg = (PC_MG*)pc->data; 1760 PC_MG_Levels **mglevels = mg->levels; 1761 Mat *mat; 1762 PetscInt l; 1763 PetscErrorCode ierr; 1764 1765 PetscFunctionBegin; 1766 if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling"); 1767 ierr = PetscMalloc1(mg->nlevels,&mat);CHKERRQ(ierr); 1768 for (l=1; l< mg->nlevels; l++) { 1769 mat[l-1] = mglevels[l]->interpolate; 1770 ierr = PetscObjectReference((PetscObject)mat[l-1]);CHKERRQ(ierr); 1771 } 1772 *num_levels = mg->nlevels; 1773 *interpolations = mat; 1774 PetscFunctionReturn(0); 1775 } 1776 1777 /* No new matrices are created, and the coarse operator matrices are the references to the original ones */ 1778 PetscErrorCode PCGetCoarseOperators_MG(PC pc,PetscInt *num_levels,Mat *coarseOperators[]) 1779 { 1780 PC_MG *mg = (PC_MG*)pc->data; 1781 PC_MG_Levels **mglevels = mg->levels; 1782 PetscInt l; 1783 Mat *mat; 1784 PetscErrorCode ierr; 1785 1786 PetscFunctionBegin; 1787 if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling"); 1788 ierr = PetscMalloc1(mg->nlevels,&mat);CHKERRQ(ierr); 1789 for (l=0; l<mg->nlevels-1; l++) { 1790 ierr = KSPGetOperators(mglevels[l]->smoothd,NULL,&(mat[l]));CHKERRQ(ierr); 1791 ierr = PetscObjectReference((PetscObject)mat[l]);CHKERRQ(ierr); 1792 } 1793 *num_levels = mg->nlevels; 1794 *coarseOperators = mat; 1795 PetscFunctionReturn(0); 1796 } 1797 1798 /*@C 1799 PCMGRegisterCoarseSpaceConstructor - Adds a method to the PCMG package for coarse space construction. 1800 1801 Not collective 1802 1803 Input Parameters: 1804 + name - name of the constructor 1805 - function - constructor routine 1806 1807 Notes: 1808 Calling sequence for the routine: 1809 $ my_csp(PC pc, PetscInt l, DM dm, KSP smooth, PetscInt Nc, const Vec initGuess[], Vec **coarseSp) 1810 $ pc - The PC object 1811 $ l - The multigrid level, 0 is the coarse level 1812 $ dm - The DM for this level 1813 $ smooth - The level smoother 1814 $ Nc - The size of the coarse space 1815 $ initGuess - Basis for an initial guess for the space 1816 $ coarseSp - A basis for the computed coarse space 1817 1818 Level: advanced 1819 1820 .seealso: PCMGGetCoarseSpaceConstructor(), PCRegister() 1821 @*/ 1822 PetscErrorCode PCMGRegisterCoarseSpaceConstructor(const char name[], PetscErrorCode (*function)(PC, PetscInt, DM, KSP, PetscInt, const Vec[], Vec **)) 1823 { 1824 PetscErrorCode ierr; 1825 1826 PetscFunctionBegin; 1827 ierr = PCInitializePackage();CHKERRQ(ierr); 1828 ierr = PetscFunctionListAdd(&PCMGCoarseList,name,function);CHKERRQ(ierr); 1829 PetscFunctionReturn(0); 1830 } 1831 1832 /*@C 1833 PCMGGetCoarseSpaceConstructor - Returns the given coarse space construction method. 1834 1835 Not collective 1836 1837 Input Parameter: 1838 . name - name of the constructor 1839 1840 Output Parameter: 1841 . function - constructor routine 1842 1843 Notes: 1844 Calling sequence for the routine: 1845 $ my_csp(PC pc, PetscInt l, DM dm, KSP smooth, PetscInt Nc, const Vec initGuess[], Vec **coarseSp) 1846 $ pc - The PC object 1847 $ l - The multigrid level, 0 is the coarse level 1848 $ dm - The DM for this level 1849 $ smooth - The level smoother 1850 $ Nc - The size of the coarse space 1851 $ initGuess - Basis for an initial guess for the space 1852 $ coarseSp - A basis for the computed coarse space 1853 1854 Level: advanced 1855 1856 .seealso: PCMGRegisterCoarseSpaceConstructor(), PCRegister() 1857 @*/ 1858 PetscErrorCode PCMGGetCoarseSpaceConstructor(const char name[], PetscErrorCode (**function)(PC, PetscInt, DM, KSP, PetscInt, const Vec[], Vec **)) 1859 { 1860 PetscErrorCode ierr; 1861 1862 PetscFunctionBegin; 1863 ierr = PetscFunctionListFind(PCMGCoarseList,name,function);CHKERRQ(ierr); 1864 PetscFunctionReturn(0); 1865 } 1866 1867 /* ----------------------------------------------------------------------------------------*/ 1868 1869 /*MC 1870 PCMG - Use multigrid preconditioning. This preconditioner requires you provide additional 1871 information about the coarser grid matrices and restriction/interpolation operators. 1872 1873 Options Database Keys: 1874 + -pc_mg_levels <nlevels> - number of levels including finest 1875 . -pc_mg_cycle_type <v,w> - provide the cycle desired 1876 . -pc_mg_type <additive,multiplicative,full,kaskade> - multiplicative is the default 1877 . -pc_mg_log - log information about time spent on each level of the solver 1878 . -pc_mg_distinct_smoothup - configure up (after interpolation) and down (before restriction) smoothers separately (with different options prefixes) 1879 . -pc_mg_galerkin <both,pmat,mat,none> - use Galerkin process to compute coarser operators, i.e. Acoarse = R A R' 1880 . -pc_mg_multiplicative_cycles - number of cycles to use as the preconditioner (defaults to 1) 1881 . -pc_mg_dump_matlab - dumps the matrices for each level and the restriction/interpolation matrices 1882 to the Socket viewer for reading from MATLAB. 1883 - -pc_mg_dump_binary - dumps the matrices for each level and the restriction/interpolation matrices 1884 to the binary output file called binaryoutput 1885 1886 Notes: 1887 If one uses a Krylov method such GMRES or CG as the smoother then one must use KSPFGMRES, KSPGCR, or KSPRICHARDSON as the outer Krylov method 1888 1889 When run with a single level the smoother options are used on that level NOT the coarse grid solver options 1890 1891 When run with KSPRICHARDSON the convergence test changes slightly if monitor is turned on. The iteration count may change slightly. This 1892 is because without monitoring the residual norm is computed WITHIN each multigrid cycle on the finest level after the pre-smoothing 1893 (because the residual has just been computed for the multigrid algorithm and is hence available for free) while with monitoring the 1894 residual is computed at the end of each cycle. 1895 1896 Level: intermediate 1897 1898 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType, PCEXOTIC, PCGAMG, PCML, PCHYPRE 1899 PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), PCMGSetCycleType(), 1900 PCMGSetDistinctSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(), 1901 PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(), 1902 PCMGSetCycleTypeOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR() 1903 M*/ 1904 1905 PETSC_EXTERN PetscErrorCode PCCreate_MG(PC pc) 1906 { 1907 PC_MG *mg; 1908 PetscErrorCode ierr; 1909 1910 PetscFunctionBegin; 1911 ierr = PetscNewLog(pc,&mg);CHKERRQ(ierr); 1912 pc->data = (void*)mg; 1913 mg->nlevels = -1; 1914 mg->am = PC_MG_MULTIPLICATIVE; 1915 mg->galerkin = PC_MG_GALERKIN_NONE; 1916 mg->adaptInterpolation = PETSC_FALSE; 1917 mg->Nc = -1; 1918 mg->eigenvalue = -1; 1919 1920 pc->useAmat = PETSC_TRUE; 1921 1922 pc->ops->apply = PCApply_MG; 1923 pc->ops->applytranspose = PCApplyTranspose_MG; 1924 pc->ops->matapply = PCMatApply_MG; 1925 pc->ops->setup = PCSetUp_MG; 1926 pc->ops->reset = PCReset_MG; 1927 pc->ops->destroy = PCDestroy_MG; 1928 pc->ops->setfromoptions = PCSetFromOptions_MG; 1929 pc->ops->view = PCView_MG; 1930 1931 ierr = PetscObjectComposedDataRegister(&mg->eigenvalue);CHKERRQ(ierr); 1932 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGSetGalerkin_C",PCMGSetGalerkin_MG);CHKERRQ(ierr); 1933 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGGetLevels_C",PCMGGetLevels_MG);CHKERRQ(ierr); 1934 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGSetLevels_C",PCMGSetLevels_MG);CHKERRQ(ierr); 1935 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGetInterpolations_C",PCGetInterpolations_MG);CHKERRQ(ierr); 1936 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGetCoarseOperators_C",PCGetCoarseOperators_MG);CHKERRQ(ierr); 1937 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGSetAdaptInterpolation_C",PCMGSetAdaptInterpolation_MG);CHKERRQ(ierr); 1938 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGGetAdaptInterpolation_C",PCMGGetAdaptInterpolation_MG);CHKERRQ(ierr); 1939 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGSetAdaptCR_C",PCMGSetAdaptCR_MG);CHKERRQ(ierr); 1940 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGGetAdaptCR_C",PCMGGetAdaptCR_MG);CHKERRQ(ierr); 1941 PetscFunctionReturn(0); 1942 } 1943