Actual source code: chowiluviennacl.cxx
petsc-3.12.5 2020-03-29
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: #define PETSC_SKIP_IMMINTRIN_H_CUDAWORKAROUND 1
11: #include <petsc/private/pcimpl.h>
12: #include <../src/mat/impls/aij/seq/aij.h>
13: #include <../src/vec/vec/impls/dvecimpl.h>
14: #include <../src/mat/impls/aij/seq/seqviennacl/viennaclmatimpl.h>
15: #include <../src/vec/vec/impls/seq/seqviennacl/viennaclvecimpl.h>
16: #include <viennacl/linalg/detail/ilu/chow_patel_ilu.hpp>
18: /*
19: Private context (data structure) for the CHOWILUVIENNACL preconditioner.
20: */
21: typedef struct {
22: viennacl::linalg::chow_patel_ilu_precond< viennacl::compressed_matrix<PetscScalar> > *CHOWILUVIENNACL;
23: } PC_CHOWILUVIENNACL;
26: /* -------------------------------------------------------------------------- */
27: /*
28: PCSetUp_CHOWILUVIENNACL - Prepares for the use of the CHOWILUVIENNACL preconditioner
29: by setting data structures and options.
31: Input Parameter:
32: . pc - the preconditioner context
34: Application Interface Routine: PCSetUp()
36: Notes:
37: The interface routine PCSetUp() is not usually called directly by
38: the user, but instead is called by PCApply() if necessary.
39: */
40: static PetscErrorCode PCSetUp_CHOWILUVIENNACL(PC pc)
41: {
42: PC_CHOWILUVIENNACL *ilu = (PC_CHOWILUVIENNACL*)pc->data;
43: PetscBool flg = PETSC_FALSE;
44: PetscErrorCode ierr;
45: Mat_SeqAIJViennaCL *gpustruct;
48: PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQAIJVIENNACL,&flg);
49: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Currently only handles ViennaCL matrices");
50: if (pc->setupcalled != 0) {
51: try {
52: delete ilu->CHOWILUVIENNACL;
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 CHOWILUVIENNACL preconditioner");
61: #else
62: MatViennaCLCopyToGPU(pc->pmat);
63: gpustruct = (Mat_SeqAIJViennaCL*)(pc->pmat->spptr);
64:
65: viennacl::linalg::chow_patel_tag ilu_tag;
66: ViennaCLAIJMatrix *mat = (ViennaCLAIJMatrix*)gpustruct->mat;
67: ilu->CHOWILUVIENNACL = new viennacl::linalg::chow_patel_ilu_precond<viennacl::compressed_matrix<PetscScalar> >(*mat, ilu_tag);
68: #endif
69: } catch(char *ex) {
70: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
71: }
72: return(0);
73: }
75: /* -------------------------------------------------------------------------- */
76: /*
77: PCApply_CHOWILUVIENNACL - Applies the CHOWILUVIENNACL preconditioner to a vector.
79: Input Parameters:
80: . pc - the preconditioner context
81: . x - input vector
83: Output Parameter:
84: . y - output vector
86: Application Interface Routine: PCApply()
87: */
88: static PetscErrorCode PCApply_CHOWILUVIENNACL(PC pc,Vec x,Vec y)
89: {
90: PC_CHOWILUVIENNACL *ilu = (PC_CHOWILUVIENNACL*)pc->data;
91: PetscErrorCode ierr;
92: PetscBool flg1,flg2;
93: viennacl::vector<PetscScalar> const *xarray=NULL;
94: viennacl::vector<PetscScalar> *yarray=NULL;
97: /*how to apply a certain fixed number of iterations?*/
98: PetscObjectTypeCompare((PetscObject)x,VECSEQVIENNACL,&flg1);
99: PetscObjectTypeCompare((PetscObject)y,VECSEQVIENNACL,&flg2);
100: if (!(flg1 && flg2)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP, "Currently only handles ViennaCL vectors");
101: if (!ilu->CHOWILUVIENNACL) {
102: PCSetUp_CHOWILUVIENNACL(pc);
103: }
104: VecSet(y,0.0);
105: VecViennaCLGetArrayRead(x,&xarray);
106: VecViennaCLGetArrayWrite(y,&yarray);
107: try {
108: #if defined(PETSC_USE_COMPLEX)
110: #else
111: *yarray = *xarray;
112: ilu->CHOWILUVIENNACL->apply(*yarray);
113: #endif
114: } catch(char * ex) {
115: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
116: }
117: VecViennaCLRestoreArrayRead(x,&xarray);
118: VecViennaCLRestoreArrayWrite(y,&yarray);
119: PetscObjectStateIncrease((PetscObject)y);
120: return(0);
121: }
122: /* -------------------------------------------------------------------------- */
123: /*
124: PCDestroy_CHOWILUVIENNACL - Destroys the private context for the CHOWILUVIENNACL preconditioner
125: that was created with PCCreate_CHOWILUVIENNACL().
127: Input Parameter:
128: . pc - the preconditioner context
130: Application Interface Routine: PCDestroy()
131: */
132: static PetscErrorCode PCDestroy_CHOWILUVIENNACL(PC pc)
133: {
134: PC_CHOWILUVIENNACL *ilu = (PC_CHOWILUVIENNACL*)pc->data;
138: if (ilu->CHOWILUVIENNACL) {
139: try {
140: delete ilu->CHOWILUVIENNACL;
141: } catch(char *ex) {
142: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex);
143: }
144: }
146: /*
147: Free the private data structure that was hanging off the PC
148: */
149: PetscFree(pc->data);
150: return(0);
151: }
153: static PetscErrorCode PCSetFromOptions_CHOWILUVIENNACL(PetscOptionItems *PetscOptionsObject,PC pc)
154: {
158: PetscOptionsHead(PetscOptionsObject,"CHOWILUVIENNACL options");
159: PetscOptionsTail();
160: return(0);
161: }
163: /* -------------------------------------------------------------------------- */
166: /*MC
167: PCCHOWILUViennaCL - A smoothed agglomeration algorithm that can be used via the CUDA, OpenCL, and OpenMP backends of ViennaCL
169: Level: advanced
171: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC
173: M*/
175: PETSC_EXTERN PetscErrorCode PCCreate_CHOWILUVIENNACL(PC pc)
176: {
177: PC_CHOWILUVIENNACL *ilu;
181: /*
182: Creates the private data structure for this preconditioner and
183: attach it to the PC object.
184: */
185: PetscNewLog(pc,&ilu);
186: pc->data = (void*)ilu;
188: /*
189: Initialize the pointer to zero
190: Initialize number of v-cycles to default (1)
191: */
192: ilu->CHOWILUVIENNACL = 0;
194: /*
195: Set the pointers for the functions that are provided above.
196: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
197: are called, they will automatically call these functions. Note we
198: choose not to provide a couple of these functions since they are
199: not needed.
200: */
201: pc->ops->apply = PCApply_CHOWILUVIENNACL;
202: pc->ops->applytranspose = 0;
203: pc->ops->setup = PCSetUp_CHOWILUVIENNACL;
204: pc->ops->destroy = PCDestroy_CHOWILUVIENNACL;
205: pc->ops->setfromoptions = PCSetFromOptions_CHOWILUVIENNACL;
206: pc->ops->view = 0;
207: pc->ops->applyrichardson = 0;
208: pc->ops->applysymmetricleft = 0;
209: pc->ops->applysymmetricright = 0;
210: return(0);
211: }