Actual source code: saviennacl.cxx

      1: /*
      2:    Include files needed for the ViennaCL Smoothed Aggregation preconditioner:
      3:      pcimpl.h - private include file intended for use by all preconditioners
      4: */
      5: #define PETSC_SKIP_IMMINTRIN_H_CUDAWORKAROUND 1
      6: #include <petsc/private/pcimpl.h>
      7: #include <../src/mat/impls/aij/seq/aij.h>
      8: #include <../src/vec/vec/impls/dvecimpl.h>
      9: #include <../src/mat/impls/aij/seq/seqviennacl/viennaclmatimpl.h>
     10: #include <../src/vec/vec/impls/seq/seqviennacl/viennaclvecimpl.h>
     11: #include <viennacl/linalg/amg.hpp>

     13: /*
     14:    Private context (data structure) for the SAVIENNACL preconditioner.
     15: */
     16: typedef struct {
     17:   viennacl::linalg::amg_precond<viennacl::compressed_matrix<PetscScalar>> *SAVIENNACL;
     18: } PC_SAVIENNACL;

     20: /*
     21:    PCSetUp_SAVIENNACL - Prepares for the use of the SAVIENNACL preconditioner
     22:                         by setting data structures and options.

     24:    Input Parameter:
     25: .  pc - the preconditioner context

     27:    Application Interface Routine: PCSetUp()

     29:    Note:
     30:    The interface routine PCSetUp() is not usually called directly by
     31:    the user, but instead is called by PCApply() if necessary.
     32: */
     33: static PetscErrorCode PCSetUp_SAVIENNACL(PC pc)
     34: {
     35:   PC_SAVIENNACL      *sa  = (PC_SAVIENNACL *)pc->data;
     36:   PetscBool           flg = PETSC_FALSE;
     37:   Mat_SeqAIJViennaCL *gpustruct;

     39:   PetscFunctionBegin;
     40:   PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATSEQAIJVIENNACL, &flg));
     41:   PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Currently only handles ViennaCL matrices");
     42:   if (pc->setupcalled) {
     43:     try {
     44:       delete sa->SAVIENNACL;
     45:     } catch (char *ex) {
     46:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex);
     47:     }
     48:   }
     49:   try {
     50: #if defined(PETSC_USE_COMPLEX)
     51:     gpustruct = NULL;
     52:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "No support for complex arithmetic in SAVIENNACL preconditioner");
     53: #else
     54:     PetscCall(MatViennaCLCopyToGPU(pc->pmat));
     55:     gpustruct = (Mat_SeqAIJViennaCL *)pc->pmat->spptr;

     57:     viennacl::linalg::amg_tag amg_tag_sa_pmis;
     58:     amg_tag_sa_pmis.set_coarsening_method(viennacl::linalg::AMG_COARSENING_METHOD_MIS2_AGGREGATION);
     59:     amg_tag_sa_pmis.set_interpolation_method(viennacl::linalg::AMG_INTERPOLATION_METHOD_SMOOTHED_AGGREGATION);
     60:     ViennaCLAIJMatrix *mat = (ViennaCLAIJMatrix *)gpustruct->mat;
     61:     sa->SAVIENNACL         = new viennacl::linalg::amg_precond<viennacl::compressed_matrix<PetscScalar>>(*mat, amg_tag_sa_pmis);
     62:     sa->SAVIENNACL->setup();
     63: #endif
     64:   } catch (char *ex) {
     65:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex);
     66:   }
     67:   PetscFunctionReturn(PETSC_SUCCESS);
     68: }

     70: /*
     71:    PCApply_SAVIENNACL - Applies the SAVIENNACL preconditioner to a vector.

     73:    Input Parameters:
     74: .  pc - the preconditioner context
     75: .  x - input vector

     77:    Output Parameter:
     78: .  y - output vector

     80:    Application Interface Routine: PCApply()
     81:  */
     82: static PetscErrorCode PCApply_SAVIENNACL(PC pc, Vec x, Vec y)
     83: {
     84:   PC_SAVIENNACL                       *sac = (PC_SAVIENNACL *)pc->data;
     85:   PetscBool                            flg1, flg2;
     86:   viennacl::vector<PetscScalar> const *xarray = NULL;
     87:   viennacl::vector<PetscScalar>       *yarray = NULL;

     89:   PetscFunctionBegin;
     90:   /*how to apply a certain fixed number of iterations?*/
     91:   PetscCall(PetscObjectTypeCompare((PetscObject)x, VECSEQVIENNACL, &flg1));
     92:   PetscCall(PetscObjectTypeCompare((PetscObject)y, VECSEQVIENNACL, &flg2));
     93:   PetscCheck(flg1 && flg2, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Currently only handles ViennaCL vectors");
     94:   if (!sac->SAVIENNACL) PetscCall(PCSetUp_SAVIENNACL(pc));
     95:   PetscCall(VecViennaCLGetArrayRead(x, &xarray));
     96:   PetscCall(VecViennaCLGetArrayWrite(y, &yarray));
     97:   try {
     98: #if !defined(PETSC_USE_COMPLEX)
     99:     *yarray = *xarray;
    100:     sac->SAVIENNACL->apply(*yarray);
    101: #endif
    102:   } catch (char *ex) {
    103:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex);
    104:   }
    105:   PetscCall(VecViennaCLRestoreArrayRead(x, &xarray));
    106:   PetscCall(VecViennaCLRestoreArrayWrite(y, &yarray));
    107:   PetscCall(PetscObjectStateIncrease((PetscObject)y));
    108:   PetscFunctionReturn(PETSC_SUCCESS);
    109: }

    111: /*
    112:    PCDestroy_SAVIENNACL - Destroys the private context for the SAVIENNACL preconditioner
    113:    that was created with PCCreate_SAVIENNACL().

    115:    Input Parameter:
    116: .  pc - the preconditioner context

    118:    Application Interface Routine: PCDestroy()
    119: */
    120: static PetscErrorCode PCDestroy_SAVIENNACL(PC pc)
    121: {
    122:   PC_SAVIENNACL *sac = (PC_SAVIENNACL *)pc->data;

    124:   PetscFunctionBegin;
    125:   if (sac->SAVIENNACL) {
    126:     try {
    127:       delete sac->SAVIENNACL;
    128:     } catch (char *ex) {
    129:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex);
    130:     }
    131:   }

    133:   /*
    134:       Free the private data structure that was hanging off the PC
    135:   */
    136:   PetscCall(PetscFree(pc->data));
    137:   PetscFunctionReturn(PETSC_SUCCESS);
    138: }

    140: static PetscErrorCode PCSetFromOptions_SAVIENNACL(PC pc, PetscOptionItems PetscOptionsObject)
    141: {
    142:   PetscFunctionBegin;
    143:   PetscOptionsHeadBegin(PetscOptionsObject, "SAVIENNACL options");
    144:   PetscOptionsHeadEnd();
    145:   PetscFunctionReturn(PETSC_SUCCESS);
    146: }

    148: /*MC
    149:   PCSAVIENNACL  - A smoothed agglomeration algorithm that can be used via the CUDA, OpenCL, and OpenMP backends of ViennaCL

    151:   Level: advanced

    153:   Developer Notes:
    154:   This `PCType` does not appear to be registered

    156: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`
    157: M*/

    159: PETSC_EXTERN PetscErrorCode PCCreate_SAVIENNACL(PC pc)
    160: {
    161:   PC_SAVIENNACL *sac;

    163:   PetscFunctionBegin;
    164:   /*
    165:      Creates the private data structure for this preconditioner and
    166:      attach it to the PC object.
    167:   */
    168:   PetscCall(PetscNew(&sac));
    169:   pc->data = (void *)sac;

    171:   /*
    172:      Initialize the pointer to zero
    173:      Initialize number of v-cycles to default (1)
    174:   */
    175:   sac->SAVIENNACL = 0;

    177:   /*
    178:       Set the pointers for the functions that are provided above.
    179:       Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
    180:       are called, they will automatically call these functions.  Note we
    181:       choose not to provide a couple of these functions since they are
    182:       not needed.
    183:   */
    184:   pc->ops->apply               = PCApply_SAVIENNACL;
    185:   pc->ops->applytranspose      = 0;
    186:   pc->ops->setup               = PCSetUp_SAVIENNACL;
    187:   pc->ops->destroy             = PCDestroy_SAVIENNACL;
    188:   pc->ops->setfromoptions      = PCSetFromOptions_SAVIENNACL;
    189:   pc->ops->view                = 0;
    190:   pc->ops->applyrichardson     = 0;
    191:   pc->ops->applysymmetricleft  = 0;
    192:   pc->ops->applysymmetricright = 0;
    193:   PetscFunctionReturn(PETSC_SUCCESS);
    194: }