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 51*292e2e67SJakub Kruzik #include <../src/ksp/pc/impls/deflation/deflation.h> /*I "petscpc.h" I*/ /* includes for fortran wrappers */ 5237eeb815SJakub Kruzik 5337eeb815SJakub Kruzik const char *const PCDeflationTypes[] = {"INIT","PRE","POST","PCDeflationType","PC_DEFLATION_",0}; 5437eeb815SJakub Kruzik 55e662bc50SJakub Kruzik static PetscErrorCode PCDeflationSetSpace_Deflation(PC pc,Mat W,PetscBool transpose) 56e662bc50SJakub Kruzik { 57e662bc50SJakub Kruzik PC_Deflation *def = (PC_Deflation*)pc->data; 58e662bc50SJakub Kruzik PetscErrorCode ierr; 59e662bc50SJakub Kruzik 60e662bc50SJakub Kruzik PetscFunctionBegin; 61e662bc50SJakub Kruzik if (transpose) { 62e662bc50SJakub Kruzik def->Wt = W; 63e662bc50SJakub Kruzik def->W = NULL; 64e662bc50SJakub Kruzik } else { 65e662bc50SJakub Kruzik def->W = W; 66e662bc50SJakub Kruzik } 67e662bc50SJakub Kruzik ierr = PetscObjectReference((PetscObject)W);CHKERRQ(ierr); 68e662bc50SJakub Kruzik PetscFunctionReturn(0); 69e662bc50SJakub Kruzik } 70e662bc50SJakub Kruzik 71e662bc50SJakub Kruzik /* TODO create PCDeflationSetSpaceTranspose? */ 72e662bc50SJakub Kruzik /*@ 73e662bc50SJakub Kruzik PCDeflationSetSpace - Set deflation space matrix (or its transpose). 74e662bc50SJakub Kruzik 75e662bc50SJakub Kruzik Logically Collective on PC 76e662bc50SJakub Kruzik 77e662bc50SJakub Kruzik Input Parameters: 78e662bc50SJakub Kruzik + pc - the preconditioner context 79e662bc50SJakub Kruzik . W - deflation matrix 80e662bc50SJakub Kruzik - tranpose - indicates that W is an explicit transpose of the deflation matrix 81e662bc50SJakub Kruzik 82e662bc50SJakub Kruzik Level: intermediate 83e662bc50SJakub Kruzik 84e662bc50SJakub Kruzik .seealso: PCDEFLATION 85e662bc50SJakub Kruzik @*/ 86e662bc50SJakub Kruzik PetscErrorCode PCDeflationSetSpace(PC pc,Mat W,PetscBool transpose) 87e662bc50SJakub Kruzik { 88e662bc50SJakub Kruzik PetscErrorCode ierr; 89e662bc50SJakub Kruzik 90e662bc50SJakub Kruzik PetscFunctionBegin; 91e662bc50SJakub Kruzik PetscValidHeaderSpecific(pc,PC_CLASSID,1); 92e662bc50SJakub Kruzik PetscValidHeaderSpecific(W,MAT_CLASSID,2); 93e662bc50SJakub Kruzik PetscValidLogicalCollectiveBool(pc,transpose,3); 94e662bc50SJakub Kruzik ierr = PetscTryMethod(pc,"PCDeflationSetSpace_C",(PC,Mat,PetscBool),(pc,W,transpose));CHKERRQ(ierr); 95e662bc50SJakub Kruzik PetscFunctionReturn(0); 96e662bc50SJakub Kruzik } 97e662bc50SJakub Kruzik 985c0c31c5SJakub Kruzik static PetscErrorCode PCDeflationSetLvl_Deflation(PC pc,PetscInt current,PetscInt max) 995c0c31c5SJakub Kruzik { 1005c0c31c5SJakub Kruzik PC_Deflation *def = (PC_Deflation*)pc->data; 1015c0c31c5SJakub Kruzik 1025c0c31c5SJakub Kruzik PetscFunctionBegin; 1035c0c31c5SJakub Kruzik def->nestedlvl = current; 1045c0c31c5SJakub Kruzik def->maxnestedlvl = max; 1055c0c31c5SJakub Kruzik PetscFunctionReturn(0); 1065c0c31c5SJakub Kruzik } 1075c0c31c5SJakub Kruzik 1085c0c31c5SJakub Kruzik /*@ 1095c0c31c5SJakub Kruzik PCDeflationSetMaxLvl - Set maximum level of deflation. 1105c0c31c5SJakub Kruzik 1115c0c31c5SJakub Kruzik Logically Collective on PC 1125c0c31c5SJakub Kruzik 1135c0c31c5SJakub Kruzik Input Parameters: 1145c0c31c5SJakub Kruzik + pc - the preconditioner context 1155c0c31c5SJakub Kruzik . max - maximum deflation level 1165c0c31c5SJakub Kruzik 1175c0c31c5SJakub Kruzik Level: intermediate 1185c0c31c5SJakub Kruzik 1195c0c31c5SJakub Kruzik .seealso: PCDEFLATION 1205c0c31c5SJakub Kruzik @*/ 1215c0c31c5SJakub Kruzik PetscErrorCode PCDeflationSetMaxLvl(PC pc,PetscInt max) 1225c0c31c5SJakub Kruzik { 1235c0c31c5SJakub Kruzik PetscErrorCode ierr; 1245c0c31c5SJakub Kruzik 1255c0c31c5SJakub Kruzik PetscFunctionBegin; 1265c0c31c5SJakub Kruzik PetscValidLogicalCollectiveInt(pc,max,2); 1275c0c31c5SJakub Kruzik ierr = PetscTryMethod(pc,"PCDeflationSetLvl_C",(PC,PetscInt,PetscInt),(pc,0,max));CHKERRQ(ierr); 1285c0c31c5SJakub Kruzik PetscFunctionReturn(0); 1295c0c31c5SJakub Kruzik } 130e662bc50SJakub Kruzik 13137eeb815SJakub Kruzik /* 132b27c8092SJakub Kruzik TODO CP corection? 133b27c8092SJakub Kruzik x <- x + W*(W'*A*W)^{-1}*W'*r = x + Q*r 134b27c8092SJakub Kruzik */ 135b27c8092SJakub Kruzik static PetscErrorCode PCPreSolve_Deflation(PC pc,KSP ksp,Vec b, Vec x) 136b27c8092SJakub Kruzik { 137b27c8092SJakub Kruzik PC_Deflation *def = (PC_Deflation*)pc->data; 138b27c8092SJakub Kruzik Mat A; 139b27c8092SJakub Kruzik Vec r,w1,w2; 140b27c8092SJakub Kruzik PetscBool nonzero; 141b27c8092SJakub Kruzik PetscErrorCode ierr; 14237eeb815SJakub Kruzik 143b27c8092SJakub Kruzik PetscFunctionBegin; 144b27c8092SJakub Kruzik w1 = def->workcoarse[0]; 145b27c8092SJakub Kruzik w2 = def->workcoarse[1]; 146b27c8092SJakub Kruzik r = def->work; 147b27c8092SJakub Kruzik ierr = PCGetOperators(pc,NULL,&A);CHKERRQ(ierr); 14837eeb815SJakub Kruzik 149b27c8092SJakub Kruzik ierr = KSPGetInitialGuessNonzero(ksp,&nonzero);CHKERRQ(ierr); 150b27c8092SJakub Kruzik ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr); 151b27c8092SJakub Kruzik if (nonzero) { 152b27c8092SJakub Kruzik ierr = MatMult(A,x,r);CHKERRQ(ierr); /* r <- b - Ax */ 153b27c8092SJakub Kruzik ierr = VecAYPX(r,-1.0,b);CHKERRQ(ierr); 154b27c8092SJakub Kruzik } else { 155b27c8092SJakub Kruzik ierr = VecCopy(b,r);CHKERRQ(ierr); /* r <- b (x is 0) */ 156b27c8092SJakub Kruzik } 157b27c8092SJakub Kruzik 158b27c8092SJakub Kruzik ierr = MatMultTranspose(def->W,r,w1);CHKERRQ(ierr); /* w1 <- W'*r */ 159b27c8092SJakub Kruzik ierr = KSPSolve(def->WtAWinv,w1,w2);CHKERRQ(ierr); /* w2 <- (W'*A*W)^{-1}*w1 */ 160b27c8092SJakub Kruzik ierr = MatMult(def->W,w2,r);CHKERRQ(ierr); /* r <- W*w2 */ 161b27c8092SJakub Kruzik ierr = VecAYPX(x,1.0,r);CHKERRQ(ierr); 162b27c8092SJakub Kruzik PetscFunctionReturn(0); 163b27c8092SJakub Kruzik } 16437eeb815SJakub Kruzik 165f8bf229cSJakub Kruzik /* 166f8bf229cSJakub Kruzik if (def->correct) { 167f8bf229cSJakub Kruzik z <- r - W*(W'*A*W)^{-1}*W'*(A*r -r) = (P-Q)*r 168f8bf229cSJakub Kruzik } else { 169f8bf229cSJakub Kruzik z <- r - W*(W'*A*W)^{-1}*W'*A*r = P*r 170f8bf229cSJakub Kruzik } 17137eeb815SJakub Kruzik */ 172f8bf229cSJakub Kruzik static PetscErrorCode PCApply_Deflation(PC pc,Vec r, Vec z) 173f8bf229cSJakub Kruzik { 174f8bf229cSJakub Kruzik PC_Deflation *def = (PC_Deflation*)pc->data; 175f8bf229cSJakub Kruzik Mat A; 176f8bf229cSJakub Kruzik Vec u,w1,w2; 177f8bf229cSJakub Kruzik PetscErrorCode ierr; 178f8bf229cSJakub Kruzik 179f8bf229cSJakub Kruzik PetscFunctionBegin; 180f8bf229cSJakub Kruzik w1 = def->workcoarse[0]; 181f8bf229cSJakub Kruzik w2 = def->workcoarse[1]; 182f8bf229cSJakub Kruzik u = def->work; 183f8bf229cSJakub Kruzik ierr = PCGetOperators(pc,NULL,&A);CHKERRQ(ierr); 184f8bf229cSJakub Kruzik 185f8bf229cSJakub Kruzik if (!def->AW) { 186f8bf229cSJakub Kruzik ierr = MatMult(A,r,u);CHKERRQ(ierr); /* u <- A*r */ 187f8bf229cSJakub Kruzik if (def->correct) ierr = VecAXPY(u,-1.0,r);CHKERRQ(ierr); /* u <- A*r -r */ 188f8bf229cSJakub Kruzik ierr = MatMultTranspose(def->W,u,w1);CHKERRQ(ierr); /* w1 <- W'*u */ 189f8bf229cSJakub Kruzik } else { 190f8bf229cSJakub Kruzik ierr = MatMultTranspose(def->AW,u,w1);CHKERRQ(ierr); /* u <- A*r */ 191f8bf229cSJakub Kruzik if (def->correct) { 192f8bf229cSJakub Kruzik ierr = MatMultTranspose(def->W,r,w2);CHKERRQ(ierr); /* w2 <- W'*u */ 193f8bf229cSJakub Kruzik ierr = VecAXPY(w1,-1.0,w2);CHKERRQ(ierr); /* w1 <- w1 - w2 */ 194f8bf229cSJakub Kruzik } 195f8bf229cSJakub Kruzik } 196f8bf229cSJakub Kruzik ierr = KSPSolve(def->WtAWinv,w1,w2);CHKERRQ(ierr); /* w2 <- (W'*A*W)^{-1}*w1 */ 197f8bf229cSJakub Kruzik ierr = MatMult(def->W,w2,u);CHKERRQ(ierr); /* u <- W*w2 */ 198f8bf229cSJakub Kruzik ierr = VecWAXPY(z,-1.0,u,r);CHKERRQ(ierr); /* z <- r - u */ 199f8bf229cSJakub Kruzik PetscFunctionReturn(0); 200f8bf229cSJakub Kruzik } 201f8bf229cSJakub Kruzik 202b27c8092SJakub Kruzik static PetscErrorCode PCDeflationSetType_Deflation(PC pc,PCDeflationType type) 203b27c8092SJakub Kruzik { 204b27c8092SJakub Kruzik PC_Deflation *def = (PC_Deflation*)pc->data; 205b27c8092SJakub Kruzik 206b27c8092SJakub Kruzik PetscFunctionBegin; 207b27c8092SJakub Kruzik def->init = PETSC_FALSE; 208b27c8092SJakub Kruzik def->pre = PETSC_FALSE; 209b27c8092SJakub Kruzik if (type == PC_DEFLATION_POST) { 210b27c8092SJakub Kruzik //pc->ops->postsolve = PCPostSolve_Deflation; 211b27c8092SJakub Kruzik pc->ops->presolve = 0; 212b27c8092SJakub Kruzik } else { 213b27c8092SJakub Kruzik pc->ops->presolve = PCPreSolve_Deflation; 214b27c8092SJakub Kruzik pc->ops->postsolve = 0; 215b27c8092SJakub Kruzik if (type == PC_DEFLATION_INIT) { 216b27c8092SJakub Kruzik def->init = PETSC_TRUE; 217b27c8092SJakub Kruzik pc->ops->apply = 0; 218b27c8092SJakub Kruzik } else { 219b27c8092SJakub Kruzik def->pre = PETSC_TRUE; 220b27c8092SJakub Kruzik } 221b27c8092SJakub Kruzik } 222b27c8092SJakub Kruzik PetscFunctionReturn(0); 223b27c8092SJakub Kruzik } 224b27c8092SJakub Kruzik 225b27c8092SJakub Kruzik /*@ 226b27c8092SJakub Kruzik PCDeflationSetType - Causes the deflation preconditioner to use only a special 227b27c8092SJakub Kruzik initial gues or pre/post solve solution update 228b27c8092SJakub Kruzik 229b27c8092SJakub Kruzik Logically Collective on PC 230b27c8092SJakub Kruzik 231b27c8092SJakub Kruzik Input Parameters: 232b27c8092SJakub Kruzik + pc - the preconditioner context 233b27c8092SJakub Kruzik - type - PC_DEFLATION_PRE, PC_DEFLATION_INIT, PC_DEFLATION_POST 234b27c8092SJakub Kruzik 235b27c8092SJakub Kruzik Options Database Key: 236b27c8092SJakub Kruzik . -pc_deflation_type <pre,init,post> 237b27c8092SJakub Kruzik 238b27c8092SJakub Kruzik Level: intermediate 239b27c8092SJakub Kruzik 240b27c8092SJakub Kruzik Concepts: Deflation preconditioner 241b27c8092SJakub Kruzik 242b27c8092SJakub Kruzik .seealso: PCDeflationGetType() 243b27c8092SJakub Kruzik @*/ 244b27c8092SJakub Kruzik PetscErrorCode PCDeflationSetType(PC pc,PCDeflationType type) 245b27c8092SJakub Kruzik { 246b27c8092SJakub Kruzik PetscErrorCode ierr; 247b27c8092SJakub Kruzik 248b27c8092SJakub Kruzik PetscFunctionBegin; 249b27c8092SJakub Kruzik PetscValidHeaderSpecific(pc,PC_CLASSID,1); 250b27c8092SJakub Kruzik ierr = PetscTryMethod(pc,"PCDeflationSetType_C",(PC,PCDeflationType),(pc,type));CHKERRQ(ierr); 251b27c8092SJakub Kruzik PetscFunctionReturn(0); 252b27c8092SJakub Kruzik } 253b27c8092SJakub Kruzik 254b27c8092SJakub Kruzik static PetscErrorCode PCDeflationGetType_Deflation(PC pc,PCDeflationType *type) 255b27c8092SJakub Kruzik { 256b27c8092SJakub Kruzik PC_Deflation *def = (PC_Deflation*)pc->data; 257b27c8092SJakub Kruzik 258b27c8092SJakub Kruzik PetscFunctionBegin; 259b27c8092SJakub Kruzik if (def->init) { 260b27c8092SJakub Kruzik *type = PC_DEFLATION_INIT; 261b27c8092SJakub Kruzik } else if (def->pre) { 262b27c8092SJakub Kruzik *type = PC_DEFLATION_PRE; 263b27c8092SJakub Kruzik } else { 264b27c8092SJakub Kruzik *type = PC_DEFLATION_POST; 265b27c8092SJakub Kruzik } 266b27c8092SJakub Kruzik PetscFunctionReturn(0); 267b27c8092SJakub Kruzik } 268b27c8092SJakub Kruzik 269b27c8092SJakub Kruzik /*@ 270b27c8092SJakub Kruzik PCDeflationGetType - Gets how the diagonal matrix is produced for the preconditioner 271b27c8092SJakub Kruzik 272b27c8092SJakub Kruzik Not Collective on PC 273b27c8092SJakub Kruzik 274b27c8092SJakub Kruzik Input Parameter: 275b27c8092SJakub Kruzik . pc - the preconditioner context 276b27c8092SJakub Kruzik 277b27c8092SJakub Kruzik Output Parameter: 278b27c8092SJakub Kruzik - type - PC_DEFLATION_PRE, PC_DEFLATION_INIT, PC_DEFLATION_POST 279b27c8092SJakub Kruzik 280b27c8092SJakub Kruzik Level: intermediate 281b27c8092SJakub Kruzik 282b27c8092SJakub Kruzik Concepts: Deflation preconditioner 283b27c8092SJakub Kruzik 284b27c8092SJakub Kruzik .seealso: PCDeflationSetType() 285b27c8092SJakub Kruzik @*/ 286b27c8092SJakub Kruzik PetscErrorCode PCDeflationGetType(PC pc,PCDeflationType *type) 287b27c8092SJakub Kruzik { 288b27c8092SJakub Kruzik PetscErrorCode ierr; 289b27c8092SJakub Kruzik 290b27c8092SJakub Kruzik PetscFunctionBegin; 291b27c8092SJakub Kruzik PetscValidHeaderSpecific(pc,PC_CLASSID,1); 292b27c8092SJakub Kruzik ierr = PetscUseMethod(pc,"PCDeflationGetType_C",(PC,PCDeflationType*),(pc,type));CHKERRQ(ierr); 293b27c8092SJakub Kruzik PetscFunctionReturn(0); 294b27c8092SJakub Kruzik } 295b27c8092SJakub Kruzik 29637eeb815SJakub Kruzik static PetscErrorCode PCSetUp_Deflation(PC pc) 29737eeb815SJakub Kruzik { 298409a357bSJakub Kruzik PC_Deflation *def = (PC_Deflation*)pc->data; 299409a357bSJakub Kruzik KSP innerksp; 300409a357bSJakub Kruzik PC pcinner; 301409a357bSJakub Kruzik Mat Amat,nextDef=NULL,*mats; 302409a357bSJakub Kruzik PetscInt i,m,red,size,commsize; 303409a357bSJakub Kruzik PetscBool match,flgspd,transp=PETSC_FALSE; 304409a357bSJakub Kruzik MatCompositeType ctype; 305409a357bSJakub Kruzik MPI_Comm comm; 306409a357bSJakub Kruzik const char *prefix; 30737eeb815SJakub Kruzik PetscErrorCode ierr; 30837eeb815SJakub Kruzik 30937eeb815SJakub Kruzik PetscFunctionBegin; 310409a357bSJakub Kruzik if (pc->setupcalled) PetscFunctionReturn(0); 311409a357bSJakub Kruzik ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 312f8bf229cSJakub Kruzik ierr = PCGetOperators(pc,NULL,&Amat);CHKERRQ(ierr); 313f8bf229cSJakub Kruzik 314f8bf229cSJakub Kruzik /* compute a deflation space */ 315409a357bSJakub Kruzik if (def->W || def->Wt) { 316409a357bSJakub Kruzik def->spacetype = PC_DEFLATION_SPACE_USER; 317409a357bSJakub Kruzik } else { 318409a357bSJakub Kruzik //ierr = KSPDCGComputeDeflationSpace(ksp);CHKERRQ(ierr); 31937eeb815SJakub Kruzik } 32037eeb815SJakub Kruzik 321409a357bSJakub Kruzik /* nested deflation */ 322409a357bSJakub Kruzik if (def->W) { 323409a357bSJakub Kruzik ierr = PetscObjectTypeCompare((PetscObject)def->W,MATCOMPOSITE,&match);CHKERRQ(ierr); 324409a357bSJakub Kruzik if (match) { 325409a357bSJakub Kruzik ierr = MatCompositeGetType(def->W,&ctype);CHKERRQ(ierr); 326409a357bSJakub Kruzik ierr = MatCompositeGetNumberMat(def->W,&size);CHKERRQ(ierr); 32737eeb815SJakub Kruzik } 328409a357bSJakub Kruzik } else { 329409a357bSJakub Kruzik ierr = MatCreateTranspose(def->Wt,&def->W);CHKERRQ(ierr); 330409a357bSJakub Kruzik ierr = PetscObjectTypeCompare((PetscObject)def->Wt,MATCOMPOSITE,&match);CHKERRQ(ierr); 331409a357bSJakub Kruzik if (match) { 332409a357bSJakub Kruzik ierr = MatCompositeGetType(def->Wt,&ctype);CHKERRQ(ierr); 333409a357bSJakub Kruzik ierr = MatCompositeGetNumberMat(def->Wt,&size);CHKERRQ(ierr); 334409a357bSJakub Kruzik } 335409a357bSJakub Kruzik transp = PETSC_TRUE; 336409a357bSJakub Kruzik } 337409a357bSJakub Kruzik if (match && ctype == MAT_COMPOSITE_MULTIPLICATIVE) { 338409a357bSJakub Kruzik ierr = PetscMalloc1(size,&mats);CHKERRQ(ierr); 339409a357bSJakub Kruzik if (!transp) { 340409a357bSJakub Kruzik for (i=0; i<size; i++) { 341409a357bSJakub Kruzik ierr = MatCompositeGetMat(def->W,i,&mats[i]);CHKERRQ(ierr); 342409a357bSJakub Kruzik //ierr = PetscObjectReference((PetscObject)mats[i]);CHKERRQ(ierr); 343409a357bSJakub Kruzik } 344409a357bSJakub Kruzik if (def->nestedlvl < def->maxnestedlvl) { 345409a357bSJakub Kruzik size -= 1; 346409a357bSJakub Kruzik ierr = MatDestroy(&def->W);CHKERRQ(ierr); 347409a357bSJakub Kruzik def->W = mats[size]; 348409a357bSJakub Kruzik ierr = PetscObjectReference((PetscObject)mats[size]);CHKERRQ(ierr); 349409a357bSJakub Kruzik if (size > 1) { 350409a357bSJakub Kruzik ierr = MatCreateComposite(comm,size,mats,&nextDef);CHKERRQ(ierr); 351409a357bSJakub Kruzik ierr = MatCompositeSetType(nextDef,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr); 352409a357bSJakub Kruzik } else { 353409a357bSJakub Kruzik nextDef = mats[0]; 354409a357bSJakub Kruzik ierr = PetscObjectReference((PetscObject)mats[0]);CHKERRQ(ierr); 355409a357bSJakub Kruzik } 356409a357bSJakub Kruzik } else { 357409a357bSJakub Kruzik /* ierr = MatCompositeSetMergeType(def->W,MAT_COMPOSITE_MERGE_LEFT);CHKERRQ(ierr); */ 358409a357bSJakub Kruzik ierr = MatCompositeMerge(def->W);CHKERRQ(ierr); 359409a357bSJakub Kruzik } 360409a357bSJakub Kruzik } 361409a357bSJakub Kruzik ierr = PetscFree(mats);CHKERRQ(ierr); 362409a357bSJakub Kruzik } 363409a357bSJakub Kruzik 364409a357bSJakub Kruzik /* setup coarse problem */ 365409a357bSJakub Kruzik if (!def->WtAWinv) { 366409a357bSJakub Kruzik ierr = MatGetSize(def->W,NULL,&m);CHKERRQ(ierr); /* TODO works for W MatTranspose? */ 367409a357bSJakub Kruzik if (!def->WtAW) { 368409a357bSJakub Kruzik /* TODO add implicit product version ? */ 369409a357bSJakub Kruzik if (!def->AW) { 370409a357bSJakub Kruzik ierr = MatPtAP(Amat,def->W,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtAW);CHKERRQ(ierr); 371409a357bSJakub Kruzik } else { 372409a357bSJakub Kruzik ierr = MatTransposeMatMult(def->W,def->AW,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtAW);CHKERRQ(ierr); 373409a357bSJakub Kruzik } 374409a357bSJakub Kruzik /* TODO create MatInheritOption(Mat,MatOption) */ 375409a357bSJakub Kruzik ierr = MatGetOption(Amat,MAT_SPD,&flgspd);CHKERRQ(ierr); 376409a357bSJakub Kruzik ierr = MatSetOption(def->WtAW,MAT_SPD,flgspd);CHKERRQ(ierr); 377409a357bSJakub Kruzik #if defined(PETSC_USE_DEBUG) 378409a357bSJakub Kruzik /* Check WtAW is not sigular */ 379409a357bSJakub Kruzik PetscReal *norms; 380409a357bSJakub Kruzik ierr = PetscMalloc1(m,&norms);CHKERRQ(ierr); 381409a357bSJakub Kruzik ierr = MatGetColumnNorms(def->WtAW,NORM_INFINITY,norms);CHKERRQ(ierr); 382409a357bSJakub Kruzik for (i=0; i<m; i++) { 383409a357bSJakub Kruzik if (norms[i] < 100*PETSC_MACHINE_EPSILON) { 384409a357bSJakub Kruzik SETERRQ1(comm,PETSC_ERR_SUP,"Column %D of W is in kernel of A.",i); 385409a357bSJakub Kruzik } 386409a357bSJakub Kruzik } 387409a357bSJakub Kruzik ierr = PetscFree(norms);CHKERRQ(ierr); 388409a357bSJakub Kruzik #endif 389409a357bSJakub Kruzik } else { 390409a357bSJakub Kruzik ierr = MatGetOption(def->WtAW,MAT_SPD,&flgspd);CHKERRQ(ierr); 391409a357bSJakub Kruzik } 392409a357bSJakub Kruzik /* TODO use MATINV */ 393ac66f006SJakub Kruzik ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr); 394409a357bSJakub Kruzik ierr = KSPCreate(comm,&def->WtAWinv);CHKERRQ(ierr); 395409a357bSJakub Kruzik ierr = KSPSetOperators(def->WtAWinv,def->WtAW,def->WtAW);CHKERRQ(ierr); 396409a357bSJakub Kruzik ierr = KSPSetType(def->WtAWinv,KSPPREONLY);CHKERRQ(ierr); 397409a357bSJakub Kruzik ierr = KSPGetPC(def->WtAWinv,&pcinner);CHKERRQ(ierr); 398409a357bSJakub Kruzik ierr = PCSetType(pcinner,PCTELESCOPE);CHKERRQ(ierr); 399ac66f006SJakub Kruzik /* ugly hack to not have overwritten PCTELESCOPE */ 400ac66f006SJakub Kruzik if (prefix) { 401ac66f006SJakub Kruzik ierr = KSPSetOptionsPrefix(def->WtAWinv,prefix);CHKERRQ(ierr); 402ac66f006SJakub Kruzik ierr = KSPAppendOptionsPrefix(def->WtAWinv,"tel_");CHKERRQ(ierr); 403ac66f006SJakub Kruzik } else { 404ac66f006SJakub Kruzik ierr = KSPSetOptionsPrefix(def->WtAWinv,"tel_");CHKERRQ(ierr); 405ac66f006SJakub Kruzik } 406ac66f006SJakub Kruzik ierr = PCSetFromOptions(pcinner);CHKERRQ(ierr); 407409a357bSJakub Kruzik /* Reduction factor choice */ 408409a357bSJakub Kruzik red = def->reductionfact; 409409a357bSJakub Kruzik if (red < 0) { 410409a357bSJakub Kruzik ierr = MPI_Comm_size(comm,&commsize);CHKERRQ(ierr); 411409a357bSJakub Kruzik red = ceil((float)commsize/ceil((float)m/commsize)); 412409a357bSJakub Kruzik ierr = PetscObjectTypeCompareAny((PetscObject)(def->WtAW),&match,MATSEQDENSE,MATMPIDENSE,MATDENSE,"");CHKERRQ(ierr); 413409a357bSJakub Kruzik if (match) red = commsize; 414409a357bSJakub Kruzik ierr = PetscInfo1(pc,"Auto choosing reduction factor %D\n",red);CHKERRQ(ierr); /* TODO add level? */ 415409a357bSJakub Kruzik } 416409a357bSJakub Kruzik ierr = PCTelescopeSetReductionFactor(pcinner,red);CHKERRQ(ierr); 417ac66f006SJakub Kruzik ierr = PCSetUp(pcinner);CHKERRQ(ierr); 418409a357bSJakub Kruzik ierr = PCTelescopeGetKSP(pcinner,&innerksp);CHKERRQ(ierr); 419ac66f006SJakub Kruzik if (innerksp) { 420409a357bSJakub Kruzik ierr = KSPGetPC(innerksp,&pcinner);CHKERRQ(ierr); 421409a357bSJakub Kruzik /* Setup KSP and PC */ 422409a357bSJakub Kruzik if (nextDef) { /* next level for multilevel deflation */ 423409a357bSJakub Kruzik /* set default KSPtype */ 424409a357bSJakub Kruzik if (!def->ksptype) { 425409a357bSJakub Kruzik def->ksptype = KSPFGMRES; 426409a357bSJakub Kruzik if (flgspd) { /* SPD system */ 427409a357bSJakub Kruzik def->ksptype = KSPFCG; 428409a357bSJakub Kruzik } 429409a357bSJakub Kruzik } 430409a357bSJakub Kruzik ierr = KSPSetType(innerksp,def->ksptype);CHKERRQ(ierr); /* TODO iherit from KSP */ 431409a357bSJakub Kruzik ierr = PCSetType(pcinner,PCDEFLATION);CHKERRQ(ierr); /* TODO create coarse preconditinoner M_c = WtMW ? */ 432409a357bSJakub Kruzik ierr = PCDeflationSetSpace(pcinner,nextDef,transp);CHKERRQ(ierr); 4335c0c31c5SJakub Kruzik ierr = PCDeflationSetLvl_Deflation(pcinner,def->nestedlvl+1,def->maxnestedlvl);CHKERRQ(ierr); 434409a357bSJakub Kruzik /* inherit options TODO if not set */ 435409a357bSJakub Kruzik ((PC_Deflation*)(pcinner))->ksptype = def->ksptype; 436409a357bSJakub Kruzik ((PC_Deflation*)(pcinner))->correct = def->correct; 437409a357bSJakub Kruzik ((PC_Deflation*)(pcinner))->adaptiveconv = def->adaptiveconv; 438409a357bSJakub Kruzik ((PC_Deflation*)(pcinner))->adaptiveconst = def->adaptiveconst; 439409a357bSJakub Kruzik ierr = MatDestroy(&nextDef);CHKERRQ(ierr); 440409a357bSJakub Kruzik } else { /* the last level */ 441ac66f006SJakub Kruzik ierr = KSPSetType(innerksp,KSPGMRES);CHKERRQ(ierr); 442409a357bSJakub Kruzik /* TODO Cholesky if flgspd? */ 443ac66f006SJakub Kruzik ierr = PCSetType(pcinner,PCLU);CHKERRQ(ierr); 444409a357bSJakub Kruzik //TODO remove explicit matSolverPackage 445409a357bSJakub Kruzik if (commsize == red) { 446ac66f006SJakub Kruzik ierr = PCFactorSetMatSolverType(pcinner,MATSOLVERSUPERLU);CHKERRQ(ierr); 447409a357bSJakub Kruzik } else { 448ac66f006SJakub Kruzik ierr = PCFactorSetMatSolverType(pcinner,MATSOLVERSUPERLU_DIST);CHKERRQ(ierr); 449409a357bSJakub Kruzik } 450409a357bSJakub Kruzik } 451409a357bSJakub Kruzik /* TODO use def_[lvl]_ if lvl > 0? */ 452409a357bSJakub Kruzik if (prefix) { 453409a357bSJakub Kruzik ierr = KSPSetOptionsPrefix(innerksp,prefix);CHKERRQ(ierr); 454409a357bSJakub Kruzik ierr = KSPAppendOptionsPrefix(innerksp,"def_");CHKERRQ(ierr); 455409a357bSJakub Kruzik } else { 456409a357bSJakub Kruzik ierr = KSPSetOptionsPrefix(innerksp,"def_");CHKERRQ(ierr); 457409a357bSJakub Kruzik } 458ac66f006SJakub Kruzik } 459409a357bSJakub Kruzik /* TODO: check if WtAWinv is KSP and move following from this if */ 460409a357bSJakub Kruzik ierr = KSPSetFromOptions(def->WtAWinv);CHKERRQ(ierr); 461409a357bSJakub Kruzik //if (def->adaptiveconv) { 462409a357bSJakub Kruzik // PetscReal *rnorm; 463409a357bSJakub Kruzik // PetscNew(&rnorm); 464409a357bSJakub Kruzik // ierr = KSPSetConvergenceTest(def->WtAWinv,KSPDCGConvergedAdaptive_DCG,rnorm,NULL);CHKERRQ(ierr); 465409a357bSJakub Kruzik //} 466409a357bSJakub Kruzik ierr = KSPSetUp(def->WtAWinv);CHKERRQ(ierr); 467409a357bSJakub Kruzik } 468409a357bSJakub Kruzik 469f8bf229cSJakub Kruzik /* create work vecs */ 470f8bf229cSJakub Kruzik ierr = MatCreateVecs(Amat,NULL,&def->work);CHKERRQ(ierr); 471f8bf229cSJakub Kruzik ierr = KSPCreateVecs(def->WtAWinv,2,&def->workcoarse,0,NULL);CHKERRQ(ierr); 47237eeb815SJakub Kruzik PetscFunctionReturn(0); 47337eeb815SJakub Kruzik } 47437eeb815SJakub Kruzik /* -------------------------------------------------------------------------- */ 47537eeb815SJakub Kruzik static PetscErrorCode PCReset_Deflation(PC pc) 47637eeb815SJakub Kruzik { 47737eeb815SJakub Kruzik PC_Deflation *jac = (PC_Deflation*)pc->data; 47837eeb815SJakub Kruzik PetscErrorCode ierr; 47937eeb815SJakub Kruzik 48037eeb815SJakub Kruzik PetscFunctionBegin; 48137eeb815SJakub Kruzik PetscFunctionReturn(0); 48237eeb815SJakub Kruzik } 48337eeb815SJakub Kruzik 48437eeb815SJakub Kruzik /* 48537eeb815SJakub Kruzik PCDestroy_Deflation - Destroys the private context for the Deflation preconditioner 48637eeb815SJakub Kruzik that was created with PCCreate_Deflation(). 48737eeb815SJakub Kruzik 48837eeb815SJakub Kruzik Input Parameter: 48937eeb815SJakub Kruzik . pc - the preconditioner context 49037eeb815SJakub Kruzik 49137eeb815SJakub Kruzik Application Interface Routine: PCDestroy() 49237eeb815SJakub Kruzik */ 49337eeb815SJakub Kruzik static PetscErrorCode PCDestroy_Deflation(PC pc) 49437eeb815SJakub Kruzik { 49537eeb815SJakub Kruzik PetscErrorCode ierr; 49637eeb815SJakub Kruzik 49737eeb815SJakub Kruzik PetscFunctionBegin; 49837eeb815SJakub Kruzik ierr = PCReset_Deflation(pc);CHKERRQ(ierr); 49937eeb815SJakub Kruzik 50037eeb815SJakub Kruzik /* 50137eeb815SJakub Kruzik Free the private data structure that was hanging off the PC 50237eeb815SJakub Kruzik */ 50337eeb815SJakub Kruzik ierr = PetscFree(pc->data);CHKERRQ(ierr); 50437eeb815SJakub Kruzik PetscFunctionReturn(0); 50537eeb815SJakub Kruzik } 50637eeb815SJakub Kruzik 50737eeb815SJakub Kruzik static PetscErrorCode PCSetFromOptions_Deflation(PetscOptionItems *PetscOptionsObject,PC pc) 50837eeb815SJakub Kruzik { 50937eeb815SJakub Kruzik PC_Deflation *jac = (PC_Deflation*)pc->data; 51037eeb815SJakub Kruzik PetscErrorCode ierr; 51137eeb815SJakub Kruzik PetscBool flg; 51237eeb815SJakub Kruzik PCDeflationType deflt,type; 51337eeb815SJakub Kruzik 51437eeb815SJakub Kruzik PetscFunctionBegin; 51537eeb815SJakub Kruzik ierr = PCDeflationGetType(pc,&deflt);CHKERRQ(ierr); 51637eeb815SJakub Kruzik ierr = PetscOptionsHead(PetscOptionsObject,"Deflation options");CHKERRQ(ierr); 51737eeb815SJakub Kruzik ierr = PetscOptionsEnum("-pc_jacobi_type","How to construct diagonal matrix","PCDeflationSetType",PCDeflationTypes,(PetscEnum)deflt,(PetscEnum*)&type,&flg);CHKERRQ(ierr); 51837eeb815SJakub Kruzik if (flg) { 51937eeb815SJakub Kruzik ierr = PCDeflationSetType(pc,type);CHKERRQ(ierr); 52037eeb815SJakub Kruzik } 52137eeb815SJakub Kruzik ierr = PetscOptionsTail();CHKERRQ(ierr); 52237eeb815SJakub Kruzik PetscFunctionReturn(0); 52337eeb815SJakub Kruzik } 52437eeb815SJakub Kruzik 52537eeb815SJakub Kruzik /*MC 526e361f8b5SJakub Kruzik PCDEFLATION - Deflation preconditioner shifts part of the spectrum to zero (deflates) 527e361f8b5SJakub Kruzik or to a predefined value 52837eeb815SJakub Kruzik 52937eeb815SJakub Kruzik Options Database Key: 530e361f8b5SJakub Kruzik + -pc_deflation_type <init,pre,post> - selects approach to deflation (default: pre) 53137eeb815SJakub Kruzik - -pc_jacobi_abs - use the absolute value of the diagonal entry 53237eeb815SJakub Kruzik 53337eeb815SJakub Kruzik Level: beginner 53437eeb815SJakub Kruzik 53537eeb815SJakub Kruzik Notes: 536e361f8b5SJakub Kruzik todo 53737eeb815SJakub Kruzik 53837eeb815SJakub Kruzik .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 539e662bc50SJakub Kruzik PCDeflationSetType(), PCDeflationSetSpace() 54037eeb815SJakub Kruzik M*/ 54137eeb815SJakub Kruzik 54237eeb815SJakub Kruzik PETSC_EXTERN PetscErrorCode PCCreate_Deflation(PC pc) 54337eeb815SJakub Kruzik { 54437eeb815SJakub Kruzik PC_Deflation *def; 54537eeb815SJakub Kruzik PetscErrorCode ierr; 54637eeb815SJakub Kruzik 54737eeb815SJakub Kruzik PetscFunctionBegin; 54837eeb815SJakub Kruzik ierr = PetscNewLog(pc,&def);CHKERRQ(ierr); 54937eeb815SJakub Kruzik pc->data = (void*)def; 55037eeb815SJakub Kruzik 551e662bc50SJakub Kruzik def->init = PETSC_FALSE; 552e662bc50SJakub Kruzik def->pre = PETSC_TRUE; 553e662bc50SJakub Kruzik def->correct = PETSC_FALSE; 554e662bc50SJakub Kruzik def->truenorm = PETSC_TRUE; 555e662bc50SJakub Kruzik def->reductionfact = -1; 556e662bc50SJakub Kruzik def->spacetype = PC_DEFLATION_SPACE_HAAR; 557e662bc50SJakub Kruzik def->spacesize = 1; 558e662bc50SJakub Kruzik def->extendsp = PETSC_FALSE; 559e662bc50SJakub Kruzik def->nestedlvl = 0; 560e662bc50SJakub Kruzik def->maxnestedlvl = 0; 561e662bc50SJakub Kruzik def->adaptiveconv = PETSC_FALSE; 562e662bc50SJakub Kruzik def->adaptiveconst = 1.0; 56337eeb815SJakub Kruzik 56437eeb815SJakub Kruzik /* 56537eeb815SJakub Kruzik Set the pointers for the functions that are provided above. 56637eeb815SJakub Kruzik Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 56737eeb815SJakub Kruzik are called, they will automatically call these functions. Note we 56837eeb815SJakub Kruzik choose not to provide a couple of these functions since they are 56937eeb815SJakub Kruzik not needed. 57037eeb815SJakub Kruzik */ 57137eeb815SJakub Kruzik pc->ops->apply = PCApply_Deflation; 57237eeb815SJakub Kruzik pc->ops->applytranspose = PCApply_Deflation; 573b27c8092SJakub Kruzik pc->ops->presolve = PCPreSolve_Deflation; 574b27c8092SJakub Kruzik pc->ops->postsolve = 0; 57537eeb815SJakub Kruzik pc->ops->setup = PCSetUp_Deflation; 57637eeb815SJakub Kruzik pc->ops->reset = PCReset_Deflation; 57737eeb815SJakub Kruzik pc->ops->destroy = PCDestroy_Deflation; 57837eeb815SJakub Kruzik pc->ops->setfromoptions = PCSetFromOptions_Deflation; 57937eeb815SJakub Kruzik pc->ops->view = 0; 58037eeb815SJakub Kruzik pc->ops->applyrichardson = 0; 58137eeb815SJakub Kruzik pc->ops->applysymmetricleft = 0; 58237eeb815SJakub Kruzik pc->ops->applysymmetricright = 0; 58337eeb815SJakub Kruzik 58437eeb815SJakub Kruzik ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetType_C",PCDeflationSetType_Deflation);CHKERRQ(ierr); 58537eeb815SJakub Kruzik ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationGetType_C",PCDeflationGetType_Deflation);CHKERRQ(ierr); 586e662bc50SJakub Kruzik ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetSpace_C",PCDeflationSetSpace_Deflation);CHKERRQ(ierr); 5875c0c31c5SJakub Kruzik ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetLvl_C",PCDeflationSetLvl_Deflation);CHKERRQ(ierr); 58837eeb815SJakub Kruzik PetscFunctionReturn(0); 58937eeb815SJakub Kruzik } 59037eeb815SJakub Kruzik 591