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