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