Actual source code: fdiff.c
1: #include "taosolver.h" /*I "taosolver.h" I*/
2: #include "private/taosolver_impl.h"
3: #include "petscsnes.h"
6: /*
7: For finited difference computations of the Hessian, we use PETSc's SNESDefaultComputeJacobian
8: */
12: static PetscErrorCode Fsnes(SNES snes ,Vec X,Vec G,void*ctx){
14: TaoSolver tao = (TaoSolver)ctx;
17: ierr=TaoComputeGradient(tao,X,G);
18: return(0);
19: }
23: /*@C
24: TaoDefaultComputeGradient - computes the gradient using finite differences.
25:
26: Collective on TaoSolver
28: Input Parameters:
29: + tao - the TaoSolver context
30: . X - compute gradient at this point
31: - dummy - not used
33: Output Parameters:
34: . G - Gradient Vector
36: Options Database Key:
37: + -tao_fd_gradient - Activates TaoDefaultComputeGradient()
38: - -tao_fd_delta <delta> - change in x used to calculate finite differences
43: Level: advanced
45: Note:
46: This routine is slow and expensive, and is not currently optimized
47: to take advantage of sparsity in the problem. Although
48: TaoAppDefaultComputeGradient is not recommended for general use
49: in large-scale applications, It can be useful in checking the
50: correctness of a user-provided gradient. Use the tao method "tao_fd_test"
51: to get an indication of whether your gradient is correct.
54: Note:
55: This finite difference gradient evaluation can be set using the routine TaoSetGradientRoutine() or by using the command line option -tao_fd_gradient
57: .seealso: TaoSetGradientRoutine()
59: @*/
60: PetscErrorCode TaoDefaultComputeGradient(TaoSolver tao,Vec X,Vec G,void *dummy)
61: {
62: PetscReal *g;
63: PetscReal f, f2;
65: PetscInt low,high,N,i;
66: PetscBool flg;
67: PetscReal h=PETSC_SQRT_MACHINE_EPSILON;
69: TaoComputeObjective(tao, X,&f);
70: PetscOptionsGetReal(PETSC_NULL,"-tao_fd_delta",&h,&flg);
71: VecGetSize(X,&N);
72: VecGetOwnershipRange(X,&low,&high);
73: VecGetArray(G,&g);
74: for (i=0;i<N;i++) {
75: VecSetValue(X,i,h,ADD_VALUES);
76: VecAssemblyBegin(X);
77: VecAssemblyEnd(X);
79: TaoComputeObjective(tao,X,&f2);
81: VecSetValue(X,i,-h,ADD_VALUES);
82: VecAssemblyBegin(X);
83: VecAssemblyEnd(X);
84:
85: if (i>=low && i<high) {
86: g[i-low]=(f2-f)/h;
87: }
88: }
89: VecRestoreArray(G,&g);
90: return(0);
91: }
95: /*@C
96: TaoDefaultComputeHessian - Computes the Hessian using finite differences.
98: Collective on TaoSolver
100: Input Parameters:
101: + tao - the TaoSolver context
102: . V - compute Hessian at this point
103: - dummy - not used
105: Output Parameters:
106: + H - Hessian matrix (not altered in this routine)
107: . B - newly computed Hessian matrix to use with preconditioner (generally the same as H)
108: - flag - flag indicating whether the matrix sparsity structure has changed
110: Options Database Key:
111: + -tao_fd - Activates TaoDefaultComputeHessian()
112: - -tao_view_hessian - view the hessian after each evaluation using PETSC_VIEWER_STDOUT_WORLD
114: Level: advanced
116: Notes:
117: This routine is slow and expensive, and is not currently optimized
118: to take advantage of sparsity in the problem. Although
119: TaoDefaultComputeHessian() is not recommended for general use
120: in large-scale applications, It can be useful in checking the
121: correctness of a user-provided Hessian.
125: .seealso: TaoSetHessianRoutine(), TaoDefaultComputeHessianColor(), SNESDefaultComputeJacobian(), TaoSetGradientRoutine(), TaoDefaultComputeGradient()
127: @*/
128: PetscErrorCode TaoDefaultComputeHessian(TaoSolver tao,Vec V,Mat *H,Mat *B,
129: MatStructure *flag,void *dummy){
130: PetscErrorCode ierr;
131: MPI_Comm comm;
132: Vec G;
133: SNES snes;
137: VecDuplicate(V,&G);
139: PetscInfo(tao,"TAO Using finite differences w/o coloring to compute Hessian matrix\n");
141: TaoComputeGradient(tao,V,G);
143: PetscObjectGetComm((PetscObject)(*H),&comm);
144: SNESCreate(comm,&snes);
146: SNESSetFunction(snes,G,Fsnes,tao);
147: SNESDefaultComputeJacobian(snes,V,H,B,flag,tao);
149: SNESDestroy(&snes);
150:
151: VecDestroy(&G);
152:
153: return(0);
154: }
161: /*@C
162: TaoDefaultComputeHessianColor - Computes the Hessian using colored finite differences.
164: Collective on TaoSolver
166: Input Parameters:
167: + tao - the TaoSolver context
168: . V - compute Hessian at this point
169: - ctx - the PetscColoring object (must be of type MatFDColoring)
171: Output Parameters:
172: + H - Hessian matrix (not altered in this routine)
173: . B - newly computed Hessian matrix to use with preconditioner (generally the same as H)
174: - flag - flag indicating whether the matrix sparsity structure has changed
176: Level: advanced
179: .seealso: TaoSetHessianRoutine(), TaoDefaultComputeHessian(),SNESDefaultComputeJacobianColor(), TaoSetGradientRoutine()
181: @*/
182: PetscErrorCode TaoDefaultComputeHessianColor(TaoSolver tao, Vec V, Mat *H,Mat *B,MatStructure *flag,void *ctx){
183: PetscErrorCode ierr;
184: MatFDColoring coloring = (MatFDColoring)ctx;
189:
190:
191: *flag = SAME_NONZERO_PATTERN;
193: ierr=PetscInfo(tao,"TAO computing matrix using finite differences Hessian and coloring\n");
194: MatFDColoringApply(*B,coloring,V,flag,ctx);
196: if (*H != *B) {
197: MatAssemblyBegin(*H, MAT_FINAL_ASSEMBLY);
198: MatAssemblyEnd(*H, MAT_FINAL_ASSEMBLY);
199: }
200: return(0);
201: }