1 2 /* -------------------------------------------------------------------- */ 3 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 10 11 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 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> 17 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; 24 25 /* -------------------------------------------------------------------------- */ 26 /* 27 PCSetUp_CHOWILUVIENNACL - Prepares for the use of the CHOWILUVIENNACL preconditioner 28 by setting data structures and options. 29 30 Input Parameter: 31 . pc - the preconditioner context 32 33 Application Interface Routine: PCSetUp() 34 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 PC_CHOWILUVIENNACL *ilu = (PC_CHOWILUVIENNACL *)pc->data; 41 PetscBool flg = PETSC_FALSE; 42 Mat_SeqAIJViennaCL *gpustruct; 43 44 PetscFunctionBegin; 45 PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATSEQAIJVIENNACL, &flg)); 46 PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Currently only handles ViennaCL matrices"); 47 if (pc->setupcalled != 0) { 48 try { 49 delete ilu->CHOWILUVIENNACL; 50 } catch (char *ex) { SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex); } 51 } 52 try { 53 #if defined(PETSC_USE_COMPLEX) 54 gpustruct = NULL; 55 SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "No support for complex arithmetic in CHOWILUVIENNACL preconditioner"); 56 #else 57 PetscCall(MatViennaCLCopyToGPU(pc->pmat)); 58 gpustruct = (Mat_SeqAIJViennaCL *)(pc->pmat->spptr); 59 60 viennacl::linalg::chow_patel_tag ilu_tag; 61 ViennaCLAIJMatrix *mat = (ViennaCLAIJMatrix *)gpustruct->mat; 62 ilu->CHOWILUVIENNACL = new viennacl::linalg::chow_patel_ilu_precond<viennacl::compressed_matrix<PetscScalar>>(*mat, ilu_tag); 63 #endif 64 } catch (char *ex) { SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex); } 65 PetscFunctionReturn(0); 66 } 67 68 /* -------------------------------------------------------------------------- */ 69 /* 70 PCApply_CHOWILUVIENNACL - Applies the CHOWILUVIENNACL preconditioner to a vector. 71 72 Input Parameters: 73 . pc - the preconditioner context 74 . x - input vector 75 76 Output Parameter: 77 . y - output vector 78 79 Application Interface Routine: PCApply() 80 */ 81 static PetscErrorCode PCApply_CHOWILUVIENNACL(PC pc, Vec x, Vec y) { 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; 86 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) 98 99 #else 100 *yarray = *xarray; 101 ilu->CHOWILUVIENNACL->apply(*yarray); 102 #endif 103 } catch (char *ex) { SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex); } 104 PetscCall(VecViennaCLRestoreArrayRead(x, &xarray)); 105 PetscCall(VecViennaCLRestoreArrayWrite(y, &yarray)); 106 PetscCall(PetscObjectStateIncrease((PetscObject)y)); 107 PetscFunctionReturn(0); 108 } 109 /* -------------------------------------------------------------------------- */ 110 /* 111 PCDestroy_CHOWILUVIENNACL - Destroys the private context for the CHOWILUVIENNACL preconditioner 112 that was created with PCCreate_CHOWILUVIENNACL(). 113 114 Input Parameter: 115 . pc - the preconditioner context 116 117 Application Interface Routine: PCDestroy() 118 */ 119 static PetscErrorCode PCDestroy_CHOWILUVIENNACL(PC pc) { 120 PC_CHOWILUVIENNACL *ilu = (PC_CHOWILUVIENNACL *)pc->data; 121 122 PetscFunctionBegin; 123 if (ilu->CHOWILUVIENNACL) { 124 try { 125 delete ilu->CHOWILUVIENNACL; 126 } catch (char *ex) { SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex); } 127 } 128 129 /* 130 Free the private data structure that was hanging off the PC 131 */ 132 PetscCall(PetscFree(pc->data)); 133 PetscFunctionReturn(0); 134 } 135 136 static PetscErrorCode PCSetFromOptions_CHOWILUVIENNACL(PC pc, PetscOptionItems *PetscOptionsObject) { 137 PetscFunctionBegin; 138 PetscOptionsHeadBegin(PetscOptionsObject, "CHOWILUVIENNACL options"); 139 PetscOptionsHeadEnd(); 140 PetscFunctionReturn(0); 141 } 142 143 /* -------------------------------------------------------------------------- */ 144 145 /*MC 146 PCCHOWILUViennaCL - A smoothed agglomeration algorithm that can be used via the CUDA, OpenCL, and OpenMP backends of ViennaCL 147 148 Level: advanced 149 150 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC` 151 152 M*/ 153 154 PETSC_EXTERN PetscErrorCode PCCreate_CHOWILUVIENNACL(PC pc) { 155 PC_CHOWILUVIENNACL *ilu; 156 157 PetscFunctionBegin; 158 /* 159 Creates the private data structure for this preconditioner and 160 attach it to the PC object. 161 */ 162 PetscCall(PetscNewLog(pc, &ilu)); 163 pc->data = (void *)ilu; 164 165 /* 166 Initialize the pointer to zero 167 Initialize number of v-cycles to default (1) 168 */ 169 ilu->CHOWILUVIENNACL = 0; 170 171 /* 172 Set the pointers for the functions that are provided above. 173 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 174 are called, they will automatically call these functions. Note we 175 choose not to provide a couple of these functions since they are 176 not needed. 177 */ 178 pc->ops->apply = PCApply_CHOWILUVIENNACL; 179 pc->ops->applytranspose = 0; 180 pc->ops->setup = PCSetUp_CHOWILUVIENNACL; 181 pc->ops->destroy = PCDestroy_CHOWILUVIENNACL; 182 pc->ops->setfromoptions = PCSetFromOptions_CHOWILUVIENNACL; 183 pc->ops->view = 0; 184 pc->ops->applyrichardson = 0; 185 pc->ops->applysymmetricleft = 0; 186 pc->ops->applysymmetricright = 0; 187 PetscFunctionReturn(0); 188 } 189