Actual source code: bicgstabcusp.cu
petsc-3.7.3 2016-08-01
2: /* -------------------------------------------------------------------- */
4: /*
5: Include files needed for the CUSP BiCGSTAB preconditioner:
6: pcimpl.h - private include file intended for use by all preconditioners
7: */
8: #define PETSC_SKIP_SPINLOCK
10: #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/
11: #include <../src/mat/impls/aij/seq/aij.h>
12: #include <cusp/monitor.h>
13: #include <cusp/krylov/bicgstab.h>
14: #include <../src/vec/vec/impls/dvecimpl.h>
15: #include <../src/mat/impls/aij/seq/seqcusp/cuspmatimpl.h>
16: #include <../src/vec/vec/impls/seq/seqcusp/cuspvecimpl.h>
18: /*
19: Private context (data structure) for the CUSP BiCGStab preconditioner.
20: */
21: typedef struct {
22: PetscInt maxits;
23: PetscReal rtol;
24: PetscBool monitorverbose;
25: CUSPMATRIX * mat;
26: } PC_BiCGStabCUSP;
30: static PetscErrorCode PCBiCGStabCUSPSetTolerance_BiCGStabCUSP(PC pc,PetscReal rtol)
31: {
32: PC_BiCGStabCUSP *bicg = (PC_BiCGStabCUSP*)pc->data;
35: bicg->rtol = rtol;
36: return(0);
37: }
41: static PetscErrorCode PCBiCGStabCUSPSetUseVerboseMonitor_BiCGStabCUSP(PC pc, PetscBool useverbose)
42: {
43: PC_BiCGStabCUSP *bicg = (PC_BiCGStabCUSP*)pc->data;
46: bicg->monitorverbose = useverbose;
47: return(0);
48: }
52: PetscErrorCode PCBiCGStabCUSPSetUseVerboseMonitor(PC pc, PetscBool useverbose)
53: {
58: PetscTryMethod(pc, "PCBiCGStabCUSPSetUseVerboseMonitors_C",(PC,PetscBool),(pc,useverbose));
59: return(0);
60: }
64: static PetscErrorCode PCBiCGStabCUSPSetIterations_BiCGStabCUSP(PC pc, PetscInt its)
65: {
66: PC_BiCGStabCUSP *bicg = (PC_BiCGStabCUSP*)pc->data;
69: bicg->maxits = its;
70: return(0);
71: }
75: PetscErrorCode PCBiCGStabCUSPSetITerations(PC pc, PetscInt its)
76: {
81: PetscTryMethod(pc, "PCBiCGStabCUSPSetIterations_C",(PC,PetscInt),(pc,its));
82: return(0);
83: }
87: PetscErrorCode PCBiCGStabCUSPSetTolerance(PC pc, PetscReal rtol)
88: {
93: PetscTryMethod(pc, "PCBiCGStabCUSPSetTolerance_C",(PC,PetscReal),(pc,rtol));
94: return(0);
95: }
97: /* -------------------------------------------------------------------------- */
98: /*
99: PCSetUp_BiCGStabCUSP - Prepares for the use of the CUSP BiCGStab preconditioner
100: by setting data structures and options.
102: Input Parameter:
103: . pc - the preconditioner context
105: Application Interface Routine: PCSetUp()
107: Notes:
108: The interface routine PCSetUp() is not usually called directly by
109: the user, but instead is called by PCApply() if necessary.
110: */
113: static PetscErrorCode PCSetUp_BiCGStabCUSP(PC pc)
114: {
115: PC_BiCGStabCUSP *bicg = (PC_BiCGStabCUSP*)pc->data;
116: PetscBool flg = PETSC_FALSE;
117: Mat_SeqAIJCUSP *gpustruct;
118: PetscErrorCode ierr;
121: PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQAIJCUSP,&flg);
122: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Currently only handles SEQAIJCUSP matrices");
123: try {
124: MatCUSPCopyToGPU(pc->pmat);
125: gpustruct = (Mat_SeqAIJCUSP*)(pc->pmat->spptr);
126: bicg->mat = (CUSPMATRIX*)gpustruct->mat;
127: } catch(char *ex) {
128: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CUSP error: %s",ex);
129: }
130: return(0);
131: }
133: /* -------------------------------------------------------------------------- */
134: /*
135: PCApply_BiCGStabCUSP - Applies the BiCGStabCUSP 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_BiCGStabCUSP(PC pc,Vec x,Vec y)
149: {
150: PC_BiCGStabCUSP *bicg = (PC_BiCGStabCUSP*)pc->data;
151: PetscErrorCode ierr;
152: PetscBool flg1,flg2;
153: CUSPARRAY *xarray=NULL,*yarray=NULL;
156: PetscObjectTypeCompare((PetscObject)x,VECSEQCUSP,&flg1);
157: PetscObjectTypeCompare((PetscObject)y,VECSEQCUSP,&flg2);
158: if (!(flg1 && flg2)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP, "Currently only handles CUSP vectors");
159: if (!bicg->mat) {
160: PCSetUp_BiCGStabCUSP(pc);
161: }
162: VecSet(y,0.0);
163: VecCUSPGetArrayRead(x,&xarray);
164: VecCUSPGetArrayWrite(y,&yarray);
165: try {
166: #if defined(CUSP_VERSION) && CUSP_VERSION >= 500
167: cusp::monitor<PetscReal> monitor(*xarray,bicg->maxits,bicg->rtol);
168: cusp::krylov::bicgstab(*bicg->mat,*yarray,*xarray,monitor);
169: #else
170: cusp::default_monitor<PetscReal> monitor(*xarray,bicg->maxits,bicg->rtol);
171: if (bicg->monitorverbose) {
172: cusp::verbose_monitor<PetscReal> verbosemonitor(*xarray,bicg->maxits,bicg->rtol);
173: cusp::krylov::bicgstab(*bicg->mat,*yarray,*xarray,verbosemonitor);
174: } else {
175: cusp::krylov::bicgstab(*bicg->mat,*yarray,*xarray,monitor);
176: }
177: #endif
178: } catch(char *ex) {
179: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CUSP error: %s", ex);
180: }
181: VecCUSPRestoreArrayRead(x,&xarray);
182: VecCUSPRestoreArrayWrite(y,&yarray);
183: PetscObjectStateIncrease((PetscObject)y);
184: return(0);
185: }
186: /* -------------------------------------------------------------------------- */
187: /*
188: PCDestroy_BiCGStabCUSP - Destroys the private context for the BiCGStabCUSP preconditioner
189: that was created with PCCreate_BiCGStabCUSP().
191: Input Parameter:
192: . pc - the preconditioner context
194: Application Interface Routine: PCDestroy()
195: */
198: static PetscErrorCode PCDestroy_BiCGStabCUSP(PC pc)
199: {
200: PetscErrorCode ierr;
203: /*
204: Free the private data structure that was hanging off the PC
205: */
206: PetscFree(pc->data);
207: return(0);
208: }
212: static PetscErrorCode PCSetFromOptions_BiCGStabCUSP(PetscOptionItems *PetscOptionsObject,PC pc)
213: {
214: PC_BiCGStabCUSP *bicg = (PC_BiCGStabCUSP*)pc->data;
215: PetscErrorCode ierr;
218: PetscOptionsHead(PetscOptionsObject,"BiCGStabCUSP options");
219: PetscOptionsReal("-pc_bicgstabcusp_rtol","relative tolerance for BiCGStabCUSP preconditioner","PCBiCGStabCUSPSetTolerance",bicg->rtol,&bicg->rtol,0);
220: PetscOptionsInt("-pc_bicgstabcusp_max_it","maximum iterations for BiCGStabCUSP preconditioner","PCBiCGStabCUSPSetIterations",bicg->maxits,&bicg->maxits,0);
221: PetscOptionsBool("-pc_bicgstabcusp_monitor_verbose","Print information about GPU BiCGStabCUSP iterations","PCBiCGStabCUSPSetUseVerboseMonitor",bicg->monitorverbose,&bicg->monitorverbose,0);
222: PetscOptionsTail();
223: return(0);
224: }
226: /* -------------------------------------------------------------------------- */
230: PETSC_EXTERN PetscErrorCode PCCreate_BiCGStabCUSP(PC pc)
231: {
232: PC_BiCGStabCUSP *bicg;
233: PetscErrorCode ierr;
236: /*
237: Creates the private data structure for this preconditioner and
238: attach it to the PC object.
239: */
240: PetscNewLog(pc,&bicg);
241: /*
242: 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.
243: */
244: bicg->maxits = 1000;
245: bicg->rtol = 1.e-1;
246: bicg->monitorverbose = PETSC_FALSE;
247: pc->data = (void*)bicg;
248: /*
249: Set the pointers for the functions that are provided above.
250: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
251: are called, they will automatically call these functions. Note we
252: choose not to provide a couple of these functions since they are
253: not needed.
254: */
255: pc->ops->apply = PCApply_BiCGStabCUSP;
256: pc->ops->applytranspose = 0;
257: pc->ops->setup = PCSetUp_BiCGStabCUSP;
258: pc->ops->destroy = PCDestroy_BiCGStabCUSP;
259: pc->ops->setfromoptions = PCSetFromOptions_BiCGStabCUSP;
260: pc->ops->view = 0;
261: pc->ops->applyrichardson = 0;
262: pc->ops->applysymmetricleft = 0;
263: pc->ops->applysymmetricright = 0;
265: PetscObjectComposeFunction((PetscObject)pc,"PCBiCGStabCUSPSetTolerance_C",PCBiCGStabCUSPSetTolerance_BiCGStabCUSP);
266: PetscObjectComposeFunction((PetscObject)pc, "PCBiCGStabCUSPSetIterations_C",PCBiCGStabCUSPSetIterations_BiCGStabCUSP);
267: PetscObjectComposeFunction((PetscObject)pc, "PCBiCGStabCUSPSetUseVerboseMonitor_C", PCBiCGStabCUSPSetUseVerboseMonitor_BiCGStabCUSP);
268: return(0);
269: }