Actual source code: chowiluviennacl.cxx


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

  4: /*
  5:    Include files needed for the ViennaCL Chow-Patel parallel ILU 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/detail/ilu/chow_patel_ilu.hpp>

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

 25: /* -------------------------------------------------------------------------- */
 26: /*
 27:    PCSetUp_CHOWILUVIENNACL - Prepares for the use of the CHOWILUVIENNACL 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_CHOWILUVIENNACL(PC pc)
 40: {
 41:   PC_CHOWILUVIENNACL *ilu = (PC_CHOWILUVIENNACL*)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 ilu->CHOWILUVIENNACL;
 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 CHOWILUVIENNACL preconditioner");
 58: #else
 59:     MatViennaCLCopyToGPU(pc->pmat);
 60:     gpustruct = (Mat_SeqAIJViennaCL*)(pc->pmat->spptr);

 62:     viennacl::linalg::chow_patel_tag ilu_tag;
 63:     ViennaCLAIJMatrix *mat = (ViennaCLAIJMatrix*)gpustruct->mat;
 64:     ilu->CHOWILUVIENNACL = new viennacl::linalg::chow_patel_ilu_precond<viennacl::compressed_matrix<PetscScalar> >(*mat, ilu_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_CHOWILUVIENNACL - Applies the CHOWILUVIENNACL 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_CHOWILUVIENNACL(PC pc,Vec x,Vec y)
 86: {
 87:   PC_CHOWILUVIENNACL            *ilu = (PC_CHOWILUVIENNACL*)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->CHOWILUVIENNACL) {
 97:     PCSetUp_CHOWILUVIENNACL(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->CHOWILUVIENNACL->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_CHOWILUVIENNACL - Destroys the private context for the CHOWILUVIENNACL preconditioner
120:    that was created with PCCreate_CHOWILUVIENNACL().

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

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

131:   if (ilu->CHOWILUVIENNACL) {
132:     try {
133:       delete ilu->CHOWILUVIENNACL;
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_CHOWILUVIENNACL(PetscOptionItems *PetscOptionsObject,PC pc)
147: {
148:   PetscOptionsHead(PetscOptionsObject,"CHOWILUVIENNACL options");
149:   PetscOptionsTail();
150:   return 0;
151: }

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

155: /*MC
156:      PCCHOWILUViennaCL  - A smoothed agglomeration algorithm 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_CHOWILUVIENNACL(PC pc)
165: {
166:   PC_CHOWILUVIENNACL  *ilu;

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

175:   /*
176:      Initialize the pointer to zero
177:      Initialize number of v-cycles to default (1)
178:   */
179:   ilu->CHOWILUVIENNACL = 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_CHOWILUVIENNACL;
189:   pc->ops->applytranspose      = 0;
190:   pc->ops->setup               = PCSetUp_CHOWILUVIENNACL;
191:   pc->ops->destroy             = PCDestroy_CHOWILUVIENNACL;
192:   pc->ops->setfromoptions      = PCSetFromOptions_CHOWILUVIENNACL;
193:   pc->ops->view                = 0;
194:   pc->ops->applyrichardson     = 0;
195:   pc->ops->applysymmetricleft  = 0;
196:   pc->ops->applysymmetricright = 0;
197:   return 0;
198: }