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