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