Actual source code: rosenbrock1_taoterm.c

      1: /*  Include "petsctao.h" so we can use TAO solvers.  */
      2: #include <petsctao.h>
      3: #include "rosenbrock1.h" // defines AppCtx, AppCtxFormFunctionGradient(), and AppCtxFormHessian()

      5: static char help[] = "This example demonstrates use of the TaoTerm\n\
      6: interface for defining problems in the Tao library.  This example\n\
      7: should be compared to rosenbrock1.c, which uses the callback interface\n\
      8: to define the Rosenbrock function.\n";

     10: static PetscErrorCode FormFunctionGradient(TaoTerm, Vec, Vec, PetscReal *, Vec);
     11: static PetscErrorCode FormHessian(TaoTerm, Vec, Vec, Mat, Mat);
     12: static PetscErrorCode CreateSolutionVec(TaoTerm, Vec *);

     14: static PetscErrorCode CtxDestroy(PetscCtxRt ctx)
     15: {
     16:   PetscFunctionBeginUser;
     17:   PetscFunctionReturn(PETSC_SUCCESS);
     18: }

     20: int main(int argc, char **argv)
     21: {
     22:   TaoTerm     objective;
     23:   Tao         tao;  /* Tao solver context */
     24:   PetscMPIInt size; /* number of processes running */
     25:   AppCtx      user; /* user-defined application context */
     26:   MPI_Comm    comm;
     27:   PetscBool   test_gradient_fd_check = PETSC_FALSE; /* test that FD delta is preserved */
     28:   PetscReal   fd_delta_set           = 1.e-6;

     30:   /* Initialize TAO and PETSc */
     31:   PetscFunctionBeginUser;
     32:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
     33:   comm = PETSC_COMM_WORLD;
     34:   PetscCallMPI(MPI_Comm_size(comm, &size));
     35:   PetscCheck(size == 1, comm, PETSC_ERR_WRONG_MPI_SIZE, "Incorrect number of processors");

     37:   PetscOptionsBegin(comm, "", help, "none");
     38:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_gradient_fd_check", &test_gradient_fd_check, NULL));
     39:   PetscOptionsEnd();
     40:   /* Initialize problem parameters */
     41:   PetscCall(AppCtxInitialize(comm, &user));

     43:   /* Define the objective function */
     44:   PetscCall(TaoTermCreateShell(comm, &user, CtxDestroy, &objective));
     45:   PetscCall(TaoTermSetParametersMode(objective, TAOTERM_PARAMETERS_NONE));
     46:   PetscCall(TaoTermShellSetCreateSolutionVec(objective, CreateSolutionVec));
     47:   PetscCall(TaoTermShellSetObjectiveAndGradient(objective, FormFunctionGradient));
     48:   PetscCall(TaoTermShellSetCreateHessianMatrices(objective, TaoTermCreateHessianMatricesDefault));
     49:   if (user.jacobi_pc) PetscCall(TaoTermSetCreateHessianMode(objective, PETSC_FALSE, MATBAIJ, MATBAIJ));
     50:   else PetscCall(TaoTermSetCreateHessianMode(objective, PETSC_TRUE /* H == Hpre */, MATBAIJ, NULL));
     51:   PetscCall(TaoTermShellSetHessian(objective, FormHessian));

     53:   /* Create TAO solver with desired solution method */
     54:   PetscCall(TaoCreate(PETSC_COMM_SELF, &tao));
     55:   PetscCall(TaoSetType(tao, TAOLMVM));

     57:   if (user.use_fd) {
     58:     PetscCall(TaoTermSetFDDelta(objective, 7.e-9));
     59:     PetscCall(TaoTermComputeGradientSetUseFD(objective, PETSC_TRUE));
     60:     PetscCall(TaoTermComputeHessianSetUseFD(objective, PETSC_TRUE));
     61:   }
     62:   /* Set routines for function, gradient, hessian evaluation */
     63:   PetscCall(TaoAddTerm(tao, NULL, 1.0, objective, NULL, NULL));

     65:   /* Check for TAO command line options */
     66:   PetscCall(TaoSetFromOptions(tao));

     68:   /* Set FD delta for testing if requested (after options processing) */
     69:   if (test_gradient_fd_check) PetscCall(TaoTermSetFDDelta(objective, fd_delta_set));

     71:   /* SOLVE THE APPLICATION */
     72:   PetscCall(TaoSolve(tao));

     74:   /* Check that FD delta is preserved if testing */
     75:   if (test_gradient_fd_check) {
     76:     PetscReal fd_delta_get;

     78:     PetscCall(TaoTermGetFDDelta(objective, &fd_delta_get));
     79:     PetscCheck(PetscAbsReal(fd_delta_get - 1.e-6) < 1.e-15, comm, PETSC_ERR_PLIB, "FD delta changed: set %g, got %g", (double)fd_delta_set, (double)fd_delta_get);
     80:   }

     82:   /* Clean up */
     83:   PetscCall(AppCtxFinalize(&user, tao));
     84:   PetscCall(TaoDestroy(&tao));
     85:   PetscCall(TaoTermDestroy(&objective));

     87:   PetscCall(PetscFinalize());
     88:   return 0;
     89: }

     91: /*
     92:   FormFunctionGradient - Evaluates the function, f(X), and gradient, G(X).

     94:   Input Parameters:
     95: + term              - the `TaoTerm` for the objective function
     96: . X                 - input vector
     97: - parameters_unused - optional vector of parameters that this rosenbrock function does not use

     99:   Output Parameters:
    100: + f - function value
    101: - G - vector containing the newly evaluated gradient

    103:   Note:
    104:   Some optimization methods ask for the function and the gradient evaluation
    105:   at the same time.  Evaluating both at once may be more efficient that
    106:   evaluating each separately.
    107: */
    108: static PetscErrorCode FormFunctionGradient(TaoTerm term, Vec X, Vec parameters_unused, PetscReal *f, Vec G)
    109: {
    110:   AppCtx *user;

    112:   PetscFunctionBeginUser;
    113:   PetscCheck(parameters_unused == NULL, PetscObjectComm((PetscObject)term), PETSC_ERR_ARG_WRONG, "Rosenbrock function does not take a parameter vector");
    114:   PetscCall(TaoTermShellGetContext(term, &user));
    115:   PetscCall(AppCtxFormFunctionGradient(user, X, f, G));
    116:   PetscFunctionReturn(PETSC_SUCCESS);
    117: }

    119: /*
    120:   FormHessian - Evaluates Hessian matrix.

    122:   Input Parameters:
    123: + tao    - the Tao context
    124: . x      - input vector
    125: . params - optional vector of parameters that this rosenbrock function does not use
    126: - ptr    - optional user-defined context, as set by TaoSetHessian()

    128:   Output Parameters:
    129: + H    - Hessian matrix
    130: - Hpre - Preconditioning matrix

    132:   Note:  Providing the Hessian may not be necessary.  Only some solvers
    133:   require this matrix.
    134: */
    135: static PetscErrorCode FormHessian(TaoTerm term, Vec X, Vec params, Mat H, Mat Hpre)
    136: {
    137:   AppCtx *user;

    139:   PetscFunctionBeginUser;
    140:   PetscCheck(params == NULL, PetscObjectComm((PetscObject)term), PETSC_ERR_ARG_WRONG, "Rosenbrock function does not take a parameter vector");
    141:   PetscCall(TaoTermShellGetContext(term, &user));
    142:   if (H) PetscCall(AppCtxFormHessian(user, X, H));
    143:   if (Hpre && Hpre != H) {
    144:     if (user->jacobi_pc) {
    145:       Vec v;

    147:       PetscCall(VecDuplicate(X, &v));
    148:       PetscCall(MatGetDiagonal(H, v));
    149:       PetscCall(MatZeroEntries(Hpre));
    150:       PetscCall(MatDiagonalSet(Hpre, v, INSERT_VALUES));
    151:       PetscCall(VecDestroy(&v));
    152:     } else {
    153:       if (H) PetscCall(MatCopy(H, Hpre, SAME_NONZERO_PATTERN));
    154:       else PetscCall(AppCtxFormHessian(user, X, Hpre));
    155:     }
    156:   }
    157:   PetscFunctionReturn(PETSC_SUCCESS);
    158: }

    160: static PetscErrorCode CreateSolutionVec(TaoTerm term, Vec *solution)
    161: {
    162:   AppCtx *user;

    164:   PetscFunctionBeginUser;
    165:   PetscCall(TaoTermShellGetContext(term, &user));
    166:   PetscCall(AppCtxCreateSolution(user, solution));
    167:   PetscFunctionReturn(PETSC_SUCCESS);
    168: }

    170: /*TEST

    172:    build:
    173:      requires: !complex !single !quad !defined(PETSC_USE_64BIT_INDICES) !__float128

    175:    test:
    176:      output_file: output/rosenbrock1_1.out
    177:      args: -tao_monitor_short -tao_type nls -tao_gatol 1.e-4

    179:    test:
    180:      suffix: test_gradient
    181:      args: -tao_monitor_short -tao_type nls -tao_gatol 1.e-4 -tao_test_gradient -tao_fd_delta 1.e-6 -n 4 -chained -tao_term_hessian_mat_type aij -alpha 49.0

    183:    test:
    184:      suffix: test_gradient_fd_check
    185:      args: -tao_monitor_short -tao_type nls -tao_gatol 1.e-4 -tao_test_gradient
    186:      args: -n 4 -chained -tao_term_hessian_mat_type aij -alpha 49.0 -test_gradient_fd_check 1

    188:    test:
    189:      suffix: fd_grad
    190:      args: -tao_monitor_short -tao_type nls -tao_term_gradient_use_fd

    192:    test:
    193:      suffix: fd_hess
    194:      args: -tao_monitor_short -tao_type nls -tao_term_hessian_use_fd

    196:    test:
    197:      suffix: fd_hess_diffpre
    198:      args: -tao_monitor_short -tao_type nls -tao_term_hessian_use_fd -tao_term_hessian_pre_is_hessian 0

    200:    test:
    201:      suffix: use_fd
    202:      args: -tao_monitor_short -tao_type nls -use_fd

    204:    test:
    205:      suffix: test_fd_hess
    206:      args: -tao_type nls -tao_fd_hessian -tao_view -tao_monitor_short

    208:    test:
    209:      suffix: test_mf_hessian
    210:      args: -tao_type nls -tao_term_hessian_mat_type mffd -tao_monitor_short -tao_view

    212:    test:
    213:      suffix: add_term
    214:      args: -tao_type nls -tao_add_terms extra_ -extra_tao_term_type halfl2squared -tao_view ::ascii_info_detail

    216:    test:
    217:      suffix: separate_hessians
    218:      requires: !single
    219:      args: -tao_monitor_short -tao_type nls -tao_gatol 1.e-4 -tao_term_hessian_pre_is_hessian 0 -tao_term_hessian_mat_type aij -tao_term_hessian_pre_mat_type sbaij -tao_view

    221: TEST*/