Actual source code: chowiluviennacl.cxx
petsc-3.11.4 2019-09-28
2: /* -------------------------------------------------------------------- */
4: /*
5: Include files needed for the ViennaCL Chow-Patel parallel ILU 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/detail/ilu/chow_patel_ilu.hpp>
16: /*
17: Private context (data structure) for the CHOWILUVIENNACL preconditioner.
18: */
19: typedef struct {
20: viennacl::linalg::chow_patel_ilu_precond< viennacl::compressed_matrix<PetscScalar> > *CHOWILUVIENNACL;
21: } PC_CHOWILUVIENNACL;
24: /* -------------------------------------------------------------------------- */
25: /*
26: PCSetUp_CHOWILUVIENNACL - Prepares for the use of the CHOWILUVIENNACL 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_CHOWILUVIENNACL(PC pc)
39: {
40: PC_CHOWILUVIENNACL *ilu = (PC_CHOWILUVIENNACL*)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 ilu->CHOWILUVIENNACL;
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 CHOWILUVIENNACL preconditioner");
59: #else
60: MatViennaCLCopyToGPU(pc->pmat);
61: gpustruct = (Mat_SeqAIJViennaCL*)(pc->pmat->spptr);
62:
63: viennacl::linalg::chow_patel_tag ilu_tag;
64: ViennaCLAIJMatrix *mat = (ViennaCLAIJMatrix*)gpustruct->mat;
65: ilu->CHOWILUVIENNACL = new viennacl::linalg::chow_patel_ilu_precond<viennacl::compressed_matrix<PetscScalar> >(*mat, ilu_tag);
66: #endif
67: } catch(char *ex) {
68: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
69: }
70: return(0);
71: }
73: /* -------------------------------------------------------------------------- */
74: /*
75: PCApply_CHOWILUVIENNACL - Applies the CHOWILUVIENNACL preconditioner to a vector.
77: Input Parameters:
78: . pc - the preconditioner context
79: . x - input vector
81: Output Parameter:
82: . y - output vector
84: Application Interface Routine: PCApply()
85: */
86: static PetscErrorCode PCApply_CHOWILUVIENNACL(PC pc,Vec x,Vec y)
87: {
88: PC_CHOWILUVIENNACL *ilu = (PC_CHOWILUVIENNACL*)pc->data;
89: PetscErrorCode ierr;
90: PetscBool flg1,flg2;
91: viennacl::vector<PetscScalar> const *xarray=NULL;
92: viennacl::vector<PetscScalar> *yarray=NULL;
95: /*how to apply a certain fixed number of iterations?*/
96: PetscObjectTypeCompare((PetscObject)x,VECSEQVIENNACL,&flg1);
97: PetscObjectTypeCompare((PetscObject)y,VECSEQVIENNACL,&flg2);
98: if (!(flg1 && flg2)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP, "Currently only handles ViennaCL vectors");
99: if (!ilu->CHOWILUVIENNACL) {
100: PCSetUp_CHOWILUVIENNACL(pc);
101: }
102: VecSet(y,0.0);
103: VecViennaCLGetArrayRead(x,&xarray);
104: VecViennaCLGetArrayWrite(y,&yarray);
105: try {
106: #if defined(PETSC_USE_COMPLEX)
108: #else
109: *yarray = *xarray;
110: ilu->CHOWILUVIENNACL->apply(*yarray);
111: #endif
112: } catch(char * ex) {
113: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
114: }
115: VecViennaCLRestoreArrayRead(x,&xarray);
116: VecViennaCLRestoreArrayWrite(y,&yarray);
117: PetscObjectStateIncrease((PetscObject)y);
118: return(0);
119: }
120: /* -------------------------------------------------------------------------- */
121: /*
122: PCDestroy_CHOWILUVIENNACL - Destroys the private context for the CHOWILUVIENNACL preconditioner
123: that was created with PCCreate_CHOWILUVIENNACL().
125: Input Parameter:
126: . pc - the preconditioner context
128: Application Interface Routine: PCDestroy()
129: */
130: static PetscErrorCode PCDestroy_CHOWILUVIENNACL(PC pc)
131: {
132: PC_CHOWILUVIENNACL *ilu = (PC_CHOWILUVIENNACL*)pc->data;
136: if (ilu->CHOWILUVIENNACL) {
137: try {
138: delete ilu->CHOWILUVIENNACL;
139: } catch(char *ex) {
140: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
141: }
142: }
144: /*
145: Free the private data structure that was hanging off the PC
146: */
147: PetscFree(pc->data);
148: return(0);
149: }
151: static PetscErrorCode PCSetFromOptions_CHOWILUVIENNACL(PetscOptionItems *PetscOptionsObject,PC pc)
152: {
156: PetscOptionsHead(PetscOptionsObject,"CHOWILUVIENNACL options");
157: PetscOptionsTail();
158: return(0);
159: }
161: /* -------------------------------------------------------------------------- */
164: /*MC
165: PCCHOWILUViennaCL - A smoothed agglomeration algorithm that can be used via the CUDA, OpenCL, and OpenMP backends of ViennaCL
167: Level: advanced
169: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC
171: M*/
173: PETSC_EXTERN PetscErrorCode PCCreate_CHOWILUVIENNACL(PC pc)
174: {
175: PC_CHOWILUVIENNACL *ilu;
179: /*
180: Creates the private data structure for this preconditioner and
181: attach it to the PC object.
182: */
183: PetscNewLog(pc,&ilu);
184: pc->data = (void*)ilu;
186: /*
187: Initialize the pointer to zero
188: Initialize number of v-cycles to default (1)
189: */
190: ilu->CHOWILUVIENNACL = 0;
192: /*
193: Set the pointers for the functions that are provided above.
194: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
195: are called, they will automatically call these functions. Note we
196: choose not to provide a couple of these functions since they are
197: not needed.
198: */
199: pc->ops->apply = PCApply_CHOWILUVIENNACL;
200: pc->ops->applytranspose = 0;
201: pc->ops->setup = PCSetUp_CHOWILUVIENNACL;
202: pc->ops->destroy = PCDestroy_CHOWILUVIENNACL;
203: pc->ops->setfromoptions = PCSetFromOptions_CHOWILUVIENNACL;
204: pc->ops->view = 0;
205: pc->ops->applyrichardson = 0;
206: pc->ops->applysymmetricleft = 0;
207: pc->ops->applysymmetricright = 0;
208: return(0);
209: }