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