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