Actual source code: saviennacl.cxx


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

  4: /*
  5:    Include files needed for the ViennaCL Smoothed Aggregation 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
 10: #include <petsc/private/pcimpl.h>
 11: #include <../src/mat/impls/aij/seq/aij.h>
 12: #include <../src/vec/vec/impls/dvecimpl.h>
 13: #include <../src/mat/impls/aij/seq/seqviennacl/viennaclmatimpl.h>
 14: #include <../src/vec/vec/impls/seq/seqviennacl/viennaclvecimpl.h>
 15: #include <viennacl/linalg/amg.hpp>

 17: /*
 18:    Private context (data structure) for the SAVIENNACL preconditioner.
 19: */
 20: typedef struct {
 21:   viennacl::linalg::amg_precond<viennacl::compressed_matrix<PetscScalar> > *SAVIENNACL;
 22: } PC_SAVIENNACL;

 24: /* -------------------------------------------------------------------------- */
 25: /*
 26:    PCSetUp_SAVIENNACL - Prepares for the use of the SAVIENNACL preconditioner
 27:                         by setting data structures and options.

 29:    Input Parameter:
 30: .  pc - the preconditioner context

 32:    Application Interface Routine: PCSetUp()

 34:    Notes:
 35:    The interface routine PCSetUp() is not usually called directly by
 36:    the user, but instead is called by PCApply() if necessary.
 37: */
 38: static PetscErrorCode PCSetUp_SAVIENNACL(PC pc)
 39: {
 40:   PC_SAVIENNACL      *sa = (PC_SAVIENNACL*)pc->data;
 41:   PetscBool          flg = PETSC_FALSE;
 42:   Mat_SeqAIJViennaCL *gpustruct;

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

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

 74: /* -------------------------------------------------------------------------- */
 75: /*
 76:    PCApply_SAVIENNACL - Applies the SAVIENNACL preconditioner to a vector.

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

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

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

 94:   /*how to apply a certain fixed number of iterations?*/
 95:   PetscObjectTypeCompare((PetscObject)x,VECSEQVIENNACL,&flg1);
 96:   PetscObjectTypeCompare((PetscObject)y,VECSEQVIENNACL,&flg2);
 98:   if (!sac->SAVIENNACL) {
 99:     PCSetUp_SAVIENNACL(pc);
100:   }
101:   VecViennaCLGetArrayRead(x,&xarray);
102:   VecViennaCLGetArrayWrite(y,&yarray);
103:   try {
104: #if !defined(PETSC_USE_COMPLEX)
105:     *yarray = *xarray;
106:     sac->SAVIENNACL->apply(*yarray);
107: #endif
108:   } catch(char * ex) {
109:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
110:   }
111:   VecViennaCLRestoreArrayRead(x,&xarray);
112:   VecViennaCLRestoreArrayWrite(y,&yarray);
113:   PetscObjectStateIncrease((PetscObject)y);
114:   return 0;
115: }
116: /* -------------------------------------------------------------------------- */
117: /*
118:    PCDestroy_SAVIENNACL - Destroys the private context for the SAVIENNACL preconditioner
119:    that was created with PCCreate_SAVIENNACL().

121:    Input Parameter:
122: .  pc - the preconditioner context

124:    Application Interface Routine: PCDestroy()
125: */
126: static PetscErrorCode PCDestroy_SAVIENNACL(PC pc)
127: {
128:   PC_SAVIENNACL  *sac = (PC_SAVIENNACL*)pc->data;

130:   if (sac->SAVIENNACL) {
131:     try {
132:       delete sac->SAVIENNACL;
133:     } catch(char * ex) {
134:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
135:     }
136:   }

138:   /*
139:       Free the private data structure that was hanging off the PC
140:   */
141:   PetscFree(pc->data);
142:   return 0;
143: }

145: static PetscErrorCode PCSetFromOptions_SAVIENNACL(PetscOptionItems *PetscOptionsObject,PC pc)
146: {
147:   PetscOptionsHead(PetscOptionsObject,"SAVIENNACL options");
148:   PetscOptionsTail();
149:   return 0;
150: }

152: /* -------------------------------------------------------------------------- */

154: /*MC
155:      PCSAViennaCL  - A smoothed agglomeration algorithm that can be used via the CUDA, OpenCL, and OpenMP backends of ViennaCL

157:    Level: advanced

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

161: M*/

163: PETSC_EXTERN PetscErrorCode PCCreate_SAVIENNACL(PC pc)
164: {
165:   PC_SAVIENNACL  *sac;

167:   /*
168:      Creates the private data structure for this preconditioner and
169:      attach it to the PC object.
170:   */
171:   PetscNewLog(pc,&sac);
172:   pc->data = (void*)sac;

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

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