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