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 finite difference computations of the Hessian, we use PETSc's SNESComputeJacobianDefault
      8: */
      9: static PetscErrorCode Fsnes(SNES snes, Vec X, Vec G, PetscCtx 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 using the command-line option `-tao_test_gradient`
     44:   This finite difference gradient evaluation can be set using the routine `TaoSetGradient()` or by using the command line option -tao_fd_gradient

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

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

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

     87:   Collective

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

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

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

    101:   Level: advanced

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

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

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

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

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

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

    148:   Collective

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

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

    159:   Level: advanced

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

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

    178: /*@
    179:   TaoDefaultComputeHessianMFFD - Computes the Hessian using finite differences with `MATMFFD`.

    181:   Collective

    183:   Input Parameters:
    184: + tao - the `Tao` context
    185: . X   - compute Hessian at this point
    186: - ctx - ignored

    188:   Output Parameters:
    189: + H - Hessian matrix of type `MATMFFD`
    190: - B - should be `NULL` or equal to `H`

    192:   Level: advanced

    194:   Note:
    195:   This can be passed to `TaoSetHessian()` to use `MATMFFD` for approximate Hessian-vector products.  The matrix `H` can originate from
    196:   `MatCreateMFFD()` or from `TaoTermCreateHessianMFFD()`.

    198: .seealso: `Tao`, `MATMFFD`, `MatCreateMFFD()`, `TaoTermCreateHessianMFFD()`
    199: @*/
    200: PetscErrorCode TaoDefaultComputeHessianMFFD(Tao tao, Vec X, Mat H, Mat B, PetscCtx ctx)
    201: {
    202:   PetscInt  n, N;
    203:   PetscBool assembled;

    205:   PetscFunctionBegin;
    206:   PetscCheck(!B || B == H, PetscObjectComm((PetscObject)tao), PETSC_ERR_SUP, "Preconditioning Hessian matrix");
    207:   PetscCall(MatAssembled(H, &assembled));
    208:   if (!assembled) {
    209:     PetscCall(VecGetSize(X, &N));
    210:     PetscCall(VecGetLocalSize(X, &n));
    211:     PetscCall(MatSetSizes(H, n, n, N, N));
    212:     PetscCall(MatSetType(H, MATMFFD));
    213:     PetscCall(MatSetUp(H));
    214:     PetscCall(MatMFFDSetFunction(H, (PetscErrorCode (*)(void *, Vec, Vec))TaoComputeGradient, tao));
    215:   }
    216:   PetscCall(MatMFFDSetBase(H, X, NULL));
    217:   PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
    218:   PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
    219:   PetscFunctionReturn(PETSC_SUCCESS);
    220: }