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;


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

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

 34:    Application Interface Routine: PCSetUp()

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

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

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

 75: /* -------------------------------------------------------------------------- */
 76: /*
 77:    PCApply_ROWSCALINGVIENNACL - Applies the ROWSCALINGVIENNACL preconditioner to a vector.

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

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

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

 97:   /*how to apply a certain fixed number of iterations?*/
 98:   PetscObjectTypeCompare((PetscObject)x,VECSEQVIENNACL,&flg1);
 99:   PetscObjectTypeCompare((PetscObject)y,VECSEQVIENNACL,&flg2);
100:   if (!(flg1 && flg2)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP, "Currently only handles ViennaCL vectors");
101:   if (!ilu->ROWSCALINGVIENNACL) {
102:     PCSetUp_ROWSCALINGVIENNACL(pc);
103:   }
104:   VecSet(y,0.0);
105:   VecViennaCLGetArrayRead(x,&xarray);
106:   VecViennaCLGetArrayWrite(y,&yarray);
107:   try {
108: #if defined(PETSC_USE_COMPLEX)

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

127:    Input Parameter:
128: .  pc - the preconditioner context

130:    Application Interface Routine: PCDestroy()
131: */
132: static PetscErrorCode PCDestroy_ROWSCALINGVIENNACL(PC pc)
133: {
134:   PC_ROWSCALINGVIENNACL  *rowscaling = (PC_ROWSCALINGVIENNACL*)pc->data;

138:   if (rowscaling->ROWSCALINGVIENNACL) {
139:     try {
140:       delete rowscaling->ROWSCALINGVIENNACL;
141:     } catch(char *ex) {
142:       SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
143:     }
144:   }

146:   /*
147:       Free the private data structure that was hanging off the PC
148:   */
149:   PetscFree(pc->data);
150:   return(0);
151: }

153: static PetscErrorCode PCSetFromOptions_ROWSCALINGVIENNACL(PetscOptionItems *PetscOptionsObject,PC pc)
154: {

158:   PetscOptionsHead(PetscOptionsObject,"ROWSCALINGVIENNACL options");
159:   PetscOptionsTail();
160:   return(0);
161: }

163: /* -------------------------------------------------------------------------- */


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

169:    Level: advanced

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

173: M*/

175: PETSC_EXTERN PetscErrorCode PCCreate_ROWSCALINGVIENNACL(PC pc)
176: {
177:   PC_ROWSCALINGVIENNACL  *rowscaling;

181:   /*
182:      Creates the private data structure for this preconditioner and
183:      attach it to the PC object.
184:   */
185:   PetscNewLog(pc,&rowscaling);
186:   pc->data = (void*)rowscaling;

188:   /*
189:      Initialize the pointer to zero
190:      Initialize number of v-cycles to default (1)
191:   */
192:   rowscaling->ROWSCALINGVIENNACL = 0;

194:   /*
195:       Set the pointers for the functions that are provided above.
196:       Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
197:       are called, they will automatically call these functions.  Note we
198:       choose not to provide a couple of these functions since they are
199:       not needed.
200:   */
201:   pc->ops->apply               = PCApply_ROWSCALINGVIENNACL;
202:   pc->ops->applytranspose      = 0;
203:   pc->ops->setup               = PCSetUp_ROWSCALINGVIENNACL;
204:   pc->ops->destroy             = PCDestroy_ROWSCALINGVIENNACL;
205:   pc->ops->setfromoptions      = PCSetFromOptions_ROWSCALINGVIENNACL;
206:   pc->ops->view                = 0;
207:   pc->ops->applyrichardson     = 0;
208:   pc->ops->applysymmetricleft  = 0;
209:   pc->ops->applysymmetricright = 0;
210:   return(0);
211: }