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