Actual source code: ainvcusp.cu
petsc-3.6.1 2015-08-06
2: /* -------------------------------------------------------------------- */
4: /*
5: Include files needed for the CUSP AINV 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: #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>
18: #define cuspainvprecondscaled cusp::precond::scaled_bridson_ainv<PetscScalar,cusp::device_memory>
19: #define cuspainvprecond cusp::precond::bridson_ainv<PetscScalar,cusp::device_memory>
21: /*
22: Private context (data structure) for the CUSP AINV preconditioner. Note that this only works on CUSP SPD matrices.
23: */
24: typedef struct {
25: void *AINVCUSP;
26: PetscBool scaled; /* Whether to use the scaled version of the Bridson AINV or not */
28: PetscInt nonzeros; /* can only use one of nonzeros, droptolerance, linparam at once */
29: PetscReal droptolerance;
30: PetscInt linparam;
31: PetscBool uselin;
32: } PC_AINVCUSP;
34: /* -------------------------------------------------------------------------- */
35: /*
36: PCSetUp_AINVCUSP - Prepares for the use of the CUSP AINV 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_AINVCUSP(PC pc)
51: {
52: PC_AINVCUSP *ainv = (PC_AINVCUSP*)pc->data;
53: PetscBool flg = PETSC_FALSE;
54: #if !defined(PETSC_USE_COMPLEX)
55: // protect these in order to avoid compiler warnings. This preconditioner does
56: // not work for complex types.
57: Mat_SeqAIJCUSP *gpustruct;
58: #endif
62: PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQAIJCUSP,&flg);
63: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Currently only handles CUSP matrices");
64: if (pc->setupcalled != 0) {
65: try {
66: if (ainv->scaled) delete (cuspainvprecondscaled*)ainv->AINVCUSP;
67: else delete (cuspainvprecond*)ainv->AINVCUSP;
68: } catch(char *ex) {
69: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CUSP error: %s", ex);
70: }
71: }
72: try {
73: MatCUSPCopyToGPU(pc->pmat);
74: #if defined(PETSC_USE_COMPLEX)
75: ainv->AINVCUSP = 0;CHKERRQ(1); /* TODO */
76: #else
77: gpustruct = (Mat_SeqAIJCUSP*)(pc->pmat->spptr);
79: if (gpustruct->format==MAT_CUSP_ELL) {
80: CUSPMATRIXELL *mat = (CUSPMATRIXELL*)gpustruct->mat;
81: if (ainv->scaled) ainv->AINVCUSP = new cuspainvprecondscaled(*mat, ainv->droptolerance,ainv->nonzeros,ainv->uselin,ainv->linparam);
82: else ainv->AINVCUSP = new cuspainvprecond(*mat, ainv->droptolerance,ainv->nonzeros,ainv->uselin,ainv->linparam);
83: } else if (gpustruct->format==MAT_CUSP_DIA) {
84: CUSPMATRIXDIA *mat = (CUSPMATRIXDIA*)gpustruct->mat;
85: if (ainv->scaled) ainv->AINVCUSP = new cuspainvprecondscaled(*mat, ainv->droptolerance,ainv->nonzeros,ainv->uselin,ainv->linparam);
86: else ainv->AINVCUSP = new cuspainvprecond(*mat, ainv->droptolerance,ainv->nonzeros,ainv->uselin,ainv->linparam);
87: } else {
88: CUSPMATRIX *mat = (CUSPMATRIX*)gpustruct->mat;
89: if (ainv->scaled) ainv->AINVCUSP = new cuspainvprecondscaled(*mat, ainv->droptolerance,ainv->nonzeros,ainv->uselin,ainv->linparam);
90: else ainv->AINVCUSP = new cuspainvprecond(*mat, ainv->droptolerance,ainv->nonzeros,ainv->uselin,ainv->linparam);
91: }
92: #endif
93: } catch(char *ex) {
94: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CUSP error: %s",ex);
95: }
96: return(0);
97: }
99: /* -------------------------------------------------------------------------- */
100: /*
101: PCApply_AINVCUSP - Applies the CUSP AINV preconditioner to a vector.
103: Input Parameters:
104: . pc - the preconditioner context
105: . x - input vector
107: Output Parameter:
108: . y - output vector
110: Application Interface Routine: PCApply()
111: */
114: static PetscErrorCode PCApply_AINVCUSP(PC pc,Vec x,Vec y)
115: {
116: PC_AINVCUSP *ainv = (PC_AINVCUSP*)pc->data;
118: PetscBool flg1,flg2;
119: CUSPARRAY *xarray=NULL,*yarray=NULL;
122: PetscObjectTypeCompare((PetscObject)x,VECSEQCUSP,&flg1);
123: PetscObjectTypeCompare((PetscObject)y,VECSEQCUSP,&flg2);
124: if (!(flg1 && flg2)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP, "Currently only handles CUSP vectors");
125: if (!ainv->AINVCUSP) {
126: PCSetUp_AINVCUSP(pc);
127: }
128: VecSet(y,0.0);
129: VecCUSPGetArrayRead(x,&xarray);
130: VecCUSPGetArrayWrite(y,&yarray);
131: try {
132: if (ainv->scaled) cusp::multiply(*(cuspainvprecondscaled*)ainv->AINVCUSP,*xarray,*yarray);
133: else cusp::multiply(*(cuspainvprecond*)ainv->AINVCUSP,*xarray,*yarray);
134: } catch(char* ex) {
135: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CUSP error: %s", ex);
136: }
137: VecCUSPRestoreArrayRead(x,&xarray);
138: VecCUSPRestoreArrayWrite(y,&yarray);
139: PetscObjectStateIncrease((PetscObject)y);
140: return(0);
141: }
142: /* -------------------------------------------------------------------------- */
146: static PetscErrorCode PCReset_AINVCUSP(PC pc)
147: {
148: PC_AINVCUSP *ainv = (PC_AINVCUSP*)pc->data;
151: if (ainv->AINVCUSP) {
152: try {
153: if (ainv->scaled) delete (cuspainvprecondscaled*)ainv->AINVCUSP;
154: else delete (cuspainvprecond*)ainv->AINVCUSP;
155: } catch(char* ex) {
156: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CUSP error: %s", ex);
157: }
158: ainv->AINVCUSP = NULL;
159: }
160: return(0);
161: }
163: /*
164: PCDestroy_AINVCUSP - Destroys the private context for the AINVCUSP preconditioner
165: that was created with PCCreate_AINVCUSP().
167: Input Parameter:
168: . pc - the preconditioner context
170: Application Interface Routine: PCDestroy()
171: */
174: static PetscErrorCode PCDestroy_AINVCUSP(PC pc)
175: {
179: PCReset_AINVCUSP(pc);
181: /*
182: Free the private data structure that was hanging off the PC
183: */
184: PetscFree(pc->data);
185: return(0);
186: }
192: static PetscErrorCode PCAINVCUSPSetDropTolerance_AINVCUSP(PC pc, PetscReal droptolerance)
193: {
194: PC_AINVCUSP *ainv = (PC_AINVCUSP*)pc->data;
197: ainv->droptolerance = droptolerance;
198: ainv->uselin = PETSC_FALSE;
199: ainv->linparam = 1;
200: ainv->nonzeros = -1;
201: return(0);
202: }
206: PetscErrorCode PCAINVCUSPSetDropTolerance(PC pc, PetscReal droptolerance)
207: {
212: PetscTryMethod(pc, "PCAINVCUSPSetDropTolerance_C",(PC,PetscReal),(pc,droptolerance));
213: return(0);
214: }
217: static PetscErrorCode PCAINVCUSPSetNonzeros_AINVCUSP(PC pc, PetscInt nonzeros)
218: {
219: PC_AINVCUSP *ainv = (PC_AINVCUSP*)pc->data;
222: ainv->droptolerance = 0;
223: ainv->uselin = PETSC_FALSE;
224: ainv->linparam = 1;
225: ainv->nonzeros = nonzeros;
226: return(0);
227: }
231: PetscErrorCode PCAINVCUSPSetNonzeros(PC pc, PetscInt nonzeros)
232: {
237: PetscTryMethod(pc, "PCAINVCUSPSetNonzeros_C",(PC,PetscInt),(pc,nonzeros));
238: return(0);
239: }
242: static PetscErrorCode PCAINVCUSPSetLinParameter_AINVCUSP(PC pc, PetscInt param)
243: {
244: PC_AINVCUSP *ainv = (PC_AINVCUSP*)pc->data;
247: ainv->droptolerance = 0;
248: ainv->uselin = PETSC_TRUE;
249: ainv->linparam = param;
250: ainv->nonzeros = -1;
251: return(0);
252: }
256: PetscErrorCode PCAINVCUSPSetLinParameter(PC pc, PetscInt param)
257: {
262: PetscTryMethod(pc, "PCAINVCUSPSetLinParameter_C",(PC,PetscInt),(pc,param));
263: return(0);
264: }
267: static PetscErrorCode PCAINVCUSPUseScaling_AINVCUSP(PC pc, PetscBool scaled)
268: {
269: PC_AINVCUSP *ainv = (PC_AINVCUSP*)pc->data;
272: ainv->scaled = scaled;
273: return(0);
274: }
278: PetscErrorCode PCAINVCUSPUseScaling(PC pc, PetscBool scaled)
279: {
284: PetscTryMethod(pc, "PCAINVCUSPUseScaling_C",(PC,PetscBool),(pc,scaled));
285: return(0);
286: }
290: static PetscErrorCode PCSetFromOptions_AINVCUSP(PetscOptions *PetscOptionsObject,PC pc)
291: {
292: PC_AINVCUSP *ainv = (PC_AINVCUSP*)pc->data;
293: PetscBool flag = PETSC_FALSE;
297: PetscOptionsHead(PetscOptionsObject,"AINVCUSP options");
298: PetscOptionsReal("-pc_ainvcusp_droptol","drop tolerance for AINVCUSP preconditioner","PCAINVCUSPSetDropTolerance",ainv->droptolerance,&ainv->droptolerance,&flag);
299: if (flag) {
300: ainv->nonzeros = -1;
301: ainv->uselin = PETSC_FALSE;
302: ainv->linparam = 1;
303: flag = PETSC_FALSE;
304: }
305: PetscOptionsInt("-pc_ainvcusp_nonzeros","nonzeros/row for AINVCUSP preconditioner","PCAINVCUSPSetNonzeros",ainv->nonzeros,&ainv->nonzeros,&flag);
306: if (flag) {
307: ainv->droptolerance = 0;
308: ainv->uselin = PETSC_FALSE;
309: ainv->linparam = 1;
310: flag = PETSC_FALSE;
311: }
312: PetscOptionsInt("-pc_ainvcusp_linparameter","Lin parameter for AINVCUSP preconditioner","PCAINVCUSPSetLinParameter",ainv->linparam,&ainv->linparam,&flag);
313: if (flag) {
314: ainv->droptolerance = 0;
315: ainv->uselin = PETSC_TRUE;
316: ainv->droptolerance = 0;
317: ainv->nonzeros = -1;
318: }
319: PetscOptionsBool("-pc_ainvcusp_scale","Whether to use scaled AINVCUSP preconditioner or not","PCAINVCUSPUseScaling",ainv->scaled,&ainv->scaled,0);
320: PetscOptionsTail();
321: return(0);
322: }
324: /* -------------------------------------------------------------------------- */
326: /*MC
327: PCAINVCUSP - A sparse approximate inverse precondition that runs on the Nvidia GPU.
330: http://docs.cusp-library.googlecode.com/hg/classcusp_1_1precond_1_1bridson__ainv.html
332: Level: advanced
334: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC
336: M*/
340: PETSC_EXTERN PetscErrorCode PCCreate_AINVCUSP(PC pc)
341: {
342: PC_AINVCUSP *ainv;
346: /*
347: Creates the private data structure for this preconditioner and
348: attach it to the PC object.
349: */
350: PetscNewLog(pc,&ainv);
351: pc->data = (void*)ainv;
352: ainv->AINVCUSP = 0;
353: ainv->droptolerance = 0.1;
354: ainv->nonzeros = -1;
355: ainv->scaled = PETSC_TRUE;
356: ainv->linparam = 1;
357: ainv->uselin = PETSC_FALSE;
358: /*
359: Set the pointers for the functions that are provided above.
360: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
361: are called, they will automatically call these functions. Note we
362: choose not to provide a couple of these functions since they are
363: not needed.
364: */
365: pc->ops->apply = PCApply_AINVCUSP;
366: pc->ops->applytranspose = 0;
367: pc->ops->setup = PCSetUp_AINVCUSP;
368: pc->ops->reset = PCReset_AINVCUSP;
369: pc->ops->destroy = PCDestroy_AINVCUSP;
370: pc->ops->setfromoptions = PCSetFromOptions_AINVCUSP;
371: pc->ops->view = 0;
372: pc->ops->applyrichardson = 0;
373: pc->ops->applysymmetricleft = 0;
374: pc->ops->applysymmetricright = 0;
376: PetscObjectComposeFunction((PetscObject)pc, "PCAINVCUSPSetDropTolerance_C", PCAINVCUSPSetDropTolerance_AINVCUSP);
377: PetscObjectComposeFunction((PetscObject)pc, "PCAINVCUSPUseScaling_C", PCAINVCUSPUseScaling_AINVCUSP);
378: PetscObjectComposeFunction((PetscObject)pc, "PCAINVCUSPSetLinParameter_C", PCAINVCUSPSetLinParameter_AINVCUSP);
379: PetscObjectComposeFunction((PetscObject)pc, "PCAINVCUSPSetNonzeros_C", PCAINVCUSPSetNonzeros_AINVCUSP);
380: return(0);
381: }