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