Actual source code: saviennacl.cxx
2: /*
3: Include files needed for the ViennaCL Smoothed Aggregation 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
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/amg.hpp>
15: /*
16: Private context (data structure) for the SAVIENNACL preconditioner.
17: */
18: typedef struct {
19: viennacl::linalg::amg_precond<viennacl::compressed_matrix<PetscScalar>> *SAVIENNACL;
20: } PC_SAVIENNACL;
22: /*
23: PCSetUp_SAVIENNACL - Prepares for the use of the SAVIENNACL 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_SAVIENNACL(PC pc)
36: {
37: PC_SAVIENNACL *sa = (PC_SAVIENNACL *)pc->data;
38: PetscBool flg = PETSC_FALSE;
39: Mat_SeqAIJViennaCL *gpustruct;
41: PetscObjectTypeCompare((PetscObject)pc->pmat, MATSEQAIJVIENNACL, &flg);
43: if (pc->setupcalled != 0) {
44: try {
45: delete sa->SAVIENNACL;
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 SAVIENNACL preconditioner");
54: #else
55: MatViennaCLCopyToGPU(pc->pmat);
56: gpustruct = (Mat_SeqAIJViennaCL *)(pc->pmat->spptr);
58: viennacl::linalg::amg_tag amg_tag_sa_pmis;
59: amg_tag_sa_pmis.set_coarsening_method(viennacl::linalg::AMG_COARSENING_METHOD_MIS2_AGGREGATION);
60: amg_tag_sa_pmis.set_interpolation_method(viennacl::linalg::AMG_INTERPOLATION_METHOD_SMOOTHED_AGGREGATION);
61: ViennaCLAIJMatrix *mat = (ViennaCLAIJMatrix *)gpustruct->mat;
62: sa->SAVIENNACL = new viennacl::linalg::amg_precond<viennacl::compressed_matrix<PetscScalar>>(*mat, amg_tag_sa_pmis);
63: sa->SAVIENNACL->setup();
64: #endif
65: } catch (char *ex) {
66: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex);
67: }
68: return 0;
69: }
71: /*
72: PCApply_SAVIENNACL - Applies the SAVIENNACL preconditioner to a vector.
74: Input Parameters:
75: . pc - the preconditioner context
76: . x - input vector
78: Output Parameter:
79: . y - output vector
81: Application Interface Routine: PCApply()
82: */
83: static PetscErrorCode PCApply_SAVIENNACL(PC pc, Vec x, Vec y)
84: {
85: PC_SAVIENNACL *sac = (PC_SAVIENNACL *)pc->data;
86: PetscBool flg1, flg2;
87: viennacl::vector<PetscScalar> const *xarray = NULL;
88: viennacl::vector<PetscScalar> *yarray = NULL;
90: /*how to apply a certain fixed number of iterations?*/
91: PetscObjectTypeCompare((PetscObject)x, VECSEQVIENNACL, &flg1);
92: PetscObjectTypeCompare((PetscObject)y, VECSEQVIENNACL, &flg2);
94: if (!sac->SAVIENNACL) PCSetUp_SAVIENNACL(pc);
95: VecViennaCLGetArrayRead(x, &xarray);
96: 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: VecViennaCLRestoreArrayRead(x, &xarray);
106: VecViennaCLRestoreArrayWrite(y, &yarray);
107: PetscObjectStateIncrease((PetscObject)y);
108: return 0;
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: if (sac->SAVIENNACL) {
125: try {
126: delete sac->SAVIENNACL;
127: } catch (char *ex) {
128: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex);
129: }
130: }
132: /*
133: Free the private data structure that was hanging off the PC
134: */
135: PetscFree(pc->data);
136: return 0;
137: }
139: static PetscErrorCode PCSetFromOptions_SAVIENNACL(PC pc, PetscOptionItems *PetscOptionsObject)
140: {
141: PetscOptionsHeadBegin(PetscOptionsObject, "SAVIENNACL options");
142: PetscOptionsHeadEnd();
143: return 0;
144: }
146: /*MC
147: PCSAViennaCL - A smoothed agglomeration algorithm that can be used via the CUDA, OpenCL, and OpenMP backends of ViennaCL
149: Level: advanced
151: Developer Note:
152: This `PCType` does not appear to be registered
154: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`
155: M*/
157: PETSC_EXTERN PetscErrorCode PCCreate_SAVIENNACL(PC pc)
158: {
159: PC_SAVIENNACL *sac;
161: /*
162: Creates the private data structure for this preconditioner and
163: attach it to the PC object.
164: */
165: PetscNew(&sac);
166: pc->data = (void *)sac;
168: /*
169: Initialize the pointer to zero
170: Initialize number of v-cycles to default (1)
171: */
172: sac->SAVIENNACL = 0;
174: /*
175: Set the pointers for the functions that are provided above.
176: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
177: are called, they will automatically call these functions. Note we
178: choose not to provide a couple of these functions since they are
179: not needed.
180: */
181: pc->ops->apply = PCApply_SAVIENNACL;
182: pc->ops->applytranspose = 0;
183: pc->ops->setup = PCSetUp_SAVIENNACL;
184: pc->ops->destroy = PCDestroy_SAVIENNACL;
185: pc->ops->setfromoptions = PCSetFromOptions_SAVIENNACL;
186: pc->ops->view = 0;
187: pc->ops->applyrichardson = 0;
188: pc->ops->applysymmetricleft = 0;
189: pc->ops->applysymmetricright = 0;
190: return 0;
191: }