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 != 0) {
 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: }