Actual source code: taocutestimpl.h
1: #pragma once
3: #include <petsctao.h>
4: EXTERN_C_BEGIN
5: #include <cutest.h>
6: EXTERN_C_END
8: typedef struct {
9: integer n, lh;
10: integer *rows, *cols;
11: PetscScalar *values;
12: Vec hessian_x;
13: logical goth;
14: } CUTEstCtx;
16: #define PetscCallCUTEst(routine, ...) \
17: do { \
18: integer status; \
19: routine(&status, __VA_ARGS__); \
20: PetscCheck(!status, PETSC_COMM_SELF, PETSC_ERR_LIB, "%s() returned CUTEst status %d", #routine, status); \
21: } while (0)
23: static PetscErrorCode CUTEstFormObjective(Tao tao, Vec X, PetscReal *f, PetscCtx ctx)
24: {
25: CUTEstCtx *user = (CUTEstCtx *)ctx;
26: const PetscScalar *x;
28: PetscFunctionBeginUser;
29: user->goth = false;
30: PetscCall(VecGetArrayRead(X, &x));
31: PetscCallCUTEst(CUTEST_ufn, &user->n, x, f);
32: PetscCall(VecRestoreArrayRead(X, &x));
33: PetscFunctionReturn(PETSC_SUCCESS);
34: }
36: static PetscErrorCode CUTEstFormGradient(Tao tao, Vec X, Vec G, PetscCtx ctx)
37: {
38: CUTEstCtx *user = (CUTEstCtx *)ctx;
39: const PetscScalar *x;
40: PetscScalar *g;
42: PetscFunctionBeginUser;
43: user->goth = false;
44: PetscCall(VecGetArrayRead(X, &x));
45: PetscCall(VecGetArrayWrite(G, &g));
46: PetscCallCUTEst(CUTEST_ugr, &user->n, x, g);
47: PetscCall(VecRestoreArrayWrite(G, &g));
48: PetscCall(VecRestoreArrayRead(X, &x));
49: PetscFunctionReturn(PETSC_SUCCESS);
50: }
52: static PetscErrorCode CUTEstFormObjectiveGradient(Tao tao, Vec X, PetscReal *f, Vec G, PetscCtx ctx)
53: {
54: CUTEstCtx *user = (CUTEstCtx *)ctx;
55: const PetscScalar *x;
56: PetscScalar *g;
57: logical grad = true;
59: PetscFunctionBeginUser;
60: user->goth = false;
61: PetscCall(VecGetArrayRead(X, &x));
62: PetscCall(VecGetArrayWrite(G, &g));
63: PetscCallCUTEst(CUTEST_uofg, &user->n, x, f, g, &grad);
64: PetscCall(VecRestoreArrayWrite(G, &g));
65: PetscCall(VecRestoreArrayRead(X, &x));
66: PetscFunctionReturn(PETSC_SUCCESS);
67: }
69: static PetscErrorCode CUTEstHessianMult(Mat H, Vec V, Vec W)
70: {
71: CUTEstCtx *user;
72: const PetscScalar *x, *v;
73: PetscScalar *w;
75: PetscFunctionBeginUser;
76: PetscCall(MatShellGetContext(H, &user));
77: PetscCall(VecGetArrayRead(user->hessian_x, &x));
78: PetscCall(VecGetArrayRead(V, &v));
79: PetscCall(VecGetArrayWrite(W, &w));
80: PetscCallCUTEst(CUTEST_uhprod, &user->n, &user->goth, x, v, w);
81: user->goth = true;
82: PetscCall(VecRestoreArrayWrite(W, &w));
83: PetscCall(VecRestoreArrayRead(V, &v));
84: PetscCall(VecRestoreArrayRead(user->hessian_x, &x));
85: PetscFunctionReturn(PETSC_SUCCESS);
86: }
88: static PetscErrorCode CUTEstFormHessian(Tao tao, Vec X, Mat H, Mat Hpre, PetscCtx ctx)
89: {
90: CUTEstCtx *user = (CUTEstCtx *)ctx;
92: PetscFunctionBeginUser;
93: user->goth = false;
94: if (user->hessian_x) { /* MATSHELL */
95: PetscCall(VecCopy(X, user->hessian_x));
96: PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
97: PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
98: } else {
99: const PetscScalar *x;
100: PetscCount nz;
101: integer nnz;
103: PetscCall(VecGetArrayRead(X, &x));
104: PetscCallCUTEst(CUTEST_ush, &user->n, x, &nnz, &user->lh, user->values, user->rows, user->cols);
105: PetscCall(VecRestoreArrayRead(X, &x));
106: /* CUTEST_ush() uses the fixed ordering from CUTEST_ushp(). Preserve its first nnz values and append the mirrored entries. */
107: nz = nnz;
108: for (PetscInt k = 0; k < nnz; ++k) {
109: if (user->rows[k] != user->cols[k]) user->values[nz++] = user->values[k];
110: }
111: PetscCall(MatSetValuesCOO(H, user->values, INSERT_VALUES));
112: }
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: static PetscErrorCode CUTEstCreateHessian(Vec X, PetscBool shell, CUTEstCtx *user, Mat *H)
117: {
118: PetscFunctionBeginUser;
119: if (shell) {
120: PetscCall(VecDuplicate(X, &user->hessian_x));
121: PetscCall(MatCreateShell(PETSC_COMM_SELF, user->n, user->n, user->n, user->n, user, H));
122: PetscCall(MatShellSetOperation(*H, MATOP_MULT, (PetscErrorCodeFn *)CUTEstHessianMult));
123: PetscCall(MatShellSetOperation(*H, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)CUTEstHessianMult));
124: } else {
125: PetscInt *rows, *cols;
126: PetscCount nz;
127: integer nnz;
129: PetscCallCUTEst(CUTEST_udimsh, &user->lh);
130: nz = 2 * (PetscCount)user->lh + user->n;
131: PetscCall(PetscCalloc3(user->lh, &user->rows, user->lh, &user->cols, nz, &user->values));
132: PetscCallCUTEst(CUTEST_ushp, &user->n, &nnz, &user->lh, user->rows, user->cols);
133: PetscCall(PetscMalloc2(nz, &rows, nz, &cols));
134: for (PetscInt k = 0; k < nnz; ++k) {
135: rows[k] = user->rows[k] - 1;
136: cols[k] = user->cols[k] - 1;
137: }
138: nz = nnz;
139: for (PetscInt k = 0; k < nnz; ++k) {
140: if (rows[k] != cols[k]) {
141: rows[nz] = cols[k];
142: cols[nz++] = rows[k];
143: }
144: }
145: /* CUTEst can omit diagonals of linear variables. Keep MatShift() from changing the COO structure. */
146: for (PetscInt i = 0; i < user->n; ++i) {
147: rows[nz] = i;
148: cols[nz++] = i;
149: }
150: PetscCall(MatCreate(PETSC_COMM_SELF, H));
151: PetscCall(MatSetSizes(*H, user->n, user->n, user->n, user->n));
152: PetscCall(MatSetType(*H, MATAIJ));
153: PetscCall(MatSetFromOptions(*H));
154: PetscCall(MatSetPreallocationCOO(*H, nz, rows, cols));
155: PetscCall(PetscFree2(rows, cols));
156: }
157: PetscCall(MatSetOption(*H, MAT_SYMMETRIC, PETSC_TRUE));
158: PetscCall(MatSetOption(*H, MAT_SYMMETRY_ETERNAL, PETSC_TRUE));
159: PetscFunctionReturn(PETSC_SUCCESS);
160: }
162: static PetscErrorCode CUTEstLoadProblem(const char library[], const char data[], CUTEstCtx *user, PetscDLLibrary *dll, Vec *X)
163: {
164: const char *symbols[] = {"elfun_", "group_", "range_"};
165: char fullpath[PETSC_MAX_PATH_LEN], name[FSTRING_LEN + 1];
166: PetscScalar *x, *lower, *upper;
167: integer *types;
168: integer input = 55, output = 6, buffer = 77, m, status;
170: PetscFunctionBeginUser;
171: PetscCheck(library[0], PETSC_COMM_SELF, PETSC_ERR_USER_INPUT, "Specify the decoded problem library with -cutest_lib");
172: PetscCall(PetscGetFullPath(library, fullpath, sizeof(fullpath)));
173: /* CUTEst's loader has no status argument, so validate its required symbols first. */
174: for (PetscInt i = 0; i < 3; ++i) {
175: void *symbol;
177: PetscCall(PetscDLLibrarySym(PETSC_COMM_SELF, dll, fullpath, symbols[i], &symbol));
178: PetscCheck(symbol, PETSC_COMM_SELF, PETSC_ERR_LIB, "Decoded problem library %s does not provide %s", library, symbols[i]);
179: }
180: CUTEST_load_routines(fullpath);
181: FORTRAN_open(&input, data, &status);
182: PetscCheck(!status, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open CUTEst data file %s: status %d", data, status);
183: PetscCallCUTEst(CUTEST_cdimen, &input, &user->n, &m);
184: PetscCheck(!m, PETSC_COMM_SELF, PETSC_ERR_SUP, "This driver supports unconstrained problems; the problem has %d constraints", m);
185: PetscCall(VecCreateSeq(PETSC_COMM_SELF, user->n, X));
186: PetscCall(PetscMalloc2(user->n, &lower, user->n, &upper));
187: PetscCall(VecGetArray(*X, &x));
188: PetscCallCUTEst(CUTEST_usetup, &input, &output, &buffer, &user->n, x, lower, upper);
189: PetscCall(VecRestoreArray(*X, &x));
190: FORTRAN_close(&input, &status);
191: PetscCheck(!status, PETSC_COMM_SELF, PETSC_ERR_LIB, "Cannot close CUTEst data file: status %d", status);
192: for (PetscInt i = 0; i < user->n; ++i) PetscCheck(lower[i] <= -CUTE_INF && upper[i] >= CUTE_INF, PETSC_COMM_SELF, PETSC_ERR_SUP, "Variable %" PetscInt_FMT " has finite bounds; this driver supports unconstrained problems", i);
193: PetscCall(PetscFree2(lower, upper));
194: PetscCall(PetscMalloc1(user->n, &types));
195: PetscCallCUTEst(CUTEST_uvartype, &user->n, types);
196: for (PetscInt i = 0; i < user->n; ++i) PetscCheck(!types[i], PETSC_COMM_SELF, PETSC_ERR_SUP, "Variable %" PetscInt_FMT " is discrete; this driver supports continuous problems", i);
197: PetscCall(PetscFree(types));
198: PetscCallCUTEst(CUTEST_probname, name);
199: name[FSTRING_LEN] = '\0';
200: for (PetscInt i = FSTRING_LEN; i > 0 && name[i - 1] == ' '; --i) name[i - 1] = '\0';
201: PetscCall(PetscPrintf(PETSC_COMM_SELF, "CUTEst problem %s: %d variables\n", name, user->n));
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
205: static PetscErrorCode CUTEstUnloadProblem(PetscDLLibrary dll)
206: {
207: integer status;
209: PetscFunctionBeginUser;
210: CUTEST_uterminate(&status);
211: PetscCheck(!status, PETSC_COMM_SELF, PETSC_ERR_LIB, "CUTEST_uterminate() returned status %d", status);
212: CUTEST_unload_routines();
213: PetscCall(PetscDLLibraryClose(dll));
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }