1 #define PETSCKSP_DLL 2 3 /* 4 This file defines a "solve the problem redundantly on each subgroup of processor" preconditioner. 5 */ 6 #include "private/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 xsub,ysub; /* vectors of a subcommunicator to hold parallel vectors of pc->comm */ 12 Vec xdup,ydup; /* parallel vector that congregates xsub or ysub facilitating vector scattering */ 13 Mat pmats; /* matrix and optional preconditioner matrix belong to a subcommunicator */ 14 VecScatter scatterin,scatterout; /* scatter used to move all values to each processor group (subcommunicator) */ 15 PetscTruth useparallelmat; 16 PetscSubcomm *psubcomm; 17 PetscInt nsubcomm; /* num of data structure PetscSubcomm */ 18 } PC_Redundant; 19 20 #undef __FUNCT__ 21 #define __FUNCT__ "PCView_Redundant" 22 static PetscErrorCode PCView_Redundant(PC pc,PetscViewer viewer) 23 { 24 PC_Redundant *red = (PC_Redundant*)pc->data; 25 PetscErrorCode ierr; 26 PetscMPIInt rank; 27 PetscTruth iascii,isstring; 28 PetscViewer sviewer,subviewer; 29 PetscInt color = red->psubcomm->color; 30 31 PetscFunctionBegin; 32 ierr = MPI_Comm_rank(pc->comm,&rank);CHKERRQ(ierr); 33 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 34 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);CHKERRQ(ierr); 35 if (iascii) { 36 ierr = PetscViewerASCIIPrintf(viewer," Redundant solver preconditioner: First PC (color=0) follows\n");CHKERRQ(ierr); 37 ierr = PetscViewerGetSubcomm(viewer,red->pc->comm,&subviewer);CHKERRQ(ierr); 38 if (!color) { /* only view first redundant pc */ 39 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 40 ierr = PCView(red->pc,subviewer);CHKERRQ(ierr); 41 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 42 } 43 ierr = PetscViewerRestoreSubcomm(viewer,red->pc->comm,&subviewer);CHKERRQ(ierr); 44 } else if (isstring) { /* not test it yet! */ 45 ierr = PetscViewerStringSPrintf(viewer," Redundant solver preconditioner");CHKERRQ(ierr); 46 ierr = PetscViewerGetSingleton(viewer,&sviewer);CHKERRQ(ierr); 47 if (!rank) { 48 ierr = PCView(red->pc,sviewer);CHKERRQ(ierr); 49 } 50 ierr = PetscViewerRestoreSingleton(viewer,&sviewer);CHKERRQ(ierr); 51 } else { 52 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PC redundant",((PetscObject)viewer)->type_name); 53 } 54 PetscFunctionReturn(0); 55 } 56 57 #include "include/private/matimpl.h" /*I "petscmat.h" I*/ 58 #undef __FUNCT__ 59 #define __FUNCT__ "PCSetUp_Redundant" 60 static PetscErrorCode PCSetUp_Redundant(PC pc) 61 { 62 PC_Redundant *red = (PC_Redundant*)pc->data; 63 PetscErrorCode ierr; 64 PetscInt mstart,mend,mlocal,m; 65 PetscMPIInt size; 66 MatReuse reuse = MAT_INITIAL_MATRIX; 67 MatStructure str = DIFFERENT_NONZERO_PATTERN; 68 MPI_Comm comm = pc->comm,subcomm; 69 Vec vec; 70 PetscInt mlocal_sub; 71 PetscMPIInt subsize,subrank; 72 PetscInt rstart_sub,rend_sub,mloc_sub; 73 const char *prefix; 74 75 PetscFunctionBegin; 76 ierr = MatGetVecs(pc->pmat,&vec,0);CHKERRQ(ierr); 77 ierr = VecGetSize(vec,&m);CHKERRQ(ierr); 78 79 if (!pc->setupcalled) { 80 ierr = PetscSubcommCreate(comm,red->nsubcomm,&red->psubcomm);CHKERRQ(ierr); 81 ierr = PetscLogObjectMemory(pc,sizeof(PetscSubcomm));CHKERRQ(ierr); 82 83 /* create a new PC that processors in each subcomm have copy of */ 84 subcomm = red->psubcomm->comm; 85 ierr = PCCreate(subcomm,&red->pc);CHKERRQ(ierr); 86 ierr = PCSetType(red->pc,PCLU);CHKERRQ(ierr); 87 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr); 88 ierr = PCSetOptionsPrefix(red->pc,prefix);CHKERRQ(ierr); 89 ierr = PCAppendOptionsPrefix(red->pc,"redundant_");CHKERRQ(ierr); 90 ierr = PCSetFromOptions(red->pc);CHKERRQ(ierr); 91 92 /* create working vectors xsub/ysub and xdup/ydup */ 93 ierr = VecGetLocalSize(vec,&mlocal);CHKERRQ(ierr); 94 ierr = VecGetOwnershipRange(vec,&mstart,&mend);CHKERRQ(ierr); 95 96 /* get local size of xsub/ysub */ 97 ierr = MPI_Comm_size(subcomm,&subsize);CHKERRQ(ierr); 98 ierr = MPI_Comm_rank(subcomm,&subrank);CHKERRQ(ierr); 99 rstart_sub = pc->pmat->rmap.range[red->psubcomm->n*subrank]; /* rstart in xsub/ysub */ 100 if (subrank+1 < subsize){ 101 rend_sub = pc->pmat->rmap.range[red->psubcomm->n*(subrank+1)]; 102 } else { 103 rend_sub = m; 104 } 105 mloc_sub = rend_sub - rstart_sub; 106 ierr = VecCreateMPI(subcomm,mloc_sub,PETSC_DECIDE,&red->ysub);CHKERRQ(ierr); 107 /* create xsub with empty local arrays, because xdup's arrays will be placed into it */ 108 ierr = VecCreateMPIWithArray(subcomm,mloc_sub,PETSC_DECIDE,PETSC_NULL,&red->xsub);CHKERRQ(ierr); 109 110 /* create xdup and ydup. ydup has empty local arrays because ysub's arrays will be place into it. 111 Note: we use communicator dupcomm, not pc->comm! */ 112 ierr = VecCreateMPI(red->psubcomm->dupparent,mloc_sub,PETSC_DECIDE,&red->xdup);CHKERRQ(ierr); 113 ierr = VecCreateMPIWithArray(red->psubcomm->dupparent,mloc_sub,PETSC_DECIDE,PETSC_NULL,&red->ydup);CHKERRQ(ierr); 114 115 /* create vec scatters */ 116 if (!red->scatterin){ 117 IS is1,is2; 118 PetscInt *idx1,*idx2,i,j,k; 119 ierr = PetscMalloc(2*red->psubcomm->n*mlocal*sizeof(PetscInt),&idx1);CHKERRQ(ierr); 120 idx2 = idx1 + red->psubcomm->n*mlocal; 121 j = 0; 122 for (k=0; k<red->psubcomm->n; k++){ 123 for (i=mstart; i<mend; i++){ 124 idx1[j] = i; 125 idx2[j++] = i + m*k; 126 } 127 } 128 ierr = ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx1,&is1);CHKERRQ(ierr); 129 ierr = ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx2,&is2);CHKERRQ(ierr); 130 ierr = VecScatterCreate(vec,is1,red->xdup,is2,&red->scatterin);CHKERRQ(ierr); 131 ierr = ISDestroy(is1);CHKERRQ(ierr); 132 ierr = ISDestroy(is2);CHKERRQ(ierr); 133 134 ierr = ISCreateStride(comm,mlocal,mstart+ red->psubcomm->color*m,1,&is1);CHKERRQ(ierr); 135 ierr = ISCreateStride(comm,mlocal,mstart,1,&is2);CHKERRQ(ierr); 136 ierr = VecScatterCreate(red->xdup,is1,vec,is2,&red->scatterout);CHKERRQ(ierr); 137 ierr = ISDestroy(is1);CHKERRQ(ierr); 138 ierr = ISDestroy(is2);CHKERRQ(ierr); 139 ierr = PetscFree(idx1);CHKERRQ(ierr); 140 } 141 } 142 ierr = VecDestroy(vec);CHKERRQ(ierr); 143 144 /* if pmatrix set by user is sequential then we do not need to gather the parallel matrix */ 145 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 146 if (size == 1) { 147 red->useparallelmat = PETSC_FALSE; 148 } 149 150 if (red->useparallelmat) { 151 if (pc->setupcalled == 1 && pc->flag == DIFFERENT_NONZERO_PATTERN) { 152 /* destroy old matrices */ 153 if (red->pmats) { 154 ierr = MatDestroy(red->pmats);CHKERRQ(ierr); 155 } 156 } else if (pc->setupcalled == 1) { 157 reuse = MAT_REUSE_MATRIX; 158 str = SAME_NONZERO_PATTERN; 159 } 160 161 /* grab the parallel matrix and put it into processors of a subcomminicator */ 162 /*--------------------------------------------------------------------------*/ 163 ierr = VecGetLocalSize(red->ysub,&mlocal_sub);CHKERRQ(ierr); 164 ierr = MatGetRedundantMatrix(pc->pmat,red->psubcomm->n,red->psubcomm->comm,mlocal_sub,reuse,&red->pmats);CHKERRQ(ierr); 165 166 /* tell PC of the subcommunicator its operators */ 167 ierr = PCSetOperators(red->pc,red->pmats,red->pmats,str);CHKERRQ(ierr); 168 } else { 169 ierr = PCSetOperators(red->pc,pc->mat,pc->pmat,pc->flag);CHKERRQ(ierr); 170 } 171 ierr = PCSetFromOptions(red->pc);CHKERRQ(ierr); 172 ierr = PCSetUp(red->pc);CHKERRQ(ierr); 173 PetscFunctionReturn(0); 174 } 175 176 #undef __FUNCT__ 177 #define __FUNCT__ "PCApply_Redundant" 178 static PetscErrorCode PCApply_Redundant(PC pc,Vec x,Vec y) 179 { 180 PC_Redundant *red = (PC_Redundant*)pc->data; 181 PetscErrorCode ierr; 182 PetscScalar *array; 183 184 PetscFunctionBegin; 185 /* scatter x to xdup */ 186 ierr = VecScatterBegin(x,red->xdup,INSERT_VALUES,SCATTER_FORWARD,red->scatterin);CHKERRQ(ierr); 187 ierr = VecScatterEnd(x,red->xdup,INSERT_VALUES,SCATTER_FORWARD,red->scatterin);CHKERRQ(ierr); 188 189 /* place xdup's local array into xsub */ 190 ierr = VecGetArray(red->xdup,&array);CHKERRQ(ierr); 191 ierr = VecPlaceArray(red->xsub,(const PetscScalar*)array);CHKERRQ(ierr); 192 193 /* apply preconditioner on each processor */ 194 ierr = PCApply(red->pc,red->xsub,red->ysub);CHKERRQ(ierr); 195 ierr = VecResetArray(red->xsub);CHKERRQ(ierr); 196 ierr = VecRestoreArray(red->xdup,&array);CHKERRQ(ierr); 197 198 /* place ysub's local array into ydup */ 199 ierr = VecGetArray(red->ysub,&array);CHKERRQ(ierr); 200 ierr = VecPlaceArray(red->ydup,(const PetscScalar*)array);CHKERRQ(ierr); 201 202 /* scatter ydup to y */ 203 ierr = VecScatterBegin(red->ydup,y,INSERT_VALUES,SCATTER_FORWARD,red->scatterout);CHKERRQ(ierr); 204 ierr = VecScatterEnd(red->ydup,y,INSERT_VALUES,SCATTER_FORWARD,red->scatterout);CHKERRQ(ierr); 205 ierr = VecResetArray(red->ydup);CHKERRQ(ierr); 206 ierr = VecRestoreArray(red->ysub,&array);CHKERRQ(ierr); 207 PetscFunctionReturn(0); 208 } 209 210 #undef __FUNCT__ 211 #define __FUNCT__ "PCDestroy_Redundant" 212 static PetscErrorCode PCDestroy_Redundant(PC pc) 213 { 214 PC_Redundant *red = (PC_Redundant*)pc->data; 215 PetscErrorCode ierr; 216 217 PetscFunctionBegin; 218 if (red->scatterin) {ierr = VecScatterDestroy(red->scatterin);CHKERRQ(ierr);} 219 if (red->scatterout) {ierr = VecScatterDestroy(red->scatterout);CHKERRQ(ierr);} 220 if (red->ysub) {ierr = VecDestroy(red->ysub);CHKERRQ(ierr);} 221 if (red->xsub) {ierr = VecDestroy(red->xsub);CHKERRQ(ierr);} 222 if (red->xdup) {ierr = VecDestroy(red->xdup);CHKERRQ(ierr);} 223 if (red->ydup) {ierr = VecDestroy(red->ydup);CHKERRQ(ierr);} 224 if (red->pmats) { 225 ierr = MatDestroy(red->pmats);CHKERRQ(ierr); 226 } 227 if (red->psubcomm) {ierr = PetscSubcommDestroy(red->psubcomm);CHKERRQ(ierr);} 228 if (red->pc) {ierr = PCDestroy(red->pc);CHKERRQ(ierr);} 229 ierr = PetscFree(red);CHKERRQ(ierr); 230 PetscFunctionReturn(0); 231 } 232 233 #undef __FUNCT__ 234 #define __FUNCT__ "PCSetFromOptions_Redundant" 235 static PetscErrorCode PCSetFromOptions_Redundant(PC pc) 236 { 237 PetscErrorCode ierr; 238 PC_Redundant *red = (PC_Redundant*)pc->data; 239 240 PetscFunctionBegin; 241 ierr = PetscOptionsHead("Redundant options");CHKERRQ(ierr); 242 ierr = PetscOptionsInt("-pc_redundant_number_comm","Number of subcommunicators","PCRedundantSetNumComm",red->nsubcomm,&red->nsubcomm,0);CHKERRQ(ierr); 243 ierr = PetscOptionsTail();CHKERRQ(ierr); 244 PetscFunctionReturn(0); 245 } 246 247 EXTERN_C_BEGIN 248 #undef __FUNCT__ 249 #define __FUNCT__ "PCRedundantSetScatter_Redundant" 250 PetscErrorCode PETSCKSP_DLLEXPORT PCRedundantSetScatter_Redundant(PC pc,VecScatter in,VecScatter out) 251 { 252 PC_Redundant *red = (PC_Redundant*)pc->data; 253 PetscErrorCode ierr; 254 255 PetscFunctionBegin; 256 red->scatterin = in; 257 red->scatterout = out; 258 ierr = PetscObjectReference((PetscObject)in);CHKERRQ(ierr); 259 ierr = PetscObjectReference((PetscObject)out);CHKERRQ(ierr); 260 PetscFunctionReturn(0); 261 } 262 EXTERN_C_END 263 264 #undef __FUNCT__ 265 #define __FUNCT__ "PCRedundantSetScatter" 266 /*@ 267 PCRedundantSetScatter - Sets the scatter used to copy values into the 268 redundant local solve and the scatter to move them back into the global 269 vector. 270 271 Collective on PC 272 273 Input Parameters: 274 + pc - the preconditioner context 275 . in - the scatter to move the values in 276 - out - the scatter to move them out 277 278 Level: advanced 279 280 .keywords: PC, redundant solve 281 @*/ 282 PetscErrorCode PETSCKSP_DLLEXPORT PCRedundantSetScatter(PC pc,VecScatter in,VecScatter out) 283 { 284 PetscErrorCode ierr,(*f)(PC,VecScatter,VecScatter); 285 286 PetscFunctionBegin; 287 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 288 PetscValidHeaderSpecific(in,VEC_SCATTER_COOKIE,2); 289 PetscValidHeaderSpecific(out,VEC_SCATTER_COOKIE,3); 290 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCRedundantSetScatter_C",(void (**)(void))&f);CHKERRQ(ierr); 291 if (f) { 292 ierr = (*f)(pc,in,out);CHKERRQ(ierr); 293 } 294 PetscFunctionReturn(0); 295 } 296 297 EXTERN_C_BEGIN 298 #undef __FUNCT__ 299 #define __FUNCT__ "PCRedundantGetPC_Redundant" 300 PetscErrorCode PETSCKSP_DLLEXPORT PCRedundantGetPC_Redundant(PC pc,PC *innerpc) 301 { 302 PC_Redundant *red = (PC_Redundant*)pc->data; 303 304 PetscFunctionBegin; 305 *innerpc = red->pc; 306 PetscFunctionReturn(0); 307 } 308 EXTERN_C_END 309 310 #undef __FUNCT__ 311 #define __FUNCT__ "PCRedundantGetPC" 312 /*@ 313 PCRedundantGetPC - Gets the sequential PC created by the redundant PC. 314 315 Not Collective 316 317 Input Parameter: 318 . pc - the preconditioner context 319 320 Output Parameter: 321 . innerpc - the sequential PC 322 323 Level: advanced 324 325 .keywords: PC, redundant solve 326 @*/ 327 PetscErrorCode PETSCKSP_DLLEXPORT PCRedundantGetPC(PC pc,PC *innerpc) 328 { 329 PetscErrorCode ierr,(*f)(PC,PC*); 330 331 PetscFunctionBegin; 332 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 333 PetscValidPointer(innerpc,2); 334 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCRedundantGetPC_C",(void (**)(void))&f);CHKERRQ(ierr); 335 if (f) { 336 ierr = (*f)(pc,innerpc);CHKERRQ(ierr); 337 } 338 PetscFunctionReturn(0); 339 } 340 341 EXTERN_C_BEGIN 342 #undef __FUNCT__ 343 #define __FUNCT__ "PCRedundantGetOperators_Redundant" 344 PetscErrorCode PETSCKSP_DLLEXPORT PCRedundantGetOperators_Redundant(PC pc,Mat *mat,Mat *pmat) 345 { 346 PC_Redundant *red = (PC_Redundant*)pc->data; 347 348 PetscFunctionBegin; 349 if (mat) *mat = red->pmats; 350 if (pmat) *pmat = red->pmats; 351 PetscFunctionReturn(0); 352 } 353 EXTERN_C_END 354 355 #undef __FUNCT__ 356 #define __FUNCT__ "PCRedundantGetOperators" 357 /*@ 358 PCRedundantGetOperators - gets the sequential matrix and preconditioner matrix 359 360 Not Collective 361 362 Input Parameter: 363 . pc - the preconditioner context 364 365 Output Parameters: 366 + mat - the matrix 367 - pmat - the (possibly different) preconditioner matrix 368 369 Level: advanced 370 371 .keywords: PC, redundant solve 372 @*/ 373 PetscErrorCode PETSCKSP_DLLEXPORT PCRedundantGetOperators(PC pc,Mat *mat,Mat *pmat) 374 { 375 PetscErrorCode ierr,(*f)(PC,Mat*,Mat*); 376 377 PetscFunctionBegin; 378 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 379 if (mat) PetscValidPointer(mat,2); 380 if (pmat) PetscValidPointer(pmat,3); 381 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCRedundantGetOperators_C",(void (**)(void))&f);CHKERRQ(ierr); 382 if (f) { 383 ierr = (*f)(pc,mat,pmat);CHKERRQ(ierr); 384 } 385 PetscFunctionReturn(0); 386 } 387 388 /* -------------------------------------------------------------------------------------*/ 389 /*MC 390 PCREDUNDANT - Runs a preconditioner for the entire problem on subgroups of processors 391 392 Options for the redundant preconditioners can be set with -redundant_pc_xxx 393 394 Options Database: 395 . -pc_redundant_number_comm - number of sub communicators to use 396 397 Level: intermediate 398 399 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCRedundantSetScatter(), 400 PCRedundantGetPC(), PCRedundantGetOperators() 401 M*/ 402 403 EXTERN_C_BEGIN 404 #undef __FUNCT__ 405 #define __FUNCT__ "PCCreate_Redundant" 406 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_Redundant(PC pc) 407 { 408 PetscErrorCode ierr; 409 PC_Redundant *red; 410 PetscMPIInt size; 411 412 PetscFunctionBegin; 413 ierr = PetscNew(PC_Redundant,&red);CHKERRQ(ierr); 414 ierr = PetscLogObjectMemory(pc,sizeof(PC_Redundant));CHKERRQ(ierr); 415 ierr = MPI_Comm_size(pc->comm,&size);CHKERRQ(ierr); 416 red->nsubcomm = size; 417 red->useparallelmat = PETSC_TRUE; 418 pc->data = (void*)red; 419 420 pc->ops->apply = PCApply_Redundant; 421 pc->ops->applytranspose = 0; 422 pc->ops->setup = PCSetUp_Redundant; 423 pc->ops->destroy = PCDestroy_Redundant; 424 pc->ops->setfromoptions = PCSetFromOptions_Redundant; 425 pc->ops->view = PCView_Redundant; 426 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantSetScatter_C","PCRedundantSetScatter_Redundant", 427 PCRedundantSetScatter_Redundant);CHKERRQ(ierr); 428 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantGetPC_C","PCRedundantGetPC_Redundant", 429 PCRedundantGetPC_Redundant);CHKERRQ(ierr); 430 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantGetOperators_C","PCRedundantGetOperators_Redundant", 431 PCRedundantGetOperators_Redundant);CHKERRQ(ierr); 432 PetscFunctionReturn(0); 433 } 434 EXTERN_C_END 435