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 if (current) 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 ierr = PetscTryMethod(pc,"PCDeflationSetLvl_C",(PC,PetscInt,PetscInt),(pc,0,max));CHKERRQ(ierr); 148 PetscFunctionReturn(0); 149 } 150 151 static PetscErrorCode PCDeflationSetPC_Deflation(PC pc,PC apc) 152 { 153 PC_Deflation *def = (PC_Deflation*)pc->data; 154 PetscErrorCode ierr; 155 156 PetscFunctionBegin; 157 ierr = PCDestroy(&def->pc);CHKERRQ(ierr); 158 def->pc = apc; 159 ierr = PetscObjectReference((PetscObject)apc);CHKERRQ(ierr); 160 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)def->pc);CHKERRQ(ierr); 161 PetscFunctionReturn(0); 162 } 163 164 /*@ 165 PCDeflationSetPC - Set additional preconditioner. 166 167 Logically Collective on PC 168 169 Input Parameters: 170 + pc - the preconditioner context 171 - apc - additional preconditioner 172 173 Level: advanced 174 175 .seealso: PCDEFLATION 176 @*/ 177 PetscErrorCode PCDeflationSetPC(PC pc,PC apc) 178 { 179 PetscErrorCode ierr; 180 181 PetscFunctionBegin; 182 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 183 PetscValidHeaderSpecific(apc,PC_CLASSID,2); 184 PetscCheckSameComm(pc,1,apc,2); 185 ierr = PetscTryMethod(pc,"PCDeflationSetPC_C",(PC,PC),(pc,apc));CHKERRQ(ierr); 186 PetscFunctionReturn(0); 187 } 188 189 /* 190 x <- x + W*(W'*A*W)^{-1}*W'*r = x + Q*r 191 */ 192 static PetscErrorCode PCPreSolve_Deflation(PC pc,KSP ksp,Vec b, Vec x) 193 { 194 PC_Deflation *def = (PC_Deflation*)pc->data; 195 Mat A; 196 Vec r,w1,w2; 197 PetscBool nonzero; 198 PetscErrorCode ierr; 199 200 PetscFunctionBegin; 201 w1 = def->workcoarse[0]; 202 w2 = def->workcoarse[1]; 203 r = def->work; 204 ierr = PCGetOperators(pc,NULL,&A);CHKERRQ(ierr); 205 206 ierr = KSPGetInitialGuessNonzero(ksp,&nonzero);CHKERRQ(ierr); 207 ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr); 208 if (nonzero) { 209 ierr = MatMult(A,x,r);CHKERRQ(ierr); /* r <- b - Ax */ 210 ierr = VecAYPX(r,-1.0,b);CHKERRQ(ierr); 211 } else { 212 ierr = VecCopy(b,r);CHKERRQ(ierr); /* r <- b (x is 0) */ 213 } 214 215 ierr = MatMultTranspose(def->W,r,w1);CHKERRQ(ierr); /* w1 <- W'*r */ 216 ierr = KSPSolve(def->WtAWinv,w1,w2);CHKERRQ(ierr); /* w2 <- (W'*A*W)^{-1}*w1 */ 217 ierr = MatMult(def->W,w2,r);CHKERRQ(ierr); /* r <- W*w2 */ 218 ierr = VecAYPX(x,1.0,r);CHKERRQ(ierr); 219 PetscFunctionReturn(0); 220 } 221 222 /* 223 if (def->correct) { 224 z <- r - W*(W'*A*W)^{-1}*W'*(A*r -r) = (P-Q)*M^{-1}*r 225 } else { 226 z <- r - W*(W'*A*W)^{-1}*W'*A*r = P*M^{-1}*r 227 } 228 */ 229 static PetscErrorCode PCApply_Deflation(PC pc,Vec r,Vec z) 230 { 231 PC_Deflation *def = (PC_Deflation*)pc->data; 232 Mat A; 233 Vec u,w1,w2; 234 PetscErrorCode ierr; 235 236 PetscFunctionBegin; 237 w1 = def->workcoarse[0]; 238 w2 = def->workcoarse[1]; 239 u = def->work; 240 ierr = PCGetOperators(pc,NULL,&A);CHKERRQ(ierr); 241 242 ierr = PCApply(def->pc,r,z);CHKERRQ(ierr); 243 if (!def->init) { 244 if (!def->AW) { 245 ierr = MatMult(A,z,u);CHKERRQ(ierr); /* u <- A*z */ 246 /* TODO correct const */ 247 if (def->correct) ierr = VecAXPY(u,-1.0,z);CHKERRQ(ierr); /* u <- A*z -z */ 248 ierr = MatMultTranspose(def->W,u,w1);CHKERRQ(ierr); /* w1 <- W'*u */ 249 } else { 250 /* ONLY if A SYM */ 251 ierr = MatMultTranspose(def->AW,z,w1);CHKERRQ(ierr); /* w1 <- W'*A*r */ 252 if (def->correct) { 253 ierr = MatMultTranspose(def->W,z,w2);CHKERRQ(ierr); /* w2 <- W'*u */ 254 ierr = VecAXPY(w1,-1.0,w2);CHKERRQ(ierr); /* w1 <- w1 - w2 */ 255 } 256 } 257 ierr = KSPSolve(def->WtAWinv,w1,w2);CHKERRQ(ierr); /* w2 <- (W'*A*W)^{-1}*w1 */ 258 ierr = MatMult(def->W,w2,u);CHKERRQ(ierr); /* u <- W*w2 */ 259 ierr = VecAXPY(z,-1.0,u);CHKERRQ(ierr); /* z <- z - u */ 260 } 261 PetscFunctionReturn(0); 262 } 263 264 static PetscErrorCode PCDeflationSetType_Deflation(PC pc,PCDeflationType type) 265 { 266 PC_Deflation *def = (PC_Deflation*)pc->data; 267 268 PetscFunctionBegin; 269 def->init = PETSC_FALSE; 270 def->pre = PETSC_FALSE; 271 if (type == PC_DEFLATION_POST) { 272 //pc->ops->postsolve = PCPostSolve_Deflation; 273 pc->ops->presolve = 0; 274 } else { 275 pc->ops->presolve = PCPreSolve_Deflation; 276 pc->ops->postsolve = 0; 277 if (type == PC_DEFLATION_INIT) { 278 def->init = PETSC_TRUE; 279 pc->ops->apply = 0; 280 } else { 281 def->pre = PETSC_TRUE; 282 } 283 } 284 PetscFunctionReturn(0); 285 } 286 287 /*@ 288 PCDeflationSetType - Causes the deflation preconditioner to use only a special 289 initial gues or pre/post solve solution update 290 291 Logically Collective on PC 292 293 Input Parameters: 294 + pc - the preconditioner context 295 - type - PC_DEFLATION_PRE, PC_DEFLATION_INIT, PC_DEFLATION_POST 296 297 Options Database Key: 298 . -pc_deflation_type <pre,init,post> 299 300 Level: intermediate 301 302 Concepts: Deflation preconditioner 303 304 .seealso: PCDeflationGetType() 305 @*/ 306 PetscErrorCode PCDeflationSetType(PC pc,PCDeflationType type) 307 { 308 PetscErrorCode ierr; 309 310 PetscFunctionBegin; 311 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 312 ierr = PetscTryMethod(pc,"PCDeflationSetType_C",(PC,PCDeflationType),(pc,type));CHKERRQ(ierr); 313 PetscFunctionReturn(0); 314 } 315 316 static PetscErrorCode PCDeflationGetType_Deflation(PC pc,PCDeflationType *type) 317 { 318 PC_Deflation *def = (PC_Deflation*)pc->data; 319 320 PetscFunctionBegin; 321 if (def->init) { 322 *type = PC_DEFLATION_INIT; 323 } else if (def->pre) { 324 *type = PC_DEFLATION_PRE; 325 } else { 326 *type = PC_DEFLATION_POST; 327 } 328 PetscFunctionReturn(0); 329 } 330 331 /*@ 332 PCDeflationGetType - Gets how the diagonal matrix is produced for the preconditioner 333 334 Not Collective on PC 335 336 Input Parameter: 337 . pc - the preconditioner context 338 339 Output Parameter: 340 - type - PC_DEFLATION_PRE, PC_DEFLATION_INIT, PC_DEFLATION_POST 341 342 Level: intermediate 343 344 Concepts: Deflation preconditioner 345 346 .seealso: PCDeflationSetType() 347 @*/ 348 PetscErrorCode PCDeflationGetType(PC pc,PCDeflationType *type) 349 { 350 PetscErrorCode ierr; 351 352 PetscFunctionBegin; 353 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 354 ierr = PetscUseMethod(pc,"PCDeflationGetType_C",(PC,PCDeflationType*),(pc,type));CHKERRQ(ierr); 355 PetscFunctionReturn(0); 356 } 357 358 static PetscErrorCode PCSetUp_Deflation(PC pc) 359 { 360 PC_Deflation *def = (PC_Deflation*)pc->data; 361 KSP innerksp; 362 PC pcinner; 363 Mat Amat,nextDef=NULL,*mats; 364 PetscInt i,m,red,size,commsize; 365 PetscBool match,flgspd,transp=PETSC_FALSE; 366 MatCompositeType ctype; 367 MPI_Comm comm; 368 const char *prefix; 369 PetscErrorCode ierr; 370 371 PetscFunctionBegin; 372 if (pc->setupcalled) PetscFunctionReturn(0); 373 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 374 ierr = PCGetOperators(pc,NULL,&Amat);CHKERRQ(ierr); 375 376 /* compute a deflation space */ 377 if (def->W || def->Wt) { 378 def->spacetype = PC_DEFLATION_SPACE_USER; 379 } else { 380 ierr = PCDeflationComputeSpace(pc);CHKERRQ(ierr); 381 } 382 383 /* nested deflation */ 384 if (def->W) { 385 ierr = PetscObjectTypeCompare((PetscObject)def->W,MATCOMPOSITE,&match);CHKERRQ(ierr); 386 if (match) { 387 ierr = MatCompositeGetType(def->W,&ctype);CHKERRQ(ierr); 388 ierr = MatCompositeGetNumberMat(def->W,&size);CHKERRQ(ierr); 389 } 390 } else { 391 ierr = MatCreateTranspose(def->Wt,&def->W);CHKERRQ(ierr); 392 ierr = PetscObjectTypeCompare((PetscObject)def->Wt,MATCOMPOSITE,&match);CHKERRQ(ierr); 393 if (match) { 394 ierr = MatCompositeGetType(def->Wt,&ctype);CHKERRQ(ierr); 395 ierr = MatCompositeGetNumberMat(def->Wt,&size);CHKERRQ(ierr); 396 } 397 transp = PETSC_TRUE; 398 } 399 if (match && ctype == MAT_COMPOSITE_MULTIPLICATIVE) { 400 if (!transp) { 401 if (def->nestedlvl < def->maxnestedlvl) { 402 ierr = PetscMalloc1(size,&mats);CHKERRQ(ierr); 403 for (i=0; i<size; i++) { 404 ierr = MatCompositeGetMat(def->W,i,&mats[i]);CHKERRQ(ierr); 405 } 406 size -= 1; 407 ierr = MatDestroy(&def->W);CHKERRQ(ierr); 408 def->W = mats[size]; 409 ierr = PetscObjectReference((PetscObject)mats[size]);CHKERRQ(ierr); 410 if (size > 1) { 411 ierr = MatCreateComposite(comm,size,mats,&nextDef);CHKERRQ(ierr); 412 ierr = MatCompositeSetType(nextDef,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr); 413 } else { 414 nextDef = mats[0]; 415 ierr = PetscObjectReference((PetscObject)mats[0]);CHKERRQ(ierr); 416 } 417 ierr = PetscFree(mats);CHKERRQ(ierr); 418 } else { 419 /* ierr = MatCompositeSetMergeType(def->W,MAT_COMPOSITE_MERGE_LEFT);CHKERRQ(ierr); */ 420 ierr = MatCompositeMerge(def->W);CHKERRQ(ierr); 421 } 422 } else { 423 if (def->nestedlvl < def->maxnestedlvl) { 424 ierr = PetscMalloc1(size,&mats);CHKERRQ(ierr); 425 for (i=0; i<size; i++) { 426 ierr = MatCompositeGetMat(def->Wt,i,&mats[i]);CHKERRQ(ierr); 427 } 428 size -= 1; 429 ierr = MatDestroy(&def->Wt);CHKERRQ(ierr); 430 def->Wt = mats[0]; 431 ierr = PetscObjectReference((PetscObject)mats[0]);CHKERRQ(ierr); 432 if (size > 1) { 433 ierr = MatCreateComposite(comm,size,&mats[1],&nextDef);CHKERRQ(ierr); 434 ierr = MatCompositeSetType(nextDef,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr); 435 } else { 436 nextDef = mats[1]; 437 ierr = PetscObjectReference((PetscObject)mats[1]);CHKERRQ(ierr); 438 } 439 ierr = PetscFree(mats);CHKERRQ(ierr); 440 } else { 441 /* ierr = MatCompositeSetMergeType(def->W,MAT_COMPOSITE_MERGE_LEFT);CHKERRQ(ierr); */ 442 ierr = MatCompositeMerge(def->Wt);CHKERRQ(ierr); 443 } 444 } 445 } 446 447 if (transp) { 448 ierr = MatDestroy(&def->W);CHKERRQ(ierr); 449 ierr = MatTranspose(def->Wt,MAT_INITIAL_MATRIX,&def->W);CHKERRQ(ierr); 450 } 451 452 453 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr); 454 455 /* setup coarse problem */ 456 if (!def->WtAWinv) { 457 ierr = MatGetSize(def->W,NULL,&m);CHKERRQ(ierr); 458 if (!def->WtAW) { 459 /* TODO add implicit product version ? */ 460 if (!def->AW) { 461 ierr = MatPtAP(Amat,def->W,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtAW);CHKERRQ(ierr); 462 } else { 463 ierr = MatTransposeMatMult(def->W,def->AW,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtAW);CHKERRQ(ierr); 464 } 465 /* TODO create MatInheritOption(Mat,MatOption) */ 466 ierr = MatGetOption(Amat,MAT_SPD,&flgspd);CHKERRQ(ierr); 467 ierr = MatSetOption(def->WtAW,MAT_SPD,flgspd);CHKERRQ(ierr); 468 #if defined(PETSC_USE_DEBUG) 469 /* Check WtAW is not sigular */ 470 PetscReal *norms; 471 ierr = PetscMalloc1(m,&norms);CHKERRQ(ierr); 472 ierr = MatGetColumnNorms(def->WtAW,NORM_INFINITY,norms);CHKERRQ(ierr); 473 for (i=0; i<m; i++) { 474 if (norms[i] < 100*PETSC_MACHINE_EPSILON) { 475 SETERRQ1(comm,PETSC_ERR_SUP,"Column %D of W is in kernel of A.",i); 476 } 477 } 478 ierr = PetscFree(norms);CHKERRQ(ierr); 479 #endif 480 } else { 481 ierr = MatGetOption(def->WtAW,MAT_SPD,&flgspd);CHKERRQ(ierr); 482 } 483 /* TODO use MATINV */ 484 ierr = KSPCreate(comm,&def->WtAWinv);CHKERRQ(ierr); 485 ierr = KSPSetOperators(def->WtAWinv,def->WtAW,def->WtAW);CHKERRQ(ierr); 486 ierr = KSPGetPC(def->WtAWinv,&pcinner);CHKERRQ(ierr); 487 /* Setup KSP and PC */ 488 if (nextDef) { /* next level for multilevel deflation */ 489 innerksp = def->WtAWinv; 490 /* set default KSPtype */ 491 if (!def->ksptype) { 492 def->ksptype = KSPFGMRES; 493 if (flgspd) { /* SPD system */ 494 def->ksptype = KSPFCG; 495 } 496 } 497 ierr = KSPSetType(innerksp,def->ksptype);CHKERRQ(ierr); /* TODO iherit from KSP */ 498 ierr = PCSetType(pcinner,PCDEFLATION);CHKERRQ(ierr); /* TODO create coarse preconditinoner M_c = WtMW ? */ 499 ierr = PCDeflationSetSpace(pcinner,nextDef,transp);CHKERRQ(ierr); 500 ierr = PCDeflationSetLvl_Deflation(pcinner,def->nestedlvl+1,def->maxnestedlvl);CHKERRQ(ierr); 501 /* inherit options TODO if not set */ 502 ((PC_Deflation*)(pcinner->data))->ksptype = def->ksptype; 503 ((PC_Deflation*)(pcinner->data))->correct = def->correct; 504 ((PC_Deflation*)(pcinner->data))->adaptiveconv = def->adaptiveconv; 505 ((PC_Deflation*)(pcinner->data))->adaptiveconst = def->adaptiveconst; 506 ierr = MatDestroy(&nextDef);CHKERRQ(ierr); 507 } else { /* the last level */ 508 ierr = KSPSetType(def->WtAWinv,KSPPREONLY);CHKERRQ(ierr); 509 ierr = PCSetType(pcinner,PCTELESCOPE);CHKERRQ(ierr); 510 /* ugly hack to not have overwritten PCTELESCOPE */ 511 if (prefix) { 512 ierr = KSPSetOptionsPrefix(def->WtAWinv,prefix);CHKERRQ(ierr); 513 } 514 ierr = KSPAppendOptionsPrefix(def->WtAWinv,"tel_");CHKERRQ(ierr); 515 ierr = PCSetFromOptions(pcinner);CHKERRQ(ierr); 516 /* Reduction factor choice */ 517 red = def->reductionfact; 518 if (red < 0) { 519 ierr = MPI_Comm_size(comm,&commsize);CHKERRQ(ierr); 520 red = ceil((float)commsize/ceil((float)m/commsize)); 521 ierr = PetscObjectTypeCompareAny((PetscObject)(def->WtAW),&match,MATSEQDENSE,MATMPIDENSE,MATDENSE,"");CHKERRQ(ierr); 522 if (match) red = commsize; 523 ierr = PetscInfo1(pc,"Auto choosing reduction factor %D\n",red);CHKERRQ(ierr); /* TODO add level? */ 524 } 525 ierr = PCTelescopeSetReductionFactor(pcinner,red);CHKERRQ(ierr); 526 ierr = PCSetUp(pcinner);CHKERRQ(ierr); 527 ierr = PCTelescopeGetKSP(pcinner,&innerksp);CHKERRQ(ierr); 528 if (innerksp) { 529 ierr = KSPGetPC(innerksp,&pcinner);CHKERRQ(ierr); 530 /* TODO Cholesky if flgspd? */ 531 ierr = PCSetType(pcinner,PCLU);CHKERRQ(ierr); 532 //TODO remove explicit matSolverPackage 533 if (commsize == red) { 534 ierr = PCFactorSetMatSolverType(pcinner,MATSOLVERSUPERLU);CHKERRQ(ierr); 535 } else { 536 ierr = PCFactorSetMatSolverType(pcinner,MATSOLVERSUPERLU_DIST);CHKERRQ(ierr); 537 } 538 } 539 } 540 541 if (innerksp) { 542 /* TODO use def_[lvl]_ if lvl > 0? */ 543 if (prefix) { 544 ierr = KSPSetOptionsPrefix(innerksp,prefix);CHKERRQ(ierr); 545 } 546 ierr = KSPAppendOptionsPrefix(innerksp,"def_");CHKERRQ(ierr); 547 ierr = KSPSetFromOptions(innerksp);CHKERRQ(ierr); 548 ierr = KSPSetUp(innerksp);CHKERRQ(ierr); 549 } 550 /* TODO: check if WtAWinv is KSP and move following from this if */ 551 //if (def->adaptiveconv) { 552 // PetscReal *rnorm; 553 // PetscNew(&rnorm); 554 // ierr = KSPSetConvergenceTest(def->WtAWinv,KSPDCGConvergedAdaptive_DCG,rnorm,NULL);CHKERRQ(ierr); 555 //} 556 } 557 ierr = KSPSetFromOptions(def->WtAWinv);CHKERRQ(ierr); 558 ierr = KSPSetUp(def->WtAWinv);CHKERRQ(ierr); 559 560 if (!def->pc) { 561 ierr = PCCreate(comm,&def->pc);CHKERRQ(ierr); 562 ierr = PCSetOperators(def->pc,Amat,Amat);CHKERRQ(ierr); 563 ierr = PCSetType(def->pc,PCNONE);CHKERRQ(ierr); 564 if (prefix) { 565 ierr = PCSetOptionsPrefix(def->pc,prefix);CHKERRQ(ierr); 566 } 567 ierr = PCAppendOptionsPrefix(def->pc,"def_pc_");CHKERRQ(ierr); 568 ierr = PCSetFromOptions(def->pc);CHKERRQ(ierr); 569 ierr = PCSetUp(def->pc);CHKERRQ(ierr); 570 } 571 572 /* create work vecs */ 573 ierr = MatCreateVecs(Amat,NULL,&def->work);CHKERRQ(ierr); 574 ierr = KSPCreateVecs(def->WtAWinv,2,&def->workcoarse,0,NULL);CHKERRQ(ierr); 575 PetscFunctionReturn(0); 576 } 577 578 static PetscErrorCode PCReset_Deflation(PC pc) 579 { 580 PC_Deflation *def = (PC_Deflation*)pc->data; 581 PetscErrorCode ierr; 582 583 PetscFunctionBegin; 584 ierr = VecDestroy(&def->work);CHKERRQ(ierr); 585 ierr = VecDestroyVecs(2,&def->workcoarse);CHKERRQ(ierr); 586 ierr = MatDestroy(&def->W);CHKERRQ(ierr); 587 ierr = MatDestroy(&def->Wt);CHKERRQ(ierr); 588 ierr = MatDestroy(&def->AW);CHKERRQ(ierr); 589 ierr = MatDestroy(&def->WtAW);CHKERRQ(ierr); 590 ierr = KSPDestroy(&def->WtAWinv);CHKERRQ(ierr); 591 ierr = PCDestroy(&def->pc);CHKERRQ(ierr); 592 PetscFunctionReturn(0); 593 } 594 595 /* 596 PCDestroy_Deflation - Destroys the private context for the Deflation preconditioner 597 that was created with PCCreate_Deflation(). 598 599 Input Parameter: 600 . pc - the preconditioner context 601 602 Application Interface Routine: PCDestroy() 603 */ 604 static PetscErrorCode PCDestroy_Deflation(PC pc) 605 { 606 PetscErrorCode ierr; 607 608 PetscFunctionBegin; 609 ierr = PCReset_Deflation(pc);CHKERRQ(ierr); 610 ierr = PetscFree(pc->data);CHKERRQ(ierr); 611 PetscFunctionReturn(0); 612 } 613 614 static PetscErrorCode PCView_Deflation(PC pc,PetscViewer viewer) 615 { 616 PC_Deflation *def = (PC_Deflation*)pc->data; 617 PetscBool iascii; 618 PetscErrorCode ierr; 619 620 PetscFunctionBegin; 621 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 622 if (iascii) { 623 //if (cg->adaptiveconv) { 624 // ierr = PetscViewerASCIIPrintf(viewer," DCG: using adaptive precision for inner solve with C=%.1e\n",cg->adaptiveconst);CHKERRQ(ierr); 625 //} 626 if (def->correct) { 627 ierr = PetscViewerASCIIPrintf(viewer," Using CP correction\n");CHKERRQ(ierr); 628 } 629 if (!def->nestedlvl) { 630 ierr = PetscViewerASCIIPrintf(viewer," Deflation space type: %s\n",PCDeflationSpaceTypes[def->spacetype]);CHKERRQ(ierr); 631 ierr = PetscViewerASCIIPrintf(viewer," DCG %s\n",def->extendsp ? "extended" : "truncated");CHKERRQ(ierr); 632 } 633 634 ierr = PetscViewerASCIIPrintf(viewer,"--- Additional PC\n");CHKERRQ(ierr); 635 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 636 ierr = PCView(def->pc,viewer);CHKERRQ(ierr); 637 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 638 639 ierr = PetscViewerASCIIPrintf(viewer,"--- Coarse solver\n");CHKERRQ(ierr); 640 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 641 ierr = KSPView(def->WtAWinv,viewer);CHKERRQ(ierr); 642 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 643 } 644 PetscFunctionReturn(0); 645 } 646 647 static PetscErrorCode PCSetFromOptions_Deflation(PetscOptionItems *PetscOptionsObject,PC pc) 648 { 649 PC_Deflation *def = (PC_Deflation*)pc->data; 650 PetscErrorCode ierr; 651 PetscBool flg; 652 PCDeflationType deflt,type; 653 654 PetscFunctionBegin; 655 ierr = PCDeflationGetType(pc,&deflt);CHKERRQ(ierr); 656 ierr = PetscOptionsHead(PetscOptionsObject,"Deflation options");CHKERRQ(ierr); 657 ierr = PetscOptionsEnum("-pc_deflation_type","Determine type of deflation","PCDeflationSetType",PCDeflationTypes,(PetscEnum)deflt,(PetscEnum*)&type,&flg);CHKERRQ(ierr); 658 if (flg) { 659 ierr = PCDeflationSetType(pc,type);CHKERRQ(ierr); 660 } 661 ierr = PetscOptionsEnum("-pc_deflation_compute_space","Compute deflation space","PCDeflationSetSpace",PCDeflationSpaceTypes,(PetscEnum)def->spacetype,(PetscEnum*)&def->spacetype,NULL);CHKERRQ(ierr); 662 ierr = PetscOptionsInt("-pc_deflation_compute_space_size","Set size of the deflation space to compute","PCDeflationSetSpace",def->spacesize,&def->spacesize,NULL);CHKERRQ(ierr); 663 ierr = PetscOptionsBool("-pc_deflation_space_extend","Extend deflation space instead of truncating (wavelets)","PCDeflation",def->extendsp,&def->extendsp,NULL);CHKERRQ(ierr); 664 //TODO add set function and fix manpages 665 ierr = PetscOptionsBool("-pc_deflation_initdef","Use only initialization step - Initdef","PCDeflation",def->init,&def->init,NULL);CHKERRQ(ierr); 666 ierr = PetscOptionsBool("-pc_deflation_correct","Add Qr to descent direction","PCDeflation",def->correct,&def->correct,NULL);CHKERRQ(ierr); 667 ierr = PetscOptionsBool("-pc_deflation_adaptive","Adaptive stopping criteria","PCDeflation",def->adaptiveconv,&def->adaptiveconv,NULL);CHKERRQ(ierr); 668 ierr = PetscOptionsReal("-pc_deflation_adaptive_const","Adaptive stopping criteria constant","PCDeflation",def->adaptiveconst,&def->adaptiveconst,NULL);CHKERRQ(ierr); 669 ierr = PetscOptionsInt("-pc_deflation_redfact","Reduction factor for coarse problem solution","PCDeflation",def->reductionfact,&def->reductionfact,NULL);CHKERRQ(ierr); 670 ierr = PetscOptionsInt("-pc_deflation_max_nested_lvl","Maximum of nested deflation levels","PCDeflation",def->maxnestedlvl,&def->maxnestedlvl,NULL);CHKERRQ(ierr); 671 ierr = PetscOptionsTail();CHKERRQ(ierr); 672 PetscFunctionReturn(0); 673 } 674 675 /*MC 676 PCDEFLATION - Deflation preconditioner shifts part of the spectrum to zero (deflates) 677 or to a predefined value 678 679 Options Database Key: 680 + -pc_deflation_type <init,pre,post> - selects approach to deflation (default: pre) 681 - -pc_jacobi_abs - use the absolute value of the diagonal entry 682 683 Level: beginner 684 685 Notes: 686 todo 687 688 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 689 PCDeflationSetType(), PCDeflationSetSpace() 690 M*/ 691 692 PETSC_EXTERN PetscErrorCode PCCreate_Deflation(PC pc) 693 { 694 PC_Deflation *def; 695 PetscErrorCode ierr; 696 697 PetscFunctionBegin; 698 ierr = PetscNewLog(pc,&def);CHKERRQ(ierr); 699 pc->data = (void*)def; 700 701 def->init = PETSC_FALSE; 702 def->pre = PETSC_TRUE; 703 def->correct = PETSC_FALSE; 704 def->truenorm = PETSC_TRUE; 705 def->reductionfact = -1; 706 def->spacetype = PC_DEFLATION_SPACE_HAAR; 707 def->spacesize = 1; 708 def->extendsp = PETSC_FALSE; 709 def->nestedlvl = 0; 710 def->maxnestedlvl = 0; 711 def->adaptiveconv = PETSC_FALSE; 712 def->adaptiveconst = 1.0; 713 714 /* 715 Set the pointers for the functions that are provided above. 716 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 717 are called, they will automatically call these functions. Note we 718 choose not to provide a couple of these functions since they are 719 not needed. 720 */ 721 pc->ops->apply = PCApply_Deflation; 722 pc->ops->applytranspose = PCApply_Deflation; 723 pc->ops->presolve = PCPreSolve_Deflation; 724 pc->ops->postsolve = 0; 725 pc->ops->setup = PCSetUp_Deflation; 726 pc->ops->reset = PCReset_Deflation; 727 pc->ops->destroy = PCDestroy_Deflation; 728 pc->ops->setfromoptions = PCSetFromOptions_Deflation; 729 pc->ops->view = PCView_Deflation; 730 pc->ops->applyrichardson = 0; 731 pc->ops->applysymmetricleft = 0; 732 pc->ops->applysymmetricright = 0; 733 734 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetType_C",PCDeflationSetType_Deflation);CHKERRQ(ierr); 735 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationGetType_C",PCDeflationGetType_Deflation);CHKERRQ(ierr); 736 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetSpace_C",PCDeflationSetSpace_Deflation);CHKERRQ(ierr); 737 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetLvl_C",PCDeflationSetLvl_Deflation);CHKERRQ(ierr); 738 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetPC_C",PCDeflationSetPC_Deflation);CHKERRQ(ierr); 739 PetscFunctionReturn(0); 740 } 741 742