Actual source code: sacusppoly.cu
petsc-3.3-p7 2013-05-11
2: /* -------------------------------------------------------------------- */
4: /*
5: Include files needed for the CUSP Smoothed Aggregation preconditioner with Chebyshev polynomial smoothing:
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: #define USE_POLY_SMOOTHER 1
14: #include <cusp/precond/smoothed_aggregation.h>
15: #undef USE_POLY_SMOOTHER
16: #define VecType char*
17: #include <../src/vec/vec/impls/dvecimpl.h>
18: #include <../src/mat/impls/aij/seq/seqcusp/cuspmatimpl.h>
20: #define cuspsaprecond cusp::precond::smoothed_aggregation<PetscInt,PetscScalar,cusp::device_memory>
22: /*
23: Private context (data structure) for the SACUSPPoly preconditioner.
24: */
25: typedef struct {
26: cuspsaprecond* SACUSPPoly;
27: /*int cycles; */
28: } PC_SACUSPPoly;
31: /* -------------------------------------------------------------------------- */
32: /*
33: PCSetUp_SACUSPPoly - Prepares for the use of the SACUSPPoly preconditioner
34: by setting data structures and options.
36: Input Parameter:
37: . pc - the preconditioner context
39: Application Interface Routine: PCSetUp()
41: Notes:
42: The interface routine PCSetUp() is not usually called directly by
43: the user, but instead is called by PCApply() if necessary.
44: */
47: static PetscErrorCode PCSetUp_SACUSPPoly(PC pc)
48: {
49: PC_SACUSPPoly *sa = (PC_SACUSPPoly*)pc->data;
50: PetscBool flg = PETSC_FALSE;
52: Mat_SeqAIJCUSP *gpustruct;
55: PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQAIJCUSP,&flg);
56: if (!flg) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Currently only handles CUSP matrices");
57: if (pc->setupcalled != 0){
58: try {
59: delete sa->SACUSPPoly;
60: } catch(char* ex) {
61: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CUSP error: %s", ex);
62: }
63: }
64: try {
65: MatCUSPCopyToGPU(pc->pmat);
66: gpustruct = (Mat_SeqAIJCUSP *)(pc->pmat->spptr);
67: sa->SACUSPPoly = new cuspsaprecond(*(CUSPMATRIX*)gpustruct->mat);
68: } catch(char* ex) {
69: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CUSP error: %s", ex);
70: }
71: /*PetscOptionsInt("-pc_sacusp_cycles","Number of v-cycles to perform","PCSACUSPSetCycles",sa->cycles,
72: &sa->cycles,PETSC_NULL);*/
73: return(0);
74: }
78: static PetscErrorCode PCApplyRichardson_SACUSPPoly(PC pc, Vec b, Vec y, Vec w,PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt its, PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
79: {
80: PC_SACUSPPoly *sac = (PC_SACUSPPoly*)pc->data;
82: CUSPARRAY *barray,*yarray;
83:
85: /* how to incorporate dtol, guesszero, w?*/
87: VecCUSPGetArrayRead(b,&barray);
88: VecCUSPGetArrayReadWrite(y,&yarray);
89: cusp::default_monitor<PetscScalar> monitor(*barray,its,rtol,abstol);
90: sac->SACUSPPoly->solve(*barray,*yarray,monitor);
91: *outits = monitor.iteration_count();
92: if (monitor.converged()){
93: /* how to discern between converging from RTOL or ATOL?*/
94: *reason = PCRICHARDSON_CONVERGED_RTOL;
95: } else{
96: *reason = PCRICHARDSON_CONVERGED_ITS;
97: }
98: PetscObjectStateIncrease((PetscObject)y);
99: VecCUSPRestoreArrayRead(b,&barray);
100: VecCUSPRestoreArrayReadWrite(y,&yarray);
101: return(0);
102: }
104: /* -------------------------------------------------------------------------- */
105: /*
106: PCApply_SACUSPPoly - Applies the SACUSPPoly preconditioner to a vector.
108: Input Parameters:
109: . pc - the preconditioner context
110: . x - input vector
112: Output Parameter:
113: . y - output vector
115: Application Interface Routine: PCApply()
116: */
119: static PetscErrorCode PCApply_SACUSPPoly(PC pc,Vec x,Vec y)
120: {
121: PC_SACUSPPoly *sac = (PC_SACUSPPoly*)pc->data;
123: PetscBool flg1,flg2;
124: CUSPARRAY *xarray,*yarray;
127: /*how to apply a certain fixed number of iterations?*/
128: PetscObjectTypeCompare((PetscObject)x,VECSEQCUSP,&flg1);
129: PetscObjectTypeCompare((PetscObject)y,VECSEQCUSP,&flg2);
130: if (!(flg1 && flg2)) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_SUP, "Currently only handles CUSP vectors");
131: if (!sac->SACUSPPoly) {
132: PCSetUp_SACUSPPoly(pc);
133: }
134: VecSet(y,0.0);
135: VecCUSPGetArrayRead(x,&xarray);
136: VecCUSPGetArrayWrite(y,&yarray);
137: try {
138: cusp::multiply(*sac->SACUSPPoly,*xarray,*yarray);
139: } catch(char* ex) {
140: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CUSP error: %s", ex);
141: }
142: VecCUSPRestoreArrayRead(x,&xarray);
143: VecCUSPRestoreArrayWrite(y,&yarray);
144: PetscObjectStateIncrease((PetscObject)y);
145: return(0);
146: }
147: /* -------------------------------------------------------------------------- */
148: /*
149: PCDestroy_SACUSPPoly - Destroys the private context for the SACUSPPoly preconditioner
150: that was created with PCCreate_SACUSPPoly().
152: Input Parameter:
153: . pc - the preconditioner context
155: Application Interface Routine: PCDestroy()
156: */
159: static PetscErrorCode PCDestroy_SACUSPPoly(PC pc)
160: {
161: PC_SACUSPPoly *sac = (PC_SACUSPPoly*)pc->data;
165: if (sac->SACUSPPoly) {
166: try {
167: delete sac->SACUSPPoly;
168: } catch(char* ex) {
169: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CUSP error: %s", ex);
170: }
171: }
173: /*
174: Free the private data structure that was hanging off the PC
175: */
176: PetscFree(sac);
177: return(0);
178: }
182: static PetscErrorCode PCSetFromOptions_SACUSPPoly(PC pc)
183: {
187: PetscOptionsHead("SACUSPPoly options");
188: PetscOptionsTail();
189: return(0);
190: }
192: /* -------------------------------------------------------------------------- */
196: EXTERN_C_BEGIN
199: PetscErrorCode PCCreate_SACUSPPoly(PC pc)
200: {
201: PC_SACUSPPoly *sac;
205: /*
206: Creates the private data structure for this preconditioner and
207: attach it to the PC object.
208: */
209: PetscNewLog(pc,PC_SACUSPPoly,&sac);
210: pc->data = (void*)sac;
212: /*
213: Initialize the pointer to zero
214: Initialize number of v-cycles to default (1)
215: */
216: sac->SACUSPPoly = 0;
217: /*sac->cycles=1;*/
220: /*
221: Set the pointers for the functions that are provided above.
222: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
223: are called, they will automatically call these functions. Note we
224: choose not to provide a couple of these functions since they are
225: not needed.
226: */
227: pc->ops->apply = PCApply_SACUSPPoly;
228: pc->ops->applytranspose = 0;
229: pc->ops->setup = PCSetUp_SACUSPPoly;
230: pc->ops->destroy = PCDestroy_SACUSPPoly;
231: pc->ops->setfromoptions = PCSetFromOptions_SACUSPPoly;
232: pc->ops->view = 0;
233: pc->ops->applyrichardson = PCApplyRichardson_SACUSPPoly;
234: pc->ops->applysymmetricleft = 0;
235: pc->ops->applysymmetricright = 0;
236: return(0);
237: }
238: EXTERN_C_END