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 */
PCSetUp_ROWSCALINGVIENNACL(PC pc)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 */
PCApply_ROWSCALINGVIENNACL(PC pc,Vec x,Vec y)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 */
PCDestroy_ROWSCALINGVIENNACL(PC pc)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
PCSetFromOptions_ROWSCALINGVIENNACL(PC pc,PetscOptionItems PetscOptionsObject)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*/
PCCreate_ROWSCALINGVIENNACL(PC pc)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