Actual source code: redundant.c

      1: #define PETSCKSP_DLL

      3: /*
      4:   This file defines a "solve the problem redundantly on each processor" preconditioner.

      6: */
     7:  #include src/ksp/pc/pcimpl.h
     8:  #include petscksp.h

     10: typedef struct {
     11:   PC         pc;                    /* actual preconditioner used on each processor */
     12:   Vec        x,b;                   /* sequential vectors to hold parallel vectors */
     13:   Mat        *pmats;                /* matrix and optional preconditioner matrix */
     14:   VecScatter scatterin,scatterout;  /* scatter used to move all values to each processor */
     15:   PetscTruth useparallelmat;
     16: } PC_Redundant;

     20: static PetscErrorCode PCView_Redundant(PC pc,PetscViewer viewer)
     21: {
     22:   PC_Redundant   *red = (PC_Redundant*)pc->data;
     24:   PetscMPIInt    rank;
     25:   PetscTruth     iascii,isstring;
     26:   PetscViewer    sviewer;

     29:   MPI_Comm_rank(pc->comm,&rank);
     30:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
     31:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
     32:   if (iascii) {
     33:     PetscViewerASCIIPrintf(viewer,"  Redundant solver preconditioner: Actual PC follows\n");
     34:     PetscViewerGetSingleton(viewer,&sviewer);
     35:     if (!rank) {
     36:       PetscViewerASCIIPushTab(viewer);
     37:       PCView(red->pc,sviewer);
     38:       PetscViewerASCIIPopTab(viewer);
     39:     }
     40:     PetscViewerRestoreSingleton(viewer,&sviewer);
     41:   } else if (isstring) {
     42:     PetscViewerStringSPrintf(viewer," Redundant solver preconditioner");
     43:     PetscViewerGetSingleton(viewer,&sviewer);
     44:     if (!rank) {
     45:       PCView(red->pc,sviewer);
     46:     }
     47:     PetscViewerRestoreSingleton(viewer,&sviewer);
     48:   } else {
     49:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PC redundant",((PetscObject)viewer)->type_name);
     50:   }
     51:   return(0);
     52: }

     56: static PetscErrorCode PCSetUp_Redundant(PC pc)
     57: {
     58:   PC_Redundant   *red  = (PC_Redundant*)pc->data;
     60:   PetscInt       mstart,mlocal,m;
     61:   PetscMPIInt    size;
     62:   IS             isl;
     63:   MatReuse       reuse = MAT_INITIAL_MATRIX;
     64:   MatStructure   str   = DIFFERENT_NONZERO_PATTERN;
     65:   MPI_Comm       comm;
     66:   Vec            vec;

     69:   MatGetVecs(pc->pmat,&vec,0);
     70:   PCSetFromOptions(red->pc);
     71:   VecGetSize(vec,&m);
     72:   if (!pc->setupcalled) {
     73:     VecGetLocalSize(vec,&mlocal);
     74:     VecCreateSeq(PETSC_COMM_SELF,m,&red->x);
     75:     VecDuplicate(red->x,&red->b);
     76:     if (!red->scatterin) {

     78:       /*
     79:          Create the vectors and vector scatter to get the entire vector onto each processor
     80:       */
     81:       VecGetOwnershipRange(vec,&mstart,PETSC_NULL);
     82:       VecScatterCreate(vec,0,red->x,0,&red->scatterin);
     83:       ISCreateStride(pc->comm,mlocal,mstart,1,&isl);
     84:       VecScatterCreate(red->x,isl,vec,isl,&red->scatterout);
     85:       ISDestroy(isl);
     86:     }
     87:   }
     88:   VecDestroy(vec);

     90:   /* if pmatrix set by user is sequential then we do not need to gather the parallel matrix*/

     92:   PetscObjectGetComm((PetscObject)pc->pmat,&comm);
     93:   MPI_Comm_size(comm,&size);
     94:   if (size == 1) {
     95:     red->useparallelmat = PETSC_FALSE;
     96:   }

     98:   if (red->useparallelmat) {
     99:     if (pc->setupcalled == 1 && pc->flag == DIFFERENT_NONZERO_PATTERN) {
    100:       /* destroy old matrices */
    101:       if (red->pmats) {
    102:         MatDestroyMatrices(1,&red->pmats);
    103:       }
    104:     } else if (pc->setupcalled == 1) {
    105:       reuse = MAT_REUSE_MATRIX;
    106:       str   = SAME_NONZERO_PATTERN;
    107:     }
    108:
    109:     /*
    110:        grab the parallel matrix and put it on each processor
    111:     */
    112:     ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);
    113:     MatGetSubMatrices(pc->pmat,1,&isl,&isl,reuse,&red->pmats);
    114:     ISDestroy(isl);

    116:     /* tell sequential PC its operators */
    117:     PCSetOperators(red->pc,red->pmats[0],red->pmats[0],str);
    118:   } else {
    119:     PCSetOperators(red->pc,pc->mat,pc->pmat,pc->flag);
    120:   }
    121:   PCSetFromOptions(red->pc);
    122:   PCSetUp(red->pc);
    123:   return(0);
    124: }

    129: static PetscErrorCode PCApply_Redundant(PC pc,Vec x,Vec y)
    130: {
    131:   PC_Redundant   *red = (PC_Redundant*)pc->data;

    135:   /* move all values to each processor */
    136:   VecScatterBegin(x,red->b,INSERT_VALUES,SCATTER_FORWARD,red->scatterin);
    137:   VecScatterEnd(x,red->b,INSERT_VALUES,SCATTER_FORWARD,red->scatterin);

    139:   /* apply preconditioner on each processor */
    140:   PCApply(red->pc,red->b,red->x);

    142:   /* move local part of values into y vector */
    143:   VecScatterBegin(red->x,y,INSERT_VALUES,SCATTER_FORWARD,red->scatterout);
    144:   VecScatterEnd(red->x,y,INSERT_VALUES,SCATTER_FORWARD,red->scatterout);
    145:   return(0);
    146: }

    151: static PetscErrorCode PCDestroy_Redundant(PC pc)
    152: {
    153:   PC_Redundant   *red = (PC_Redundant*)pc->data;

    157:   if (red->scatterin)  {VecScatterDestroy(red->scatterin);}
    158:   if (red->scatterout) {VecScatterDestroy(red->scatterout);}
    159:   if (red->x)          {VecDestroy(red->x);}
    160:   if (red->b)          {VecDestroy(red->b);}
    161:   if (red->pmats) {
    162:     MatDestroyMatrices(1,&red->pmats);
    163:   }
    164:   PCDestroy(red->pc);
    165:   PetscFree(red);
    166:   return(0);
    167: }

    171: static PetscErrorCode PCSetFromOptions_Redundant(PC pc)
    172: {
    174:   return(0);
    175: }

    180: PetscErrorCode PETSCKSP_DLLEXPORT PCRedundantSetScatter_Redundant(PC pc,VecScatter in,VecScatter out)
    181: {
    182:   PC_Redundant   *red = (PC_Redundant*)pc->data;

    186:   red->scatterin  = in;
    187:   red->scatterout = out;
    188:   PetscObjectReference((PetscObject)in);
    189:   PetscObjectReference((PetscObject)out);
    190:   return(0);
    191: }

    196: /*@
    197:    PCRedundantSetScatter - Sets the scatter used to copy values into the
    198:      redundant local solve and the scatter to move them back into the global
    199:      vector.

    201:    Collective on PC

    203:    Input Parameters:
    204: +  pc - the preconditioner context
    205: .  in - the scatter to move the values in
    206: -  out - the scatter to move them out

    208:    Level: advanced

    210: .keywords: PC, redundant solve
    211: @*/
    212: PetscErrorCode PETSCKSP_DLLEXPORT PCRedundantSetScatter(PC pc,VecScatter in,VecScatter out)
    213: {
    214:   PetscErrorCode ierr,(*f)(PC,VecScatter,VecScatter);

    220:   PetscObjectQueryFunction((PetscObject)pc,"PCRedundantSetScatter_C",(void (**)(void))&f);
    221:   if (f) {
    222:     (*f)(pc,in,out);
    223:   }
    224:   return(0);
    225: }

    230: PetscErrorCode PETSCKSP_DLLEXPORT PCRedundantGetPC_Redundant(PC pc,PC *innerpc)
    231: {
    232:   PC_Redundant *red = (PC_Redundant*)pc->data;

    235:   *innerpc = red->pc;
    236:   return(0);
    237: }

    242: /*@
    243:    PCRedundantGetPC - Gets the sequential PC created by the redundant PC.

    245:    Not Collective

    247:    Input Parameter:
    248: .  pc - the preconditioner context

    250:    Output Parameter:
    251: .  innerpc - the sequential PC

    253:    Level: advanced

    255: .keywords: PC, redundant solve
    256: @*/
    257: PetscErrorCode PETSCKSP_DLLEXPORT PCRedundantGetPC(PC pc,PC *innerpc)
    258: {
    259:   PetscErrorCode ierr,(*f)(PC,PC*);

    264:   PetscObjectQueryFunction((PetscObject)pc,"PCRedundantGetPC_C",(void (**)(void))&f);
    265:   if (f) {
    266:     (*f)(pc,innerpc);
    267:   }
    268:   return(0);
    269: }

    274: PetscErrorCode PETSCKSP_DLLEXPORT PCRedundantGetOperators_Redundant(PC pc,Mat *mat,Mat *pmat)
    275: {
    276:   PC_Redundant *red = (PC_Redundant*)pc->data;

    279:   if (mat)  *mat  = red->pmats[0];
    280:   if (pmat) *pmat = red->pmats[0];
    281:   return(0);
    282: }

    287: /*@
    288:    PCRedundantGetOperators - gets the sequential matrix and preconditioner matrix

    290:    Not Collective

    292:    Input Parameter:
    293: .  pc - the preconditioner context

    295:    Output Parameters:
    296: +  mat - the matrix
    297: -  pmat - the (possibly different) preconditioner matrix

    299:    Level: advanced

    301: .keywords: PC, redundant solve
    302: @*/
    303: PetscErrorCode PETSCKSP_DLLEXPORT PCRedundantGetOperators(PC pc,Mat *mat,Mat *pmat)
    304: {
    305:   PetscErrorCode ierr,(*f)(PC,Mat*,Mat*);

    311:   PetscObjectQueryFunction((PetscObject)pc,"PCRedundantGetOperators_C",(void (**)(void))&f);
    312:   if (f) {
    313:     (*f)(pc,mat,pmat);
    314:   }
    315:   return(0);
    316: }

    318: /* -------------------------------------------------------------------------------------*/
    319: /*MC
    320:      PCREDUNDANT - Runs a preconditioner for the entire problem on each processor

    323:      Options for the redundant preconditioners can be set with -redundant_pc_xxx

    325:    Level: intermediate

    328: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PCRedundantSetScatter(),
    329:            PCRedundantGetPC(), PCRedundantGetOperators()

    331: M*/

    336: PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_Redundant(PC pc)
    337: {
    339:   PC_Redundant   *red;
    340:   const char     *prefix;

    343:   PetscNew(PC_Redundant,&red);
    344:   PetscLogObjectMemory(pc,sizeof(PC_Redundant));
    345:   red->useparallelmat   = PETSC_TRUE;

    347:   /* create the sequential PC that each processor has copy of */
    348:   PCCreate(PETSC_COMM_SELF,&red->pc);
    349:   PCSetType(red->pc,PCLU);
    350:   PCGetOptionsPrefix(pc,&prefix);
    351:   PCSetOptionsPrefix(red->pc,prefix);
    352:   PCAppendOptionsPrefix(red->pc,"redundant_");

    354:   pc->ops->apply             = PCApply_Redundant;
    355:   pc->ops->applytranspose    = 0;
    356:   pc->ops->setup             = PCSetUp_Redundant;
    357:   pc->ops->destroy           = PCDestroy_Redundant;
    358:   pc->ops->setfromoptions    = PCSetFromOptions_Redundant;
    359:   pc->ops->setuponblocks     = 0;
    360:   pc->ops->view              = PCView_Redundant;
    361:   pc->ops->applyrichardson   = 0;

    363:   pc->data              = (void*)red;

    365:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantSetScatter_C","PCRedundantSetScatter_Redundant",
    366:                     PCRedundantSetScatter_Redundant);
    367:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantGetPC_C","PCRedundantGetPC_Redundant",
    368:                     PCRedundantGetPC_Redundant);
    369:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCRedundantGetOperators_C","PCRedundantGetOperators_Redundant",
    370:                     PCRedundantGetOperators_Redundant);

    372:   return(0);
    373: }