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