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