Actual source code: saviennacl.cxx
2: /* -------------------------------------------------------------------- */
4: /*
5: Include files needed for the ViennaCL Smoothed Aggregation preconditioner:
6: pcimpl.h - private include file intended for use by all preconditioners
7: */
8: #define PETSC_SKIP_SPINLOCK
9: #define PETSC_SKIP_IMMINTRIN_H_CUDAWORKAROUND 1
10: #include <petsc/private/pcimpl.h>
11: #include <../src/mat/impls/aij/seq/aij.h>
12: #include <../src/vec/vec/impls/dvecimpl.h>
13: #include <../src/mat/impls/aij/seq/seqviennacl/viennaclmatimpl.h>
14: #include <../src/vec/vec/impls/seq/seqviennacl/viennaclvecimpl.h>
15: #include <viennacl/linalg/amg.hpp>
17: /*
18: Private context (data structure) for the SAVIENNACL preconditioner.
19: */
20: typedef struct {
21: viennacl::linalg::amg_precond<viennacl::compressed_matrix<PetscScalar> > *SAVIENNACL;
22: } PC_SAVIENNACL;
24: /* -------------------------------------------------------------------------- */
25: /*
26: PCSetUp_SAVIENNACL - Prepares for the use of the SAVIENNACL preconditioner
27: by setting data structures and options.
29: Input Parameter:
30: . pc - the preconditioner context
32: Application Interface Routine: PCSetUp()
34: Notes:
35: The interface routine PCSetUp() is not usually called directly by
36: the user, but instead is called by PCApply() if necessary.
37: */
38: static PetscErrorCode PCSetUp_SAVIENNACL(PC pc)
39: {
40: PC_SAVIENNACL *sa = (PC_SAVIENNACL*)pc->data;
41: PetscBool flg = PETSC_FALSE;
42: Mat_SeqAIJViennaCL *gpustruct;
44: PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQAIJVIENNACL,&flg);
46: if (pc->setupcalled != 0) {
47: try {
48: delete sa->SAVIENNACL;
49: } catch(char *ex) {
50: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
51: }
52: }
53: try {
54: #if defined(PETSC_USE_COMPLEX)
55: gpustruct = NULL;
56: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"No support for complex arithmetic in SAVIENNACL preconditioner");
57: #else
58: MatViennaCLCopyToGPU(pc->pmat);
59: gpustruct = (Mat_SeqAIJViennaCL*)(pc->pmat->spptr);
61: viennacl::linalg::amg_tag amg_tag_sa_pmis;
62: amg_tag_sa_pmis.set_coarsening_method(viennacl::linalg::AMG_COARSENING_METHOD_MIS2_AGGREGATION);
63: amg_tag_sa_pmis.set_interpolation_method(viennacl::linalg::AMG_INTERPOLATION_METHOD_SMOOTHED_AGGREGATION);
64: ViennaCLAIJMatrix *mat = (ViennaCLAIJMatrix*)gpustruct->mat;
65: sa->SAVIENNACL = new viennacl::linalg::amg_precond<viennacl::compressed_matrix<PetscScalar> >(*mat, amg_tag_sa_pmis);
66: sa->SAVIENNACL->setup();
67: #endif
68: } catch(char *ex) {
69: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
70: }
71: return 0;
72: }
74: /* -------------------------------------------------------------------------- */
75: /*
76: PCApply_SAVIENNACL - Applies the SAVIENNACL preconditioner to a vector.
78: Input Parameters:
79: . pc - the preconditioner context
80: . x - input vector
82: Output Parameter:
83: . y - output vector
85: Application Interface Routine: PCApply()
86: */
87: static PetscErrorCode PCApply_SAVIENNACL(PC pc,Vec x,Vec y)
88: {
89: PC_SAVIENNACL *sac = (PC_SAVIENNACL*)pc->data;
90: PetscBool flg1,flg2;
91: viennacl::vector<PetscScalar> const *xarray=NULL;
92: viennacl::vector<PetscScalar> *yarray=NULL;
94: /*how to apply a certain fixed number of iterations?*/
95: PetscObjectTypeCompare((PetscObject)x,VECSEQVIENNACL,&flg1);
96: PetscObjectTypeCompare((PetscObject)y,VECSEQVIENNACL,&flg2);
98: if (!sac->SAVIENNACL) {
99: PCSetUp_SAVIENNACL(pc);
100: }
101: VecViennaCLGetArrayRead(x,&xarray);
102: VecViennaCLGetArrayWrite(y,&yarray);
103: try {
104: #if !defined(PETSC_USE_COMPLEX)
105: *yarray = *xarray;
106: sac->SAVIENNACL->apply(*yarray);
107: #endif
108: } catch(char * ex) {
109: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
110: }
111: VecViennaCLRestoreArrayRead(x,&xarray);
112: VecViennaCLRestoreArrayWrite(y,&yarray);
113: PetscObjectStateIncrease((PetscObject)y);
114: return 0;
115: }
116: /* -------------------------------------------------------------------------- */
117: /*
118: PCDestroy_SAVIENNACL - Destroys the private context for the SAVIENNACL preconditioner
119: that was created with PCCreate_SAVIENNACL().
121: Input Parameter:
122: . pc - the preconditioner context
124: Application Interface Routine: PCDestroy()
125: */
126: static PetscErrorCode PCDestroy_SAVIENNACL(PC pc)
127: {
128: PC_SAVIENNACL *sac = (PC_SAVIENNACL*)pc->data;
130: if (sac->SAVIENNACL) {
131: try {
132: delete sac->SAVIENNACL;
133: } catch(char * ex) {
134: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
135: }
136: }
138: /*
139: Free the private data structure that was hanging off the PC
140: */
141: PetscFree(pc->data);
142: return 0;
143: }
145: static PetscErrorCode PCSetFromOptions_SAVIENNACL(PetscOptionItems *PetscOptionsObject,PC pc)
146: {
147: PetscOptionsHead(PetscOptionsObject,"SAVIENNACL options");
148: PetscOptionsTail();
149: return 0;
150: }
152: /* -------------------------------------------------------------------------- */
154: /*MC
155: PCSAViennaCL - A smoothed agglomeration algorithm that can be used via the CUDA, OpenCL, and OpenMP backends of ViennaCL
157: Level: advanced
159: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC
161: M*/
163: PETSC_EXTERN PetscErrorCode PCCreate_SAVIENNACL(PC pc)
164: {
165: PC_SAVIENNACL *sac;
167: /*
168: Creates the private data structure for this preconditioner and
169: attach it to the PC object.
170: */
171: PetscNewLog(pc,&sac);
172: pc->data = (void*)sac;
174: /*
175: Initialize the pointer to zero
176: Initialize number of v-cycles to default (1)
177: */
178: sac->SAVIENNACL = 0;
180: /*
181: Set the pointers for the functions that are provided above.
182: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
183: are called, they will automatically call these functions. Note we
184: choose not to provide a couple of these functions since they are
185: not needed.
186: */
187: pc->ops->apply = PCApply_SAVIENNACL;
188: pc->ops->applytranspose = 0;
189: pc->ops->setup = PCSetUp_SAVIENNACL;
190: pc->ops->destroy = PCDestroy_SAVIENNACL;
191: pc->ops->setfromoptions = PCSetFromOptions_SAVIENNACL;
192: pc->ops->view = 0;
193: pc->ops->applyrichardson = 0;
194: pc->ops->applysymmetricleft = 0;
195: pc->ops->applysymmetricright = 0;
196: return 0;
197: }