Actual source code: cutest.c
1: static char help[] = "Compare assembled and sampled shell CUTEst Hessians.\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_check_type (matrix|mult) Compare full matrices or forward and transpose products (default matrix)\n\
5: -cutest_check_vectors n Number of products for each multiplication check (default 8)\n";
7: #include <petsc/private/taocutestimpl.h>
9: static PetscErrorCode CheckHessian(Mat H, Mat reference, PetscBool mult, PetscInt nvec, PetscInt point, const char stage[])
10: {
11: Mat computed, difference;
12: PetscReal norm, error;
13: PetscBool match;
15: PetscFunctionBeginUser;
16: for (PetscInt transpose = 0; transpose < 2; ++transpose) {
17: if (mult) {
18: if (transpose) PetscCall(MatMultTransposeEqual(reference, H, nvec, &match));
19: else PetscCall(MatMultEqual(reference, H, nvec, &match));
20: PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Hessian%s mismatch at point %" PetscInt_FMT " (%s); rerun with -cutest_check_type matrix to display the matrices", transpose ? " transpose" : "", point, stage);
21: continue;
22: }
23: if (transpose) PetscCall(MatComputeOperatorTranspose(H, MATDENSE, &computed));
24: else PetscCall(MatComputeOperator(H, MATDENSE, &computed));
25: PetscCall(MatDuplicate(computed, MAT_COPY_VALUES, &difference));
26: PetscCall(MatNorm(computed, NORM_FROBENIUS, &norm));
27: PetscCall(MatAXPY(difference, -1.0, reference, DIFFERENT_NONZERO_PATTERN));
28: PetscCall(MatNorm(difference, NORM_FROBENIUS, &error));
29: match = (PetscBool)(!PetscIsInfOrNanReal(norm) && !PetscIsInfOrNanReal(error) && error <= 100.0 * PETSC_MACHINE_EPSILON * PetscMax(norm, 1.0));
30: if (!match) {
31: PetscViewer viewer = PETSC_VIEWER_STDOUT_SELF;
33: PetscCall(PetscViewerASCIIPrintf(viewer, "Hessian%s mismatch at point %" PetscInt_FMT " (%s)\nComputed Hessian:\n", transpose ? " transpose" : "", point, stage));
34: PetscCall(MatView(computed, viewer));
35: PetscCall(PetscViewerASCIIPrintf(viewer, "Reference Hessian:\n"));
36: PetscCall(MatView(reference, viewer));
37: PetscCall(PetscViewerASCIIPrintf(viewer, "Difference (computed - reference):\n"));
38: PetscCall(MatView(difference, viewer));
39: PetscCall(PetscViewerFlush(viewer));
40: }
41: PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Hessian%s mismatch at point %" PetscInt_FMT " (%s): error %g, computed norm %g", transpose ? " transpose" : "", point, stage, (double)error, (double)norm);
42: PetscCall(PetscInfo(H, "Point %" PetscInt_FMT " (%s): relative Hessian%s difference %g\n", point, stage, transpose ? " transpose" : "", (double)(error / PetscMax(norm, 1.0))));
43: PetscCall(MatDestroy(&computed));
44: PetscCall(MatDestroy(&difference));
45: }
46: PetscFunctionReturn(PETSC_SUCCESS);
47: }
49: int main(int argc, char **argv)
50: {
51: const char *comparisons[] = {"matrix", "mult"};
52: char library[PETSC_MAX_PATH_LEN] = "", data[PETSC_MAX_PATH_LEN] = "OUTSDIF.d";
53: PetscDLLibrary dll = NULL;
54: CUTEstCtx assembled = {0}, shell = {0};
55: Vec X, initial, trial, G;
56: Mat H, S;
57: PetscInt comparison = 0, nvec = 8;
58: PetscMPIInt size;
60: PetscFunctionBeginUser;
61: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
62: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
63: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This CUTEst test runs on one MPI process");
64: PetscCall(PetscOptionsGetString(NULL, NULL, "-cutest_lib", library, sizeof(library), NULL));
65: PetscCall(PetscOptionsGetString(NULL, NULL, "-cutest_data", data, sizeof(data), NULL));
66: PetscOptionsBegin(PETSC_COMM_SELF, NULL, "CUTEst Hessian comparison", NULL);
67: PetscCall(PetscOptionsEList("-cutest_check_type", "Hessian comparison method", NULL, comparisons, 2, comparisons[comparison], &comparison, NULL));
68: PetscCall(PetscOptionsInt("-cutest_check_vectors", "Number of random products", "MatMultEqual", nvec, &nvec, NULL));
69: PetscOptionsEnd();
70: PetscCheck(nvec > 0, PETSC_COMM_SELF, PETSC_ERR_USER_INPUT, "The number of random products must be positive");
71: PetscCall(CUTEstLoadProblem(library, data, &assembled, &dll, &X));
72: shell.n = assembled.n;
73: PetscCall(CUTEstCreateHessian(X, PETSC_FALSE, &assembled, &H));
74: PetscCall(CUTEstCreateHessian(X, PETSC_TRUE, &shell, &S));
75: PetscCall(VecDuplicate(X, &initial));
76: PetscCall(VecCopy(X, initial));
77: PetscCall(VecDuplicate(X, &trial));
78: PetscCall(VecDuplicate(X, &G));
79: for (PetscInt point = 0; point < 5; ++point) {
80: Mat reference;
81: PetscScalar *x;
82: PetscReal f;
84: if (!point || point == 4) PetscCall(VecCopy(initial, X));
85: else {
86: PetscCall(VecGetArrayWrite(X, &x));
87: for (PetscInt i = 0; i < assembled.n; ++i) x[i] = point == 1 ? 0.0 : (point == 2 ? 0.125 : -0.25) * (i + 1);
88: PetscCall(VecRestoreArrayWrite(X, &x));
89: }
90: PetscCall(CUTEstFormHessian(NULL, X, H, H, &assembled));
91: PetscCall(CUTEstFormHessian(NULL, X, S, S, &shell));
92: /* Keep an immutable reference across trial evaluations and solver shifts. */
93: if (comparison) PetscCall(MatDuplicate(H, MAT_COPY_VALUES, &reference));
94: else PetscCall(MatComputeOperator(S, MATDENSE, &reference));
95: PetscCall(CheckHessian(H, reference, (PetscBool)comparison, nvec, point, "assembled"));
96: PetscCall(CheckHessian(S, reference, (PetscBool)comparison, nvec, point, "cached shell"));
98: PetscCall(VecCopy(X, trial));
99: PetscCall(VecShift(trial, 0.375));
100: PetscCall(CUTEstFormObjective(NULL, trial, &f, &shell));
101: PetscCall(CheckHessian(S, reference, (PetscBool)comparison, nvec, point, "objective trial"));
102: PetscCall(CUTEstFormGradient(NULL, trial, G, &shell));
103: PetscCall(CheckHessian(S, reference, (PetscBool)comparison, nvec, point, "gradient trial"));
104: PetscCall(CUTEstFormObjectiveGradient(NULL, trial, &f, G, &shell));
105: PetscCall(CheckHessian(S, reference, (PetscBool)comparison, nvec, point, "objective/gradient trial"));
107: PetscCall(MatShift(H, 1.0));
108: PetscCall(MatShift(S, 2.0));
109: PetscCall(CUTEstFormHessian(NULL, X, H, H, &assembled));
110: PetscCall(CUTEstFormHessian(NULL, X, S, S, &shell));
111: PetscCall(CheckHessian(H, reference, (PetscBool)comparison, nvec, point, "assembled after shift"));
112: PetscCall(CheckHessian(S, reference, (PetscBool)comparison, nvec, point, "shell after shift"));
113: PetscCall(MatDestroy(&reference));
114: }
115: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Assembled and shell Hessians and their transposes agree at all 5 points\n"));
116: PetscCall(MatDestroy(&H));
117: PetscCall(MatDestroy(&S));
118: PetscCall(VecDestroy(&shell.hessian_x));
119: PetscCall(PetscFree3(assembled.rows, assembled.cols, assembled.values));
120: PetscCall(VecDestroy(&initial));
121: PetscCall(VecDestroy(&trial));
122: PetscCall(VecDestroy(&G));
123: PetscCall(VecDestroy(&X));
124: PetscCall(CUTEstUnloadProblem(dll));
125: PetscCall(PetscFinalize());
126: return 0;
127: }
129: /*TEST
131: build:
132: requires: cutest double !complex defined(PETSC_HAVE_DYNAMIC_LIBRARIES)
134: testset:
135: requires: sifdecode
136: args: -cutest_check_type {{matrix mult}}
137: 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/${CUTEST_SIF}.SIF -- ${args} @SUBARGS@
138: test:
139: suffix: rosenbrock
140: output_file: output/cutest_rosenbrock.out
141: env: CUTEST_SIF=ROSENBR
142: test:
143: suffix: polynomial
144: output_file: output/cutest_polynomial.out
145: env: CUTEST_SIF=COOHESS
147: TEST*/