Actual source code: bicgstabcusp.cu
petsc-3.6.4 2016-04-12
2: /* -------------------------------------------------------------------- */
4: /*
5: Include files needed for the CUSP BiCGSTAB preconditioner:
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/krylov/bicgstab.h>
13: #include <../src/vec/vec/impls/dvecimpl.h>
14: #include <../src/mat/impls/aij/seq/seqcusp/cuspmatimpl.h>
17: /*
18: Private context (data structure) for the CUSP BiCGStab preconditioner.
19: */
20: typedef struct {
21: PetscInt maxits;
22: PetscReal rtol;
23: PetscBool monitorverbose;
24: CUSPMATRIX * mat;
25: } PC_BiCGStabCUSP;
29: static PetscErrorCode PCBiCGStabCUSPSetTolerance_BiCGStabCUSP(PC pc,PetscReal rtol)
30: {
31: PC_BiCGStabCUSP *bicg = (PC_BiCGStabCUSP*)pc->data;
34: bicg->rtol = rtol;
35: return(0);
36: }
40: static PetscErrorCode PCBiCGStabCUSPSetUseVerboseMonitor_BiCGStabCUSP(PC pc, PetscBool useverbose)
41: {
42: PC_BiCGStabCUSP *bicg = (PC_BiCGStabCUSP*)pc->data;
45: bicg->monitorverbose = useverbose;
46: return(0);
47: }
51: PetscErrorCode PCBiCGStabCUSPSetUseVerboseMonitor(PC pc, PetscBool useverbose)
52: {
57: PetscTryMethod(pc, "PCBiCGStabCUSPSetUseVerboseMonitors_C",(PC,PetscBool),(pc,useverbose));
58: return(0);
59: }
63: static PetscErrorCode PCBiCGStabCUSPSetIterations_BiCGStabCUSP(PC pc, PetscInt its)
64: {
65: PC_BiCGStabCUSP *bicg = (PC_BiCGStabCUSP*)pc->data;
68: bicg->maxits = its;
69: return(0);
70: }
74: PetscErrorCode PCBiCGStabCUSPSetITerations(PC pc, PetscInt its)
75: {
80: PetscTryMethod(pc, "PCBiCGStabCUSPSetIterations_C",(PC,PetscInt),(pc,its));
81: return(0);
82: }
86: PetscErrorCode PCBiCGStabCUSPSetTolerance(PC pc, PetscReal rtol)
87: {
92: PetscTryMethod(pc, "PCBiCGStabCUSPSetTolerance_C",(PC,PetscReal),(pc,rtol));
93: return(0);
94: }
96: /* -------------------------------------------------------------------------- */
97: /*
98: PCSetUp_BiCGStabCUSP - Prepares for the use of the CUSP BiCGStab preconditioner
99: by setting data structures and options.
101: Input Parameter:
102: . pc - the preconditioner context
104: Application Interface Routine: PCSetUp()
106: Notes:
107: The interface routine PCSetUp() is not usually called directly by
108: the user, but instead is called by PCApply() if necessary.
109: */
112: static PetscErrorCode PCSetUp_BiCGStabCUSP(PC pc)
113: {
114: PC_BiCGStabCUSP *bicg = (PC_BiCGStabCUSP*)pc->data;
115: PetscBool flg = PETSC_FALSE;
116: Mat_SeqAIJCUSP *gpustruct;
117: PetscErrorCode ierr;
120: PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQAIJCUSP,&flg);
121: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Currently only handles CUSP matrices");
122: try {
123: MatCUSPCopyToGPU(pc->pmat);
124: gpustruct = (Mat_SeqAIJCUSP*)(pc->pmat->spptr);
125: bicg->mat = (CUSPMATRIX*)gpustruct->mat;
126: } catch(char *ex) {
127: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CUSP error: %s",ex);
128: }
129: return(0);
130: }
132: /* -------------------------------------------------------------------------- */
133: /*
134: PCApply_BiCGStabCUSP - Applies the BiCGStabCUSP preconditioner to a vector.
136: Input Parameters:
137: . pc - the preconditioner context
138: . x - input vector
140: Output Parameter:
141: . y - output vector
143: Application Interface Routine: PCApply()
144: */
147: static PetscErrorCode PCApply_BiCGStabCUSP(PC pc,Vec x,Vec y)
148: {
149: PC_BiCGStabCUSP *bicg = (PC_BiCGStabCUSP*)pc->data;
150: PetscErrorCode ierr;
151: PetscBool flg1,flg2;
152: CUSPARRAY *xarray=NULL,*yarray=NULL;
155: PetscObjectTypeCompare((PetscObject)x,VECSEQCUSP,&flg1);
156: PetscObjectTypeCompare((PetscObject)y,VECSEQCUSP,&flg2);
157: if (!(flg1 && flg2)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP, "Currently only handles CUSP vectors");
158: if (!bicg->mat) {
159: PCSetUp_BiCGStabCUSP(pc);
160: }
161: VecSet(y,0.0);
162: VecCUSPGetArrayRead(x,&xarray);
163: VecCUSPGetArrayWrite(y,&yarray);
164: try {
165: #if defined(CUSP_VERSION) && CUSP_VERSION >= 500
166: cusp::monitor<PetscReal> monitor(*xarray,bicg->maxits,bicg->rtol);
167: cusp::krylov::bicgstab(*bicg->mat,*yarray,*xarray,monitor);
168: #else
169: cusp::default_monitor<PetscReal> monitor(*xarray,bicg->maxits,bicg->rtol);
170: if (bicg->monitorverbose) {
171: cusp::verbose_monitor<PetscReal> verbosemonitor(*xarray,bicg->maxits,bicg->rtol);
172: cusp::krylov::bicgstab(*bicg->mat,*yarray,*xarray,verbosemonitor);
173: } else {
174: cusp::krylov::bicgstab(*bicg->mat,*yarray,*xarray,monitor);
175: }
176: #endif
177: } catch(char *ex) {
178: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CUSP error: %s", ex);
179: }
180: VecCUSPRestoreArrayRead(x,&xarray);
181: VecCUSPRestoreArrayWrite(y,&yarray);
182: PetscObjectStateIncrease((PetscObject)y);
183: return(0);
184: }
185: /* -------------------------------------------------------------------------- */
186: /*
187: PCDestroy_BiCGStabCUSP - Destroys the private context for the BiCGStabCUSP preconditioner
188: that was created with PCCreate_BiCGStabCUSP().
190: Input Parameter:
191: . pc - the preconditioner context
193: Application Interface Routine: PCDestroy()
194: */
197: static PetscErrorCode PCDestroy_BiCGStabCUSP(PC pc)
198: {
199: PetscErrorCode ierr;
202: /*
203: Free the private data structure that was hanging off the PC
204: */
205: PetscFree(pc->data);
206: return(0);
207: }
211: static PetscErrorCode PCSetFromOptions_BiCGStabCUSP(PetscOptions *PetscOptionsObject,PC pc)
212: {
213: PC_BiCGStabCUSP *bicg = (PC_BiCGStabCUSP*)pc->data;
214: PetscErrorCode ierr;
217: PetscOptionsHead(PetscOptionsObject,"BiCGStabCUSP options");
218: PetscOptionsReal("-pc_bicgstabcusp_rtol","relative tolerance for BiCGStabCUSP preconditioner","PCBiCGStabCUSPSetTolerance",bicg->rtol,&bicg->rtol,0);
219: PetscOptionsInt("-pc_bicgstabcusp_max_it","maximum iterations for BiCGStabCUSP preconditioner","PCBiCGStabCUSPSetIterations",bicg->maxits,&bicg->maxits,0);
220: PetscOptionsBool("-pc_bicgstabcusp_monitor_verbose","Print information about GPU BiCGStabCUSP iterations","PCBiCGStabCUSPSetUseVerboseMonitor",bicg->monitorverbose,&bicg->monitorverbose,0);
221: PetscOptionsTail();
222: return(0);
223: }
225: /* -------------------------------------------------------------------------- */
229: PETSC_EXTERN PetscErrorCode PCCreate_BiCGStabCUSP(PC pc)
230: {
231: PC_BiCGStabCUSP *bicg;
232: PetscErrorCode ierr;
235: /*
236: Creates the private data structure for this preconditioner and
237: attach it to the PC object.
238: */
239: PetscNewLog(pc,&bicg);
240: /*
241: Set default values. We don't actually want to set max iterations as far as I know, but the Cusp monitor requires them so we use a large number.
242: */
243: bicg->maxits = 1000;
244: bicg->rtol = 1.e-1;
245: bicg->monitorverbose = PETSC_FALSE;
246: pc->data = (void*)bicg;
247: /*
248: Set the pointers for the functions that are provided above.
249: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
250: are called, they will automatically call these functions. Note we
251: choose not to provide a couple of these functions since they are
252: not needed.
253: */
254: pc->ops->apply = PCApply_BiCGStabCUSP;
255: pc->ops->applytranspose = 0;
256: pc->ops->setup = PCSetUp_BiCGStabCUSP;
257: pc->ops->destroy = PCDestroy_BiCGStabCUSP;
258: pc->ops->setfromoptions = PCSetFromOptions_BiCGStabCUSP;
259: pc->ops->view = 0;
260: pc->ops->applyrichardson = 0;
261: pc->ops->applysymmetricleft = 0;
262: pc->ops->applysymmetricright = 0;
264: PetscObjectComposeFunction((PetscObject)pc,"PCBiCGStabCUSPSetTolerance_C",PCBiCGStabCUSPSetTolerance_BiCGStabCUSP);
265: PetscObjectComposeFunction((PetscObject)pc, "PCBiCGStabCUSPSetIterations_C",PCBiCGStabCUSPSetIterations_BiCGStabCUSP);
266: PetscObjectComposeFunction((PetscObject)pc, "PCBiCGStabCUSPSetUseVerboseMonitor_C", PCBiCGStabCUSPSetUseVerboseMonitor_BiCGStabCUSP);
267: return(0);
268: }