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