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