1 2 /* -------------------------------------------------------------------- 3 4 This file implements a Deflation preconditioner in PETSc as part of PC. 5 You can use this as a starting point for implementing your own 6 preconditioner that is not provided with PETSc. (You might also consider 7 just using PCSHELL) 8 9 The following basic routines are required for each preconditioner. 10 PCCreate_XXX() - Creates a preconditioner context 11 PCSetFromOptions_XXX() - Sets runtime options 12 PCApply_XXX() - Applies the preconditioner 13 PCDestroy_XXX() - Destroys the preconditioner context 14 where the suffix "_XXX" denotes a particular implementation, in 15 this case we use _Deflation (e.g., PCCreate_Deflation, PCApply_Deflation). 16 These routines are actually called via the common user interface 17 routines PCCreate(), PCSetFromOptions(), PCApply(), and PCDestroy(), 18 so the application code interface remains identical for all 19 preconditioners. 20 21 Another key routine is: 22 PCSetUp_XXX() - Prepares for the use of a preconditioner 23 by setting data structures and options. The interface routine PCSetUp() 24 is not usually called directly by the user, but instead is called by 25 PCApply() if necessary. 26 27 Additional basic routines are: 28 PCView_XXX() - Prints details of runtime options that 29 have actually been used. 30 These are called by application codes via the interface routines 31 PCView(). 32 33 The various types of solvers (preconditioners, Krylov subspace methods, 34 nonlinear solvers, timesteppers) are all organized similarly, so the 35 above description applies to these categories also. One exception is 36 that the analogues of PCApply() for these components are KSPSolve(), 37 SNESSolve(), and TSSolve(). 38 39 Additional optional functionality unique to preconditioners is left and 40 right symmetric preconditioner application via PCApplySymmetricLeft() 41 and PCApplySymmetricRight(). The Deflation implementation is 42 PCApplySymmetricLeftOrRight_Deflation(). 43 44 -------------------------------------------------------------------- */ 45 46 /* 47 Include files needed for the Deflation preconditioner: 48 pcimpl.h - private include file intended for use by all preconditioners 49 */ 50 51 #include <../src/ksp/pc/impls/deflation/deflation.h> /*I "petscpc.h" I*/ /* includes for fortran wrappers */ 52 53 const char *const PCDeflationTypes[] = {"INIT","PRE","POST","PCDeflationType","PC_DEFLATION_",0}; 54 55 const char *const PCDeflationSpaceTypes[] = { 56 "haar", 57 "jacket-haar", 58 "db2", 59 "db4", 60 "db8", 61 "db16", 62 "biorth22", 63 "meyer", 64 "aggregation", 65 "slepc", 66 "slepc-cheap", 67 "user", 68 "DdefSpaceType", 69 "Ddef_SPACE_", 70 0 71 }; 72 73 74 static PetscErrorCode PCDeflationSetSpace_Deflation(PC pc,Mat W,PetscBool transpose) 75 { 76 PC_Deflation *def = (PC_Deflation*)pc->data; 77 PetscErrorCode ierr; 78 79 PetscFunctionBegin; 80 if (transpose) { 81 def->Wt = W; 82 def->W = NULL; 83 } else { 84 def->W = W; 85 } 86 ierr = PetscObjectReference((PetscObject)W);CHKERRQ(ierr); 87 PetscFunctionReturn(0); 88 } 89 90 /* TODO create PCDeflationSetSpaceTranspose? */ 91 /*@ 92 PCDeflationSetSpace - Set deflation space matrix (or its transpose). 93 94 Logically Collective on PC 95 96 Input Parameters: 97 + pc - the preconditioner context 98 . W - deflation matrix 99 - tranpose - indicates that W is an explicit transpose of the deflation matrix 100 101 Level: intermediate 102 103 .seealso: PCDEFLATION 104 @*/ 105 PetscErrorCode PCDeflationSetSpace(PC pc,Mat W,PetscBool transpose) 106 { 107 PetscErrorCode ierr; 108 109 PetscFunctionBegin; 110 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 111 PetscValidHeaderSpecific(W,MAT_CLASSID,2); 112 PetscValidLogicalCollectiveBool(pc,transpose,3); 113 ierr = PetscTryMethod(pc,"PCDeflationSetSpace_C",(PC,Mat,PetscBool),(pc,W,transpose));CHKERRQ(ierr); 114 PetscFunctionReturn(0); 115 } 116 117 static PetscErrorCode PCDeflationSetLvl_Deflation(PC pc,PetscInt current,PetscInt max) 118 { 119 PC_Deflation *def = (PC_Deflation*)pc->data; 120 121 PetscFunctionBegin; 122 def->nestedlvl = current; 123 def->maxnestedlvl = max; 124 PetscFunctionReturn(0); 125 } 126 127 /*@ 128 PCDeflationSetMaxLvl - Set maximum level of deflation. 129 130 Logically Collective on PC 131 132 Input Parameters: 133 + pc - the preconditioner context 134 - max - maximum deflation level 135 136 Level: intermediate 137 138 .seealso: PCDEFLATION 139 @*/ 140 PetscErrorCode PCDeflationSetMaxLvl(PC pc,PetscInt max) 141 { 142 PetscErrorCode ierr; 143 144 PetscFunctionBegin; 145 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 146 PetscValidLogicalCollectiveInt(pc,max,2); 147 /* TODO allow setting only for the top level */ 148 ierr = PetscTryMethod(pc,"PCDeflationSetLvl_C",(PC,PetscInt,PetscInt),(pc,0,max));CHKERRQ(ierr); 149 PetscFunctionReturn(0); 150 } 151 152 static PetscErrorCode PCDeflationSetPC_Deflation(PC pc,PC apc) 153 { 154 PC_Deflation *def = (PC_Deflation*)pc->data; 155 PetscErrorCode ierr; 156 157 PetscFunctionBegin; 158 ierr = PCDestroy(&def->pc);CHKERRQ(ierr); 159 def->pc = apc; 160 ierr = PetscObjectReference((PetscObject)apc);CHKERRQ(ierr); 161 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)def->pc);CHKERRQ(ierr); 162 PetscFunctionReturn(0); 163 } 164 165 /*@ 166 PCDeflationSetPC - Set additional preconditioner. 167 168 Logically Collective on PC 169 170 Input Parameters: 171 + pc - the preconditioner context 172 - apc - additional preconditioner 173 174 Level: advanced 175 176 .seealso: PCDEFLATION 177 @*/ 178 PetscErrorCode PCDeflationSetPC(PC pc,PC apc) 179 { 180 PetscErrorCode ierr; 181 182 PetscFunctionBegin; 183 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 184 PetscValidHeaderSpecific(apc,PC_CLASSID,2); 185 PetscCheckSameComm(pc,1,apc,2); 186 ierr = PetscTryMethod(pc,"PCDeflationSetPC_C",(PC,PC),(pc,apc));CHKERRQ(ierr); 187 PetscFunctionReturn(0); 188 } 189 190 /* 191 x <- x + W*(W'*A*W)^{-1}*W'*r = x + Q*r 192 */ 193 static PetscErrorCode PCPreSolve_Deflation(PC pc,KSP ksp,Vec b, Vec x) 194 { 195 PC_Deflation *def = (PC_Deflation*)pc->data; 196 Mat A; 197 Vec r,w1,w2; 198 PetscBool nonzero; 199 PetscErrorCode ierr; 200 201 PetscFunctionBegin; 202 w1 = def->workcoarse[0]; 203 w2 = def->workcoarse[1]; 204 r = def->work; 205 ierr = PCGetOperators(pc,NULL,&A);CHKERRQ(ierr); 206 207 ierr = KSPGetInitialGuessNonzero(ksp,&nonzero);CHKERRQ(ierr); 208 ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr); 209 if (nonzero) { 210 ierr = MatMult(A,x,r);CHKERRQ(ierr); /* r <- b - Ax */ 211 ierr = VecAYPX(r,-1.0,b);CHKERRQ(ierr); 212 } else { 213 ierr = VecCopy(b,r);CHKERRQ(ierr); /* r <- b (x is 0) */ 214 } 215 216 ierr = MatMultTranspose(def->W,r,w1);CHKERRQ(ierr); /* w1 <- W'*r */ 217 ierr = KSPSolve(def->WtAWinv,w1,w2);CHKERRQ(ierr); /* w2 <- (W'*A*W)^{-1}*w1 */ 218 ierr = MatMult(def->W,w2,r);CHKERRQ(ierr); /* r <- W*w2 */ 219 ierr = VecAYPX(x,1.0,r);CHKERRQ(ierr); 220 PetscFunctionReturn(0); 221 } 222 223 /* 224 if (def->correct) { 225 z <- r - W*(W'*A*W)^{-1}*W'*(A*r -r) = (P-Q)*M^{-1}*r 226 } else { 227 z <- r - W*(W'*A*W)^{-1}*W'*A*r = P*M^{-1}*r 228 } 229 */ 230 static PetscErrorCode PCApply_Deflation(PC pc,Vec r,Vec z) 231 { 232 PC_Deflation *def = (PC_Deflation*)pc->data; 233 Mat A; 234 Vec u,w1,w2; 235 PetscErrorCode ierr; 236 237 PetscFunctionBegin; 238 w1 = def->workcoarse[0]; 239 w2 = def->workcoarse[1]; 240 u = def->work; 241 ierr = PCGetOperators(pc,NULL,&A);CHKERRQ(ierr); 242 243 ierr = PCApply(def->pc,r,z);CHKERRQ(ierr); 244 if (!def->init) { 245 if (!def->AW) { 246 ierr = MatMult(A,z,u);CHKERRQ(ierr); /* u <- A*z */ 247 /* TODO correct const */ 248 if (def->correct) ierr = VecAXPY(u,-1.0,z);CHKERRQ(ierr); /* u <- A*z -z */ 249 ierr = MatMultTranspose(def->W,u,w1);CHKERRQ(ierr); /* w1 <- W'*u */ 250 } else { 251 /* ONLY if A SYM */ 252 ierr = MatMultTranspose(def->AW,z,w1);CHKERRQ(ierr); /* w1 <- W'*A*r */ 253 if (def->correct) { 254 ierr = MatMultTranspose(def->W,z,w2);CHKERRQ(ierr); /* w2 <- W'*u */ 255 ierr = VecAXPY(w1,-1.0,w2);CHKERRQ(ierr); /* w1 <- w1 - w2 */ 256 } 257 } 258 ierr = KSPSolve(def->WtAWinv,w1,w2);CHKERRQ(ierr); /* w2 <- (W'*A*W)^{-1}*w1 */ 259 ierr = MatMult(def->W,w2,u);CHKERRQ(ierr); /* u <- W*w2 */ 260 ierr = VecAXPY(z,-1.0,u);CHKERRQ(ierr); /* z <- z - u */ 261 } 262 PetscFunctionReturn(0); 263 } 264 265 static PetscErrorCode PCDeflationSetType_Deflation(PC pc,PCDeflationType type) 266 { 267 PC_Deflation *def = (PC_Deflation*)pc->data; 268 269 PetscFunctionBegin; 270 def->init = PETSC_FALSE; 271 def->pre = PETSC_FALSE; 272 if (type == PC_DEFLATION_POST) { 273 //pc->ops->postsolve = PCPostSolve_Deflation; 274 pc->ops->presolve = 0; 275 } else { 276 pc->ops->presolve = PCPreSolve_Deflation; 277 pc->ops->postsolve = 0; 278 if (type == PC_DEFLATION_INIT) { 279 def->init = PETSC_TRUE; 280 pc->ops->apply = 0; 281 } else { 282 def->pre = PETSC_TRUE; 283 } 284 } 285 PetscFunctionReturn(0); 286 } 287 288 /*@ 289 PCDeflationSetType - Causes the deflation preconditioner to use only a special 290 initial gues or pre/post solve solution update 291 292 Logically Collective on PC 293 294 Input Parameters: 295 + pc - the preconditioner context 296 - type - PC_DEFLATION_PRE, PC_DEFLATION_INIT, PC_DEFLATION_POST 297 298 Options Database Key: 299 . -pc_deflation_type <pre,init,post> 300 301 Level: intermediate 302 303 Concepts: Deflation preconditioner 304 305 .seealso: PCDeflationGetType() 306 @*/ 307 PetscErrorCode PCDeflationSetType(PC pc,PCDeflationType type) 308 { 309 PetscErrorCode ierr; 310 311 PetscFunctionBegin; 312 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 313 ierr = PetscTryMethod(pc,"PCDeflationSetType_C",(PC,PCDeflationType),(pc,type));CHKERRQ(ierr); 314 PetscFunctionReturn(0); 315 } 316 317 static PetscErrorCode PCDeflationGetType_Deflation(PC pc,PCDeflationType *type) 318 { 319 PC_Deflation *def = (PC_Deflation*)pc->data; 320 321 PetscFunctionBegin; 322 if (def->init) { 323 *type = PC_DEFLATION_INIT; 324 } else if (def->pre) { 325 *type = PC_DEFLATION_PRE; 326 } else { 327 *type = PC_DEFLATION_POST; 328 } 329 PetscFunctionReturn(0); 330 } 331 332 /*@ 333 PCDeflationGetType - Gets how the diagonal matrix is produced for the preconditioner 334 335 Not Collective on PC 336 337 Input Parameter: 338 . pc - the preconditioner context 339 340 Output Parameter: 341 - type - PC_DEFLATION_PRE, PC_DEFLATION_INIT, PC_DEFLATION_POST 342 343 Level: intermediate 344 345 Concepts: Deflation preconditioner 346 347 .seealso: PCDeflationSetType() 348 @*/ 349 PetscErrorCode PCDeflationGetType(PC pc,PCDeflationType *type) 350 { 351 PetscErrorCode ierr; 352 353 PetscFunctionBegin; 354 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 355 ierr = PetscUseMethod(pc,"PCDeflationGetType_C",(PC,PCDeflationType*),(pc,type));CHKERRQ(ierr); 356 PetscFunctionReturn(0); 357 } 358 359 static PetscErrorCode PCSetUp_Deflation(PC pc) 360 { 361 PC_Deflation *def = (PC_Deflation*)pc->data; 362 KSP innerksp; 363 PC pcinner; 364 Mat Amat,nextDef=NULL,*mats; 365 PetscInt i,m,red,size,commsize; 366 PetscBool match,flgspd,transp=PETSC_FALSE; 367 MatCompositeType ctype; 368 MPI_Comm comm; 369 const char *prefix; 370 PetscErrorCode ierr; 371 372 PetscFunctionBegin; 373 if (pc->setupcalled) PetscFunctionReturn(0); 374 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 375 ierr = PCGetOperators(pc,NULL,&Amat);CHKERRQ(ierr); 376 377 /* compute a deflation space */ 378 if (def->W || def->Wt) { 379 def->spacetype = PC_DEFLATION_SPACE_USER; 380 } else { 381 ierr = PCDeflationComputeSpace(pc);CHKERRQ(ierr); 382 } 383 384 /* nested deflation */ 385 if (def->W) { 386 ierr = PetscObjectTypeCompare((PetscObject)def->W,MATCOMPOSITE,&match);CHKERRQ(ierr); 387 if (match) { 388 ierr = MatCompositeGetType(def->W,&ctype);CHKERRQ(ierr); 389 ierr = MatCompositeGetNumberMat(def->W,&size);CHKERRQ(ierr); 390 } 391 } else { 392 ierr = MatCreateTranspose(def->Wt,&def->W);CHKERRQ(ierr); 393 ierr = PetscObjectTypeCompare((PetscObject)def->Wt,MATCOMPOSITE,&match);CHKERRQ(ierr); 394 if (match) { 395 ierr = MatCompositeGetType(def->Wt,&ctype);CHKERRQ(ierr); 396 ierr = MatCompositeGetNumberMat(def->Wt,&size);CHKERRQ(ierr); 397 } 398 transp = PETSC_TRUE; 399 } 400 if (match && ctype == MAT_COMPOSITE_MULTIPLICATIVE) { 401 if (!transp) { 402 if (def->nestedlvl < def->maxnestedlvl) { 403 ierr = PetscMalloc1(size,&mats);CHKERRQ(ierr); 404 for (i=0; i<size; i++) { 405 ierr = MatCompositeGetMat(def->W,i,&mats[i]);CHKERRQ(ierr); 406 } 407 size -= 1; 408 ierr = MatDestroy(&def->W);CHKERRQ(ierr); 409 def->W = mats[size]; 410 ierr = PetscObjectReference((PetscObject)mats[size]);CHKERRQ(ierr); 411 if (size > 1) { 412 ierr = MatCreateComposite(comm,size,mats,&nextDef);CHKERRQ(ierr); 413 ierr = MatCompositeSetType(nextDef,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr); 414 } else { 415 nextDef = mats[0]; 416 ierr = PetscObjectReference((PetscObject)mats[0]);CHKERRQ(ierr); 417 } 418 ierr = PetscFree(mats);CHKERRQ(ierr); 419 } else { 420 /* ierr = MatCompositeSetMergeType(def->W,MAT_COMPOSITE_MERGE_LEFT);CHKERRQ(ierr); */ 421 ierr = MatCompositeMerge(def->W);CHKERRQ(ierr); 422 } 423 } else { 424 if (def->nestedlvl < def->maxnestedlvl) { 425 ierr = PetscMalloc1(size,&mats);CHKERRQ(ierr); 426 for (i=0; i<size; i++) { 427 ierr = MatCompositeGetMat(def->Wt,i,&mats[i]);CHKERRQ(ierr); 428 } 429 size -= 1; 430 ierr = MatDestroy(&def->Wt);CHKERRQ(ierr); 431 def->Wt = mats[0]; 432 ierr = PetscObjectReference((PetscObject)mats[0]);CHKERRQ(ierr); 433 if (size > 1) { 434 ierr = MatCreateComposite(comm,size,&mats[1],&nextDef);CHKERRQ(ierr); 435 ierr = MatCompositeSetType(nextDef,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr); 436 } else { 437 nextDef = mats[1]; 438 ierr = PetscObjectReference((PetscObject)mats[1]);CHKERRQ(ierr); 439 } 440 ierr = PetscFree(mats);CHKERRQ(ierr); 441 } else { 442 /* ierr = MatCompositeSetMergeType(def->W,MAT_COMPOSITE_MERGE_LEFT);CHKERRQ(ierr); */ 443 ierr = MatCompositeMerge(def->Wt);CHKERRQ(ierr); 444 } 445 } 446 } 447 448 if (transp) { 449 ierr = MatDestroy(&def->W);CHKERRQ(ierr); 450 ierr = MatTranspose(def->Wt,MAT_INITIAL_MATRIX,&def->W);CHKERRQ(ierr); 451 } 452 453 454 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr); 455 456 /* setup coarse problem */ 457 if (!def->WtAWinv) { 458 ierr = MatGetSize(def->W,NULL,&m);CHKERRQ(ierr); 459 if (!def->WtAW) { 460 /* TODO add implicit product version ? */ 461 if (!def->AW) { 462 ierr = MatPtAP(Amat,def->W,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtAW);CHKERRQ(ierr); 463 } else { 464 ierr = MatTransposeMatMult(def->W,def->AW,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtAW);CHKERRQ(ierr); 465 } 466 /* TODO create MatInheritOption(Mat,MatOption) */ 467 ierr = MatGetOption(Amat,MAT_SPD,&flgspd);CHKERRQ(ierr); 468 ierr = MatSetOption(def->WtAW,MAT_SPD,flgspd);CHKERRQ(ierr); 469 #if defined(PETSC_USE_DEBUG) 470 /* Check WtAW is not sigular */ 471 PetscReal *norms; 472 ierr = PetscMalloc1(m,&norms);CHKERRQ(ierr); 473 ierr = MatGetColumnNorms(def->WtAW,NORM_INFINITY,norms);CHKERRQ(ierr); 474 for (i=0; i<m; i++) { 475 if (norms[i] < 100*PETSC_MACHINE_EPSILON) { 476 SETERRQ1(comm,PETSC_ERR_SUP,"Column %D of W is in kernel of A.",i); 477 } 478 } 479 ierr = PetscFree(norms);CHKERRQ(ierr); 480 #endif 481 } else { 482 ierr = MatGetOption(def->WtAW,MAT_SPD,&flgspd);CHKERRQ(ierr); 483 } 484 /* TODO use MATINV */ 485 ierr = KSPCreate(comm,&def->WtAWinv);CHKERRQ(ierr); 486 ierr = KSPSetOperators(def->WtAWinv,def->WtAW,def->WtAW);CHKERRQ(ierr); 487 ierr = KSPGetPC(def->WtAWinv,&pcinner);CHKERRQ(ierr); 488 /* Setup KSP and PC */ 489 if (nextDef) { /* next level for multilevel deflation */ 490 innerksp = def->WtAWinv; 491 /* set default KSPtype */ 492 if (!def->ksptype) { 493 def->ksptype = KSPFGMRES; 494 if (flgspd) { /* SPD system */ 495 def->ksptype = KSPFCG; 496 } 497 } 498 ierr = KSPSetType(innerksp,def->ksptype);CHKERRQ(ierr); /* TODO iherit from KSP */ 499 ierr = PCSetType(pcinner,PCDEFLATION);CHKERRQ(ierr); /* TODO create coarse preconditinoner M_c = WtMW ? */ 500 ierr = PCDeflationSetSpace(pcinner,nextDef,transp);CHKERRQ(ierr); 501 ierr = PCDeflationSetLvl_Deflation(pcinner,def->nestedlvl+1,def->maxnestedlvl);CHKERRQ(ierr); 502 /* inherit options TODO if not set */ 503 ((PC_Deflation*)(pcinner->data))->ksptype = def->ksptype; 504 ((PC_Deflation*)(pcinner->data))->correct = def->correct; 505 ((PC_Deflation*)(pcinner->data))->adaptiveconv = def->adaptiveconv; 506 ((PC_Deflation*)(pcinner->data))->adaptiveconst = def->adaptiveconst; 507 ierr = MatDestroy(&nextDef);CHKERRQ(ierr); 508 } else { /* the last level */ 509 ierr = KSPSetType(def->WtAWinv,KSPPREONLY);CHKERRQ(ierr); 510 ierr = PCSetType(pcinner,PCTELESCOPE);CHKERRQ(ierr); 511 /* ugly hack to not have overwritten PCTELESCOPE */ 512 if (prefix) { 513 ierr = KSPSetOptionsPrefix(def->WtAWinv,prefix);CHKERRQ(ierr); 514 } 515 ierr = KSPAppendOptionsPrefix(def->WtAWinv,"tel_");CHKERRQ(ierr); 516 ierr = PCSetFromOptions(pcinner);CHKERRQ(ierr); 517 /* Reduction factor choice */ 518 red = def->reductionfact; 519 if (red < 0) { 520 ierr = MPI_Comm_size(comm,&commsize);CHKERRQ(ierr); 521 red = ceil((float)commsize/ceil((float)m/commsize)); 522 ierr = PetscObjectTypeCompareAny((PetscObject)(def->WtAW),&match,MATSEQDENSE,MATMPIDENSE,MATDENSE,"");CHKERRQ(ierr); 523 if (match) red = commsize; 524 ierr = PetscInfo1(pc,"Auto choosing reduction factor %D\n",red);CHKERRQ(ierr); /* TODO add level? */ 525 } 526 ierr = PCTelescopeSetReductionFactor(pcinner,red);CHKERRQ(ierr); 527 ierr = PCSetUp(pcinner);CHKERRQ(ierr); 528 ierr = PCTelescopeGetKSP(pcinner,&innerksp);CHKERRQ(ierr); 529 if (innerksp) { 530 ierr = KSPGetPC(innerksp,&pcinner);CHKERRQ(ierr); 531 /* TODO Cholesky if flgspd? */ 532 ierr = PCSetType(pcinner,PCLU);CHKERRQ(ierr); 533 //TODO remove explicit matSolverPackage 534 if (commsize == red) { 535 ierr = PCFactorSetMatSolverType(pcinner,MATSOLVERSUPERLU);CHKERRQ(ierr); 536 } else { 537 ierr = PCFactorSetMatSolverType(pcinner,MATSOLVERSUPERLU_DIST);CHKERRQ(ierr); 538 } 539 } 540 } 541 542 if (innerksp) { 543 /* TODO use def_[lvl]_ if lvl > 0? */ 544 if (prefix) { 545 ierr = KSPSetOptionsPrefix(innerksp,prefix);CHKERRQ(ierr); 546 } 547 ierr = KSPAppendOptionsPrefix(innerksp,"def_");CHKERRQ(ierr); 548 ierr = KSPSetFromOptions(innerksp);CHKERRQ(ierr); 549 ierr = KSPSetUp(innerksp);CHKERRQ(ierr); 550 } 551 /* TODO: check if WtAWinv is KSP and move following from this if */ 552 //if (def->adaptiveconv) { 553 // PetscReal *rnorm; 554 // PetscNew(&rnorm); 555 // ierr = KSPSetConvergenceTest(def->WtAWinv,KSPDCGConvergedAdaptive_DCG,rnorm,NULL);CHKERRQ(ierr); 556 //} 557 } 558 ierr = KSPSetFromOptions(def->WtAWinv);CHKERRQ(ierr); 559 ierr = KSPSetUp(def->WtAWinv);CHKERRQ(ierr); 560 561 if (!def->pc) { 562 ierr = PCCreate(comm,&def->pc);CHKERRQ(ierr); 563 ierr = PCSetOperators(def->pc,Amat,Amat);CHKERRQ(ierr); 564 ierr = PCSetType(def->pc,PCNONE);CHKERRQ(ierr); 565 if (prefix) { 566 ierr = PCSetOptionsPrefix(def->pc,prefix);CHKERRQ(ierr); 567 } 568 ierr = PCAppendOptionsPrefix(def->pc,"def_pc_");CHKERRQ(ierr); 569 ierr = PCSetFromOptions(def->pc);CHKERRQ(ierr); 570 ierr = PCSetUp(def->pc);CHKERRQ(ierr); 571 } 572 573 /* create work vecs */ 574 ierr = MatCreateVecs(Amat,NULL,&def->work);CHKERRQ(ierr); 575 ierr = KSPCreateVecs(def->WtAWinv,2,&def->workcoarse,0,NULL);CHKERRQ(ierr); 576 PetscFunctionReturn(0); 577 } 578 579 static PetscErrorCode PCReset_Deflation(PC pc) 580 { 581 PC_Deflation *def = (PC_Deflation*)pc->data; 582 PetscErrorCode ierr; 583 584 PetscFunctionBegin; 585 ierr = VecDestroy(&def->work);CHKERRQ(ierr); 586 ierr = VecDestroyVecs(2,&def->workcoarse);CHKERRQ(ierr); 587 ierr = MatDestroy(&def->W);CHKERRQ(ierr); 588 ierr = MatDestroy(&def->Wt);CHKERRQ(ierr); 589 ierr = MatDestroy(&def->AW);CHKERRQ(ierr); 590 ierr = MatDestroy(&def->WtAW);CHKERRQ(ierr); 591 ierr = KSPDestroy(&def->WtAWinv);CHKERRQ(ierr); 592 ierr = PCDestroy(&def->pc);CHKERRQ(ierr); 593 PetscFunctionReturn(0); 594 } 595 596 /* 597 PCDestroy_Deflation - Destroys the private context for the Deflation preconditioner 598 that was created with PCCreate_Deflation(). 599 600 Input Parameter: 601 . pc - the preconditioner context 602 603 Application Interface Routine: PCDestroy() 604 */ 605 static PetscErrorCode PCDestroy_Deflation(PC pc) 606 { 607 PetscErrorCode ierr; 608 609 PetscFunctionBegin; 610 ierr = PCReset_Deflation(pc);CHKERRQ(ierr); 611 ierr = PetscFree(pc->data);CHKERRQ(ierr); 612 PetscFunctionReturn(0); 613 } 614 615 static PetscErrorCode PCView_Deflation(PC pc,PetscViewer viewer) 616 { 617 PC_Deflation *def = (PC_Deflation*)pc->data; 618 PetscBool iascii; 619 PetscErrorCode ierr; 620 621 PetscFunctionBegin; 622 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 623 if (iascii) { 624 //if (cg->adaptiveconv) { 625 // ierr = PetscViewerASCIIPrintf(viewer," DCG: using adaptive precision for inner solve with C=%.1e\n",cg->adaptiveconst);CHKERRQ(ierr); 626 //} 627 if (def->correct) { 628 ierr = PetscViewerASCIIPrintf(viewer," Using CP correction\n");CHKERRQ(ierr); 629 } 630 if (!def->nestedlvl) { 631 ierr = PetscViewerASCIIPrintf(viewer," Deflation space type: %s\n",PCDeflationSpaceTypes[def->spacetype]);CHKERRQ(ierr); 632 ierr = PetscViewerASCIIPrintf(viewer," DCG %s\n",def->extendsp ? "extended" : "truncated");CHKERRQ(ierr); 633 } 634 635 ierr = PetscViewerASCIIPrintf(viewer,"--- Additional PC\n");CHKERRQ(ierr); 636 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 637 ierr = PCView(def->pc,viewer);CHKERRQ(ierr); 638 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 639 640 ierr = PetscViewerASCIIPrintf(viewer,"--- Coarse solver\n");CHKERRQ(ierr); 641 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 642 ierr = KSPView(def->WtAWinv,viewer);CHKERRQ(ierr); 643 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 644 } 645 PetscFunctionReturn(0); 646 } 647 648 static PetscErrorCode PCSetFromOptions_Deflation(PetscOptionItems *PetscOptionsObject,PC pc) 649 { 650 PC_Deflation *def = (PC_Deflation*)pc->data; 651 PetscErrorCode ierr; 652 PetscBool flg; 653 PCDeflationType deflt,type; 654 655 PetscFunctionBegin; 656 ierr = PCDeflationGetType(pc,&deflt);CHKERRQ(ierr); 657 ierr = PetscOptionsHead(PetscOptionsObject,"Deflation options");CHKERRQ(ierr); 658 ierr = PetscOptionsEnum("-pc_deflation_type","Determine type of deflation","PCDeflationSetType",PCDeflationTypes,(PetscEnum)deflt,(PetscEnum*)&type,&flg);CHKERRQ(ierr); 659 if (flg) { 660 ierr = PCDeflationSetType(pc,type);CHKERRQ(ierr); 661 } 662 ierr = PetscOptionsEnum("-pc_deflation_compute_space","Compute deflation space","PCDeflationSetSpace",PCDeflationSpaceTypes,(PetscEnum)def->spacetype,(PetscEnum*)&def->spacetype,NULL);CHKERRQ(ierr); 663 ierr = PetscOptionsInt("-pc_deflation_compute_space_size","Set size of the deflation space to compute","PCDeflationSetSpace",def->spacesize,&def->spacesize,NULL);CHKERRQ(ierr); 664 ierr = PetscOptionsBool("-pc_deflation_space_extend","Extend deflation space instead of truncating (wavelets)","PCDeflation",def->extendsp,&def->extendsp,NULL);CHKERRQ(ierr); 665 //TODO add set function and fix manpages 666 ierr = PetscOptionsBool("-pc_deflation_initdef","Use only initialization step - Initdef","PCDeflation",def->init,&def->init,NULL);CHKERRQ(ierr); 667 ierr = PetscOptionsBool("-pc_deflation_correct","Add Qr to descent direction","PCDeflation",def->correct,&def->correct,NULL);CHKERRQ(ierr); 668 ierr = PetscOptionsBool("-pc_deflation_adaptive","Adaptive stopping criteria","PCDeflation",def->adaptiveconv,&def->adaptiveconv,NULL);CHKERRQ(ierr); 669 ierr = PetscOptionsReal("-pc_deflation_adaptive_const","Adaptive stopping criteria constant","PCDeflation",def->adaptiveconst,&def->adaptiveconst,NULL);CHKERRQ(ierr); 670 ierr = PetscOptionsInt("-pc_deflation_redfact","Reduction factor for coarse problem solution","PCDeflation",def->reductionfact,&def->reductionfact,NULL);CHKERRQ(ierr); 671 ierr = PetscOptionsInt("-pc_deflation_max_nested_lvl","Maximum of nested deflation levels","PCDeflation",def->maxnestedlvl,&def->maxnestedlvl,NULL);CHKERRQ(ierr); 672 ierr = PetscOptionsTail();CHKERRQ(ierr); 673 PetscFunctionReturn(0); 674 } 675 676 /*MC 677 PCDEFLATION - Deflation preconditioner shifts part of the spectrum to zero (deflates) 678 or to a predefined value 679 680 Options Database Key: 681 + -pc_deflation_type <init,pre,post> - selects approach to deflation (default: pre) 682 - -pc_jacobi_abs - use the absolute value of the diagonal entry 683 684 Level: beginner 685 686 Notes: 687 todo 688 689 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 690 PCDeflationSetType(), PCDeflationSetSpace() 691 M*/ 692 693 PETSC_EXTERN PetscErrorCode PCCreate_Deflation(PC pc) 694 { 695 PC_Deflation *def; 696 PetscErrorCode ierr; 697 698 PetscFunctionBegin; 699 ierr = PetscNewLog(pc,&def);CHKERRQ(ierr); 700 pc->data = (void*)def; 701 702 def->init = PETSC_FALSE; 703 def->pre = PETSC_TRUE; 704 def->correct = PETSC_FALSE; 705 def->truenorm = PETSC_TRUE; 706 def->reductionfact = -1; 707 def->spacetype = PC_DEFLATION_SPACE_HAAR; 708 def->spacesize = 1; 709 def->extendsp = PETSC_FALSE; 710 def->nestedlvl = 0; 711 def->maxnestedlvl = 0; 712 def->adaptiveconv = PETSC_FALSE; 713 def->adaptiveconst = 1.0; 714 715 /* 716 Set the pointers for the functions that are provided above. 717 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 718 are called, they will automatically call these functions. Note we 719 choose not to provide a couple of these functions since they are 720 not needed. 721 */ 722 pc->ops->apply = PCApply_Deflation; 723 pc->ops->applytranspose = PCApply_Deflation; 724 pc->ops->presolve = PCPreSolve_Deflation; 725 pc->ops->postsolve = 0; 726 pc->ops->setup = PCSetUp_Deflation; 727 pc->ops->reset = PCReset_Deflation; 728 pc->ops->destroy = PCDestroy_Deflation; 729 pc->ops->setfromoptions = PCSetFromOptions_Deflation; 730 pc->ops->view = PCView_Deflation; 731 pc->ops->applyrichardson = 0; 732 pc->ops->applysymmetricleft = 0; 733 pc->ops->applysymmetricright = 0; 734 735 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetType_C",PCDeflationSetType_Deflation);CHKERRQ(ierr); 736 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationGetType_C",PCDeflationGetType_Deflation);CHKERRQ(ierr); 737 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetSpace_C",PCDeflationSetSpace_Deflation);CHKERRQ(ierr); 738 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetLvl_C",PCDeflationSetLvl_Deflation);CHKERRQ(ierr); 739 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetPC_C",PCDeflationSetPC_Deflation);CHKERRQ(ierr); 740 PetscFunctionReturn(0); 741 } 742 743