Actual source code: fdiff.c

  1: #include "taoapp.h"         /*I  "taoapp.h"  I*/
  2: #include "src/tao_impl.h"      /*I "src/tao_impl.h"  I*/
  3: #include "src/petsctao/include/taopetsc.h"          /*I "src/petsctao/include/taopetsc.h" I*/
  4: #include "petscsnes.h"



 12: static int Ftemp(SNES snes ,Vec X,Vec G,void*ctx){
 13:   int info;
 14:   TAO_APPLICATION taoapp = (TAO_APPLICATION)ctx;
 15:   TAO_SOLVER tao;
 18:   info = TaoAppGetTaoSolver(taoapp,&tao);
 19:   info=TaoAppComputeGradient(taoapp,X,G); CHKERRQ(info);
 20:   tao->ngrads++;
 21:   return(0);
 22: }

 27: /*@C
 28:   TaoAppDefaultComputeHessian - computes the gradient using finite differences.
 29:  
 30:   Collective on TAO_APPLICATION

 32:   Input Parameters:
 33: + taoapp - the TAO_APPLICATION context
 34: . X - compute gradient at this point
 35: - ctx - the TAO_APPLICATION structure, cast to (void*)

 37:   Output Parameters:
 38: . G - Gradient Vector

 40:    Options Database Key:
 41: +  -tao_fd_gradient - Activates TaoAppDefaultComputeGradient()
 42: -  -tao_fd_delta <delta> - change in x used to calculate finite differences




 47:    Level: intermediate

 49:    Note:
 50:    This routine is slow and expensive, and is not currently optimized
 51:    to take advantage of sparsity in the problem.  Although
 52:    TaoAppDefaultComputeGradient is not recommended for general use
 53:    in large-scale applications, It can be useful in checking the
 54:    correctness of a user-provided gradient.  Use the tao method "tao_fd_test"
 55:    to get an indication of whether your gradient is correct.


 58:    Note:
 59:    The gradient evaluation must be set using the routine TaoAppSetGradientRoutine().

 61: .keywords: TAO_APPLICATION, finite differences, Hessian

 63: .seealso: TaoAppDefaultComputeGradient(),  TaoAppSetGradientRoutine()

 65: @*/
 66: int TaoAppDefaultComputeGradient(TAO_APPLICATION taoapp,Vec X,Vec G,void*) 
 67: {
 68:   Vec TempX;
 69:   double *g;
 70:   double f, f2;
 71:   int info;
 72:   PetscInt low,high,N,i;
 73:   PetscTruth flg;
 74:   double h=1.0e-6;
 75:   TAO_SOLVER tao;
 77:   info = TaoAppGetTaoSolver(taoapp, &tao);
 78:   //PetscStackPush("TAO Finite Difference gradient");
 79:   //info = PetscLogEventBegin(Tao_FiniteDifferenceGradient,taoapp,X,0,0);
 80:   info = TaoAppComputeObjective(taoapp, X,&f); CHKERRQ(info);
 81:   tao->nfuncs++;
 82:   info = PetscOptionsGetReal(PETSC_NULL,"-tao_fd_delta",&h,&flg);
 83:   info = VecDuplicate(X,&TempX); CHKERRQ(info);
 84:   info = VecCopy(X,TempX); CHKERRQ(info);
 85:   info = VecGetSize(X,&N); CHKERRQ(info);
 86:   info = VecGetOwnershipRange(TempX,&low,&high); CHKERRQ(info);
 87:   info = VecGetArray(G,&g); CHKERRQ(info);
 88:   for (i=0;i<N;i++) {
 89:     info = VecSetValue(TempX,i,h,ADD_VALUES);
 90:     info = TaoAppComputeObjective(taoapp,TempX,&f2); CHKERRQ(info);
 91:     tao->nfuncs++;
 92:     info = VecSetValue(TempX,i,-h,ADD_VALUES);
 93:     if (i>=low && i<high) {
 94:       g[i-low]=(f2-f)/h;
 95:     }

 97:   }
 98:   info = VecRestoreArray(G,&g); CHKERRQ(info);
 99:   
100:   //info = PetscLogEventEnd(Tao_FiniteDifferenceGradient, taoapp, X, 0,0); CHKERRQ(info);
101:   //PetscStackPop;
102:   
103:   return(0);
104: }

108: /*@C
109:    TaoAppDefaultComputeHessian - Computes the Hessian using finite differences. 

111:    Collective on TAO_APPLICATION

113:    Input Parameters:
114: +  taoapp - the TAO_APPLICATION context 
115: .  V - compute Hessian at this point
116: -  ctx - the TAO_APPLICATION structure, cast to (void*)

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

123:    Options Database Key:
124: +  -tao_fd - Activates TaoAppDefaultComputeHessian()
125: -  -tao_view_hessian - view the hessian after each evaluation using PETSC_VIEWER_STDOUT_WORLD

127:    Level: intermediate

129:    Notes:
130:    This routine is slow and expensive, and is not currently optimized
131:    to take advantage of sparsity in the problem.  Although
132:    TaoAppDefaultComputeHessian() is not recommended for general use
133:    in large-scale applications, It can be useful in checking the
134:    correctness of a user-provided Hessian.

136:    Note:
137:    The gradient evaluation must be set using the routine TaoAppSetGradientRoutine().

139: .keywords: TAO_APPLICATION, finite differences, Hessian

141: .seealso: TaoAppSetHessianRoutine(), TaoAppDefaultComputeHessianColor(), SNESDefaultComputeJacobian(),
142:           TaoAppSetGradientRoutine(), TaoAppDefaultComputeGradient()

144: @*/
145: int TaoAppDefaultComputeHessian(TAO_APPLICATION taoapp,Vec V,Mat *H,Mat *B,
146:                              MatStructure *flag,void *ctx){
147:   int                  info;
148:   MPI_Comm             comm;
149:   Vec                  G;
150:   SNES                 snes;
151:   TAO_SOLVER tao;


155:   info = TaoAppGetTaoSolver(taoapp, &tao);

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

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

162:   info = TaoAppComputeGradient(taoapp,V,G); CHKERRQ(info);
163:   tao->ngrads++;

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

168:   info = SNESSetFunction(snes,G,Ftemp,taoapp);CHKERRQ(info);
169:   info = SNESDefaultComputeJacobian(snes,V,H,B,flag,taoapp);CHKERRQ(info);

171:   info = SNESDestroy(snes);CHKERRQ(info);
172:   
173:   info = VecDestroy(G);CHKERRQ(info);
174:   
175:   return(0);
176: }




185: /*@C
186:    TaoAppDefaultComputeHessianColor - Computes the Hessian using colored finite differences. 

188:    Collective on TAO_APPLICATION

190:    Input Parameters:
191: +  tao - the TAO_APPLICATION context
192: .  V - compute Hessian at this point
193: -  ctx - the TAO_APPLICATION structure, cast to (void*)

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

200:    Options Database Keys:
201: +  -mat_fd_coloring_freq <freq>
202: -  -tao_view_hessian - view the hessian after each evaluation using PETSC_VIEWER_STDOUT_WORLD

204:    Level: intermediate

206:    Note:
207:    The gradient evaluation must be set using the routine TaoSetPetscGradient().

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

211: .seealso: TaoAppSetHessianRoutine(), TaoAppDefaultComputeHessian(),SNESDefaultComputeJacobianColor(), 
212:           TaoAppSetGradientRoutine(), TaoAppSetColoring()

214: @*/
215: int TaoAppDefaultComputeHessianColor(TAO_APPLICATION taoapp, Vec V, Mat *HH,Mat *BB,
216:                                   MatStructure *flag,void *ctx){
217:   int                 info;
218:   MPI_Comm            comm;
219:   Vec                 G=0;
220:   Mat                 H=*HH,B=*BB;
221:   SNES                snes;
222:   ISColoring          iscoloring;
223:   MatFDColoring       matfdcoloring;
224:   TAO_SOLVER          tao;


232:   info = TaoAppGetTaoSolver(taoapp,&tao); CHKERRQ(info);
233:   info = TaoAppGetColoring(taoapp,&iscoloring); CHKERRQ(info);
234:   if (!iscoloring){
235:     SETERRQ(1,"Must set coloring before using this routine.  Try Finite Differences without coloring\n");
236:   }
237:   info = VecDuplicate(V,&G);CHKERRQ(info);

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

241:   info=TaoAppComputeGradient(taoapp,V,G); CHKERRQ(info);
242:   tao->ngrads++;

244:   info = PetscObjectGetComm((PetscObject)(H),&comm);CHKERRQ(info);
245:   info = SNESCreate(comm,&snes);CHKERRQ(info);

247:   info = MatFDColoringCreate(H,iscoloring,&matfdcoloring);CHKERRQ(info);
248:   info = MatFDColoringSetFunction(matfdcoloring,(int (*)(void)) Ftemp,taoapp);CHKERRQ(info);
249:   info = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(info);

251:   info = SNESSetFunction(snes,G,Ftemp,taoapp);CHKERRQ(info);
252:   info = SNESSetJacobian(snes,H,B,SNESDefaultComputeJacobianColor,(void*)matfdcoloring);CHKERRQ(info);
253:   info = SNESDefaultComputeJacobianColor(snes,V,HH,BB,flag,matfdcoloring);CHKERRQ(info);

255:   info = MatFDColoringDestroy(matfdcoloring);CHKERRQ(info);
256:   info = SNESDestroy(snes);CHKERRQ(info);
257:   
258:   info = VecDestroy(G);CHKERRQ(info);
259:   return(0);
260: }

265: /* @
266:   TaoAppSetFiniteDifferencesOptions - Sets various TAO parameters from user options

268:    Collective on TAO_APPLICATION

270:    Input Parameters:
271: +  taoapp - the TAO Application (optional)

273:    Level: beginner

275: .keywords:  options, finite differences

277: .seealso: TaoSolveApplication();

279: @ */
280: int TaoAppSetFiniteDifferencesOptions(TAO_APPLICATION taoapp){
281:   int info;
282:   PetscTruth flg;



288:   flg=PETSC_FALSE;
289:   info = PetscOptionsName("-tao_fd","use finite differences for Hessian","TaoAppDefaultComputeHessian",&flg);CHKERRQ(info);
290:   if (flg) {
291:     info = TaoAppSetHessianRoutine(taoapp,TaoAppDefaultComputeHessian,(void*)taoapp);CHKERRQ(info);
292:     info = PetscInfo(taoapp,"Setting default finite difference Hessian matrix\n"); CHKERRQ(info);
293:   }

295:   flg=PETSC_FALSE;
296:   info = PetscOptionsName("-tao_fdgrad","use finite differences for gradient","TaoAppDefaultComputeGradient",&flg);CHKERRQ(info);
297:   if (flg) {
298:     info = TaoAppSetGradientRoutine(taoapp,TaoAppDefaultComputeGradient,(void*)taoapp);CHKERRQ(info);
299:     info = PetscInfo(taoapp,"Setting default finite difference gradient routine\n"); CHKERRQ(info);
300:   }


303:   flg=PETSC_FALSE;
304:   info = PetscOptionsName("-tao_fd_coloring","use finite differences with coloring to compute Hessian","TaoAppDefaultComputeHessianColor",&flg);CHKERRQ(info);
305:   if (flg) {
306:     info = TaoAppSetHessianRoutine(taoapp,TaoAppDefaultComputeHessianColor,(void*)taoapp);CHKERRQ(info);
307:     info = PetscInfo(taoapp,"Use finite differencing with coloring to compute Hessian \n"); CHKERRQ(info);
308:   }
309:     
310:   return(0);
311: }


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

316: typedef struct {
317:   ISColoring coloring;
318: } TaoAppColorit;

322: static int TaoAppDestroyColoringXXX(void*ctx){
323:   int info;
324:   TaoAppColorit *mctx=(TaoAppColorit*)ctx;
326:   if (mctx){
327:     /*
328:     if (mctx->coloring){  
329:       info = ISColoringDestroy(mctx->coloring);CHKERRQ(info);
330:     }
331:     */
332:     info = PetscFree(mctx); CHKERRQ(info);
333:   }
334:   return(0);
335: }

339: /*@
340:    TaoAppSetColoring - Set the matrix coloring to be used when computing the
341:    Hessian by finite differences.

343:    Collective on TAO_APPLICATION

345:    Input Parameters:
346: +  app - the TAO_APPLICATION context
347: -  coloring - the coloring to be used

349:    Level: intermediate

351:    Note:
352:    The gradient evaluation must be set using the routine TaoSetPetscGradient().

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

356: .seealso: TaoAppSetHessianRoutine(), TaoAppDefaultComputeHessianColor(),
357:           TaoAppSetGradientRoutine()

359: @*/
360: int TaoAppSetColoring(TAO_APPLICATION taoapp, ISColoring coloring){
361:   int info;
362:   PetscInt ii;
363:   TaoAppColorit *mctx;
366:   
367:   info = TaoAppQueryForObject(taoapp,TaoAppColoringXXX,(void**)&mctx); CHKERRQ(info);
368:   if (mctx==0){
369:     info=PetscNew(TaoAppColorit,(void**)&mctx); CHKERRQ(info);
370:     info=TaoAppAddObject(taoapp,TaoAppColoringXXX,(void*)mctx,&ii); CHKERRQ(info);
371:     info=TaoAppSetDestroyRoutine(taoapp,TaoAppDestroyColoringXXX, (void*)mctx); CHKERRQ(info);
372:   }
373:   /*
374:   if (coloring){
375:   info=PetscObjectReference((PetscObject)coloring); CHKERRQ(info);
376:   }
377:   if (mctx->coloring){
378:        info=ISColoringDestroy(mctx->coloring); CHKERRQ(info);
379:   }
380:   */
381:   mctx->coloring=coloring;  
382:   return(0);
383: }

387: /*@C
388:    TaoAppGetColoring - Set the matrix coloring to be used when computing the
389:    Hessian by finite differences.

391:    Collective on TAO_APPLICATION

393:    Input Parameters:
394: +  app - the TAO_APPLICATION context
395: -  coloring - the coloring to be used

397:    Level: advanced

399:    Note:
400:    The gradient evaluation must be set using the routine TaoAppSetGradientRoutine().

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

404: .seealso: TaoAppSetHessianRoutine(), TaoAppDefaultComputeHessianColor(),
405:           TaoAppSetGradientRoutine()

407: @*/
408: int TaoAppGetColoring(TAO_APPLICATION taoapp, ISColoring *coloring){
409:   int info;
410:   TaoAppColorit *mctx;
413:   if (coloring){
414:     info = TaoAppQueryForObject(taoapp,TaoAppColoringXXX,(void**)&mctx); CHKERRQ(info);
415:     if (mctx){
416:       *coloring=mctx->coloring;
417:     } else {
418:       *coloring=0;
419:     }
420:   }
421:   return(0);
422: }


427: int TaoAppAddFiniteDifferences(TAO_APPLICATION taoapp){
428:   int info;
430:   info = TaoAppSetOptionsRoutine(taoapp,TaoAppSetFiniteDifferencesOptions); CHKERRQ(info);
431:   return(0);
432: }

434: int MatTAOMFSetGradient(Mat mat,Vec v,int (*func)(TAO_APPLICATION,Vec,Vec,void *),void *funcctx){
436:   return(0);
437: }