Actual source code: redundant.c

petsc-3.3-p7 2013-05-11
  2: /*
  3:   This file defines a "solve the problem redundantly on each subgroup of processor" preconditioner.
  4: */
  5: #include <petsc-private/pcimpl.h>     /*I "petscpc.h" I*/
  6: #include <petscksp.h>           /*I "petscksp.h" I*/

  8: typedef struct {
  9:   KSP          ksp;
 10:   PC           pc;                   /* actual preconditioner used on each processor */
 11:   Vec          xsub,ysub;            /* vectors of a subcommunicator to hold parallel vectors of ((PetscObject)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:   PetscBool    useparallelmat;
 16:   PetscSubcomm psubcomm;
 17:   PetscInt     nsubcomm;           /* num of data structure PetscSubcomm */
 18: } PC_Redundant;

 22: static PetscErrorCode PCView_Redundant(PC pc,PetscViewer viewer)
 23: {
 24:   PC_Redundant   *red = (PC_Redundant*)pc->data;
 26:   PetscBool      iascii,isstring;
 27:   PetscViewer    subviewer;

 30:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
 31:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
 32:   if (iascii) {
 33:     if (!red->psubcomm) {
 34:       PetscViewerASCIIPrintf(viewer,"  Redundant preconditioner: Not yet setup\n");
 35:     } else {
 36:       PetscViewerASCIIPrintf(viewer,"  Redundant preconditioner: First (color=0) of %D PCs follows\n",red->nsubcomm);
 37:       PetscViewerGetSubcomm(viewer,((PetscObject)red->pc)->comm,&subviewer);
 38:       if (!red->psubcomm->color) { /* only view first redundant pc */
 39:         PetscViewerASCIIPushTab(viewer);
 40:         KSPView(red->ksp,subviewer);
 41:         PetscViewerASCIIPopTab(viewer);
 42:       }
 43:       PetscViewerRestoreSubcomm(viewer,((PetscObject)red->pc)->comm,&subviewer);
 44:     }
 45:   } else if (isstring) {
 46:     PetscViewerStringSPrintf(viewer," Redundant solver preconditioner");
 47:   } else {
 48:     SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for PC redundant",((PetscObject)viewer)->type_name);
 49:   }
 50:   return(0);
 51: }

 55: static PetscErrorCode PCSetUp_Redundant(PC pc)
 56: {
 57:   PC_Redundant   *red = (PC_Redundant*)pc->data;
 59:   PetscInt       mstart,mend,mlocal,m,mlocal_sub,rstart_sub,rend_sub,mloc_sub;
 60:   PetscMPIInt    size;
 61:   MatReuse       reuse = MAT_INITIAL_MATRIX;
 62:   MatStructure   str = DIFFERENT_NONZERO_PATTERN;
 63:   MPI_Comm       comm = ((PetscObject)pc)->comm,subcomm;
 64:   Vec            vec;
 65:   PetscMPIInt    subsize,subrank;
 66:   const char     *prefix;
 67:   const PetscInt *range;

 70:   MatGetVecs(pc->pmat,&vec,0);
 71:   VecGetSize(vec,&m);

 73:   if (!pc->setupcalled) {
 74:     if (!red->psubcomm) {
 75:       PetscSubcommCreate(comm,&red->psubcomm);
 76:       PetscSubcommSetNumber(red->psubcomm,red->nsubcomm);
 77:       PetscSubcommSetType(red->psubcomm,PETSC_SUBCOMM_INTERLACED);
 78:       PetscLogObjectMemory(pc,sizeof(PetscSubcomm));

 80:       /* create a new PC that processors in each subcomm have copy of */
 81:       subcomm = red->psubcomm->comm;
 82:       KSPCreate(subcomm,&red->ksp);
 83:       PetscObjectIncrementTabLevel((PetscObject)red->ksp,(PetscObject)pc,1);
 84:       PetscLogObjectParent(pc,red->ksp);
 85:       KSPSetType(red->ksp,KSPPREONLY);
 86:       KSPGetPC(red->ksp,&red->pc);
 87:       PCSetType(red->pc,PCLU);

 89:       PCGetOptionsPrefix(pc,&prefix);
 90:       KSPSetOptionsPrefix(red->ksp,prefix);
 91:       KSPAppendOptionsPrefix(red->ksp,"redundant_");
 92:     } else {
 93:        subcomm = red->psubcomm->comm;
 94:     }

 96:     /* create working vectors xsub/ysub and xdup/ydup */
 97:     VecGetLocalSize(vec,&mlocal);
 98:     VecGetOwnershipRange(vec,&mstart,&mend);

100:     /* get local size of xsub/ysub */
101:     MPI_Comm_size(subcomm,&subsize);
102:     MPI_Comm_rank(subcomm,&subrank);
103:     MatGetOwnershipRanges(pc->pmat,&range);
104:     rstart_sub = range[red->psubcomm->n*subrank]; /* rstart in xsub/ysub */
105:     if (subrank+1 < subsize){
106:       rend_sub = range[red->psubcomm->n*(subrank+1)];
107:     } else {
108:       rend_sub = m;
109:     }
110:     mloc_sub = rend_sub - rstart_sub;
111:     VecCreateMPI(subcomm,mloc_sub,PETSC_DECIDE,&red->ysub);
112:     /* create xsub with empty local arrays, because xdup's arrays will be placed into it */
113:     VecCreateMPIWithArray(subcomm,1,mloc_sub,PETSC_DECIDE,PETSC_NULL,&red->xsub);

115:     /* create xdup and ydup. ydup has empty local arrays because ysub's arrays will be place into it. 
116:        Note: we use communicator dupcomm, not ((PetscObject)pc)->comm! */
117:     VecCreateMPI(red->psubcomm->dupparent,mloc_sub,PETSC_DECIDE,&red->xdup);
118:     VecCreateMPIWithArray(red->psubcomm->dupparent,1,mloc_sub,PETSC_DECIDE,PETSC_NULL,&red->ydup);
119: 
120:     /* create vec scatters */
121:     if (!red->scatterin){
122:       IS       is1,is2;
123:       PetscInt *idx1,*idx2,i,j,k;

125:       PetscMalloc2(red->psubcomm->n*mlocal,PetscInt,&idx1,red->psubcomm->n*mlocal,PetscInt,&idx2);
126:       j = 0;
127:       for (k=0; k<red->psubcomm->n; k++){
128:         for (i=mstart; i<mend; i++){
129:           idx1[j]   = i;
130:           idx2[j++] = i + m*k;
131:         }
132:       }
133:       ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx1,PETSC_COPY_VALUES,&is1);
134:       ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx2,PETSC_COPY_VALUES,&is2);
135:       VecScatterCreate(vec,is1,red->xdup,is2,&red->scatterin);
136:       ISDestroy(&is1);
137:       ISDestroy(&is2);

139:       ISCreateStride(comm,mlocal,mstart+ red->psubcomm->color*m,1,&is1);
140:       ISCreateStride(comm,mlocal,mstart,1,&is2);
141:       VecScatterCreate(red->xdup,is1,vec,is2,&red->scatterout);
142:       ISDestroy(&is1);
143:       ISDestroy(&is2);
144:       PetscFree2(idx1,idx2);
145:     }
146:   }
147:   VecDestroy(&vec);

149:   /* if pmatrix set by user is sequential then we do not need to gather the parallel matrix */
150:   MPI_Comm_size(comm,&size);
151:   if (size == 1) {
152:     red->useparallelmat = PETSC_FALSE;
153:   }

155:   if (red->useparallelmat) {
156:     if (pc->setupcalled == 1 && pc->flag == DIFFERENT_NONZERO_PATTERN) {
157:       /* destroy old matrices */
158:       MatDestroy(&red->pmats);
159:     } else if (pc->setupcalled == 1) {
160:       reuse = MAT_REUSE_MATRIX;
161:       str   = SAME_NONZERO_PATTERN;
162:     }
163: 
164:     /* grab the parallel matrix and put it into processors of a subcomminicator */
165:     /*--------------------------------------------------------------------------*/
166:     VecGetLocalSize(red->ysub,&mlocal_sub);
167:     MatGetRedundantMatrix(pc->pmat,red->psubcomm->n,red->psubcomm->comm,mlocal_sub,reuse,&red->pmats);
168:     /* tell PC of the subcommunicator its operators */
169:     KSPSetOperators(red->ksp,red->pmats,red->pmats,str);
170:   } else {
171:     KSPSetOperators(red->ksp,pc->mat,pc->pmat,pc->flag);
172:   }
173:   if (pc->setfromoptionscalled){
174:     KSPSetFromOptions(red->ksp);
175:   }
176:   KSPSetUp(red->ksp);
177:   return(0);
178: }

182: static PetscErrorCode PCApply_Redundant(PC pc,Vec x,Vec y)
183: {
184:   PC_Redundant   *red = (PC_Redundant*)pc->data;
186:   PetscScalar    *array;

189:   /* scatter x to xdup */
190:   VecScatterBegin(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);
191:   VecScatterEnd(red->scatterin,x,red->xdup,INSERT_VALUES,SCATTER_FORWARD);
192: 
193:   /* place xdup's local array into xsub */
194:   VecGetArray(red->xdup,&array);
195:   VecPlaceArray(red->xsub,(const PetscScalar*)array);

197:   /* apply preconditioner on each processor */
198:   KSPSolve(red->ksp,red->xsub,red->ysub);
199:   VecResetArray(red->xsub);
200:   VecRestoreArray(red->xdup,&array);
201: 
202:   /* place ysub's local array into ydup */
203:   VecGetArray(red->ysub,&array);
204:   VecPlaceArray(red->ydup,(const PetscScalar*)array);

206:   /* scatter ydup to y */
207:   VecScatterBegin(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);
208:   VecScatterEnd(red->scatterout,red->ydup,y,INSERT_VALUES,SCATTER_FORWARD);
209:   VecResetArray(red->ydup);
210:   VecRestoreArray(red->ysub,&array);
211:   return(0);
212: }

216: static PetscErrorCode PCReset_Redundant(PC pc)
217: {
218:   PC_Redundant   *red = (PC_Redundant*)pc->data;

222:   VecScatterDestroy(&red->scatterin);
223:   VecScatterDestroy(&red->scatterout);
224:   VecDestroy(&red->ysub);
225:   VecDestroy(&red->xsub);
226:   VecDestroy(&red->xdup);
227:   VecDestroy(&red->ydup);
228:   MatDestroy(&red->pmats);
229:   if (red->ksp) {KSPReset(red->ksp);}
230:   return(0);
231: }

235: static PetscErrorCode PCDestroy_Redundant(PC pc)
236: {
237:   PC_Redundant   *red = (PC_Redundant*)pc->data;

241:   PCReset_Redundant(pc);
242:   KSPDestroy(&red->ksp);
243:   PetscSubcommDestroy(&red->psubcomm);
244:   PetscFree(pc->data);
245:   return(0);
246: }

250: static PetscErrorCode PCSetFromOptions_Redundant(PC pc)
251: {
253:   PC_Redundant   *red = (PC_Redundant*)pc->data;

256:   PetscOptionsHead("Redundant options");
257:   PetscOptionsInt("-pc_redundant_number","Number of redundant pc","PCRedundantSetNumber",red->nsubcomm,&red->nsubcomm,0);
258:   PetscOptionsTail();
259:   return(0);
260: }

262: EXTERN_C_BEGIN
265: PetscErrorCode  PCRedundantSetNumber_Redundant(PC pc,PetscInt nreds)
266: {
267:   PC_Redundant *red = (PC_Redundant*)pc->data;

270:   red->nsubcomm = nreds;
271:   return(0);
272: }
273: EXTERN_C_END

277: /*@
278:    PCRedundantSetNumber - Sets the number of redundant preconditioner contexts.

280:    Logically Collective on PC

282:    Input Parameters:
283: +  pc - the preconditioner context
284: -  nredundant - number of redundant preconditioner contexts; for example if you are using 64 MPI processes and
285:                               use an nredundant of 4 there will be 4 parallel solves each on 16 = 64/4 processes.

287:    Level: advanced

289: .keywords: PC, redundant solve
290: @*/
291: PetscErrorCode  PCRedundantSetNumber(PC pc,PetscInt nredundant)
292: {

297:   if (nredundant <= 0) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "num of redundant pc %D must be positive",nredundant);
298:   PetscTryMethod(pc,"PCRedundantSetNumber_C",(PC,PetscInt),(pc,nredundant));
299:   return(0);
300: }

302: EXTERN_C_BEGIN
305: PetscErrorCode  PCRedundantSetScatter_Redundant(PC pc,VecScatter in,VecScatter out)
306: {
307:   PC_Redundant   *red = (PC_Redundant*)pc->data;

311:   PetscObjectReference((PetscObject)in);
312:   VecScatterDestroy(&red->scatterin);
313:   red->scatterin  = in;
314:   PetscObjectReference((PetscObject)out);
315:   VecScatterDestroy(&red->scatterout);
316:   red->scatterout = out;
317:   return(0);
318: }
319: EXTERN_C_END

323: /*@
324:    PCRedundantSetScatter - Sets the scatter used to copy values into the
325:      redundant local solve and the scatter to move them back into the global
326:      vector.

328:    Logically Collective on PC

330:    Input Parameters:
331: +  pc - the preconditioner context
332: .  in - the scatter to move the values in
333: -  out - the scatter to move them out

335:    Level: advanced

337: .keywords: PC, redundant solve
338: @*/
339: PetscErrorCode  PCRedundantSetScatter(PC pc,VecScatter in,VecScatter out)
340: {

347:   PetscTryMethod(pc,"PCRedundantSetScatter_C",(PC,VecScatter,VecScatter),(pc,in,out));
348:   return(0);
349: }

351: EXTERN_C_BEGIN
354: PetscErrorCode  PCRedundantGetKSP_Redundant(PC pc,KSP *innerksp)
355: {
357:   PC_Redundant   *red = (PC_Redundant*)pc->data;
358:   MPI_Comm       comm,subcomm;
359:   const char     *prefix;

362:   if (!red->psubcomm) {
363:     PetscObjectGetComm((PetscObject)pc,&comm);
364:     PetscSubcommCreate(comm,&red->psubcomm);
365:     PetscSubcommSetNumber(red->psubcomm,red->nsubcomm);
366:     PetscSubcommSetType(red->psubcomm,PETSC_SUBCOMM_INTERLACED);
367:     PetscLogObjectMemory(pc,sizeof(PetscSubcomm));

369:     /* create a new PC that processors in each subcomm have copy of */
370:     subcomm = red->psubcomm->comm;
371:     KSPCreate(subcomm,&red->ksp);
372:     PetscObjectIncrementTabLevel((PetscObject)red->ksp,(PetscObject)pc,1);
373:     PetscLogObjectParent(pc,red->ksp);
374:     KSPSetType(red->ksp,KSPPREONLY);
375:     KSPGetPC(red->ksp,&red->pc);
376:     PCSetType(red->pc,PCLU);

378:     PCGetOptionsPrefix(pc,&prefix);
379:     KSPSetOptionsPrefix(red->ksp,prefix);
380:     KSPAppendOptionsPrefix(red->ksp,"redundant_");
381:   }
382:   *innerksp = red->ksp;
383:   return(0);
384: }
385: EXTERN_C_END

389: /*@
390:    PCRedundantGetKSP - Gets the less parallel KSP created by the redundant PC.

392:    Not Collective

394:    Input Parameter:
395: .  pc - the preconditioner context

397:    Output Parameter:
398: .  innerksp - the KSP on the smaller set of processes

400:    Level: advanced

402: .keywords: PC, redundant solve
403: @*/
404: PetscErrorCode  PCRedundantGetKSP(PC pc,KSP *innerksp)
405: {

411:   PetscTryMethod(pc,"PCRedundantGetKSP_C",(PC,KSP*),(pc,innerksp));
412:   return(0);
413: }

415: EXTERN_C_BEGIN
418: PetscErrorCode  PCRedundantGetOperators_Redundant(PC pc,Mat *mat,Mat *pmat)
419: {
420:   PC_Redundant *red = (PC_Redundant*)pc->data;

423:   if (mat)  *mat  = red->pmats;
424:   if (pmat) *pmat = red->pmats;
425:   return(0);
426: }
427: EXTERN_C_END

431: /*@
432:    PCRedundantGetOperators - gets the sequential matrix and preconditioner matrix

434:    Not Collective

436:    Input Parameter:
437: .  pc - the preconditioner context

439:    Output Parameters:
440: +  mat - the matrix
441: -  pmat - the (possibly different) preconditioner matrix

443:    Level: advanced

445: .keywords: PC, redundant solve
446: @*/
447: PetscErrorCode  PCRedundantGetOperators(PC pc,Mat *mat,Mat *pmat)
448: {

455:   PetscTryMethod(pc,"PCRedundantGetOperators_C",(PC,Mat*,Mat*),(pc,mat,pmat));
456:   return(0);
457: }

459: /* -------------------------------------------------------------------------------------*/
460: /*MC
461:      PCREDUNDANT - Runs a KSP solver with preconditioner for the entire problem on subgroups of processors

463:      Options for the redundant preconditioners can be set with -redundant_pc_xxx for the redundant KSP with -redundant_ksp_xxx

465:   Options Database:
466: .  -pc_redundant_number <n> - number of redundant solves, for example if you are using 64 MPI processes and
467:                               use an n of 4 there will be 4 parallel solves each on 16 = 64/4 processes.

469:    Level: intermediate

471:    Notes: The default KSP is preonly and the default PC is LU.

473:    Developer Notes: Note that PCSetInitialGuessNonzero()  is not used by this class but likely should be.

475: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PCRedundantSetScatter(),
476:            PCRedundantGetKSP(), PCRedundantGetOperators(), PCRedundantSetNumber()
477: M*/

479: EXTERN_C_BEGIN
482: PetscErrorCode  PCCreate_Redundant(PC pc)
483: {
485:   PC_Redundant   *red;
486:   PetscMPIInt    size;
487: 
489:   PetscNewLog(pc,PC_Redundant,&red);
490:   MPI_Comm_size(((PetscObject)pc)->comm,&size);
491:   red->nsubcomm       = size;
492:   red->useparallelmat = PETSC_TRUE;
493:   pc->data            = (void*)red;

495:   pc->ops->apply           = PCApply_Redundant;
496:   pc->ops->applytranspose  = 0;
497:   pc->ops->setup           = PCSetUp_Redundant;
498:   pc->ops->destroy         = PCDestroy_Redundant;
499:   pc->ops->reset           = PCReset_Redundant;
500:   pc->ops->setfromoptions  = PCSetFromOptions_Redundant;
501:   pc->ops->view            = PCView_Redundant;
502:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantSetScatter_C","PCRedundantSetScatter_Redundant",
503:                     PCRedundantSetScatter_Redundant);
504:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantSetNumber_C","PCRedundantSetNumber_Redundant",
505:                     PCRedundantSetNumber_Redundant);
506:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantGetKSP_C","PCRedundantGetKSP_Redundant",
507:                     PCRedundantGetKSP_Redundant);
508:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantGetOperators_C","PCRedundantGetOperators_Redundant",
509:                     PCRedundantGetOperators_Redundant);
510:   return(0);
511: }
512: EXTERN_C_END