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; 69 Vec vec; 70 PetscInt mlocal_sub; 71 PetscMPIInt subsize,subrank; 72 PetscInt rstart_sub,rend_sub,mloc_sub,nsubcomm; 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 MPI_Comm 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 ierr = PetscSubcommDestroy(red->psubcomm);CHKERRQ(ierr); 228 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 PetscMPIInt size; 240 241 PetscFunctionBegin; 242 ierr = PetscOptionsHead("Redundant options");CHKERRQ(ierr); 243 ierr = PetscOptionsInt("-pc_redundant_number_comm","Number of subcommunicators","PCRedundantSetNumComm",red->nsubcomm,&red->nsubcomm,0);CHKERRQ(ierr); 244 ierr = PetscOptionsTail();CHKERRQ(ierr); 245 PetscFunctionReturn(0); 246 } 247 248 EXTERN_C_BEGIN 249 #undef __FUNCT__ 250 #define __FUNCT__ "PCRedundantSetScatter_Redundant" 251 PetscErrorCode PETSCKSP_DLLEXPORT PCRedundantSetScatter_Redundant(PC pc,VecScatter in,VecScatter out) 252 { 253 PC_Redundant *red = (PC_Redundant*)pc->data; 254 PetscErrorCode ierr; 255 256 PetscFunctionBegin; 257 red->scatterin = in; 258 red->scatterout = out; 259 ierr = PetscObjectReference((PetscObject)in);CHKERRQ(ierr); 260 ierr = PetscObjectReference((PetscObject)out);CHKERRQ(ierr); 261 PetscFunctionReturn(0); 262 } 263 EXTERN_C_END 264 265 #undef __FUNCT__ 266 #define __FUNCT__ "PCRedundantSetScatter" 267 /*@ 268 PCRedundantSetScatter - Sets the scatter used to copy values into the 269 redundant local solve and the scatter to move them back into the global 270 vector. 271 272 Collective on PC 273 274 Input Parameters: 275 + pc - the preconditioner context 276 . in - the scatter to move the values in 277 - out - the scatter to move them out 278 279 Level: advanced 280 281 .keywords: PC, redundant solve 282 @*/ 283 PetscErrorCode PETSCKSP_DLLEXPORT PCRedundantSetScatter(PC pc,VecScatter in,VecScatter out) 284 { 285 PetscErrorCode ierr,(*f)(PC,VecScatter,VecScatter); 286 287 PetscFunctionBegin; 288 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 289 PetscValidHeaderSpecific(in,VEC_SCATTER_COOKIE,2); 290 PetscValidHeaderSpecific(out,VEC_SCATTER_COOKIE,3); 291 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCRedundantSetScatter_C",(void (**)(void))&f);CHKERRQ(ierr); 292 if (f) { 293 ierr = (*f)(pc,in,out);CHKERRQ(ierr); 294 } 295 PetscFunctionReturn(0); 296 } 297 298 EXTERN_C_BEGIN 299 #undef __FUNCT__ 300 #define __FUNCT__ "PCRedundantGetPC_Redundant" 301 PetscErrorCode PETSCKSP_DLLEXPORT PCRedundantGetPC_Redundant(PC pc,PC *innerpc) 302 { 303 PC_Redundant *red = (PC_Redundant*)pc->data; 304 305 PetscFunctionBegin; 306 *innerpc = red->pc; 307 PetscFunctionReturn(0); 308 } 309 EXTERN_C_END 310 311 #undef __FUNCT__ 312 #define __FUNCT__ "PCRedundantGetPC" 313 /*@ 314 PCRedundantGetPC - Gets the sequential PC created by the redundant PC. 315 316 Not Collective 317 318 Input Parameter: 319 . pc - the preconditioner context 320 321 Output Parameter: 322 . innerpc - the sequential PC 323 324 Level: advanced 325 326 .keywords: PC, redundant solve 327 @*/ 328 PetscErrorCode PETSCKSP_DLLEXPORT PCRedundantGetPC(PC pc,PC *innerpc) 329 { 330 PetscErrorCode ierr,(*f)(PC,PC*); 331 332 PetscFunctionBegin; 333 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 334 PetscValidPointer(innerpc,2); 335 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCRedundantGetPC_C",(void (**)(void))&f);CHKERRQ(ierr); 336 if (f) { 337 ierr = (*f)(pc,innerpc);CHKERRQ(ierr); 338 } 339 PetscFunctionReturn(0); 340 } 341 342 EXTERN_C_BEGIN 343 #undef __FUNCT__ 344 #define __FUNCT__ "PCRedundantGetOperators_Redundant" 345 PetscErrorCode PETSCKSP_DLLEXPORT PCRedundantGetOperators_Redundant(PC pc,Mat *mat,Mat *pmat) 346 { 347 PC_Redundant *red = (PC_Redundant*)pc->data; 348 349 PetscFunctionBegin; 350 if (mat) *mat = red->pmats; 351 if (pmat) *pmat = red->pmats; 352 PetscFunctionReturn(0); 353 } 354 EXTERN_C_END 355 356 #undef __FUNCT__ 357 #define __FUNCT__ "PCRedundantGetOperators" 358 /*@ 359 PCRedundantGetOperators - gets the sequential matrix and preconditioner matrix 360 361 Not Collective 362 363 Input Parameter: 364 . pc - the preconditioner context 365 366 Output Parameters: 367 + mat - the matrix 368 - pmat - the (possibly different) preconditioner matrix 369 370 Level: advanced 371 372 .keywords: PC, redundant solve 373 @*/ 374 PetscErrorCode PETSCKSP_DLLEXPORT PCRedundantGetOperators(PC pc,Mat *mat,Mat *pmat) 375 { 376 PetscErrorCode ierr,(*f)(PC,Mat*,Mat*); 377 378 PetscFunctionBegin; 379 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 380 if (mat) PetscValidPointer(mat,2); 381 if (pmat) PetscValidPointer(pmat,3); 382 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCRedundantGetOperators_C",(void (**)(void))&f);CHKERRQ(ierr); 383 if (f) { 384 ierr = (*f)(pc,mat,pmat);CHKERRQ(ierr); 385 } 386 PetscFunctionReturn(0); 387 } 388 389 /* -------------------------------------------------------------------------------------*/ 390 /*MC 391 PCREDUNDANT - Runs a preconditioner for the entire problem on subgroups of processors 392 393 Options for the redundant preconditioners can be set with -redundant_pc_xxx 394 395 Options Database: 396 . -pc_redundant_number_comm - number of sub communicators to use 397 398 Level: intermediate 399 400 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCRedundantSetScatter(), 401 PCRedundantGetPC(), PCRedundantGetOperators() 402 M*/ 403 404 EXTERN_C_BEGIN 405 #undef __FUNCT__ 406 #define __FUNCT__ "PCCreate_Redundant" 407 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_Redundant(PC pc) 408 { 409 PetscErrorCode ierr; 410 PC_Redundant *red; 411 PetscMPIInt size; 412 413 PetscFunctionBegin; 414 ierr = PetscNew(PC_Redundant,&red);CHKERRQ(ierr); 415 ierr = PetscLogObjectMemory(pc,sizeof(PC_Redundant));CHKERRQ(ierr); 416 ierr = MPI_Comm_size(pc->comm,&size);CHKERRQ(ierr); 417 red->nsubcomm = size; 418 red->useparallelmat = PETSC_TRUE; 419 pc->data = (void*)red; 420 421 pc->ops->apply = PCApply_Redundant; 422 pc->ops->applytranspose = 0; 423 pc->ops->setup = PCSetUp_Redundant; 424 pc->ops->destroy = PCDestroy_Redundant; 425 pc->ops->setfromoptions = PCSetFromOptions_Redundant; 426 pc->ops->view = PCView_Redundant; 427 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantSetScatter_C","PCRedundantSetScatter_Redundant", 428 PCRedundantSetScatter_Redundant);CHKERRQ(ierr); 429 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantGetPC_C","PCRedundantGetPC_Redundant", 430 PCRedundantGetPC_Redundant);CHKERRQ(ierr); 431 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantGetOperators_C","PCRedundantGetOperators_Redundant", 432 PCRedundantGetOperators_Redundant);CHKERRQ(ierr); 433 PetscFunctionReturn(0); 434 } 435 EXTERN_C_END 436