Actual source code: fdiff.c

      1: #include "taoapp.h"         /*I  "tao.h"  I*/
      2: #include "petscsnes.h"

     10: static int Ftemp(SNES snes ,Vec X,Vec G,void*ctx){
     11:   int info;
     12:   TAO_APPLICATION taoapp = (TAO_APPLICATION)ctx;
     14:   info=TaoAppComputeGradient(taoapp,X,G); CHKERRQ(info);
     15:   return(0);
     16: }

     23: /*@C
     24:    TaoAppDefaultComputeHessian - Computes the Hessian using finite differences.

     26:    Collective on TAO_APPLICATION

     28:    Input Parameters:
     29: +  taoapp - the TAO_APPLICATION context
     30: .  V - compute Hessian at this point
     31: -  ctx - the TAO_APPLICATION structure, cast to (void*)

     33:    Output Parameters:
     34: +  H - Hessian matrix (not altered in this routine)
     35: .  B - newly computed Hessian matrix to use with preconditioner (generally the same as H)
     36: -  flag - flag indicating whether the matrix sparsity structure has changed

     38:    Options Database Key:
     39: +  -tao_fd - Activates TaoAppDefaultComputeHessian()
     40: -  -tao_view_hessian - view the hessian after each evaluation using PETSC_VIEWER_STDOUT_WORLD

     42:    Level: intermediate

     44:    Notes:
     45:    This routine is slow and expensive, and is not currently optimized
     46:    to take advantage of sparsity in the problem.  Although
     47:    TaoAppDefaultComputeHessian() is not recommended for general use
     48:    in large-scale applications, It can be useful in checking the
     49:    correctness of a user-provided Hessian.

     51:    Note:
     52:    The gradient evaluation must be set using the routine TaoAppSetGradientRoutine().

     54: .keywords: TAO_APPLICATION, finite differences, Hessian

     56: .seealso: TaoAppSetHessianRoutine(), TaoAppDefaultComputeHessianColor(), SNESDefaultComputeJacobian(),
     57:           TaoAppSetGradientRoutine()

     59: @*/
     60: int TaoAppDefaultComputeHessian(TAO_APPLICATION taoapp,Vec V,Mat *H,Mat *B,
     61:                              MatStructure *flag,void *ctx){
     62:   int                  info;
     63:   MPI_Comm             comm;
     64:   Vec                  G;
     65:   SNES                 snes;

     70:   info = VecDuplicate(V,&G);CHKERRQ(info);

     72:   info = PetscLogInfo((G,"TAO Using finite differences w/o coloring to compute matrix\n")); CHKERRQ(info);

     74:   info = TaoAppComputeGradient(taoapp,V,G); CHKERRQ(info);

     76:   info = PetscObjectGetComm((PetscObject)(*H),&comm);CHKERRQ(info);
     77:   info = SNESCreate(comm,&snes);CHKERRQ(info);

     79:   info = SNESSetFunction(snes,G,Ftemp,ctx);CHKERRQ(info);
     80:   info = SNESSetJacobian(snes,*H,*B,SNESDefaultComputeJacobian,ctx);CHKERRQ(info);
     81:   //  info = SNESComputeJacobian(snes,V,H,B,flag);CHKERRQ(info);
     82:   info = SNESDefaultComputeJacobian(snes,V,H,B,flag,ctx);CHKERRQ(info);

     84:   info = SNESDestroy(snes);CHKERRQ(info);
     85:
     86:   info = VecDestroy(G);CHKERRQ(info);
     87:
     88:   return(0);
     89: }

     98: /*@C
     99:    TaoAppDefaultComputeHessianColor - Computes the Hessian using colored finite differences.

    101:    Collective on TAO_APPLICATION

    103:    Input Parameters:
    104: +  tao - the TAO_APPLICATION context
    105: .  V - compute Hessian at this point
    106: -  ctx - the TAO_APPLICATION structure, cast to (void*)

    108:    Output Parameters:
    109: +  H - Hessian matrix (not altered in this routine)
    110: .  B - newly computed Hessian matrix to use with preconditioner (generally the same as H)
    111: -  flag - flag indicating whether the matrix sparsity structure has changed

    113:    Options Database Keys:
    114: +  -mat_fd_coloring_freq <freq>
    115: -  -tao_view_hessian - view the hessian after each evaluation using PETSC_VIEWER_STDOUT_WORLD

    117:    Level: intermediate

    119:    Note:
    120:    The gradient evaluation must be set using the routine TaoSetPetscGradient().

    122:  .keywords: TAO_APPLICATION, finite differences, Hessian, coloring, sparse

    124: .seealso: TaoAppSetHessianRoutine(), TaoAppDefaultComputeHessian(),SNESDefaultComputeJacobianColor(),
    125:           TaoAppSetGradientRoutine(), TaoAppSetColoring()

    127: @*/
    128: int TaoAppDefaultComputeHessianColor(TAO_APPLICATION taoapp, Vec V, Mat *HH,Mat *BB,
    129:                                   MatStructure *flag,void *ctx){
    130:   int                 info;
    131:   MPI_Comm            comm;
    132:   Vec                 G=0;
    133:   Mat                 H=*HH,B=*BB;
    134:   SNES                snes;
    135:   ISColoring          iscoloring;
    136:   MatFDColoring       matfdcoloring;

    144:   info = TaoAppGetColoring(taoapp,&iscoloring); CHKERRQ(info);
    145:   if (!iscoloring){
    146:     SETERRQ(1,"Must set coloring before using this routine.  Try Finite Differences without coloring\n");
    147:   }
    148:   info = VecDuplicate(V,&G);CHKERRQ(info);

    150:   info=PetscLogInfo((G,"TAO computing matrix using finite differences and coloring\n")); CHKERRQ(info);

    152:   info=TaoAppComputeGradient(taoapp,V,G); CHKERRQ(info);

    154:   info = PetscObjectGetComm((PetscObject)(H),&comm);CHKERRQ(info);
    155:   info = SNESCreate(comm,&snes);CHKERRQ(info);

    157:   info = MatFDColoringCreate(H,iscoloring,&matfdcoloring);CHKERRQ(info);
    158:   info = MatFDColoringSetFunction(matfdcoloring,(int (*)(void)) Ftemp,ctx);CHKERRQ(info);
    159:   info = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(info);

    161:   info = SNESSetFunction(snes,G,Ftemp,ctx);CHKERRQ(info);
    162:   info = SNESSetJacobian(snes,H,B,SNESDefaultComputeJacobianColor,(void*)matfdcoloring);CHKERRQ(info);
    163:   info = SNESDefaultComputeJacobianColor(snes,V,HH,BB,flag,matfdcoloring);CHKERRQ(info);

    165:   info = MatFDColoringDestroy(matfdcoloring);CHKERRQ(info);
    166:   info = SNESDestroy(snes);CHKERRQ(info);
    167:
    168:   info = VecDestroy(G);CHKERRQ(info);
    169:   return(0);
    170: }

    175: /* @
    176:   TaoAppSetFiniteDifferencesOptions - Sets various TAO parameters from user options

    178:    Collective on TAO_APPLICATION

    180:    Input Parameters:
    181: +  taoapp - the TAO Application (optional)

    183:    Level: beginner

    185: .keywords:  options, finite differences

    187: .seealso: TaoSolveApplication();

    189: @ */
    190: int TaoAppSetFiniteDifferencesOptions(TAO_APPLICATION taoapp){
    191:   int info;
    192:   PetscTruth flg;

    198:   flg=PETSC_FALSE;
    199:   info = PetscOptionsName("-tao_fd","use finite differences for Hessian","TaoAppDefaultComputeHessian",&flg);CHKERRQ(info);
    200:   if (flg) {
    201:     info = TaoAppSetHessianRoutine(taoapp,TaoAppDefaultComputeHessian,(void*)taoapp);CHKERRQ(info);
    202:     info = PetscLogInfo((taoapp,"TaoSetOptions: Setting default finite difference Hessian matrix\n")); CHKERRQ(info);
    203:   }

    205:   flg=PETSC_FALSE;
    206:   info = PetscOptionsName("-tao_fd_coloring","use finite differences with coloring to compute Hessian","TaoAppDefaultComputeHessianColor",&flg);CHKERRQ(info);
    207:   if (flg) {
    208:     info = TaoAppSetHessianRoutine(taoapp,TaoAppDefaultComputeHessianColor,(void*)taoapp);CHKERRQ(info);
    209:     info = PetscLogInfo((taoapp,"TaoAppSetFromOptions: Use finite differencing with coloring to compute Hessian \n")); CHKERRQ(info);
    210:   }
    211:
    212:   return(0);
    213: }

    216: static char TaoAppColoringXXX[] = "TaoColoring";

    218: typedef struct {
    219:   ISColoring coloring;
    220: } TaoAppColorit;

    224: static int TaoAppDestroyColoringXXX(void*ctx){
    225:   //  int info;
    226:   TaoAppColorit *mctx=(TaoAppColorit*)ctx;
    228:   if (mctx){
    229:     /*
    230:     if (mctx->coloring){
    231:       info = ISColoringDestroy(mctx->coloring);CHKERRQ(info);
    232:     }
    233:     */
    234:     PetscFree(mctx);
    235:   }
    236:   return(0);
    237: }

    241: /*@
    242:    TaoAppSetColoring - Set the matrix coloring to be used when computing the
    243:    Hessian by finite differences.

    245:    Collective on TAO_APPLICATION

    247:    Input Parameters:
    248: +  app - the TAO_APPLICATION context
    249: -  coloring - the coloring to be used

    251:    Level: intermediate

    253:    Note:
    254:    The gradient evaluation must be set using the routine TaoSetPetscGradient().

    256:  .keywords: TAO_APPLICATION, finite differences, Hessian, coloring, sparse

    258: .seealso: TaoAppSetHessianRoutine(), TaoAppDefaultComputeHessianColor(),
    259:           TaoAppSetGradientRoutine()

    261: @*/
    262: int TaoAppSetColoring(TAO_APPLICATION taoapp, ISColoring coloring){
    263:   int ii,info;
    264:   TaoAppColorit *mctx;
    267:
    268:   info = TaoAppQueryForObject(taoapp,TaoAppColoringXXX,(void**)&mctx); CHKERRQ(info);
    269:   if (mctx==0){
    270:     info=PetscNew(TaoAppColorit,(void**)&mctx); CHKERRQ(info);
    271:     info=TaoAppAddObject(taoapp,TaoAppColoringXXX,(void*)mctx,&ii); CHKERRQ(info);
    272:     info=TaoAppSetDestroyRoutine(taoapp,TaoAppDestroyColoringXXX, (void*)mctx); CHKERRQ(info);
    273:   }
    274:   /*
    275:   if (coloring){
    276:   info=PetscObjectReference((PetscObject)coloring); CHKERRQ(info);
    277:   }
    278:   if (mctx->coloring){
    279:        info=ISColoringDestroy(mctx->coloring); CHKERRQ(info);
    280:   }
    281:   */
    282:   mctx->coloring=coloring;
    283:   return(0);
    284: }

    288: /*@C
    289:    TaoAppGetColoring - Set the matrix coloring to be used when computing the
    290:    Hessian by finite differences.

    292:    Collective on TAO_APPLICATION

    294:    Input Parameters:
    295: +  app - the TAO_APPLICATION context
    296: -  coloring - the coloring to be used

    298:    Level: advanced

    300:    Note:
    301:    The gradient evaluation must be set using the routine TaoAppSetGradientRoutine().

    303:  .keywords: TAO_APPLICATION, finite differences, Hessian, coloring, sparse

    305: .seealso: TaoAppSetHessianRoutine(), TaoAppDefaultComputeHessianColor(),
    306:           TaoAppSetGradientRoutine()

    308: @*/
    309: int TaoAppGetColoring(TAO_APPLICATION taoapp, ISColoring *coloring){
    310:   int info;
    311:   TaoAppColorit *mctx;
    314:   if (coloring){
    315:     info = TaoAppQueryForObject(taoapp,TaoAppColoringXXX,(void**)&mctx); CHKERRQ(info);
    316:     if (mctx){
    317:       *coloring=mctx->coloring;
    318:     } else {
    319:       *coloring=0;
    320:     }
    321:   }
    322:   return(0);
    323: }

    328: int TaoAppAddFiniteDifferences(TAO_APPLICATION taoapp){
    329:   int info;
    331:   info = TaoAppSetOptionsRoutine(taoapp,TaoAppSetFiniteDifferencesOptions); CHKERRQ(info);
    332:   return(0);
    333: }

    335: int MatTAOMFSetGradient(Mat mat,Vec v,int (*func)(TAO_APPLICATION,Vec,Vec,void *),void *funcctx){
    337:   return(0);
    338: }