Actual source code: fdiff.c

  1: #include "taosolver.h"         /*I  "taosolver.h"  I*/
  2: #include "private/taosolver_impl.h"    
  3: #include "petscsnes.h"


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

 12: static PetscErrorCode Fsnes(SNES snes ,Vec X,Vec G,void*ctx){
 14:   TaoSolver tao = (TaoSolver)ctx;
 17:   ierr=TaoComputeGradient(tao,X,G); 
 18:   return(0);
 19: }

 23: /*@C
 24:   TaoDefaultComputeGradient - computes the gradient using finite differences.
 25:  
 26:   Collective on TaoSolver

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

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

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




 43:    Level: advanced

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


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

 57: .seealso: TaoSetGradientRoutine()

 59: @*/
 60: PetscErrorCode TaoDefaultComputeGradient(TaoSolver tao,Vec X,Vec G,void *dummy) 
 61: {
 62:   PetscReal *g;
 63:   PetscReal f, f2;
 65:   PetscInt low,high,N,i;
 66:   PetscBool flg;
 67:   PetscReal h=PETSC_SQRT_MACHINE_EPSILON;
 69:   TaoComputeObjective(tao, X,&f); 
 70:   PetscOptionsGetReal(PETSC_NULL,"-tao_fd_delta",&h,&flg); 
 71:   VecGetSize(X,&N); 
 72:   VecGetOwnershipRange(X,&low,&high); 
 73:   VecGetArray(G,&g); 
 74:   for (i=0;i<N;i++) {
 75:       VecSetValue(X,i,h,ADD_VALUES); 
 76:       VecAssemblyBegin(X); 
 77:       VecAssemblyEnd(X); 

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

 81:       VecSetValue(X,i,-h,ADD_VALUES);
 82:       VecAssemblyBegin(X); 
 83:       VecAssemblyEnd(X); 
 84:       
 85:       if (i>=low && i<high) {
 86:           g[i-low]=(f2-f)/h;
 87:       }
 88:   }
 89:   VecRestoreArray(G,&g); 
 90:   return(0);
 91: }

 95: /*@C
 96:    TaoDefaultComputeHessian - Computes the Hessian using finite differences. 

 98:    Collective on TaoSolver

100:    Input Parameters:
101: +  tao - the TaoSolver context 
102: .  V - compute Hessian at this point
103: -  dummy - not used

105:    Output Parameters:
106: +  H - Hessian matrix (not altered in this routine)
107: .  B - newly computed Hessian matrix to use with preconditioner (generally the same as H)
108: -  flag - flag indicating whether the matrix sparsity structure has changed

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

114:    Level: advanced

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



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

127: @*/
128: PetscErrorCode TaoDefaultComputeHessian(TaoSolver tao,Vec V,Mat *H,Mat *B,
129:                              MatStructure *flag,void *dummy){
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:   SNESDefaultComputeJacobian(snes,V,H,B,flag,tao);

149:   SNESDestroy(&snes);
150:   
151:   VecDestroy(&G);
152:   
153:   return(0);
154: }




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

164:    Collective on TaoSolver

166:    Input Parameters:
167: +  tao - the TaoSolver 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)
174: -  flag - flag indicating whether the matrix sparsity structure has changed

176:    Level: advanced


179: .seealso: TaoSetHessianRoutine(), TaoDefaultComputeHessian(),SNESDefaultComputeJacobianColor(), TaoSetGradientRoutine()

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


189:   
190:   
191:   *flag = SAME_NONZERO_PATTERN;

193:   ierr=PetscInfo(tao,"TAO computing matrix using finite differences Hessian and coloring\n"); 
194:   MatFDColoringApply(*B,coloring,V,flag,ctx); 

196:   if (*H != *B) {
197:       MatAssemblyBegin(*H, MAT_FINAL_ASSEMBLY); 
198:       MatAssemblyEnd(*H, MAT_FINAL_ASSEMBLY); 
199:   }
200:   return(0);
201: }