Actual source code: chowiluviennacl.cxx


  2: /*
  3:    Include files needed for the ViennaCL Chow-Patel parallel ILU preconditioner:
  4:      pcimpl.h - private include file intended for use by all preconditioners
  5: */
  6: #define PETSC_SKIP_SPINLOCK
  7: #define PETSC_SKIP_IMMINTRIN_H_CUDAWORKAROUND 1

  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/detail/ilu/chow_patel_ilu.hpp>

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

 23: /*
 24:    PCSetUp_CHOWILUVIENNACL - Prepares for the use of the CHOWILUVIENNACL preconditioner
 25:                              by setting data structures and options.

 27:    Input Parameter:
 28: .  pc - the preconditioner context

 30:    Application Interface Routine: PCSetUp()

 32:    Note:
 33:    The interface routine PCSetUp() is not usually called directly by
 34:    the user, but instead is called by PCApply() if necessary.
 35: */
 36: static PetscErrorCode PCSetUp_CHOWILUVIENNACL(PC pc)
 37: {
 38:   PC_CHOWILUVIENNACL *ilu = (PC_CHOWILUVIENNACL *)pc->data;
 39:   PetscBool           flg = PETSC_FALSE;
 40:   Mat_SeqAIJViennaCL *gpustruct;

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

 59:     viennacl::linalg::chow_patel_tag ilu_tag;
 60:     ViennaCLAIJMatrix               *mat = (ViennaCLAIJMatrix *)gpustruct->mat;
 61:     ilu->CHOWILUVIENNACL                 = new viennacl::linalg::chow_patel_ilu_precond<viennacl::compressed_matrix<PetscScalar>>(*mat, ilu_tag);
 62: #endif
 63:   } catch (char *ex) {
 64:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex);
 65:   }
 66:   return 0;
 67: }

 69: /*
 70:    PCApply_CHOWILUVIENNACL - Applies the CHOWILUVIENNACL preconditioner to a vector.

 72:    Input Parameters:
 73: .  pc - the preconditioner context
 74: .  x - input vector

 76:    Output Parameter:
 77: .  y - output vector

 79:    Application Interface Routine: PCApply()
 80:  */
 81: static PetscErrorCode PCApply_CHOWILUVIENNACL(PC pc, Vec x, Vec y)
 82: {
 83:   PC_CHOWILUVIENNACL                  *ilu = (PC_CHOWILUVIENNACL *)pc->data;
 84:   PetscBool                            flg1, flg2;
 85:   viennacl::vector<PetscScalar> const *xarray = NULL;
 86:   viennacl::vector<PetscScalar>       *yarray = NULL;

 88:   /*how to apply a certain fixed number of iterations?*/
 89:   PetscObjectTypeCompare((PetscObject)x, VECSEQVIENNACL, &flg1);
 90:   PetscObjectTypeCompare((PetscObject)y, VECSEQVIENNACL, &flg2);
 92:   if (!ilu->CHOWILUVIENNACL) PCSetUp_CHOWILUVIENNACL(pc);
 93:   VecSet(y, 0.0);
 94:   VecViennaCLGetArrayRead(x, &xarray);
 95:   VecViennaCLGetArrayWrite(y, &yarray);
 96:   try {
 97: #if defined(PETSC_USE_COMPLEX)

 99: #else
100:     *yarray                              = *xarray;
101:     ilu->CHOWILUVIENNACL->apply(*yarray);
102: #endif
103:   } catch (char *ex) {
104:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex);
105:   }
106:   VecViennaCLRestoreArrayRead(x, &xarray);
107:   VecViennaCLRestoreArrayWrite(y, &yarray);
108:   PetscObjectStateIncrease((PetscObject)y);
109:   return 0;
110: }

112: /*
113:    PCDestroy_CHOWILUVIENNACL - Destroys the private context for the CHOWILUVIENNACL preconditioner
114:    that was created with PCCreate_CHOWILUVIENNACL().

116:    Input Parameter:
117: .  pc - the preconditioner context

119:    Application Interface Routine: PCDestroy()
120: */
121: static PetscErrorCode PCDestroy_CHOWILUVIENNACL(PC pc)
122: {
123:   PC_CHOWILUVIENNACL *ilu = (PC_CHOWILUVIENNACL *)pc->data;

125:   if (ilu->CHOWILUVIENNACL) {
126:     try {
127:       delete ilu->CHOWILUVIENNACL;
128:     } catch (char *ex) {
129:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex);
130:     }
131:   }

133:   /*
134:       Free the private data structure that was hanging off the PC
135:   */
136:   PetscFree(pc->data);
137:   return 0;
138: }

140: static PetscErrorCode PCSetFromOptions_CHOWILUVIENNACL(PC pc, PetscOptionItems *PetscOptionsObject)
141: {
142:   PetscOptionsHeadBegin(PetscOptionsObject, "CHOWILUVIENNACL options");
143:   PetscOptionsHeadEnd();
144:   return 0;
145: }

147: /*MC
148:      PCCHOWILUViennaCL  - A smoothed agglomeration algorithm that can be used via the CUDA, OpenCL, and OpenMP backends of ViennaCL

150:    Level: advanced

152:    Developer Note:
153:    This does not appear to be wired up with `PCRegisterType()`

155: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`
156: M*/

158: PETSC_EXTERN PetscErrorCode PCCreate_CHOWILUVIENNACL(PC pc)
159: {
160:   PC_CHOWILUVIENNACL *ilu;

162:   /*
163:      Creates the private data structure for this preconditioner and
164:      attach it to the PC object.
165:   */
166:   PetscNew(&ilu);
167:   pc->data = (void *)ilu;

169:   /*
170:      Initialize the pointer to zero
171:      Initialize number of v-cycles to default (1)
172:   */
173:   ilu->CHOWILUVIENNACL = 0;

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