1 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 8 9 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 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> 15 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; 22 23 /* 24 PCSetUp_CHOWILUVIENNACL - Prepares for the use of the CHOWILUVIENNACL preconditioner 25 by setting data structures and options. 26 27 Input Parameter: 28 . pc - the preconditioner context 29 30 Application Interface Routine: PCSetUp() 31 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 PC_CHOWILUVIENNACL *ilu = (PC_CHOWILUVIENNACL *)pc->data; 38 PetscBool flg = PETSC_FALSE; 39 Mat_SeqAIJViennaCL *gpustruct; 40 41 PetscFunctionBegin; 42 PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATSEQAIJVIENNACL, &flg)); 43 PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Currently only handles ViennaCL matrices"); 44 if (pc->setupcalled != 0) { 45 try { 46 delete ilu->CHOWILUVIENNACL; 47 } catch (char *ex) { SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex); } 48 } 49 try { 50 #if defined(PETSC_USE_COMPLEX) 51 gpustruct = NULL; 52 SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "No support for complex arithmetic in CHOWILUVIENNACL preconditioner"); 53 #else 54 PetscCall(MatViennaCLCopyToGPU(pc->pmat)); 55 gpustruct = (Mat_SeqAIJViennaCL *)(pc->pmat->spptr); 56 57 viennacl::linalg::chow_patel_tag ilu_tag; 58 ViennaCLAIJMatrix *mat = (ViennaCLAIJMatrix *)gpustruct->mat; 59 ilu->CHOWILUVIENNACL = new viennacl::linalg::chow_patel_ilu_precond<viennacl::compressed_matrix<PetscScalar>>(*mat, ilu_tag); 60 #endif 61 } catch (char *ex) { SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex); } 62 PetscFunctionReturn(0); 63 } 64 65 /* 66 PCApply_CHOWILUVIENNACL - Applies the CHOWILUVIENNACL preconditioner to a vector. 67 68 Input Parameters: 69 . pc - the preconditioner context 70 . x - input vector 71 72 Output Parameter: 73 . y - output vector 74 75 Application Interface Routine: PCApply() 76 */ 77 static PetscErrorCode PCApply_CHOWILUVIENNACL(PC pc, Vec x, Vec y) { 78 PC_CHOWILUVIENNACL *ilu = (PC_CHOWILUVIENNACL *)pc->data; 79 PetscBool flg1, flg2; 80 viennacl::vector<PetscScalar> const *xarray = NULL; 81 viennacl::vector<PetscScalar> *yarray = NULL; 82 83 PetscFunctionBegin; 84 /*how to apply a certain fixed number of iterations?*/ 85 PetscCall(PetscObjectTypeCompare((PetscObject)x, VECSEQVIENNACL, &flg1)); 86 PetscCall(PetscObjectTypeCompare((PetscObject)y, VECSEQVIENNACL, &flg2)); 87 PetscCheck((flg1 && flg2), PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Currently only handles ViennaCL vectors"); 88 if (!ilu->CHOWILUVIENNACL) PetscCall(PCSetUp_CHOWILUVIENNACL(pc)); 89 PetscCall(VecSet(y, 0.0)); 90 PetscCall(VecViennaCLGetArrayRead(x, &xarray)); 91 PetscCall(VecViennaCLGetArrayWrite(y, &yarray)); 92 try { 93 #if defined(PETSC_USE_COMPLEX) 94 95 #else 96 *yarray = *xarray; 97 ilu->CHOWILUVIENNACL->apply(*yarray); 98 #endif 99 } catch (char *ex) { SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex); } 100 PetscCall(VecViennaCLRestoreArrayRead(x, &xarray)); 101 PetscCall(VecViennaCLRestoreArrayWrite(y, &yarray)); 102 PetscCall(PetscObjectStateIncrease((PetscObject)y)); 103 PetscFunctionReturn(0); 104 } 105 106 /* 107 PCDestroy_CHOWILUVIENNACL - Destroys the private context for the CHOWILUVIENNACL preconditioner 108 that was created with PCCreate_CHOWILUVIENNACL(). 109 110 Input Parameter: 111 . pc - the preconditioner context 112 113 Application Interface Routine: PCDestroy() 114 */ 115 static PetscErrorCode PCDestroy_CHOWILUVIENNACL(PC pc) { 116 PC_CHOWILUVIENNACL *ilu = (PC_CHOWILUVIENNACL *)pc->data; 117 118 PetscFunctionBegin; 119 if (ilu->CHOWILUVIENNACL) { 120 try { 121 delete ilu->CHOWILUVIENNACL; 122 } catch (char *ex) { SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex); } 123 } 124 125 /* 126 Free the private data structure that was hanging off the PC 127 */ 128 PetscCall(PetscFree(pc->data)); 129 PetscFunctionReturn(0); 130 } 131 132 static PetscErrorCode PCSetFromOptions_CHOWILUVIENNACL(PC pc, PetscOptionItems *PetscOptionsObject) { 133 PetscFunctionBegin; 134 PetscOptionsHeadBegin(PetscOptionsObject, "CHOWILUVIENNACL options"); 135 PetscOptionsHeadEnd(); 136 PetscFunctionReturn(0); 137 } 138 139 /*MC 140 PCCHOWILUViennaCL - A smoothed agglomeration algorithm that can be used via the CUDA, OpenCL, and OpenMP backends of ViennaCL 141 142 Level: advanced 143 144 Developer Note: 145 This does not appear to be wired up with `PCRegisterType()` 146 147 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC` 148 M*/ 149 150 PETSC_EXTERN PetscErrorCode PCCreate_CHOWILUVIENNACL(PC pc) { 151 PC_CHOWILUVIENNACL *ilu; 152 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; 160 161 /* 162 Initialize the pointer to zero 163 Initialize number of v-cycles to default (1) 164 */ 165 ilu->CHOWILUVIENNACL = 0; 166 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(0); 184 } 185