1 /*
2 Include files needed for the ViennaCL Chow-Patel parallel ILU 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/detail/ilu/chow_patel_ilu.hpp>
13
14 /*
15 Private context (data structure) for the CHOWILUVIENNACL preconditioner.
16 */
17 typedef struct {
18 viennacl::linalg::chow_patel_ilu_precond<viennacl::compressed_matrix<PetscScalar>> *CHOWILUVIENNACL;
19 } PC_CHOWILUVIENNACL;
20
21 /*
22 PCSetUp_CHOWILUVIENNACL - Prepares for the use of the CHOWILUVIENNACL 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_CHOWILUVIENNACL(PC pc)34 static PetscErrorCode PCSetUp_CHOWILUVIENNACL(PC pc)
35 {
36 PC_CHOWILUVIENNACL *ilu = (PC_CHOWILUVIENNACL *)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 ilu->CHOWILUVIENNACL;
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 CHOWILUVIENNACL preconditioner");
54 #else
55 PetscCall(MatViennaCLCopyToGPU(pc->pmat));
56 gpustruct = (Mat_SeqAIJViennaCL *)pc->pmat->spptr;
57
58 viennacl::linalg::chow_patel_tag ilu_tag;
59 ViennaCLAIJMatrix *mat = (ViennaCLAIJMatrix *)gpustruct->mat;
60 ilu->CHOWILUVIENNACL = new viennacl::linalg::chow_patel_ilu_precond<viennacl::compressed_matrix<PetscScalar>>(*mat, ilu_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_CHOWILUVIENNACL - Applies the CHOWILUVIENNACL 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_CHOWILUVIENNACL(PC pc,Vec x,Vec y)80 static PetscErrorCode PCApply_CHOWILUVIENNACL(PC pc, Vec x, Vec y)
81 {
82 PC_CHOWILUVIENNACL *ilu = (PC_CHOWILUVIENNACL *)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->CHOWILUVIENNACL) PetscCall(PCSetUp_CHOWILUVIENNACL(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->CHOWILUVIENNACL->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_CHOWILUVIENNACL - Destroys the private context for the CHOWILUVIENNACL preconditioner
114 that was created with PCCreate_CHOWILUVIENNACL().
115
116 Input Parameter:
117 . pc - the preconditioner context
118
119 Application Interface Routine: PCDestroy()
120 */
PCDestroy_CHOWILUVIENNACL(PC pc)121 static PetscErrorCode PCDestroy_CHOWILUVIENNACL(PC pc)
122 {
123 PC_CHOWILUVIENNACL *ilu = (PC_CHOWILUVIENNACL *)pc->data;
124
125 PetscFunctionBegin;
126 if (ilu->CHOWILUVIENNACL) {
127 try {
128 delete ilu->CHOWILUVIENNACL;
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_CHOWILUVIENNACL(PC pc,PetscOptionItems PetscOptionsObject)141 static PetscErrorCode PCSetFromOptions_CHOWILUVIENNACL(PC pc, PetscOptionItems PetscOptionsObject)
142 {
143 PetscFunctionBegin;
144 PetscOptionsHeadBegin(PetscOptionsObject, "CHOWILUVIENNACL options");
145 PetscOptionsHeadEnd();
146 PetscFunctionReturn(PETSC_SUCCESS);
147 }
148
PCCreate_CHOWILUVIENNACL(PC pc)149 PETSC_EXTERN PetscErrorCode PCCreate_CHOWILUVIENNACL(PC pc)
150 {
151 PC_CHOWILUVIENNACL *ilu;
152
153 PetscFunctionBegin;
154 /*
155 Creates the private data structure for this preconditioner and
156 attach it to the PC object.
157 */
158 PetscCall(PetscNew(&ilu));
159 pc->data = (void *)ilu;
160
161 /*
162 Initialize the pointer to zero
163 Initialize number of v-cycles to default (1)
164 */
165 ilu->CHOWILUVIENNACL = 0;
166
167 /*
168 Set the pointers for the functions that are provided above.
169 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
170 are called, they will automatically call these functions. Note we
171 choose not to provide a couple of these functions since they are
172 not needed.
173 */
174 pc->ops->apply = PCApply_CHOWILUVIENNACL;
175 pc->ops->applytranspose = 0;
176 pc->ops->setup = PCSetUp_CHOWILUVIENNACL;
177 pc->ops->destroy = PCDestroy_CHOWILUVIENNACL;
178 pc->ops->setfromoptions = PCSetFromOptions_CHOWILUVIENNACL;
179 pc->ops->view = 0;
180 pc->ops->applyrichardson = 0;
181 pc->ops->applysymmetricleft = 0;
182 pc->ops->applysymmetricright = 0;
183 PetscFunctionReturn(PETSC_SUCCESS);
184 }
185