Actual source code: fdiff.c
1: #include <petsctao.h>
2: #include <petsc/private/taoimpl.h>
3: #include <petscsnes.h>
4: #include <petscdmshell.h>
6: /*
7: For finited difference computations of the Hessian, we use PETSc's SNESComputeJacobianDefault
8: */
9: static PetscErrorCode Fsnes(SNES snes,Vec X,Vec G,void* ctx)
10: {
11: Tao tao = (Tao)ctx;
14: TaoComputeGradient(tao,X,G);
15: return 0;
16: }
18: /*@C
19: TaoDefaultComputeGradient - computes the gradient using finite differences.
21: Collective on Tao
23: Input Parameters:
24: + tao - the Tao context
25: . X - compute gradient at this point
26: - dummy - not used
28: Output Parameter:
29: . G - Gradient Vector
31: Options Database Key:
32: + -tao_fd_gradient - activates TaoDefaultComputeGradient()
33: - -tao_fd_delta <delta> - change in X used to calculate finite differences
35: Level: advanced
37: Notes:
38: This routine is slow and expensive, and is not currently optimized
39: to take advantage of sparsity in the problem. Although
40: TaoDefaultComputeGradient is not recommended for general use
41: in large-scale applications, It can be useful in checking the
42: correctness of a user-provided gradient. Use the tao method TAOTEST
43: to get an indication of whether your gradient is correct.
44: This finite difference gradient evaluation can be set using the routine TaoSetGradient() or by using the command line option -tao_fd_gradient
46: .seealso: TaoSetGradient()
47: @*/
48: PetscErrorCode TaoDefaultComputeGradient(Tao tao,Vec Xin,Vec G,void *dummy)
49: {
50: Vec X;
51: PetscScalar *g;
52: PetscReal f, f2;
53: PetscInt low,high,N,i;
54: PetscBool flg;
55: PetscReal h=.5*PETSC_SQRT_MACHINE_EPSILON;
57: PetscOptionsGetReal(((PetscObject)tao)->options,((PetscObject)tao)->prefix,"-tao_fd_delta",&h,&flg);
58: VecDuplicate(Xin,&X);
59: VecCopy(Xin,X);
60: VecGetSize(X,&N);
61: VecGetOwnershipRange(X,&low,&high);
62: VecSetOption(X,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
63: VecGetArray(G,&g);
64: for (i=0;i<N;i++) {
65: VecSetValue(X,i,-h,ADD_VALUES);
66: VecAssemblyBegin(X);
67: VecAssemblyEnd(X);
68: TaoComputeObjective(tao,X,&f);
69: VecSetValue(X,i,2.0*h,ADD_VALUES);
70: VecAssemblyBegin(X);
71: VecAssemblyEnd(X);
72: TaoComputeObjective(tao,X,&f2);
73: VecSetValue(X,i,-h,ADD_VALUES);
74: VecAssemblyBegin(X);
75: VecAssemblyEnd(X);
76: if (i>=low && i<high) {
77: g[i-low]=(f2-f)/(2.0*h);
78: }
79: }
80: VecRestoreArray(G,&g);
81: VecDestroy(&X);
82: return 0;
83: }
85: /*@C
86: TaoDefaultComputeHessian - Computes the Hessian using finite differences.
88: Collective on Tao
90: Input Parameters:
91: + tao - the Tao context
92: . V - compute Hessian at this point
93: - dummy - not used
95: Output Parameters:
96: + H - Hessian matrix (not altered in this routine)
97: - B - newly computed Hessian matrix to use with preconditioner (generally the same as H)
99: Options Database Key:
100: . -tao_fd_hessian - activates TaoDefaultComputeHessian()
102: Level: advanced
104: Notes:
105: This routine is slow and expensive, and is not currently optimized
106: to take advantage of sparsity in the problem. Although
107: TaoDefaultComputeHessian() is not recommended for general use
108: in large-scale applications, It can be useful in checking the
109: correctness of a user-provided Hessian.
111: .seealso: TaoSetHessian(), TaoDefaultComputeHessianColor(), SNESComputeJacobianDefault(), TaoSetGradient(), TaoDefaultComputeGradient()
112: @*/
113: PetscErrorCode TaoDefaultComputeHessian(Tao tao,Vec V,Mat H,Mat B,void *dummy)
114: {
115: SNES snes;
116: DM dm;
118: PetscInfo(tao,"TAO Using finite differences w/o coloring to compute Hessian matrix\n");
119: SNESCreate(PetscObjectComm((PetscObject)H),&snes);
120: SNESSetFunction(snes,NULL,Fsnes,tao);
121: SNESGetDM(snes,&dm);
122: DMShellSetGlobalVector(dm,V);
123: SNESSetUp(snes);
124: if (H) {
125: PetscInt n,N;
127: VecGetSize(V,&N);
128: VecGetLocalSize(V,&n);
129: MatSetSizes(H,n,n,N,N);
130: MatSetUp(H);
131: }
132: if (B && B != H) {
133: PetscInt n,N;
135: VecGetSize(V,&N);
136: VecGetLocalSize(V,&n);
137: MatSetSizes(B,n,n,N,N);
138: MatSetUp(B);
139: }
140: SNESComputeJacobianDefault(snes,V,H,B,NULL);
141: SNESDestroy(&snes);
142: return 0;
143: }
145: /*@C
146: TaoDefaultComputeHessianColor - Computes the Hessian using colored finite differences.
148: Collective on Tao
150: Input Parameters:
151: + tao - the Tao context
152: . V - compute Hessian at this point
153: - ctx - the PetscColoring object (must be of type MatFDColoring)
155: Output Parameters:
156: + H - Hessian matrix (not altered in this routine)
157: - B - newly computed Hessian matrix to use with preconditioner (generally the same as H)
159: Level: advanced
161: .seealso: TaoSetHessian(), TaoDefaultComputeHessian(),SNESComputeJacobianDefaultColor(), TaoSetGradient()
162: @*/
163: PetscErrorCode TaoDefaultComputeHessianColor(Tao tao,Vec V,Mat H,Mat B,void *ctx)
164: {
165: MatFDColoring coloring = (MatFDColoring)ctx;
168: PetscInfo(tao,"TAO computing matrix using finite differences Hessian and coloring\n");
169: MatFDColoringApply(B,coloring,V,ctx);
170: if (H != B) {
171: MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
172: MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
173: }
174: return 0;
175: }
177: PetscErrorCode TaoDefaultComputeHessianMFFD(Tao tao,Vec X,Mat H,Mat B,void *ctx)
178: {
179: PetscInt n,N;
180: PetscBool assembled;
183: MatAssembled(H, &assembled);
184: if (!assembled) {
185: VecGetSize(X,&N);
186: VecGetLocalSize(X,&n);
187: MatSetSizes(H,n,n,N,N);
188: MatSetType(H,MATMFFD);
189: MatSetUp(H);
190: MatMFFDSetFunction(H,(PetscErrorCode (*)(void*,Vec,Vec))TaoComputeGradient,tao);
191: }
192: MatMFFDSetBase(H,X,NULL);
193: MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
194: MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
195: return 0;
196: }