Actual source code: sacusppoly.cu
petsc-3.6.1 2015-08-06
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: */
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>
25: /*
26: Private context (data structure) for the SACUSPPoly preconditioner.
27: */
28: typedef struct {
29: cuspsaprecond * SACUSPPoly;
30: /*int cycles; */
31: } PC_SACUSPPoly;
34: /* -------------------------------------------------------------------------- */
35: /*
36: PCSetUp_SACUSPPoly - Prepares for the use of the SACUSPPoly preconditioner
37: by setting data structures and options.
39: Input Parameter:
40: . pc - the preconditioner context
42: Application Interface Routine: PCSetUp()
44: Notes:
45: The interface routine PCSetUp() is not usually called directly by
46: the user, but instead is called by PCApply() if necessary.
47: */
50: static PetscErrorCode PCSetUp_SACUSPPoly(PC pc)
51: {
52: PC_SACUSPPoly *sa = (PC_SACUSPPoly*)pc->data;
53: PetscBool flg = PETSC_FALSE;
55: #if !defined(PETSC_USE_COMPLEX)
56: // protect these in order to avoid compiler warnings. This preconditioner does
57: // not work for complex types.
58: Mat_SeqAIJCUSP *gpustruct;
59: #endif
61: PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQAIJCUSP,&flg);
62: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Currently only handles CUSP matrices");
63: if (pc->setupcalled != 0) {
64: try {
65: delete sa->SACUSPPoly;
66: } catch(char *ex) {
67: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CUSP error: %s", ex);
68: }
69: }
70: try {
71: #if defined(PETSC_USE_COMPLEX)
72: sa->SACUSPPoly = 0;CHKERRQ(1); /* TODO */
73: #else
74: MatCUSPCopyToGPU(pc->pmat);
75: gpustruct = (Mat_SeqAIJCUSP*)(pc->pmat->spptr);
76:
77: if (gpustruct->format==MAT_CUSP_ELL) {
78: CUSPMATRIXELL *mat = (CUSPMATRIXELL*)gpustruct->mat;
79: sa->SACUSPPoly = new cuspsaprecond(*mat);
80: } else if (gpustruct->format==MAT_CUSP_DIA) {
81: CUSPMATRIXDIA *mat = (CUSPMATRIXDIA*)gpustruct->mat;
82: sa->SACUSPPoly = new cuspsaprecond(*mat);
83: } else {
84: CUSPMATRIX *mat = (CUSPMATRIX*)gpustruct->mat;
85: sa->SACUSPPoly = new cuspsaprecond(*mat);
86: }
87: #endif
88: } catch(char *ex) {
89: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CUSP error: %s", ex);
90: }
91: /*PetscOptionsInt("-pc_sacusp_cycles","Number of v-cycles to perform","PCSACUSPSetCycles",sa->cycles,
92: &sa->cycles,NULL);*/
93: return(0);
94: }
98: 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)
99: {
100: #if !defined(PETSC_USE_COMPLEX)
101: // protect these in order to avoid compiler warnings. This preconditioner does
102: // not work for complex types.
103: PC_SACUSPPoly *sac = (PC_SACUSPPoly*)pc->data;
104: #endif
106: CUSPARRAY *barray,*yarray;
109: /* how to incorporate dtol, guesszero, w?*/
111: VecCUSPGetArrayRead(b,&barray);
112: VecCUSPGetArrayReadWrite(y,&yarray);
113: #if defined(CUSP_VERSION) && CUSP_VERSION >= 500
114: cusp::monitor<PetscReal> monitor(*barray,its,rtol,abstol);
115: #else
116: cusp::default_monitor<PetscReal> monitor(*barray,its,rtol,abstol);
117: #endif
118: #if defined(PETSC_USE_COMPLEX)
119: CHKERRQ(1); /* TODO */
120: #else
121: sac->SACUSPPoly->solve(*barray,*yarray,monitor);
122: #endif
123: *outits = monitor.iteration_count();
124: if (monitor.converged()) *reason = PCRICHARDSON_CONVERGED_RTOL; /* how to discern between converging from RTOL or ATOL?*/
125: else *reason = PCRICHARDSON_CONVERGED_ITS;
127: PetscObjectStateIncrease((PetscObject)y);
128: VecCUSPRestoreArrayRead(b,&barray);
129: VecCUSPRestoreArrayReadWrite(y,&yarray);
130: return(0);
131: }
133: /* -------------------------------------------------------------------------- */
134: /*
135: PCApply_SACUSPPoly - Applies the SACUSPPoly preconditioner to a vector.
137: Input Parameters:
138: . pc - the preconditioner context
139: . x - input vector
141: Output Parameter:
142: . y - output vector
144: Application Interface Routine: PCApply()
145: */
148: static PetscErrorCode PCApply_SACUSPPoly(PC pc,Vec x,Vec y)
149: {
150: PC_SACUSPPoly *sac = (PC_SACUSPPoly*)pc->data;
152: PetscBool flg1,flg2;
153: CUSPARRAY *xarray=NULL,*yarray=NULL;
156: /*how to apply a certain fixed number of iterations?*/
157: PetscObjectTypeCompare((PetscObject)x,VECSEQCUSP,&flg1);
158: PetscObjectTypeCompare((PetscObject)y,VECSEQCUSP,&flg2);
159: if (!(flg1 && flg2)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP, "Currently only handles CUSP vectors");
160: if (!sac->SACUSPPoly) {
161: PCSetUp_SACUSPPoly(pc);
162: }
163: VecSet(y,0.0);
164: VecCUSPGetArrayRead(x,&xarray);
165: VecCUSPGetArrayWrite(y,&yarray);
166: try {
167: #if defined(PETSC_USE_COMPLEX)
168: CHKERRQ(1); /* TODO */
169: #else
170: cusp::multiply(*sac->SACUSPPoly,*xarray,*yarray);
171: #endif
172: } catch(char *ex) {
173: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CUSP error: %s", ex);
174: }
175: VecCUSPRestoreArrayRead(x,&xarray);
176: VecCUSPRestoreArrayWrite(y,&yarray);
177: PetscObjectStateIncrease((PetscObject)y);
178: return(0);
179: }
180: /* -------------------------------------------------------------------------- */
181: /*
182: PCDestroy_SACUSPPoly - Destroys the private context for the SACUSPPoly preconditioner
183: that was created with PCCreate_SACUSPPoly().
185: Input Parameter:
186: . pc - the preconditioner context
188: Application Interface Routine: PCDestroy()
189: */
192: static PetscErrorCode PCDestroy_SACUSPPoly(PC pc)
193: {
194: PC_SACUSPPoly *sac = (PC_SACUSPPoly*)pc->data;
198: if (sac->SACUSPPoly) {
199: try {
200: delete sac->SACUSPPoly;
201: } catch(char *ex) {
202: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CUSP error: %s", ex);
203: }
204: }
206: /*
207: Free the private data structure that was hanging off the PC
208: */
209: PetscFree(sac);
210: return(0);
211: }
215: static PetscErrorCode PCSetFromOptions_SACUSPPoly(PetscOptions *PetscOptionsObject,PC pc)
216: {
220: PetscOptionsHead(PetscOptionsObject,"SACUSPPoly options");
221: PetscOptionsTail();
222: return(0);
223: }
225: /* -------------------------------------------------------------------------- */
229: PETSC_EXTERN PetscErrorCode PCCreate_SACUSPPoly(PC pc)
230: {
231: PC_SACUSPPoly *sac;
235: /*
236: Creates the private data structure for this preconditioner and
237: attach it to the PC object.
238: */
239: PetscNewLog(pc,&sac);
240: pc->data = (void*)sac;
242: /*
243: Initialize the pointer to zero
244: Initialize number of v-cycles to default (1)
245: */
246: sac->SACUSPPoly = 0;
247: /*sac->cycles=1;*/
250: /*
251: Set the pointers for the functions that are provided above.
252: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
253: are called, they will automatically call these functions. Note we
254: choose not to provide a couple of these functions since they are
255: not needed.
256: */
257: pc->ops->apply = PCApply_SACUSPPoly;
258: pc->ops->applytranspose = 0;
259: pc->ops->setup = PCSetUp_SACUSPPoly;
260: pc->ops->destroy = PCDestroy_SACUSPPoly;
261: pc->ops->setfromoptions = PCSetFromOptions_SACUSPPoly;
262: pc->ops->view = 0;
263: pc->ops->applyrichardson = PCApplyRichardson_SACUSPPoly;
264: pc->ops->applysymmetricleft = 0;
265: pc->ops->applysymmetricright = 0;
266: return(0);
267: }