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