1 2 /* 3 Include files needed for the ViennaCL Smoothed Aggregation 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 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 9 #include <../src/mat/impls/aij/seq/aij.h> 10 #include <../src/vec/vec/impls/dvecimpl.h> 11 #include <../src/mat/impls/aij/seq/seqviennacl/viennaclmatimpl.h> 12 #include <../src/vec/vec/impls/seq/seqviennacl/viennaclvecimpl.h> 13 #include <viennacl/linalg/amg.hpp> 14 15 /* 16 Private context (data structure) for the SAVIENNACL preconditioner. 17 */ 18 typedef struct { 19 viennacl::linalg::amg_precond<viennacl::compressed_matrix<PetscScalar>> *SAVIENNACL; 20 } PC_SAVIENNACL; 21 22 /* 23 PCSetUp_SAVIENNACL - Prepares for the use of the SAVIENNACL preconditioner 24 by setting data structures and options. 25 26 Input Parameter: 27 . pc - the preconditioner context 28 29 Application Interface Routine: PCSetUp() 30 31 Note: 32 The interface routine PCSetUp() is not usually called directly by 33 the user, but instead is called by PCApply() if necessary. 34 */ 35 static PetscErrorCode PCSetUp_SAVIENNACL(PC pc) { 36 PC_SAVIENNACL *sa = (PC_SAVIENNACL *)pc->data; 37 PetscBool flg = PETSC_FALSE; 38 Mat_SeqAIJViennaCL *gpustruct; 39 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 sa->SAVIENNACL; 46 } catch (char *ex) { SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex); } 47 } 48 try { 49 #if defined(PETSC_USE_COMPLEX) 50 gpustruct = NULL; 51 SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "No support for complex arithmetic in SAVIENNACL preconditioner"); 52 #else 53 PetscCall(MatViennaCLCopyToGPU(pc->pmat)); 54 gpustruct = (Mat_SeqAIJViennaCL *)(pc->pmat->spptr); 55 56 viennacl::linalg::amg_tag amg_tag_sa_pmis; 57 amg_tag_sa_pmis.set_coarsening_method(viennacl::linalg::AMG_COARSENING_METHOD_MIS2_AGGREGATION); 58 amg_tag_sa_pmis.set_interpolation_method(viennacl::linalg::AMG_INTERPOLATION_METHOD_SMOOTHED_AGGREGATION); 59 ViennaCLAIJMatrix *mat = (ViennaCLAIJMatrix *)gpustruct->mat; 60 sa->SAVIENNACL = new viennacl::linalg::amg_precond<viennacl::compressed_matrix<PetscScalar>>(*mat, amg_tag_sa_pmis); 61 sa->SAVIENNACL->setup(); 62 #endif 63 } catch (char *ex) { SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex); } 64 PetscFunctionReturn(0); 65 } 66 67 /* 68 PCApply_SAVIENNACL - Applies the SAVIENNACL preconditioner to a vector. 69 70 Input Parameters: 71 . pc - the preconditioner context 72 . x - input vector 73 74 Output Parameter: 75 . y - output vector 76 77 Application Interface Routine: PCApply() 78 */ 79 static PetscErrorCode PCApply_SAVIENNACL(PC pc, Vec x, Vec y) { 80 PC_SAVIENNACL *sac = (PC_SAVIENNACL *)pc->data; 81 PetscBool flg1, flg2; 82 viennacl::vector<PetscScalar> const *xarray = NULL; 83 viennacl::vector<PetscScalar> *yarray = NULL; 84 85 PetscFunctionBegin; 86 /*how to apply a certain fixed number of iterations?*/ 87 PetscCall(PetscObjectTypeCompare((PetscObject)x, VECSEQVIENNACL, &flg1)); 88 PetscCall(PetscObjectTypeCompare((PetscObject)y, VECSEQVIENNACL, &flg2)); 89 PetscCheck((flg1 && flg2), PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Currently only handles ViennaCL vectors"); 90 if (!sac->SAVIENNACL) PetscCall(PCSetUp_SAVIENNACL(pc)); 91 PetscCall(VecViennaCLGetArrayRead(x, &xarray)); 92 PetscCall(VecViennaCLGetArrayWrite(y, &yarray)); 93 try { 94 #if !defined(PETSC_USE_COMPLEX) 95 *yarray = *xarray; 96 sac->SAVIENNACL->apply(*yarray); 97 #endif 98 } catch (char *ex) { SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex); } 99 PetscCall(VecViennaCLRestoreArrayRead(x, &xarray)); 100 PetscCall(VecViennaCLRestoreArrayWrite(y, &yarray)); 101 PetscCall(PetscObjectStateIncrease((PetscObject)y)); 102 PetscFunctionReturn(0); 103 } 104 105 /* 106 PCDestroy_SAVIENNACL - Destroys the private context for the SAVIENNACL preconditioner 107 that was created with PCCreate_SAVIENNACL(). 108 109 Input Parameter: 110 . pc - the preconditioner context 111 112 Application Interface Routine: PCDestroy() 113 */ 114 static PetscErrorCode PCDestroy_SAVIENNACL(PC pc) { 115 PC_SAVIENNACL *sac = (PC_SAVIENNACL *)pc->data; 116 117 PetscFunctionBegin; 118 if (sac->SAVIENNACL) { 119 try { 120 delete sac->SAVIENNACL; 121 } catch (char *ex) { SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex); } 122 } 123 124 /* 125 Free the private data structure that was hanging off the PC 126 */ 127 PetscCall(PetscFree(pc->data)); 128 PetscFunctionReturn(0); 129 } 130 131 static PetscErrorCode PCSetFromOptions_SAVIENNACL(PC pc, PetscOptionItems *PetscOptionsObject) { 132 PetscFunctionBegin; 133 PetscOptionsHeadBegin(PetscOptionsObject, "SAVIENNACL options"); 134 PetscOptionsHeadEnd(); 135 PetscFunctionReturn(0); 136 } 137 138 /*MC 139 PCSAViennaCL - A smoothed agglomeration algorithm that can be used via the CUDA, OpenCL, and OpenMP backends of ViennaCL 140 141 Level: advanced 142 143 Developer Note: 144 This `PCType` does not appear to be registered 145 146 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC` 147 M*/ 148 149 PETSC_EXTERN PetscErrorCode PCCreate_SAVIENNACL(PC pc) { 150 PC_SAVIENNACL *sac; 151 152 PetscFunctionBegin; 153 /* 154 Creates the private data structure for this preconditioner and 155 attach it to the PC object. 156 */ 157 PetscCall(PetscNew(&sac)); 158 pc->data = (void *)sac; 159 160 /* 161 Initialize the pointer to zero 162 Initialize number of v-cycles to default (1) 163 */ 164 sac->SAVIENNACL = 0; 165 166 /* 167 Set the pointers for the functions that are provided above. 168 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 169 are called, they will automatically call these functions. Note we 170 choose not to provide a couple of these functions since they are 171 not needed. 172 */ 173 pc->ops->apply = PCApply_SAVIENNACL; 174 pc->ops->applytranspose = 0; 175 pc->ops->setup = PCSetUp_SAVIENNACL; 176 pc->ops->destroy = PCDestroy_SAVIENNACL; 177 pc->ops->setfromoptions = PCSetFromOptions_SAVIENNACL; 178 pc->ops->view = 0; 179 pc->ops->applyrichardson = 0; 180 pc->ops->applysymmetricleft = 0; 181 pc->ops->applysymmetricright = 0; 182 PetscFunctionReturn(0); 183 } 184