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: {
     11:   Tao tao = (Tao)ctx;

     13:   PetscFunctionBegin;
     15:   PetscCall(TaoComputeGradient(tao, X, G));
     16:   PetscFunctionReturn(PETSC_SUCCESS);
     17: }

     19: /*@C
     20:   TaoDefaultComputeGradient - computes the gradient using finite differences.

     22:   Collective

     24:   Input Parameters:
     25: + tao   - the Tao context
     26: . Xin   - compute gradient at this point
     27: - dummy - not used

     29:   Output Parameter:
     30: . G - Gradient Vector

     32:   Options Database Key:
     33: + -tao_fd_gradient      - activates TaoDefaultComputeGradient()
     34: - -tao_fd_delta <delta> - change in X used to calculate finite differences

     36:   Level: advanced

     38:   Notes:
     39:   This routine is slow and expensive, and is not optimized
     40:   to take advantage of sparsity in the problem.  Although
     41:   not recommended for general use
     42:   in large-scale applications, it can be useful in checking the
     43:   correctness of a user-provided gradient.  Use the tao method TAOTEST
     44:   to get an indication of whether your gradient is correct.
     45:   This finite difference gradient evaluation can be set using the routine `TaoSetGradient()` or by using the command line option -tao_fd_gradient

     47: .seealso: `Tao`, `TaoSetGradient()`
     48: @*/
     49: PetscErrorCode TaoDefaultComputeGradient(Tao tao, Vec Xin, Vec G, void *dummy)
     50: {
     51:   Vec          X;
     52:   PetscScalar *g;
     53:   PetscReal    f, f2;
     54:   PetscInt     low, high, N, i;
     55:   PetscBool    flg;
     56:   PetscReal    h = .5 * PETSC_SQRT_MACHINE_EPSILON;

     58:   PetscFunctionBegin;
     59:   PetscCall(PetscOptionsGetReal(((PetscObject)tao)->options, ((PetscObject)tao)->prefix, "-tao_fd_delta", &h, &flg));
     60:   PetscCall(VecDuplicate(Xin, &X));
     61:   PetscCall(VecCopy(Xin, X));
     62:   PetscCall(VecGetSize(X, &N));
     63:   PetscCall(VecGetOwnershipRange(X, &low, &high));
     64:   PetscCall(VecSetOption(X, VEC_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE));
     65:   PetscCall(VecGetArray(G, &g));
     66:   for (i = 0; i < N; i++) {
     67:     PetscCall(VecSetValue(X, i, -h, ADD_VALUES));
     68:     PetscCall(VecAssemblyBegin(X));
     69:     PetscCall(VecAssemblyEnd(X));
     70:     PetscCall(TaoComputeObjective(tao, X, &f));
     71:     PetscCall(VecSetValue(X, i, 2.0 * h, ADD_VALUES));
     72:     PetscCall(VecAssemblyBegin(X));
     73:     PetscCall(VecAssemblyEnd(X));
     74:     PetscCall(TaoComputeObjective(tao, X, &f2));
     75:     PetscCall(VecSetValue(X, i, -h, ADD_VALUES));
     76:     PetscCall(VecAssemblyBegin(X));
     77:     PetscCall(VecAssemblyEnd(X));
     78:     if (i >= low && i < high) g[i - low] = (f2 - f) / (2.0 * h);
     79:   }
     80:   PetscCall(VecRestoreArray(G, &g));
     81:   PetscCall(VecDestroy(&X));
     82:   PetscFunctionReturn(PETSC_SUCCESS);
     83: }

     85: /*@C
     86:   TaoDefaultComputeHessian - Computes the Hessian using finite differences.

     88:   Collective

     90:   Input Parameters:
     91: + tao   - the Tao context
     92: . V     - compute Hessian at this point
     93: - dummy - not used

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

     99:   Options Database Key:
    100: . -tao_fd_hessian - activates TaoDefaultComputeHessian()

    102:   Level: advanced

    104:   Notes:
    105:   This routine is slow and expensive, and is not optimized
    106:   to take advantage of sparsity in the problem.  Although
    107:   it is not recommended for general use
    108:   in large-scale applications, It can be useful in checking the
    109:   correctness of a user-provided Hessian.

    111: .seealso: `Tao`, `TaoSetHessian()`, `TaoDefaultComputeHessianColor()`, `SNESComputeJacobianDefault()`, `TaoSetGradient()`, `TaoDefaultComputeGradient()`
    112: @*/
    113: PetscErrorCode TaoDefaultComputeHessian(Tao tao, Vec V, Mat H, Mat B, void *dummy)
    114: {
    115:   SNES snes;
    116:   DM   dm;

    118:   PetscFunctionBegin;
    119:   PetscCall(PetscInfo(tao, "TAO Using finite differences w/o coloring to compute Hessian matrix\n"));
    120:   PetscCall(SNESCreate(PetscObjectComm((PetscObject)H), &snes));
    121:   PetscCall(SNESSetFunction(snes, NULL, Fsnes, tao));
    122:   PetscCall(SNESGetDM(snes, &dm));
    123:   PetscCall(DMShellSetGlobalVector(dm, V));
    124:   PetscCall(SNESSetUp(snes));
    125:   if (H) {
    126:     PetscInt n, N;

    128:     PetscCall(VecGetSize(V, &N));
    129:     PetscCall(VecGetLocalSize(V, &n));
    130:     PetscCall(MatSetSizes(H, n, n, N, N));
    131:     PetscCall(MatSetUp(H));
    132:   }
    133:   if (B && B != H) {
    134:     PetscInt n, N;

    136:     PetscCall(VecGetSize(V, &N));
    137:     PetscCall(VecGetLocalSize(V, &n));
    138:     PetscCall(MatSetSizes(B, n, n, N, N));
    139:     PetscCall(MatSetUp(B));
    140:   }
    141:   PetscCall(SNESComputeJacobianDefault(snes, V, H, B, NULL));
    142:   PetscCall(SNESDestroy(&snes));
    143:   PetscFunctionReturn(PETSC_SUCCESS);
    144: }

    146: /*@C
    147:   TaoDefaultComputeHessianColor - Computes the Hessian using colored finite differences.

    149:   Collective

    151:   Input Parameters:
    152: + tao - the Tao context
    153: . V   - compute Hessian at this point
    154: - ctx - the color object of type `MatFDColoring`

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

    160:   Level: advanced

    162: .seealso: `Tao`, `MatColoring`, `TaoSetHessian()`, `TaoDefaultComputeHessian()`, `SNESComputeJacobianDefaultColor()`, `TaoSetGradient()`
    163: @*/
    164: PetscErrorCode TaoDefaultComputeHessianColor(Tao tao, Vec V, Mat H, Mat B, void *ctx)
    165: {
    166:   MatFDColoring coloring = (MatFDColoring)ctx;

    168:   PetscFunctionBegin;
    170:   PetscCall(PetscInfo(tao, "TAO computing matrix using finite differences Hessian and coloring\n"));
    171:   PetscCall(MatFDColoringApply(B, coloring, V, ctx));
    172:   if (H != B) {
    173:     PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
    174:     PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
    175:   }
    176:   PetscFunctionReturn(PETSC_SUCCESS);
    177: }

    179: PetscErrorCode TaoDefaultComputeHessianMFFD(Tao tao, Vec X, Mat H, Mat B, void *ctx)
    180: {
    181:   PetscInt  n, N;
    182:   PetscBool assembled;

    184:   PetscFunctionBegin;
    185:   PetscCheck(!B || B == H, PetscObjectComm((PetscObject)tao), PETSC_ERR_SUP, "Preconditioning Hessian matrix");
    186:   PetscCall(MatAssembled(H, &assembled));
    187:   if (!assembled) {
    188:     PetscCall(VecGetSize(X, &N));
    189:     PetscCall(VecGetLocalSize(X, &n));
    190:     PetscCall(MatSetSizes(H, n, n, N, N));
    191:     PetscCall(MatSetType(H, MATMFFD));
    192:     PetscCall(MatSetUp(H));
    193:     PetscCall(MatMFFDSetFunction(H, (PetscErrorCode (*)(void *, Vec, Vec))TaoComputeGradient, tao));
    194:   }
    195:   PetscCall(MatMFFDSetBase(H, X, NULL));
    196:   PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
    197:   PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
    198:   PetscFunctionReturn(PETSC_SUCCESS);
    199: }