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