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