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 { 40 PC_SAVIENNACL *sa = (PC_SAVIENNACL*)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 sa->SAVIENNACL; 50 } catch(char *ex) { 51 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex); 52 } 53 } 54 try { 55 #if defined(PETSC_USE_COMPLEX) 56 gpustruct = NULL; 57 SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"No support for complex arithmetic in SAVIENNACL preconditioner"); 58 #else 59 PetscCall(MatViennaCLCopyToGPU(pc->pmat)); 60 gpustruct = (Mat_SeqAIJViennaCL*)(pc->pmat->spptr); 61 62 viennacl::linalg::amg_tag amg_tag_sa_pmis; 63 amg_tag_sa_pmis.set_coarsening_method(viennacl::linalg::AMG_COARSENING_METHOD_MIS2_AGGREGATION); 64 amg_tag_sa_pmis.set_interpolation_method(viennacl::linalg::AMG_INTERPOLATION_METHOD_SMOOTHED_AGGREGATION); 65 ViennaCLAIJMatrix *mat = (ViennaCLAIJMatrix*)gpustruct->mat; 66 sa->SAVIENNACL = new viennacl::linalg::amg_precond<viennacl::compressed_matrix<PetscScalar> >(*mat, amg_tag_sa_pmis); 67 sa->SAVIENNACL->setup(); 68 #endif 69 } catch(char *ex) { 70 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex); 71 } 72 PetscFunctionReturn(0); 73 } 74 75 /* -------------------------------------------------------------------------- */ 76 /* 77 PCApply_SAVIENNACL - Applies the SAVIENNACL preconditioner to a vector. 78 79 Input Parameters: 80 . pc - the preconditioner context 81 . x - input vector 82 83 Output Parameter: 84 . y - output vector 85 86 Application Interface Routine: PCApply() 87 */ 88 static PetscErrorCode PCApply_SAVIENNACL(PC pc,Vec x,Vec y) 89 { 90 PC_SAVIENNACL *sac = (PC_SAVIENNACL*)pc->data; 91 PetscBool flg1,flg2; 92 viennacl::vector<PetscScalar> const *xarray=NULL; 93 viennacl::vector<PetscScalar> *yarray=NULL; 94 95 PetscFunctionBegin; 96 /*how to apply a certain fixed number of iterations?*/ 97 PetscCall(PetscObjectTypeCompare((PetscObject)x,VECSEQVIENNACL,&flg1)); 98 PetscCall(PetscObjectTypeCompare((PetscObject)y,VECSEQVIENNACL,&flg2)); 99 PetscCheck((flg1 && flg2),PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP, "Currently only handles ViennaCL vectors"); 100 if (!sac->SAVIENNACL) { 101 PetscCall(PCSetUp_SAVIENNACL(pc)); 102 } 103 PetscCall(VecViennaCLGetArrayRead(x,&xarray)); 104 PetscCall(VecViennaCLGetArrayWrite(y,&yarray)); 105 try { 106 #if !defined(PETSC_USE_COMPLEX) 107 *yarray = *xarray; 108 sac->SAVIENNACL->apply(*yarray); 109 #endif 110 } catch(char * ex) { 111 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex); 112 } 113 PetscCall(VecViennaCLRestoreArrayRead(x,&xarray)); 114 PetscCall(VecViennaCLRestoreArrayWrite(y,&yarray)); 115 PetscCall(PetscObjectStateIncrease((PetscObject)y)); 116 PetscFunctionReturn(0); 117 } 118 /* -------------------------------------------------------------------------- */ 119 /* 120 PCDestroy_SAVIENNACL - Destroys the private context for the SAVIENNACL preconditioner 121 that was created with PCCreate_SAVIENNACL(). 122 123 Input Parameter: 124 . pc - the preconditioner context 125 126 Application Interface Routine: PCDestroy() 127 */ 128 static PetscErrorCode PCDestroy_SAVIENNACL(PC pc) 129 { 130 PC_SAVIENNACL *sac = (PC_SAVIENNACL*)pc->data; 131 132 PetscFunctionBegin; 133 if (sac->SAVIENNACL) { 134 try { 135 delete sac->SAVIENNACL; 136 } catch(char * ex) { 137 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex); 138 } 139 } 140 141 /* 142 Free the private data structure that was hanging off the PC 143 */ 144 PetscCall(PetscFree(pc->data)); 145 PetscFunctionReturn(0); 146 } 147 148 static PetscErrorCode PCSetFromOptions_SAVIENNACL(PetscOptionItems *PetscOptionsObject,PC pc) 149 { 150 PetscFunctionBegin; 151 PetscOptionsHeadBegin(PetscOptionsObject,"SAVIENNACL options"); 152 PetscOptionsHeadEnd(); 153 PetscFunctionReturn(0); 154 } 155 156 /* -------------------------------------------------------------------------- */ 157 158 /*MC 159 PCSAViennaCL - A smoothed agglomeration algorithm that can be used via the CUDA, OpenCL, and OpenMP backends of ViennaCL 160 161 Level: advanced 162 163 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC 164 165 M*/ 166 167 PETSC_EXTERN PetscErrorCode PCCreate_SAVIENNACL(PC pc) 168 { 169 PC_SAVIENNACL *sac; 170 171 PetscFunctionBegin; 172 /* 173 Creates the private data structure for this preconditioner and 174 attach it to the PC object. 175 */ 176 PetscCall(PetscNewLog(pc,&sac)); 177 pc->data = (void*)sac; 178 179 /* 180 Initialize the pointer to zero 181 Initialize number of v-cycles to default (1) 182 */ 183 sac->SAVIENNACL = 0; 184 185 /* 186 Set the pointers for the functions that are provided above. 187 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 188 are called, they will automatically call these functions. Note we 189 choose not to provide a couple of these functions since they are 190 not needed. 191 */ 192 pc->ops->apply = PCApply_SAVIENNACL; 193 pc->ops->applytranspose = 0; 194 pc->ops->setup = PCSetUp_SAVIENNACL; 195 pc->ops->destroy = PCDestroy_SAVIENNACL; 196 pc->ops->setfromoptions = PCSetFromOptions_SAVIENNACL; 197 pc->ops->view = 0; 198 pc->ops->applyrichardson = 0; 199 pc->ops->applysymmetricleft = 0; 200 pc->ops->applysymmetricright = 0; 201 PetscFunctionReturn(0); 202 } 203