Actual source code: daapp.c

  1: #include "taodaapplication.h"     /*I "taodaapplication.h" I*/
  2: #include "taoapp.h"
  3: #include "petsc.h"
  4: #include "daapp_impl.h"

  6: int DAAPP_COOKIE=0;
  7: static char PPetscDAAppXXX[] = "PPetscDAApp";

  9: int TaoAppDAApp(TAO_APPLICATION, DA_APPLICATION *);
 10: int DAAppDestroy(DA_APPLICATION);
 11: int DAAppExtend(TAO_APPLICATION);

 15: /* @C
 16:   DAApplicationCreate - Creates a TaoApplication defined
 17: a PETSc Distributed Array (DA) object.  The routines used to evaluate
 18: the objective function, constraint, and derivative information are
 19: designed specifically for these problems.

 21:    Input Parameters:
 22: +  comm - an MPI communicator
 23: .  tda - an array of Distributed Array objects
 24: -  nnda - the number of Distibuted Array objects

 26:    Output Parameters:
 27: .  newdaapplication - the TaoApplication structure

 29: .seealso TaoDAAppSolve(), TaoApplicationCreate(), TaoAppDestroy();

 31:    Level: beginner

 33: .keywords: Application, DA
 34: @ */
 35: int DAApplicationCreate(MPI_Comm comm, DA* tda, PetscInt nnda, TAO_APPLICATION* newdaapplication){
 36:   int info;
 37:   TAO_APPLICATION ttapp;
 39:   info = TaoApplicationCreate(comm,&ttapp); CHKERRQ(info);
 40:   info = TaoAppSetDAApp(ttapp,tda,nnda); CHKERRQ(info);
 41:   info = PetscInfo(tda[0],"Creating a DA_APPLICATION object.\n"); CHKERRQ(info);
 42:   *newdaapplication=ttapp;
 43:   return(0);
 44: }

 48: int DAAppExtend(TAO_APPLICATION daapplication){
 49:   int info;
 50:   MPI_Comm comm;
 51:   DA_APPLICATION daapp;

 54:   if (DAAPP_COOKIE==0){
 55:       info=PetscCookieRegister("TAO DA Application",&DAAPP_COOKIE);CHKERRQ(info);
 56:   }
 57:   info=PetscObjectGetComm((PetscObject)daapplication,&comm); CHKERRQ(info);
 58:   info = PetscHeaderCreate(daapp,_p_DA_APPLICATION,PetscInt,DAAPP_COOKIE,-1,"DA Application",comm,DAAppDestroy,0); CHKERRQ(info);
 59:   
 60:   daapp->nda=0;
 61:   daapp->ndamax=PETSCDAAPPMAXGRIDS;
 62:   daapp->currentlevel=0;
 63:   daapp->IsComplementarity=PETSC_FALSE;
 64:   daapp->kspflag=SAME_NONZERO_PATTERN;
 65:   daapp->computedafunction=0;
 66:   daapp->computedagradient=0;
 67:   daapp->computedafunctiongradient=0;
 68:   daapp->computedahessian=0;
 69:   daapp->usrdafctx=0; daapp->usrdafgctx=0; daapp->usrdagctx=0; daapp->usrdahctx=0;
 70:   daapp->computedabounds=0;
 71:   daapp->bounddactx=0;
 72:   daapp->nbeforemonitors=0;
 73:   daapp->naftermonitors=0;
 74:   PetscInt ii;
 75:   info=TaoAppAddObject(daapplication,PPetscDAAppXXX,(void*)daapp,&ii); CHKERRQ(info);
 76:   info=TaoAppSetDestroyRoutine(daapplication,(int(*)(void*))(TaoAppDestroyDAApp), (void*)daapplication); CHKERRQ(info);
 77:   info=TaoAppSetOptionsRoutine(daapplication,DAAppSetOptions); CHKERRQ(info);

 79:   return(0);
 80: }

 84: /*@
 85:   TaoAppSetDAApp - Extend the functionality of a TaoApplication by 
 86: adding support for applications defined
 87: a PETSc Distributed Array (DA) object.  The routines used to evaluate
 88: the objective function, constraint, and derivative information are
 89: designed specifically for these problems.

 91:    Input Parameters:
 92: +  daappliation - an existing application object
 93: .  tda - an array of Distributed Array objects
 94: -  nnda - the number of Distibuted Array objects

 96: .seealso TaoDAAppSolve(), TaoApplicationCreate();

 98:    Level: beginner

100: .keywords: Application, DA
101: @*/
102: int TaoAppSetDAApp(TAO_APPLICATION daapplication, DA* tda, PetscInt nnda){
103:   int info;
104:   PetscInt i;
105:   DA_APPLICATION daapp;
106:   PetscScalar zero=0.0;
108:   if (nnda > PETSCDAAPPMAXGRIDS){
109:     SETERRQ1(1,"Number of grids cannot exceed %d.\n",PETSCDAAPPMAXGRIDS); }
110:   info = DAAppExtend(daapplication); CHKERRQ(info);
111:   info = TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
112:   info = PetscInfo(daapplication,"TAO Extending TAO_APPLICATION for DA_APPLICATION object.\n"); CHKERRQ(info);
113:   info = DAAppSetMatType(daapplication,(const MatType)MATAIJ); CHKERRQ(info);
114:   daapp->nda=nnda;
115:   daapp->currentlevel=0;
116:   for (i=0;i<nnda; i++){
118:     info = PetscObjectReference((PetscObject)(tda[i])); CHKERRQ(info);
119:     daapp->grid[i].da=tda[i];
120:     info = DACreateGlobalVector(daapp->grid[i].da,&daapp->grid[i].X); CHKERRQ(info);
121:     info = VecDuplicate(daapp->grid[i].X,&daapp->grid[i].R); CHKERRQ(info);
122:     info = VecDuplicate(daapp->grid[i].X,&daapp->grid[i].RHS); CHKERRQ(info);
123:     daapp->grid[i].XL=0;
124:     daapp->grid[i].XU=0;
125:     daapp->grid[i].H=0;
126:     daapp->grid[i].Interpolate=0;
127:     daapp->grid[i].CScale=0;
128:     daapp->grid[i].mgrid=0;
129:     if (i>0){
130:       info = DAGetInterpolation(daapp->grid[i-1].da,daapp->grid[i].da,&daapp->grid[i].Interpolate,&daapp->grid[i].CScale); CHKERRQ(info);
131:     }
132:   }
133:   info = VecSet(daapp->grid[0].X, zero); CHKERRQ(info);
134:   info = TaoAppSetInitialSolutionVec(daapplication,daapp->grid[0].X); CHKERRQ(info);
135:   info = PetscInfo(daapp,"Create work vectors for DA_APPLICATION object.\n"); CHKERRQ(info);

137:   return(0);
138: }

142: int TaoAppDAApp(TAO_APPLICATION daapplication, DA_APPLICATION *daapp){
143:   int info;
144:   DA_APPLICATION daapp2;
147:   info=TaoAppQueryForObject(daapplication,PPetscDAAppXXX,(void**)&daapp2); CHKERRQ(info);
148:   if (daapp2){
150:   } else {
151:     SETERRQ(1,"TAO ERROR: Must call TaoAppSetDAApp() first");
152:   }
153:   *daapp=daapp2;
154:   return(0);
155: }

159: int TaoAppDestroyDAApp(TAO_APPLICATION daapplication){
160:   int info;
161:   DA_APPLICATION daapp;
163:   info=TaoAppQueryForObject(daapplication,PPetscDAAppXXX,(void**)&daapp); CHKERRQ(info);
164:   if (daapp){
165:     info=DAAppDestroy(daapp); CHKERRQ(info);
166:     //    info=TaoAppQueryRemoveObject(daapplication,PPetscDAAppXXX); CHKERRQ(info);
167:   }

169:   return(0);
170: }

174: int DAAppDestroy(DA_APPLICATION daapp){
175:   PetscInt i,nnda=daapp->nda;
176:   int info;


181:   if (--((PetscObject)daapp)->refct > 0) return(0);

183:   info = PetscInfo(daapp,"Destroying work vectors for DA_APPLICATION object.\n"); CHKERRQ(info);
184:   for (i=0;i<nnda; i++){
185:     if (daapp->grid[i].coloring){
186:       info = ISColoringDestroy(daapp->grid[i].coloring); CHKERRQ(info); daapp->grid[i].coloring=0;
187:     }
188:     if (daapp->grid[i].H){
189:       info = MatDestroy(daapp->grid[i].H); CHKERRQ(info);daapp->grid[i].H=0;
190:     }
191:     if (daapp->grid[i].XL && daapp->grid[i].XU){
192:       info = VecDestroy(daapp->grid[i].XL); CHKERRQ(info);daapp->grid[i].XL=0;
193:       info = VecDestroy(daapp->grid[i].XU); CHKERRQ(info);daapp->grid[i].XU=0;
194:     }
195:     if (i>0){
196:       info = MatDestroy(daapp->grid[i].Interpolate); CHKERRQ(info);daapp->grid[i].Interpolate=0;
197:       info = VecDestroy(daapp->grid[i].CScale); CHKERRQ(info);daapp->grid[i].CScale=0;
198:     }
199:     info = VecDestroy(daapp->grid[i].RHS); CHKERRQ(info);daapp->grid[i].RHS=0;
200:     info = VecDestroy(daapp->grid[i].R); CHKERRQ(info);daapp->grid[i].R=0;
201:     info = VecDestroy(daapp->grid[i].X); CHKERRQ(info);daapp->grid[i].X=0;
202:     info = DADestroy(daapp->grid[i].da); CHKERRQ(info);daapp->grid[i].da=0;
203:   }
204:   daapp->nda=0;
205:   daapp->currentlevel=0;
206:   PetscHeaderDestroy(daapp); 
207:   return(0);
208: }

212: /*@
213:   DAAppGetDA - Get the DA on a the specified level.
214:   

216:    Input Parameters:
217: +  daapplication - the TAO DAApplication structure
218: -  n - the level

220:    Output Parameter:
221: .  da - address of the pointer to the DA.


224: .seealso TaoAppSetDAApp();

226:    Level: intermediate

228: .keywords: Application, DA
229: @*/
230: int DAAppGetDA(TAO_APPLICATION daapplication, PetscInt n, DA *da){
231:   int info;
232:   DA_APPLICATION daapp;
234:   info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
235:   *da=daapp->grid[n].da;
236:   return(0);
237: }

241: /*@
242:   DAAppGetNumberOfDAGrids - Get the number of grids specified on the application
243:   

245:    Input Parameters:
246: .  daapplication - the TAO DAApplication structure

248:    Output Parameter:
249: .  n - number of DA grids specified

251: .seealso TaoAppSetDAApp();

253:    Level: intermediate

255: .keywords: Application, DA
256: @*/
257: int DAAppGetNumberOfDAGrids(TAO_APPLICATION daapplication, PetscInt *n){
258:   int info;
259:   DA_APPLICATION daapp;
261:   info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
262:   *n=daapp->nda;
263:   return(0);
264: }

268: /*@
269:   DAAppGetCurrentLevel - Get the number of the current DA grid

271:    Input Parameters:
272: .  daapplication - the TAO DAApplication structure

274:    Output Parameter:
275: .  n - number of DA grids specified

277: .seealso TaoAppSetDAApp();

279:    Level: intermediate

281: .keywords: Application, DA
282: @*/
283: int DAAppGetCurrentLevel(TAO_APPLICATION daapplication, PetscInt *n){
284:   int info;
285:   DA_APPLICATION daapp;
287:   info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
288:   *n=daapp->currentlevel;
289:   return(0);
290: }


295: static int DAAppPetscHessian(TAO_APPLICATION daapplication , Vec X , Mat *M, Mat *MPre, MatStructure *flag, void*ctx){
296:   int     info;
297:   PetscInt clevel;
298:   Mat     H=*M;
299:   DA_APPLICATION daapp=(DA_APPLICATION)ctx;

304:   PetscStackPush("TAO User DA Hessian");
305:   clevel=daapp->currentlevel;
306:   
307:   //    daapp->FDGrad=daapp->grid[i].RHS;   /* Needed for finite differences gradient vector */
308:   info = PetscInfo1(daapp,"TAO hessian evaluation at DA_APPLICATION object, level %d.\n",clevel); CHKERRQ(info);
309:   info = (*daapp->computedahessian)(daapplication,daapp->grid[clevel].da,X,H,daapp->usrdahctx);CHKERRQ(info);
310:   *flag=daapp->kspflag;

312:   PetscStackPop;

314:   return(0);
315: }


320: /*@
321:   DAAppSetHessianMat - Creates the matrices used to store the Hessian at finer grid levels

323:    Collective on TAO_APPLICATION

325:    Input Parameters:
326: .  daapplication - the DA Application object

328:    Level: advanced


331: .keywords:  DA, hessian

333: .seealso: DAAppSetHessianRoutine(), DAAppSetMatType();

335: @*/
336: int DAAppSetHessianMat(TAO_APPLICATION daapplication){
337:   int info;
338:   PetscInt i,nnda;
339:   DA_APPLICATION daapp;

342:   info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
343:   nnda=daapp->nda;
344:   for (i=0;i<nnda; i++){
345:     if (daapp->grid[i].H==0){
346:       info = DAGetMatrix(daapp->grid[i].da,daapp->HessianMatrixType,&daapp->grid[i].H); CHKERRQ(info);
347:       info = MatSetOption(daapp->grid[i].H,MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE); CHKERRQ(info);
348:       info = MatSetOption(daapp->grid[i].H,MAT_STRUCTURALLY_SYMMETRIC, PETSC_TRUE); CHKERRQ(info);
349:       info = MatSetOption(daapp->grid[i].H,MAT_SYMMETRIC, PETSC_TRUE); CHKERRQ(info);
350:       info = PetscInfo1(daapp->grid[i].da,"TAO create Hessian for DA_APPLICATION object, level %d.\n",i); CHKERRQ(info);
351:       info = DAGetColoring(daapp->grid[i].da,IS_COLORING_GLOBAL,daapp->HessianMatrixType,&(daapp->grid[i].coloring));CHKERRQ(info);
352:     }
353:   }
354:   info = TaoAppSetHessianMat(daapplication,daapp->grid[0].H,daapp->grid[0].H);CHKERRQ(info);
355:   
356:   return(0);
357: }
360: /*@C
361:   DAAppSetHessianRoutine - Set a routine that will evaluate the hessian
362:   on the given DA at the given point.

364:    Collective on TAO_APPLICATION

366:    Input Parameters:
367: +  daapplication - the DA Application object
368: .  hess - the function pointer for the hessian evaluation routine
369: -  ctx - the hessian context

371:    Calling sequence of hess:
372: $     hess(TAO_APPLICATION daapplication,DA da, Vec x,Mat H,void *ctx);

374: +  daapplication - the TAO_APPLICATION context
375: .  da - the Distributed Array
376: .  x - input vector
377: .  H - hessian matrix
378: -  ctx - user-defined hessian context set from DAAppSetHessianRoutine()

380:    Level: beginner

382:    Options Database Key:
383: .  -tao_view_hessian - view the hessian after each evaluation using PETSC_VIEWER_STDOUT_WORLD

385: .keywords:  DA, hessian

387: .seealso: DAAppSetObjectiveAndGradientRoutine();
388: @*/
389: int DAAppSetHessianRoutine(TAO_APPLICATION daapplication, int (*hess)(TAO_APPLICATION,DA,Vec,Mat,void*),void *ctx){
390:   int info;
391:   DA_APPLICATION daapp;

394:   info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
395:   daapp->computedahessian=hess;
396:   daapp->usrdahctx=ctx;
397:   if (hess){
398:     info = TaoAppSetHessianRoutine(daapplication,DAAppPetscHessian,(void*)daapp);CHKERRQ(info);
399:   } else {
400:     info = TaoAppSetHessianRoutine(daapplication,0,0);CHKERRQ(info);
401:   }
402:   info = DAAppSetHessianMat(daapplication);CHKERRQ(info);
403:   return(0);
404: }

408: static int DAAppPetscFunctionGradient(TAO_APPLICATION daapplication , Vec X ,double *ff, Vec G, void*ctx){
409:   int info;
410:   DA_APPLICATION daapp=(DA_APPLICATION)ctx;

415:   PetscStackPush("TAO User DA Objective and gradient function");
416:   info = (*daapp->computedafunctiongradient)(daapplication,daapp->grid[daapp->currentlevel].da,X,ff,G,daapp->usrdafgctx);
417:   CHKERRQ(info);
418:   info = PetscInfo1(daapplication,"TAO function evaluation at DA_APPLICATION object, level %d.\n",daapp->currentlevel); CHKERRQ(info);
419:   PetscStackPop;

421:   return(0);
422: }


427: /*@C
428:   DAAppSetObjectiveAndGradientRoutine - Set a routine that will evaluate the objective and
429:   gradient functions on the given DA at the given point.

431:    Collective on TAO_APPLICATION

433:    Input Parameters:
434: +  daapplication - the DA Application object
435: .  grad - the function pointer for the gradient evaluation routine
436: -  ctx - the function-gradient context

438:    Calling sequence of funcgrad:
439: $     funcgrad(TAO_APPLICATION daapplication,DA da, Vec x,double *f,Vec g,void *ctx);

441: +  daapplication - the TAO_APPLICATION daapplication context
442: .  da - the Distributed Array
443: .  x - input vector
444: .  f - objective value 
445: .  g - gradient vector 
446: -  ctx - user defined function-gradient context set from DAAppSetObjectiveAndGradientRoutine()

448:    Fortran Note:
449:    If your Fortran compiler does not recognize symbols over 31 characters in length, then
450:    use the identical routine with the shortened name DAAppSetObjectiveAndGradientRou()

452:    Level: beginner

454:    Options Database Key:
455: .  -tao_view_gradient - view the gradient after each evaluation using PETSC_VIEWER_STDOUT_SELF

457: .keywords: DA, Gradient, Objective Function

459: .seealso: DAAppSetObjectiveRoutine(), DAAppSetGradientRoutine();

461: @*/
462: int DAAppSetObjectiveAndGradientRoutine(TAO_APPLICATION daapplication, int (*funcgrad)(TAO_APPLICATION,DA,Vec,double*,Vec, void*),void *ctx){
463:   int info;
464:   DA_APPLICATION daapp;
466:   info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
467:   if (funcgrad){ 
468:     info=TaoAppSetObjectiveAndGradientRoutine(daapplication,
469:                                               DAAppPetscFunctionGradient,(void*)daapp);CHKERRQ(info);
470:   } else {
471:     info=TaoAppSetObjectiveAndGradientRoutine(daapplication,0,0);CHKERRQ(info);
472:   }
473:   info=TaoAppSetInitialSolutionVec(daapplication,daapp->grid[0].X);CHKERRQ(info);
474:   daapp->computedafunctiongradient=funcgrad;
475:   daapp->usrdafgctx=ctx;
476:   info = PetscInfo(daapp,"Set objective function pointer for DA_APPLICATION object.\n"); CHKERRQ(info);
477:   return(0);
478: }

482: static int DAAppPetscGradient(TAO_APPLICATION daapplication , Vec X ,Vec G, void*ctx){
483:   int info;
484:   DA_APPLICATION daapp=(DA_APPLICATION)ctx;

489:   PetscStackPush("TAO User DA Gradient Evaluation");
490:   info = (*daapp->computedagradient)(daapplication,daapp->grid[daapp->currentlevel].da,X,G,daapp->usrdafgctx);
491:   CHKERRQ(info);
492:   info = PetscInfo1(daapplication,"TAO gradient evaluation at DA_APPLICATION object, level %d.\n",daapp->currentlevel); CHKERRQ(info);
493:   PetscStackPop;

495:   return(0);
496: }


501: /*@C
502:   DAAppSetGradientRoutine - Set a routine that will evaluate the gradient 
503:   function on the given DA at the given point.

505:    Collective on TAO_APPLICATION

507:    Input Parameters:
508: +  daapplication - the DA Application object
509: .  grad - the function pointer for the gradient evaluation routine
510: -  ctx - the gradient context

512:    Calling sequence of grad:
513: $     grad(TAO_APPLICATION daapplication,DA da, Vec x,Vec g,void *ctx);

515: +  daapplication - the TAO_APPLICATION context
516: .  da - the Distributed Array
517: .  x - input vector
518: .  g - gradient vector 
519: -  ctx - user defined gradient context set from DAAppSetGradientRoutine()

521:    Level: intermediate

523:    Options Database Key:
524: .  -tao_view_gradient - view the gradient after each evaluation using PETSC_VIEWER_STDOUT_SELF

526: .keywords:  DA, gradient

528: .seealso: DAAppSetObjectiveRoutine(), DAAppSetObjectiveAndGradientRoutine();

530: @*/
531: int DAAppSetGradientRoutine(TAO_APPLICATION daapplication, int (*grad)(TAO_APPLICATION,DA,Vec,Vec, void*),void *ctx){
532:   int info;
533:   DA_APPLICATION daapp;
535:   info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
536:   if (grad){ 
537:     info=TaoAppSetGradientRoutine(daapplication, DAAppPetscGradient, (void*)daapp);CHKERRQ(info);
538:   } else {
539:     info=TaoAppSetGradientRoutine(daapplication,0,0);CHKERRQ(info);
540:   }
541:   daapp->computedagradient=grad;
542:   daapp->usrdagctx=ctx;
543:   info = PetscInfo(daapp,"Set gradient function pointer for DA_APPLICATION object.\n"); CHKERRQ(info);
544:   return(0);
545: }

549: static int DAAppPetscObjectiveFunction(TAO_APPLICATION daapplication , Vec X ,double* f, void*ctx){
550:   int info;
551:   DA_APPLICATION daapp=(DA_APPLICATION)ctx;

556:   PetscStackPush("TAO User DA Objective function");
557:   info = (*daapp->computedafunction)(daapplication,daapp->grid[daapp->currentlevel].da,X,f,daapp->usrdafgctx);
558:   CHKERRQ(info);
559:   info = PetscInfo1(daapp,"TAO gradient evaluation at DA_APPLICATION object, level %d.\n",daapp->currentlevel); CHKERRQ(info);
560:   PetscStackPop;

562:   return(0);
563: }


568: /*@C
569:   DAAppSetObjectiveRoutine - Set a routine that will evaluate the objective 
570:   function on the given DA at the given point.

572:    Collective on TAO_APPLICATION

574:    Input Parameters:
575: +  daapplication - the DA Application object
576: .  func - the function pointer for the objecive evaluation routine
577: -  ctx - the monitor context

579:    Calling sequence of func:
580: $     func(TAO_APPLICATION daapplication,DA da,Vec x,double *f,void *ctx);

582: +  daapplication - the TAO_APPLICATION context
583: .  da - the Distributed Array
584: .  x - input vector
585: .  f - application sets equal to the function value 
586: -  ctx - user-defined function context set from DAAppSetObjectiveRoutine()

588:    Level: beginner

590: .keywords:  DA, objective

592: .seealso: DAAppSetObjectiveAndGradientRoutine();

594: @*/
595: int DAAppSetObjectiveRoutine(TAO_APPLICATION daapplication, int (*func)(TAO_APPLICATION,DA,Vec,double*, void*),void *ctx){
596:   int info;
597:   DA_APPLICATION daapp;
599:   info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
600:   if (func){
601:     info=TaoAppSetObjectiveRoutine(daapplication, DAAppPetscObjectiveFunction, (void*)daapp);CHKERRQ(info);
602:   } else {
603:     info=TaoAppSetObjectiveRoutine(daapplication,0,0);CHKERRQ(info);
604:   }
605:   info=TaoAppSetInitialSolutionVec(daapplication, daapp->grid[0].X);CHKERRQ(info);
606:   daapp->computedafunction=func;
607:   daapp->usrdafctx=ctx;
608:   info = PetscInfo(daapp,"Set objective function pointer for DA_APPLICATION object.\n"); CHKERRQ(info);
609:   return(0);
610: }

614: int DAApplicationEvaluateVariableBounds(TAO_APPLICATION daapplication, Vec XL, Vec XU, void*ctx){
615:   DA_APPLICATION daapp=(DA_APPLICATION)ctx;
616:   PetscInt level=daapp->currentlevel;
617:   int  info;

621:   info = PetscInfo1(daapp->grid[level].da,"Compute variable bounds in level %d of DA_APPLICATION object.\n",level); CHKERRQ(info);
622:   info = PetscObjectReference((PetscObject)XL);
623:   info = PetscObjectReference((PetscObject)XU);
624:   daapp->grid[level].XL=XL;
625:   daapp->grid[level].XU=XU;
626:   if (daapp->computedabounds){
627:     info = (*daapp->computedabounds)(daapplication,daapp->grid[level].da,XL,XU,daapp->bounddactx); CHKERRQ(info);
628:   }
629:   return(0);
630: }


635: /*@C
636:   DAAppSetVariableBoundsRoutine - Set a routine that will evaluate the bounds
637:   on the variables.

639:    Collective on TAO_APPLICATION

641:    Input Parameters:
642: +  daapplication - the DA Application object
643: .  bound - the function pointer for the bound evaluation routine
644: -  ctx - the monitor context

646:    Calling sequence of bound:
647: $     bound(TAO_APPLICATION daapplication, DA da,Vec xl, Vec xu, void *ctx);

649: +  daapplication - the TAO_APPLICATION context
650: .  da - the Distributed Array
651: .  xl - vector of upper bounds
652: .  xu - vector of lower bounds
653: -  ctx - user-defined monitor context set from DAAppSetVariableBoundsRoutine()

655:    Level: beginner

657: .keywords:  DA, bounds

659: .seealso: TaoDAAppSolve();

661: @*/
662: int DAAppSetVariableBoundsRoutine(TAO_APPLICATION daapplication, int (*bound)(TAO_APPLICATION,DA,Vec,Vec, void*),void *ctx){
663:   int info;
664:   DA_APPLICATION daapp;
666:   info = TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
667:   info=TaoAppSetVariableBoundsRoutine(daapplication,DAApplicationEvaluateVariableBounds,(void*)daapp);CHKERRQ(info);
668:   daapp->computedabounds=bound;
669:   daapp->bounddactx=ctx;
670:   info = PetscInfo(daapp,"Set variable bounds in DA_APPLICATION object.\n"); CHKERRQ(info);
671:   return(0);
672: }



678: /*@C
679:   DAAppSetConstraintRoutine - Set a routine that will evaluate the constraint
680:   function on the given DA at the given point.

682:    Collective on TAO_APPLICATION

684:    Input Parameters:
685: +  daapplication - the TAO Application object
686: .  grad - the function pointer for the gradient evaluation routine
687: -  ctx - the gradient context

689:    Calling sequence of grad:
690: $     f(TAO_APPLICATION daapplication,DA da, Vec x,Vec r,void *ctx);

692: +  daapplication - the DA_APPLICATION context
693: .  da - the Distributed Array
694: .  x - input vector
695: .  r - constraint vector 
696: -  ctx - user defined gradient context set from DAAppSetGradientRoutine()

698:    Level: intermediate

700:    Options Database Key:
701: .  -tao_view - view the gradient after each evaluation using PETSC_VIEWER_STDOUT_SELF

703: .keywords:  DA, gradient

705: .seealso: DAAppSetJacobianRoutine();
706: @*/
707: int DAAppSetConstraintRoutine(TAO_APPLICATION daapplication, int (*f)(TAO_APPLICATION,DA,Vec,Vec, void*),void *ctx){
708:   int info;
709:   DA_APPLICATION daapp;
711:   info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
712:   if (f){ 
713:     info=TaoAppSetConstraintRoutine(daapplication, DAAppPetscGradient, (void*)daapp);CHKERRQ(info);
714:   } else {
715:     info=TaoAppSetConstraintRoutine(daapplication,0,0);CHKERRQ(info);
716:   }
717:   daapp->IsComplementarity=PETSC_TRUE;
718:   info=DAAppSetGradientRoutine(daapplication,f,ctx); CHKERRQ(info);
719:   daapp->computedagradient=f;
720:   daapp->usrdagctx=ctx;
721:   daapp->usrdafgctx=ctx;

723:   info = PetscInfo(daapp,"Set constraint function pointer for DA_APPLICATION object.\n"); CHKERRQ(info);
724:   return(0);
725: }


730: /*@C
731:   DAAppSetJacobianRoutine - Set a routine that will evaluate the Jacobian
732:   on the given DA at the given point.

734:    Collective on TAO_APPLICATION

736:    Input Parameters:
737: +  daapplication - the DA Application object
738: .  jac - the function pointer for the Jacobian evaluation routine
739: -  ctx - the jacobian context

741:    Calling sequence of hess:
742: $     jac(TAO_APPLICATION daapplication, DA da, Vec x,Mat J,void *ctx);

744: +  daapplication - the TAO_APPLICATION context
745: .  da - the Distributed Array
746: .  x - input vector
747: .  J - Jacobian matrix
748: -  ctx - user-defined hessian context set from DAAppSetHessianRoutine()

750:    Level: intermediate

752:    Options Database Key:
753: .  -tao_view_jacobian - view the jacobian after each evaluation using PETSC_VIEWER_STDOUT_WORLD

755: .keywords:  DA, hessian

757: .seealso: DAAppSetConstraintRoutine();

759: @*/
760: int DAAppSetJacobianRoutine(TAO_APPLICATION daapplication, int (*jac)(TAO_APPLICATION,DA,Vec,Mat,void*),void *ctx){
761:   int info;
762:   PetscInt i,nnda;
763:   DA_APPLICATION daapp;
765:   info = TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
766:   if (!jac){ SETERRQ(1,"No DA Hessian routine provided\n");}
767:   daapp->computedahessian=jac;
768:   daapp->usrdahctx=ctx;
769:   nnda=daapp->nda;;
770:   for (i=0;i<nnda; i++){
771:     if (daapp->grid[i].H==0){
772:       info = DAGetMatrix(daapp->grid[i].da,daapp->HessianMatrixType,&daapp->grid[i].H); CHKERRQ(info);
773:       info = MatSetOption(daapp->grid[i].H,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE); CHKERRQ(info);
774:       info = MatSetOption(daapp->grid[i].H,MAT_STRUCTURALLY_SYMMETRIC,PETSC_TRUE); CHKERRQ(info);
775:       info = PetscInfo1(daapp->grid[i].da,"TAO create Jacobian for DA_APPLICATION object, level %d.\n",i); CHKERRQ(info);
776:     }
777:   }
778:   
779:   info = TaoAppSetJacobianRoutine(daapplication,DAAppPetscHessian,(void*)daapp);CHKERRQ(info);
780:   info = TaoAppSetJacobianMat(daapplication,daapp->grid[0].H,daapp->grid[0].H);CHKERRQ(info);
781:   
782:   return(0);
783: }

787: /*@C
788:   DAAppSetBeforeMonitor - Set a routine that will be called before the
789:   optimization on each grid.

791:    Collective on TAO_APPLICATION

793:    Input Parameters:
794: +  daapplication - the DA Application object
795: .  beforemonitor - a monitor routine called before solving the application on each DA 
796: -  ctx - the monitor context

798:    Calling sequence of monitor:
799: $     beforemonitor(TAO_APPLICATION daapplication, DA da,int level, void *ctx);

801: +  daapplication - this TAO_APPLICATION context
802: .  da - the Distributed Array
803: .  level - the grid level that will be solved next (level 0 is coarsest)
804: -  ctx - user-defined function context set from DAAppSetBeforeMonitor()

806: Note:
807:    These monitors are different from the monitors that can be called after
808:    each iteration of the optimization algorithm.

810: Note:
811:    The beforemonitor and aftermonitor are good for setting up and destroying the application data.

813:    Level: intermediate

815: .keywords:  DA, monitor

817: .seealso: DAAppSetAfterMonitor(), TaoSetMonitor(), TaoAppSetMonitor();

819: @*/
820: int DAAppSetBeforeMonitor(TAO_APPLICATION daapplication, int (*beforemonitor)(TAO_APPLICATION,DA,PetscInt, void*), void *ctx){
821:   PetscInt n;
822:   int info;
823:   DA_APPLICATION daapp;
825:   info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
826:   if (beforemonitor){
827:     n=daapp->nbeforemonitors;
828:     if (n>=MAX_DAAP_MONITORS){
829:       SETERRQ(1,"TAO ERROR: Too many beforemonitors set in DAAPP");
830:     }
831:     daapp->beforemonitor[n]=beforemonitor;
832:     daapp->beforemonitorctx[n]=ctx;
833:     daapp->nbeforemonitors++;
834:     info = PetscInfo(daapplication,"TAO: Set user beforemonitor for DA_APPLICATION object.\n"); CHKERRQ(info);
835:   }
836:   return(0);
837: }

841: /*@C
842:   DAAppSetAfterMonitor - Set a routine that will be called after the
843:   optimization on each grid.

845:    Collective on TAO_APPLICATION

847:    Input Parameters:
848: +  daapplication - the DA Application object
849: .  aftermonitor - a monitor routine called after solving the application on each DA 
850: -  ctx - the monitor context

852:    Calling sequence of monitor:
853: $     aftermonitor(TAO_APPLICATION daapplication, DA da, int level, void *ctx);

855: +  daapplication - this TAO_APPLICATION context
856: .  da - the Distributed Array
857: .  level - the grid level that will be solved next (level 0 is coarsest)
858: -  ctx - user-defined function context set from DAAppSetAfterMonitor()

860: Note:
861:    These monitors are different from the monitors that can be called after
862:    each iteration of the optimization algorithm.

864: Note:
865:    The beforemonitor and aftermonitor are good for setting up and destroying the application data.

867:    Level: intermediate

869: .keywords:  DA, monitor

871: .seealso: DAAppSetBeforeMonitor(), TaoSetMonitor(), TaoAppSetMonitor();

873: @*/
874: int DAAppSetAfterMonitor(TAO_APPLICATION daapplication, int (*aftermonitor)(TAO_APPLICATION,DA,PetscInt, void*), void *ctx){
875:   int info;
876:   PetscInt n;
877:   DA_APPLICATION daapp;
879:   info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
880:   if (aftermonitor){
881:     n=daapp->naftermonitors;
882:     if (n>=MAX_DAAP_MONITORS){
883:       SETERRQ(1,"TAO ERROR: Too many aftermonitors set in DAAPP");
884:     }
885:     daapp->aftermonitor[n]=aftermonitor;
886:     daapp->aftermonitorctx[n]=ctx;
887:     daapp->naftermonitors++;
888:     info = PetscInfo(daapp,"TAO: Set user aftermonitor for DA_APPLICATION object.\n"); CHKERRQ(info);
889:   }
890:   return(0);
891: }



897: /*@
898:   DAAppGetSolution - Sets a pointer to an existing vector that contains
899:   the solution.

901:    Collective on TAO_APPLICATION

903:    Input Parameters:
904: +  daapplication - the DA Application object
905: -  level - the DA on which the current solution is desired (Level 0 is coarsest)

907:    Output Parameters:
908: .  X - the current solution

910:    Level: beginner

912: .seealso: DAAppSetBeforeMonitor(), DAAppSetAfterMonitor();

914: .keywords: Solution, DA
915: @*/
916: int DAAppGetSolution(TAO_APPLICATION daapplication, PetscInt level, Vec *X){
917:   int info;
918:   DA_APPLICATION daapp;
920:   info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
921:   if (level<0 || level>=daapp->nda){
922:     SETERRQ(1,"No Solution available.  This grid does not exist");
923:   } else if (X && daapp && daapp->grid && level<daapp->nda){
924:     *X=daapp->grid[level].X;
925:   }
926:   return(0);
927: }


932: /*@C
933:   DAAppGetHessianMat - Sets a pointer to an existing matrix that contains
934:   the Hessian at the specified level.

936:    Collective on TAO_APPLICATION

938:    Input Parameters:
939: +  daapplication - the DA Application object
940: -  level - the DA on which the current Hessian is desired (Level 0 is coarsest)

942:    Output Parameters:
943: .  H - the current Hessian

945:    Level: intermediate

947: .seealso: DAAppSetHessianRoutine;

949: .keywords: Hessian, DA
950: @*/
951: int DAAppGetHessianMat(TAO_APPLICATION daapplication, PetscInt level, Mat *H){
952:   int info;
953:   DA_APPLICATION daapp;
955:   info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
956:   if (level<0 || level>=daapp->nda){
957:     SETERRQ(1,"No Hessian available.  This grid does not exist");
958:   } else if (H && daapp && daapp->grid && level<daapp->nda){
959:     *H=daapp->grid[level].H;
960:   }
961:   return(0);
962: }

966: /*@
967:   DAAppGetInterpolationMatrix - Sets pointers to existing vectors
968:   containg upper and lower bounds on the variables.

970:    Collective on TAO_APPLICATION

972:    Input Parameters:
973: +  daapplication - the DA Application object
974: -  level - the DA on which the current solution is desired (Level 0 is coarsest)

976:    Output Parameters:
977: +  Interpolate - the interpolating matrix
978: -  CScale - the scaling matrix for restriction Matrix

980:    Level: intermediate

982: .seealso: DAAppSetBeforeMonitor(), DAAppSetAfterMonitor();

984: .keywords: Interpolate, DA
985: @*/
986: int DAAppGetInterpolationMatrix(TAO_APPLICATION daapplication, PetscInt level, Mat *Interpolate, Vec *CScale){
987:   int info;
988:   DA_APPLICATION daapp;
990:   info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
991:   if ( !daapp || level<1 || level>=daapp->nda){
992:     SETERRQ1(1,"No Interpolator available for level %d.  This grid does not exist",level);
993:   } 
994:   if (Interpolate){
995:     *Interpolate=daapp->grid[level].Interpolate;
996:   }
997:   if (CScale){
998:     *CScale=daapp->grid[level].CScale;
999:   }
1000:   return(0);
1001: }


1006: /*@
1007:   DAAppGetVariableBounds - Sets pointers to existing vectors
1008:   containg upper and lower bounds on the variables.

1010:    Collective on TAO_APPLICATION

1012:    Input Parameters:
1013: +  daapplication - the DA Application object
1014: -  level - the DA on which the current solution is desired (Level 0 is coarsest)

1016:    Output Parameters:
1017: +  XL - the current solution
1018: -  XU - the current solution

1020:    Level: intermediate

1022: .seealso: DAAppSetBeforeMonitor(), DAAppSetAfterMonitor();

1024: .keywords: Solution, DA
1025: @*/
1026: int DAAppGetVariableBounds(TAO_APPLICATION daapplication, PetscInt level, Vec *XL, Vec *XU){
1027:   int info;
1028:   DA_APPLICATION daapp;
1030:   info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
1031:   if ( !daapp || level<0 || level>=daapp->nda){
1032:     SETERRQ(1,"No Solution available.  This grid does not exist.");
1033:   } 
1034:   if (XL){
1035:     *XL=daapp->grid[level].XL;
1036:   }
1037:   if (XU){
1038:     *XU=daapp->grid[level].XU;
1039:   }
1040:   return(0);
1041: }



1047: /*@
1048:   DAAppSetInitialSolution - Sets the initial solution for
1049:   the grid application.

1051:    Collective on TAO_APPLICATION

1053:    Input Parameters:
1054: +  daapplication - the DA Application object
1055: -  X0 - the initial solution vector

1057:    Level: beginner

1059: .seealso: DAAppGetSolution();

1061: .keywords: Initial Solution, DA
1062: @*/
1063: int DAAppSetInitialSolution(TAO_APPLICATION daapplication, Vec X0){
1064:   int info;
1065:   DA_APPLICATION daapp;
1068:   info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
1069:   info=VecCopy(X0,daapp->grid[0].X); CHKERRQ(info);
1070:   return(0);
1071: }



1077: /*
1078:   DAAppSetMatStructure - Set the MatStructure flag to be used in KSPSetOperators:

1080:   Collective on TAO_APPLICATION
1081:   
1082:   Input Parameters:
1083: +  daapplication - the DA Application object
1084: -  flag - indicates the KSP flag

1086:    Level: intermediate

1088:    Note:
1089:    The default value is SAME_NONZERO_PATTERN

1091: .seealso: KSPSetOperators(), TaoAppGetKSP();

1093: .keywords: Linear Solver, Multigrid, DA, KSP

1095: @*/
1096: int DAAppSetMatStructure(TAO_APPLICATION daapplication, MatStructure pflag){
1097:   int info;
1098:   DA_APPLICATION daapp;

1101:   info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
1102:   daapp->kspflag=pflag;
1103:   info = PetscInfo(daapp,"Set MatStructure flag in TAO solver.\n"); CHKERRQ(info);

1105:   return(0);
1106: }

1110: /*@
1111:   DAAppSetMatType - set the matrix type to be used for the Hessian matrix

1113:   Collective on TAO_APPLICATION
1114:   
1115:   Input Parameters:
1116: +  daapplication - the DA Application object
1117: -  mattype - the matrix type

1119:    Level: intermediate

1121:    Options Database Key:
1122: .  -tao_da_mattype - select matrix type

1124:    Note:
1125:    The default value is MATMPIAIJ.

1127:    Note:
1128:    Call before DAAppSetHessianRoutine()

1130: .seealso: DAAppSetMatStructure(),  DAAppSetHessianRoutine(), DAGetMatrix();

1132: .keywords: DA, Hessian

1134: @*/
1135: int DAAppSetMatType(TAO_APPLICATION daapplication, const MatType mattype){
1136:   int info;
1137:   DA_APPLICATION daapp;

1140:   info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
1141:   info=PetscStrcpy(daapp->HessianMatrixType,mattype); CHKERRQ(info);

1143:   return(0);
1144: }


1149: /*@
1150:   DAAppSetOptions - Sets various parameters to be used in this application
1151:   and the TAO solver.

1153:    Collective on TAO_APPLICATION

1155:    Input Parameters:
1156: .  daapplication - the DA Application object

1158:    Level: intermediate

1160: Note:  This routine will be called by the TaoAppSetFromOptions().

1162: .keywords:  DA, options

1164: .seealso: TaoDAAppSolve();

1166: @*/
1167: int DAAppSetOptions(TAO_APPLICATION daapplication){
1168:   int info;
1169:   //  char *prefix;
1170:   PetscTruth flg1=PETSC_FALSE,flg2=PETSC_FALSE;
1171:   DA_APPLICATION daapp;
1172:   //  MPI_Comm comm;

1176:   info = PetscInfo(daapplication,"TaoDAAppSetOptions(): Reading command line for options\n"); CHKERRQ(info);
1177:   info = TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
1178:   /*
1179:   info = PetscObjectGetOptionsPrefix((PetscObject)daapplication,&prefix); CHKERRQ(info);
1180:   info = PetscObjectGetComm((PetscObject)daapplication,&comm); CHKERRQ(info);
1181:   info = PetscOptionsBegin(comm,prefix,"TAO PETSC DA APPLICATIONS","solver");CHKERRQ(info);
1182:   */
1183:   info = DAAppSetMultiGridOptions(daapplication);CHKERRQ(info);


1186:   flg1=PETSC_FALSE,flg2=PETSC_FALSE;
1187:   info = PetscOptionsTruth("-tao_da_samepreconditioner","Use same preconditioner in KSP","DAAppSetMatStructure",
1188:                              PETSC_FALSE,&flg2,&flg1);CHKERRQ(info);
1189:   if (flg1 && flg2) {
1190:     info = DAAppSetMatStructure(daapplication, SAME_PRECONDITIONER);CHKERRQ(info); 
1191:   } else if (flg1){
1192:     info = DAAppSetMatStructure(daapplication, SAME_NONZERO_PATTERN);CHKERRQ(info); 
1193:   }

1195:   info = PetscOptionsString("-tao_da_mattype","the hessian should have matrix type","DAAppSetMatType",
1196:                             daapp->HessianMatrixType,daapp->HessianMatrixType,18,&flg1);CHKERRQ(info);

1198:   /*
1199:   info = PetscOptionsEnd();CHKERRQ(info);
1200:   */
1201:   return(0);
1202: }