Actual source code: daapp_mgrid.c
1: #include "taodaapplication.h" /* taodaapplication.h */
2: //#include "taoapp.h"
3: #include src/petsctao/include/taopetsc.h
5: #include "../../application/taoapp/taoapp_petsc.h"
6: #include "../interface/daapp_impl.h"
7: #include multigridmat.h
9: #include "petscmg.h"
11: class TaoMultiGridApplication: public TaoPetscApplication {
13: protected:
15: public:
17: TaoMultiGridApplication(MPI_Comm);
18: TaoMultiGridApplication();
20: int coarselevel;
21: int EvaluateHessian(TaoVec *xx, TaoMat *HH);
23: };
28: TaoMultiGridApplication :: TaoMultiGridApplication(MPI_Comm mpicomm) : TaoPetscApplication(mpicomm) {
29: Mat M=0;
30: this->coarselevel=0;
31: this->taoh=new TaoMultiGridMatPetsc(M);
32: this->taoj=new TaoMultiGridMatPetsc(M);
33: return;
34: }
38: TaoMultiGridApplication::TaoMultiGridApplication() : TaoPetscApplication(PETSC_COMM_WORLD){
39:
40: Mat M=0;
41: this->coarselevel=0;
42: this->taoh=new TaoMultiGridMatPetsc(M);
43: this->taoj=new TaoMultiGridMatPetsc(M);
44: return;
45: }
50: int TaoMultiGridApplication::EvaluateHessian(TaoVec *xx, TaoMat *HH){
51: int info,i,clevel,cclevel;
52: Vec X = (Vec) xx->VecObject;
53: TaoMatPetsc *MM=(TaoMatPetsc*)HH;
54: Mat H, HPre;
55: MatStructure pflag;
56: DA_APPLICATION daapp;
59: info=TaoAppDAApp(this->papp,&daapp); CHKERRQ(info);
61: PetscStackPush("TAO User DA Hessian");
62: info=MM->GetMatrix(&H,&HPre,&pflag);CHKERRQ(info);
63: info=DAAppGetCurrentLevel(this->papp,&clevel); CHKERRQ(info);
64: cclevel=PetscMin(clevel,this->coarselevel);
65: cclevel=PetscMax(cclevel,0);
67: for (i=clevel; i>=cclevel; i--){
68: daapp->currentlevel=i;
69: info=TaoAppSetColoring(this->papp,daapp->grid[i].coloring); CHKERRQ(info);
70: info = PetscLogInfo((daapp,"TAO hessian evaluation at DA_APPLICATION object, level %d.\n",i)); CHKERRQ(info);
71: info = TaoAppComputeHessian(this->papp, X, &H, &HPre, &pflag); CHKERRQ(info);
72: if (i==clevel){
73: info = MM->SetMatrix(H,HPre,pflag);CHKERRQ(info);
74: }
76: if (i>cclevel){
77: info=MatMultTranspose(daapp->grid[i].Interpolate,X,daapp->grid[i-1].R);CHKERRQ(info);
78: info=VecPointwiseMult(daapp->grid[i-1].R,daapp->grid[i].CScale,daapp->grid[i-1].R);CHKERRQ(info);
79: X=daapp->grid[i-1].R;
80: H=daapp->grid[i-1].H;
81: HPre=H;
82: }
83: }
84: daapp->currentlevel=clevel;
85: info=TaoAppSetColoring(this->papp,daapp->grid[clevel].coloring); CHKERRQ(info);
86: PetscStackPop;
88: return(0);
89: }
91: int DAAppSetMultigridKSP(GridCtx *, int, KSP);
92: int DAAppSetupMultigridMonitor(TAO_APPLICATION, DA, int, void *);
93: int TaoAppGetMultiGridApplication(TAO_APPLICATION, TaoMultiGridApplication *);
94: static char TaoPetscPCMGDAAppXXX[] = "TaoPetscPCMGDAApp";
95: static char TaoPetscAppXXX[] = "TaoPetscApp";
100: static int DAAppDestroyTaoAppXXX(void*ctx){
101: TaoMultiGridApplication *mctx=(TaoMultiGridApplication*)ctx;
103: if (mctx){ delete mctx;}
104: return(0);
105: }
109: int TaoAppGetMultiGridApplication(TAO_APPLICATION daapplication, TaoMultiGridApplication **mgdaapp){
110: int ii,clevel,info;
111: Mat H,HH;
112: MatStructure pflag=SAME_NONZERO_PATTERN;
113: DA_APPLICATION daapp;
114: TaoMultiGridApplication* daapp2;
115: TaoPetscApplication *pppappp;
119: info=TaoAppQueryForObject(daapplication,TaoPetscPCMGDAAppXXX,(void**)&daapp2); CHKERRQ(info);
120: if (daapp2==0){
121: daapp2=new TaoMultiGridApplication(PETSC_COMM_WORLD);
122: daapp2->papp=daapplication;
123: info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
125: info=TaoAppQueryForObject(daapplication,TaoPetscAppXXX,(void**)&pppappp); CHKERRQ(info);
126: if (pppappp && pppappp->taoh && 0==1){
127: info=pppappp->taoh->GetMatrix(&H,&HH,&pflag); CHKERRQ(info);
128: info=daapp2->taoh->SetMatrix(H,HH,pflag);CHKERRQ(info);
129: } else {
130: info=DAAppGetCurrentLevel(daapplication,&clevel); CHKERRQ(info);
131: info=daapp2->taoh->SetMatrix(daapp->grid[clevel].H,daapp->grid[clevel].H,pflag);CHKERRQ(info);
132: }
133: if (pppappp){
134: daapp2->tao=pppappp->tao;
135: }
136: info=TaoAppAddObject(daapplication,TaoPetscPCMGDAAppXXX,(void*)daapp2,&ii); CHKERRQ(info);
137: info=TaoAppQueryRemoveObject(daapplication,TaoPetscAppXXX); CHKERRQ(info);
138: info=TaoAppAddObject(daapplication,TaoPetscAppXXX,(void*)daapp2,&ii); CHKERRQ(info);
139: info=TaoAppSetDestroyRoutine(daapplication,DAAppDestroyTaoAppXXX, (void*)daapp2); CHKERRQ(info);
140: }
141: *mgdaapp=daapp2;
142: return(0);
143: }
145: static int UsePCMGGG=0;
148: /*@
149: DAAppUseMultigrid - Set the preconditioner for the linear solver to be an algebraic multigrid.
151: Collective on TAO_APPLICATION
152:
153: Input Parameters:
154: + daapplication - the DA Application object
155: - coarselevel - the coarsest grid to be used in the multigrid preconditioner. (Grid 0 is the coarsest grid.
157: Level: intermediate
159: Options Database Key:
160: + -tao_da_multigrid - use multigrid linear solver
161: - -ksp_view - view the linear solver
163: Note:
164: This function should be called after DAAppSetHessianRoutine();
166: Note:
167: This function should be called before each optimization solver as part of the DAAppMonitor
169: Note:
170: Multigrid functionality is still under developement for good performance.
172: .seealso: TaoGetKSP(), DAAppSetupMultigrid()
174: Options Database Key:
175: . -tao_da_multigrid
177: .keywords: Linear Solver, Multigrid, DA, KSP
179: @*/
180: int DAAppUseMultigrid(TAO_APPLICATION daapplication, int coarselevel){
181: int info;
183: UsePCMGGG=coarselevel;
184: info=DAAppSetBeforeMonitor(daapplication,DAAppSetupMultigridMonitor,0);
185: CHKERRQ(info);
186: return(0);
187: }
192: int DAAppSetupMultigridMonitor(TAO_APPLICATION daapplication, DA da, int level, void *ctx){
193: int info;
195: info = DAAppSetupMultigrid(daapplication,UsePCMGGG); CHKERRQ(info);
196: return(0);
197: }
202: /*@
203: DAAppSetupMultigrid - Sets up the multigrid preconditioner.
205: Collective on TAO_APPLICATION
207: Input Parameters:
208: + daapplication - the TAO_APPLICATION context
209: - coarselevel - the coarsest grid that should be used in the multigrid (>=0)
211: Note:
212: Usually the coarselevel is set to 0;
214: Level: intermediate
216: .seealso: DAAppUseMultigrid(), TaoAppSetMonitor()
217: @*/
218: int DAAppSetupMultigrid(TAO_APPLICATION daapplication, int coarselevel){
219: int nn,level,info;
220: TaoMultiGridMatPetsc *MMM;
221: DA_APPLICATION daapp;
222: TaoMultiGridApplication *mgdaapp;
223: TAO_SOLVER tao;
224: KSP ksp;
225:
227: info=DAAppGetCurrentLevel(daapplication,&level); CHKERRQ(info);
228: if (coarselevel<level){
229: info=TaoAppGetMultiGridApplication(daapplication,&mgdaapp); CHKERRQ(info);
230: info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
231: mgdaapp->coarselevel=PetscMax(0,coarselevel);
232: nn=level-coarselevel+1;
233: if ( mgdaapp->taoh ){
234: MMM=(TaoMultiGridMatPetsc*)(mgdaapp->taoh);
235: info=MMM->SetUpMultiGrid(daapp->grid+coarselevel,nn); CHKERRQ(info);
236: } else if ( mgdaapp->taoj ){
237: MMM=(TaoMultiGridMatPetsc*)(mgdaapp->taoj);
238: info=MMM->SetUpMultiGrid(daapp->grid+coarselevel,nn); CHKERRQ(info);
239: mgdaapp->taoj=MMM;
240: }
242: info=TaoAppGetTaoSolver(daapplication,&tao); CHKERRQ(info);
243: info=TaoSetupApplicationSolver(daapplication,tao); CHKERRQ(info);
244: info=TaoAppGetKSP(daapplication,&ksp); CHKERRQ(info);
245: if (ksp){
247: info = DAAppSetMultigridKSP(daapp->grid+coarselevel,nn,ksp);CHKERRQ(info);
248: }
249: }
250: return(0);
251: }
254:
257: int DAAppSetMultigridKSP(GridCtx* dagrid, int nlevels, KSP ksp){
258: int i,info;
259: PC pc;
260: KSP sksp;
261: KSPType ksptype;
264: if (ksp==0 || nlevels<=1){
265: return(0);
266: }
269: info = PetscLogInfo((ksp,"Set multigrid precondition in optimization solver.\n")); CHKERRQ(info);
270:
271: info = KSPGetType(ksp,&ksptype); CHKERRQ(info);
273: info = KSPGetPC(ksp,&pc); CHKERRQ(info);
274: info = PCSetType(pc,PCMG); CHKERRQ(info);
275:
276: info = PCMGSetLevels(pc,nlevels,PETSC_NULL); CHKERRQ(info);
277:
278: for (i=0; i<nlevels; i++) {
279:
280: info = PCMGGetSmoother(pc,i,&sksp); CHKERRQ(info);
281: info = KSPSetType(sksp,ksptype); CHKERRQ(info);
283: info = KSPSetOperators(sksp,dagrid[i].H,dagrid[i].H,SAME_NONZERO_PATTERN); CHKERRQ(info);
284:
285: info = PCMGSetX(pc,i,dagrid[i].X); CHKERRQ(info);
286: info = PCMGSetRhs(pc,i,dagrid[i].RHS); CHKERRQ(info);
287: info = PCMGSetR(pc,i,dagrid[i].R); CHKERRQ(info);
288: info = PCMGSetResidual(pc,i,PCMGDefaultResidual,dagrid[i].H); CHKERRQ(info);
289: }
290: for (i=1; i<nlevels; i++) {
291: info = PCMGSetInterpolate(pc,i,dagrid[i].Interpolate); CHKERRQ(info);
292: info = PCMGSetRestriction(pc,i,dagrid[i].Interpolate); CHKERRQ(info);
293: }
294:
295: return(0);
296: }
300: /* @
301: DAAppSetQuadraticObjective - identify the objective function as quadratic or not.
303: Collective on TAO_APPLICATION
305: Input Parameters:
306: + daapplication - the DA Application object
307: - flag - indicates whether the objective is quadratic or not
309: Level: intermediate
311: Options Database Key:
312: . -tao_da_quadratic - use multigrid linear solver
314: Note:
315: The default value is PETSC_FALSE.
317: Note:
318: If quadratic, consider setting the flag used in KSP to SAME_PRECONDITIONER
320: Note:
321: If quadratic, this routine reduces the number of Hessian evaluations done when using
322: the multigrid preconditioner.
324: .seealso: DAAppSetMatStructure(), DAAppUseMultigrid()
327: .keywords: DA, Objective Function
329: @ */
330: int DAAppSetQuadraticObjective(TAO_APPLICATION daapplication, PetscTruth pflag){
331: int info;
332: TaoMultiGridApplication *mgdaapp;
335: info=TaoAppGetMultiGridApplication(daapplication,&mgdaapp); CHKERRQ(info);
336: // mgdaapp->IsQuadratic=pflag;
337: if (pflag==PETSC_TRUE){
338: info = PetscLogInfo((daapplication,"User labels the objective function as quadratic.\n")); CHKERRQ(info);
339: info = PetscLogInfo((daapplication,"Set KSP MatStructure to SAME_PRECONDITIONER.\n")); CHKERRQ(info);
340: info=DAAppSetMatStructure(daapplication, SAME_PRECONDITIONER);CHKERRQ(info);
341: } else {
342: info = PetscLogInfo((daapplication,"User labels the objective function as NOT quadratic.\n")); CHKERRQ(info);
343: info = PetscLogInfo((daapplication,"Set KSP MatStructure to SAME_NONZERO_PATTERN.\n")); CHKERRQ(info);
344: info=DAAppSetMatStructure(daapplication, SAME_NONZERO_PATTERN);CHKERRQ(info);
345: }
346: return(0);
347: }
352: /* @
353: DAAppSetMultiGridOptions - Sets various multigrid options to be used in this application
354: and the TAO solver.
356: Collective on TAO_APPLICATION
358: Input Parameters:
359: . daapplication - the DA Application object
361: Level: beginner
363: .keywords: options
365: .seealso: TaoDAAppSetOptions();
367: @ */
368: int DAAppSetMultiGridOptions(TAO_APPLICATION daapplication){
369: int info,nlevel;
371: PetscTruth flg1=PETSC_FALSE;
374: info = PetscLogInfo((daapplication,"TaoDAAppSetMultiGridOptions(): Reading command line for options\n")); CHKERRQ(info);
376: flg1=PETSC_FALSE;
377: info = PetscOptionsInt("-tao_da_multigrid","use multigrid linear solver","DAAppUseMultigrid",
378: PETSC_FALSE,&nlevel,&flg1);CHKERRQ(info);
379: if (flg1) {
380: info=DAAppUseMultigrid(daapplication,nlevel);CHKERRQ(info);
381: info = PetscLogInfo((daapplication,"TaoDAAppSetOptions: Use Multigrid linear solver \n")); CHKERRQ(info);
382: }
383: /*
384: flg1=PETSC_FALSE,flg2=PETSC_FALSE;
385: info = PetscOptionsTruth("-tao_da_quadratic","Identify the objective function as quadratic",
386: "DAAppSetQuadraticObjective",PETSC_FALSE,&flg2,&flg1);CHKERRQ(info);
387: if (flg1) {
388: info = DAAppSetQuadraticObjective(daapplication, flg2);CHKERRQ(info);
389: }
390: */
392: return(0);
393: }