1: #include <petsctao.h>         /*I  "petsctao.h"  I*/
  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: */

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

 18:   ierr=TaoComputeGradient(tao,X,G);
 19:   return(0);
 20: }

 24: /*@C
 25:   TaoDefaultComputeGradient - computes the gradient using finite differences.

 27:   Collective on Tao

 29:   Input Parameters:
 30: + tao - the Tao context
 31: . X - compute gradient at this point
 32: - dummy - not used

 34:   Output Parameters:
 35: . G - Gradient Vector

 37:    Options Database Key:
 38: +  -tao_fd_gradient - Activates TaoDefaultComputeGradient()
 39: -  -tao_fd_delta <delta> - change in x used to calculate finite differences

 41:    Level: advanced

 43:    Note:
 44:    This routine is slow and expensive, and is not currently optimized
 45:    to take advantage of sparsity in the problem.  Although
 46:    TaoAppDefaultComputeGradient is not recommended for general use
 47:    in large-scale applications, It can be useful in checking the
 48:    correctness of a user-provided gradient.  Use the tao method TAOTEST
 49:    to get an indication of whether your gradient is correct.


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

 55: .seealso: TaoSetGradientRoutine()

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

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

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

 81:     if (i>=low && i<high) {
 82:       VecGetArray(X,&x);
 83:       x[i-low] += 2*h;
 84:       VecRestoreArray(X,&x);
 85:     }

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

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

104: /*@C
105:    TaoDefaultComputeHessian - Computes the Hessian using finite differences.

107:    Collective on Tao

109:    Input Parameters:
110: +  tao - the Tao context
111: .  V - compute Hessian at this point
112: -  dummy - not used

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

118:    Options Database Key:
119: +  -tao_fd - Activates TaoDefaultComputeHessian()
120: -  -tao_view_hessian - view the hessian after each evaluation using PETSC_VIEWER_STDOUT_WORLD

122:    Level: advanced

124:    Notes:
125:    This routine is slow and expensive, and is not currently optimized
126:    to take advantage of sparsity in the problem.  Although
127:    TaoDefaultComputeHessian() is not recommended for general use
128:    in large-scale applications, It can be useful in checking the
129:    correctness of a user-provided Hessian.

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

133: @*/
134: PetscErrorCode TaoDefaultComputeHessian(Tao tao,Vec V,Mat H,Mat B,void *dummy)
135: {
136:   PetscErrorCode       ierr;
137:   MPI_Comm             comm;
138:   Vec                  G;
139:   SNES                 snes;

143:   VecDuplicate(V,&G);

145:   PetscInfo(tao,"TAO Using finite differences w/o coloring to compute Hessian matrix\n");

147:   TaoComputeGradient(tao,V,G);

149:   PetscObjectGetComm((PetscObject)H,&comm);
150:   SNESCreate(comm,&snes);

152:   SNESSetFunction(snes,G,Fsnes,tao);
153:   SNESComputeJacobianDefault(snes,V,H,B,tao);
154:   SNESDestroy(&snes);
155:   VecDestroy(&G);
156:   return(0);
157: }

161: /*@C
162:    TaoDefaultComputeHessianColor - Computes the Hessian using colored finite differences.

164:    Collective on Tao

166:    Input Parameters:
167: +  tao - the Tao 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)

175:    Level: advanced


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

180: @*/
181: PetscErrorCode TaoDefaultComputeHessianColor(Tao tao, Vec V, Mat H,Mat B,void *ctx)
182: {
183:   PetscErrorCode      ierr;
184:   MatFDColoring       coloring = (MatFDColoring)ctx;

188:   ierr=PetscInfo(tao,"TAO computing matrix using finite differences Hessian and coloring\n");
189:   MatFDColoringApply(B,coloring,V,ctx);
190:   if (H != B) {
191:     MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
192:     MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
193:   }
194:   return(0);
195: }