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_SPINLOCK
7: #define PETSC_SKIP_IMMINTRIN_H_CUDAWORKAROUND 1
9: #include <petsc/private/pcimpl.h>
10: #include <../src/mat/impls/aij/seq/aij.h>
11: #include <../src/vec/vec/impls/dvecimpl.h>
12: #include <../src/mat/impls/aij/seq/seqviennacl/viennaclmatimpl.h>
13: #include <../src/vec/vec/impls/seq/seqviennacl/viennaclvecimpl.h>
14: #include <viennacl/linalg/detail/ilu/chow_patel_ilu.hpp>
16: /*
17: Private context (data structure) for the CHOWILUVIENNACL preconditioner.
18: */
19: typedef struct {
20: viennacl::linalg::chow_patel_ilu_precond<viennacl::compressed_matrix<PetscScalar>> *CHOWILUVIENNACL;
21: } PC_CHOWILUVIENNACL;
23: /*
24: PCSetUp_CHOWILUVIENNACL - Prepares for the use of the CHOWILUVIENNACL preconditioner
25: by setting data structures and options.
27: Input Parameter:
28: . pc - the preconditioner context
30: Application Interface Routine: PCSetUp()
32: Note:
33: The interface routine PCSetUp() is not usually called directly by
34: the user, but instead is called by PCApply() if necessary.
35: */
36: static PetscErrorCode PCSetUp_CHOWILUVIENNACL(PC pc)
37: {
38: PC_CHOWILUVIENNACL *ilu = (PC_CHOWILUVIENNACL *)pc->data;
39: PetscBool flg = PETSC_FALSE;
40: Mat_SeqAIJViennaCL *gpustruct;
42: PetscObjectTypeCompare((PetscObject)pc->pmat, MATSEQAIJVIENNACL, &flg);
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: 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: return 0;
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: /*how to apply a certain fixed number of iterations?*/
89: PetscObjectTypeCompare((PetscObject)x, VECSEQVIENNACL, &flg1);
90: PetscObjectTypeCompare((PetscObject)y, VECSEQVIENNACL, &flg2);
92: if (!ilu->CHOWILUVIENNACL) PCSetUp_CHOWILUVIENNACL(pc);
93: VecSet(y, 0.0);
94: VecViennaCLGetArrayRead(x, &xarray);
95: 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: VecViennaCLRestoreArrayRead(x, &xarray);
107: VecViennaCLRestoreArrayWrite(y, &yarray);
108: PetscObjectStateIncrease((PetscObject)y);
109: return 0;
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: if (ilu->CHOWILUVIENNACL) {
126: try {
127: delete ilu->CHOWILUVIENNACL;
128: } catch (char *ex) {
129: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex);
130: }
131: }
133: /*
134: Free the private data structure that was hanging off the PC
135: */
136: PetscFree(pc->data);
137: return 0;
138: }
140: static PetscErrorCode PCSetFromOptions_CHOWILUVIENNACL(PC pc, PetscOptionItems *PetscOptionsObject)
141: {
142: PetscOptionsHeadBegin(PetscOptionsObject, "CHOWILUVIENNACL options");
143: PetscOptionsHeadEnd();
144: return 0;
145: }
147: /*MC
148: PCCHOWILUViennaCL - A smoothed agglomeration algorithm that can be used via the CUDA, OpenCL, and OpenMP backends of ViennaCL
150: Level: advanced
152: Developer Note:
153: This does not appear to be wired up with `PCRegisterType()`
155: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`
156: M*/
158: PETSC_EXTERN PetscErrorCode PCCreate_CHOWILUVIENNACL(PC pc)
159: {
160: PC_CHOWILUVIENNACL *ilu;
162: /*
163: Creates the private data structure for this preconditioner and
164: attach it to the PC object.
165: */
166: PetscNew(&ilu);
167: pc->data = (void *)ilu;
169: /*
170: Initialize the pointer to zero
171: Initialize number of v-cycles to default (1)
172: */
173: ilu->CHOWILUVIENNACL = 0;
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_CHOWILUVIENNACL;
183: pc->ops->applytranspose = 0;
184: pc->ops->setup = PCSetUp_CHOWILUVIENNACL;
185: pc->ops->destroy = PCDestroy_CHOWILUVIENNACL;
186: pc->ops->setfromoptions = PCSetFromOptions_CHOWILUVIENNACL;
187: pc->ops->view = 0;
188: pc->ops->applyrichardson = 0;
189: pc->ops->applysymmetricleft = 0;
190: pc->ops->applysymmetricright = 0;
191: return 0;
192: }