Actual source code: taotermtest1.c
1: const char help[] = "TaoTerm coverage test comparing TaoTerm interface with traditional callbacks for Rosenbrock problem.\n\
2: Tests different TaoTerm configurations for L1, and HALFL2SQUARED types with various matrix and parameter options.\n";
4: #include <petsctao.h>
5: #include "../unconstrained/tutorials/rosenbrock4.h"
7: typedef struct {
8: AppCtx user; /* Note: AppCtx is a pointer type in rosenbrock4.h */
9: AppCtx user2; /* Second user context for callback version */
10: PetscBool test_print;
11: PetscBool test_print_map;
12: PetscBool use_term1;
13: PetscBool term1_has_A;
14: PetscBool term1_has_params;
15: PetscReal term1_scale;
16: PetscBool use_term2;
17: PetscBool term2_has_A;
18: PetscBool term2_has_params;
19: PetscReal term2_scale;
20: PetscReal term1_scale_callback; /* Callback-only scale options (for testing equivalence with TaoTerm scales) */
21: PetscReal term2_scale_callback;
22: PetscInt map_row_size;
23: PetscBool print_debug;
24: TaoTerm term1, term2;
25: Vec term1_params, term2_params;
26: Vec X_mapped1, X_mapped2, G_mapped1, G_mapped2, G_work;
27: Mat term1_A_callback; /* AIJ version for callbacks */
28: Mat term2_A_callback; /* AIJ version for callbacks */
29: } TestCtx;
31: /* Forward declarations */
32: static PetscErrorCode TestCtxInitialize(MPI_Comm, TestCtx *);
33: static PetscErrorCode TestCtxFinalize(TestCtx *);
34: static PetscErrorCode CreateTaoTermWithOptions(TestCtx *, TaoTerm *, Vec *, Mat *, const char *, const char *, PetscBool, PetscBool);
35: static PetscErrorCode FormFunctionGradient_TaoTerm(Tao, Vec, PetscReal *, Vec, void *);
36: static PetscErrorCode FormHessian_TaoTerm(Tao, Vec, Mat, Mat, void *);
37: static PetscErrorCode FormFunctionGradient_Callbacks(Tao, Vec, PetscReal *, Vec, void *);
38: static PetscErrorCode FormHessian_Callbacks(Tao, Vec, Mat, Mat, void *);
39: static PetscErrorCode CompareSolutions(Tao, Tao, TestCtx *);
41: int main(int argc, char **argv)
42: {
43: TestCtx ctx;
44: Tao tao_term, tao_callback;
45: Vec x_term, x_callback;
46: Mat H_term, H_callback;
47: TaoTerm term1, term2;
48: Vec term1_params = NULL;
49: Vec term2_params = NULL;
50: Vec term1_params_callback = NULL;
51: Vec term2_params_callback = NULL;
52: Mat term1_A, term2_A;
53: MPI_Comm comm;
55: PetscFunctionBeginUser;
56: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
57: comm = PETSC_COMM_WORLD;
59: PetscCall(TestCtxInitialize(comm, &ctx));
60: PetscCall(TaoCreate(comm, &tao_term));
61: PetscCall(TaoSetType(tao_term, TAOLMVM));
63: /* Create Rosenbrock objective using traditional TaoSet interface with user context */
64: PetscCall(CreateHessian(ctx.user, &H_term));
65: PetscCall(CreateVectors(ctx.user, H_term, &x_term, NULL));
66: PetscCall(VecZeroEntries(x_term));
67: PetscCall(TaoSetSolution(tao_term, x_term));
68: PetscCall(TaoSetObjectiveAndGradient(tao_term, NULL, FormFunctionGradient_TaoTerm, &ctx));
69: PetscCall(TaoSetHessian(tao_term, H_term, H_term, FormHessian_TaoTerm, &ctx));
71: /* Add term 1 if requested */
72: if (ctx.use_term1) {
73: PetscCall(CreateTaoTermWithOptions(&ctx, &term1, &term1_params, &term1_A, "reg1_", "A1_", ctx.term1_has_A, ctx.term1_has_params));
74: if (ctx.test_print) PetscCall(PetscObjectSetName((PetscObject)term1, "Regularizer TaoTerm"));
75: if (ctx.test_print_map) PetscCall(PetscObjectSetName((PetscObject)term1_A, "Regularizer TaoTerm Map"));
76: PetscCall(TaoAddTerm(tao_term, "reg1_", ctx.term1_scale, term1, term1_params, term1_A));
77: PetscCall(TaoTermDestroy(&term1));
78: }
79: /* Add term 2 if requested */
80: if (ctx.use_term2) {
81: PetscCall(CreateTaoTermWithOptions(&ctx, &term2, &term2_params, &term2_A, "reg2_", "A2_", ctx.term2_has_A, ctx.term2_has_params));
82: PetscCall(TaoAddTerm(tao_term, "reg2_", ctx.term2_scale, term2, term2_params, term2_A));
83: PetscCall(TaoTermDestroy(&term2));
84: }
86: PetscCall(TaoSetFromOptions(tao_term));
87: PetscCall(TaoSolve(tao_term));
89: /* Setup Tao with traditional callbacks */
90: PetscCall(TaoCreate(comm, &tao_callback));
91: PetscCall(TaoSetType(tao_callback, TAOLMVM));
93: /* Create Rosenbrock objective using callbacks that manually add term evaluations with user2 context */
94: PetscCall(CreateHessian(ctx.user2, &H_callback));
95: PetscCall(CreateVectors(ctx.user2, H_callback, &x_callback, NULL));
96: PetscCall(VecZeroEntries(x_callback));
97: PetscCall(TaoSetSolution(tao_callback, x_callback));
98: PetscCall(TaoSetObjectiveAndGradient(tao_callback, NULL, FormFunctionGradient_Callbacks, &ctx));
99: PetscCall(TaoSetHessian(tao_callback, H_callback, H_callback, FormHessian_Callbacks, &ctx));
101: /* creating duplicate term1 for callback if requested */
102: if (ctx.use_term1) {
103: if (term1_params) {
104: PetscCall(VecDuplicate(term1_params, &term1_params_callback));
105: PetscCall(VecCopy(term1_params, term1_params_callback));
106: }
107: PetscCall(TaoTermCreate(comm, &term1));
108: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)term1, "reg1_"));
109: PetscCall(TaoTermSetSolutionSizes(term1, PETSC_DECIDE, ctx.user->n, 1));
110: PetscCall(TaoTermSetFromOptions(term1));
111: /* Store for callback implementation */
112: ctx.term1 = term1;
113: ctx.term1_params = term1_params_callback;
115: if (term1_A) PetscCall(MatDuplicate(term1_A, MAT_COPY_VALUES, &ctx.term1_A_callback));
116: }
117: /* creating duplicate term2 for callback if requested */
118: if (ctx.use_term2) {
119: if (term2_params) {
120: PetscCall(VecDuplicate(term2_params, &term2_params_callback));
121: PetscCall(VecCopy(term2_params, term2_params_callback));
122: }
123: PetscCall(TaoTermCreate(comm, &term2));
124: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)term2, "reg2_"));
125: PetscCall(TaoTermSetSolutionSizes(term2, PETSC_DECIDE, ctx.user->n, 1));
126: PetscCall(TaoTermSetFromOptions(term2));
127: /* Store for callback implementation */
128: ctx.term2 = term2;
129: ctx.term2_params = term2_params_callback;
131: if (term2_A) PetscCall(MatDuplicate(term2_A, MAT_COPY_VALUES, &ctx.term2_A_callback));
132: }
134: if (ctx.use_term1) {
135: if (ctx.term1_has_A) PetscCall(MatCreateVecs(ctx.term1_A_callback, NULL, &ctx.X_mapped1));
136: else PetscCall(VecDuplicate(x_callback, &ctx.X_mapped1));
137: PetscCall(VecDuplicate(ctx.X_mapped1, &ctx.G_mapped1));
138: }
139: if (ctx.use_term2) {
140: if (ctx.term2_has_A) PetscCall(MatCreateVecs(ctx.term2_A_callback, NULL, &ctx.X_mapped2));
141: else PetscCall(VecDuplicate(x_callback, &ctx.X_mapped2));
142: PetscCall(VecDuplicate(ctx.X_mapped2, &ctx.G_mapped2));
143: }
145: PetscCall(VecDuplicate(x_callback, &ctx.G_work));
146: PetscCall(TaoSetFromOptions(tao_callback));
147: if (ctx.print_debug) PetscCall(PetscPrintf(comm, "Solving Callback version \n"));
148: PetscCall(TaoSolve(tao_callback));
150: /* Compare solutions */
151: PetscCall(CompareSolutions(tao_term, tao_callback, &ctx));
153: if (ctx.use_term1) {
154: PetscCall(VecDestroy(&term1_params));
155: PetscCall(MatDestroy(&term1_A));
156: PetscCall(MatDestroy(&ctx.term1_A_callback));
157: PetscCall(VecDestroy(&ctx.term1_params));
158: PetscCall(TaoTermDestroy(&term1));
159: }
160: if (ctx.use_term2) {
161: PetscCall(VecDestroy(&term2_params));
162: PetscCall(MatDestroy(&term2_A));
163: PetscCall(MatDestroy(&ctx.term2_A_callback));
164: PetscCall(VecDestroy(&ctx.term2_params));
165: PetscCall(TaoTermDestroy(&term2));
166: }
167: PetscCall(TaoDestroy(&tao_term));
168: PetscCall(TaoDestroy(&tao_callback));
169: PetscCall(VecDestroy(&x_term));
170: PetscCall(VecDestroy(&x_callback));
171: PetscCall(MatDestroy(&H_term));
172: PetscCall(MatDestroy(&H_callback));
173: PetscCall(TestCtxFinalize(&ctx));
174: PetscCall(PetscFinalize());
175: return 0;
176: }
178: static PetscErrorCode TestCtxInitialize(MPI_Comm comm, TestCtx *ctx)
179: {
180: PetscFunctionBeginUser;
181: PetscCall(PetscMemzero(ctx, sizeof(TestCtx)));
183: /* Initialize Rosenbrock contexts */
184: PetscCall(AppCtxCreate(comm, &ctx->user));
185: PetscCall(AppCtxCreate(comm, &ctx->user2));
187: /* Default configuration */
188: ctx->test_print = PETSC_FALSE;
189: ctx->test_print_map = PETSC_FALSE;
190: ctx->use_term1 = PETSC_FALSE;
191: ctx->use_term2 = PETSC_FALSE;
192: ctx->term1_has_A = PETSC_FALSE;
193: ctx->term1_has_params = PETSC_FALSE;
194: ctx->term2_has_A = PETSC_FALSE;
195: ctx->term2_has_params = PETSC_FALSE;
196: ctx->term1_scale = 0.1;
197: ctx->term2_scale = 0.05;
198: ctx->term1_scale_callback = 0;
199: ctx->term2_scale_callback = 0;
200: ctx->print_debug = PETSC_FALSE;
201: ctx->map_row_size = ctx->user->n - 1;
203: PetscOptionsBegin(comm, "", "TaoTerm Coverage Test Options", "TAO");
204: PetscCall(PetscOptionsBool("-test_print", "Test TaoView of term with name", "", ctx->test_print, &ctx->test_print, NULL));
205: PetscCall(PetscOptionsBool("-test_print_map", "Test TaoView of map of term with name", "", ctx->test_print_map, &ctx->test_print_map, NULL));
206: PetscCall(PetscOptionsBool("-use_term1", "Use first additional term", "", ctx->use_term1, &ctx->use_term1, NULL));
207: PetscCall(PetscOptionsBool("-use_term2", "Use second additional term", "", ctx->use_term2, &ctx->use_term2, NULL));
208: PetscCall(PetscOptionsBool("-term1_has_A", "Term 1 has a map matrix A", "", ctx->term1_has_A, &ctx->term1_has_A, NULL));
209: PetscCall(PetscOptionsBool("-term1_has_params", "Term 1 has parameters", "", ctx->term1_has_params, &ctx->term1_has_params, NULL));
210: PetscCall(PetscOptionsBool("-term2_has_A", "Term 2 has a map matrix A", "", ctx->term2_has_A, &ctx->term2_has_A, NULL));
211: PetscCall(PetscOptionsBool("-term2_has_params", "Term 2 has parameters", "", ctx->term2_has_params, &ctx->term2_has_params, NULL));
212: PetscCall(PetscOptionsReal("-term1_scale", "Scaling for term 1", "", ctx->term1_scale, &ctx->term1_scale, NULL));
213: PetscCall(PetscOptionsReal("-term2_scale", "Scaling for term 2", "", ctx->term2_scale, &ctx->term2_scale, NULL));
214: PetscCall(PetscOptionsReal("-term1_scale_callback", "Scaling for term 1 in callback version", "", ctx->term1_scale_callback, &ctx->term1_scale_callback, NULL));
215: PetscCall(PetscOptionsReal("-term2_scale_callback", "Scaling for term 2 in callback version", "", ctx->term2_scale_callback, &ctx->term2_scale_callback, NULL));
216: PetscCall(PetscOptionsInt("-map_row_size", "Row size of mapping matrix", "", ctx->map_row_size, &ctx->map_row_size, NULL));
217: PetscCall(PetscOptionsBool("-print_debug", "Print extra floating point for comparison", "", ctx->print_debug, &ctx->print_debug, NULL));
218: PetscOptionsEnd();
220: if (ctx->term1_scale_callback == 0) ctx->term1_scale_callback = ctx->term1_scale;
221: if (ctx->term2_scale_callback == 0) ctx->term2_scale_callback = ctx->term2_scale;
223: ctx->term1 = NULL;
224: ctx->term2 = NULL;
225: ctx->term1_params = NULL;
226: ctx->term2_params = NULL;
227: ctx->term1_A_callback = NULL;
228: ctx->term2_A_callback = NULL;
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
232: /* Finalize test context */
233: static PetscErrorCode TestCtxFinalize(TestCtx *ctx)
234: {
235: PetscFunctionBeginUser;
236: PetscCall(VecDestroy(&ctx->X_mapped1));
237: PetscCall(VecDestroy(&ctx->X_mapped2));
238: PetscCall(VecDestroy(&ctx->G_mapped1));
239: PetscCall(VecDestroy(&ctx->G_mapped2));
240: PetscCall(VecDestroy(&ctx->G_work));
241: PetscCall(AppCtxDestroy(&ctx->user));
242: PetscCall(AppCtxDestroy(&ctx->user2));
243: PetscFunctionReturn(PETSC_SUCCESS);
244: }
246: static PetscErrorCode CreateTaoTermWithOptions(TestCtx *ctx, TaoTerm *term, Vec *params, Mat *A, const char *term_prefix, const char *A_prefix, PetscBool has_A, PetscBool has_params)
247: {
248: MPI_Comm comm = ctx->user->comm;
249: PetscMPIInt size;
251: PetscFunctionBeginUser;
252: *term = NULL;
253: *params = NULL;
254: *A = NULL;
256: PetscCallMPI(MPI_Comm_size(comm, &size));
257: /* Create parameters if requested */
258: if (has_params) {
259: PetscCall(VecCreate(comm, params));
260: PetscCall(VecSetSizes(*params, PETSC_DECIDE, has_A ? ctx->map_row_size : ctx->user->n));
261: PetscCall(VecSetFromOptions(*params));
262: PetscCall(VecSetRandom(*params, NULL));
263: }
265: /* Create map matrix A if requested */
266: if (has_A) {
267: PetscCall(MatCreate(comm, A));
268: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*A, A_prefix));
269: PetscCall(MatSetSizes(*A, PETSC_DECIDE, PETSC_DECIDE, ctx->map_row_size, ctx->user->n));
270: PetscCall(MatSetType(*A, MATAIJ)); /* Set default type before SetFromOptions */
271: PetscCall(MatSetFromOptions(*A));
272: /* Check matrix type and set up accordingly */
273: if (size == 1) PetscCall(MatSeqAIJSetPreallocation(*A, PETSC_DEFAULT, NULL));
274: else PetscCall(MatMPIAIJSetPreallocation(*A, 5, NULL, 5, NULL));
275: PetscCall(MatSetUp(*A));
276: PetscCall(MatSetRandom(*A, NULL));
277: PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
278: PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
279: }
280: /* Create TaoTerm, set prefix, and configure from options */
281: PetscCall(TaoTermCreate(comm, term));
282: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*term, term_prefix));
283: PetscCall(TaoTermSetSolutionSizes(*term, PETSC_DECIDE, has_A ? ctx->map_row_size : ctx->user->n, 1));
284: PetscCall(TaoTermSetFromOptions(*term));
285: PetscFunctionReturn(PETSC_SUCCESS);
286: }
288: static PetscErrorCode FormFunctionGradient_TaoTerm(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
289: {
290: TestCtx *ctx = (TestCtx *)ptr;
292: PetscFunctionBeginUser;
293: PetscCall(FormObjectiveGradient(tao, X, f, G, ctx->user));
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: static PetscErrorCode FormHessian_TaoTerm(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr)
298: {
299: TestCtx *ctx = (TestCtx *)ptr;
301: PetscFunctionBeginUser;
302: PetscCall(FormHessian(tao, X, H, Hpre, ctx->user));
303: PetscFunctionReturn(PETSC_SUCCESS);
304: }
306: /* Form function and gradient for callback version (Rosenbrock + terms manually) */
307: static PetscErrorCode FormFunctionGradient_Callbacks(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
308: {
309: TestCtx *ctx = (TestCtx *)ptr;
310: PetscReal f_term;
312: PetscFunctionBeginUser;
313: /* Compute Rosenbrock part */
314: PetscCall(FormObjectiveGradient(tao, X, f, G, ctx->user2));
315: /* Add term 1 contribution */
316: if (ctx->use_term1 && ctx->term1) {
317: /* Map X if needed */
318: if (ctx->term1_A_callback) PetscCall(MatMult(ctx->term1_A_callback, X, ctx->X_mapped1));
319: else PetscCall(VecCopy(X, ctx->X_mapped1));
321: /* Compute term objective and gradient */
322: if (ctx->term1_A_callback) {
323: PetscCall(TaoTermComputeObjectiveAndGradient(ctx->term1, ctx->X_mapped1, ctx->term1_params, &f_term, ctx->G_mapped1));
324: PetscCall(MatMultTranspose(ctx->term1_A_callback, ctx->G_mapped1, ctx->G_work));
325: PetscCall(VecAXPY(G, ctx->term1_scale_callback, ctx->G_work));
326: } else {
327: PetscCall(TaoTermComputeObjectiveAndGradient(ctx->term1, ctx->X_mapped1, ctx->term1_params, &f_term, ctx->G_work));
328: PetscCall(VecAXPY(G, ctx->term1_scale_callback, ctx->G_work));
329: }
330: *f += ctx->term1_scale_callback * f_term;
331: }
333: /* Add term 2 contribution */
334: if (ctx->use_term2 && ctx->term2) {
335: /* Map X if needed */
336: if (ctx->term2_A_callback) PetscCall(MatMult(ctx->term2_A_callback, X, ctx->X_mapped2));
337: else PetscCall(VecCopy(X, ctx->X_mapped2));
339: /* Compute term objective and gradient */
340: if (ctx->term2_A_callback) {
341: PetscCall(TaoTermComputeObjectiveAndGradient(ctx->term2, ctx->X_mapped2, ctx->term2_params, &f_term, ctx->G_mapped2));
342: /* Map gradient back and add to G */
343: PetscCall(MatMultTranspose(ctx->term2_A_callback, ctx->G_mapped2, ctx->G_work));
344: PetscCall(VecAXPY(G, ctx->term2_scale_callback, ctx->G_work));
345: } else {
346: PetscCall(TaoTermComputeObjectiveAndGradient(ctx->term2, ctx->X_mapped2, ctx->term2_params, &f_term, ctx->G_work));
347: PetscCall(VecAXPY(G, ctx->term2_scale_callback, ctx->G_work));
348: }
349: *f += ctx->term2_scale_callback * f_term;
350: }
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }
354: /* Form Hessian for callback version (Rosenbrock + terms manually) */
355: static PetscErrorCode FormHessian_Callbacks(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr)
356: {
357: TestCtx *ctx = (TestCtx *)ptr;
358: Mat H_term;
359: PetscInt m, n;
360: PetscBool is_assembled;
362: PetscFunctionBeginUser;
363: /* Compute Rosenbrock Hessian */
364: PetscCall(FormHessian(tao, X, H, Hpre, ctx->user2));
365: /* Add term 1 Hessian contribution */
366: if (ctx->use_term1 && ctx->term1) {
367: if (ctx->term1_A_callback) {
368: PetscCall(MatMult(ctx->term1_A_callback, X, ctx->X_mapped1));
369: PetscCall(VecGetSize(ctx->X_mapped1, &m));
370: PetscCall(MatCreate(ctx->user2->comm, &H_term));
371: PetscCall(MatSetSizes(H_term, PETSC_DECIDE, PETSC_DECIDE, m, m));
372: PetscCall(MatSetType(H_term, MATAIJ));
373: PetscCall(MatSetUp(H_term));
374: } else {
375: PetscCall(VecCopy(X, ctx->X_mapped1));
376: PetscCall(VecGetSize(X, &n));
377: PetscCall(MatCreate(ctx->user2->comm, &H_term));
378: PetscCall(MatSetSizes(H_term, PETSC_DECIDE, PETSC_DECIDE, n, n));
379: PetscCall(MatSetType(H_term, MATAIJ));
380: PetscCall(MatSetUp(H_term));
381: }
382: PetscCall(MatAssembled(H_term, &is_assembled));
383: if (!is_assembled) {
384: PetscCall(MatAssemblyBegin(H_term, MAT_FINAL_ASSEMBLY));
385: PetscCall(MatAssemblyEnd(H_term, MAT_FINAL_ASSEMBLY));
386: }
388: PetscCall(TaoTermComputeHessian(ctx->term1, ctx->X_mapped1, ctx->term1_params, H_term, NULL));
389: if (ctx->term1_A_callback) {
390: Mat H_mapped;
391: /* H = H + scale * A^T * H_term * A */
392: PetscCall(MatMatMult(H_term, ctx->term1_A_callback, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &H_mapped));
393: PetscCall(MatDestroy(&H_term));
394: PetscCall(MatTransposeMatMult(ctx->term1_A_callback, H_mapped, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &H_term));
395: PetscCall(MatDestroy(&H_mapped));
396: PetscCall(MatScale(H_term, ctx->term1_scale_callback));
397: PetscCall(MatAXPY(H, 1.0, H_term, DIFFERENT_NONZERO_PATTERN));
398: PetscCall(MatDestroy(&H_term));
399: } else {
400: PetscCall(MatScale(H_term, ctx->term1_scale_callback));
401: PetscCall(MatAXPY(H, 1.0, H_term, DIFFERENT_NONZERO_PATTERN));
402: PetscCall(MatDestroy(&H_term));
403: }
404: }
406: /* Add term 2 Hessian contribution */
407: if (ctx->use_term2 && ctx->term2) {
408: if (ctx->term2_A_callback) {
409: PetscCall(MatMult(ctx->term2_A_callback, X, ctx->X_mapped2));
410: PetscCall(VecGetSize(ctx->X_mapped2, &m));
411: PetscCall(MatCreate(ctx->user2->comm, &H_term));
412: PetscCall(MatSetSizes(H_term, PETSC_DECIDE, PETSC_DECIDE, m, m));
413: PetscCall(MatSetType(H_term, MATAIJ));
414: PetscCall(MatSetUp(H_term));
415: } else {
416: PetscCall(VecCopy(X, ctx->X_mapped2));
417: PetscCall(VecGetSize(X, &n));
418: PetscCall(MatCreate(ctx->user2->comm, &H_term));
419: PetscCall(MatSetSizes(H_term, PETSC_DECIDE, PETSC_DECIDE, n, n));
420: PetscCall(MatSetType(H_term, MATAIJ));
421: PetscCall(MatSetUp(H_term));
422: }
423: PetscCall(MatAssembled(H_term, &is_assembled));
424: if (!is_assembled) {
425: PetscCall(MatAssemblyBegin(H_term, MAT_FINAL_ASSEMBLY));
426: PetscCall(MatAssemblyEnd(H_term, MAT_FINAL_ASSEMBLY));
427: }
429: PetscCall(TaoTermComputeHessian(ctx->term2, ctx->X_mapped2, ctx->term2_params, H_term, NULL));
431: if (ctx->term2_A_callback) {
432: Mat H_mapped;
433: /* H = H + scale * A^T * H_term * A */
434: PetscCall(MatMatMult(H_term, ctx->term2_A_callback, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &H_mapped));
435: PetscCall(MatDestroy(&H_term));
436: PetscCall(MatTransposeMatMult(ctx->term2_A_callback, H_mapped, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &H_term));
437: PetscCall(MatDestroy(&H_mapped));
438: PetscCall(MatScale(H_term, ctx->term2_scale_callback));
439: PetscCall(MatAXPY(H, 1.0, H_term, DIFFERENT_NONZERO_PATTERN));
440: PetscCall(MatDestroy(&H_term));
441: } else {
442: PetscCall(MatScale(H_term, ctx->term2_scale_callback));
443: PetscCall(MatAXPY(H, 1.0, H_term, DIFFERENT_NONZERO_PATTERN));
444: PetscCall(MatDestroy(&H_term));
445: }
446: }
447: PetscFunctionReturn(PETSC_SUCCESS);
448: }
450: /* Compare solutions from both methods */
451: static PetscErrorCode CompareSolutions(Tao tao_term, Tao tao_callback, TestCtx *ctx)
452: {
453: Vec x_term, x_callback, diff;
454: PetscReal norm_term, norm_callback, norm_diff, rel_diff;
455: PetscReal f_term, f_callback;
456: TaoConvergedReason reason_term, reason_callback;
458: PetscFunctionBeginUser;
459: PetscCall(TaoGetSolution(tao_term, &x_term));
460: PetscCall(TaoGetSolution(tao_callback, &x_callback));
462: /* Compute difference */
463: PetscCall(VecDuplicate(x_term, &diff));
464: PetscCall(VecCopy(x_term, diff));
465: PetscCall(VecAXPY(diff, -1.0, x_callback));
467: PetscCall(VecNorm(x_term, NORM_2, &norm_term));
468: PetscCall(VecNorm(x_callback, NORM_2, &norm_callback));
469: PetscCall(VecNorm(diff, NORM_2, &norm_diff));
471: rel_diff = norm_diff / PetscMax(norm_term, 1.0e-10);
473: /* Get objective values */
474: PetscCall(TaoGetSolutionStatus(tao_term, NULL, &f_term, NULL, NULL, NULL, NULL));
475: PetscCall(TaoGetSolutionStatus(tao_callback, NULL, &f_callback, NULL, NULL, NULL, NULL));
477: /* Get convergence reasons */
478: PetscCall(TaoGetConvergedReason(tao_term, &reason_term));
479: PetscCall(TaoGetConvergedReason(tao_callback, &reason_callback));
481: if (rel_diff < 1.0e-12) {
482: if (ctx->print_debug) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)tao_term), "Solutions match (relative difference < 1e-12), %6.10e\n", (double)rel_diff));
483: else PetscCall(PetscPrintf(PetscObjectComm((PetscObject)tao_term), "Solutions match (relative difference < 1e-12)\n"));
484: } else if (rel_diff < 1.0e-6) {
485: if (ctx->print_debug) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)tao_term), "Solutions approximately match (relative difference < 1e-6), %6.10e\n", (double)rel_diff));
486: else PetscCall(PetscPrintf(PetscObjectComm((PetscObject)tao_term), "Solutions approximately match (relative difference < 1e-6)\n"));
487: } else {
488: if (ctx->print_debug) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)tao_term), "Solutions differ significantly (relative difference >= 1e-6), %6.10e\n", (double)rel_diff));
489: else PetscCall(PetscPrintf(PetscObjectComm((PetscObject)tao_term), "Solutions differ significantly (relative difference >= 1e-6)\n"));
490: }
491: PetscCall(VecDestroy(&diff));
492: PetscFunctionReturn(PETSC_SUCCESS);
493: }
495: /*TEST
497: build:
498: requires: !complex !single !quad !defined(PETSC_USE_64BIT_INDICES) !__float128
500: test:
501: suffix: test_print_noparam
502: args: -tao_view -tao_type nls -use_term1 -term1_has_params 0 -term1_has_A
503: args: -reg1_tao_term_type l1 -test_print -test_print_map -reg1_tao_term_parameters_mode none
505: test:
506: suffix: test_print_p_req
507: args: -tao_view -tao_type nls -use_term1 -term1_has_params 1 -term1_has_A
508: args: -reg1_tao_term_type l1 -test_print -test_print_map -reg1_tao_term_parameters_mode required
510: # Single term tests
511: # L1 with A
512: testset:
513: suffix: term1_l1_A
514: output_file: output/taotermtest1.out
515: args: -tao_type nls -use_term1 -term1_has_params {{0 1}} -term1_has_A
516: args: -reg1_tao_term_type l1
517: test:
518: args: -term1_scale 0.012
519: test:
520: args: -tao_term_sum_reg1_scale 0.012 -term1_scale_callback 0.012
522: # L1 without A
523: testset:
524: suffix: term1_l1_no_A
525: output_file: output/taotermtest1.out
526: args: -tao_type nls -use_term1 -term1_has_params {{0 1}}
527: args: -term1_has_A 0 -reg1_tao_term_type l1
528: test:
529: args: -term1_scale 0.012
530: test:
531: args: -tao_term_sum_reg1_scale 0.012 -term1_scale_callback 0.012
533: # L2 with A
534: testset:
535: suffix: term1_l2_A
536: output_file: output/taotermtest1.out
537: args: -tao_type nls -use_term1 -term1_has_params {{0 1}} -term1_has_A
538: args: -reg1_tao_term_type halfl2squared
539: test:
540: args: -term1_scale 0.12
541: test:
542: args: -tao_term_sum_reg1_scale 0.12 -term1_scale_callback 0.12
544: # L2 without A
545: testset:
546: suffix: term1_l2_no_A
547: output_file: output/taotermtest1.out
548: args: -tao_type nls -use_term1 -term1_has_params {{0 1}}
549: args: -term1_has_A 0
550: args: -reg1_tao_term_type halfl2squared
551: test:
552: args: -term1_scale 0.12
553: test:
554: args: -tao_term_sum_reg1_scale 0.12 -term1_scale_callback 0.12
556: # Two terms, with no mapping
557: testset:
558: suffix: l1_l2_no_A1_no_A2
559: output_file: output/taotermtest1.out
560: args: -tao_type nls -use_term1 -use_term2 -term1_has_params {{0 1}} -term2_has_params {{0 1}}
561: args: -term1_has_A 0 -term2_has_A 0 -reg1_tao_term_type l1 -reg2_tao_term_type halfl2squared
562: test:
563: args: -term1_scale 0.123 -term2_scale 0.098
564: test:
565: args: -term1_scale 0.123 -tao_term_sum_reg2_scale 0.098 -term2_scale_callback 0.098
566: test:
567: args: -tao_term_sum_reg1_scale 0.123 -term1_scale_callback 0.123 -term2_scale 0.098
568: test:
569: args: -tao_term_sum_reg1_scale 0.123 -term1_scale_callback 0.123 -tao_term_sum_reg2_scale 0.098 -term2_scale_callback 0.098
571: # Two terms, with one mapping
572: testset:
573: suffix: l1_l2_no_A1_A2
574: output_file: output/taotermtest1.out
575: args: -tao_type nls -use_term1 -use_term2 -term1_has_params {{0 1}} -term2_has_params {{0 1}}
576: args: -term1_has_A 0 -term2_has_A 1 -reg1_tao_term_type l1 -reg2_tao_term_type halfl2squared
577: test:
578: args: -term1_scale 0.123 -term2_scale 0.098
579: test:
580: args: -term1_scale 0.123 -tao_term_sum_reg2_scale 0.098 -term2_scale_callback 0.098
581: test:
582: args: -tao_term_sum_reg1_scale 0.123 -term1_scale_callback 0.123 -term2_scale 0.098
583: test:
584: args: -tao_term_sum_reg1_scale 0.123 -term1_scale_callback 0.123 -tao_term_sum_reg2_scale 0.098 -term2_scale_callback 0.098
586: testset:
587: suffix: l2_l1_no_A1_A2
588: output_file: output/taotermtest1.out
589: args: -tao_type nls -use_term1 -use_term2 -term1_has_params {{0 1}} -term2_has_params {{0 1}}
590: args: -term1_has_A 0 -term2_has_A 1 -reg1_tao_term_type halfl2squared -reg2_tao_term_type l1
591: test:
592: args: -term1_scale 0.123 -term2_scale 0.098
593: test:
594: args: -term1_scale 0.123 -tao_term_sum_reg2_scale 0.098 -term2_scale_callback 0.098
595: test:
596: args: -tao_term_sum_reg1_scale 0.123 -term1_scale_callback 0.123 -term2_scale 0.098
597: test:
598: args: -tao_term_sum_reg1_scale 0.123 -term1_scale_callback 0.123 -tao_term_sum_reg2_scale 0.098 -term2_scale_callback 0.098
600: # Two terms: term1 has A, term2 has A
601: testset:
602: suffix: l1_l2_A1_A2
603: output_file: output/taotermtest1.out
604: args: -tao_type nls -use_term1 -use_term2 -term1_has_params {{0 1}} -term2_has_params {{0 1}}
605: args: -term1_has_A 1 -term2_has_A 1 -reg1_tao_term_type l1 -reg2_tao_term_type halfl2squared
606: test:
607: args: -term1_scale 0.075 -term2_scale 0.098
608: test:
609: args: -term1_scale 0.075 -tao_term_sum_reg2_scale 0.098 -term2_scale_callback 0.098
610: test:
611: args: -tao_term_sum_reg1_scale 0.075 -term1_scale_callback 0.075 -term2_scale 0.098
612: test:
613: args: -tao_term_sum_reg1_scale 0.075 -term1_scale_callback 0.075 -tao_term_sum_reg2_scale 0.098 -term2_scale_callback 0.098
615: TEST*/