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