xref: /petsc/src/ksp/pc/impls/saviennacl/saviennacl.cxx (revision 2bf68e3e0f2a61f71e7c65bee250bfa1c8ce0cdb)
1 
2 /*  -------------------------------------------------------------------- */
3 
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 #include <petsc/private/pcimpl.h>   /*I "petscpc.h" I*/
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/amg.hpp>
15 
16 /*
17    Private context (data structure) for the SAVIENNACL preconditioner.
18 */
19 typedef struct {
20   viennacl::linalg::amg_precond<viennacl::compressed_matrix<PetscScalar> > *SAVIENNACL;
21 } PC_SAVIENNACL;
22 
23 
24 /* -------------------------------------------------------------------------- */
25 /*
26    PCSetUp_SAVIENNACL - Prepares for the use of the SAVIENNACL preconditioner
27                         by setting data structures and options.
28 
29    Input Parameter:
30 .  pc - the preconditioner context
31 
32    Application Interface Routine: PCSetUp()
33 
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 #undef __FUNCT__
39 #define __FUNCT__ "PCSetUp_SAVIENNACL"
40 static PetscErrorCode PCSetUp_SAVIENNACL(PC pc)
41 {
42   PC_SAVIENNACL      *sa = (PC_SAVIENNACL*)pc->data;
43   PetscBool          flg = PETSC_FALSE;
44   PetscErrorCode     ierr;
45   Mat_SeqAIJViennaCL *gpustruct;
46 
47   PetscFunctionBegin;
48   ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQAIJVIENNACL,&flg);CHKERRQ(ierr);
49   if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Currently only handles ViennaCL matrices");
50   if (pc->setupcalled != 0) {
51     try {
52       delete sa->SAVIENNACL;
53     } catch(char *ex) {
54       SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
55     }
56   }
57   try {
58 #if defined(PETSC_USE_COMPLEX)
59     gpustruct = NULL;
60     SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"No support for complex arithmetic in SAVIENNACL preconditioner");
61 #else
62     ierr      = MatViennaCLCopyToGPU(pc->pmat);CHKERRQ(ierr);
63     gpustruct = (Mat_SeqAIJViennaCL*)(pc->pmat->spptr);
64 
65     viennacl::linalg::amg_tag amg_tag_sa_pmis;
66     amg_tag_sa_pmis.set_coarsening_method(viennacl::linalg::AMG_COARSENING_METHOD_MIS2_AGGREGATION);
67     amg_tag_sa_pmis.set_interpolation_method(viennacl::linalg::AMG_INTERPOLATION_METHOD_SMOOTHED_AGGREGATION);
68     ViennaCLAIJMatrix *mat = (ViennaCLAIJMatrix*)gpustruct->mat;
69     sa->SAVIENNACL = new viennacl::linalg::amg_precond<viennacl::compressed_matrix<PetscScalar> >(*mat, amg_tag_sa_pmis);
70     sa->SAVIENNACL->setup();
71 #endif
72   } catch(char *ex) {
73     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
74   }
75   PetscFunctionReturn(0);
76 }
77 
78 /* -------------------------------------------------------------------------- */
79 /*
80    PCApply_SAVIENNACL - Applies the SAVIENNACL preconditioner to a vector.
81 
82    Input Parameters:
83 .  pc - the preconditioner context
84 .  x - input vector
85 
86    Output Parameter:
87 .  y - output vector
88 
89    Application Interface Routine: PCApply()
90  */
91 #undef __FUNCT__
92 #define __FUNCT__ "PCApply_SAVIENNACL"
93 static PetscErrorCode PCApply_SAVIENNACL(PC pc,Vec x,Vec y)
94 {
95   PC_SAVIENNACL                 *sac = (PC_SAVIENNACL*)pc->data;
96   PetscErrorCode                ierr;
97   PetscBool                     flg1,flg2;
98   viennacl::vector<PetscScalar> const *xarray=NULL;
99   viennacl::vector<PetscScalar> *yarray=NULL;
100 
101   PetscFunctionBegin;
102   /*how to apply a certain fixed number of iterations?*/
103   ierr = PetscObjectTypeCompare((PetscObject)x,VECSEQVIENNACL,&flg1);CHKERRQ(ierr);
104   ierr = PetscObjectTypeCompare((PetscObject)y,VECSEQVIENNACL,&flg2);CHKERRQ(ierr);
105   if (!(flg1 && flg2)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP, "Currently only handles ViennaCL vectors");
106   if (!sac->SAVIENNACL) {
107     ierr = PCSetUp_SAVIENNACL(pc);CHKERRQ(ierr);
108   }
109   ierr = VecSet(y,0.0);CHKERRQ(ierr);
110   ierr = VecViennaCLGetArrayRead(x,&xarray);CHKERRQ(ierr);
111   ierr = VecViennaCLGetArrayWrite(y,&yarray);CHKERRQ(ierr);
112   try {
113 #if defined(PETSC_USE_COMPLEX)
114 
115 #else
116     *yarray = *xarray;
117     sac->SAVIENNACL->apply(*yarray);
118 #endif
119   } catch(char * ex) {
120     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
121   }
122   ierr = VecViennaCLRestoreArrayRead(x,&xarray);CHKERRQ(ierr);
123   ierr = VecViennaCLRestoreArrayWrite(y,&yarray);CHKERRQ(ierr);
124   ierr = PetscObjectStateIncrease((PetscObject)y);CHKERRQ(ierr);
125   PetscFunctionReturn(0);
126 }
127 /* -------------------------------------------------------------------------- */
128 /*
129    PCDestroy_SAVIENNACL - Destroys the private context for the SAVIENNACL preconditioner
130    that was created with PCCreate_SAVIENNACL().
131 
132    Input Parameter:
133 .  pc - the preconditioner context
134 
135    Application Interface Routine: PCDestroy()
136 */
137 #undef __FUNCT__
138 #define __FUNCT__ "PCDestroy_SAVIENNACL"
139 static PetscErrorCode PCDestroy_SAVIENNACL(PC pc)
140 {
141   PC_SAVIENNACL  *sac = (PC_SAVIENNACL*)pc->data;
142   PetscErrorCode ierr;
143 
144   PetscFunctionBegin;
145   if (sac->SAVIENNACL) {
146     try {
147       delete sac->SAVIENNACL;
148     } catch(char * ex) {
149       SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
150     }
151   }
152 
153   /*
154       Free the private data structure that was hanging off the PC
155   */
156   ierr = PetscFree(pc->data);CHKERRQ(ierr);
157   PetscFunctionReturn(0);
158 }
159 
160 #undef __FUNCT__
161 #define __FUNCT__ "PCSetFromOptions_SAVIENNACL"
162 static PetscErrorCode PCSetFromOptions_SAVIENNACL(PetscOptionItems *PetscOptionsObject,PC pc)
163 {
164   PetscErrorCode ierr;
165 
166   PetscFunctionBegin;
167   ierr = PetscOptionsHead(PetscOptionsObject,"SAVIENNACL options");CHKERRQ(ierr);
168   ierr = PetscOptionsTail();CHKERRQ(ierr);
169   PetscFunctionReturn(0);
170 }
171 
172 /* -------------------------------------------------------------------------- */
173 
174 
175 /*MC
176      PCSAViennaCL  - A smoothed agglomeration algorithm that can be used via the CUDA, OpenCL, and OpenMP backends of ViennaCL
177 
178    Level: advanced
179 
180 .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC
181 
182 M*/
183 
184 #undef __FUNCT__
185 #define __FUNCT__ "PCCreate_SAVIENNACL"
186 PETSC_EXTERN PetscErrorCode PCCreate_SAVIENNACL(PC pc)
187 {
188   PC_SAVIENNACL  *sac;
189   PetscErrorCode ierr;
190 
191   PetscFunctionBegin;
192   /*
193      Creates the private data structure for this preconditioner and
194      attach it to the PC object.
195   */
196   ierr     = PetscNewLog(pc,&sac);CHKERRQ(ierr);
197   pc->data = (void*)sac;
198 
199   /*
200      Initialize the pointer to zero
201      Initialize number of v-cycles to default (1)
202   */
203   sac->SAVIENNACL = 0;
204 
205   /*
206       Set the pointers for the functions that are provided above.
207       Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
208       are called, they will automatically call these functions.  Note we
209       choose not to provide a couple of these functions since they are
210       not needed.
211   */
212   pc->ops->apply               = PCApply_SAVIENNACL;
213   pc->ops->applytranspose      = 0;
214   pc->ops->setup               = PCSetUp_SAVIENNACL;
215   pc->ops->destroy             = PCDestroy_SAVIENNACL;
216   pc->ops->setfromoptions      = PCSetFromOptions_SAVIENNACL;
217   pc->ops->view                = 0;
218   pc->ops->applyrichardson     = 0;
219   pc->ops->applysymmetricleft  = 0;
220   pc->ops->applysymmetricright = 0;
221   PetscFunctionReturn(0);
222 }
223 
224