Actual source code: daapp_mgrid.c
1: #include "tao_general.h"
2: #include "taodaapplication.h" /* taodaapplication.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: PetscInt 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: {
52: TaoVecPetsc *px = dynamic_cast <TaoVecPetsc *> (xx);
53: TaoMatPetsc *pH = dynamic_cast <TaoMatPetsc *> (HH);
55: int info;
56: PetscInt i,clevel,cclevel;
57: Vec X = px->GetVec();
58: Mat H, HPre;
59: MatStructure pflag;
60: DA_APPLICATION daapp;
63: info=TaoAppDAApp(this->papp,&daapp); CHKERRQ(info);
65: PetscStackPush("TAO User DA Hessian");
66: info=pH->GetMatrix(&H,&HPre,&pflag);CHKERRQ(info);
67: info=DAAppGetCurrentLevel(this->papp,&clevel); CHKERRQ(info);
68: cclevel=PetscMin(clevel,this->coarselevel);
69: cclevel=PetscMax(cclevel,0);
71: for (i=clevel; i>=cclevel; i--){
72: daapp->currentlevel=i;
73: info=TaoAppSetColoring(this->papp,daapp->grid[i].coloring); CHKERRQ(info);
74: info = PetscInfo1(daapp,"TAO hessian evaluation at DA_APPLICATION object, level %d.\n",i); CHKERRQ(info);
75: info = TaoAppComputeHessian(this->papp, X, &H, &HPre, &pflag); CHKERRQ(info);
76: if (i==clevel){
77: info = pH->SetMatrix(H,HPre,pflag);CHKERRQ(info);
78: }
80: if (i>cclevel){
81: info=MatMultTranspose(daapp->grid[i].Interpolate,X,daapp->grid[i-1].R);CHKERRQ(info);
82: info=VecPointwiseMult(daapp->grid[i-1].R,daapp->grid[i].CScale,daapp->grid[i-1].R);CHKERRQ(info);
83: X=daapp->grid[i-1].R;
84: H=daapp->grid[i-1].H;
85: HPre=H;
86: }
87: }
88: daapp->currentlevel=clevel;
89: info=TaoAppSetColoring(this->papp,daapp->grid[clevel].coloring); CHKERRQ(info);
90: PetscStackPop;
92: return(0);
93: }
95: int DAAppSetMultigridKSP(GridCtx *, PetscInt, KSP);
96: int DAAppSetupMultigridMonitor(TAO_APPLICATION, DA, PetscInt, void *);
97: int TaoAppGetMultiGridApplication(TAO_APPLICATION, TaoMultiGridApplication *);
98: static char TaoPetscPCMGDAAppXXX[] = "TaoPetscPCMGDAApp";
99: static char TaoPetscAppXXX[] = "TaoPetscApp";
104: static int DAAppDestroyTaoAppXXX(void*ctx){
105: TaoMultiGridApplication *mctx=(TaoMultiGridApplication*)ctx;
107: if (mctx){ delete mctx;}
108: return(0);
109: }
113: int TaoAppGetMultiGridApplication(TAO_APPLICATION daapplication, TaoMultiGridApplication **mgdaapp){
114: int info;
115: PetscInt clevel,ii;
116: Mat H,HH;
117: MatStructure pflag=SAME_NONZERO_PATTERN;
118: DA_APPLICATION daapp;
119: TaoMultiGridApplication* daapp2;
120: TaoPetscApplication *pppappp;
124: info=TaoAppQueryForObject(daapplication,TaoPetscPCMGDAAppXXX,(void**)&daapp2); CHKERRQ(info);
125: if (daapp2==0){
126: daapp2=new TaoMultiGridApplication(PETSC_COMM_WORLD);
127: daapp2->papp=daapplication;
128: info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
130: info=TaoAppQueryForObject(daapplication,TaoPetscAppXXX,(void**)&pppappp); CHKERRQ(info);
131: if (pppappp && pppappp->taoh && 0==1){
132: info=pppappp->taoh->GetMatrix(&H,&HH,&pflag); CHKERRQ(info);
133: info=daapp2->taoh->SetMatrix(H,HH,pflag);CHKERRQ(info);
134: } else {
135: info=DAAppGetCurrentLevel(daapplication,&clevel); CHKERRQ(info);
136: info=daapp2->taoh->SetMatrix(daapp->grid[clevel].H,daapp->grid[clevel].H,pflag);CHKERRQ(info);
137: }
138: if (pppappp){
139: daapp2->tao=pppappp->tao;
140: }
141: info=TaoAppAddObject(daapplication,TaoPetscPCMGDAAppXXX,(void*)daapp2,&ii); CHKERRQ(info);
142: info=TaoAppQueryRemoveObject(daapplication,TaoPetscAppXXX); CHKERRQ(info);
143: info=TaoAppAddObject(daapplication,TaoPetscAppXXX,(void*)daapp2,&ii); CHKERRQ(info);
144: info=TaoAppSetDestroyRoutine(daapplication,DAAppDestroyTaoAppXXX, (void*)daapp2); CHKERRQ(info);
145: }
146: *mgdaapp=daapp2;
147: return(0);
148: }
150: static int UsePCMGGG=0;
153: /*@
154: DAAppUseMultigrid - Set the preconditioner for the linear solver to be an algebraic multigrid.
156: Collective on TAO_APPLICATION
157:
158: Input Parameters:
159: + daapplication - the DA Application object
160: - coarselevel - the coarsest grid to be used in the multigrid preconditioner. (Grid 0 is the coarsest grid.
162: Level: intermediate
164: Options Database Key:
165: + -tao_da_multigrid - use multigrid linear solver
166: - -ksp_view - view the linear solver
168: Note:
169: This function should be called after DAAppSetHessianRoutine();
171: Note:
172: This function should be called before each optimization solver as part of the DAAppMonitor
174: Note:
175: Multigrid functionality is still under developement for good performance.
177: .seealso: TaoAppGetKSP(), DAAppSetupMultigrid()
179: Options Database Key:
180: . -tao_da_multigrid
182: .keywords: Linear Solver, Multigrid, DA, KSP
184: @*/
185: int DAAppUseMultigrid(TAO_APPLICATION daapplication, PetscInt coarselevel){
186: int info;
188: UsePCMGGG=coarselevel;
189: info=DAAppSetBeforeMonitor(daapplication,DAAppSetupMultigridMonitor,0);
190: CHKERRQ(info);
191: return(0);
192: }
197: int DAAppSetupMultigridMonitor(TAO_APPLICATION daapplication, DA da, PetscInt level, void *ctx){
198: int info;
200: info = DAAppSetupMultigrid(daapplication,UsePCMGGG); CHKERRQ(info);
201: return(0);
202: }
207: /*@
208: DAAppSetupMultigrid - Sets up the multigrid preconditioner.
210: Collective on TAO_APPLICATION
212: Input Parameters:
213: + daapplication - the TAO_APPLICATION context
214: - coarselevel - the coarsest grid that should be used in the multigrid (>=0)
216: Note:
217: Usually the coarselevel is set to 0;
219: Level: intermediate
221: .seealso: DAAppUseMultigrid(), TaoAppSetMonitor()
222: @*/
223: int DAAppSetupMultigrid(TAO_APPLICATION daapplication, PetscInt coarselevel){
224: int info;
225: PetscInt nn,level;
226: TaoMultiGridMatPetsc *MMM;
227: DA_APPLICATION daapp;
228: TaoMultiGridApplication *mgdaapp;
229: TAO_SOLVER tao;
230: KSP ksp;
231:
233: info=DAAppGetCurrentLevel(daapplication,&level); CHKERRQ(info);
234: info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
235: if (daapp->grid[level].mgrid==0 && coarselevel<level){
236: info=TaoAppGetMultiGridApplication(daapplication,&mgdaapp); CHKERRQ(info);
237: mgdaapp->coarselevel=PetscMax(0,coarselevel);
238: nn=level-coarselevel+1;
239: if ( mgdaapp->taoh ){
240: MMM=(TaoMultiGridMatPetsc*)(mgdaapp->taoh);
241: info=MMM->SetUpMultiGrid(daapp->grid+coarselevel,nn); CHKERRQ(info);
242: } else if ( mgdaapp->taoj ){
243: MMM=(TaoMultiGridMatPetsc*)(mgdaapp->taoj);
244: info=MMM->SetUpMultiGrid(daapp->grid+coarselevel,nn); CHKERRQ(info);
245: mgdaapp->taoj=MMM;
246: }
248: info=TaoAppGetTaoSolver(daapplication,&tao); CHKERRQ(info);
249: info=TaoSetupApplicationSolver(daapplication,tao); CHKERRQ(info);
250: info=TaoAppGetKSP(daapplication,&ksp); CHKERRQ(info);
251: if (ksp){
253: info = DAAppSetMultigridKSP(daapp->grid+coarselevel,nn,ksp);CHKERRQ(info);
254: }
255: daapp->grid[level].mgrid=1;
256: }
257: return(0);
258: }
261:
264: int DAAppSetMultigridKSP(GridCtx* dagrid, PetscInt nlevels, KSP ksp){
265: int info;
266: PetscInt i;
267: PC pc,pc2;
268: KSP sksp;
269: const KSPType ksptype;
270: const PCType pctype;
271: PCMGType form;
274: if (ksp==0 || nlevels<=1){
275: return(0);
276: }
279: info = PetscInfo(ksp,"Set multigrid precondition in optimization solver.\n"); CHKERRQ(info);
280:
281: info = KSPGetType(ksp,&ksptype); CHKERRQ(info);
282:
283: info = KSPGetPC(ksp,&pc); CHKERRQ(info);
284: info = PCGetType(pc,&pctype); CHKERRQ(info);
285: info = PCSetType(pc,PCMG); CHKERRQ(info);
286:
287: info = PCMGSetLevels(pc,nlevels,PETSC_NULL); CHKERRQ(info);
289: for (i=0; i<nlevels; i++) {
290:
291: info = PCMGGetSmoother(pc,i,&sksp); CHKERRQ(info);
292: info = KSPSetType(sksp,ksptype); CHKERRQ(info);
294: info = KSPGetPC(sksp,&pc2); CHKERRQ(info);
295: info = PCSetType(pc2,PCJACOBI); CHKERRQ(info);
297: form=PC_MG_MULTIPLICATIVE;
298: info=PCMGSetType(pc,form);CHKERRQ(info);
300: info = KSPSetOperators(sksp,dagrid[i].H,dagrid[i].H,SAME_NONZERO_PATTERN); CHKERRQ(info);
301:
302: if (i<nlevels-1){
303: info = PCMGSetX(pc,i,dagrid[i].X); CHKERRQ(info);
304: info = PCMGSetRhs(pc,i,dagrid[i].RHS); CHKERRQ(info);
305: }
306: if (i>0){
307: info = PCMGSetR(pc,i,dagrid[i].R); CHKERRQ(info);
308: }
309: info = PCMGSetResidual(pc,i,PCMGDefaultResidual,dagrid[i].H); CHKERRQ(info);
310: }
311: for (i=1; i<nlevels; i++) {
312: info = PCMGSetInterpolation(pc,i,dagrid[i].Interpolate); CHKERRQ(info);
313: info = PCMGSetRestriction(pc,i,dagrid[i].Interpolate); CHKERRQ(info);
314: }
315:
316: return(0);
317: }
321: /* @
322: DAAppSetQuadraticObjective - identify the objective function as quadratic or not.
324: Collective on TAO_APPLICATION
326: Input Parameters:
327: + daapplication - the DA Application object
328: - flag - indicates whether the objective is quadratic or not
330: Level: intermediate
332: Options Database Key:
333: . -tao_da_quadratic - use multigrid linear solver
335: Note:
336: The default value is PETSC_FALSE.
338: Note:
339: If quadratic, consider setting the flag used in KSP to SAME_PRECONDITIONER
341: Note:
342: If quadratic, this routine reduces the number of Hessian evaluations done when using
343: the multigrid preconditioner.
345: .seealso: DAAppSetMatStructure(), DAAppUseMultigrid()
348: .keywords: DA, Objective Function
350: @ */
351: int DAAppSetQuadraticObjective(TAO_APPLICATION daapplication, PetscTruth pflag){
352: int info;
353: TaoMultiGridApplication *mgdaapp;
356: info=TaoAppGetMultiGridApplication(daapplication,&mgdaapp); CHKERRQ(info);
357: // mgdaapp->IsQuadratic=pflag;
358: if (pflag==PETSC_TRUE){
359: info = PetscInfo(daapplication,"User labels the objective function as quadratic.\n"); CHKERRQ(info);
360: info = PetscInfo(daapplication,"Set KSP MatStructure to SAME_PRECONDITIONER.\n"); CHKERRQ(info);
361: info=DAAppSetMatStructure(daapplication, SAME_PRECONDITIONER);CHKERRQ(info);
362: } else {
363: info = PetscInfo(daapplication,"User labels the objective function as NOT quadratic.\n"); CHKERRQ(info);
364: info = PetscInfo(daapplication,"Set KSP MatStructure to SAME_NONZERO_PATTERN.\n"); CHKERRQ(info);
365: info=DAAppSetMatStructure(daapplication, SAME_NONZERO_PATTERN);CHKERRQ(info);
366: }
367: return(0);
368: }
373: /* @
374: DAAppSetMultiGridOptions - Sets various multigrid options to be used in this application
375: and the TAO solver.
377: Collective on TAO_APPLICATION
379: Input Parameters:
380: . daapplication - the DA Application object
382: Level: beginner
384: .keywords: options
386: .seealso: TaoDAAppSetOptions();
388: @ */
389: int DAAppSetMultiGridOptions(TAO_APPLICATION daapplication){
390: int info;
391: PetscInt nlevel;
393: PetscTruth flg1=PETSC_FALSE;
396: info = PetscInfo(daapplication,"TaoDAAppSetMultiGridOptions(): Reading command line for options\n"); CHKERRQ(info);
398: flg1=PETSC_FALSE;
399: info = PetscOptionsInt("-tao_da_multigrid","use multigrid linear solver","DAAppUseMultigrid",
400: PETSC_FALSE,&nlevel,&flg1);CHKERRQ(info);
401: if (flg1) {
402: info=DAAppUseMultigrid(daapplication,nlevel);CHKERRQ(info);
403: info = PetscInfo(daapplication,"TaoDAAppSetOptions: Use Multigrid linear solver \n"); CHKERRQ(info);
404: }
405: /*
406: flg1=PETSC_FALSE,flg2=PETSC_FALSE;
407: info = PetscOptionsTruth("-tao_da_quadratic","Identify the objective function as quadratic",
408: "DAAppSetQuadraticObjective",PETSC_FALSE,&flg2,&flg1);CHKERRQ(info);
409: if (flg1) {
410: info = DAAppSetQuadraticObjective(daapplication, flg2);CHKERRQ(info);
411: }
412: */
414: return(0);
415: }