Actual source code: rowscalingviennacl.cxx


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

  4: /*
  5:    Include files needed for the ViennaCL row-scaling 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

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

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

 25: /* -------------------------------------------------------------------------- */
 26: /*
 27:    PCSetUp_ROWSCALINGVIENNACL - Prepares for the use of the ROWSCALINGVIENNACL preconditioner
 28:                                 by setting data structures and options.

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

 33:    Application Interface Routine: PCSetUp()

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

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

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

 72: /* -------------------------------------------------------------------------- */
 73: /*
 74:    PCApply_ROWSCALINGVIENNACL - Applies the ROWSCALINGVIENNACL preconditioner to a vector.

 76:    Input Parameters:
 77: .  pc - the preconditioner context
 78: .  x - input vector

 80:    Output Parameter:
 81: .  y - output vector

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

 92:   /*how to apply a certain fixed number of iterations?*/
 93:   PetscObjectTypeCompare((PetscObject)x,VECSEQVIENNACL,&flg1);
 94:   PetscObjectTypeCompare((PetscObject)y,VECSEQVIENNACL,&flg2);
 96:   if (!ilu->ROWSCALINGVIENNACL) {
 97:     PCSetUp_ROWSCALINGVIENNACL(pc);
 98:   }
 99:   VecSet(y,0.0);
100:   VecViennaCLGetArrayRead(x,&xarray);
101:   VecViennaCLGetArrayWrite(y,&yarray);
102:   try {
103: #if defined(PETSC_USE_COMPLEX)

105: #else
106:     *yarray = *xarray;
107:     ilu->ROWSCALINGVIENNACL->apply(*yarray);
108: #endif
109:   } catch(char * ex) {
110:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
111:   }
112:   VecViennaCLRestoreArrayRead(x,&xarray);
113:   VecViennaCLRestoreArrayWrite(y,&yarray);
114:   PetscObjectStateIncrease((PetscObject)y);
115:   return 0;
116: }
117: /* -------------------------------------------------------------------------- */
118: /*
119:    PCDestroy_ROWSCALINGVIENNACL - Destroys the private context for the ROWSCALINGVIENNACL preconditioner
120:    that was created with PCCreate_ROWSCALINGVIENNACL().

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

125:    Application Interface Routine: PCDestroy()
126: */
127: static PetscErrorCode PCDestroy_ROWSCALINGVIENNACL(PC pc)
128: {
129:   PC_ROWSCALINGVIENNACL  *rowscaling = (PC_ROWSCALINGVIENNACL*)pc->data;

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

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

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

153: /* -------------------------------------------------------------------------- */

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

158:    Level: advanced

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

162: M*/

164: PETSC_EXTERN PetscErrorCode PCCreate_ROWSCALINGVIENNACL(PC pc)
165: {
166:   PC_ROWSCALINGVIENNACL  *rowscaling;

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

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

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