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, int 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 = PetscLogInfo((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=PetscLogClassRegister(&DAAPP_COOKIE,"TAO DA Application");CHKERRQ(info);
 56:   }
 57:   info=PetscObjectGetComm((PetscObject)daapplication,&comm); CHKERRQ(info);
 58:   info = PetscHeaderCreate(daapp,_p_DA_APPLICATION,int,DAAPP_COOKIE,-1,"DA Application",comm,DAAppDestroy,0); CHKERRQ(info);
 59:   info = PetscLogObjectCreate(daapp); CHKERRQ(info);
 60:   info = PetscLogObjectMemory(daapp, sizeof(struct _p_DA_APPLICATION)); CHKERRQ(info);
 61: 
 62:   daapp->nda=0;
 63:   daapp->ndamax=PETSCDAAPPMAXGRIDS;
 64:   daapp->currentlevel=0;
 65:   daapp->IsComplementarity=PETSC_FALSE;
 66:   daapp->kspflag=SAME_NONZERO_PATTERN;
 67:   daapp->computedafunction=0;
 68:   daapp->computedagradient=0;
 69:   daapp->computedafunctiongradient=0;
 70:   daapp->computedahessian=0;
 71:   daapp->usrdafctx=0; daapp->usrdafgctx=0; daapp->usrdagctx=0; daapp->usrdahctx=0;
 72:   daapp->computedabounds=0;
 73:   daapp->bounddactx=0;
 74:   daapp->nbeforemonitors=0;
 75:   daapp->naftermonitors=0;
 76:   int ii;
 77:   info=TaoAppAddObject(daapplication,PPetscDAAppXXX,(void*)daapp,&ii); CHKERRQ(info);
 78:   info=TaoAppSetDestroyRoutine(daapplication,(int(*)(void*))(TaoAppDestroyDAApp), (void*)daapplication); CHKERRQ(info);
 79:   info=TaoAppSetOptionsRoutine(daapplication,DAAppSetOptions); CHKERRQ(info);

 81:   return(0);
 82: }

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

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

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

100:    Level: beginner

102: .keywords: Application, DA
103: @*/
104: int TaoAppSetDAApp(TAO_APPLICATION daapplication, DA* tda, int nnda){
105:   int i,info;
106:   DA_APPLICATION daapp;
107:   PetscScalar zero=0.0;
109:   if (nnda > PETSCDAAPPMAXGRIDS){
110:     SETERRQ1(1,"Number of grids cannot exceed %d.\n",PETSCDAAPPMAXGRIDS); }
111:   info = DAAppExtend(daapplication); CHKERRQ(info);
112:   info = TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
113:   info = PetscLogInfo((daapplication,"TAO Extending TAO_APPLICATION for DA_APPLICATION object.\n")); CHKERRQ(info);
114:   info = DAAppSetMatType(daapplication,MATAIJ); CHKERRQ(info);
115:   daapp->nda=nnda;
116:   daapp->currentlevel=0;
117:   for (i=0;i<nnda; i++){
119:     info = PetscObjectReference((PetscObject)(tda[i])); CHKERRQ(info);
120:     daapp->grid[i].da=tda[i];
121:     info = DACreateGlobalVector(daapp->grid[i].da,&daapp->grid[i].X); CHKERRQ(info);
122:     info = VecDuplicate(daapp->grid[i].X,&daapp->grid[i].R); CHKERRQ(info);
123:     info = VecDuplicate(daapp->grid[i].X,&daapp->grid[i].RHS); CHKERRQ(info);
124:     daapp->grid[i].XL=0;
125:     daapp->grid[i].XU=0;
126:     daapp->grid[i].H=0;
127:     daapp->grid[i].Interpolate=0;
128:     daapp->grid[i].CScale=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 = PetscLogInfo((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:   int i,info,nnda=daapp->nda;


180:   if (--daapp->refct > 0) return(0);

182:   info = PetscLogInfo((daapp,"Destroying work vectors for DA_APPLICATION object.\n")); CHKERRQ(info);
183:   for (i=0;i<nnda; i++){
184:     if (daapp->grid[i].coloring){
185:       info = ISColoringDestroy(daapp->grid[i].coloring); CHKERRQ(info) daapp->grid[i].coloring=0;
186:     }
187:     if (daapp->grid[i].H){
188:       info = MatDestroy(daapp->grid[i].H); CHKERRQ(info);daapp->grid[i].H=0;
189:     }
190:     if (daapp->grid[i].XL && daapp->grid[i].XU){
191:       info = VecDestroy(daapp->grid[i].XL); CHKERRQ(info);daapp->grid[i].XL=0;
192:       info = VecDestroy(daapp->grid[i].XU); CHKERRQ(info);daapp->grid[i].XU=0;
193:     }
194:     if (i>0){
195:       info = MatDestroy(daapp->grid[i].Interpolate); CHKERRQ(info);daapp->grid[i].Interpolate=0;
196:       info = VecDestroy(daapp->grid[i].CScale); CHKERRQ(info);daapp->grid[i].CScale=0;
197:     }
198:     info = VecDestroy(daapp->grid[i].RHS); CHKERRQ(info);daapp->grid[i].RHS=0;
199:     info = VecDestroy(daapp->grid[i].R); CHKERRQ(info);daapp->grid[i].R=0;
200:     info = VecDestroy(daapp->grid[i].X); CHKERRQ(info);daapp->grid[i].X=0;
201:     info = DADestroy(daapp->grid[i].da); CHKERRQ(info);daapp->grid[i].da=0;
202:   }
203:   daapp->nda=0;
204:   daapp->currentlevel=0;
205:   PetscLogObjectDestroy(daapp);
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, int 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, int *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, int *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, clevel;
297:   Mat     H=*M;
298:   DA_APPLICATION daapp=(DA_APPLICATION)ctx;

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

311:   PetscStackPop;

313:   return(0);
314: }


319: /*@
320:   DAAppSetHessianMat - Sets the matrices used to store the Hessian

322:    Collective on TAO_APPLICATION

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

327:    Level: intermediate

329: Note: This routine should be used when finite differences are used to compute
330: the Hessian.

332: .keywords:  DA, hessian

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

336: @*/
337: int DAAppSetHessianMat(TAO_APPLICATION daapplication){
338:   int i,info,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_NO_NEW_NONZERO_LOCATIONS); CHKERRQ(info);
348:       info = MatSetOption(daapp->grid[i].H,MAT_STRUCTURALLY_SYMMETRIC); CHKERRQ(info);
349:       info = MatSetOption(daapp->grid[i].H,MAT_SYMMETRIC); CHKERRQ(info);
350:       info = PetscLogInfo((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_LOCAL,&(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 = PetscLogInfo((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 = PetscLogInfo((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 = PetscLogInfo((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 = PetscLogInfo((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 = PetscLogInfo((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 = PetscLogInfo((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:   int level=daapp->currentlevel, info;

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


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

638:    Collective on TAO_APPLICATION

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

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

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

654:    Level: beginner

656: .keywords:  DA, bounds

658: .seealso: TaoDAAppSolve();

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



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

681:    Collective on TAO_APPLICATION

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

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

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

697:    Level: intermediate

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

702: .keywords:  DA, gradient

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

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


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

733:    Collective on TAO_APPLICATION

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

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

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

749:    Level: intermediate

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

754: .keywords:  DA, hessian

756: .seealso: DAAppSetConstraintRoutine();

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

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

789:    Collective on TAO_APPLICATION

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

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

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

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

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

811:    Level: intermediate

813: .keywords:  DA, monitor

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

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

838: /*@C
839:   DAAppSetAfterMonitor - Set a routine that will be called after the
840:   optimization on each grid.

842:    Collective on TAO_APPLICATION

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

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

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

857: Note:
858:    These monitors are different from the monitors that can be called after
859:    each iteration of the optimization algorithm.

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

864:    Level: intermediate

866: .keywords:  DA, monitor

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

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



893: /*@
894:   DAAppGetSolution - Sets a pointer to an existing vector that contains
895:   the solution.

897:    Collective on TAO_APPLICATION

899:    Input Parameters:
900: +  daapplication - the DA Application object
901: -  level - the DA on which the current solution is desired (Level 0 is coarsest)

903:    Output Parameters:
904: .  X - the current solution

906:    Level: beginner

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

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


928: /*@C
929:   DAAppGetHessianMat - Sets a pointer to an existing matrix that contains
930:   the Hessian at the specified level.

932:    Collective on TAO_APPLICATION

934:    Input Parameters:
935: +  daapplication - the DA Application object
936: -  level - the DA on which the current Hessian is desired (Level 0 is coarsest)

938:    Output Parameters:
939: .  H - the current Hessian

941:    Level: intermediate

943: .seealso: DAAppSetHessianRoutine;

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

962: /*@
963:   DAAppGetInterpolationMatrix - Sets pointers to existing vectors
964:   containg upper and lower bounds on the variables.

966:    Collective on TAO_APPLICATION

968:    Input Parameters:
969: +  daapplication - the DA Application object
970: -  level - the DA on which the current solution is desired (Level 0 is coarsest)

972:    Output Parameters:
973: +  Interpolate - the interpolating matrix
974: -  CScale - the scaling matrix for restriction Matrix

976:    Level: intermediate

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

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


1002: /*@
1003:   DAAppGetVariableBounds - Sets pointers to existing vectors
1004:   containg upper and lower bounds on the variables.

1006:    Collective on TAO_APPLICATION

1008:    Input Parameters:
1009: +  daapplication - the DA Application object
1010: -  level - the DA on which the current solution is desired (Level 0 is coarsest)

1012:    Output Parameters:
1013: +  XL - the current solution
1014: -  XU - the current solution

1016:    Level: intermediate

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

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



1043: /*@
1044:   DAAppSetInitialSolution - Sets the initial solution for
1045:   the grid application.

1047:    Collective on TAO_APPLICATION

1049:    Input Parameters:
1050: +  daapplication - the DA Application object
1051: -  X0 - the initial solution vector

1053:    Level: beginner

1055: .seealso: DAAppGetSolution();

1057: .keywords: Initial Solution, DA
1058: @*/
1059: int DAAppSetInitialSolution(TAO_APPLICATION daapplication, Vec X0){
1060:   int info;
1061:   DA_APPLICATION daapp;
1064:   info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
1065:   info=VecCopy(X0,daapp->grid[0].X); CHKERRQ(info);
1066:   return(0);
1067: }



1073: /*
1074:   DAAppSetMatStructure - Set the MatStructure flag to be used in KSPSetOperators:

1076:   Collective on TAO_APPLICATION
1077:   
1078:   Input Parameters:
1079: +  daapplication - the DA Application object
1080: -  flag - indicates the KSP flag

1082:    Level: intermediate

1084:    Note:
1085:    The default value is SAME_NONZERO_PATTERN

1087: .seealso: KSPSetOperators(), TaoGetKSP();

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

1091: @*/
1092: int DAAppSetMatStructure(TAO_APPLICATION daapplication, MatStructure pflag){
1093:   int info;
1094:   DA_APPLICATION daapp;
1096:   info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
1097:   daapp->kspflag=pflag;
1098:   info = PetscLogInfo((daapp,"Set MatStructure flag in TAO solver.\n")); CHKERRQ(info);
1099:   return(0);
1100: }

1104: /*@
1105:   DAAppSetMatType - set the matrix type to be used for the Hessian matrix

1107:   Collective on TAO_APPLICATION
1108:   
1109:   Input Parameters:
1110: +  daapplication - the DA Application object
1111: -  mattype - the matrix type

1113:    Level: intermediate

1115:    Options Database Key:
1116: .  -tao_da_mattype - select matrix type

1118:    Note:
1119:    The default value is MATMPIAIJ.

1121:    Note:
1122:    Call before DAAppSetHessianRoutine()

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

1126: .keywords: DA, Hessian

1128: @*/
1129: int DAAppSetMatType(TAO_APPLICATION daapplication, MatType mattype){
1130:   int info;
1131:   DA_APPLICATION daapp;
1133:   info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
1134:   info=PetscStrcpy(daapp->HessianMatrixType,mattype); CHKERRQ(info);
1135:   return(0);
1136: }


1141: /*@
1142:   DAAppSetOptions - Sets various parameters to be used in this application
1143:   and the TAO solver.

1145:    Collective on TAO_APPLICATION

1147:    Input Parameters:
1148: .  daapplication - the DA Application object

1150:    Level: intermediate

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

1154: .keywords:  DA, options

1156: .seealso: TaoDAAppSolve();

1158: @*/
1159: int DAAppSetOptions(TAO_APPLICATION daapplication){
1160:   int info;
1161:   //  char *prefix;
1162:   PetscTruth flg1=PETSC_FALSE,flg2=PETSC_FALSE;
1163:   DA_APPLICATION daapp;
1164:   //  MPI_Comm comm;

1168:   info = PetscLogInfo((daapplication,"TaoDAAppSetOptions(): Reading command line for options\n")); CHKERRQ(info);
1169:   info = TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
1170:   /*
1171:   info = PetscObjectGetOptionsPrefix((PetscObject)daapplication,&prefix); CHKERRQ(info);
1172:   info = PetscObjectGetComm((PetscObject)daapplication,&comm); CHKERRQ(info);
1173:   info = PetscOptionsBegin(comm,prefix,"TAO PETSC DA APPLICATIONS","solver");CHKERRQ(info);
1174:   */
1175:   info = DAAppSetMultiGridOptions(daapplication);CHKERRQ(info);

1177:   flg1=PETSC_FALSE,flg2=PETSC_FALSE;
1178:   info = PetscOptionsTruth("-tao_da_samepreconditioner","Use same preconditioner in KSP","DAAppSetMatStructure",
1179:                              PETSC_FALSE,&flg2,&flg1);CHKERRQ(info);
1180:   if (flg1 && flg2) {
1181:     info = DAAppSetMatStructure(daapplication, SAME_PRECONDITIONER);CHKERRQ(info);
1182:   } else if (flg1){
1183:     info = DAAppSetMatStructure(daapplication, SAME_NONZERO_PATTERN);CHKERRQ(info);
1184:   }

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

1189:   /*
1190:   info = PetscOptionsEnd();CHKERRQ(info);
1191:   */
1192:   return(0);
1193: }