Actual source code: bicgstabcusp.cu

petsc-3.6.1 2015-08-06
Report Typos and Errors
  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: }