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