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: {
     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:   SNES           snes;
    121:   DM             dm;

    124:   PetscInfo(tao,"TAO Using finite differences w/o coloring to compute Hessian matrix\n");
    125:   SNESCreate(PetscObjectComm((PetscObject)H),&snes);
    126:   SNESSetFunction(snes,NULL,Fsnes,tao);
    127:   SNESGetDM(snes,&dm);
    128:   DMShellSetGlobalVector(dm,V);
    129:   SNESSetUp(snes);
    130:   if (H) {
    131:     PetscInt n,N;

    133:     VecGetSize(V,&N);
    134:     VecGetLocalSize(V,&n);
    135:     MatSetSizes(H,n,n,N,N);
    136:     MatSetUp(H);
    137:   }
    138:   if (B && B != H) {
    139:     PetscInt n,N;

    141:     VecGetSize(V,&N);
    142:     VecGetLocalSize(V,&n);
    143:     MatSetSizes(B,n,n,N,N);
    144:     MatSetUp(B);
    145:   }
    146:   SNESComputeJacobianDefault(snes,V,H,B,NULL);
    147:   SNESDestroy(&snes);
    148:   return(0);
    149: }

    151: /*@C
    152:    TaoDefaultComputeHessianColor - Computes the Hessian using colored finite differences.

    154:    Collective on Tao

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

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

    165:    Level: advanced

    168: .seealso: TaoSetHessianRoutine(), TaoDefaultComputeHessian(),SNESComputeJacobianDefaultColor(), TaoSetGradientRoutine()
    169: @*/
    170: PetscErrorCode TaoDefaultComputeHessianColor(Tao tao,Vec V,Mat H,Mat B,void *ctx)
    171: {
    172:   PetscErrorCode      ierr;
    173:   MatFDColoring       coloring = (MatFDColoring)ctx;

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

    186: PetscErrorCode TaoDefaultComputeHessianMFFD(Tao tao,Vec X,Mat H,Mat B,void *ctx)
    187: {
    188:   PetscInt       n,N;

    192:   if (B && B != H) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_SUP,"Preconditioning Hessian matrix");
    193:   VecGetSize(X,&N);
    194:   VecGetLocalSize(X,&n);
    195:   MatSetSizes(H,n,n,N,N);
    196:   MatSetType(H,MATMFFD);
    197:   MatSetUp(H);
    198:   MatMFFDSetBase(H,X,NULL);
    199:   MatMFFDSetFunction(H,(PetscErrorCode (*)(void*,Vec,Vec))TaoComputeGradient,tao);
    200:   MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
    201:   MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
    202:   return(0);
    203: }