Actual source code: chowiluviennacl.cxx

      2: /*  -------------------------------------------------------------------- */

      4: /*
      5:    Include files needed for the ViennaCL Chow-Patel parallel ILU preconditioner:
      6:      pcimpl.h - private include file intended for use by all preconditioners
      7: */
      8: #define PETSC_SKIP_SPINLOCK
      9: #define PETSC_SKIP_IMMINTRIN_H_CUDAWORKAROUND 1

     11: #include <petsc/private/pcimpl.h>
     12: #include <../src/mat/impls/aij/seq/aij.h>
     13: #include <../src/vec/vec/impls/dvecimpl.h>
     14: #include <../src/mat/impls/aij/seq/seqviennacl/viennaclmatimpl.h>
     15: #include <../src/vec/vec/impls/seq/seqviennacl/viennaclvecimpl.h>
     16: #include <viennacl/linalg/detail/ilu/chow_patel_ilu.hpp>

     18: /*
     19:    Private context (data structure) for the CHOWILUVIENNACL preconditioner.
     20: */
     21: typedef struct {
     22:   viennacl::linalg::chow_patel_ilu_precond< viennacl::compressed_matrix<PetscScalar> > *CHOWILUVIENNACL;
     23: } PC_CHOWILUVIENNACL;

     26: /* -------------------------------------------------------------------------- */
     27: /*
     28:    PCSetUp_CHOWILUVIENNACL - Prepares for the use of the CHOWILUVIENNACL preconditioner
     29:                              by setting data structures and options.

     31:    Input Parameter:
     32: .  pc - the preconditioner context

     34:    Application Interface Routine: PCSetUp()

     36:    Notes:
     37:    The interface routine PCSetUp() is not usually called directly by
     38:    the user, but instead is called by PCApply() if necessary.
     39: */
     40: static PetscErrorCode PCSetUp_CHOWILUVIENNACL(PC pc)
     41: {
     42:   PC_CHOWILUVIENNACL *ilu = (PC_CHOWILUVIENNACL*)pc->data;
     43:   PetscBool           flg = PETSC_FALSE;
     44:   PetscErrorCode      ierr;
     45:   Mat_SeqAIJViennaCL  *gpustruct;

     48:   PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQAIJVIENNACL,&flg);
     49:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Currently only handles ViennaCL matrices");
     50:   if (pc->setupcalled != 0) {
     51:     try {
     52:       delete ilu->CHOWILUVIENNACL;
     53:     } catch(char *ex) {
     54:       SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
     55:     }
     56:   }
     57:   try {
     58: #if defined(PETSC_USE_COMPLEX)
     59:     gpustruct = NULL;
     60:     SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"No support for complex arithmetic in CHOWILUVIENNACL preconditioner");
     61: #else
     62:     MatViennaCLCopyToGPU(pc->pmat);
     63:     gpustruct = (Mat_SeqAIJViennaCL*)(pc->pmat->spptr);

     65:     viennacl::linalg::chow_patel_tag ilu_tag;
     66:     ViennaCLAIJMatrix *mat = (ViennaCLAIJMatrix*)gpustruct->mat;
     67:     ilu->CHOWILUVIENNACL = new viennacl::linalg::chow_patel_ilu_precond<viennacl::compressed_matrix<PetscScalar> >(*mat, ilu_tag);
     68: #endif
     69:   } catch(char *ex) {
     70:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
     71:   }
     72:   return(0);
     73: }

     75: /* -------------------------------------------------------------------------- */
     76: /*
     77:    PCApply_CHOWILUVIENNACL - Applies the CHOWILUVIENNACL preconditioner to a vector.

     79:    Input Parameters:
     80: .  pc - the preconditioner context
     81: .  x - input vector

     83:    Output Parameter:
     84: .  y - output vector

     86:    Application Interface Routine: PCApply()
     87:  */
     88: static PetscErrorCode PCApply_CHOWILUVIENNACL(PC pc,Vec x,Vec y)
     89: {
     90:   PC_CHOWILUVIENNACL            *ilu = (PC_CHOWILUVIENNACL*)pc->data;
     91:   PetscErrorCode                ierr;
     92:   PetscBool                     flg1,flg2;
     93:   viennacl::vector<PetscScalar> const *xarray=NULL;
     94:   viennacl::vector<PetscScalar> *yarray=NULL;

     97:   /*how to apply a certain fixed number of iterations?*/
     98:   PetscObjectTypeCompare((PetscObject)x,VECSEQVIENNACL,&flg1);
     99:   PetscObjectTypeCompare((PetscObject)y,VECSEQVIENNACL,&flg2);
    100:   if (!(flg1 && flg2)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP, "Currently only handles ViennaCL vectors");
    101:   if (!ilu->CHOWILUVIENNACL) {
    102:     PCSetUp_CHOWILUVIENNACL(pc);
    103:   }
    104:   VecSet(y,0.0);
    105:   VecViennaCLGetArrayRead(x,&xarray);
    106:   VecViennaCLGetArrayWrite(y,&yarray);
    107:   try {
    108: #if defined(PETSC_USE_COMPLEX)

    110: #else
    111:     *yarray = *xarray;
    112:     ilu->CHOWILUVIENNACL->apply(*yarray);
    113: #endif
    114:   } catch(char * ex) {
    115:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
    116:   }
    117:   VecViennaCLRestoreArrayRead(x,&xarray);
    118:   VecViennaCLRestoreArrayWrite(y,&yarray);
    119:   PetscObjectStateIncrease((PetscObject)y);
    120:   return(0);
    121: }
    122: /* -------------------------------------------------------------------------- */
    123: /*
    124:    PCDestroy_CHOWILUVIENNACL - Destroys the private context for the CHOWILUVIENNACL preconditioner
    125:    that was created with PCCreate_CHOWILUVIENNACL().

    127:    Input Parameter:
    128: .  pc - the preconditioner context

    130:    Application Interface Routine: PCDestroy()
    131: */
    132: static PetscErrorCode PCDestroy_CHOWILUVIENNACL(PC pc)
    133: {
    134:   PC_CHOWILUVIENNACL  *ilu = (PC_CHOWILUVIENNACL*)pc->data;

    138:   if (ilu->CHOWILUVIENNACL) {
    139:     try {
    140:       delete ilu->CHOWILUVIENNACL;
    141:     } catch(char *ex) {
    142:       SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
    143:     }
    144:   }

    146:   /*
    147:       Free the private data structure that was hanging off the PC
    148:   */
    149:   PetscFree(pc->data);
    150:   return(0);
    151: }

    153: static PetscErrorCode PCSetFromOptions_CHOWILUVIENNACL(PetscOptionItems *PetscOptionsObject,PC pc)
    154: {

    158:   PetscOptionsHead(PetscOptionsObject,"CHOWILUVIENNACL options");
    159:   PetscOptionsTail();
    160:   return(0);
    161: }

    163: /* -------------------------------------------------------------------------- */

    166: /*MC
    167:      PCCHOWILUViennaCL  - A smoothed agglomeration algorithm that can be used via the CUDA, OpenCL, and OpenMP backends of ViennaCL

    169:    Level: advanced

    171: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC

    173: M*/

    175: PETSC_EXTERN PetscErrorCode PCCreate_CHOWILUVIENNACL(PC pc)
    176: {
    177:   PC_CHOWILUVIENNACL  *ilu;

    181:   /*
    182:      Creates the private data structure for this preconditioner and
    183:      attach it to the PC object.
    184:   */
    185:   PetscNewLog(pc,&ilu);
    186:   pc->data = (void*)ilu;

    188:   /*
    189:      Initialize the pointer to zero
    190:      Initialize number of v-cycles to default (1)
    191:   */
    192:   ilu->CHOWILUVIENNACL = 0;

    194:   /*
    195:       Set the pointers for the functions that are provided above.
    196:       Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
    197:       are called, they will automatically call these functions.  Note we
    198:       choose not to provide a couple of these functions since they are
    199:       not needed.
    200:   */
    201:   pc->ops->apply               = PCApply_CHOWILUVIENNACL;
    202:   pc->ops->applytranspose      = 0;
    203:   pc->ops->setup               = PCSetUp_CHOWILUVIENNACL;
    204:   pc->ops->destroy             = PCDestroy_CHOWILUVIENNACL;
    205:   pc->ops->setfromoptions      = PCSetFromOptions_CHOWILUVIENNACL;
    206:   pc->ops->view                = 0;
    207:   pc->ops->applyrichardson     = 0;
    208:   pc->ops->applysymmetricleft  = 0;
    209:   pc->ops->applysymmetricright = 0;
    210:   return(0);
    211: }