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 { 38 PC_CHOWILUVIENNACL *ilu = (PC_CHOWILUVIENNACL *)pc->data; 39 PetscBool flg = PETSC_FALSE; 40 Mat_SeqAIJViennaCL *gpustruct; 41 42 PetscFunctionBegin; 43 PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATSEQAIJVIENNACL, &flg)); 44 PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Currently only handles ViennaCL matrices"); 45 if (pc->setupcalled != 0) { 46 try { 47 delete ilu->CHOWILUVIENNACL; 48 } catch (char *ex) { 49 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex); 50 } 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) { 65 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex); 66 } 67 PetscFunctionReturn(0); 68 } 69 70 /* 71 PCApply_CHOWILUVIENNACL - Applies the CHOWILUVIENNACL preconditioner to a vector. 72 73 Input Parameters: 74 . pc - the preconditioner context 75 . x - input vector 76 77 Output Parameter: 78 . y - output vector 79 80 Application Interface Routine: PCApply() 81 */ 82 static PetscErrorCode PCApply_CHOWILUVIENNACL(PC pc, Vec x, Vec y) 83 { 84 PC_CHOWILUVIENNACL *ilu = (PC_CHOWILUVIENNACL *)pc->data; 85 PetscBool flg1, flg2; 86 viennacl::vector<PetscScalar> const *xarray = NULL; 87 viennacl::vector<PetscScalar> *yarray = NULL; 88 89 PetscFunctionBegin; 90 /*how to apply a certain fixed number of iterations?*/ 91 PetscCall(PetscObjectTypeCompare((PetscObject)x, VECSEQVIENNACL, &flg1)); 92 PetscCall(PetscObjectTypeCompare((PetscObject)y, VECSEQVIENNACL, &flg2)); 93 PetscCheck((flg1 && flg2), PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Currently only handles ViennaCL vectors"); 94 if (!ilu->CHOWILUVIENNACL) PetscCall(PCSetUp_CHOWILUVIENNACL(pc)); 95 PetscCall(VecSet(y, 0.0)); 96 PetscCall(VecViennaCLGetArrayRead(x, &xarray)); 97 PetscCall(VecViennaCLGetArrayWrite(y, &yarray)); 98 try { 99 #if defined(PETSC_USE_COMPLEX) 100 101 #else 102 *yarray = *xarray; 103 ilu->CHOWILUVIENNACL->apply(*yarray); 104 #endif 105 } catch (char *ex) { 106 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex); 107 } 108 PetscCall(VecViennaCLRestoreArrayRead(x, &xarray)); 109 PetscCall(VecViennaCLRestoreArrayWrite(y, &yarray)); 110 PetscCall(PetscObjectStateIncrease((PetscObject)y)); 111 PetscFunctionReturn(0); 112 } 113 114 /* 115 PCDestroy_CHOWILUVIENNACL - Destroys the private context for the CHOWILUVIENNACL preconditioner 116 that was created with PCCreate_CHOWILUVIENNACL(). 117 118 Input Parameter: 119 . pc - the preconditioner context 120 121 Application Interface Routine: PCDestroy() 122 */ 123 static PetscErrorCode PCDestroy_CHOWILUVIENNACL(PC pc) 124 { 125 PC_CHOWILUVIENNACL *ilu = (PC_CHOWILUVIENNACL *)pc->data; 126 127 PetscFunctionBegin; 128 if (ilu->CHOWILUVIENNACL) { 129 try { 130 delete ilu->CHOWILUVIENNACL; 131 } catch (char *ex) { 132 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex); 133 } 134 } 135 136 /* 137 Free the private data structure that was hanging off the PC 138 */ 139 PetscCall(PetscFree(pc->data)); 140 PetscFunctionReturn(0); 141 } 142 143 static PetscErrorCode PCSetFromOptions_CHOWILUVIENNACL(PC pc, PetscOptionItems *PetscOptionsObject) 144 { 145 PetscFunctionBegin; 146 PetscOptionsHeadBegin(PetscOptionsObject, "CHOWILUVIENNACL options"); 147 PetscOptionsHeadEnd(); 148 PetscFunctionReturn(0); 149 } 150 151 /*MC 152 PCCHOWILUViennaCL - A smoothed agglomeration algorithm that can be used via the CUDA, OpenCL, and OpenMP backends of ViennaCL 153 154 Level: advanced 155 156 Developer Note: 157 This does not appear to be wired up with `PCRegisterType()` 158 159 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC` 160 M*/ 161 162 PETSC_EXTERN PetscErrorCode PCCreate_CHOWILUVIENNACL(PC pc) 163 { 164 PC_CHOWILUVIENNACL *ilu; 165 166 PetscFunctionBegin; 167 /* 168 Creates the private data structure for this preconditioner and 169 attach it to the PC object. 170 */ 171 PetscCall(PetscNew(&ilu)); 172 pc->data = (void *)ilu; 173 174 /* 175 Initialize the pointer to zero 176 Initialize number of v-cycles to default (1) 177 */ 178 ilu->CHOWILUVIENNACL = 0; 179 180 /* 181 Set the pointers for the functions that are provided above. 182 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 183 are called, they will automatically call these functions. Note we 184 choose not to provide a couple of these functions since they are 185 not needed. 186 */ 187 pc->ops->apply = PCApply_CHOWILUVIENNACL; 188 pc->ops->applytranspose = 0; 189 pc->ops->setup = PCSetUp_CHOWILUVIENNACL; 190 pc->ops->destroy = PCDestroy_CHOWILUVIENNACL; 191 pc->ops->setfromoptions = PCSetFromOptions_CHOWILUVIENNACL; 192 pc->ops->view = 0; 193 pc->ops->applyrichardson = 0; 194 pc->ops->applysymmetricleft = 0; 195 pc->ops->applysymmetricright = 0; 196 PetscFunctionReturn(0); 197 } 198