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: }