Actual source code: fdiff.c
petsc-3.11.4 2019-09-28
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: {
12: Tao tao = (Tao)ctx;
16: TaoComputeGradient(tao,X,G);
17: return(0);
18: }
20: /*@C
21: TaoDefaultComputeGradient - computes the gradient using finite differences.
23: Collective on Tao
25: Input Parameters:
26: + tao - the Tao context
27: . X - compute gradient at this point
28: - dummy - not used
30: Output Parameters:
31: . G - Gradient Vector
33: Options Database Key:
34: + -tao_fd_gradient - activates TaoDefaultComputeGradient()
35: - -tao_fd_delta <delta> - change in X used to calculate finite differences
37: Level: advanced
39: Notes:
40: This routine is slow and expensive, and is not currently optimized
41: to take advantage of sparsity in the problem. Although
42: TaoDefaultComputeGradient is not recommended for general use
43: in large-scale applications, It can be useful in checking the
44: correctness of a user-provided gradient. Use the tao method TAOTEST
45: to get an indication of whether your gradient is correct.
46: This finite difference gradient evaluation can be set using the routine TaoSetGradientRoutine() or by using the command line option -tao_fd_gradient
48: .seealso: TaoSetGradientRoutine()
49: @*/
50: PetscErrorCode TaoDefaultComputeGradient(Tao tao,Vec Xin,Vec G,void *dummy)
51: {
52: Vec X;
53: PetscScalar *g;
54: PetscReal f, f2;
56: PetscInt low,high,N,i;
57: PetscBool flg;
58: PetscReal h=.5*PETSC_SQRT_MACHINE_EPSILON;
61: PetscOptionsGetReal(((PetscObject)tao)->options,((PetscObject)tao)->prefix,"-tao_fd_delta",&h,&flg);
62: VecDuplicate(Xin,&X);
63: VecCopy(Xin,X);
64: VecGetSize(X,&N);
65: VecGetOwnershipRange(X,&low,&high);
66: VecSetOption(X,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
67: VecGetArray(G,&g);
68: for (i=0;i<N;i++) {
69: VecSetValue(X,i,-h,ADD_VALUES);
70: VecAssemblyBegin(X);
71: VecAssemblyEnd(X);
72: TaoComputeObjective(tao,X,&f);
73: VecSetValue(X,i,2.0*h,ADD_VALUES);
74: VecAssemblyBegin(X);
75: VecAssemblyEnd(X);
76: TaoComputeObjective(tao,X,&f2);
77: VecSetValue(X,i,-h,ADD_VALUES);
78: VecAssemblyBegin(X);
79: VecAssemblyEnd(X);
80: if (i>=low && i<high) {
81: g[i-low]=(f2-f)/(2.0*h);
82: }
83: }
84: VecRestoreArray(G,&g);
85: VecDestroy(&X);
86: return(0);
87: }
89: /*@C
90: TaoDefaultComputeHessian - Computes the Hessian using finite differences.
92: Collective on Tao
94: Input Parameters:
95: + tao - the Tao context
96: . V - compute Hessian at this point
97: - dummy - not used
99: Output Parameters:
100: + H - Hessian matrix (not altered in this routine)
101: - B - newly computed Hessian matrix to use with preconditioner (generally the same as H)
103: Options Database Key:
104: . -tao_fd_hessian - activates TaoDefaultComputeHessian()
106: Level: advanced
108: Notes:
109: This routine is slow and expensive, and is not currently optimized
110: to take advantage of sparsity in the problem. Although
111: TaoDefaultComputeHessian() is not recommended for general use
112: in large-scale applications, It can be useful in checking the
113: correctness of a user-provided Hessian.
115: .seealso: TaoSetHessianRoutine(), TaoDefaultComputeHessianColor(), SNESComputeJacobianDefault(), TaoSetGradientRoutine(), TaoDefaultComputeGradient()
116: @*/
117: PetscErrorCode TaoDefaultComputeHessian(Tao tao,Vec V,Mat H,Mat B,void *dummy)
118: {
120: Vec G;
121: SNES snes;
122: DM dm;
125: VecDuplicate(V,&G);
126: PetscInfo(tao,"TAO Using finite differences w/o coloring to compute Hessian matrix\n");
127: TaoComputeGradient(tao,V,G);
128: SNESCreate(PetscObjectComm((PetscObject)H),&snes);
129: SNESSetFunction(snes,G,Fsnes,tao);
130: SNESGetDM(snes,&dm);
131: DMShellSetGlobalVector(dm,V);
132: SNESSetUp(snes);
133: if (H) {
134: PetscInt n,N;
136: VecGetSize(V,&N);
137: VecGetLocalSize(V,&n);
138: MatSetSizes(H,n,n,N,N);
139: MatSetUp(H);
140: }
141: if (B && B != H) {
142: PetscInt n,N;
144: VecGetSize(V,&N);
145: VecGetLocalSize(V,&n);
146: MatSetSizes(B,n,n,N,N);
147: MatSetUp(B);
148: }
149: SNESComputeJacobianDefault(snes,V,H,B,NULL);
150: SNESDestroy(&snes);
151: VecDestroy(&G);
152: return(0);
153: }
155: /*@C
156: TaoDefaultComputeHessianColor - Computes the Hessian using colored finite differences.
158: Collective on Tao
160: Input Parameters:
161: + tao - the Tao context
162: . V - compute Hessian at this point
163: - ctx - the PetscColoring object (must be of type MatFDColoring)
165: Output Parameters:
166: + H - Hessian matrix (not altered in this routine)
167: - B - newly computed Hessian matrix to use with preconditioner (generally the same as H)
169: Level: advanced
172: .seealso: TaoSetHessianRoutine(), TaoDefaultComputeHessian(),SNESComputeJacobianDefaultColor(), TaoSetGradientRoutine()
173: @*/
174: PetscErrorCode TaoDefaultComputeHessianColor(Tao tao,Vec V,Mat H,Mat B,void *ctx)
175: {
176: PetscErrorCode ierr;
177: MatFDColoring coloring = (MatFDColoring)ctx;
181: PetscInfo(tao,"TAO computing matrix using finite differences Hessian and coloring\n");
182: MatFDColoringApply(B,coloring,V,ctx);
183: if (H != B) {
184: MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
185: MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
186: }
187: return(0);
188: }
190: PetscErrorCode TaoDefaultComputeHessianMFFD(Tao tao,Vec X,Mat H,Mat B,void *ctx)
191: {
192: PetscInt n,N;
196: if (B && B != H) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_SUP,"Preconditioning Hessian matrix");
197: VecGetSize(X,&N);
198: VecGetLocalSize(X,&n);
199: MatSetSizes(H,n,n,N,N);
200: MatSetType(H,MATMFFD);
201: MatSetUp(H);
202: MatMFFDSetBase(H,X,NULL);
203: MatMFFDSetFunction(H,(PetscErrorCode (*)(void*,Vec,Vec))TaoComputeGradient,tao);
204: MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
205: MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
206: return(0);
207: }