Actual source code: fdiff.c

petsc-master 2019-12-13
Report Typos and Errors
  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: }