Actual source code: sacusppoly.cu
petsc-3.7.3 2016-08-01
2: /* -------------------------------------------------------------------- */
4: /*
5: Include files needed for the CUSP Smoothed Aggregation preconditioner with Chebyshev polynomial smoothing:
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 <cusp/monitor.h>
12: #include <cusp/version.h>
13: #define USE_POLY_SMOOTHER 1
14: #if CUSP_VERSION >= 400
15: #include <cusp/precond/aggregation/smoothed_aggregation.h>
16: #define cuspsaprecond cusp::precond::aggregation::smoothed_aggregation<PetscInt,PetscScalar,cusp::device_memory>
17: #else
18: #include <cusp/precond/smoothed_aggregation.h>
19: #define cuspsaprecond cusp::precond::smoothed_aggregation<PetscInt,PetscScalar,cusp::device_memory>
20: #endif
21: #undef USE_POLY_SMOOTHER
22: #include <../src/vec/vec/impls/dvecimpl.h>
23: #include <../src/mat/impls/aij/seq/seqcusp/cuspmatimpl.h>
24: #include <../src/vec/vec/impls/seq/seqcusp/cuspvecimpl.h>
26: /*
27: Private context (data structure) for the SACUSPPoly preconditioner.
28: */
29: typedef struct {
30: cuspsaprecond * SACUSPPoly;
31: /*int cycles; */
32: } PC_SACUSPPoly;
35: /* -------------------------------------------------------------------------- */
36: /*
37: PCSetUp_SACUSPPoly - Prepares for the use of the SACUSPPoly preconditioner
38: by setting data structures and options.
40: Input Parameter:
41: . pc - the preconditioner context
43: Application Interface Routine: PCSetUp()
45: Notes:
46: The interface routine PCSetUp() is not usually called directly by
47: the user, but instead is called by PCApply() if necessary.
48: */
51: static PetscErrorCode PCSetUp_SACUSPPoly(PC pc)
52: {
53: PC_SACUSPPoly *sa = (PC_SACUSPPoly*)pc->data;
54: PetscBool flg = PETSC_FALSE;
56: #if !defined(PETSC_USE_COMPLEX)
57: // protect these in order to avoid compiler warnings. This preconditioner does
58: // not work for complex types.
59: Mat_SeqAIJCUSP *gpustruct;
60: #endif
62: PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQAIJCUSP,&flg);
63: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Currently only handles CUSP matrices");
64: if (pc->setupcalled != 0) {
65: try {
66: delete sa->SACUSPPoly;
67: } catch(char *ex) {
68: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CUSP error: %s", ex);
69: }
70: }
71: try {
72: #if defined(PETSC_USE_COMPLEX)
73: sa->SACUSPPoly = 0;CHKERRQ(1); /* TODO */
74: #else
75: MatCUSPCopyToGPU(pc->pmat);
76: gpustruct = (Mat_SeqAIJCUSP*)(pc->pmat->spptr);
77:
78: if (gpustruct->format==MAT_CUSP_ELL) {
79: CUSPMATRIXELL *mat = (CUSPMATRIXELL*)gpustruct->mat;
80: sa->SACUSPPoly = new cuspsaprecond(*mat);
81: } else if (gpustruct->format==MAT_CUSP_DIA) {
82: CUSPMATRIXDIA *mat = (CUSPMATRIXDIA*)gpustruct->mat;
83: sa->SACUSPPoly = new cuspsaprecond(*mat);
84: } else {
85: CUSPMATRIX *mat = (CUSPMATRIX*)gpustruct->mat;
86: sa->SACUSPPoly = new cuspsaprecond(*mat);
87: }
88: #endif
89: } catch(char *ex) {
90: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CUSP error: %s", ex);
91: }
92: /*PetscOptionsInt("-pc_sacusp_cycles","Number of v-cycles to perform","PCSACUSPSetCycles",sa->cycles,
93: &sa->cycles,NULL);*/
94: return(0);
95: }
99: static PetscErrorCode PCApplyRichardson_SACUSPPoly(PC pc, Vec b, Vec y, Vec w,PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt its, PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
100: {
101: #if !defined(PETSC_USE_COMPLEX)
102: // protect these in order to avoid compiler warnings. This preconditioner does
103: // not work for complex types.
104: PC_SACUSPPoly *sac = (PC_SACUSPPoly*)pc->data;
105: #endif
107: CUSPARRAY *barray,*yarray;
110: /* how to incorporate dtol, guesszero, w?*/
112: VecCUSPGetArrayRead(b,&barray);
113: VecCUSPGetArrayReadWrite(y,&yarray);
114: #if defined(CUSP_VERSION) && CUSP_VERSION >= 500
115: cusp::monitor<PetscReal> monitor(*barray,its,rtol,abstol);
116: #else
117: cusp::default_monitor<PetscReal> monitor(*barray,its,rtol,abstol);
118: #endif
119: #if defined(PETSC_USE_COMPLEX)
120: CHKERRQ(1); /* TODO */
121: #else
122: sac->SACUSPPoly->solve(*barray,*yarray,monitor);
123: #endif
124: *outits = monitor.iteration_count();
125: if (monitor.converged()) *reason = PCRICHARDSON_CONVERGED_RTOL; /* how to discern between converging from RTOL or ATOL?*/
126: else *reason = PCRICHARDSON_CONVERGED_ITS;
128: PetscObjectStateIncrease((PetscObject)y);
129: VecCUSPRestoreArrayRead(b,&barray);
130: VecCUSPRestoreArrayReadWrite(y,&yarray);
131: return(0);
132: }
134: /* -------------------------------------------------------------------------- */
135: /*
136: PCApply_SACUSPPoly - Applies the SACUSPPoly preconditioner to a vector.
138: Input Parameters:
139: . pc - the preconditioner context
140: . x - input vector
142: Output Parameter:
143: . y - output vector
145: Application Interface Routine: PCApply()
146: */
149: static PetscErrorCode PCApply_SACUSPPoly(PC pc,Vec x,Vec y)
150: {
151: PC_SACUSPPoly *sac = (PC_SACUSPPoly*)pc->data;
153: PetscBool flg1,flg2;
154: CUSPARRAY *xarray=NULL,*yarray=NULL;
157: /*how to apply a certain fixed number of iterations?*/
158: PetscObjectTypeCompare((PetscObject)x,VECSEQCUSP,&flg1);
159: PetscObjectTypeCompare((PetscObject)y,VECSEQCUSP,&flg2);
160: if (!(flg1 && flg2)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP, "Currently only handles CUSP vectors");
161: if (!sac->SACUSPPoly) {
162: PCSetUp_SACUSPPoly(pc);
163: }
164: VecSet(y,0.0);
165: VecCUSPGetArrayRead(x,&xarray);
166: VecCUSPGetArrayWrite(y,&yarray);
167: try {
168: #if defined(PETSC_USE_COMPLEX)
169: CHKERRQ(1); /* TODO */
170: #else
171: cusp::multiply(*sac->SACUSPPoly,*xarray,*yarray);
172: #endif
173: } catch(char *ex) {
174: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CUSP error: %s", ex);
175: }
176: VecCUSPRestoreArrayRead(x,&xarray);
177: VecCUSPRestoreArrayWrite(y,&yarray);
178: PetscObjectStateIncrease((PetscObject)y);
179: return(0);
180: }
181: /* -------------------------------------------------------------------------- */
182: /*
183: PCDestroy_SACUSPPoly - Destroys the private context for the SACUSPPoly preconditioner
184: that was created with PCCreate_SACUSPPoly().
186: Input Parameter:
187: . pc - the preconditioner context
189: Application Interface Routine: PCDestroy()
190: */
193: static PetscErrorCode PCDestroy_SACUSPPoly(PC pc)
194: {
195: PC_SACUSPPoly *sac = (PC_SACUSPPoly*)pc->data;
199: if (sac->SACUSPPoly) {
200: try {
201: delete sac->SACUSPPoly;
202: } catch(char *ex) {
203: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CUSP error: %s", ex);
204: }
205: }
207: /*
208: Free the private data structure that was hanging off the PC
209: */
210: PetscFree(sac);
211: return(0);
212: }
216: static PetscErrorCode PCSetFromOptions_SACUSPPoly(PetscOptionItems *PetscOptionsObject,PC pc)
217: {
221: PetscOptionsHead(PetscOptionsObject,"SACUSPPoly options");
222: PetscOptionsTail();
223: return(0);
224: }
226: /* -------------------------------------------------------------------------- */
230: PETSC_EXTERN PetscErrorCode PCCreate_SACUSPPoly(PC pc)
231: {
232: PC_SACUSPPoly *sac;
236: /*
237: Creates the private data structure for this preconditioner and
238: attach it to the PC object.
239: */
240: PetscNewLog(pc,&sac);
241: pc->data = (void*)sac;
243: /*
244: Initialize the pointer to zero
245: Initialize number of v-cycles to default (1)
246: */
247: sac->SACUSPPoly = 0;
248: /*sac->cycles=1;*/
251: /*
252: Set the pointers for the functions that are provided above.
253: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
254: are called, they will automatically call these functions. Note we
255: choose not to provide a couple of these functions since they are
256: not needed.
257: */
258: pc->ops->apply = PCApply_SACUSPPoly;
259: pc->ops->applytranspose = 0;
260: pc->ops->setup = PCSetUp_SACUSPPoly;
261: pc->ops->destroy = PCDestroy_SACUSPPoly;
262: pc->ops->setfromoptions = PCSetFromOptions_SACUSPPoly;
263: pc->ops->view = 0;
264: pc->ops->applyrichardson = PCApplyRichardson_SACUSPPoly;
265: pc->ops->applysymmetricleft = 0;
266: pc->ops->applysymmetricright = 0;
267: return(0);
268: }