1 /* 2 This file defines a "solve the problem redundantly on each processor" preconditioner. 3 4 */ 5 #include "src/ksp/pc/pcimpl.h" /*I "petscpc.h" I*/ 6 #include "petscksp.h" 7 8 typedef struct { 9 PC pc; /* actual preconditioner used on each processor */ 10 Vec x,b; /* sequential vectors to hold parallel vectors */ 11 Mat *pmats; /* matrix and optional preconditioner matrix */ 12 VecScatter scatterin,scatterout; /* scatter used to move all values to each processor */ 13 PetscTruth useparallelmat; 14 } PC_Redundant; 15 16 #undef __FUNCT__ 17 #define __FUNCT__ "PCView_Redundant" 18 static PetscErrorCode PCView_Redundant(PC pc,PetscViewer viewer) 19 { 20 PC_Redundant *red = (PC_Redundant*)pc->data; 21 PetscErrorCode ierr; 22 int rank; 23 PetscTruth iascii,isstring; 24 PetscViewer sviewer; 25 26 PetscFunctionBegin; 27 ierr = MPI_Comm_rank(pc->comm,&rank);CHKERRQ(ierr); 28 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 29 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);CHKERRQ(ierr); 30 if (iascii) { 31 ierr = PetscViewerASCIIPrintf(viewer," Redundant solver preconditioner: Actual PC follows\n");CHKERRQ(ierr); 32 ierr = PetscViewerGetSingleton(viewer,&sviewer);CHKERRQ(ierr); 33 if (!rank) { 34 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 35 ierr = PCView(red->pc,sviewer);CHKERRQ(ierr); 36 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 37 } 38 ierr = PetscViewerRestoreSingleton(viewer,&sviewer);CHKERRQ(ierr); 39 } else if (isstring) { 40 ierr = PetscViewerStringSPrintf(viewer," Redundant solver preconditioner");CHKERRQ(ierr); 41 ierr = PetscViewerGetSingleton(viewer,&sviewer);CHKERRQ(ierr); 42 if (!rank) { 43 ierr = PCView(red->pc,sviewer);CHKERRQ(ierr); 44 } 45 ierr = PetscViewerRestoreSingleton(viewer,&sviewer);CHKERRQ(ierr); 46 } else { 47 SETERRQ1(1,"Viewer type %s not supported for PC redundant",((PetscObject)viewer)->type_name); 48 } 49 PetscFunctionReturn(0); 50 } 51 52 #undef __FUNCT__ 53 #define __FUNCT__ "PCSetUp_Redundant" 54 static PetscErrorCode PCSetUp_Redundant(PC pc) 55 { 56 PC_Redundant *red = (PC_Redundant*)pc->data; 57 PetscErrorCode ierr; 58 int mstart,mlocal,m,size; 59 IS isl; 60 MatReuse reuse = MAT_INITIAL_MATRIX; 61 MatStructure str = DIFFERENT_NONZERO_PATTERN; 62 MPI_Comm comm; 63 Vec vec; 64 65 PetscFunctionBegin; 66 ierr = MatGetVecs(pc->pmat,&vec,0);CHKERRQ(ierr); 67 ierr = PCSetFromOptions(red->pc);CHKERRQ(ierr); 68 ierr = VecGetSize(vec,&m);CHKERRQ(ierr); 69 if (!pc->setupcalled) { 70 ierr = VecGetLocalSize(vec,&mlocal);CHKERRQ(ierr); 71 ierr = VecCreateSeq(PETSC_COMM_SELF,m,&red->x);CHKERRQ(ierr); 72 ierr = VecDuplicate(red->x,&red->b);CHKERRQ(ierr); 73 if (!red->scatterin) { 74 75 /* 76 Create the vectors and vector scatter to get the entire vector onto each processor 77 */ 78 ierr = VecGetOwnershipRange(vec,&mstart,PETSC_NULL);CHKERRQ(ierr); 79 ierr = VecScatterCreate(vec,0,red->x,0,&red->scatterin);CHKERRQ(ierr); 80 ierr = ISCreateStride(pc->comm,mlocal,mstart,1,&isl);CHKERRQ(ierr); 81 ierr = VecScatterCreate(red->x,isl,vec,isl,&red->scatterout);CHKERRQ(ierr); 82 ierr = ISDestroy(isl);CHKERRQ(ierr); 83 } 84 } 85 ierr = VecDestroy(vec);CHKERRQ(ierr); 86 87 /* if pmatrix set by user is sequential then we do not need to gather the parallel matrix*/ 88 89 ierr = PetscObjectGetComm((PetscObject)pc->pmat,&comm);CHKERRQ(ierr); 90 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);CHKERRQ(ierr); 91 if (size == 1) { 92 red->useparallelmat = PETSC_FALSE; 93 } 94 95 if (red->useparallelmat) { 96 if (pc->setupcalled == 1 && pc->flag == DIFFERENT_NONZERO_PATTERN) { 97 /* destroy old matrices */ 98 if (red->pmats) { 99 ierr = MatDestroyMatrices(1,&red->pmats);CHKERRQ(ierr); 100 } 101 } else if (pc->setupcalled == 1) { 102 reuse = MAT_REUSE_MATRIX; 103 str = SAME_NONZERO_PATTERN; 104 } 105 106 /* 107 grab the parallel matrix and put it on each processor 108 */ 109 ierr = ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);CHKERRQ(ierr); 110 ierr = MatGetSubMatrices(pc->pmat,1,&isl,&isl,reuse,&red->pmats);CHKERRQ(ierr); 111 ierr = ISDestroy(isl);CHKERRQ(ierr); 112 113 /* tell sequential PC its operators */ 114 ierr = PCSetOperators(red->pc,red->pmats[0],red->pmats[0],str);CHKERRQ(ierr); 115 } else { 116 ierr = PCSetOperators(red->pc,pc->mat,pc->pmat,pc->flag);CHKERRQ(ierr); 117 } 118 ierr = PCSetFromOptions(red->pc);CHKERRQ(ierr); 119 ierr = PCSetUp(red->pc);CHKERRQ(ierr); 120 PetscFunctionReturn(0); 121 } 122 123 124 #undef __FUNCT__ 125 #define __FUNCT__ "PCApply_Redundant" 126 static PetscErrorCode PCApply_Redundant(PC pc,Vec x,Vec y) 127 { 128 PC_Redundant *red = (PC_Redundant*)pc->data; 129 PetscErrorCode ierr; 130 131 PetscFunctionBegin; 132 /* move all values to each processor */ 133 ierr = VecScatterBegin(x,red->b,INSERT_VALUES,SCATTER_FORWARD,red->scatterin);CHKERRQ(ierr); 134 ierr = VecScatterEnd(x,red->b,INSERT_VALUES,SCATTER_FORWARD,red->scatterin);CHKERRQ(ierr); 135 136 /* apply preconditioner on each processor */ 137 ierr = PCApply(red->pc,red->b,red->x);CHKERRQ(ierr); 138 139 /* move local part of values into y vector */ 140 ierr = VecScatterBegin(red->x,y,INSERT_VALUES,SCATTER_FORWARD,red->scatterout);CHKERRQ(ierr); 141 ierr = VecScatterEnd(red->x,y,INSERT_VALUES,SCATTER_FORWARD,red->scatterout);CHKERRQ(ierr); 142 PetscFunctionReturn(0); 143 } 144 145 146 #undef __FUNCT__ 147 #define __FUNCT__ "PCDestroy_Redundant" 148 static PetscErrorCode PCDestroy_Redundant(PC pc) 149 { 150 PC_Redundant *red = (PC_Redundant*)pc->data; 151 PetscErrorCode ierr; 152 153 PetscFunctionBegin; 154 if (red->scatterin) {ierr = VecScatterDestroy(red->scatterin);CHKERRQ(ierr);} 155 if (red->scatterout) {ierr = VecScatterDestroy(red->scatterout);CHKERRQ(ierr);} 156 if (red->x) {ierr = VecDestroy(red->x);CHKERRQ(ierr);} 157 if (red->b) {ierr = VecDestroy(red->b);CHKERRQ(ierr);} 158 if (red->pmats) { 159 ierr = MatDestroyMatrices(1,&red->pmats);CHKERRQ(ierr); 160 } 161 ierr = PCDestroy(red->pc);CHKERRQ(ierr); 162 ierr = PetscFree(red);CHKERRQ(ierr); 163 PetscFunctionReturn(0); 164 } 165 166 #undef __FUNCT__ 167 #define __FUNCT__ "PCSetFromOptions_Redundant" 168 static PetscErrorCode PCSetFromOptions_Redundant(PC pc) 169 { 170 PetscFunctionBegin; 171 PetscFunctionReturn(0); 172 } 173 174 EXTERN_C_BEGIN 175 #undef __FUNCT__ 176 #define __FUNCT__ "PCRedundantSetScatter_Redundant" 177 PetscErrorCode PCRedundantSetScatter_Redundant(PC pc,VecScatter in,VecScatter out) 178 { 179 PC_Redundant *red = (PC_Redundant*)pc->data; 180 PetscErrorCode ierr; 181 182 PetscFunctionBegin; 183 red->scatterin = in; 184 red->scatterout = out; 185 ierr = PetscObjectReference((PetscObject)in);CHKERRQ(ierr); 186 ierr = PetscObjectReference((PetscObject)out);CHKERRQ(ierr); 187 PetscFunctionReturn(0); 188 } 189 EXTERN_C_END 190 191 #undef __FUNCT__ 192 #define __FUNCT__ "PCRedundantSetScatter" 193 /*@ 194 PCRedundantSetScatter - Sets the scatter used to copy values into the 195 redundant local solve and the scatter to move them back into the global 196 vector. 197 198 Collective on PC 199 200 Input Parameters: 201 + pc - the preconditioner context 202 . in - the scatter to move the values in 203 - out - the scatter to move them out 204 205 Level: advanced 206 207 .keywords: PC, redundant solve 208 @*/ 209 PetscErrorCode PCRedundantSetScatter(PC pc,VecScatter in,VecScatter out) 210 { 211 PetscErrorCode ierr,(*f)(PC,VecScatter,VecScatter); 212 213 PetscFunctionBegin; 214 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 215 PetscValidHeaderSpecific(in,VEC_SCATTER_COOKIE,2); 216 PetscValidHeaderSpecific(out,VEC_SCATTER_COOKIE,3); 217 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCRedundantSetScatter_C",(void (**)(void))&f);CHKERRQ(ierr); 218 if (f) { 219 ierr = (*f)(pc,in,out);CHKERRQ(ierr); 220 } 221 PetscFunctionReturn(0); 222 } 223 224 EXTERN_C_BEGIN 225 #undef __FUNCT__ 226 #define __FUNCT__ "PCRedundantGetPC_Redundant" 227 PetscErrorCode PCRedundantGetPC_Redundant(PC pc,PC *innerpc) 228 { 229 PC_Redundant *red = (PC_Redundant*)pc->data; 230 231 PetscFunctionBegin; 232 *innerpc = red->pc; 233 PetscFunctionReturn(0); 234 } 235 EXTERN_C_END 236 237 #undef __FUNCT__ 238 #define __FUNCT__ "PCRedundantGetPC" 239 /*@ 240 PCRedundantGetPC - Gets the sequential PC created by the redundant PC. 241 242 Not Collective 243 244 Input Parameter: 245 . pc - the preconditioner context 246 247 Output Parameter: 248 . innerpc - the sequential PC 249 250 Level: advanced 251 252 .keywords: PC, redundant solve 253 @*/ 254 PetscErrorCode PCRedundantGetPC(PC pc,PC *innerpc) 255 { 256 PetscErrorCode ierr,(*f)(PC,PC*); 257 258 PetscFunctionBegin; 259 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 260 PetscValidPointer(innerpc,2); 261 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCRedundantGetPC_C",(void (**)(void))&f);CHKERRQ(ierr); 262 if (f) { 263 ierr = (*f)(pc,innerpc);CHKERRQ(ierr); 264 } 265 PetscFunctionReturn(0); 266 } 267 268 EXTERN_C_BEGIN 269 #undef __FUNCT__ 270 #define __FUNCT__ "PCRedundantGetOperators_Redundant" 271 PetscErrorCode PCRedundantGetOperators_Redundant(PC pc,Mat *mat,Mat *pmat) 272 { 273 PC_Redundant *red = (PC_Redundant*)pc->data; 274 275 PetscFunctionBegin; 276 if (mat) *mat = red->pmats[0]; 277 if (pmat) *pmat = red->pmats[0]; 278 PetscFunctionReturn(0); 279 } 280 EXTERN_C_END 281 282 #undef __FUNCT__ 283 #define __FUNCT__ "PCRedundantGetOperators" 284 /*@ 285 PCRedundantGetOperators - gets the sequential matrix and preconditioner matrix 286 287 Not Collective 288 289 Input Parameter: 290 . pc - the preconditioner context 291 292 Output Parameters: 293 + mat - the matrix 294 - pmat - the (possibly different) preconditioner matrix 295 296 Level: advanced 297 298 .keywords: PC, redundant solve 299 @*/ 300 PetscErrorCode PCRedundantGetOperators(PC pc,Mat *mat,Mat *pmat) 301 { 302 PetscErrorCode ierr,(*f)(PC,Mat*,Mat*); 303 304 PetscFunctionBegin; 305 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 306 if (mat) PetscValidPointer(mat,2); 307 if (pmat) PetscValidPointer(pmat,3); 308 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCRedundantGetOperators_C",(void (**)(void))&f);CHKERRQ(ierr); 309 if (f) { 310 ierr = (*f)(pc,mat,pmat);CHKERRQ(ierr); 311 } 312 PetscFunctionReturn(0); 313 } 314 315 /* -------------------------------------------------------------------------------------*/ 316 /*MC 317 PCREDUNDANT - Runs a preconditioner for the entire problem on each processor 318 319 320 Options for the redundant preconditioners can be set with -redundant_pc_xxx 321 322 Level: intermediate 323 324 325 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCRedundantSetScatter(), 326 PCRedundantGetPC(), PCRedundantGetOperators() 327 328 M*/ 329 330 EXTERN_C_BEGIN 331 #undef __FUNCT__ 332 #define __FUNCT__ "PCCreate_Redundant" 333 PetscErrorCode PCCreate_Redundant(PC pc) 334 { 335 PetscErrorCode ierr; 336 PC_Redundant *red; 337 char *prefix; 338 339 PetscFunctionBegin; 340 ierr = PetscNew(PC_Redundant,&red);CHKERRQ(ierr); 341 PetscLogObjectMemory(pc,sizeof(PC_Redundant)); 342 ierr = PetscMemzero(red,sizeof(PC_Redundant));CHKERRQ(ierr); 343 red->useparallelmat = PETSC_TRUE; 344 345 /* create the sequential PC that each processor has copy of */ 346 ierr = PCCreate(PETSC_COMM_SELF,&red->pc);CHKERRQ(ierr); 347 ierr = PCSetType(red->pc,PCLU);CHKERRQ(ierr); 348 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr); 349 ierr = PCSetOptionsPrefix(red->pc,prefix);CHKERRQ(ierr); 350 ierr = PCAppendOptionsPrefix(red->pc,"redundant_");CHKERRQ(ierr); 351 352 pc->ops->apply = PCApply_Redundant; 353 pc->ops->applytranspose = 0; 354 pc->ops->setup = PCSetUp_Redundant; 355 pc->ops->destroy = PCDestroy_Redundant; 356 pc->ops->setfromoptions = PCSetFromOptions_Redundant; 357 pc->ops->setuponblocks = 0; 358 pc->ops->view = PCView_Redundant; 359 pc->ops->applyrichardson = 0; 360 361 pc->data = (void*)red; 362 363 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantSetScatter_C","PCRedundantSetScatter_Redundant", 364 PCRedundantSetScatter_Redundant);CHKERRQ(ierr); 365 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantGetPC_C","PCRedundantGetPC_Redundant", 366 PCRedundantGetPC_Redundant);CHKERRQ(ierr); 367 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantGetOperators_C","PCRedundantGetOperators_Redundant", 368 PCRedundantGetOperators_Redundant);CHKERRQ(ierr); 369 370 PetscFunctionReturn(0); 371 } 372 EXTERN_C_END 373