1:  #include <petsctao.h>
  2:  #include <petsc/private/taoimpl.h>
  3:  #include <petscsnes.h>

  5: /*
  6:    For finited difference computations of the Hessian, we use PETSc's SNESComputeJacobianDefault
  7: */

  9: static PetscErrorCode Fsnes(SNES snes ,Vec X,Vec G,void*ctx)
 10: {
 12:   Tao            tao = (Tao)ctx;

 16:   ierr=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:    Note:
 40:    This routine is slow and expensive, and is not currently optimized
 41:    to take advantage of sparsity in the problem.  Although
 42:    TaoAppDefaultComputeGradient 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.


 48:    Note:
 49:    This finite difference gradient evaluation can be set using the routine TaoSetGradientRoutine() or by using the command line option -tao_fd_gradient

 51: .seealso: TaoSetGradientRoutine()

 53: @*/
 54: PetscErrorCode TaoDefaultComputeGradient(Tao tao,Vec X,Vec G,void *dummy)
 55: {
 56:   PetscScalar    *x,*g;
 57:   PetscReal      f, f2;
 59:   PetscInt       low,high,N,i;
 60:   PetscBool      flg;
 61:   PetscReal      h=.5*PETSC_SQRT_MACHINE_EPSILON;

 64:   PetscOptionsGetReal(((PetscObject)tao)->options,((PetscObject)tao)->prefix,"-tao_fd_delta",&h,&flg);
 65:   VecGetSize(X,&N);
 66:   VecGetOwnershipRange(X,&low,&high);
 67:   VecGetArray(G,&g);
 68:   for (i=0;i<N;i++) {
 69:     if (i>=low && i<high) {
 70:       VecGetArray(X,&x);
 71:       x[i-low] -= h;
 72:       VecRestoreArray(X,&x);
 73:     }

 75:     TaoComputeObjective(tao, X,&f);

 77:     if (i>=low && i<high) {
 78:       VecGetArray(X,&x);
 79:       x[i-low] += 2*h;
 80:       VecRestoreArray(X,&x);
 81:     }

 83:     TaoComputeObjective(tao,X,&f2);

 85:     if (i>=low && i<high) {
 86:       VecGetArray(X,&x);
 87:       x[i-low] -= h;
 88:       VecRestoreArray(X,&x);
 89:     }
 90:     if (i>=low && i<high) {
 91:       g[i-low]=(f2-f)/(2.0*h);
 92:     }
 93:   }
 94:   VecRestoreArray(G,&g);
 95:   return(0);
 96: }

 98: /*@C
 99:    TaoDefaultComputeHessian - Computes the Hessian using finite differences.

101:    Collective on Tao

103:    Input Parameters:
104: +  tao - the Tao context
105: .  V - compute Hessian at this point
106: -  dummy - not used

108:    Output Parameters:
109: +  H - Hessian matrix (not altered in this routine)
110: -  B - newly computed Hessian matrix to use with preconditioner (generally the same as H)

112:    Options Database Key:
113: +  -tao_fd - Activates TaoDefaultComputeHessian()
114: -  -tao_view_hessian - view the hessian after each evaluation using PETSC_VIEWER_STDOUT_WORLD

116:    Level: advanced

118:    Notes:
119:    This routine is slow and expensive, and is not currently optimized
120:    to take advantage of sparsity in the problem.  Although
121:    TaoDefaultComputeHessian() is not recommended for general use
122:    in large-scale applications, It can be useful in checking the
123:    correctness of a user-provided Hessian.

125: .seealso: TaoSetHessianRoutine(), TaoDefaultComputeHessianColor(), SNESComputeJacobianDefault(), TaoSetGradientRoutine(), TaoDefaultComputeGradient()

127: @*/
128: PetscErrorCode TaoDefaultComputeHessian(Tao tao,Vec V,Mat H,Mat B,void *dummy)
129: {
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:   SNESComputeJacobianDefault(snes,V,H,B,tao);
148:   SNESDestroy(&snes);
149:   VecDestroy(&G);
150:   return(0);
151: }

153: /*@C
154:    TaoDefaultComputeHessianColor - Computes the Hessian using colored finite differences.

156:    Collective on Tao

158:    Input Parameters:
159: +  tao - the Tao context
160: .  V - compute Hessian at this point
161: -  ctx - the PetscColoring object (must be of type MatFDColoring)

163:    Output Parameters:
164: +  H - Hessian matrix (not altered in this routine)
165: -  B - newly computed Hessian matrix to use with preconditioner (generally the same as H)

167:    Level: advanced


170: .seealso: TaoSetHessianRoutine(), TaoDefaultComputeHessian(),SNESComputeJacobianDefaultColor(), TaoSetGradientRoutine()

172: @*/
173: PetscErrorCode TaoDefaultComputeHessianColor(Tao tao, Vec V, Mat H,Mat B,void *ctx)
174: {
175:   PetscErrorCode      ierr;
176:   MatFDColoring       coloring = (MatFDColoring)ctx;

180:   ierr=PetscInfo(tao,"TAO computing matrix using finite differences Hessian and coloring\n");
181:   MatFDColoringApply(B,coloring,V,ctx);
182:   if (H != B) {
183:     MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
184:     MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
185:   }
186:   return(0);
187: }