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:  #include <petsc/private/pcimpl.h>
 10:  #include <../src/mat/impls/aij/seq/aij.h>
 11:  #include <../src/vec/vec/impls/dvecimpl.h>
 12:  #include <../src/mat/impls/aij/seq/seqviennacl/viennaclmatimpl.h>
 13:  #include <../src/vec/vec/impls/seq/seqviennacl/viennaclvecimpl.h>
 14: #include <viennacl/linalg/row_scaling.hpp>

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


 24: /* -------------------------------------------------------------------------- */
 25: /*
 26:    PCSetUp_ROWSCALINGVIENNACL - Prepares for the use of the ROWSCALINGVIENNACL 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_ROWSCALINGVIENNACL(PC pc)
 39: {
 40:   PC_ROWSCALINGVIENNACL  *rowscaling = (PC_ROWSCALINGVIENNACL*)pc->data;
 41:   PetscBool              flg = PETSC_FALSE;
 42:   PetscErrorCode         ierr;
 43:   Mat_SeqAIJViennaCL     *gpustruct;

 46:   PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQAIJVIENNACL,&flg);
 47:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Currently only handles ViennaCL matrices");
 48:   if (pc->setupcalled != 0) {
 49:     try {
 50:       delete rowscaling->ROWSCALINGVIENNACL;
 51:     } catch(char *ex) {
 52:       SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
 53:     }
 54:   }
 55:   try {
 56: #if defined(PETSC_USE_COMPLEX)
 57:     gpustruct = NULL;
 58:     SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"No support for complex arithmetic in ROWSCALINGVIENNACL preconditioner");
 59: #else
 60:     MatViennaCLCopyToGPU(pc->pmat);
 61:     gpustruct = (Mat_SeqAIJViennaCL*)(pc->pmat->spptr);
 62: 
 63:     viennacl::linalg::row_scaling_tag pc_tag(1);
 64:     ViennaCLAIJMatrix *mat = (ViennaCLAIJMatrix*)gpustruct->mat;
 65:     rowscaling->ROWSCALINGVIENNACL = new viennacl::linalg::row_scaling<viennacl::compressed_matrix<PetscScalar> >(*mat, pc_tag);
 66: #endif
 67:   } catch(char *ex) {
 68:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
 69:   }
 70:   return(0);
 71: }

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

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

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

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

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

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

125:    Input Parameter:
126: .  pc - the preconditioner context

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

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

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

151: static PetscErrorCode PCSetFromOptions_ROWSCALINGVIENNACL(PetscOptionItems *PetscOptionsObject,PC pc)
152: {

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

161: /* -------------------------------------------------------------------------- */


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

167:    Level: advanced

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

171: M*/

173: PETSC_EXTERN PetscErrorCode PCCreate_ROWSCALINGVIENNACL(PC pc)
174: {
175:   PC_ROWSCALINGVIENNACL  *rowscaling;

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

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

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