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*/