1 2 /* 3 Defines a preconditioner defined by R^T S R 4 */ 5 #include <petsc/private/pcimpl.h> 6 #include <petscksp.h> /*I "petscksp.h" I*/ 7 8 typedef struct { 9 KSP ksp; 10 Mat R,P; 11 Vec b,x; 12 } PC_Galerkin; 13 14 static PetscErrorCode PCApply_Galerkin(PC pc,Vec x,Vec y) 15 { 16 PetscErrorCode ierr; 17 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 18 19 PetscFunctionBegin; 20 if (jac->R) { 21 ierr = MatRestrict(jac->R,x,jac->b);CHKERRQ(ierr); 22 } else { 23 ierr = MatRestrict(jac->P,x,jac->b);CHKERRQ(ierr); 24 } 25 ierr = KSPSolve(jac->ksp,jac->b,jac->x);CHKERRQ(ierr); 26 if (jac->P) { 27 ierr = MatInterpolate(jac->P,jac->x,y);CHKERRQ(ierr); 28 } else { 29 ierr = MatInterpolate(jac->R,jac->x,y);CHKERRQ(ierr); 30 } 31 PetscFunctionReturn(0); 32 } 33 34 static PetscErrorCode PCSetUp_Galerkin(PC pc) 35 { 36 PetscErrorCode ierr; 37 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 38 PetscBool a; 39 Vec *xx,*yy; 40 41 PetscFunctionBegin; 42 if (!jac->x) { 43 ierr = KSPGetOperatorsSet(jac->ksp,&a,NULL);CHKERRQ(ierr); 44 if (!a) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set operator of PCGALERKIN KSP with PCGalerkinGetKSP()/KSPSetOperators()"); 45 ierr = KSPCreateVecs(jac->ksp,1,&xx,1,&yy);CHKERRQ(ierr); 46 jac->x = *xx; 47 jac->b = *yy; 48 ierr = PetscFree(xx);CHKERRQ(ierr); 49 ierr = PetscFree(yy);CHKERRQ(ierr); 50 } 51 if (!jac->R && !jac->P) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set restriction or interpolation of PCGALERKIN with PCGalerkinSetRestriction/Interpolation()"); 52 /* should check here that sizes of R/P match size of a */ 53 PetscFunctionReturn(0); 54 } 55 56 static PetscErrorCode PCReset_Galerkin(PC pc) 57 { 58 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 59 PetscErrorCode ierr; 60 61 PetscFunctionBegin; 62 ierr = MatDestroy(&jac->R);CHKERRQ(ierr); 63 ierr = MatDestroy(&jac->P);CHKERRQ(ierr); 64 ierr = VecDestroy(&jac->x);CHKERRQ(ierr); 65 ierr = VecDestroy(&jac->b);CHKERRQ(ierr); 66 ierr = KSPReset(jac->ksp);CHKERRQ(ierr); 67 PetscFunctionReturn(0); 68 } 69 70 static PetscErrorCode PCDestroy_Galerkin(PC pc) 71 { 72 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 73 PetscErrorCode ierr; 74 75 PetscFunctionBegin; 76 ierr = PCReset_Galerkin(pc);CHKERRQ(ierr); 77 ierr = KSPDestroy(&jac->ksp);CHKERRQ(ierr); 78 ierr = PetscFree(pc->data);CHKERRQ(ierr); 79 PetscFunctionReturn(0); 80 } 81 82 static PetscErrorCode PCView_Galerkin(PC pc,PetscViewer viewer) 83 { 84 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 85 PetscErrorCode ierr; 86 PetscBool iascii; 87 88 PetscFunctionBegin; 89 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 90 if (iascii) { 91 ierr = PetscViewerASCIIPrintf(viewer,"Galerkin PC\n");CHKERRQ(ierr); 92 ierr = PetscViewerASCIIPrintf(viewer,"KSP on Galerkin follow\n");CHKERRQ(ierr); 93 ierr = PetscViewerASCIIPrintf(viewer,"---------------------------------\n");CHKERRQ(ierr); 94 } 95 ierr = KSPView(jac->ksp,viewer);CHKERRQ(ierr); 96 PetscFunctionReturn(0); 97 } 98 99 static PetscErrorCode PCGalerkinGetKSP_Galerkin(PC pc,KSP *ksp) 100 { 101 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 102 103 PetscFunctionBegin; 104 *ksp = jac->ksp; 105 PetscFunctionReturn(0); 106 } 107 108 static PetscErrorCode PCGalerkinSetRestriction_Galerkin(PC pc,Mat R) 109 { 110 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 111 PetscErrorCode ierr; 112 113 PetscFunctionBegin; 114 ierr = PetscObjectReference((PetscObject)R);CHKERRQ(ierr); 115 ierr = MatDestroy(&jac->R);CHKERRQ(ierr); 116 jac->R = R; 117 PetscFunctionReturn(0); 118 } 119 120 static PetscErrorCode PCGalerkinSetInterpolation_Galerkin(PC pc,Mat P) 121 { 122 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 123 PetscErrorCode ierr; 124 125 PetscFunctionBegin; 126 ierr = PetscObjectReference((PetscObject)P);CHKERRQ(ierr); 127 ierr = MatDestroy(&jac->P);CHKERRQ(ierr); 128 jac->P = P; 129 PetscFunctionReturn(0); 130 } 131 132 /* -------------------------------------------------------------------------------- */ 133 /*@ 134 PCGalerkinSetRestriction - Sets the restriction operator for the "Galerkin-type" preconditioner 135 136 Logically Collective on PC 137 138 Input Parameter: 139 + pc - the preconditioner context 140 - R - the restriction operator 141 142 Notes: Either this or PCGalerkinSetInterpolation() or both must be called 143 144 Level: Intermediate 145 146 .keywords: PC, set, Galerkin preconditioner 147 148 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN, 149 PCGalerkinSetInterpolation(), PCGalerkinGetKSP() 150 151 @*/ 152 PetscErrorCode PCGalerkinSetRestriction(PC pc,Mat R) 153 { 154 PetscErrorCode ierr; 155 156 PetscFunctionBegin; 157 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 158 ierr = PetscTryMethod(pc,"PCGalerkinSetRestriction_C",(PC,Mat),(pc,R));CHKERRQ(ierr); 159 PetscFunctionReturn(0); 160 } 161 162 /*@ 163 PCGalerkinSetInterpolation - Sets the interpolation operator for the "Galerkin-type" preconditioner 164 165 Logically Collective on PC 166 167 Input Parameter: 168 + pc - the preconditioner context 169 - R - the interpolation operator 170 171 Notes: Either this or PCGalerkinSetRestriction() or both must be called 172 173 Level: Intermediate 174 175 .keywords: PC, set, Galerkin preconditioner 176 177 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN, 178 PCGalerkinSetRestriction(), PCGalerkinGetKSP() 179 180 @*/ 181 PetscErrorCode PCGalerkinSetInterpolation(PC pc,Mat P) 182 { 183 PetscErrorCode ierr; 184 185 PetscFunctionBegin; 186 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 187 ierr = PetscTryMethod(pc,"PCGalerkinSetInterpolation_C",(PC,Mat),(pc,P));CHKERRQ(ierr); 188 PetscFunctionReturn(0); 189 } 190 191 /*@ 192 PCGalerkinGetKSP - Gets the KSP object in the Galerkin PC. 193 194 Not Collective 195 196 Input Parameter: 197 . pc - the preconditioner context 198 199 Output Parameters: 200 . ksp - the KSP object 201 202 Level: Intermediate 203 204 .keywords: PC, get, Galerkin preconditioner, sub preconditioner 205 206 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN, 207 PCGalerkinSetRestriction(), PCGalerkinSetInterpolation() 208 209 @*/ 210 PetscErrorCode PCGalerkinGetKSP(PC pc,KSP *ksp) 211 { 212 PetscErrorCode ierr; 213 214 PetscFunctionBegin; 215 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 216 PetscValidPointer(ksp,2); 217 ierr = PetscUseMethod(pc,"PCGalerkinGetKSP_C",(PC,KSP*),(pc,ksp));CHKERRQ(ierr); 218 PetscFunctionReturn(0); 219 } 220 221 222 /* -------------------------------------------------------------------------------------------*/ 223 224 /*MC 225 PCGALERKIN - Build (part of) a preconditioner by P S R (where P is often R^T) 226 227 $ Use PCGalerkinSetRestriction(pc,R) and/or PCGalerkinSetInterpolation(pc,P) followed by 228 $ PCGalerkinGetKSP(pc,&ksp); KSPSetOperators(ksp,A,....) 229 230 Level: intermediate 231 232 Developer Note: If KSPSetOperators() has not been called on the inner KSP then PCGALERKIN could use MatRARt() or MatPtAP() to compute 233 the operators automatically. 234 Should there be a prefix for the inner KSP. 235 There is no KSPSetFromOptions_Galerkin() that calls KSPSetFromOptions() on the inner KSP 236 237 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 238 PCSHELL, PCKSP, PCGalerkinSetRestriction(), PCGalerkinSetInterpolation(), PCGalerkinGetKSP() 239 240 M*/ 241 242 PETSC_EXTERN PetscErrorCode PCCreate_Galerkin(PC pc) 243 { 244 PetscErrorCode ierr; 245 PC_Galerkin *jac; 246 247 PetscFunctionBegin; 248 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 249 250 pc->ops->apply = PCApply_Galerkin; 251 pc->ops->setup = PCSetUp_Galerkin; 252 pc->ops->reset = PCReset_Galerkin; 253 pc->ops->destroy = PCDestroy_Galerkin; 254 pc->ops->view = PCView_Galerkin; 255 pc->ops->applyrichardson = 0; 256 257 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&jac->ksp);CHKERRQ(ierr); 258 ierr = KSPSetErrorIfNotConverged(jac->ksp,pc->erroriffailure);CHKERRQ(ierr); 259 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 260 261 pc->data = (void*)jac; 262 263 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinSetRestriction_C",PCGalerkinSetRestriction_Galerkin);CHKERRQ(ierr); 264 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinSetInterpolation_C",PCGalerkinSetInterpolation_Galerkin);CHKERRQ(ierr); 265 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinGetKSP_C",PCGalerkinGetKSP_Galerkin);CHKERRQ(ierr); 266 PetscFunctionReturn(0); 267 } 268 269