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; 69*f8bf229cSJakub Kruzik Vec work; 70*f8bf229cSJakub 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 7937eeb815SJakub Kruzik static PetscErrorCode PCDeflationSetType_Deflation(PC pc,PCDeflationType type) 8037eeb815SJakub Kruzik { 8137eeb815SJakub Kruzik PC_Deflation *def = (PC_Deflation*)pc->data; 8237eeb815SJakub Kruzik 8337eeb815SJakub Kruzik PetscFunctionBegin; 8437eeb815SJakub Kruzik def->init = PETSC_FALSE; 8537eeb815SJakub Kruzik def->pre = PETSC_FALSE; 8637eeb815SJakub Kruzik if (type == PC_DEFLATION_INIT) { 8737eeb815SJakub Kruzik def->init = PETSC_TRUE; 8837eeb815SJakub Kruzik def->pre = PETSC_TRUE; 8937eeb815SJakub Kruzik } else if (type == PC_DEFLATION_PRE) { 9037eeb815SJakub Kruzik def->pre = PETSC_TRUE; 9137eeb815SJakub Kruzik } 9237eeb815SJakub Kruzik PetscFunctionReturn(0); 9337eeb815SJakub Kruzik } 9437eeb815SJakub Kruzik 95bb9edbfaSJakub Kruzik /*@ 96bb9edbfaSJakub Kruzik PCDeflationSetType - Causes the deflation preconditioner to use only a special 97bb9edbfaSJakub Kruzik initial gues or pre/post solve solution update 98bb9edbfaSJakub Kruzik 99bb9edbfaSJakub Kruzik Logically Collective on PC 100bb9edbfaSJakub Kruzik 101bb9edbfaSJakub Kruzik Input Parameters: 102bb9edbfaSJakub Kruzik + pc - the preconditioner context 103bb9edbfaSJakub Kruzik - type - PC_DEFLATION_PRE, PC_DEFLATION_INIT, PC_DEFLATION_POST 104bb9edbfaSJakub Kruzik 105bb9edbfaSJakub Kruzik Options Database Key: 106bb9edbfaSJakub Kruzik . -pc_deflation_type <pre,init,post> 107bb9edbfaSJakub Kruzik 108bb9edbfaSJakub Kruzik Level: intermediate 109bb9edbfaSJakub Kruzik 110bb9edbfaSJakub Kruzik Concepts: Deflation preconditioner 111bb9edbfaSJakub Kruzik 112bb9edbfaSJakub Kruzik .seealso: PCDeflationGetType() 113bb9edbfaSJakub Kruzik @*/ 114bb9edbfaSJakub Kruzik PetscErrorCode PCDeflationSetType(PC pc,PCDeflationType type) 115bb9edbfaSJakub Kruzik { 116bb9edbfaSJakub Kruzik PetscErrorCode ierr; 117bb9edbfaSJakub Kruzik 118bb9edbfaSJakub Kruzik PetscFunctionBegin; 119bb9edbfaSJakub Kruzik PetscValidHeaderSpecific(pc,PC_CLASSID,1); 120bb9edbfaSJakub Kruzik ierr = PetscTryMethod(pc,"PCDeflationSetType_C",(PC,PCDeflationType),(pc,type));CHKERRQ(ierr); 121bb9edbfaSJakub Kruzik PetscFunctionReturn(0); 122bb9edbfaSJakub Kruzik } 123bb9edbfaSJakub Kruzik 12437eeb815SJakub Kruzik static PetscErrorCode PCDeflationGetType_Deflation(PC pc,PCDeflationType *type) 12537eeb815SJakub Kruzik { 12637eeb815SJakub Kruzik PC_Deflation *def = (PC_Deflation*)pc->data; 12737eeb815SJakub Kruzik 12837eeb815SJakub Kruzik PetscFunctionBegin; 12937eeb815SJakub Kruzik if (def->init) { 13037eeb815SJakub Kruzik *type = PC_DEFLATION_INIT; 13137eeb815SJakub Kruzik } else if (def->pre) { 13237eeb815SJakub Kruzik *type = PC_DEFLATION_PRE; 13337eeb815SJakub Kruzik } else { 13437eeb815SJakub Kruzik *type = PC_DEFLATION_POST; 13537eeb815SJakub Kruzik } 13637eeb815SJakub Kruzik PetscFunctionReturn(0); 13737eeb815SJakub Kruzik } 13837eeb815SJakub Kruzik 139bb9edbfaSJakub Kruzik /*@ 140bb9edbfaSJakub Kruzik PCDeflationGetType - Gets how the diagonal matrix is produced for the preconditioner 141bb9edbfaSJakub Kruzik 142bb9edbfaSJakub Kruzik Not Collective on PC 143bb9edbfaSJakub Kruzik 144bb9edbfaSJakub Kruzik Input Parameter: 145bb9edbfaSJakub Kruzik . pc - the preconditioner context 146bb9edbfaSJakub Kruzik 147bb9edbfaSJakub Kruzik Output Parameter: 148bb9edbfaSJakub Kruzik - type - PC_DEFLATION_PRE, PC_DEFLATION_INIT, PC_DEFLATION_POST 149bb9edbfaSJakub Kruzik 150bb9edbfaSJakub Kruzik Level: intermediate 151bb9edbfaSJakub Kruzik 152bb9edbfaSJakub Kruzik Concepts: Deflation preconditioner 153bb9edbfaSJakub Kruzik 154bb9edbfaSJakub Kruzik .seealso: PCDeflationSetType() 155bb9edbfaSJakub Kruzik @*/ 156bb9edbfaSJakub Kruzik PetscErrorCode PCDeflationGetType(PC pc,PCDeflationType *type) 157bb9edbfaSJakub Kruzik { 158bb9edbfaSJakub Kruzik PetscErrorCode ierr; 159bb9edbfaSJakub Kruzik 160bb9edbfaSJakub Kruzik PetscFunctionBegin; 161bb9edbfaSJakub Kruzik PetscValidHeaderSpecific(pc,PC_CLASSID,1); 162bb9edbfaSJakub Kruzik ierr = PetscUseMethod(pc,"PCDeflationGetType_C",(PC,PCDeflationType*),(pc,type));CHKERRQ(ierr); 163bb9edbfaSJakub Kruzik PetscFunctionReturn(0); 164bb9edbfaSJakub Kruzik } 165bb9edbfaSJakub Kruzik 166e662bc50SJakub Kruzik static PetscErrorCode PCDeflationSetSpace_Deflation(PC pc,Mat W,PetscBool transpose) 167e662bc50SJakub Kruzik { 168e662bc50SJakub Kruzik PC_Deflation *def = (PC_Deflation*)pc->data; 169e662bc50SJakub Kruzik PetscErrorCode ierr; 170e662bc50SJakub Kruzik 171e662bc50SJakub Kruzik PetscFunctionBegin; 172e662bc50SJakub Kruzik if (transpose) { 173e662bc50SJakub Kruzik def->Wt = W; 174e662bc50SJakub Kruzik def->W = NULL; 175e662bc50SJakub Kruzik } else { 176e662bc50SJakub Kruzik def->W = W; 177e662bc50SJakub Kruzik } 178e662bc50SJakub Kruzik ierr = PetscObjectReference((PetscObject)W);CHKERRQ(ierr); 179e662bc50SJakub Kruzik PetscFunctionReturn(0); 180e662bc50SJakub Kruzik } 181e662bc50SJakub Kruzik 182e662bc50SJakub Kruzik /* TODO create PCDeflationSetSpaceTranspose? */ 183e662bc50SJakub Kruzik /*@ 184e662bc50SJakub Kruzik PCDeflationSetSpace - Set deflation space matrix (or its transpose). 185e662bc50SJakub Kruzik 186e662bc50SJakub Kruzik Logically Collective on PC 187e662bc50SJakub Kruzik 188e662bc50SJakub Kruzik Input Parameters: 189e662bc50SJakub Kruzik + pc - the preconditioner context 190e662bc50SJakub Kruzik . W - deflation matrix 191e662bc50SJakub Kruzik - tranpose - indicates that W is an explicit transpose of the deflation matrix 192e662bc50SJakub Kruzik 193e662bc50SJakub Kruzik Level: intermediate 194e662bc50SJakub Kruzik 195e662bc50SJakub Kruzik .seealso: PCDEFLATION 196e662bc50SJakub Kruzik @*/ 197e662bc50SJakub Kruzik PetscErrorCode PCDeflationSetSpace(PC pc,Mat W,PetscBool transpose) 198e662bc50SJakub Kruzik { 199e662bc50SJakub Kruzik PetscErrorCode ierr; 200e662bc50SJakub Kruzik 201e662bc50SJakub Kruzik PetscFunctionBegin; 202e662bc50SJakub Kruzik PetscValidHeaderSpecific(pc,PC_CLASSID,1); 203e662bc50SJakub Kruzik PetscValidHeaderSpecific(W,MAT_CLASSID,2); 204e662bc50SJakub Kruzik PetscValidLogicalCollectiveBool(pc,transpose,3); 205e662bc50SJakub Kruzik ierr = PetscTryMethod(pc,"PCDeflationSetSpace_C",(PC,Mat,PetscBool),(pc,W,transpose));CHKERRQ(ierr); 206e662bc50SJakub Kruzik PetscFunctionReturn(0); 207e662bc50SJakub Kruzik } 208e662bc50SJakub Kruzik 2095c0c31c5SJakub Kruzik static PetscErrorCode PCDeflationSetLvl_Deflation(PC pc,PetscInt current,PetscInt max) 2105c0c31c5SJakub Kruzik { 2115c0c31c5SJakub Kruzik PC_Deflation *def = (PC_Deflation*)pc->data; 2125c0c31c5SJakub Kruzik 2135c0c31c5SJakub Kruzik PetscFunctionBegin; 2145c0c31c5SJakub Kruzik def->nestedlvl = current; 2155c0c31c5SJakub Kruzik def->maxnestedlvl = max; 2165c0c31c5SJakub Kruzik PetscFunctionReturn(0); 2175c0c31c5SJakub Kruzik } 2185c0c31c5SJakub Kruzik 2195c0c31c5SJakub Kruzik /*@ 2205c0c31c5SJakub Kruzik PCDeflationSetMaxLvl - Set maximum level of deflation. 2215c0c31c5SJakub Kruzik 2225c0c31c5SJakub Kruzik Logically Collective on PC 2235c0c31c5SJakub Kruzik 2245c0c31c5SJakub Kruzik Input Parameters: 2255c0c31c5SJakub Kruzik + pc - the preconditioner context 2265c0c31c5SJakub Kruzik . max - maximum deflation level 2275c0c31c5SJakub Kruzik 2285c0c31c5SJakub Kruzik Level: intermediate 2295c0c31c5SJakub Kruzik 2305c0c31c5SJakub Kruzik .seealso: PCDEFLATION 2315c0c31c5SJakub Kruzik @*/ 2325c0c31c5SJakub Kruzik PetscErrorCode PCDeflationSetMaxLvl(PC pc,PetscInt max) 2335c0c31c5SJakub Kruzik { 2345c0c31c5SJakub Kruzik PetscErrorCode ierr; 2355c0c31c5SJakub Kruzik 2365c0c31c5SJakub Kruzik PetscFunctionBegin; 2375c0c31c5SJakub Kruzik PetscValidLogicalCollectiveInt(pc,max,2); 2385c0c31c5SJakub Kruzik ierr = PetscTryMethod(pc,"PCDeflationSetLvl_C",(PC,PetscInt,PetscInt),(pc,0,max));CHKERRQ(ierr); 2395c0c31c5SJakub Kruzik PetscFunctionReturn(0); 2405c0c31c5SJakub Kruzik } 241e662bc50SJakub Kruzik 24237eeb815SJakub Kruzik /* 24337eeb815SJakub Kruzik PCSetUp_Deflation - Prepares for the use of the Deflation preconditioner 24437eeb815SJakub Kruzik by setting data structures and options. 24537eeb815SJakub Kruzik 24637eeb815SJakub Kruzik Input Parameter: 24737eeb815SJakub Kruzik . pc - the preconditioner context 24837eeb815SJakub Kruzik 24937eeb815SJakub Kruzik Application Interface Routine: PCSetUp() 25037eeb815SJakub Kruzik 251*f8bf229cSJakub Kruzik /* 252*f8bf229cSJakub Kruzik if (def->correct) { 253*f8bf229cSJakub Kruzik z <- r - W*(W'*A*W)^{-1}*W'*(A*r -r) = (P-Q)*r 254*f8bf229cSJakub Kruzik } else { 255*f8bf229cSJakub Kruzik z <- r - W*(W'*A*W)^{-1}*W'*A*r = P*r 256*f8bf229cSJakub Kruzik } 25737eeb815SJakub Kruzik */ 258*f8bf229cSJakub Kruzik static PetscErrorCode PCApply_Deflation(PC pc,Vec r, Vec z) 259*f8bf229cSJakub Kruzik { 260*f8bf229cSJakub Kruzik PC_Deflation *def = (PC_Deflation*)pc->data; 261*f8bf229cSJakub Kruzik Mat A; 262*f8bf229cSJakub Kruzik Vec u,w1,w2; 263*f8bf229cSJakub Kruzik PetscErrorCode ierr; 264*f8bf229cSJakub Kruzik 265*f8bf229cSJakub Kruzik PetscFunctionBegin; 266*f8bf229cSJakub Kruzik w1 = def->workcoarse[0]; 267*f8bf229cSJakub Kruzik w2 = def->workcoarse[1]; 268*f8bf229cSJakub Kruzik u = def->work; 269*f8bf229cSJakub Kruzik ierr = PCGetOperators(pc,NULL,&A);CHKERRQ(ierr); 270*f8bf229cSJakub Kruzik 271*f8bf229cSJakub Kruzik if (!def->AW) { 272*f8bf229cSJakub Kruzik ierr = MatMult(A,r,u);CHKERRQ(ierr); /* u <- A*r */ 273*f8bf229cSJakub Kruzik if (def->correct) ierr = VecAXPY(u,-1.0,r);CHKERRQ(ierr); /* u <- A*r -r */ 274*f8bf229cSJakub Kruzik ierr = MatMultTranspose(def->W,u,w1);CHKERRQ(ierr); /* w1 <- W'*u */ 275*f8bf229cSJakub Kruzik } else { 276*f8bf229cSJakub Kruzik ierr = MatMultTranspose(def->AW,u,w1);CHKERRQ(ierr); /* u <- A*r */ 277*f8bf229cSJakub Kruzik if (def->correct) { 278*f8bf229cSJakub Kruzik ierr = MatMultTranspose(def->W,r,w2);CHKERRQ(ierr); /* w2 <- W'*u */ 279*f8bf229cSJakub Kruzik ierr = VecAXPY(w1,-1.0,w2);CHKERRQ(ierr); /* w1 <- w1 - w2 */ 280*f8bf229cSJakub Kruzik } 281*f8bf229cSJakub Kruzik } 282*f8bf229cSJakub Kruzik ierr = KSPSolve(def->WtAWinv,w1,w2);CHKERRQ(ierr); /* w2 <- (W'*A*W)^{-1}*w1 */ 283*f8bf229cSJakub Kruzik ierr = MatMult(def->W,w2,u);CHKERRQ(ierr); /* u <- W*w2 */ 284*f8bf229cSJakub Kruzik ierr = VecWAXPY(z,-1.0,u,r);CHKERRQ(ierr); /* z <- r - u */ 285*f8bf229cSJakub Kruzik PetscFunctionReturn(0); 286*f8bf229cSJakub Kruzik } 287*f8bf229cSJakub Kruzik 28837eeb815SJakub Kruzik static PetscErrorCode PCSetUp_Deflation(PC pc) 28937eeb815SJakub Kruzik { 290409a357bSJakub Kruzik PC_Deflation *def = (PC_Deflation*)pc->data; 291409a357bSJakub Kruzik KSP innerksp; 292409a357bSJakub Kruzik PC pcinner; 293409a357bSJakub Kruzik Mat Amat,nextDef=NULL,*mats; 294409a357bSJakub Kruzik PetscInt i,m,red,size,commsize; 295409a357bSJakub Kruzik PetscBool match,flgspd,transp=PETSC_FALSE; 296409a357bSJakub Kruzik MatCompositeType ctype; 297409a357bSJakub Kruzik MPI_Comm comm; 298409a357bSJakub Kruzik const char *prefix; 29937eeb815SJakub Kruzik PetscErrorCode ierr; 30037eeb815SJakub Kruzik 30137eeb815SJakub Kruzik PetscFunctionBegin; 302409a357bSJakub Kruzik if (pc->setupcalled) PetscFunctionReturn(0); 303409a357bSJakub Kruzik ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 304*f8bf229cSJakub Kruzik ierr = PCGetOperators(pc,NULL,&Amat);CHKERRQ(ierr); 305*f8bf229cSJakub Kruzik 306*f8bf229cSJakub Kruzik /* compute a deflation space */ 307409a357bSJakub Kruzik if (def->W || def->Wt) { 308409a357bSJakub Kruzik def->spacetype = PC_DEFLATION_SPACE_USER; 309409a357bSJakub Kruzik } else { 310409a357bSJakub Kruzik //ierr = KSPDCGComputeDeflationSpace(ksp);CHKERRQ(ierr); 31137eeb815SJakub Kruzik } 31237eeb815SJakub Kruzik 313409a357bSJakub Kruzik /* nested deflation */ 314409a357bSJakub Kruzik if (def->W) { 315409a357bSJakub Kruzik ierr = PetscObjectTypeCompare((PetscObject)def->W,MATCOMPOSITE,&match);CHKERRQ(ierr); 316409a357bSJakub Kruzik if (match) { 317409a357bSJakub Kruzik ierr = MatCompositeGetType(def->W,&ctype);CHKERRQ(ierr); 318409a357bSJakub Kruzik ierr = MatCompositeGetNumberMat(def->W,&size);CHKERRQ(ierr); 31937eeb815SJakub Kruzik } 320409a357bSJakub Kruzik } else { 321409a357bSJakub Kruzik ierr = MatCreateTranspose(def->Wt,&def->W);CHKERRQ(ierr); 322409a357bSJakub Kruzik ierr = PetscObjectTypeCompare((PetscObject)def->Wt,MATCOMPOSITE,&match);CHKERRQ(ierr); 323409a357bSJakub Kruzik if (match) { 324409a357bSJakub Kruzik ierr = MatCompositeGetType(def->Wt,&ctype);CHKERRQ(ierr); 325409a357bSJakub Kruzik ierr = MatCompositeGetNumberMat(def->Wt,&size);CHKERRQ(ierr); 326409a357bSJakub Kruzik } 327409a357bSJakub Kruzik transp = PETSC_TRUE; 328409a357bSJakub Kruzik } 329409a357bSJakub Kruzik if (match && ctype == MAT_COMPOSITE_MULTIPLICATIVE) { 330409a357bSJakub Kruzik ierr = PetscMalloc1(size,&mats);CHKERRQ(ierr); 331409a357bSJakub Kruzik if (!transp) { 332409a357bSJakub Kruzik for (i=0; i<size; i++) { 333409a357bSJakub Kruzik ierr = MatCompositeGetMat(def->W,i,&mats[i]);CHKERRQ(ierr); 334409a357bSJakub Kruzik //ierr = PetscObjectReference((PetscObject)mats[i]);CHKERRQ(ierr); 335409a357bSJakub Kruzik } 336409a357bSJakub Kruzik if (def->nestedlvl < def->maxnestedlvl) { 337409a357bSJakub Kruzik size -= 1; 338409a357bSJakub Kruzik ierr = MatDestroy(&def->W);CHKERRQ(ierr); 339409a357bSJakub Kruzik def->W = mats[size]; 340409a357bSJakub Kruzik ierr = PetscObjectReference((PetscObject)mats[size]);CHKERRQ(ierr); 341409a357bSJakub Kruzik if (size > 1) { 342409a357bSJakub Kruzik ierr = MatCreateComposite(comm,size,mats,&nextDef);CHKERRQ(ierr); 343409a357bSJakub Kruzik ierr = MatCompositeSetType(nextDef,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr); 344409a357bSJakub Kruzik } else { 345409a357bSJakub Kruzik nextDef = mats[0]; 346409a357bSJakub Kruzik ierr = PetscObjectReference((PetscObject)mats[0]);CHKERRQ(ierr); 347409a357bSJakub Kruzik } 348409a357bSJakub Kruzik } else { 349409a357bSJakub Kruzik /* ierr = MatCompositeSetMergeType(def->W,MAT_COMPOSITE_MERGE_LEFT);CHKERRQ(ierr); */ 350409a357bSJakub Kruzik ierr = MatCompositeMerge(def->W);CHKERRQ(ierr); 351409a357bSJakub Kruzik } 352409a357bSJakub Kruzik } 353409a357bSJakub Kruzik ierr = PetscFree(mats);CHKERRQ(ierr); 354409a357bSJakub Kruzik } 355409a357bSJakub Kruzik 356409a357bSJakub Kruzik /* setup coarse problem */ 357409a357bSJakub Kruzik if (!def->WtAWinv) { 358409a357bSJakub Kruzik ierr = MatGetSize(def->W,NULL,&m);CHKERRQ(ierr); /* TODO works for W MatTranspose? */ 359409a357bSJakub Kruzik if (!def->WtAW) { 360409a357bSJakub Kruzik /* TODO add implicit product version ? */ 361409a357bSJakub Kruzik if (!def->AW) { 362409a357bSJakub Kruzik ierr = MatPtAP(Amat,def->W,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtAW);CHKERRQ(ierr); 363409a357bSJakub Kruzik } else { 364409a357bSJakub Kruzik ierr = MatTransposeMatMult(def->W,def->AW,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtAW);CHKERRQ(ierr); 365409a357bSJakub Kruzik } 366409a357bSJakub Kruzik /* TODO create MatInheritOption(Mat,MatOption) */ 367409a357bSJakub Kruzik ierr = MatGetOption(Amat,MAT_SPD,&flgspd);CHKERRQ(ierr); 368409a357bSJakub Kruzik ierr = MatSetOption(def->WtAW,MAT_SPD,flgspd);CHKERRQ(ierr); 369409a357bSJakub Kruzik #if defined(PETSC_USE_DEBUG) 370409a357bSJakub Kruzik /* Check WtAW is not sigular */ 371409a357bSJakub Kruzik PetscReal *norms; 372409a357bSJakub Kruzik ierr = PetscMalloc1(m,&norms);CHKERRQ(ierr); 373409a357bSJakub Kruzik ierr = MatGetColumnNorms(def->WtAW,NORM_INFINITY,norms);CHKERRQ(ierr); 374409a357bSJakub Kruzik for (i=0; i<m; i++) { 375409a357bSJakub Kruzik if (norms[i] < 100*PETSC_MACHINE_EPSILON) { 376409a357bSJakub Kruzik SETERRQ1(comm,PETSC_ERR_SUP,"Column %D of W is in kernel of A.",i); 377409a357bSJakub Kruzik } 378409a357bSJakub Kruzik } 379409a357bSJakub Kruzik ierr = PetscFree(norms);CHKERRQ(ierr); 380409a357bSJakub Kruzik #endif 381409a357bSJakub Kruzik } else { 382409a357bSJakub Kruzik ierr = MatGetOption(def->WtAW,MAT_SPD,&flgspd);CHKERRQ(ierr); 383409a357bSJakub Kruzik } 384409a357bSJakub Kruzik /* TODO use MATINV */ 385409a357bSJakub Kruzik ierr = KSPCreate(comm,&def->WtAWinv);CHKERRQ(ierr); 386409a357bSJakub Kruzik ierr = KSPSetOperators(def->WtAWinv,def->WtAW,def->WtAW);CHKERRQ(ierr); 387409a357bSJakub Kruzik ierr = KSPSetType(def->WtAWinv,KSPPREONLY);CHKERRQ(ierr); 388409a357bSJakub Kruzik ierr = KSPGetPC(def->WtAWinv,&pcinner);CHKERRQ(ierr); 389409a357bSJakub Kruzik ierr = PCSetType(pcinner,PCTELESCOPE);CHKERRQ(ierr); 390409a357bSJakub Kruzik /* Reduction factor choice */ 391409a357bSJakub Kruzik red = def->reductionfact; 392409a357bSJakub Kruzik if (red < 0) { 393409a357bSJakub Kruzik ierr = MPI_Comm_size(comm,&commsize);CHKERRQ(ierr); 394409a357bSJakub Kruzik red = ceil((float)commsize/ceil((float)m/commsize)); 395409a357bSJakub Kruzik ierr = PetscObjectTypeCompareAny((PetscObject)(def->WtAW),&match,MATSEQDENSE,MATMPIDENSE,MATDENSE,"");CHKERRQ(ierr); 396409a357bSJakub Kruzik if (match) red = commsize; 397409a357bSJakub Kruzik ierr = PetscInfo1(pc,"Auto choosing reduction factor %D\n",red);CHKERRQ(ierr); /* TODO add level? */ 398409a357bSJakub Kruzik } 399409a357bSJakub Kruzik ierr = PCTelescopeSetReductionFactor(pcinner,red);CHKERRQ(ierr); 400409a357bSJakub Kruzik ierr = PCTelescopeGetKSP(pcinner,&innerksp);CHKERRQ(ierr); 401409a357bSJakub Kruzik ierr = KSPGetPC(innerksp,&pcinner);CHKERRQ(ierr); 402409a357bSJakub Kruzik /* Setup KSP and PC */ 403409a357bSJakub Kruzik if (nextDef) { /* next level for multilevel deflation */ 404409a357bSJakub Kruzik /* set default KSPtype */ 405409a357bSJakub Kruzik if (!def->ksptype) { 406409a357bSJakub Kruzik def->ksptype = KSPFGMRES; 407409a357bSJakub Kruzik if (flgspd) { /* SPD system */ 408409a357bSJakub Kruzik def->ksptype = KSPFCG; 409409a357bSJakub Kruzik } 410409a357bSJakub Kruzik } 411409a357bSJakub Kruzik ierr = KSPSetType(innerksp,def->ksptype);CHKERRQ(ierr); /* TODO iherit from KSP */ 412409a357bSJakub Kruzik ierr = PCSetType(pcinner,PCDEFLATION);CHKERRQ(ierr); /* TODO create coarse preconditinoner M_c = WtMW ? */ 413409a357bSJakub Kruzik ierr = PCDeflationSetSpace(pcinner,nextDef,transp);CHKERRQ(ierr); 4145c0c31c5SJakub Kruzik ierr = PCDeflationSetLvl_Deflation(pcinner,def->nestedlvl+1,def->maxnestedlvl);CHKERRQ(ierr); 415409a357bSJakub Kruzik /* inherit options TODO if not set */ 416409a357bSJakub Kruzik ((PC_Deflation*)(pcinner))->ksptype = def->ksptype; 417409a357bSJakub Kruzik ((PC_Deflation*)(pcinner))->correct = def->correct; 418409a357bSJakub Kruzik ((PC_Deflation*)(pcinner))->adaptiveconv = def->adaptiveconv; 419409a357bSJakub Kruzik ((PC_Deflation*)(pcinner))->adaptiveconst = def->adaptiveconst; 420409a357bSJakub Kruzik ierr = MatDestroy(&nextDef);CHKERRQ(ierr); 421409a357bSJakub Kruzik } else { /* the last level */ 422409a357bSJakub Kruzik ierr = KSPSetType(innerksp,KSPPREONLY);CHKERRQ(ierr); 423409a357bSJakub Kruzik /* TODO Cholesky if flgspd? */ 424409a357bSJakub Kruzik ierr = PCSetType(pc,PCLU);CHKERRQ(ierr); 425409a357bSJakub Kruzik //TODO remove explicit matSolverPackage 426409a357bSJakub Kruzik if (commsize == red) { 427409a357bSJakub Kruzik ierr = PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU);CHKERRQ(ierr); 428409a357bSJakub Kruzik } else { 429409a357bSJakub Kruzik ierr = PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU_DIST);CHKERRQ(ierr); 430409a357bSJakub Kruzik } 431409a357bSJakub Kruzik } 432409a357bSJakub Kruzik /* TODO use def_[lvl]_ if lvl > 0? */ 433409a357bSJakub Kruzik ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr); 434409a357bSJakub Kruzik if (prefix) { 435409a357bSJakub Kruzik ierr = KSPSetOptionsPrefix(innerksp,prefix);CHKERRQ(ierr); 436409a357bSJakub Kruzik ierr = KSPAppendOptionsPrefix(innerksp,"def_");CHKERRQ(ierr); 437409a357bSJakub Kruzik } else { 438409a357bSJakub Kruzik ierr = KSPSetOptionsPrefix(innerksp,"def_");CHKERRQ(ierr); 439409a357bSJakub Kruzik } 440409a357bSJakub Kruzik /* TODO: check if WtAWinv is KSP and move following from this if */ 441409a357bSJakub Kruzik ierr = KSPSetFromOptions(def->WtAWinv);CHKERRQ(ierr); 442409a357bSJakub Kruzik //if (def->adaptiveconv) { 443409a357bSJakub Kruzik // PetscReal *rnorm; 444409a357bSJakub Kruzik // PetscNew(&rnorm); 445409a357bSJakub Kruzik // ierr = KSPSetConvergenceTest(def->WtAWinv,KSPDCGConvergedAdaptive_DCG,rnorm,NULL);CHKERRQ(ierr); 446409a357bSJakub Kruzik //} 447409a357bSJakub Kruzik ierr = KSPSetUp(def->WtAWinv);CHKERRQ(ierr); 448409a357bSJakub Kruzik } 449409a357bSJakub Kruzik 450*f8bf229cSJakub Kruzik /* create work vecs */ 451*f8bf229cSJakub Kruzik ierr = MatCreateVecs(Amat,NULL,&def->work);CHKERRQ(ierr); 452*f8bf229cSJakub Kruzik ierr = KSPCreateVecs(def->WtAWinv,2,&def->workcoarse,0,NULL);CHKERRQ(ierr); 45337eeb815SJakub Kruzik PetscFunctionReturn(0); 45437eeb815SJakub Kruzik } 45537eeb815SJakub Kruzik /* -------------------------------------------------------------------------- */ 45637eeb815SJakub Kruzik static PetscErrorCode PCReset_Deflation(PC pc) 45737eeb815SJakub Kruzik { 45837eeb815SJakub Kruzik PC_Deflation *jac = (PC_Deflation*)pc->data; 45937eeb815SJakub Kruzik PetscErrorCode ierr; 46037eeb815SJakub Kruzik 46137eeb815SJakub Kruzik PetscFunctionBegin; 46237eeb815SJakub Kruzik PetscFunctionReturn(0); 46337eeb815SJakub Kruzik } 46437eeb815SJakub Kruzik 46537eeb815SJakub Kruzik /* 46637eeb815SJakub Kruzik PCDestroy_Deflation - Destroys the private context for the Deflation preconditioner 46737eeb815SJakub Kruzik that was created with PCCreate_Deflation(). 46837eeb815SJakub Kruzik 46937eeb815SJakub Kruzik Input Parameter: 47037eeb815SJakub Kruzik . pc - the preconditioner context 47137eeb815SJakub Kruzik 47237eeb815SJakub Kruzik Application Interface Routine: PCDestroy() 47337eeb815SJakub Kruzik */ 47437eeb815SJakub Kruzik static PetscErrorCode PCDestroy_Deflation(PC pc) 47537eeb815SJakub Kruzik { 47637eeb815SJakub Kruzik PetscErrorCode ierr; 47737eeb815SJakub Kruzik 47837eeb815SJakub Kruzik PetscFunctionBegin; 47937eeb815SJakub Kruzik ierr = PCReset_Deflation(pc);CHKERRQ(ierr); 48037eeb815SJakub Kruzik 48137eeb815SJakub Kruzik /* 48237eeb815SJakub Kruzik Free the private data structure that was hanging off the PC 48337eeb815SJakub Kruzik */ 48437eeb815SJakub Kruzik ierr = PetscFree(pc->data);CHKERRQ(ierr); 48537eeb815SJakub Kruzik PetscFunctionReturn(0); 48637eeb815SJakub Kruzik } 48737eeb815SJakub Kruzik 48837eeb815SJakub Kruzik static PetscErrorCode PCSetFromOptions_Deflation(PetscOptionItems *PetscOptionsObject,PC pc) 48937eeb815SJakub Kruzik { 49037eeb815SJakub Kruzik PC_Deflation *jac = (PC_Deflation*)pc->data; 49137eeb815SJakub Kruzik PetscErrorCode ierr; 49237eeb815SJakub Kruzik PetscBool flg; 49337eeb815SJakub Kruzik PCDeflationType deflt,type; 49437eeb815SJakub Kruzik 49537eeb815SJakub Kruzik PetscFunctionBegin; 49637eeb815SJakub Kruzik ierr = PCDeflationGetType(pc,&deflt);CHKERRQ(ierr); 49737eeb815SJakub Kruzik ierr = PetscOptionsHead(PetscOptionsObject,"Deflation options");CHKERRQ(ierr); 49837eeb815SJakub Kruzik ierr = PetscOptionsEnum("-pc_jacobi_type","How to construct diagonal matrix","PCDeflationSetType",PCDeflationTypes,(PetscEnum)deflt,(PetscEnum*)&type,&flg);CHKERRQ(ierr); 49937eeb815SJakub Kruzik if (flg) { 50037eeb815SJakub Kruzik ierr = PCDeflationSetType(pc,type);CHKERRQ(ierr); 50137eeb815SJakub Kruzik } 50237eeb815SJakub Kruzik ierr = PetscOptionsTail();CHKERRQ(ierr); 50337eeb815SJakub Kruzik PetscFunctionReturn(0); 50437eeb815SJakub Kruzik } 50537eeb815SJakub Kruzik 50637eeb815SJakub Kruzik /*MC 507e361f8b5SJakub Kruzik PCDEFLATION - Deflation preconditioner shifts part of the spectrum to zero (deflates) 508e361f8b5SJakub Kruzik or to a predefined value 50937eeb815SJakub Kruzik 51037eeb815SJakub Kruzik Options Database Key: 511e361f8b5SJakub Kruzik + -pc_deflation_type <init,pre,post> - selects approach to deflation (default: pre) 51237eeb815SJakub Kruzik - -pc_jacobi_abs - use the absolute value of the diagonal entry 51337eeb815SJakub Kruzik 51437eeb815SJakub Kruzik Level: beginner 51537eeb815SJakub Kruzik 51637eeb815SJakub Kruzik Notes: 517e361f8b5SJakub Kruzik todo 51837eeb815SJakub Kruzik 51937eeb815SJakub Kruzik .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 520e662bc50SJakub Kruzik PCDeflationSetType(), PCDeflationSetSpace() 52137eeb815SJakub Kruzik M*/ 52237eeb815SJakub Kruzik 52337eeb815SJakub Kruzik PETSC_EXTERN PetscErrorCode PCCreate_Deflation(PC pc) 52437eeb815SJakub Kruzik { 52537eeb815SJakub Kruzik PC_Deflation *def; 52637eeb815SJakub Kruzik PetscErrorCode ierr; 52737eeb815SJakub Kruzik 52837eeb815SJakub Kruzik PetscFunctionBegin; 52937eeb815SJakub Kruzik ierr = PetscNewLog(pc,&def);CHKERRQ(ierr); 53037eeb815SJakub Kruzik pc->data = (void*)def; 53137eeb815SJakub Kruzik 532e662bc50SJakub Kruzik def->init = PETSC_FALSE; 533e662bc50SJakub Kruzik def->pre = PETSC_TRUE; 534e662bc50SJakub Kruzik def->correct = PETSC_FALSE; 535e662bc50SJakub Kruzik def->truenorm = PETSC_TRUE; 536e662bc50SJakub Kruzik def->reductionfact = -1; 537e662bc50SJakub Kruzik def->spacetype = PC_DEFLATION_SPACE_HAAR; 538e662bc50SJakub Kruzik def->spacesize = 1; 539e662bc50SJakub Kruzik def->extendsp = PETSC_FALSE; 540e662bc50SJakub Kruzik def->nestedlvl = 0; 541e662bc50SJakub Kruzik def->maxnestedlvl = 0; 542e662bc50SJakub Kruzik def->adaptiveconv = PETSC_FALSE; 543e662bc50SJakub Kruzik def->adaptiveconst = 1.0; 54437eeb815SJakub Kruzik 54537eeb815SJakub Kruzik /* 54637eeb815SJakub Kruzik Set the pointers for the functions that are provided above. 54737eeb815SJakub Kruzik Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 54837eeb815SJakub Kruzik are called, they will automatically call these functions. Note we 54937eeb815SJakub Kruzik choose not to provide a couple of these functions since they are 55037eeb815SJakub Kruzik not needed. 55137eeb815SJakub Kruzik */ 55237eeb815SJakub Kruzik pc->ops->apply = PCApply_Deflation; 55337eeb815SJakub Kruzik pc->ops->applytranspose = PCApply_Deflation; 55437eeb815SJakub Kruzik pc->ops->setup = PCSetUp_Deflation; 55537eeb815SJakub Kruzik pc->ops->reset = PCReset_Deflation; 55637eeb815SJakub Kruzik pc->ops->destroy = PCDestroy_Deflation; 55737eeb815SJakub Kruzik pc->ops->setfromoptions = PCSetFromOptions_Deflation; 55837eeb815SJakub Kruzik pc->ops->view = 0; 55937eeb815SJakub Kruzik pc->ops->applyrichardson = 0; 56037eeb815SJakub Kruzik pc->ops->applysymmetricleft = 0; 56137eeb815SJakub Kruzik pc->ops->applysymmetricright = 0; 56237eeb815SJakub Kruzik 56337eeb815SJakub Kruzik ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetType_C",PCDeflationSetType_Deflation);CHKERRQ(ierr); 56437eeb815SJakub Kruzik ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationGetType_C",PCDeflationGetType_Deflation);CHKERRQ(ierr); 565e662bc50SJakub Kruzik ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetSpace_C",PCDeflationSetSpace_Deflation);CHKERRQ(ierr); 5665c0c31c5SJakub Kruzik ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetLvl_C",PCDeflationSetLvl_Deflation);CHKERRQ(ierr); 56737eeb815SJakub Kruzik PetscFunctionReturn(0); 56837eeb815SJakub Kruzik } 56937eeb815SJakub Kruzik 570