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 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 10 #include <../src/mat/impls/aij/seq/aij.h> 11 #include <../src/vec/vec/impls/dvecimpl.h> 12 #include <../src/mat/impls/aij/seq/seqviennacl/viennaclmatimpl.h> 13 #include <../src/vec/vec/impls/seq/seqviennacl/viennaclvecimpl.h> 14 #include <viennacl/linalg/row_scaling.hpp> 15 16 /* 17 Private context (data structure) for the ROWSCALINGVIENNACL preconditioner. 18 */ 19 typedef struct { 20 viennacl::linalg::row_scaling< viennacl::compressed_matrix<PetscScalar> > *ROWSCALINGVIENNACL; 21 } PC_ROWSCALINGVIENNACL; 22 23 24 /* -------------------------------------------------------------------------- */ 25 /* 26 PCSetUp_ROWSCALINGVIENNACL - Prepares for the use of the ROWSCALINGVIENNACL 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_ROWSCALINGVIENNACL(PC pc) 39 { 40 PC_ROWSCALINGVIENNACL *rowscaling = (PC_ROWSCALINGVIENNACL*)pc->data; 41 PetscBool flg = PETSC_FALSE; 42 PetscErrorCode ierr; 43 Mat_SeqAIJViennaCL *gpustruct; 44 45 PetscFunctionBegin; 46 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQAIJVIENNACL,&flg);CHKERRQ(ierr); 47 if (!flg) SETERRQ(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 SETERRQ1(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 ierr = MatViennaCLCopyToGPU(pc->pmat);CHKERRQ(ierr); 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 SETERRQ1(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 PetscErrorCode ierr; 90 PetscBool flg1,flg2; 91 viennacl::vector<PetscScalar> const *xarray=NULL; 92 viennacl::vector<PetscScalar> *yarray=NULL; 93 94 PetscFunctionBegin; 95 /*how to apply a certain fixed number of iterations?*/ 96 ierr = PetscObjectTypeCompare((PetscObject)x,VECSEQVIENNACL,&flg1);CHKERRQ(ierr); 97 ierr = PetscObjectTypeCompare((PetscObject)y,VECSEQVIENNACL,&flg2);CHKERRQ(ierr); 98 if (!(flg1 && flg2)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP, "Currently only handles ViennaCL vectors"); 99 if (!ilu->ROWSCALINGVIENNACL) { 100 ierr = PCSetUp_ROWSCALINGVIENNACL(pc);CHKERRQ(ierr); 101 } 102 ierr = VecSet(y,0.0);CHKERRQ(ierr); 103 ierr = VecViennaCLGetArrayRead(x,&xarray);CHKERRQ(ierr); 104 ierr = VecViennaCLGetArrayWrite(y,&yarray);CHKERRQ(ierr); 105 try { 106 #if defined(PETSC_USE_COMPLEX) 107 108 #else 109 *yarray = *xarray; 110 ilu->ROWSCALINGVIENNACL->apply(*yarray); 111 #endif 112 } catch(char * ex) { 113 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex); 114 } 115 ierr = VecViennaCLRestoreArrayRead(x,&xarray);CHKERRQ(ierr); 116 ierr = VecViennaCLRestoreArrayWrite(y,&yarray);CHKERRQ(ierr); 117 ierr = PetscObjectStateIncrease((PetscObject)y);CHKERRQ(ierr); 118 PetscFunctionReturn(0); 119 } 120 /* -------------------------------------------------------------------------- */ 121 /* 122 PCDestroy_ROWSCALINGVIENNACL - Destroys the private context for the ROWSCALINGVIENNACL preconditioner 123 that was created with PCCreate_ROWSCALINGVIENNACL(). 124 125 Input Parameter: 126 . pc - the preconditioner context 127 128 Application Interface Routine: PCDestroy() 129 */ 130 static PetscErrorCode PCDestroy_ROWSCALINGVIENNACL(PC pc) 131 { 132 PC_ROWSCALINGVIENNACL *rowscaling = (PC_ROWSCALINGVIENNACL*)pc->data; 133 PetscErrorCode ierr; 134 135 PetscFunctionBegin; 136 if (rowscaling->ROWSCALINGVIENNACL) { 137 try { 138 delete rowscaling->ROWSCALINGVIENNACL; 139 } catch(char *ex) { 140 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex); 141 } 142 } 143 144 /* 145 Free the private data structure that was hanging off the PC 146 */ 147 ierr = PetscFree(pc->data);CHKERRQ(ierr); 148 PetscFunctionReturn(0); 149 } 150 151 static PetscErrorCode PCSetFromOptions_ROWSCALINGVIENNACL(PetscOptionItems *PetscOptionsObject,PC pc) 152 { 153 PetscErrorCode ierr; 154 155 PetscFunctionBegin; 156 ierr = PetscOptionsHead(PetscOptionsObject,"ROWSCALINGVIENNACL options");CHKERRQ(ierr); 157 ierr = PetscOptionsTail();CHKERRQ(ierr); 158 PetscFunctionReturn(0); 159 } 160 161 /* -------------------------------------------------------------------------- */ 162 163 164 /*MC 165 PCRowScalingViennaCL - A diagonal preconditioner (scaling rows of matrices by their norm) that can be used via the CUDA, OpenCL, and OpenMP backends of ViennaCL 166 167 Level: advanced 168 169 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC 170 171 M*/ 172 173 PETSC_EXTERN PetscErrorCode PCCreate_ROWSCALINGVIENNACL(PC pc) 174 { 175 PC_ROWSCALINGVIENNACL *rowscaling; 176 PetscErrorCode ierr; 177 178 PetscFunctionBegin; 179 /* 180 Creates the private data structure for this preconditioner and 181 attach it to the PC object. 182 */ 183 ierr = PetscNewLog(pc,&rowscaling);CHKERRQ(ierr); 184 pc->data = (void*)rowscaling; 185 186 /* 187 Initialize the pointer to zero 188 Initialize number of v-cycles to default (1) 189 */ 190 rowscaling->ROWSCALINGVIENNACL = 0; 191 192 /* 193 Set the pointers for the functions that are provided above. 194 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 195 are called, they will automatically call these functions. Note we 196 choose not to provide a couple of these functions since they are 197 not needed. 198 */ 199 pc->ops->apply = PCApply_ROWSCALINGVIENNACL; 200 pc->ops->applytranspose = 0; 201 pc->ops->setup = PCSetUp_ROWSCALINGVIENNACL; 202 pc->ops->destroy = PCDestroy_ROWSCALINGVIENNACL; 203 pc->ops->setfromoptions = PCSetFromOptions_ROWSCALINGVIENNACL; 204 pc->ops->view = 0; 205 pc->ops->applyrichardson = 0; 206 pc->ops->applysymmetricleft = 0; 207 pc->ops->applysymmetricright = 0; 208 PetscFunctionReturn(0); 209 } 210 211