Actual source code: rowscalingviennacl.cxx

      1: /*
      2:    Include files needed for the ViennaCL row-scaling preconditioner:
      3:      pcimpl.h - private include file intended for use by all preconditioners
      4: */
      5: #define PETSC_SKIP_IMMINTRIN_H_CUDAWORKAROUND 1

      7: #include <petsc/private/pcimpl.h>
      8: #include <../src/mat/impls/aij/seq/aij.h>
      9: #include <../src/vec/vec/impls/dvecimpl.h>
     10: #include <../src/mat/impls/aij/seq/seqviennacl/viennaclmatimpl.h>
     11: #include <../src/vec/vec/impls/seq/seqviennacl/viennaclvecimpl.h>
     12: #include <viennacl/linalg/row_scaling.hpp>

     14: /*
     15:    Private context (data structure) for the ROWSCALINGVIENNACL preconditioner.
     16: */
     17: typedef struct {
     18:   viennacl::linalg::row_scaling<viennacl::compressed_matrix<PetscScalar>> *ROWSCALINGVIENNACL;
     19: } PC_ROWSCALINGVIENNACL;

     21: /*
     22:    PCSetUp_ROWSCALINGVIENNACL - Prepares for the use of the ROWSCALINGVIENNACL preconditioner
     23:                                 by setting data structures and options.

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

     28:    Application Interface Routine: PCSetUp()

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

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

     58:     viennacl::linalg::row_scaling_tag pc_tag(1);
     59:     ViennaCLAIJMatrix                *mat = (ViennaCLAIJMatrix *)gpustruct->mat;
     60:     rowscaling->ROWSCALINGVIENNACL        = new viennacl::linalg::row_scaling<viennacl::compressed_matrix<PetscScalar>>(*mat, pc_tag);
     61: #endif
     62:   } catch (char *ex) {
     63:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex);
     64:   }
     65:   PetscFunctionReturn(PETSC_SUCCESS);
     66: }

     68: /*
     69:    PCApply_ROWSCALINGVIENNACL - Applies the ROWSCALINGVIENNACL preconditioner to a vector.

     71:    Input Parameters:
     72: .  pc - the preconditioner context
     73: .  x - input vector

     75:    Output Parameter:
     76: .  y - output vector

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

     87:   PetscFunctionBegin;
     88:   /*how to apply a certain fixed number of iterations?*/
     89:   PetscCall(PetscObjectTypeCompare((PetscObject)x, VECSEQVIENNACL, &flg1));
     90:   PetscCall(PetscObjectTypeCompare((PetscObject)y, VECSEQVIENNACL, &flg2));
     91:   PetscCheck(flg1 && flg2, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Currently only handles ViennaCL vectors");
     92:   if (!ilu->ROWSCALINGVIENNACL) PetscCall(PCSetUp_ROWSCALINGVIENNACL(pc));
     93:   PetscCall(VecSet(y, 0.0));
     94:   PetscCall(VecViennaCLGetArrayRead(x, &xarray));
     95:   PetscCall(VecViennaCLGetArrayWrite(y, &yarray));
     96:   try {
     97: #if defined(PETSC_USE_COMPLEX)

     99: #else
    100:     *yarray = *xarray;
    101:     ilu->ROWSCALINGVIENNACL->apply(*yarray);
    102: #endif
    103:   } catch (char *ex) {
    104:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex);
    105:   }
    106:   PetscCall(VecViennaCLRestoreArrayRead(x, &xarray));
    107:   PetscCall(VecViennaCLRestoreArrayWrite(y, &yarray));
    108:   PetscCall(PetscObjectStateIncrease((PetscObject)y));
    109:   PetscFunctionReturn(PETSC_SUCCESS);
    110: }

    112: /*
    113:    PCDestroy_ROWSCALINGVIENNACL - Destroys the private context for the ROWSCALINGVIENNACL preconditioner
    114:    that was created with PCCreate_ROWSCALINGVIENNACL().

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

    119:    Application Interface Routine: PCDestroy()
    120: */
    121: static PetscErrorCode PCDestroy_ROWSCALINGVIENNACL(PC pc)
    122: {
    123:   PC_ROWSCALINGVIENNACL *rowscaling = (PC_ROWSCALINGVIENNACL *)pc->data;

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

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

    141: static PetscErrorCode PCSetFromOptions_ROWSCALINGVIENNACL(PC pc, PetscOptionItems *PetscOptionsObject)
    142: {
    143:   PetscFunctionBegin;
    144:   PetscOptionsHeadBegin(PetscOptionsObject, "ROWSCALINGVIENNACL options");
    145:   PetscOptionsHeadEnd();
    146:   PetscFunctionReturn(PETSC_SUCCESS);
    147: }

    149: /*MC
    150:   PCROWSCALINGVIENNACLL  - A diagonal preconditioner (scaling rows of matrices by their norm) that can be used via the CUDA,
    151:   OpenCL, and OpenMP backends of ViennaCL

    153:   Level: advanced

    155: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`
    156: M*/
    157: PETSC_EXTERN PetscErrorCode PCCreate_ROWSCALINGVIENNACL(PC pc)
    158: {
    159:   PC_ROWSCALINGVIENNACL *rowscaling;

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

    169:   /*
    170:      Initialize the pointer to zero
    171:      Initialize number of v-cycles to default (1)
    172:   */
    173:   rowscaling->ROWSCALINGVIENNACL = 0;

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