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 #undef __FUNCT__ 39 #define __FUNCT__ "PCSetUp_ROWSCALINGVIENNACL" 40 static PetscErrorCode PCSetUp_ROWSCALINGVIENNACL(PC pc) 41 { 42 PC_ROWSCALINGVIENNACL *rowscaling = (PC_ROWSCALINGVIENNACL*)pc->data; 43 PetscBool flg = PETSC_FALSE; 44 PetscErrorCode ierr; 45 Mat_SeqAIJViennaCL *gpustruct; 46 47 PetscFunctionBegin; 48 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQAIJVIENNACL,&flg);CHKERRQ(ierr); 49 if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Currently only handles ViennaCL matrices"); 50 if (pc->setupcalled != 0) { 51 try { 52 delete rowscaling->ROWSCALINGVIENNACL; 53 } catch(char *ex) { 54 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex); 55 } 56 } 57 try { 58 #if defined(PETSC_USE_COMPLEX) 59 gpustruct = NULL; 60 SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"No support for complex arithmetic in ROWSCALINGVIENNACL preconditioner"); 61 #else 62 ierr = MatViennaCLCopyToGPU(pc->pmat);CHKERRQ(ierr); 63 gpustruct = (Mat_SeqAIJViennaCL*)(pc->pmat->spptr); 64 65 viennacl::linalg::row_scaling_tag pc_tag(1); 66 ViennaCLAIJMatrix *mat = (ViennaCLAIJMatrix*)gpustruct->mat; 67 rowscaling->ROWSCALINGVIENNACL = new viennacl::linalg::row_scaling<viennacl::compressed_matrix<PetscScalar> >(*mat, pc_tag); 68 #endif 69 } catch(char *ex) { 70 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex); 71 } 72 PetscFunctionReturn(0); 73 } 74 75 /* -------------------------------------------------------------------------- */ 76 /* 77 PCApply_ROWSCALINGVIENNACL - Applies the ROWSCALINGVIENNACL 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 #undef __FUNCT__ 89 #define __FUNCT__ "PCApply_ROWSCALINGVIENNACL" 90 static PetscErrorCode PCApply_ROWSCALINGVIENNACL(PC pc,Vec x,Vec y) 91 { 92 PC_ROWSCALINGVIENNACL *ilu = (PC_ROWSCALINGVIENNACL*)pc->data; 93 PetscErrorCode ierr; 94 PetscBool flg1,flg2; 95 viennacl::vector<PetscScalar> const *xarray=NULL; 96 viennacl::vector<PetscScalar> *yarray=NULL; 97 98 PetscFunctionBegin; 99 /*how to apply a certain fixed number of iterations?*/ 100 ierr = PetscObjectTypeCompare((PetscObject)x,VECSEQVIENNACL,&flg1);CHKERRQ(ierr); 101 ierr = PetscObjectTypeCompare((PetscObject)y,VECSEQVIENNACL,&flg2);CHKERRQ(ierr); 102 if (!(flg1 && flg2)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP, "Currently only handles ViennaCL vectors"); 103 if (!ilu->ROWSCALINGVIENNACL) { 104 ierr = PCSetUp_ROWSCALINGVIENNACL(pc);CHKERRQ(ierr); 105 } 106 ierr = VecSet(y,0.0);CHKERRQ(ierr); 107 ierr = VecViennaCLGetArrayRead(x,&xarray);CHKERRQ(ierr); 108 ierr = VecViennaCLGetArrayWrite(y,&yarray);CHKERRQ(ierr); 109 try { 110 #if defined(PETSC_USE_COMPLEX) 111 112 #else 113 *yarray = *xarray; 114 ilu->ROWSCALINGVIENNACL->apply(*yarray); 115 #endif 116 } catch(char * ex) { 117 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex); 118 } 119 ierr = VecViennaCLRestoreArrayRead(x,&xarray);CHKERRQ(ierr); 120 ierr = VecViennaCLRestoreArrayWrite(y,&yarray);CHKERRQ(ierr); 121 ierr = PetscObjectStateIncrease((PetscObject)y);CHKERRQ(ierr); 122 PetscFunctionReturn(0); 123 } 124 /* -------------------------------------------------------------------------- */ 125 /* 126 PCDestroy_ROWSCALINGVIENNACL - Destroys the private context for the ROWSCALINGVIENNACL preconditioner 127 that was created with PCCreate_ROWSCALINGVIENNACL(). 128 129 Input Parameter: 130 . pc - the preconditioner context 131 132 Application Interface Routine: PCDestroy() 133 */ 134 #undef __FUNCT__ 135 #define __FUNCT__ "PCDestroy_ROWSCALINGVIENNACL" 136 static PetscErrorCode PCDestroy_ROWSCALINGVIENNACL(PC pc) 137 { 138 PC_ROWSCALINGVIENNACL *rowscaling = (PC_ROWSCALINGVIENNACL*)pc->data; 139 PetscErrorCode ierr; 140 141 PetscFunctionBegin; 142 if (rowscaling->ROWSCALINGVIENNACL) { 143 try { 144 delete rowscaling->ROWSCALINGVIENNACL; 145 } catch(char *ex) { 146 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex); 147 } 148 } 149 150 /* 151 Free the private data structure that was hanging off the PC 152 */ 153 ierr = PetscFree(pc->data);CHKERRQ(ierr); 154 PetscFunctionReturn(0); 155 } 156 157 #undef __FUNCT__ 158 #define __FUNCT__ "PCSetFromOptions_ROWSCALINGVIENNACL" 159 static PetscErrorCode PCSetFromOptions_ROWSCALINGVIENNACL(PetscOptionItems *PetscOptionsObject,PC pc) 160 { 161 PetscErrorCode ierr; 162 163 PetscFunctionBegin; 164 ierr = PetscOptionsHead(PetscOptionsObject,"ROWSCALINGVIENNACL options");CHKERRQ(ierr); 165 ierr = PetscOptionsTail();CHKERRQ(ierr); 166 PetscFunctionReturn(0); 167 } 168 169 /* -------------------------------------------------------------------------- */ 170 171 172 /*MC 173 PCRowScalingViennaCL - A diagonal preconditioner (scaling rows of matrices by their norm) that can be used via the CUDA, OpenCL, and OpenMP backends of ViennaCL 174 175 Level: advanced 176 177 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC 178 179 M*/ 180 181 #undef __FUNCT__ 182 #define __FUNCT__ "PCCreate_ROWSCALINGVIENNACL" 183 PETSC_EXTERN PetscErrorCode PCCreate_ROWSCALINGVIENNACL(PC pc) 184 { 185 PC_ROWSCALINGVIENNACL *rowscaling; 186 PetscErrorCode ierr; 187 188 PetscFunctionBegin; 189 /* 190 Creates the private data structure for this preconditioner and 191 attach it to the PC object. 192 */ 193 ierr = PetscNewLog(pc,&rowscaling);CHKERRQ(ierr); 194 pc->data = (void*)rowscaling; 195 196 /* 197 Initialize the pointer to zero 198 Initialize number of v-cycles to default (1) 199 */ 200 rowscaling->ROWSCALINGVIENNACL = 0; 201 202 /* 203 Set the pointers for the functions that are provided above. 204 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 205 are called, they will automatically call these functions. Note we 206 choose not to provide a couple of these functions since they are 207 not needed. 208 */ 209 pc->ops->apply = PCApply_ROWSCALINGVIENNACL; 210 pc->ops->applytranspose = 0; 211 pc->ops->setup = PCSetUp_ROWSCALINGVIENNACL; 212 pc->ops->destroy = PCDestroy_ROWSCALINGVIENNACL; 213 pc->ops->setfromoptions = PCSetFromOptions_ROWSCALINGVIENNACL; 214 pc->ops->view = 0; 215 pc->ops->applyrichardson = 0; 216 pc->ops->applysymmetricleft = 0; 217 pc->ops->applysymmetricright = 0; 218 PetscFunctionReturn(0); 219 } 220 221