Actual source code: rowscalingviennacl.cxx
2: /* -------------------------------------------------------------------- */
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
11: #include <petsc/private/pcimpl.h>
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>
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;
25: /* -------------------------------------------------------------------------- */
26: /*
27: PCSetUp_ROWSCALINGVIENNACL - Prepares for the use of the ROWSCALINGVIENNACL preconditioner
28: by setting data structures and options.
30: Input Parameter:
31: . pc - the preconditioner context
33: Application Interface Routine: PCSetUp()
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: PetscErrorCode ierr;
44: Mat_SeqAIJViennaCL *gpustruct;
47: PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQAIJVIENNACL,&flg);
48: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Currently only handles ViennaCL matrices");
49: if (pc->setupcalled != 0) {
50: try {
51: delete rowscaling->ROWSCALINGVIENNACL;
52: } catch(char *ex) {
53: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
54: }
55: }
56: try {
57: #if defined(PETSC_USE_COMPLEX)
58: gpustruct = NULL;
59: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"No support for complex arithmetic in ROWSCALINGVIENNACL preconditioner");
60: #else
61: MatViennaCLCopyToGPU(pc->pmat);
62: gpustruct = (Mat_SeqAIJViennaCL*)(pc->pmat->spptr);
64: viennacl::linalg::row_scaling_tag pc_tag(1);
65: ViennaCLAIJMatrix *mat = (ViennaCLAIJMatrix*)gpustruct->mat;
66: rowscaling->ROWSCALINGVIENNACL = new viennacl::linalg::row_scaling<viennacl::compressed_matrix<PetscScalar> >(*mat, pc_tag);
67: #endif
68: } catch(char *ex) {
69: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
70: }
71: return(0);
72: }
74: /* -------------------------------------------------------------------------- */
75: /*
76: PCApply_ROWSCALINGVIENNACL - Applies the ROWSCALINGVIENNACL preconditioner to a vector.
78: Input Parameters:
79: . pc - the preconditioner context
80: . x - input vector
82: Output Parameter:
83: . y - output vector
85: Application Interface Routine: PCApply()
86: */
87: static PetscErrorCode PCApply_ROWSCALINGVIENNACL(PC pc,Vec x,Vec y)
88: {
89: PC_ROWSCALINGVIENNACL *ilu = (PC_ROWSCALINGVIENNACL*)pc->data;
90: PetscErrorCode ierr;
91: PetscBool flg1,flg2;
92: viennacl::vector<PetscScalar> const *xarray=NULL;
93: viennacl::vector<PetscScalar> *yarray=NULL;
96: /*how to apply a certain fixed number of iterations?*/
97: PetscObjectTypeCompare((PetscObject)x,VECSEQVIENNACL,&flg1);
98: PetscObjectTypeCompare((PetscObject)y,VECSEQVIENNACL,&flg2);
99: if (!(flg1 && flg2)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP, "Currently only handles ViennaCL vectors");
100: if (!ilu->ROWSCALINGVIENNACL) {
101: PCSetUp_ROWSCALINGVIENNACL(pc);
102: }
103: VecSet(y,0.0);
104: VecViennaCLGetArrayRead(x,&xarray);
105: VecViennaCLGetArrayWrite(y,&yarray);
106: try {
107: #if defined(PETSC_USE_COMPLEX)
109: #else
110: *yarray = *xarray;
111: ilu->ROWSCALINGVIENNACL->apply(*yarray);
112: #endif
113: } catch(char * ex) {
114: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
115: }
116: VecViennaCLRestoreArrayRead(x,&xarray);
117: VecViennaCLRestoreArrayWrite(y,&yarray);
118: PetscObjectStateIncrease((PetscObject)y);
119: return(0);
120: }
121: /* -------------------------------------------------------------------------- */
122: /*
123: PCDestroy_ROWSCALINGVIENNACL - Destroys the private context for the ROWSCALINGVIENNACL preconditioner
124: that was created with PCCreate_ROWSCALINGVIENNACL().
126: Input Parameter:
127: . pc - the preconditioner context
129: Application Interface Routine: PCDestroy()
130: */
131: static PetscErrorCode PCDestroy_ROWSCALINGVIENNACL(PC pc)
132: {
133: PC_ROWSCALINGVIENNACL *rowscaling = (PC_ROWSCALINGVIENNACL*)pc->data;
137: if (rowscaling->ROWSCALINGVIENNACL) {
138: try {
139: delete rowscaling->ROWSCALINGVIENNACL;
140: } catch(char *ex) {
141: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
142: }
143: }
145: /*
146: Free the private data structure that was hanging off the PC
147: */
148: PetscFree(pc->data);
149: return(0);
150: }
152: static PetscErrorCode PCSetFromOptions_ROWSCALINGVIENNACL(PetscOptionItems *PetscOptionsObject,PC pc)
153: {
157: PetscOptionsHead(PetscOptionsObject,"ROWSCALINGVIENNACL options");
158: PetscOptionsTail();
159: return(0);
160: }
162: /* -------------------------------------------------------------------------- */
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
167: Level: advanced
169: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC
171: M*/
173: PETSC_EXTERN PetscErrorCode PCCreate_ROWSCALINGVIENNACL(PC pc)
174: {
175: PC_ROWSCALINGVIENNACL *rowscaling;
179: /*
180: Creates the private data structure for this preconditioner and
181: attach it to the PC object.
182: */
183: PetscNewLog(pc,&rowscaling);
184: pc->data = (void*)rowscaling;
186: /*
187: Initialize the pointer to zero
188: Initialize number of v-cycles to default (1)
189: */
190: rowscaling->ROWSCALINGVIENNACL = 0;
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: return(0);
209: }