Actual source code: cutest.c
1: static char help[] = "Solve a decoded unconstrained CUTEst problem with TAO.\n\
2: -cutest_lib filename Shared library containing the decoded problem\n\
3: -cutest_data filename Decoded data file (default OUTSDIF.d)\n\
4: -cutest_hessian (sparse|shell) Sparse Hessian or exact Hessian-vector products\n\
5: -cutest_solution_view View the final solution\n\
6: -cutest_result filename Write solve statistics as CSV instead of a screen summary\n\
7: See the TAO users manual for configuration, decoding, and running examples.\n";
9: #include <petsc/private/taocutestimpl.h>
11: typedef struct {
12: PetscErrorCode (*mult)(Mat, Vec, Vec);
13: PetscErrorCode (*multtranspose)(Mat, Vec, Vec);
14: PetscCount products;
15: } HessianCounter;
17: static PetscErrorCode CountHessianMult(Mat H, Vec X, Vec Y)
18: {
19: HessianCounter *counter;
21: PetscFunctionBeginUser;
22: PetscCall(PetscObjectContainerQuery((PetscObject)H, "CUTEstHessianProducts", &counter));
23: PetscCall((*counter->mult)(H, X, Y));
24: ++counter->products;
25: PetscFunctionReturn(PETSC_SUCCESS);
26: }
28: static PetscErrorCode CountHessianMultTranspose(Mat H, Vec X, Vec Y)
29: {
30: HessianCounter *counter;
32: PetscFunctionBeginUser;
33: PetscCall(PetscObjectContainerQuery((PetscObject)H, "CUTEstHessianProducts", &counter));
34: PetscCall((*counter->multtranspose)(H, X, Y));
35: ++counter->products;
36: PetscFunctionReturn(PETSC_SUCCESS);
37: }
39: static PetscErrorCode WriteResult(Tao tao, Vec X, CUTEstCtx *user, const char filename[], PetscReal initial_f, PetscReal initial_gnorm, const PetscReal initial_calls[], PetscLogDouble elapsed, PetscCount products)
40: {
41: PetscViewer viewer;
42: TaoConvergedReason reason;
43: SNESConvergedReason snes_reason = SNES_CONVERGED_ITERATING;
44: TaoType type;
45: const char *solver, *snes_reason_name = "";
46: char name[256];
47: Vec G;
48: PetscReal f, gnorm, calls[4], time[4];
49: PetscInt its, rejected = -1;
50: PetscBool is_snes, is_python;
52: PetscFunctionBeginUser;
53: /* Exclude the final verification evaluation from CUTEst's solve counters. */
54: PetscCallCUTEst(CUTEST_ureport, calls, time);
55: PetscCheck(!user->hessian_x || (PetscReal)products == calls[3] - initial_calls[3], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Hessian multiplication count disagrees with CUTEst");
56: PetscCall(VecDuplicate(X, &G));
57: PetscCall(CUTEstFormObjectiveGradient(tao, X, &f, G, user));
58: PetscCall(VecNorm(G, NORM_2, &gnorm));
59: PetscCall(VecDestroy(&G));
60: PetscCall(TaoGetSolutionStatus(tao, &its, NULL, NULL, NULL, NULL, &reason));
61: PetscCall(TaoGetType(tao, &type));
62: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOSNES, &is_snes));
63: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOPYTHON, &is_python));
64: solver = name;
65: if (is_snes) {
66: SNES snes;
67: SNESType snes_type;
68: PetscBool newton;
70: PetscCall(TaoSNESGetSNES(tao, &snes));
71: PetscCall(SNESGetType(snes, &snes_type));
72: PetscCall(SNESGetConvergedReason(snes, &snes_reason));
73: PetscCall(SNESGetNonlinearStepFailures(snes, &rejected));
74: snes_reason_name = SNESConvergedReasons[snes_reason];
75: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESPYTHON, &is_python));
76: if (is_python) PetscCall(SNESPythonGetType(snes, &solver));
77: else {
78: PetscCall(PetscStrncmp(snes_type, "newton", 6, &newton));
79: PetscCall(PetscSNPrintf(name, sizeof(name), "snes%s", newton ? snes_type + 6 : snes_type));
80: }
81: } else if (is_python) PetscCall(TaoPythonGetType(tao, &solver));
82: else PetscCall(PetscSNPrintf(name, sizeof(name), "tao%s", type));
83: if (filename[0]) {
84: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_SELF, filename, &viewer));
85: PetscCall(PetscViewerASCIIPrintf(viewer, "n,solver,reason_code,reason,snes_reason_code,snes_reason,iterations,initial_objective,initial_gradient_norm,objective,gradient_norm,objective_evaluations,gradient_evaluations,hessian_"
86: "evaluations,hessian_products,solve_seconds,rejected_steps,cutest_hessian_products\n"));
87: PetscCall(PetscViewerASCIIPrintf(viewer, "%d,%s,%d,%s,%d,%s,%" PetscInt_FMT ",%.17g,%.17g,%.17g,%.17g,%.0f,%.0f,%.0f,%" PetscCount_FMT ",%.17g,%" PetscInt_FMT ",%.0f\n", user->n, solver, (int)reason, TaoConvergedReasons[reason], (int)snes_reason, snes_reason_name, its, (double)initial_f, (double)initial_gnorm, (double)f, (double)gnorm, (double)(calls[0] - initial_calls[0]), (double)(calls[1] - initial_calls[1]), (double)(calls[2] - initial_calls[2]), products, (double)elapsed, rejected, (double)(calls[3] - initial_calls[3])));
88: PetscCall(PetscViewerDestroy(&viewer));
89: } else {
90: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\nCUTEst solve summary\n"));
91: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Solver: %s\n", solver));
92: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Variables: %d\n", user->n));
93: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Termination reason: %s\n", TaoConvergedReasons[reason]));
94: if (is_snes) PetscCall(PetscPrintf(PETSC_COMM_SELF, " SNES termination reason: %s\n", snes_reason_name));
95: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Iterations: %" PetscInt_FMT "\n", its));
96: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n %-26s %-16s %s\n", "", "Initial", "Final"));
97: PetscCall(PetscPrintf(PETSC_COMM_SELF, " %-26s %-16.8e %.8e\n", "Objective:", (double)initial_f, (double)f));
98: PetscCall(PetscPrintf(PETSC_COMM_SELF, " %-26s %-16.8e %.8e\n", "Gradient norm (2-norm):", (double)initial_gnorm, (double)gnorm));
99: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n Objective evaluations: %.0f\n", (double)(calls[0] - initial_calls[0])));
100: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Gradient evaluations: %.0f\n", (double)(calls[1] - initial_calls[1])));
101: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Hessian evaluations: %.0f\n", (double)(calls[2] - initial_calls[2])));
102: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Hessian-vector products: %" PetscCount_FMT "\n", products));
103: PetscCall(PetscPrintf(PETSC_COMM_SELF, " CUTEst Hessian products: %.0f\n", (double)(calls[3] - initial_calls[3])));
104: if (rejected >= 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Rejected steps: %" PetscInt_FMT "\n", rejected));
105: else PetscCall(PetscPrintf(PETSC_COMM_SELF, " Rejected steps: unavailable\n"));
106: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Solve time (seconds): %.6g\n", (double)elapsed));
107: }
108: PetscFunctionReturn(PETSC_SUCCESS);
109: }
111: int main(int argc, char **argv)
112: {
113: const char *hessians[] = {"sparse", "shell"};
114: char library[PETSC_MAX_PATH_LEN] = "", data[PETSC_MAX_PATH_LEN] = "OUTSDIF.d", result[PETSC_MAX_PATH_LEN] = "";
115: PetscDLLibrary dll = NULL;
116: CUTEstCtx user = {0};
117: HessianCounter counter = {0};
118: Tao tao;
119: Vec X;
120: Mat H;
121: PetscInt hessian = 0;
122: PetscMPIInt size;
123: PetscReal initial_f = 0.0, initial_gnorm = 0.0, initial_calls[4] = {0.0};
124: PetscLogDouble start, end;
126: PetscFunctionBeginUser;
127: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
128: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
129: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This CUTEst driver runs on one MPI process");
130: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "CUTEst problem options", "Tao");
131: PetscCall(PetscOptionsString("-cutest_lib", "Decoded problem library", NULL, library, library, sizeof(library), NULL));
132: PetscCall(PetscOptionsString("-cutest_data", "Decoded problem data", NULL, data, data, sizeof(data), NULL));
133: PetscCall(PetscOptionsEList("-cutest_hessian", "Hessian representation", NULL, hessians, 2, hessians[hessian], &hessian, NULL));
134: PetscCall(PetscOptionsString("-cutest_result", "CSV solve statistics", NULL, result, result, sizeof(result), NULL));
135: PetscOptionsEnd();
136: PetscCall(CUTEstLoadProblem(library, data, &user, &dll, &X));
138: PetscCall(CUTEstCreateHessian(X, (PetscBool)hessian, &user, &H));
139: PetscCall(TaoCreate(PETSC_COMM_SELF, &tao));
140: PetscCall(TaoSetType(tao, TAOLMVM));
141: PetscCall(TaoSetSolution(tao, X));
142: PetscCall(TaoSetObjective(tao, CUTEstFormObjective, &user));
143: PetscCall(TaoSetGradient(tao, NULL, CUTEstFormGradient, &user));
144: PetscCall(TaoSetObjectiveAndGradient(tao, NULL, CUTEstFormObjectiveGradient, &user));
145: PetscCall(PetscObjectContainerCompose((PetscObject)H, "CUTEstHessianProducts", &counter, NULL));
146: PetscCall(MatGetOperation(H, MATOP_MULT, (PetscErrorCodeFn **)&counter.mult));
147: PetscCall(MatSetOperation(H, MATOP_MULT, (PetscErrorCodeFn *)CountHessianMult));
148: PetscCall(MatGetOperation(H, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn **)&counter.multtranspose));
149: PetscCall(MatSetOperation(H, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)CountHessianMultTranspose));
150: PetscCall(TaoSetHessian(tao, H, H, CUTEstFormHessian, &user));
151: PetscCall(TaoSetFromOptions(tao));
152: {
153: Vec G;
154: PetscReal time[4];
156: PetscCall(VecDuplicate(X, &G));
157: PetscCall(CUTEstFormObjectiveGradient(tao, X, &initial_f, G, &user));
158: PetscCall(VecNorm(G, NORM_2, &initial_gnorm));
159: PetscCall(VecDestroy(&G));
160: PetscCallCUTEst(CUTEST_ureport, initial_calls, time);
161: }
162: counter.products = 0;
163: PetscCall(PetscTime(&start));
164: PetscCall(TaoSolve(tao));
165: PetscCall(PetscTime(&end));
166: PetscCall(WriteResult(tao, X, &user, result, initial_f, initial_gnorm, initial_calls, end - start, counter.products));
167: PetscCall(VecViewFromOptions(X, NULL, "-cutest_solution_view"));
168: PetscCall(TaoDestroy(&tao));
169: PetscCall(MatDestroy(&H));
170: PetscCall(VecDestroy(&user.hessian_x));
171: PetscCall(PetscFree3(user.rows, user.cols, user.values));
172: PetscCall(VecDestroy(&X));
173: PetscCall(CUTEstUnloadProblem(dll));
174: PetscCall(PetscFinalize());
175: return 0;
176: }
178: /*TEST
180: build:
181: requires: cutest double !complex defined(PETSC_HAVE_DYNAMIC_LIBRARIES)
183: testset:
184: requires: sifdecode
185: command: @PYTHON@ ${wPETSC_DIR}/share/petsc/cutest/run.py --petsc-dir ${wPETSC_DIR} --petsc-arch=${petsc_arch} --driver ../cutest --sif ${wPETSC_DIR}/share/petsc/datafiles/cutest/ROSENBR.SIF -- ${args} @SUBARGS@
186: args: -tao_gatol 1.e-6 -tao_grtol 0 -tao_gttol 0 -tao_converged_reason -cutest_result result.csv
187: temporaries: result.csv
188: filter: sed -E "s/ iterations [0-9]+//"
189: output_file: output/cutest_1.out
190: test:
191: suffix: 1
192: args: -tao_type {{lmvm nls ntr}} -cutest_hessian {{sparse shell}}
193: test:
194: suffix: snes
195: args: -tao_type snes -snes_type newtontr -snes_atol 1.e-6 -ksp_type cg -pc_type none -cutest_hessian {{sparse shell}}
196: output_file: output/cutest_snes.out
198: testset:
199: requires: sifdecode
200: command: @PYTHON@ ${wPETSC_DIR}/share/petsc/cutest/run.py --petsc-dir ${wPETSC_DIR} --petsc-arch=${petsc_arch} --driver ../cutest --sif ${wPETSC_DIR}/share/petsc/datafiles/cutest/ROSENBR.SIF -- ${args} @SUBARGS@
201: args: -tao_gatol 1.e-6 -tao_grtol 0 -tao_gttol 0
202: filter: sed -E "s,Solve time [(]seconds[)]:.*,Solve time (seconds):,"
203: test:
204: suffix: summary_ntr
205: args: -tao_type ntr
206: test:
207: suffix: summary_snes_sparse
208: args: -tao_type snes -snes_type newtonls -snes_atol 1.e-6 -ksp_type bicg -pc_type none
209: test:
210: suffix: summary_snes_shell
211: args: -tao_type snes -snes_type newtonls -snes_atol 1.e-6 -ksp_type bicg -pc_type none -cutest_hessian shell
213: TEST*/