Actual source code: chowiluviennacl.cxx

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

  7: #include <petsc/private/pcimpl.h>
  8: #include <../src/mat/impls/aij/seq/aij.h>
  9: #include <../src/vec/vec/impls/dvecimpl.h>
 10: #include <../src/mat/impls/aij/seq/seqviennacl/viennaclmatimpl.h>
 11: #include <../src/vec/vec/impls/seq/seqviennacl/viennaclvecimpl.h>
 12: #include <viennacl/linalg/detail/ilu/chow_patel_ilu.hpp>

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

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

 25:    Input Parameter:
 26: .  pc - the preconditioner context

 28:    Application Interface Routine: PCSetUp()

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

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

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

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

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

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

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

 87:   PetscFunctionBegin;
 88:   /*how to apply a certain fixed number of iterations?*/
 89:   PetscCall(PetscObjectTypeCompare((PetscObject)x, VECSEQVIENNACL, &flg1));
 90:   PetscCall(PetscObjectTypeCompare((PetscObject)y, VECSEQVIENNACL, &flg2));
 91:   PetscCheck((flg1 && flg2), PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Currently only handles ViennaCL vectors");
 92:   if (!ilu->CHOWILUVIENNACL) PetscCall(PCSetUp_CHOWILUVIENNACL(pc));
 93:   PetscCall(VecSet(y, 0.0));
 94:   PetscCall(VecViennaCLGetArrayRead(x, &xarray));
 95:   PetscCall(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:   PetscCall(VecViennaCLRestoreArrayRead(x, &xarray));
107:   PetscCall(VecViennaCLRestoreArrayWrite(y, &yarray));
108:   PetscCall(PetscObjectStateIncrease((PetscObject)y));
109:   PetscFunctionReturn(PETSC_SUCCESS);
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:   PetscFunctionBegin;
126:   if (ilu->CHOWILUVIENNACL) {
127:     try {
128:       delete ilu->CHOWILUVIENNACL;
129:     } catch (char *ex) {
130:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex);
131:     }
132:   }

134:   /*
135:       Free the private data structure that was hanging off the PC
136:   */
137:   PetscCall(PetscFree(pc->data));
138:   PetscFunctionReturn(PETSC_SUCCESS);
139: }

141: static PetscErrorCode PCSetFromOptions_CHOWILUVIENNACL(PC pc, PetscOptionItems *PetscOptionsObject)
142: {
143:   PetscFunctionBegin;
144:   PetscOptionsHeadBegin(PetscOptionsObject, "CHOWILUVIENNACL options");
145:   PetscOptionsHeadEnd();
146:   PetscFunctionReturn(PETSC_SUCCESS);
147: }

149: PETSC_EXTERN PetscErrorCode PCCreate_CHOWILUVIENNACL(PC pc)
150: {
151:   PC_CHOWILUVIENNACL *ilu;

153:   PetscFunctionBegin;
154:   /*
155:      Creates the private data structure for this preconditioner and
156:      attach it to the PC object.
157:   */
158:   PetscCall(PetscNew(&ilu));
159:   pc->data = (void *)ilu;

161:   /*
162:      Initialize the pointer to zero
163:      Initialize number of v-cycles to default (1)
164:   */
165:   ilu->CHOWILUVIENNACL = 0;

167:   /*
168:       Set the pointers for the functions that are provided above.
169:       Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
170:       are called, they will automatically call these functions.  Note we
171:       choose not to provide a couple of these functions since they are
172:       not needed.
173:   */
174:   pc->ops->apply               = PCApply_CHOWILUVIENNACL;
175:   pc->ops->applytranspose      = 0;
176:   pc->ops->setup               = PCSetUp_CHOWILUVIENNACL;
177:   pc->ops->destroy             = PCDestroy_CHOWILUVIENNACL;
178:   pc->ops->setfromoptions      = PCSetFromOptions_CHOWILUVIENNACL;
179:   pc->ops->view                = 0;
180:   pc->ops->applyrichardson     = 0;
181:   pc->ops->applysymmetricleft  = 0;
182:   pc->ops->applysymmetricright = 0;
183:   PetscFunctionReturn(PETSC_SUCCESS);
184: }