137eeb815SJakub Kruzik 237eeb815SJakub Kruzik /* -------------------------------------------------------------------- 337eeb815SJakub Kruzik 437eeb815SJakub Kruzik This file implements a Deflation preconditioner in PETSc as part of PC. 537eeb815SJakub Kruzik You can use this as a starting point for implementing your own 637eeb815SJakub Kruzik preconditioner that is not provided with PETSc. (You might also consider 737eeb815SJakub Kruzik just using PCSHELL) 837eeb815SJakub Kruzik 937eeb815SJakub Kruzik The following basic routines are required for each preconditioner. 1037eeb815SJakub Kruzik PCCreate_XXX() - Creates a preconditioner context 1137eeb815SJakub Kruzik PCSetFromOptions_XXX() - Sets runtime options 1237eeb815SJakub Kruzik PCApply_XXX() - Applies the preconditioner 1337eeb815SJakub Kruzik PCDestroy_XXX() - Destroys the preconditioner context 1437eeb815SJakub Kruzik where the suffix "_XXX" denotes a particular implementation, in 1537eeb815SJakub Kruzik this case we use _Deflation (e.g., PCCreate_Deflation, PCApply_Deflation). 1637eeb815SJakub Kruzik These routines are actually called via the common user interface 1737eeb815SJakub Kruzik routines PCCreate(), PCSetFromOptions(), PCApply(), and PCDestroy(), 1837eeb815SJakub Kruzik so the application code interface remains identical for all 1937eeb815SJakub Kruzik preconditioners. 2037eeb815SJakub Kruzik 2137eeb815SJakub Kruzik Another key routine is: 2237eeb815SJakub Kruzik PCSetUp_XXX() - Prepares for the use of a preconditioner 2337eeb815SJakub Kruzik by setting data structures and options. The interface routine PCSetUp() 2437eeb815SJakub Kruzik is not usually called directly by the user, but instead is called by 2537eeb815SJakub Kruzik PCApply() if necessary. 2637eeb815SJakub Kruzik 2737eeb815SJakub Kruzik Additional basic routines are: 2837eeb815SJakub Kruzik PCView_XXX() - Prints details of runtime options that 2937eeb815SJakub Kruzik have actually been used. 3037eeb815SJakub Kruzik These are called by application codes via the interface routines 3137eeb815SJakub Kruzik PCView(). 3237eeb815SJakub Kruzik 3337eeb815SJakub Kruzik The various types of solvers (preconditioners, Krylov subspace methods, 3437eeb815SJakub Kruzik nonlinear solvers, timesteppers) are all organized similarly, so the 3537eeb815SJakub Kruzik above description applies to these categories also. One exception is 3637eeb815SJakub Kruzik that the analogues of PCApply() for these components are KSPSolve(), 3737eeb815SJakub Kruzik SNESSolve(), and TSSolve(). 3837eeb815SJakub Kruzik 3937eeb815SJakub Kruzik Additional optional functionality unique to preconditioners is left and 4037eeb815SJakub Kruzik right symmetric preconditioner application via PCApplySymmetricLeft() 4137eeb815SJakub Kruzik and PCApplySymmetricRight(). The Deflation implementation is 4237eeb815SJakub Kruzik PCApplySymmetricLeftOrRight_Deflation(). 4337eeb815SJakub Kruzik 4437eeb815SJakub Kruzik -------------------------------------------------------------------- */ 4537eeb815SJakub Kruzik 4637eeb815SJakub Kruzik /* 4737eeb815SJakub Kruzik Include files needed for the Deflation preconditioner: 4837eeb815SJakub Kruzik pcimpl.h - private include file intended for use by all preconditioners 4937eeb815SJakub Kruzik */ 5037eeb815SJakub Kruzik 5137eeb815SJakub Kruzik #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 5237eeb815SJakub Kruzik 5337eeb815SJakub Kruzik const char *const PCDeflationTypes[] = {"INIT","PRE","POST","PCDeflationType","PC_DEFLATION_",0}; 5437eeb815SJakub Kruzik 5537eeb815SJakub Kruzik /* 5637eeb815SJakub Kruzik Private context (data structure) for the deflation preconditioner. 5737eeb815SJakub Kruzik */ 5837eeb815SJakub Kruzik typedef struct { 5937eeb815SJakub Kruzik PetscBool init; /* do only init step - error correction of direction is omitted */ 6037eeb815SJakub Kruzik PetscBool pre; /* start with x0 being the solution in the deflation space */ 6137eeb815SJakub Kruzik PetscBool correct; /* add CP (Qr) correction to descent direction */ 6237eeb815SJakub Kruzik PetscBool truenorm; 6337eeb815SJakub Kruzik PetscBool adaptiveconv; 6437eeb815SJakub Kruzik PetscReal adaptiveconst; 6537eeb815SJakub Kruzik PetscInt reductionfact; 6637eeb815SJakub Kruzik Mat W,Wt,AW,WtAW; /* deflation space, coarse problem mats */ 6737eeb815SJakub Kruzik KSP WtAWinv; /* deflation coarse problem */ 68409a357bSJakub Kruzik KSPType ksptype; 69f8bf229cSJakub Kruzik Vec work; 70f8bf229cSJakub Kruzik Vec *workcoarse; 7137eeb815SJakub Kruzik 72409a357bSJakub Kruzik PCDeflationSpaceType spacetype; 7337eeb815SJakub Kruzik PetscInt spacesize; 7437eeb815SJakub Kruzik PetscInt nestedlvl; 7537eeb815SJakub Kruzik PetscInt maxnestedlvl; 7637eeb815SJakub Kruzik PetscBool extendsp; 7737eeb815SJakub Kruzik } PC_Deflation; 7837eeb815SJakub Kruzik 79e662bc50SJakub Kruzik static PetscErrorCode PCDeflationSetSpace_Deflation(PC pc,Mat W,PetscBool transpose) 80e662bc50SJakub Kruzik { 81e662bc50SJakub Kruzik PC_Deflation *def = (PC_Deflation*)pc->data; 82e662bc50SJakub Kruzik PetscErrorCode ierr; 83e662bc50SJakub Kruzik 84e662bc50SJakub Kruzik PetscFunctionBegin; 85e662bc50SJakub Kruzik if (transpose) { 86e662bc50SJakub Kruzik def->Wt = W; 87e662bc50SJakub Kruzik def->W = NULL; 88e662bc50SJakub Kruzik } else { 89e662bc50SJakub Kruzik def->W = W; 90e662bc50SJakub Kruzik } 91e662bc50SJakub Kruzik ierr = PetscObjectReference((PetscObject)W);CHKERRQ(ierr); 92e662bc50SJakub Kruzik PetscFunctionReturn(0); 93e662bc50SJakub Kruzik } 94e662bc50SJakub Kruzik 95e662bc50SJakub Kruzik /* TODO create PCDeflationSetSpaceTranspose? */ 96e662bc50SJakub Kruzik /*@ 97e662bc50SJakub Kruzik PCDeflationSetSpace - Set deflation space matrix (or its transpose). 98e662bc50SJakub Kruzik 99e662bc50SJakub Kruzik Logically Collective on PC 100e662bc50SJakub Kruzik 101e662bc50SJakub Kruzik Input Parameters: 102e662bc50SJakub Kruzik + pc - the preconditioner context 103e662bc50SJakub Kruzik . W - deflation matrix 104e662bc50SJakub Kruzik - tranpose - indicates that W is an explicit transpose of the deflation matrix 105e662bc50SJakub Kruzik 106e662bc50SJakub Kruzik Level: intermediate 107e662bc50SJakub Kruzik 108e662bc50SJakub Kruzik .seealso: PCDEFLATION 109e662bc50SJakub Kruzik @*/ 110e662bc50SJakub Kruzik PetscErrorCode PCDeflationSetSpace(PC pc,Mat W,PetscBool transpose) 111e662bc50SJakub Kruzik { 112e662bc50SJakub Kruzik PetscErrorCode ierr; 113e662bc50SJakub Kruzik 114e662bc50SJakub Kruzik PetscFunctionBegin; 115e662bc50SJakub Kruzik PetscValidHeaderSpecific(pc,PC_CLASSID,1); 116e662bc50SJakub Kruzik PetscValidHeaderSpecific(W,MAT_CLASSID,2); 117e662bc50SJakub Kruzik PetscValidLogicalCollectiveBool(pc,transpose,3); 118e662bc50SJakub Kruzik ierr = PetscTryMethod(pc,"PCDeflationSetSpace_C",(PC,Mat,PetscBool),(pc,W,transpose));CHKERRQ(ierr); 119e662bc50SJakub Kruzik PetscFunctionReturn(0); 120e662bc50SJakub Kruzik } 121e662bc50SJakub Kruzik 1225c0c31c5SJakub Kruzik static PetscErrorCode PCDeflationSetLvl_Deflation(PC pc,PetscInt current,PetscInt max) 1235c0c31c5SJakub Kruzik { 1245c0c31c5SJakub Kruzik PC_Deflation *def = (PC_Deflation*)pc->data; 1255c0c31c5SJakub Kruzik 1265c0c31c5SJakub Kruzik PetscFunctionBegin; 1275c0c31c5SJakub Kruzik def->nestedlvl = current; 1285c0c31c5SJakub Kruzik def->maxnestedlvl = max; 1295c0c31c5SJakub Kruzik PetscFunctionReturn(0); 1305c0c31c5SJakub Kruzik } 1315c0c31c5SJakub Kruzik 1325c0c31c5SJakub Kruzik /*@ 1335c0c31c5SJakub Kruzik PCDeflationSetMaxLvl - Set maximum level of deflation. 1345c0c31c5SJakub Kruzik 1355c0c31c5SJakub Kruzik Logically Collective on PC 1365c0c31c5SJakub Kruzik 1375c0c31c5SJakub Kruzik Input Parameters: 1385c0c31c5SJakub Kruzik + pc - the preconditioner context 1395c0c31c5SJakub Kruzik . max - maximum deflation level 1405c0c31c5SJakub Kruzik 1415c0c31c5SJakub Kruzik Level: intermediate 1425c0c31c5SJakub Kruzik 1435c0c31c5SJakub Kruzik .seealso: PCDEFLATION 1445c0c31c5SJakub Kruzik @*/ 1455c0c31c5SJakub Kruzik PetscErrorCode PCDeflationSetMaxLvl(PC pc,PetscInt max) 1465c0c31c5SJakub Kruzik { 1475c0c31c5SJakub Kruzik PetscErrorCode ierr; 1485c0c31c5SJakub Kruzik 1495c0c31c5SJakub Kruzik PetscFunctionBegin; 1505c0c31c5SJakub Kruzik PetscValidLogicalCollectiveInt(pc,max,2); 1515c0c31c5SJakub Kruzik ierr = PetscTryMethod(pc,"PCDeflationSetLvl_C",(PC,PetscInt,PetscInt),(pc,0,max));CHKERRQ(ierr); 1525c0c31c5SJakub Kruzik PetscFunctionReturn(0); 1535c0c31c5SJakub Kruzik } 154e662bc50SJakub Kruzik 15537eeb815SJakub Kruzik /* 156*b27c8092SJakub Kruzik TODO CP corection? 157*b27c8092SJakub Kruzik x <- x + W*(W'*A*W)^{-1}*W'*r = x + Q*r 158*b27c8092SJakub Kruzik */ 159*b27c8092SJakub Kruzik static PetscErrorCode PCPreSolve_Deflation(PC pc,KSP ksp,Vec b, Vec x) 160*b27c8092SJakub Kruzik { 161*b27c8092SJakub Kruzik PC_Deflation *def = (PC_Deflation*)pc->data; 162*b27c8092SJakub Kruzik Mat A; 163*b27c8092SJakub Kruzik Vec r,w1,w2; 164*b27c8092SJakub Kruzik PetscBool nonzero; 165*b27c8092SJakub Kruzik PetscErrorCode ierr; 16637eeb815SJakub Kruzik 167*b27c8092SJakub Kruzik PetscFunctionBegin; 168*b27c8092SJakub Kruzik w1 = def->workcoarse[0]; 169*b27c8092SJakub Kruzik w2 = def->workcoarse[1]; 170*b27c8092SJakub Kruzik r = def->work; 171*b27c8092SJakub Kruzik ierr = PCGetOperators(pc,NULL,&A);CHKERRQ(ierr); 17237eeb815SJakub Kruzik 173*b27c8092SJakub Kruzik ierr = KSPGetInitialGuessNonzero(ksp,&nonzero);CHKERRQ(ierr); 174*b27c8092SJakub Kruzik ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr); 175*b27c8092SJakub Kruzik if (nonzero) { 176*b27c8092SJakub Kruzik ierr = MatMult(A,x,r);CHKERRQ(ierr); /* r <- b - Ax */ 177*b27c8092SJakub Kruzik ierr = VecAYPX(r,-1.0,b);CHKERRQ(ierr); 178*b27c8092SJakub Kruzik } else { 179*b27c8092SJakub Kruzik ierr = VecCopy(b,r);CHKERRQ(ierr); /* r <- b (x is 0) */ 180*b27c8092SJakub Kruzik } 181*b27c8092SJakub Kruzik 182*b27c8092SJakub Kruzik ierr = MatMultTranspose(def->W,r,w1);CHKERRQ(ierr); /* w1 <- W'*r */ 183*b27c8092SJakub Kruzik ierr = KSPSolve(def->WtAWinv,w1,w2);CHKERRQ(ierr); /* w2 <- (W'*A*W)^{-1}*w1 */ 184*b27c8092SJakub Kruzik ierr = MatMult(def->W,w2,r);CHKERRQ(ierr); /* r <- W*w2 */ 185*b27c8092SJakub Kruzik ierr = VecAYPX(x,1.0,r);CHKERRQ(ierr); 186*b27c8092SJakub Kruzik PetscFunctionReturn(0); 187*b27c8092SJakub Kruzik } 18837eeb815SJakub Kruzik 189f8bf229cSJakub Kruzik /* 190f8bf229cSJakub Kruzik if (def->correct) { 191f8bf229cSJakub Kruzik z <- r - W*(W'*A*W)^{-1}*W'*(A*r -r) = (P-Q)*r 192f8bf229cSJakub Kruzik } else { 193f8bf229cSJakub Kruzik z <- r - W*(W'*A*W)^{-1}*W'*A*r = P*r 194f8bf229cSJakub Kruzik } 19537eeb815SJakub Kruzik */ 196f8bf229cSJakub Kruzik static PetscErrorCode PCApply_Deflation(PC pc,Vec r, Vec z) 197f8bf229cSJakub Kruzik { 198f8bf229cSJakub Kruzik PC_Deflation *def = (PC_Deflation*)pc->data; 199f8bf229cSJakub Kruzik Mat A; 200f8bf229cSJakub Kruzik Vec u,w1,w2; 201f8bf229cSJakub Kruzik PetscErrorCode ierr; 202f8bf229cSJakub Kruzik 203f8bf229cSJakub Kruzik PetscFunctionBegin; 204f8bf229cSJakub Kruzik w1 = def->workcoarse[0]; 205f8bf229cSJakub Kruzik w2 = def->workcoarse[1]; 206f8bf229cSJakub Kruzik u = def->work; 207f8bf229cSJakub Kruzik ierr = PCGetOperators(pc,NULL,&A);CHKERRQ(ierr); 208f8bf229cSJakub Kruzik 209f8bf229cSJakub Kruzik if (!def->AW) { 210f8bf229cSJakub Kruzik ierr = MatMult(A,r,u);CHKERRQ(ierr); /* u <- A*r */ 211f8bf229cSJakub Kruzik if (def->correct) ierr = VecAXPY(u,-1.0,r);CHKERRQ(ierr); /* u <- A*r -r */ 212f8bf229cSJakub Kruzik ierr = MatMultTranspose(def->W,u,w1);CHKERRQ(ierr); /* w1 <- W'*u */ 213f8bf229cSJakub Kruzik } else { 214f8bf229cSJakub Kruzik ierr = MatMultTranspose(def->AW,u,w1);CHKERRQ(ierr); /* u <- A*r */ 215f8bf229cSJakub Kruzik if (def->correct) { 216f8bf229cSJakub Kruzik ierr = MatMultTranspose(def->W,r,w2);CHKERRQ(ierr); /* w2 <- W'*u */ 217f8bf229cSJakub Kruzik ierr = VecAXPY(w1,-1.0,w2);CHKERRQ(ierr); /* w1 <- w1 - w2 */ 218f8bf229cSJakub Kruzik } 219f8bf229cSJakub Kruzik } 220f8bf229cSJakub Kruzik ierr = KSPSolve(def->WtAWinv,w1,w2);CHKERRQ(ierr); /* w2 <- (W'*A*W)^{-1}*w1 */ 221f8bf229cSJakub Kruzik ierr = MatMult(def->W,w2,u);CHKERRQ(ierr); /* u <- W*w2 */ 222f8bf229cSJakub Kruzik ierr = VecWAXPY(z,-1.0,u,r);CHKERRQ(ierr); /* z <- r - u */ 223f8bf229cSJakub Kruzik PetscFunctionReturn(0); 224f8bf229cSJakub Kruzik } 225f8bf229cSJakub Kruzik 226*b27c8092SJakub Kruzik static PetscErrorCode PCDeflationSetType_Deflation(PC pc,PCDeflationType type) 227*b27c8092SJakub Kruzik { 228*b27c8092SJakub Kruzik PC_Deflation *def = (PC_Deflation*)pc->data; 229*b27c8092SJakub Kruzik 230*b27c8092SJakub Kruzik PetscFunctionBegin; 231*b27c8092SJakub Kruzik def->init = PETSC_FALSE; 232*b27c8092SJakub Kruzik def->pre = PETSC_FALSE; 233*b27c8092SJakub Kruzik if (type == PC_DEFLATION_POST) { 234*b27c8092SJakub Kruzik //pc->ops->postsolve = PCPostSolve_Deflation; 235*b27c8092SJakub Kruzik pc->ops->presolve = 0; 236*b27c8092SJakub Kruzik } else { 237*b27c8092SJakub Kruzik pc->ops->presolve = PCPreSolve_Deflation; 238*b27c8092SJakub Kruzik pc->ops->postsolve = 0; 239*b27c8092SJakub Kruzik if (type == PC_DEFLATION_INIT) { 240*b27c8092SJakub Kruzik def->init = PETSC_TRUE; 241*b27c8092SJakub Kruzik pc->ops->apply = 0; 242*b27c8092SJakub Kruzik } else { 243*b27c8092SJakub Kruzik def->pre = PETSC_TRUE; 244*b27c8092SJakub Kruzik } 245*b27c8092SJakub Kruzik } 246*b27c8092SJakub Kruzik PetscFunctionReturn(0); 247*b27c8092SJakub Kruzik } 248*b27c8092SJakub Kruzik 249*b27c8092SJakub Kruzik /*@ 250*b27c8092SJakub Kruzik PCDeflationSetType - Causes the deflation preconditioner to use only a special 251*b27c8092SJakub Kruzik initial gues or pre/post solve solution update 252*b27c8092SJakub Kruzik 253*b27c8092SJakub Kruzik Logically Collective on PC 254*b27c8092SJakub Kruzik 255*b27c8092SJakub Kruzik Input Parameters: 256*b27c8092SJakub Kruzik + pc - the preconditioner context 257*b27c8092SJakub Kruzik - type - PC_DEFLATION_PRE, PC_DEFLATION_INIT, PC_DEFLATION_POST 258*b27c8092SJakub Kruzik 259*b27c8092SJakub Kruzik Options Database Key: 260*b27c8092SJakub Kruzik . -pc_deflation_type <pre,init,post> 261*b27c8092SJakub Kruzik 262*b27c8092SJakub Kruzik Level: intermediate 263*b27c8092SJakub Kruzik 264*b27c8092SJakub Kruzik Concepts: Deflation preconditioner 265*b27c8092SJakub Kruzik 266*b27c8092SJakub Kruzik .seealso: PCDeflationGetType() 267*b27c8092SJakub Kruzik @*/ 268*b27c8092SJakub Kruzik PetscErrorCode PCDeflationSetType(PC pc,PCDeflationType type) 269*b27c8092SJakub Kruzik { 270*b27c8092SJakub Kruzik PetscErrorCode ierr; 271*b27c8092SJakub Kruzik 272*b27c8092SJakub Kruzik PetscFunctionBegin; 273*b27c8092SJakub Kruzik PetscValidHeaderSpecific(pc,PC_CLASSID,1); 274*b27c8092SJakub Kruzik ierr = PetscTryMethod(pc,"PCDeflationSetType_C",(PC,PCDeflationType),(pc,type));CHKERRQ(ierr); 275*b27c8092SJakub Kruzik PetscFunctionReturn(0); 276*b27c8092SJakub Kruzik } 277*b27c8092SJakub Kruzik 278*b27c8092SJakub Kruzik static PetscErrorCode PCDeflationGetType_Deflation(PC pc,PCDeflationType *type) 279*b27c8092SJakub Kruzik { 280*b27c8092SJakub Kruzik PC_Deflation *def = (PC_Deflation*)pc->data; 281*b27c8092SJakub Kruzik 282*b27c8092SJakub Kruzik PetscFunctionBegin; 283*b27c8092SJakub Kruzik if (def->init) { 284*b27c8092SJakub Kruzik *type = PC_DEFLATION_INIT; 285*b27c8092SJakub Kruzik } else if (def->pre) { 286*b27c8092SJakub Kruzik *type = PC_DEFLATION_PRE; 287*b27c8092SJakub Kruzik } else { 288*b27c8092SJakub Kruzik *type = PC_DEFLATION_POST; 289*b27c8092SJakub Kruzik } 290*b27c8092SJakub Kruzik PetscFunctionReturn(0); 291*b27c8092SJakub Kruzik } 292*b27c8092SJakub Kruzik 293*b27c8092SJakub Kruzik /*@ 294*b27c8092SJakub Kruzik PCDeflationGetType - Gets how the diagonal matrix is produced for the preconditioner 295*b27c8092SJakub Kruzik 296*b27c8092SJakub Kruzik Not Collective on PC 297*b27c8092SJakub Kruzik 298*b27c8092SJakub Kruzik Input Parameter: 299*b27c8092SJakub Kruzik . pc - the preconditioner context 300*b27c8092SJakub Kruzik 301*b27c8092SJakub Kruzik Output Parameter: 302*b27c8092SJakub Kruzik - type - PC_DEFLATION_PRE, PC_DEFLATION_INIT, PC_DEFLATION_POST 303*b27c8092SJakub Kruzik 304*b27c8092SJakub Kruzik Level: intermediate 305*b27c8092SJakub Kruzik 306*b27c8092SJakub Kruzik Concepts: Deflation preconditioner 307*b27c8092SJakub Kruzik 308*b27c8092SJakub Kruzik .seealso: PCDeflationSetType() 309*b27c8092SJakub Kruzik @*/ 310*b27c8092SJakub Kruzik PetscErrorCode PCDeflationGetType(PC pc,PCDeflationType *type) 311*b27c8092SJakub Kruzik { 312*b27c8092SJakub Kruzik PetscErrorCode ierr; 313*b27c8092SJakub Kruzik 314*b27c8092SJakub Kruzik PetscFunctionBegin; 315*b27c8092SJakub Kruzik PetscValidHeaderSpecific(pc,PC_CLASSID,1); 316*b27c8092SJakub Kruzik ierr = PetscUseMethod(pc,"PCDeflationGetType_C",(PC,PCDeflationType*),(pc,type));CHKERRQ(ierr); 317*b27c8092SJakub Kruzik PetscFunctionReturn(0); 318*b27c8092SJakub Kruzik } 319*b27c8092SJakub Kruzik 32037eeb815SJakub Kruzik static PetscErrorCode PCSetUp_Deflation(PC pc) 32137eeb815SJakub Kruzik { 322409a357bSJakub Kruzik PC_Deflation *def = (PC_Deflation*)pc->data; 323409a357bSJakub Kruzik KSP innerksp; 324409a357bSJakub Kruzik PC pcinner; 325409a357bSJakub Kruzik Mat Amat,nextDef=NULL,*mats; 326409a357bSJakub Kruzik PetscInt i,m,red,size,commsize; 327409a357bSJakub Kruzik PetscBool match,flgspd,transp=PETSC_FALSE; 328409a357bSJakub Kruzik MatCompositeType ctype; 329409a357bSJakub Kruzik MPI_Comm comm; 330409a357bSJakub Kruzik const char *prefix; 33137eeb815SJakub Kruzik PetscErrorCode ierr; 33237eeb815SJakub Kruzik 33337eeb815SJakub Kruzik PetscFunctionBegin; 334409a357bSJakub Kruzik if (pc->setupcalled) PetscFunctionReturn(0); 335409a357bSJakub Kruzik ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 336f8bf229cSJakub Kruzik ierr = PCGetOperators(pc,NULL,&Amat);CHKERRQ(ierr); 337f8bf229cSJakub Kruzik 338f8bf229cSJakub Kruzik /* compute a deflation space */ 339409a357bSJakub Kruzik if (def->W || def->Wt) { 340409a357bSJakub Kruzik def->spacetype = PC_DEFLATION_SPACE_USER; 341409a357bSJakub Kruzik } else { 342409a357bSJakub Kruzik //ierr = KSPDCGComputeDeflationSpace(ksp);CHKERRQ(ierr); 34337eeb815SJakub Kruzik } 34437eeb815SJakub Kruzik 345409a357bSJakub Kruzik /* nested deflation */ 346409a357bSJakub Kruzik if (def->W) { 347409a357bSJakub Kruzik ierr = PetscObjectTypeCompare((PetscObject)def->W,MATCOMPOSITE,&match);CHKERRQ(ierr); 348409a357bSJakub Kruzik if (match) { 349409a357bSJakub Kruzik ierr = MatCompositeGetType(def->W,&ctype);CHKERRQ(ierr); 350409a357bSJakub Kruzik ierr = MatCompositeGetNumberMat(def->W,&size);CHKERRQ(ierr); 35137eeb815SJakub Kruzik } 352409a357bSJakub Kruzik } else { 353409a357bSJakub Kruzik ierr = MatCreateTranspose(def->Wt,&def->W);CHKERRQ(ierr); 354409a357bSJakub Kruzik ierr = PetscObjectTypeCompare((PetscObject)def->Wt,MATCOMPOSITE,&match);CHKERRQ(ierr); 355409a357bSJakub Kruzik if (match) { 356409a357bSJakub Kruzik ierr = MatCompositeGetType(def->Wt,&ctype);CHKERRQ(ierr); 357409a357bSJakub Kruzik ierr = MatCompositeGetNumberMat(def->Wt,&size);CHKERRQ(ierr); 358409a357bSJakub Kruzik } 359409a357bSJakub Kruzik transp = PETSC_TRUE; 360409a357bSJakub Kruzik } 361409a357bSJakub Kruzik if (match && ctype == MAT_COMPOSITE_MULTIPLICATIVE) { 362409a357bSJakub Kruzik ierr = PetscMalloc1(size,&mats);CHKERRQ(ierr); 363409a357bSJakub Kruzik if (!transp) { 364409a357bSJakub Kruzik for (i=0; i<size; i++) { 365409a357bSJakub Kruzik ierr = MatCompositeGetMat(def->W,i,&mats[i]);CHKERRQ(ierr); 366409a357bSJakub Kruzik //ierr = PetscObjectReference((PetscObject)mats[i]);CHKERRQ(ierr); 367409a357bSJakub Kruzik } 368409a357bSJakub Kruzik if (def->nestedlvl < def->maxnestedlvl) { 369409a357bSJakub Kruzik size -= 1; 370409a357bSJakub Kruzik ierr = MatDestroy(&def->W);CHKERRQ(ierr); 371409a357bSJakub Kruzik def->W = mats[size]; 372409a357bSJakub Kruzik ierr = PetscObjectReference((PetscObject)mats[size]);CHKERRQ(ierr); 373409a357bSJakub Kruzik if (size > 1) { 374409a357bSJakub Kruzik ierr = MatCreateComposite(comm,size,mats,&nextDef);CHKERRQ(ierr); 375409a357bSJakub Kruzik ierr = MatCompositeSetType(nextDef,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr); 376409a357bSJakub Kruzik } else { 377409a357bSJakub Kruzik nextDef = mats[0]; 378409a357bSJakub Kruzik ierr = PetscObjectReference((PetscObject)mats[0]);CHKERRQ(ierr); 379409a357bSJakub Kruzik } 380409a357bSJakub Kruzik } else { 381409a357bSJakub Kruzik /* ierr = MatCompositeSetMergeType(def->W,MAT_COMPOSITE_MERGE_LEFT);CHKERRQ(ierr); */ 382409a357bSJakub Kruzik ierr = MatCompositeMerge(def->W);CHKERRQ(ierr); 383409a357bSJakub Kruzik } 384409a357bSJakub Kruzik } 385409a357bSJakub Kruzik ierr = PetscFree(mats);CHKERRQ(ierr); 386409a357bSJakub Kruzik } 387409a357bSJakub Kruzik 388409a357bSJakub Kruzik /* setup coarse problem */ 389409a357bSJakub Kruzik if (!def->WtAWinv) { 390409a357bSJakub Kruzik ierr = MatGetSize(def->W,NULL,&m);CHKERRQ(ierr); /* TODO works for W MatTranspose? */ 391409a357bSJakub Kruzik if (!def->WtAW) { 392409a357bSJakub Kruzik /* TODO add implicit product version ? */ 393409a357bSJakub Kruzik if (!def->AW) { 394409a357bSJakub Kruzik ierr = MatPtAP(Amat,def->W,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtAW);CHKERRQ(ierr); 395409a357bSJakub Kruzik } else { 396409a357bSJakub Kruzik ierr = MatTransposeMatMult(def->W,def->AW,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtAW);CHKERRQ(ierr); 397409a357bSJakub Kruzik } 398409a357bSJakub Kruzik /* TODO create MatInheritOption(Mat,MatOption) */ 399409a357bSJakub Kruzik ierr = MatGetOption(Amat,MAT_SPD,&flgspd);CHKERRQ(ierr); 400409a357bSJakub Kruzik ierr = MatSetOption(def->WtAW,MAT_SPD,flgspd);CHKERRQ(ierr); 401409a357bSJakub Kruzik #if defined(PETSC_USE_DEBUG) 402409a357bSJakub Kruzik /* Check WtAW is not sigular */ 403409a357bSJakub Kruzik PetscReal *norms; 404409a357bSJakub Kruzik ierr = PetscMalloc1(m,&norms);CHKERRQ(ierr); 405409a357bSJakub Kruzik ierr = MatGetColumnNorms(def->WtAW,NORM_INFINITY,norms);CHKERRQ(ierr); 406409a357bSJakub Kruzik for (i=0; i<m; i++) { 407409a357bSJakub Kruzik if (norms[i] < 100*PETSC_MACHINE_EPSILON) { 408409a357bSJakub Kruzik SETERRQ1(comm,PETSC_ERR_SUP,"Column %D of W is in kernel of A.",i); 409409a357bSJakub Kruzik } 410409a357bSJakub Kruzik } 411409a357bSJakub Kruzik ierr = PetscFree(norms);CHKERRQ(ierr); 412409a357bSJakub Kruzik #endif 413409a357bSJakub Kruzik } else { 414409a357bSJakub Kruzik ierr = MatGetOption(def->WtAW,MAT_SPD,&flgspd);CHKERRQ(ierr); 415409a357bSJakub Kruzik } 416409a357bSJakub Kruzik /* TODO use MATINV */ 417409a357bSJakub Kruzik ierr = KSPCreate(comm,&def->WtAWinv);CHKERRQ(ierr); 418409a357bSJakub Kruzik ierr = KSPSetOperators(def->WtAWinv,def->WtAW,def->WtAW);CHKERRQ(ierr); 419409a357bSJakub Kruzik ierr = KSPSetType(def->WtAWinv,KSPPREONLY);CHKERRQ(ierr); 420409a357bSJakub Kruzik ierr = KSPGetPC(def->WtAWinv,&pcinner);CHKERRQ(ierr); 421409a357bSJakub Kruzik ierr = PCSetType(pcinner,PCTELESCOPE);CHKERRQ(ierr); 422409a357bSJakub Kruzik /* Reduction factor choice */ 423409a357bSJakub Kruzik red = def->reductionfact; 424409a357bSJakub Kruzik if (red < 0) { 425409a357bSJakub Kruzik ierr = MPI_Comm_size(comm,&commsize);CHKERRQ(ierr); 426409a357bSJakub Kruzik red = ceil((float)commsize/ceil((float)m/commsize)); 427409a357bSJakub Kruzik ierr = PetscObjectTypeCompareAny((PetscObject)(def->WtAW),&match,MATSEQDENSE,MATMPIDENSE,MATDENSE,"");CHKERRQ(ierr); 428409a357bSJakub Kruzik if (match) red = commsize; 429409a357bSJakub Kruzik ierr = PetscInfo1(pc,"Auto choosing reduction factor %D\n",red);CHKERRQ(ierr); /* TODO add level? */ 430409a357bSJakub Kruzik } 431409a357bSJakub Kruzik ierr = PCTelescopeSetReductionFactor(pcinner,red);CHKERRQ(ierr); 432409a357bSJakub Kruzik ierr = PCTelescopeGetKSP(pcinner,&innerksp);CHKERRQ(ierr); 433409a357bSJakub Kruzik ierr = KSPGetPC(innerksp,&pcinner);CHKERRQ(ierr); 434409a357bSJakub Kruzik /* Setup KSP and PC */ 435409a357bSJakub Kruzik if (nextDef) { /* next level for multilevel deflation */ 436409a357bSJakub Kruzik /* set default KSPtype */ 437409a357bSJakub Kruzik if (!def->ksptype) { 438409a357bSJakub Kruzik def->ksptype = KSPFGMRES; 439409a357bSJakub Kruzik if (flgspd) { /* SPD system */ 440409a357bSJakub Kruzik def->ksptype = KSPFCG; 441409a357bSJakub Kruzik } 442409a357bSJakub Kruzik } 443409a357bSJakub Kruzik ierr = KSPSetType(innerksp,def->ksptype);CHKERRQ(ierr); /* TODO iherit from KSP */ 444409a357bSJakub Kruzik ierr = PCSetType(pcinner,PCDEFLATION);CHKERRQ(ierr); /* TODO create coarse preconditinoner M_c = WtMW ? */ 445409a357bSJakub Kruzik ierr = PCDeflationSetSpace(pcinner,nextDef,transp);CHKERRQ(ierr); 4465c0c31c5SJakub Kruzik ierr = PCDeflationSetLvl_Deflation(pcinner,def->nestedlvl+1,def->maxnestedlvl);CHKERRQ(ierr); 447409a357bSJakub Kruzik /* inherit options TODO if not set */ 448409a357bSJakub Kruzik ((PC_Deflation*)(pcinner))->ksptype = def->ksptype; 449409a357bSJakub Kruzik ((PC_Deflation*)(pcinner))->correct = def->correct; 450409a357bSJakub Kruzik ((PC_Deflation*)(pcinner))->adaptiveconv = def->adaptiveconv; 451409a357bSJakub Kruzik ((PC_Deflation*)(pcinner))->adaptiveconst = def->adaptiveconst; 452409a357bSJakub Kruzik ierr = MatDestroy(&nextDef);CHKERRQ(ierr); 453409a357bSJakub Kruzik } else { /* the last level */ 454409a357bSJakub Kruzik ierr = KSPSetType(innerksp,KSPPREONLY);CHKERRQ(ierr); 455409a357bSJakub Kruzik /* TODO Cholesky if flgspd? */ 456409a357bSJakub Kruzik ierr = PCSetType(pc,PCLU);CHKERRQ(ierr); 457409a357bSJakub Kruzik //TODO remove explicit matSolverPackage 458409a357bSJakub Kruzik if (commsize == red) { 459409a357bSJakub Kruzik ierr = PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU);CHKERRQ(ierr); 460409a357bSJakub Kruzik } else { 461409a357bSJakub Kruzik ierr = PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU_DIST);CHKERRQ(ierr); 462409a357bSJakub Kruzik } 463409a357bSJakub Kruzik } 464409a357bSJakub Kruzik /* TODO use def_[lvl]_ if lvl > 0? */ 465409a357bSJakub Kruzik ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr); 466409a357bSJakub Kruzik if (prefix) { 467409a357bSJakub Kruzik ierr = KSPSetOptionsPrefix(innerksp,prefix);CHKERRQ(ierr); 468409a357bSJakub Kruzik ierr = KSPAppendOptionsPrefix(innerksp,"def_");CHKERRQ(ierr); 469409a357bSJakub Kruzik } else { 470409a357bSJakub Kruzik ierr = KSPSetOptionsPrefix(innerksp,"def_");CHKERRQ(ierr); 471409a357bSJakub Kruzik } 472409a357bSJakub Kruzik /* TODO: check if WtAWinv is KSP and move following from this if */ 473409a357bSJakub Kruzik ierr = KSPSetFromOptions(def->WtAWinv);CHKERRQ(ierr); 474409a357bSJakub Kruzik //if (def->adaptiveconv) { 475409a357bSJakub Kruzik // PetscReal *rnorm; 476409a357bSJakub Kruzik // PetscNew(&rnorm); 477409a357bSJakub Kruzik // ierr = KSPSetConvergenceTest(def->WtAWinv,KSPDCGConvergedAdaptive_DCG,rnorm,NULL);CHKERRQ(ierr); 478409a357bSJakub Kruzik //} 479409a357bSJakub Kruzik ierr = KSPSetUp(def->WtAWinv);CHKERRQ(ierr); 480409a357bSJakub Kruzik } 481409a357bSJakub Kruzik 482f8bf229cSJakub Kruzik /* create work vecs */ 483f8bf229cSJakub Kruzik ierr = MatCreateVecs(Amat,NULL,&def->work);CHKERRQ(ierr); 484f8bf229cSJakub Kruzik ierr = KSPCreateVecs(def->WtAWinv,2,&def->workcoarse,0,NULL);CHKERRQ(ierr); 48537eeb815SJakub Kruzik PetscFunctionReturn(0); 48637eeb815SJakub Kruzik } 48737eeb815SJakub Kruzik /* -------------------------------------------------------------------------- */ 48837eeb815SJakub Kruzik static PetscErrorCode PCReset_Deflation(PC pc) 48937eeb815SJakub Kruzik { 49037eeb815SJakub Kruzik PC_Deflation *jac = (PC_Deflation*)pc->data; 49137eeb815SJakub Kruzik PetscErrorCode ierr; 49237eeb815SJakub Kruzik 49337eeb815SJakub Kruzik PetscFunctionBegin; 49437eeb815SJakub Kruzik PetscFunctionReturn(0); 49537eeb815SJakub Kruzik } 49637eeb815SJakub Kruzik 49737eeb815SJakub Kruzik /* 49837eeb815SJakub Kruzik PCDestroy_Deflation - Destroys the private context for the Deflation preconditioner 49937eeb815SJakub Kruzik that was created with PCCreate_Deflation(). 50037eeb815SJakub Kruzik 50137eeb815SJakub Kruzik Input Parameter: 50237eeb815SJakub Kruzik . pc - the preconditioner context 50337eeb815SJakub Kruzik 50437eeb815SJakub Kruzik Application Interface Routine: PCDestroy() 50537eeb815SJakub Kruzik */ 50637eeb815SJakub Kruzik static PetscErrorCode PCDestroy_Deflation(PC pc) 50737eeb815SJakub Kruzik { 50837eeb815SJakub Kruzik PetscErrorCode ierr; 50937eeb815SJakub Kruzik 51037eeb815SJakub Kruzik PetscFunctionBegin; 51137eeb815SJakub Kruzik ierr = PCReset_Deflation(pc);CHKERRQ(ierr); 51237eeb815SJakub Kruzik 51337eeb815SJakub Kruzik /* 51437eeb815SJakub Kruzik Free the private data structure that was hanging off the PC 51537eeb815SJakub Kruzik */ 51637eeb815SJakub Kruzik ierr = PetscFree(pc->data);CHKERRQ(ierr); 51737eeb815SJakub Kruzik PetscFunctionReturn(0); 51837eeb815SJakub Kruzik } 51937eeb815SJakub Kruzik 52037eeb815SJakub Kruzik static PetscErrorCode PCSetFromOptions_Deflation(PetscOptionItems *PetscOptionsObject,PC pc) 52137eeb815SJakub Kruzik { 52237eeb815SJakub Kruzik PC_Deflation *jac = (PC_Deflation*)pc->data; 52337eeb815SJakub Kruzik PetscErrorCode ierr; 52437eeb815SJakub Kruzik PetscBool flg; 52537eeb815SJakub Kruzik PCDeflationType deflt,type; 52637eeb815SJakub Kruzik 52737eeb815SJakub Kruzik PetscFunctionBegin; 52837eeb815SJakub Kruzik ierr = PCDeflationGetType(pc,&deflt);CHKERRQ(ierr); 52937eeb815SJakub Kruzik ierr = PetscOptionsHead(PetscOptionsObject,"Deflation options");CHKERRQ(ierr); 53037eeb815SJakub Kruzik ierr = PetscOptionsEnum("-pc_jacobi_type","How to construct diagonal matrix","PCDeflationSetType",PCDeflationTypes,(PetscEnum)deflt,(PetscEnum*)&type,&flg);CHKERRQ(ierr); 53137eeb815SJakub Kruzik if (flg) { 53237eeb815SJakub Kruzik ierr = PCDeflationSetType(pc,type);CHKERRQ(ierr); 53337eeb815SJakub Kruzik } 53437eeb815SJakub Kruzik ierr = PetscOptionsTail();CHKERRQ(ierr); 53537eeb815SJakub Kruzik PetscFunctionReturn(0); 53637eeb815SJakub Kruzik } 53737eeb815SJakub Kruzik 53837eeb815SJakub Kruzik /*MC 539e361f8b5SJakub Kruzik PCDEFLATION - Deflation preconditioner shifts part of the spectrum to zero (deflates) 540e361f8b5SJakub Kruzik or to a predefined value 54137eeb815SJakub Kruzik 54237eeb815SJakub Kruzik Options Database Key: 543e361f8b5SJakub Kruzik + -pc_deflation_type <init,pre,post> - selects approach to deflation (default: pre) 54437eeb815SJakub Kruzik - -pc_jacobi_abs - use the absolute value of the diagonal entry 54537eeb815SJakub Kruzik 54637eeb815SJakub Kruzik Level: beginner 54737eeb815SJakub Kruzik 54837eeb815SJakub Kruzik Notes: 549e361f8b5SJakub Kruzik todo 55037eeb815SJakub Kruzik 55137eeb815SJakub Kruzik .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 552e662bc50SJakub Kruzik PCDeflationSetType(), PCDeflationSetSpace() 55337eeb815SJakub Kruzik M*/ 55437eeb815SJakub Kruzik 55537eeb815SJakub Kruzik PETSC_EXTERN PetscErrorCode PCCreate_Deflation(PC pc) 55637eeb815SJakub Kruzik { 55737eeb815SJakub Kruzik PC_Deflation *def; 55837eeb815SJakub Kruzik PetscErrorCode ierr; 55937eeb815SJakub Kruzik 56037eeb815SJakub Kruzik PetscFunctionBegin; 56137eeb815SJakub Kruzik ierr = PetscNewLog(pc,&def);CHKERRQ(ierr); 56237eeb815SJakub Kruzik pc->data = (void*)def; 56337eeb815SJakub Kruzik 564e662bc50SJakub Kruzik def->init = PETSC_FALSE; 565e662bc50SJakub Kruzik def->pre = PETSC_TRUE; 566e662bc50SJakub Kruzik def->correct = PETSC_FALSE; 567e662bc50SJakub Kruzik def->truenorm = PETSC_TRUE; 568e662bc50SJakub Kruzik def->reductionfact = -1; 569e662bc50SJakub Kruzik def->spacetype = PC_DEFLATION_SPACE_HAAR; 570e662bc50SJakub Kruzik def->spacesize = 1; 571e662bc50SJakub Kruzik def->extendsp = PETSC_FALSE; 572e662bc50SJakub Kruzik def->nestedlvl = 0; 573e662bc50SJakub Kruzik def->maxnestedlvl = 0; 574e662bc50SJakub Kruzik def->adaptiveconv = PETSC_FALSE; 575e662bc50SJakub Kruzik def->adaptiveconst = 1.0; 57637eeb815SJakub Kruzik 57737eeb815SJakub Kruzik /* 57837eeb815SJakub Kruzik Set the pointers for the functions that are provided above. 57937eeb815SJakub Kruzik Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 58037eeb815SJakub Kruzik are called, they will automatically call these functions. Note we 58137eeb815SJakub Kruzik choose not to provide a couple of these functions since they are 58237eeb815SJakub Kruzik not needed. 58337eeb815SJakub Kruzik */ 58437eeb815SJakub Kruzik pc->ops->apply = PCApply_Deflation; 58537eeb815SJakub Kruzik pc->ops->applytranspose = PCApply_Deflation; 586*b27c8092SJakub Kruzik pc->ops->presolve = PCPreSolve_Deflation; 587*b27c8092SJakub Kruzik pc->ops->postsolve = 0; 58837eeb815SJakub Kruzik pc->ops->setup = PCSetUp_Deflation; 58937eeb815SJakub Kruzik pc->ops->reset = PCReset_Deflation; 59037eeb815SJakub Kruzik pc->ops->destroy = PCDestroy_Deflation; 59137eeb815SJakub Kruzik pc->ops->setfromoptions = PCSetFromOptions_Deflation; 59237eeb815SJakub Kruzik pc->ops->view = 0; 59337eeb815SJakub Kruzik pc->ops->applyrichardson = 0; 59437eeb815SJakub Kruzik pc->ops->applysymmetricleft = 0; 59537eeb815SJakub Kruzik pc->ops->applysymmetricright = 0; 59637eeb815SJakub Kruzik 59737eeb815SJakub Kruzik ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetType_C",PCDeflationSetType_Deflation);CHKERRQ(ierr); 59837eeb815SJakub Kruzik ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationGetType_C",PCDeflationGetType_Deflation);CHKERRQ(ierr); 599e662bc50SJakub Kruzik ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetSpace_C",PCDeflationSetSpace_Deflation);CHKERRQ(ierr); 6005c0c31c5SJakub Kruzik ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetLvl_C",PCDeflationSetLvl_Deflation);CHKERRQ(ierr); 60137eeb815SJakub Kruzik PetscFunctionReturn(0); 60237eeb815SJakub Kruzik } 60337eeb815SJakub Kruzik 604