Actual source code: ainvcusp.cu
petsc-3.7.3 2016-08-01
2: /* -------------------------------------------------------------------- */
4: /*
5: Include files needed for the CUSP AINV preconditioner:
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: #undef VecType
13: #include <cusp/precond/ainv.h>
14: #define VecType char*
15: #include <../src/vec/vec/impls/dvecimpl.h>
16: #include <../src/mat/impls/aij/seq/seqcusp/cuspmatimpl.h>
17: #include <../src/vec/vec/impls/seq/seqcusp/cuspvecimpl.h>
19: #define cuspainvprecondscaled cusp::precond::scaled_bridson_ainv<PetscScalar,cusp::device_memory>
20: #define cuspainvprecond cusp::precond::bridson_ainv<PetscScalar,cusp::device_memory>
22: /*
23: Private context (data structure) for the CUSP AINV preconditioner. Note that this only works on CUSP SPD matrices.
24: */
25: typedef struct {
26: void *AINVCUSP;
27: PetscBool scaled; /* Whether to use the scaled version of the Bridson AINV or not */
29: PetscInt nonzeros; /* can only use one of nonzeros, droptolerance, linparam at once */
30: PetscReal droptolerance;
31: PetscInt linparam;
32: PetscBool uselin;
33: } PC_AINVCUSP;
35: /* -------------------------------------------------------------------------- */
36: /*
37: PCSetUp_AINVCUSP - Prepares for the use of the CUSP AINV 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_AINVCUSP(PC pc)
52: {
53: PC_AINVCUSP *ainv = (PC_AINVCUSP*)pc->data;
54: 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
63: PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQAIJCUSP,&flg);
64: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Currently only handles CUSP matrices");
65: if (pc->setupcalled != 0) {
66: try {
67: if (ainv->scaled) delete (cuspainvprecondscaled*)ainv->AINVCUSP;
68: else delete (cuspainvprecond*)ainv->AINVCUSP;
69: } catch(char *ex) {
70: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CUSP error: %s", ex);
71: }
72: }
73: try {
74: MatCUSPCopyToGPU(pc->pmat);
75: #if defined(PETSC_USE_COMPLEX)
76: ainv->AINVCUSP = 0;CHKERRQ(1); /* TODO */
77: #else
78: gpustruct = (Mat_SeqAIJCUSP*)(pc->pmat->spptr);
80: if (gpustruct->format==MAT_CUSP_ELL) {
81: CUSPMATRIXELL *mat = (CUSPMATRIXELL*)gpustruct->mat;
82: if (ainv->scaled) ainv->AINVCUSP = new cuspainvprecondscaled(*mat, ainv->droptolerance,ainv->nonzeros,ainv->uselin,ainv->linparam);
83: else ainv->AINVCUSP = new cuspainvprecond(*mat, ainv->droptolerance,ainv->nonzeros,ainv->uselin,ainv->linparam);
84: } else if (gpustruct->format==MAT_CUSP_DIA) {
85: CUSPMATRIXDIA *mat = (CUSPMATRIXDIA*)gpustruct->mat;
86: if (ainv->scaled) ainv->AINVCUSP = new cuspainvprecondscaled(*mat, ainv->droptolerance,ainv->nonzeros,ainv->uselin,ainv->linparam);
87: else ainv->AINVCUSP = new cuspainvprecond(*mat, ainv->droptolerance,ainv->nonzeros,ainv->uselin,ainv->linparam);
88: } else {
89: CUSPMATRIX *mat = (CUSPMATRIX*)gpustruct->mat;
90: if (ainv->scaled) ainv->AINVCUSP = new cuspainvprecondscaled(*mat, ainv->droptolerance,ainv->nonzeros,ainv->uselin,ainv->linparam);
91: else ainv->AINVCUSP = new cuspainvprecond(*mat, ainv->droptolerance,ainv->nonzeros,ainv->uselin,ainv->linparam);
92: }
93: #endif
94: } catch(char *ex) {
95: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CUSP error: %s",ex);
96: }
97: return(0);
98: }
100: /* -------------------------------------------------------------------------- */
101: /*
102: PCApply_AINVCUSP - Applies the CUSP AINV preconditioner to a vector.
104: Input Parameters:
105: . pc - the preconditioner context
106: . x - input vector
108: Output Parameter:
109: . y - output vector
111: Application Interface Routine: PCApply()
112: */
115: static PetscErrorCode PCApply_AINVCUSP(PC pc,Vec x,Vec y)
116: {
117: PC_AINVCUSP *ainv = (PC_AINVCUSP*)pc->data;
119: PetscBool flg1,flg2;
120: CUSPARRAY *xarray=NULL,*yarray=NULL;
123: PetscObjectTypeCompare((PetscObject)x,VECSEQCUSP,&flg1);
124: PetscObjectTypeCompare((PetscObject)y,VECSEQCUSP,&flg2);
125: if (!(flg1 && flg2)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP, "Currently only handles CUSP vectors");
126: if (!ainv->AINVCUSP) {
127: PCSetUp_AINVCUSP(pc);
128: }
129: VecSet(y,0.0);
130: VecCUSPGetArrayRead(x,&xarray);
131: VecCUSPGetArrayWrite(y,&yarray);
132: try {
133: if (ainv->scaled) cusp::multiply(*(cuspainvprecondscaled*)ainv->AINVCUSP,*xarray,*yarray);
134: else cusp::multiply(*(cuspainvprecond*)ainv->AINVCUSP,*xarray,*yarray);
135: } catch(char* ex) {
136: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CUSP error: %s", ex);
137: }
138: VecCUSPRestoreArrayRead(x,&xarray);
139: VecCUSPRestoreArrayWrite(y,&yarray);
140: PetscObjectStateIncrease((PetscObject)y);
141: return(0);
142: }
143: /* -------------------------------------------------------------------------- */
147: static PetscErrorCode PCReset_AINVCUSP(PC pc)
148: {
149: PC_AINVCUSP *ainv = (PC_AINVCUSP*)pc->data;
152: if (ainv->AINVCUSP) {
153: try {
154: if (ainv->scaled) delete (cuspainvprecondscaled*)ainv->AINVCUSP;
155: else delete (cuspainvprecond*)ainv->AINVCUSP;
156: } catch(char* ex) {
157: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CUSP error: %s", ex);
158: }
159: ainv->AINVCUSP = NULL;
160: }
161: return(0);
162: }
164: /*
165: PCDestroy_AINVCUSP - Destroys the private context for the AINVCUSP preconditioner
166: that was created with PCCreate_AINVCUSP().
168: Input Parameter:
169: . pc - the preconditioner context
171: Application Interface Routine: PCDestroy()
172: */
175: static PetscErrorCode PCDestroy_AINVCUSP(PC pc)
176: {
180: PCReset_AINVCUSP(pc);
182: /*
183: Free the private data structure that was hanging off the PC
184: */
185: PetscFree(pc->data);
186: return(0);
187: }
193: static PetscErrorCode PCAINVCUSPSetDropTolerance_AINVCUSP(PC pc, PetscReal droptolerance)
194: {
195: PC_AINVCUSP *ainv = (PC_AINVCUSP*)pc->data;
198: ainv->droptolerance = droptolerance;
199: ainv->uselin = PETSC_FALSE;
200: ainv->linparam = 1;
201: ainv->nonzeros = -1;
202: return(0);
203: }
207: PetscErrorCode PCAINVCUSPSetDropTolerance(PC pc, PetscReal droptolerance)
208: {
213: PetscTryMethod(pc, "PCAINVCUSPSetDropTolerance_C",(PC,PetscReal),(pc,droptolerance));
214: return(0);
215: }
218: static PetscErrorCode PCAINVCUSPSetNonzeros_AINVCUSP(PC pc, PetscInt nonzeros)
219: {
220: PC_AINVCUSP *ainv = (PC_AINVCUSP*)pc->data;
223: ainv->droptolerance = 0;
224: ainv->uselin = PETSC_FALSE;
225: ainv->linparam = 1;
226: ainv->nonzeros = nonzeros;
227: return(0);
228: }
232: PetscErrorCode PCAINVCUSPSetNonzeros(PC pc, PetscInt nonzeros)
233: {
238: PetscTryMethod(pc, "PCAINVCUSPSetNonzeros_C",(PC,PetscInt),(pc,nonzeros));
239: return(0);
240: }
243: static PetscErrorCode PCAINVCUSPSetLinParameter_AINVCUSP(PC pc, PetscInt param)
244: {
245: PC_AINVCUSP *ainv = (PC_AINVCUSP*)pc->data;
248: ainv->droptolerance = 0;
249: ainv->uselin = PETSC_TRUE;
250: ainv->linparam = param;
251: ainv->nonzeros = -1;
252: return(0);
253: }
257: PetscErrorCode PCAINVCUSPSetLinParameter(PC pc, PetscInt param)
258: {
263: PetscTryMethod(pc, "PCAINVCUSPSetLinParameter_C",(PC,PetscInt),(pc,param));
264: return(0);
265: }
268: static PetscErrorCode PCAINVCUSPUseScaling_AINVCUSP(PC pc, PetscBool scaled)
269: {
270: PC_AINVCUSP *ainv = (PC_AINVCUSP*)pc->data;
273: ainv->scaled = scaled;
274: return(0);
275: }
279: PetscErrorCode PCAINVCUSPUseScaling(PC pc, PetscBool scaled)
280: {
285: PetscTryMethod(pc, "PCAINVCUSPUseScaling_C",(PC,PetscBool),(pc,scaled));
286: return(0);
287: }
291: static PetscErrorCode PCSetFromOptions_AINVCUSP(PetscOptionItems *PetscOptionsObject,PC pc)
292: {
293: PC_AINVCUSP *ainv = (PC_AINVCUSP*)pc->data;
294: PetscBool flag = PETSC_FALSE;
298: PetscOptionsHead(PetscOptionsObject,"AINVCUSP options");
299: PetscOptionsReal("-pc_ainvcusp_droptol","drop tolerance for AINVCUSP preconditioner","PCAINVCUSPSetDropTolerance",ainv->droptolerance,&ainv->droptolerance,&flag);
300: if (flag) {
301: ainv->nonzeros = -1;
302: ainv->uselin = PETSC_FALSE;
303: ainv->linparam = 1;
304: flag = PETSC_FALSE;
305: }
306: PetscOptionsInt("-pc_ainvcusp_nonzeros","nonzeros/row for AINVCUSP preconditioner","PCAINVCUSPSetNonzeros",ainv->nonzeros,&ainv->nonzeros,&flag);
307: if (flag) {
308: ainv->droptolerance = 0;
309: ainv->uselin = PETSC_FALSE;
310: ainv->linparam = 1;
311: flag = PETSC_FALSE;
312: }
313: PetscOptionsInt("-pc_ainvcusp_linparameter","Lin parameter for AINVCUSP preconditioner","PCAINVCUSPSetLinParameter",ainv->linparam,&ainv->linparam,&flag);
314: if (flag) {
315: ainv->droptolerance = 0;
316: ainv->uselin = PETSC_TRUE;
317: ainv->droptolerance = 0;
318: ainv->nonzeros = -1;
319: }
320: PetscOptionsBool("-pc_ainvcusp_scale","Whether to use scaled AINVCUSP preconditioner or not","PCAINVCUSPUseScaling",ainv->scaled,&ainv->scaled,0);
321: PetscOptionsTail();
322: return(0);
323: }
325: /* -------------------------------------------------------------------------- */
327: /*MC
328: PCAINVCUSP - A sparse approximate inverse precondition that runs on the Nvidia GPU.
331: http://docs.cusp-library.googlecode.com/hg/classcusp_1_1precond_1_1bridson__ainv.html
333: Level: advanced
335: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC
337: M*/
341: PETSC_EXTERN PetscErrorCode PCCreate_AINVCUSP(PC pc)
342: {
343: PC_AINVCUSP *ainv;
347: /*
348: Creates the private data structure for this preconditioner and
349: attach it to the PC object.
350: */
351: PetscNewLog(pc,&ainv);
352: pc->data = (void*)ainv;
353: ainv->AINVCUSP = 0;
354: ainv->droptolerance = 0.1;
355: ainv->nonzeros = -1;
356: ainv->scaled = PETSC_TRUE;
357: ainv->linparam = 1;
358: ainv->uselin = PETSC_FALSE;
359: /*
360: Set the pointers for the functions that are provided above.
361: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
362: are called, they will automatically call these functions. Note we
363: choose not to provide a couple of these functions since they are
364: not needed.
365: */
366: pc->ops->apply = PCApply_AINVCUSP;
367: pc->ops->applytranspose = 0;
368: pc->ops->setup = PCSetUp_AINVCUSP;
369: pc->ops->reset = PCReset_AINVCUSP;
370: pc->ops->destroy = PCDestroy_AINVCUSP;
371: pc->ops->setfromoptions = PCSetFromOptions_AINVCUSP;
372: pc->ops->view = 0;
373: pc->ops->applyrichardson = 0;
374: pc->ops->applysymmetricleft = 0;
375: pc->ops->applysymmetricright = 0;
377: PetscObjectComposeFunction((PetscObject)pc, "PCAINVCUSPSetDropTolerance_C", PCAINVCUSPSetDropTolerance_AINVCUSP);
378: PetscObjectComposeFunction((PetscObject)pc, "PCAINVCUSPUseScaling_C", PCAINVCUSPUseScaling_AINVCUSP);
379: PetscObjectComposeFunction((PetscObject)pc, "PCAINVCUSPSetLinParameter_C", PCAINVCUSPSetLinParameter_AINVCUSP);
380: PetscObjectComposeFunction((PetscObject)pc, "PCAINVCUSPSetNonzeros_C", PCAINVCUSPSetNonzeros_AINVCUSP);
381: return(0);
382: }