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); 213 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCRedundantSetScatter_C",(void (**)(void))&f);CHKERRQ(ierr); 214 if (f) { 215 ierr = (*f)(pc,in,out);CHKERRQ(ierr); 216 } 217 PetscFunctionReturn(0); 218 } 219 220 EXTERN_C_BEGIN 221 #undef __FUNCT__ 222 #define __FUNCT__ "PCRedundantGetPC_Redundant" 223 int PCRedundantGetPC_Redundant(PC pc,PC *innerpc) 224 { 225 PC_Redundant *red = (PC_Redundant*)pc->data; 226 227 PetscFunctionBegin; 228 *innerpc = red->pc; 229 PetscFunctionReturn(0); 230 } 231 EXTERN_C_END 232 233 #undef __FUNCT__ 234 #define __FUNCT__ "PCRedundantGetPC" 235 /*@ 236 PCRedundantGetPC - Gets the sequential PC created by the redundant PC. 237 238 Not Collective 239 240 Input Parameter: 241 . pc - the preconditioner context 242 243 Output Parameter: 244 . innerpc - the sequential PC 245 246 Level: advanced 247 248 .keywords: PC, redundant solve 249 @*/ 250 int PCRedundantGetPC(PC pc,PC *innerpc) 251 { 252 int ierr,(*f)(PC,PC*); 253 254 PetscFunctionBegin; 255 PetscValidHeaderSpecific(pc,PC_COOKIE); 256 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCRedundantGetPC_C",(void (**)(void))&f);CHKERRQ(ierr); 257 if (f) { 258 ierr = (*f)(pc,innerpc);CHKERRQ(ierr); 259 } 260 PetscFunctionReturn(0); 261 } 262 263 EXTERN_C_BEGIN 264 #undef __FUNCT__ 265 #define __FUNCT__ "PCRedundantGetOperators_Redundant" 266 int PCRedundantGetOperators_Redundant(PC pc,Mat *mat,Mat *pmat) 267 { 268 PC_Redundant *red = (PC_Redundant*)pc->data; 269 270 PetscFunctionBegin; 271 if (mat) *mat = red->pmats[0]; 272 if (pmat) *pmat = red->pmats[0]; 273 PetscFunctionReturn(0); 274 } 275 EXTERN_C_END 276 277 #undef __FUNCT__ 278 #define __FUNCT__ "PCRedundantGetOperators" 279 /*@ 280 PCRedundantGetOperators - gets the sequential matrix and preconditioner matrix 281 282 Not Collective 283 284 Input Parameter: 285 . pc - the preconditioner context 286 287 Output Parameters: 288 + mat - the matrix 289 - pmat - the (possibly different) preconditioner matrix 290 291 Level: advanced 292 293 .keywords: PC, redundant solve 294 @*/ 295 int PCRedundantGetOperators(PC pc,Mat *mat,Mat *pmat) 296 { 297 int ierr,(*f)(PC,Mat*,Mat*); 298 299 PetscFunctionBegin; 300 PetscValidHeaderSpecific(pc,PC_COOKIE); 301 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCRedundantGetOperators_C",(void (**)(void))&f);CHKERRQ(ierr); 302 if (f) { 303 ierr = (*f)(pc,mat,pmat);CHKERRQ(ierr); 304 } 305 PetscFunctionReturn(0); 306 } 307 308 /* -------------------------------------------------------------------------------------*/ 309 EXTERN_C_BEGIN 310 #undef __FUNCT__ 311 #define __FUNCT__ "PCCreate_Redundant" 312 int PCCreate_Redundant(PC pc) 313 { 314 int ierr; 315 PC_Redundant *red; 316 char *prefix; 317 318 PetscFunctionBegin; 319 ierr = PetscNew(PC_Redundant,&red);CHKERRQ(ierr); 320 PetscLogObjectMemory(pc,sizeof(PC_Redundant)); 321 ierr = PetscMemzero(red,sizeof(PC_Redundant));CHKERRQ(ierr); 322 red->useparallelmat = PETSC_TRUE; 323 324 /* create the sequential PC that each processor has copy of */ 325 ierr = PCCreate(PETSC_COMM_SELF,&red->pc);CHKERRQ(ierr); 326 ierr = PCSetType(red->pc,PCLU);CHKERRQ(ierr); 327 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr); 328 ierr = PCSetOptionsPrefix(red->pc,prefix);CHKERRQ(ierr); 329 ierr = PCAppendOptionsPrefix(red->pc,"redundant_");CHKERRQ(ierr); 330 331 pc->ops->apply = PCApply_Redundant; 332 pc->ops->applytranspose = 0; 333 pc->ops->setup = PCSetUp_Redundant; 334 pc->ops->destroy = PCDestroy_Redundant; 335 pc->ops->setfromoptions = PCSetFromOptions_Redundant; 336 pc->ops->setuponblocks = 0; 337 pc->ops->view = PCView_Redundant; 338 pc->ops->applyrichardson = 0; 339 340 pc->data = (void*)red; 341 342 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantSetScatter_C","PCRedundantSetScatter_Redundant", 343 PCRedundantSetScatter_Redundant);CHKERRQ(ierr); 344 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantGetPC_C","PCRedundantGetPC_Redundant", 345 PCRedundantGetPC_Redundant);CHKERRQ(ierr); 346 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantGetOperators_C","PCRedundantGetOperators_Redundant", 347 PCRedundantGetOperators_Redundant);CHKERRQ(ierr); 348 349 PetscFunctionReturn(0); 350 } 351 EXTERN_C_END 352