Actual source code: rowscalingviennacl.cxx


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

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

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

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

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

 29:    Application Interface Routine: PCSetUp()

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

153:    Level: advanced

155:    Developer Note:
156:    This `PCType` does not appear to be registered

158: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`
159: M*/

161: PETSC_EXTERN PetscErrorCode PCCreate_ROWSCALINGVIENNACL(PC pc)
162: {
163:   PC_ROWSCALINGVIENNACL *rowscaling;

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

173:   /*
174:      Initialize the pointer to zero
175:      Initialize number of v-cycles to default (1)
176:   */
177:   rowscaling->ROWSCALINGVIENNACL = 0;

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