Actual source code: saviennacl.cxx
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>
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>
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;
20: /*
21: PCSetUp_SAVIENNACL - Prepares for the use of the SAVIENNACL preconditioner
22: by setting data structures and options.
24: Input Parameter:
25: . pc - the preconditioner context
27: Application Interface Routine: PCSetUp()
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: */
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;
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 != 0) {
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;
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: }
70: /*
71: PCApply_SAVIENNACL - Applies the SAVIENNACL preconditioner to a vector.
73: Input Parameters:
74: . pc - the preconditioner context
75: . x - input vector
77: Output Parameter:
78: . y - output vector
80: Application Interface Routine: PCApply()
81: */
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;
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: }
111: /*
112: PCDestroy_SAVIENNACL - Destroys the private context for the SAVIENNACL preconditioner
113: that was created with PCCreate_SAVIENNACL().
115: Input Parameter:
116: . pc - the preconditioner context
118: Application Interface Routine: PCDestroy()
119: */
120: static PetscErrorCode PCDestroy_SAVIENNACL(PC pc)
121: {
122: PC_SAVIENNACL *sac = (PC_SAVIENNACL *)pc->data;
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: }
133: /*
134: Free the private data structure that was hanging off the PC
135: */
136: PetscCall(PetscFree(pc->data));
137: PetscFunctionReturn(PETSC_SUCCESS);
138: }
140: static PetscErrorCode PCSetFromOptions_SAVIENNACL(PC pc, PetscOptionItems *PetscOptionsObject)
141: {
142: PetscFunctionBegin;
143: PetscOptionsHeadBegin(PetscOptionsObject, "SAVIENNACL options");
144: PetscOptionsHeadEnd();
145: PetscFunctionReturn(PETSC_SUCCESS);
146: }
148: /*MC
149: PCSAVIENNACL - A smoothed agglomeration algorithm that can be used via the CUDA, OpenCL, and OpenMP backends of ViennaCL
151: Level: advanced
153: Developer Notes:
154: This `PCType` does not appear to be registered
156: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`
157: M*/
159: PETSC_EXTERN PetscErrorCode PCCreate_SAVIENNACL(PC pc)
160: {
161: PC_SAVIENNACL *sac;
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;
171: /*
172: Initialize the pointer to zero
173: Initialize number of v-cycles to default (1)
174: */
175: sac->SAVIENNACL = 0;
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: }