Actual source code: taotermfdiff.c

      1: #include <petsc/private/taoimpl.h>
      2: #include <petscsnes.h>
      3: #include <petscdmshell.h>

      5: typedef struct _n_TaoTermWithParameters {
      6:   TaoTerm term; // weak-reference
      7:   Vec     params;
      8: } TaoTermWithParameters;

     10: PETSC_INTERN PetscErrorCode TaoTermWithParametersDestroy(PetscCtxRt ctx)
     11: {
     12:   TaoTermWithParameters *t = (TaoTermWithParameters *)*(void **)ctx;

     14:   PetscFunctionBegin;
     15:   PetscCall(VecDestroy(&t->params));
     16:   PetscCall(PetscFree(t));
     17:   PetscFunctionReturn(PETSC_SUCCESS);
     18: }

     20: static PetscErrorCode SNESFunction_TaoTerm(SNES snes, Vec X, Vec G, void *ctx)
     21: {
     22:   TaoTermWithParameters *t = (TaoTermWithParameters *)ctx;

     24:   PetscFunctionBegin;
     26:   PetscCall(TaoTermComputeGradient(t->term, X, t->params, G));
     27:   PetscFunctionReturn(PETSC_SUCCESS);
     28: }

     30: /*@
     31:   TaoTermComputeGradientFD - Approximate the gradient of a `TaoTerm` using finite differences

     33:   Collective

     35:   Input Parameters:
     36: + term   - a `TaoTerm`
     37: . x      - a solution vector
     38: - params - parameters vector (may be `NULL`, see `TaoTermParametersMode`)

     40:   Output Parameter:
     41: . g - the computed finite difference approximation to the gradient

     43:   Options Database Keys:
     44: + -tao_term_fd_delta <delta>       - change in `x` used to calculate finite differences
     45: - -tao_term_gradient_use_fd <bool> - Use `TaoTermComputeGradientFD()` in `TaoTermComputeGradient()`

     47:   Level: advanced

     49:   Notes:
     50:   This routine is slow and expensive, and is not optimized to take advantage of
     51:   sparsity in the problem.  Although not recommended for general use in
     52:   large-scale applications, it can be useful in checking the correctness of a
     53:   user-provided gradient.  Call `TaoTermComputeGradientSetUseFD()` to start using
     54:   this routine in `TaoTermComputeGradient()`.

     56: .seealso: [](sec_tao_term),
     57:           `TaoTerm`,
     58:           `TaoTermGetFDDelta()`,
     59:           `TaoTermSetFDDelta()`,
     60:           `TaoTermComputeGradientSetUseFD()`,
     61:           `TaoTermComputeGradientGetUseFD()`,
     62:           `TaoTermComputeHessianFD()`
     63: @*/
     64: PetscErrorCode TaoTermComputeGradientFD(TaoTerm term, Vec x, Vec params, Vec g)
     65: {
     66:   Vec          x_perturbed;
     67:   PetscScalar *_g;
     68:   PetscReal    f, f2;
     69:   PetscInt     low, high, N, i;
     70:   PetscReal    h;

     72:   PetscFunctionBegin;
     77:   h = term->fd_delta;

     79:   PetscCall(VecDuplicate(x, &x_perturbed));
     80:   PetscCall(VecCopy(x, x_perturbed));
     81:   PetscCall(VecGetSize(x_perturbed, &N));
     82:   PetscCall(VecGetOwnershipRange(x_perturbed, &low, &high));
     83:   PetscCall(VecSetOption(x_perturbed, VEC_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE));
     84:   PetscCall(VecGetArray(g, &_g));
     85:   for (i = 0; i < N; i++) {
     86:     PetscCall(VecSetValue(x_perturbed, i, -h, ADD_VALUES));
     87:     PetscCall(VecAssemblyBegin(x_perturbed));
     88:     PetscCall(VecAssemblyEnd(x_perturbed));
     89:     PetscCall(TaoTermComputeObjective(term, x_perturbed, params, &f));
     90:     PetscCall(VecSetValue(x_perturbed, i, 2.0 * h, ADD_VALUES));
     91:     PetscCall(VecAssemblyBegin(x_perturbed));
     92:     PetscCall(VecAssemblyEnd(x_perturbed));
     93:     PetscCall(TaoTermComputeObjective(term, x_perturbed, params, &f2));
     94:     PetscCall(VecSetValue(x_perturbed, i, -h, ADD_VALUES));
     95:     PetscCall(VecAssemblyBegin(x_perturbed));
     96:     PetscCall(VecAssemblyEnd(x_perturbed));
     97:     if (i >= low && i < high) _g[i - low] = (f2 - f) / (2.0 * h);
     98:   }
     99:   PetscCall(VecRestoreArray(g, &_g));
    100:   PetscCall(VecDestroy(&x_perturbed));
    101:   PetscFunctionReturn(PETSC_SUCCESS);
    102: }

    104: /*@
    105:   TaoTermComputeHessianFD - Use finite difference to compute Hessian matrix.

    107:   Collective

    109:   Input Parameters:
    110: + term   - a `TaoTerm`
    111: . x      - a solution vector
    112: - params - parameters vector (may be `NULL`, see `TaoTermParametersMode`)

    114:   Output Parameters:
    115: + H    - (optional) Hessian matrix
    116: - Hpre - (optional) Hessian preconditioning matrix

    118:   Options Database Keys:
    119: + -tao_term_fd_delta <delta>      - change in X used to calculate finite differences
    120: - -tao_term_hessian_use_fd <bool> - Use `TaoTermComputeHessianFD()` in `TaoTermComputeHessian()`

    122:   Level: advanced

    124:   Notes:
    125:   This routine is slow and expensive, and is not optimized to take advantage of
    126:   sparsity in the problem.  Although not recommended for general use in
    127:   large-scale applications, it can be useful in checking the correctness of a
    128:   user-provided Hessian.  Call `TaoTermComputeHessianSetUseFD()` to start using
    129:   this routine in `TaoTermComputeHessian()`.

    131: .seealso: [](sec_tao_term),
    132:           `TaoTerm`,
    133:           `TaoTermComputeHessian()`,
    134:           `TaoTermGetFDDelta()`,
    135:           `TaoTermSetFDDelta()`,
    136:           `TaoTermComputeHessianSetUseFD()`,
    137:           `TaoTermComputeHessianGetUseFD()`
    138: @*/
    139: PetscErrorCode TaoTermComputeHessianFD(TaoTerm term, Vec x, Vec params, Mat H, Mat Hpre)
    140: {
    141:   SNES                  snes;
    142:   DM                    dm;
    143:   TaoTermWithParameters t;

    145:   PetscFunctionBegin;
    151:   PetscCheck(H || Hpre, PetscObjectComm((PetscObject)term), PETSC_ERR_ARG_NULL, "At least one of H or Hpre must be non-NULL");
    152:   // Note: Request for FD takes higher precedence over MATMFFD. Ignore MFFD
    153:   PetscCall(PetscInfo(term, "%s: TaoTerm using finite differences w/o coloring to compute Hessian matrix.\n", ((PetscObject)term)->prefix));
    154:   // Note: same routine as in fdiff.c
    155:   PetscCall(SNESCreate(PetscObjectComm((PetscObject)term), &snes));

    157:   t.term   = term;
    158:   t.params = params;
    159:   PetscCall(SNESSetFunction(snes, NULL, SNESFunction_TaoTerm, (void *)&t));
    160:   PetscCall(SNESGetDM(snes, &dm));
    161:   PetscCall(DMShellSetGlobalVector(dm, x));
    162:   PetscCall(SNESSetUp(snes));
    163:   if (H) {
    164:     PetscInt n, N;

    166:     PetscCall(VecGetSize(x, &N));
    167:     PetscCall(VecGetLocalSize(x, &n));
    168:     PetscCall(MatSetSizes(H, n, n, N, N));
    169:     PetscCall(MatSetUp(H));
    170:   }
    171:   if (Hpre && Hpre != H) {
    172:     PetscInt n, N;

    174:     PetscCall(VecGetSize(x, &N));
    175:     PetscCall(VecGetLocalSize(x, &n));
    176:     PetscCall(MatSetSizes(Hpre, n, n, N, N));
    177:     PetscCall(MatSetUp(Hpre));
    178:   }
    179:   PetscCall(SNESComputeJacobianDefault(snes, x, H ? H : Hpre, Hpre ? Hpre : H, NULL));
    180:   PetscCall(SNESDestroy(&snes));
    181:   PetscFunctionReturn(PETSC_SUCCESS);
    182: }

    184: static PetscErrorCode MatMFFDFunction_TaoTermHessianShell(void *ctx, Vec x, Vec g)
    185: {
    186:   TaoTermWithParameters *tp = (TaoTermWithParameters *)ctx;

    188:   PetscFunctionBegin;
    189:   // we expect the solution to move around in a finite difference method, but not the parameters
    190:   // TODO but not checking for it now
    191:   PetscCall(TaoTermComputeGradient(tp->term, x, tp->params, g));
    192:   tp->term->ngrad_mffd++;
    193:   PetscFunctionReturn(PETSC_SUCCESS);
    194: }

    196: static PetscErrorCode TaoTermInitializeHessianMFFD(TaoTerm term, Mat mffd)
    197: {
    198:   TaoTermWithParameters *tp;
    199:   PetscLayout            sol_layout;
    200:   VecType                sol_vec_type;
    201:   PetscContainer         container;

    203:   PetscFunctionBegin;
    204:   PetscCall(TaoTermGetSolutionLayout(term, &sol_layout));
    205:   PetscCall(MatSetLayouts(mffd, sol_layout, sol_layout));
    206:   PetscCall(TaoTermGetSolutionVecType(term, &sol_vec_type));
    207:   PetscCall(MatSetVecType(mffd, sol_vec_type));
    208:   PetscCall(MatSetType(mffd, MATMFFD));
    209:   PetscCall(PetscNew(&tp));
    210:   tp->term = term;
    211:   PetscCall(MatMFFDSetFunction(mffd, MatMFFDFunction_TaoTermHessianShell, tp));
    212:   PetscCall(MatSetOption(mffd, MAT_SYMMETRIC, PETSC_TRUE));
    213:   PetscCall(MatSetOption(mffd, MAT_SYMMETRY_ETERNAL, PETSC_TRUE));
    214:   PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)term), &container));
    215:   PetscCall(PetscContainerSetPointer(container, (void *)tp));
    216:   PetscCall(PetscContainerSetCtxDestroy(container, TaoTermWithParametersDestroy));
    217:   PetscCall(PetscObjectCompose((PetscObject)mffd, "__TaoTermWithParameters", (PetscObject)container));
    218:   PetscCall(PetscContainerDestroy(&container));
    219:   PetscCall(PetscFree(term->Hpre_mattype));
    220:   PetscCall(PetscFree(term->H_mattype));
    221:   PetscCall(PetscStrallocpy(MATMFFD, (char **)&term->H_mattype));
    222:   PetscCall(PetscStrallocpy(MATMFFD, (char **)&term->Hpre_mattype));
    223:   PetscFunctionReturn(PETSC_SUCCESS);
    224: }

    226: /*@
    227:   TaoTermCreateHessianMFFD - Create a `MATMFFD` for a matrix-free finite-difference approximation of the Hessian of a `TaoTerm`

    229:   Collective

    231:   Input Parameter:
    232: . term - a `TaoTerm`

    234:   Output Parameter:
    235: . mffd - a `Mat` of type `MATMFFD`

    237:   Level: advanced

    239: .seealso: [](sec_tao_term), `TaoTerm`, `TaoTermComputeHessianFD()`
    240: @*/
    241: PetscErrorCode TaoTermCreateHessianMFFD(TaoTerm term, Mat *mffd)
    242: {
    243:   PetscFunctionBegin;
    245:   PetscAssertPointer(mffd, 2);
    246:   PetscCall(MatCreate(PetscObjectComm((PetscObject)term), mffd));
    247:   PetscCall(TaoTermInitializeHessianMFFD(term, *mffd));
    248:   PetscFunctionReturn(PETSC_SUCCESS);
    249: }

    251: PetscErrorCode TaoTermComputeHessianMFFD(TaoTerm term, Vec x, Vec params, Mat H, Mat B)
    252: {
    253:   PetscContainer         container;
    254:   TaoTermWithParameters *tp;

    256:   PetscFunctionBegin;
    262:   PetscCall(PetscObjectQuery((PetscObject)H, "__TaoTermWithParameters", (PetscObject *)&container));
    263:   if (!container) {
    264:     PetscCall(TaoTermInitializeHessianMFFD(term, H));
    265:     PetscCall(PetscObjectQuery((PetscObject)H, "__TaoTermWithParameters", (PetscObject *)&container));
    266:     PetscCheck(container, PetscObjectComm((PetscObject)term), PETSC_ERR_PLIB, "failed to initialize mffd matrix");
    267:   }
    268:   PetscCall(PetscContainerGetPointer(container, (void **)&tp));
    269:   PetscCheck(tp->term == term, PetscObjectComm((PetscObject)term), PETSC_ERR_ARG_INCOMP, "Hessian shell matrix does not come from this TaoTerm");
    270:   PetscCall(PetscObjectReference((PetscObject)params));
    271:   PetscCall(VecDestroy(&tp->params));
    272:   tp->params = params;
    273:   PetscCall(MatMFFDSetBase(H, x, NULL));
    274:   PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
    275:   PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
    276:   PetscFunctionReturn(PETSC_SUCCESS);
    277: }