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 <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 52 53 const char *const PCDeflationTypes[] = {"INIT","PRE","POST","PCDeflationType","PC_DEFLATION_",0}; 54 55 /* 56 Private context (data structure) for the deflation preconditioner. 57 */ 58 typedef struct { 59 PetscBool init; /* do only init step - error correction of direction is omitted */ 60 PetscBool pre; /* start with x0 being the solution in the deflation space */ 61 PetscBool correct; /* add CP (Qr) correction to descent direction */ 62 PetscBool truenorm; 63 PetscBool adaptiveconv; 64 PetscReal adaptiveconst; 65 PetscInt reductionfact; 66 Mat W,Wt,AW,WtAW; /* deflation space, coarse problem mats */ 67 KSP WtAWinv; /* deflation coarse problem */ 68 Vec *work; 69 70 PCDEFLATIONSpaceType spacetype; 71 PetscInt spacesize; 72 PetscInt nestedlvl; 73 PetscInt maxnestedlvl; 74 PetscBool extendsp; 75 } PC_Deflation; 76 77 static PetscErrorCode PCDeflationSetType_Deflation(PC pc,PCDeflationType type) 78 { 79 PC_Deflation *def = (PC_Deflation*)pc->data; 80 81 PetscFunctionBegin; 82 def->init = PETSC_FALSE; 83 def->pre = PETSC_FALSE; 84 if (type == PC_DEFLATION_INIT) { 85 def->init = PETSC_TRUE; 86 def->pre = PETSC_TRUE; 87 } else if (type == PC_DEFLATION_PRE) { 88 def->pre = PETSC_TRUE; 89 } 90 PetscFunctionReturn(0); 91 } 92 93 static PetscErrorCode PCDeflationGetType_Deflation(PC pc,PCDeflationType *type) 94 { 95 PC_Deflation *def = (PC_Deflation*)pc->data; 96 97 PetscFunctionBegin; 98 if (def->init) { 99 *type = PC_DEFLATION_INIT; 100 } else if (def->pre) { 101 *type = PC_DEFLATION_PRE; 102 } else { 103 *type = PC_DEFLATION_POST; 104 } 105 PetscFunctionReturn(0); 106 } 107 108 static PetscErrorCode PCDeflationSetSpace_Deflation(PC pc,Mat W,PetscBool transpose) 109 { 110 PC_Deflation *def = (PC_Deflation*)pc->data; 111 PetscErrorCode ierr; 112 113 PetscFunctionBegin; 114 if (transpose) { 115 def->Wt = W; 116 def->W = NULL; 117 } else { 118 def->W = W; 119 } 120 ierr = PetscObjectReference((PetscObject)W);CHKERRQ(ierr); 121 PetscFunctionReturn(0); 122 } 123 124 /* TODO create PCDeflationSetSpaceTranspose? */ 125 /*@ 126 PCDeflationSetSpace - Set deflation space matrix (or its transpose). 127 128 Logically Collective on PC 129 130 Input Parameters: 131 + pc - the preconditioner context 132 . W - deflation matrix 133 - tranpose - indicates that W is an explicit transpose of the deflation matrix 134 135 Level: intermediate 136 137 .seealso: PCDEFLATION 138 @*/ 139 PetscErrorCode PCDeflationSetSpace(PC pc,Mat W,PetscBool transpose) 140 { 141 PetscErrorCode ierr; 142 143 PetscFunctionBegin; 144 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 145 PetscValidHeaderSpecific(W,MAT_CLASSID,2); 146 PetscValidLogicalCollectiveBool(pc,transpose,3); 147 ierr = PetscTryMethod(pc,"PCDeflationSetSpace_C",(PC,Mat,PetscBool),(pc,W,transpose));CHKERRQ(ierr); 148 PetscFunctionReturn(0); 149 } 150 151 152 /* -------------------------------------------------------------------------- */ 153 /* 154 PCSetUp_Deflation - Prepares for the use of the Deflation preconditioner 155 by setting data structures and options. 156 157 Input Parameter: 158 . pc - the preconditioner context 159 160 Application Interface Routine: PCSetUp() 161 162 Notes: 163 The interface routine PCSetUp() is not usually called directly by 164 the user, but instead is called by PCApply() if necessary. 165 */ 166 static PetscErrorCode PCSetUp_Deflation(PC pc) 167 { 168 PC_Deflation *jac = (PC_Deflation*)pc->data; 169 Vec diag,diagsqrt; 170 PetscErrorCode ierr; 171 PetscInt n,i; 172 PetscScalar *x; 173 PetscBool zeroflag = PETSC_FALSE; 174 175 PetscFunctionBegin; 176 /* 177 For most preconditioners the code would begin here something like 178 179 if (pc->setupcalled == 0) { allocate space the first time this is ever called 180 ierr = MatCreateVecs(pc->mat,&jac->diag);CHKERRQ(ierr); 181 PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->diag); 182 } 183 184 But for this preconditioner we want to support use of both the matrix' diagonal 185 elements (for left or right preconditioning) and square root of diagonal elements 186 (for symmetric preconditioning). Hence we do not allocate space here, since we 187 don't know at this point which will be needed (diag and/or diagsqrt) until the user 188 applies the preconditioner, and we don't want to allocate BOTH unless we need 189 them both. Thus, the diag and diagsqrt are allocated in PCSetUp_Deflation_NonSymmetric() 190 and PCSetUp_Deflation_Symmetric(), respectively. 191 */ 192 193 /* 194 Here we set up the preconditioner; that is, we copy the diagonal values from 195 the matrix and put them into a format to make them quick to apply as a preconditioner. 196 */ 197 if (zeroflag) { 198 ierr = PetscInfo(pc,"Zero detected in diagonal of matrix, using 1 at those locations\n");CHKERRQ(ierr); 199 } 200 PetscFunctionReturn(0); 201 } 202 /* -------------------------------------------------------------------------- */ 203 /* 204 PCApply_Deflation - Applies the Deflation preconditioner to a vector. 205 206 Input Parameters: 207 . pc - the preconditioner context 208 . x - input vector 209 210 Output Parameter: 211 . y - output vector 212 213 Application Interface Routine: PCApply() 214 */ 215 static PetscErrorCode PCApply_Deflation(PC pc,Vec x,Vec y) 216 { 217 PC_Deflation *jac = (PC_Deflation*)pc->data; 218 PetscErrorCode ierr; 219 220 PetscFunctionBegin; 221 PetscFunctionReturn(0); 222 } 223 /* -------------------------------------------------------------------------- */ 224 static PetscErrorCode PCReset_Deflation(PC pc) 225 { 226 PC_Deflation *jac = (PC_Deflation*)pc->data; 227 PetscErrorCode ierr; 228 229 PetscFunctionBegin; 230 PetscFunctionReturn(0); 231 } 232 233 /* 234 PCDestroy_Deflation - Destroys the private context for the Deflation preconditioner 235 that was created with PCCreate_Deflation(). 236 237 Input Parameter: 238 . pc - the preconditioner context 239 240 Application Interface Routine: PCDestroy() 241 */ 242 static PetscErrorCode PCDestroy_Deflation(PC pc) 243 { 244 PetscErrorCode ierr; 245 246 PetscFunctionBegin; 247 ierr = PCReset_Deflation(pc);CHKERRQ(ierr); 248 249 /* 250 Free the private data structure that was hanging off the PC 251 */ 252 ierr = PetscFree(pc->data);CHKERRQ(ierr); 253 PetscFunctionReturn(0); 254 } 255 256 static PetscErrorCode PCSetFromOptions_Deflation(PetscOptionItems *PetscOptionsObject,PC pc) 257 { 258 PC_Deflation *jac = (PC_Deflation*)pc->data; 259 PetscErrorCode ierr; 260 PetscBool flg; 261 PCDeflationType deflt,type; 262 263 PetscFunctionBegin; 264 ierr = PCDeflationGetType(pc,&deflt);CHKERRQ(ierr); 265 ierr = PetscOptionsHead(PetscOptionsObject,"Deflation options");CHKERRQ(ierr); 266 ierr = PetscOptionsEnum("-pc_jacobi_type","How to construct diagonal matrix","PCDeflationSetType",PCDeflationTypes,(PetscEnum)deflt,(PetscEnum*)&type,&flg);CHKERRQ(ierr); 267 if (flg) { 268 ierr = PCDeflationSetType(pc,type);CHKERRQ(ierr); 269 } 270 ierr = PetscOptionsTail();CHKERRQ(ierr); 271 PetscFunctionReturn(0); 272 } 273 274 /*MC 275 PCDEFLATION - Deflation preconditioner shifts part of the spectrum to zero (deflates) 276 or to a predefined value 277 278 Options Database Key: 279 + -pc_deflation_type <init,pre,post> - selects approach to deflation (default: pre) 280 - -pc_jacobi_abs - use the absolute value of the diagonal entry 281 282 Level: beginner 283 284 Concepts: Deflation, preconditioners 285 286 Notes: 287 todo 288 289 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 290 PCDeflationSetType(), PCDeflationSetSpace() 291 M*/ 292 293 PETSC_EXTERN PetscErrorCode PCCreate_Deflation(PC pc) 294 { 295 PC_Deflation *def; 296 PetscErrorCode ierr; 297 298 PetscFunctionBegin; 299 ierr = PetscNewLog(pc,&def);CHKERRQ(ierr); 300 pc->data = (void*)def; 301 302 def->init = PETSC_FALSE; 303 def->pre = PETSC_TRUE; 304 def->correct = PETSC_FALSE; 305 def->truenorm = PETSC_TRUE; 306 def->reductionfact = -1; 307 def->spacetype = PC_DEFLATION_SPACE_HAAR; 308 def->spacesize = 1; 309 def->extendsp = PETSC_FALSE; 310 def->nestedlvl = 0; 311 def->maxnestedlvl = 0; 312 def->adaptiveconv = PETSC_FALSE; 313 def->adaptiveconst = 1.0; 314 315 /* 316 Set the pointers for the functions that are provided above. 317 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 318 are called, they will automatically call these functions. Note we 319 choose not to provide a couple of these functions since they are 320 not needed. 321 */ 322 pc->ops->apply = PCApply_Deflation; 323 pc->ops->applytranspose = PCApply_Deflation; 324 pc->ops->setup = PCSetUp_Deflation; 325 pc->ops->reset = PCReset_Deflation; 326 pc->ops->destroy = PCDestroy_Deflation; 327 pc->ops->setfromoptions = PCSetFromOptions_Deflation; 328 pc->ops->view = 0; 329 pc->ops->applyrichardson = 0; 330 pc->ops->applysymmetricleft = 0; 331 pc->ops->applysymmetricright = 0; 332 333 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetType_C",PCDeflationSetType_Deflation);CHKERRQ(ierr); 334 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationGetType_C",PCDeflationGetType_Deflation);CHKERRQ(ierr); 335 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetSpace_C",PCDeflationSetSpace_Deflation);CHKERRQ(ierr); 336 PetscFunctionReturn(0); 337 } 338 339