Actual source code: chowiluviennacl.cxx
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
7: #include <petsc/private/pcimpl.h>
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>
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;
21: /*
22: PCSetUp_CHOWILUVIENNACL - Prepares for the use of the CHOWILUVIENNACL preconditioner
23: by setting data structures and options.
25: Input Parameter:
26: . pc - the preconditioner context
28: Application Interface Routine: PCSetUp()
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: */
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;
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 != 0) {
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;
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: }
68: /*
69: PCApply_CHOWILUVIENNACL - Applies the CHOWILUVIENNACL preconditioner to a vector.
71: Input Parameters:
72: . pc - the preconditioner context
73: . x - input vector
75: Output Parameter:
76: . y - output vector
78: Application Interface Routine: PCApply()
79: */
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;
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)
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: }
112: /*
113: PCDestroy_CHOWILUVIENNACL - Destroys the private context for the CHOWILUVIENNACL preconditioner
114: that was created with PCCreate_CHOWILUVIENNACL().
116: Input Parameter:
117: . pc - the preconditioner context
119: Application Interface Routine: PCDestroy()
120: */
121: static PetscErrorCode PCDestroy_CHOWILUVIENNACL(PC pc)
122: {
123: PC_CHOWILUVIENNACL *ilu = (PC_CHOWILUVIENNACL *)pc->data;
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: }
134: /*
135: Free the private data structure that was hanging off the PC
136: */
137: PetscCall(PetscFree(pc->data));
138: PetscFunctionReturn(PETSC_SUCCESS);
139: }
141: static PetscErrorCode PCSetFromOptions_CHOWILUVIENNACL(PC pc, PetscOptionItems *PetscOptionsObject)
142: {
143: PetscFunctionBegin;
144: PetscOptionsHeadBegin(PetscOptionsObject, "CHOWILUVIENNACL options");
145: PetscOptionsHeadEnd();
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: PETSC_EXTERN PetscErrorCode PCCreate_CHOWILUVIENNACL(PC pc)
150: {
151: PC_CHOWILUVIENNACL *ilu;
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;
161: /*
162: Initialize the pointer to zero
163: Initialize number of v-cycles to default (1)
164: */
165: ilu->CHOWILUVIENNACL = 0;
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: }