Actual source code: taotermsum.c
1: #include <petsc/private/taoimpl.h>
2: #include <../src/tao/term/impls/sum/taotermsum.h>
3: #include <ctype.h>
5: static const char *const TaoTermMasks[] = {"none", "objective", "gradient", "hessian", "TaoTermMask", "TAOTERM_MASK_", NULL};
7: typedef struct _n_TaoTerm_Sum TaoTerm_Sum;
9: typedef struct _n_TaoTermSumHessCache {
10: PetscObjectId x_id;
11: PetscObjectId p_id;
12: PetscObjectState x_state;
13: PetscObjectState p_state;
14: PetscInt n_terms;
15: Mat *hessians;
16: Vec *Axs;
17: } TaoTermSumHessCache;
19: struct _n_TaoTerm_Sum {
20: PetscInt n_terms;
21: TaoTermMapping *terms;
22: PetscReal *subterm_values;
23: TaoTermSumHessCache hessian_cache;
24: };
26: PETSC_INTERN PetscErrorCode TaoTermSumVecNestGetSubVecsRead(Vec params, PetscInt *n, Vec **subparams, PetscBool **is_dummy)
27: {
28: PetscContainer is_dummy_container = NULL;
30: PetscFunctionBegin;
31: *is_dummy = NULL;
32: PetscCall(VecNestGetSubVecsRead(params, n, subparams));
33: PetscCall(PetscObjectQuery((PetscObject)params, "__TaoTermSumParametersPack", (PetscObject *)&is_dummy_container));
34: if (is_dummy_container) PetscCall(PetscContainerGetPointer(is_dummy_container, (void **)is_dummy));
35: PetscFunctionReturn(PETSC_SUCCESS);
36: }
38: PETSC_INTERN PetscErrorCode TaoTermSumVecNestRestoreSubVecsRead(Vec params, PetscInt *n, Vec **subparams, PetscBool **is_dummy)
39: {
40: PetscFunctionBegin;
41: PetscCall(VecNestRestoreSubVecsRead(params, n, subparams));
42: PetscFunctionReturn(PETSC_SUCCESS);
43: }
45: static PetscErrorCode TaoTermSumHessCacheReset(TaoTermSumHessCache *cache)
46: {
47: PetscFunctionBegin;
48: for (PetscInt i = 0; i < cache->n_terms; i++) PetscCall(MatDestroy(&cache->hessians[i]));
49: PetscCall(PetscFree(cache->hessians));
50: for (PetscInt i = 0; i < cache->n_terms; i++) PetscCall(VecDestroy(&cache->Axs[i]));
51: PetscCall(PetscFree(cache->Axs));
52: cache->n_terms = 0;
53: cache->x_id = 0;
54: cache->p_id = 0;
55: cache->x_state = 0;
56: cache->p_state = 0;
57: PetscFunctionReturn(PETSC_SUCCESS);
58: }
60: static PetscErrorCode TaoTermSumIsDummyDestroy(PetscCtxRt ctx)
61: {
62: PetscFunctionBegin;
63: PetscCall(PetscFree(*(void **)ctx));
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
67: /*@
68: TaoTermSumParametersPack - Concatenate the parameters for terms into a `VECNEST` parameter vector for a `TAOTERMSUM`
70: Collective
72: Input Parameters:
73: + term - a `TaoTerm` of type `TAOTERMSUM`
74: - p_arr - an array of parameters `Vec`s, one for each term in the sum. An entry can be `NULL` for a term that doesn't take parameters.
76: Output Parameter:
77: . params - a `Vec` of type `VECNEST` that concatenates all of the parameters
79: Level: developer
81: Note:
82: This is a wrapper around `VecCreateNest()`, but that function does not allow `NULL` for any of the `Vec`s in the array. A 0-length
83: vector will be created for each `NULL` `Vec` that will be internally ignored by `TAOTERMSUM`.
85: .seealso: [](sec_tao_term),
86: `TaoTerm`,
87: `TAOTERMSUM`,
88: `TaoTermSumParametersUnpack()`,
89: `VECNEST`,
90: `VecNestGetTaoTermSumParameters()`,
91: `VecCreateNest()`
92: @*/
93: PetscErrorCode TaoTermSumParametersPack(TaoTerm term, Vec p_arr[], Vec *params)
94: {
95: PetscInt n_terms;
96: Vec *p;
97: PetscBool *is_dummy;
98: PetscContainer is_dummy_container;
100: PetscFunctionBegin;
102: PetscAssertPointer(p_arr, 2);
103: PetscAssertPointer(params, 3);
104: PetscCall(TaoTermSumGetNumberTerms(term, &n_terms));
105: PetscCall(PetscMalloc1(n_terms, &p));
106: PetscCall(PetscMalloc1(n_terms, &is_dummy));
107: for (PetscInt i = 0; i < n_terms; i++) {
108: if (p_arr[i]) {
110: p[i] = p_arr[i];
111: is_dummy[i] = PETSC_FALSE;
112: } else {
113: TaoTerm subterm;
114: Vec dummy_vec;
115: TaoTermParametersMode mode;
116: VecType vec_type = VECSTANDARD;
117: PetscLayout layout = NULL;
119: PetscCall(TaoTermSumGetTerm(term, i, NULL, NULL, &subterm, NULL));
120: PetscCall(TaoTermGetParametersMode(subterm, &mode));
121: if (mode != TAOTERM_PARAMETERS_NONE) {
122: PetscCall(TaoTermGetParametersVecType(subterm, &vec_type));
123: PetscCall(TaoTermGetParametersLayout(subterm, &layout));
124: layout->refcnt++;
125: } else {
126: PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)term), &layout));
127: PetscCall(PetscLayoutSetLocalSize(layout, 0));
128: PetscCall(PetscLayoutSetSize(layout, 0));
129: }
130: PetscCall(VecCreate(PetscObjectComm((PetscObject)term), &dummy_vec));
131: PetscCall(VecSetLayout(dummy_vec, layout));
132: PetscCall(PetscLayoutDestroy(&layout));
133: PetscCall(VecSetType(dummy_vec, vec_type));
134: is_dummy[i] = PETSC_TRUE;
135: p[i] = dummy_vec;
136: }
137: }
138: PetscCall(VecCreateNest(PetscObjectComm((PetscObject)term), n_terms, NULL, p, params));
139: for (PetscInt i = 0; i < n_terms; i++) {
140: if (!p_arr[i]) PetscCall(VecDestroy(&p[i]));
141: }
142: PetscCall(PetscFree(p));
143: PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)term), &is_dummy_container));
144: PetscCall(PetscContainerSetPointer(is_dummy_container, (void *)is_dummy));
145: PetscCall(PetscContainerSetCtxDestroy(is_dummy_container, TaoTermSumIsDummyDestroy));
146: PetscCall(PetscObjectCompose((PetscObject)*params, "__TaoTermSumParametersPack", (PetscObject)is_dummy_container));
147: PetscCall(PetscContainerDestroy(&is_dummy_container));
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: /*@
152: TaoTermSumParametersUnpack - Unpack the concatenated parameters created by `TaoTermSumParametersPack()` and destroy the `VECNEST`
154: Collective
156: Input Parameters:
157: + term - a `TaoTerm` of type `TAOTERMSUM`
158: - params - a `Vec` created by `TaoTermSumParametersPack()`
160: Output Parameter:
161: . p_arr - an array of parameters `Vec`s, one for each term in the sum. An entry will be `NULL` if `NULL` was passed in the same position of `TaoTermSumParametersPack()`
163: Level: intermediate
165: .seealso: [](sec_tao_term),
166: `TaoTerm`,
167: `TAOTERMSUM`,
168: `TaoTermSumParametersPack()`,
169: `VecNestGetTaoTermSumParameters()`
170: @*/
171: PetscErrorCode TaoTermSumParametersUnpack(TaoTerm term, Vec *params, Vec p_arr[])
172: {
173: PetscInt n_terms;
174: PetscBool *is_dummy = NULL;
175: PetscContainer is_dummy_container = NULL;
177: PetscFunctionBegin;
180: PetscAssertPointer(p_arr, 3);
181: PetscCall(TaoTermSumGetNumberTerms(term, &n_terms));
182: PetscCall(PetscObjectQuery((PetscObject)*params, "__TaoTermSumParametersPack", (PetscObject *)&is_dummy_container));
183: if (is_dummy_container) PetscCall(PetscContainerGetPointer(is_dummy_container, (void **)&is_dummy));
184: for (PetscInt i = 0; i < n_terms; i++) {
185: Vec subparam;
187: PetscCall(VecNestGetSubVec(*params, i, &subparam));
188: if (is_dummy && is_dummy[i]) {
189: p_arr[i] = NULL;
190: } else {
191: PetscCall(PetscObjectReference((PetscObject)subparam));
192: p_arr[i] = subparam;
193: }
194: }
195: PetscCall(VecDestroy(params));
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
199: /*@
200: VecNestGetTaoTermSumParameters - A wrapper around `VecNestGetSubVec()` for `TAOTERMSUM`.
202: Not collective
204: Input Parameters:
205: + params - a `VECNEST` that has one nested vector for each term of a `TAOTERMSUM`
206: - index - the index of a term
208: Output Parameter:
209: . subparams - the parameters of the internal terms of `TAOTERMSUM`. (may be `NULL`)
211: Level: intermediate
213: Note:
214: `VecNestGetSubVec()` cannot return `NULL` for the subvec. If `params` was
215: created by `TaoTermSumParametersPack()`, then any `NULL` subvecs that were passed
216: to that function will be returned `NULL` by this function.
218: .seealso: [](sec_tao_term),
219: `TaoTerm`,
220: `TAOTERMSUM`,
221: `TaoTermSumParametersPack()`,
222: `TaoTermSumParametersUnpack()`,
223: `VECNEST`,
224: `VecNestGetSubVec()`
225: @*/
226: PetscErrorCode VecNestGetTaoTermSumParameters(Vec params, PetscInt index, Vec *subparams)
227: {
228: PetscBool *is_dummy = NULL;
229: PetscContainer is_dummy_container = NULL;
231: PetscFunctionBegin;
233: PetscAssertPointer(subparams, 3);
234: PetscCall(PetscObjectQuery((PetscObject)params, "__TaoTermSumParametersPack", (PetscObject *)&is_dummy_container));
235: if (is_dummy_container) PetscCall(PetscContainerGetPointer(is_dummy_container, (void **)&is_dummy));
236: if (is_dummy && is_dummy[index]) {
237: *subparams = NULL;
238: } else {
239: PetscCall(VecNestGetSubVec(params, index, subparams));
240: }
241: PetscFunctionReturn(PETSC_SUCCESS);
242: }
244: static PetscErrorCode TaoTermDestroy_Sum(TaoTerm term)
245: {
246: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
248: PetscFunctionBegin;
249: for (PetscInt i = 0; i < sum->n_terms; i++) PetscCall(TaoTermMappingReset(&sum->terms[i]));
250: PetscCall(TaoTermSumHessCacheReset(&sum->hessian_cache));
251: PetscCall(PetscFree(sum->terms));
252: PetscCall(PetscFree(sum->subterm_values));
253: PetscCall(PetscFree(sum));
254: term->data = NULL;
255: PetscCall(PetscObjectReference((PetscObject)term->parameters_factory_orig));
256: PetscCall(MatDestroy(&term->parameters_factory));
257: term->parameters_factory = term->parameters_factory_orig;
258: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumGetNumberTerms_C", NULL));
259: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumSetNumberTerms_C", NULL));
260: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumGetTerm_C", NULL));
261: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumSetTerm_C", NULL));
262: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumAddTerm_C", NULL));
263: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumGetTermHessianMatrices_C", NULL));
264: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumSetTermHessianMatrices_C", NULL));
265: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumGetTermMask_C", NULL));
266: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumSetTermMask_C", NULL));
267: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumGetLastTermObjectives_C", NULL));
268: PetscFunctionReturn(PETSC_SUCCESS);
269: }
271: static PetscErrorCode TaoTermView_Sum_NameHasSpaces(const char name[], PetscBool *has_spaces)
272: {
273: size_t n;
275: PetscFunctionBegin;
276: PetscCall(PetscStrlen(name, &n));
277: for (size_t i = 0; i < n; i++) {
278: if (isspace((unsigned char)name[i])) {
279: *has_spaces = PETSC_TRUE;
280: PetscFunctionReturn(PETSC_SUCCESS);
281: }
282: }
283: *has_spaces = PETSC_FALSE;
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
287: static PetscErrorCode TaoTermViewSumPrintSubtermName(PetscViewer viewer, TaoTerm subterm, PetscInt i, const char f[], PetscBool colon_newline)
288: {
289: const char *subterm_prefix;
290: const char *subterm_name = NULL;
292: PetscFunctionBegin;
293: if (((PetscObject)subterm)->name) PetscCall(PetscObjectGetName((PetscObject)subterm, &subterm_name));
294: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)subterm, &subterm_prefix));
295: if (subterm_name) {
296: PetscBool same;
298: PetscCall(PetscStrncmp(subterm_name, "TaoTerm_", 8, &same));
299: if (same == PETSC_FALSE) {
300: PetscBool has_spaces;
302: PetscCall(TaoTermView_Sum_NameHasSpaces(subterm_name, &has_spaces));
303: if (has_spaces) PetscCall(PetscViewerASCIIPrintf(viewer, "%s_{%s}%s", f, subterm_name, colon_newline ? ":\n" : ""));
304: else PetscCall(PetscViewerASCIIPrintf(viewer, "%s%s", subterm_name, colon_newline ? ":\n" : ""));
305: PetscFunctionReturn(PETSC_SUCCESS);
306: }
307: }
308: if (subterm_prefix) PetscCall(PetscViewerASCIIPrintf(viewer, "%s_{%s}%s", f, subterm_prefix, colon_newline ? ":\n" : ""));
309: else PetscCall(PetscViewerASCIIPrintf(viewer, "%s_%" PetscInt_FMT "%s", f, i, colon_newline ? ":\n" : ""));
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: PETSC_INTERN PetscErrorCode TaoTermViewSumPrintMapName(PetscViewer viewer, Mat map, PetscInt i, const char A[], PetscBool colon_newline)
314: {
315: const char *map_prefix;
316: const char *map_name = NULL;
318: PetscFunctionBegin;
319: if (((PetscObject)map)->name) PetscCall(PetscObjectGetName((PetscObject)map, &map_name));
320: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)map, &map_prefix));
321: if (map_name) {
322: PetscBool same;
324: PetscCall(PetscStrncmp(map_name, "Mat_", 4, &same));
325: if (same == PETSC_FALSE) {
326: PetscBool has_spaces;
328: PetscCall(TaoTermView_Sum_NameHasSpaces(map_name, &has_spaces));
329: if (has_spaces) PetscCall(PetscViewerASCIIPrintf(viewer, "%s_{%s}%s", A, map_name, colon_newline ? ":\n" : ""));
330: else PetscCall(PetscViewerASCIIPrintf(viewer, "%s%s", map_name, colon_newline ? ":\n" : ""));
331: PetscFunctionReturn(PETSC_SUCCESS);
332: }
333: }
334: if (map_prefix) PetscCall(PetscViewerASCIIPrintf(viewer, "%s_{%s}%s", A, map_prefix, colon_newline ? ":\n" : ""));
335: else PetscCall(PetscViewerASCIIPrintf(viewer, "%s_%" PetscInt_FMT "%s", A, i, colon_newline ? ":\n" : ""));
336: PetscFunctionReturn(PETSC_SUCCESS);
337: }
339: PETSC_INTERN PetscErrorCode TaoTermViewSumPrintSubterm(TaoTerm term, PetscViewer viewer, Vec params, PetscInt i, PetscBool initial, PetscBool print_map, const char f[], const char A[], const char x[], const char p[])
340: {
341: PetscReal scale;
342: TaoTerm subterm;
343: Mat map;
344: TaoTermParametersMode pmode;
346: PetscFunctionBegin;
347: PetscCall(TaoTermSumGetTerm(term, i, NULL, &scale, &subterm, &map));
348: if (scale == 1.0) PetscCall(PetscViewerASCIIPrintf(viewer, "%s", initial ? "" : " + "));
349: else if (initial) PetscCall(PetscViewerASCIIPrintf(viewer, "%g ", (double)scale));
350: else PetscCall(PetscViewerASCIIPrintf(viewer, " %s %g ", scale >= 0.0 ? "+" : "-", (double)PetscAbsReal(scale)));
351: PetscCall(TaoTermViewSumPrintSubtermName(viewer, subterm, i, f, PETSC_FALSE));
352: PetscCall(PetscViewerASCIIPrintf(viewer, "("));
353: if (print_map && map) {
354: PetscCall(TaoTermViewSumPrintMapName(viewer, map, i, A, PETSC_FALSE));
355: PetscCall(PetscViewerASCIIPrintf(viewer, " "));
356: }
357: PetscCall(PetscViewerASCIIPrintf(viewer, "%s", x));
358: PetscCall(TaoTermGetParametersMode(subterm, &pmode));
359: switch (pmode) {
360: case TAOTERM_PARAMETERS_NONE:
361: break;
362: case TAOTERM_PARAMETERS_OPTIONAL:
363: PetscCall(PetscViewerASCIIPrintf(viewer, "; [p_%" PetscInt_FMT "]", i));
364: break;
365: case TAOTERM_PARAMETERS_REQUIRED:
366: PetscCall(PetscViewerASCIIPrintf(viewer, "; p_%" PetscInt_FMT, i));
367: break;
368: }
369: PetscCall(PetscViewerASCIIPrintf(viewer, ")"));
370: PetscFunctionReturn(PETSC_SUCCESS);
371: }
373: static PetscErrorCode TaoTermView_Sum_ASCII_INFO(TaoTerm term, PetscViewer viewer)
374: {
375: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
377: PetscFunctionBegin;
378: PetscCall(PetscViewerASCIIPrintf(viewer, "Sum of %" PetscInt_FMT " terms:%s", sum->n_terms, sum->n_terms > 0 ? " " : ""));
379: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
380: for (PetscInt i = 0; i < sum->n_terms; i++) PetscCall(TaoTermViewSumPrintSubterm(term, viewer, NULL, i, (i == 0) ? PETSC_TRUE : PETSC_FALSE, PETSC_TRUE, "f", "A", "x", "p"));
381: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
382: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
383: for (PetscInt i = 0; i < sum->n_terms; i++) {
384: Mat map;
385: TaoTerm subterm;
387: PetscCall(TaoTermSumGetTerm(term, i, NULL, NULL, &subterm, &map));
388: PetscCall(TaoTermViewSumPrintSubtermName(viewer, subterm, i, "f", PETSC_TRUE));
389: PetscCall(PetscViewerASCIIPushTab(viewer));
390: PetscCall(TaoTermView(subterm, viewer));
391: PetscCall(PetscViewerASCIIPopTab(viewer));
392: if (map == NULL) continue;
393: PetscCall(TaoTermViewSumPrintMapName(viewer, map, i, "A", PETSC_TRUE));
394: PetscCall(PetscViewerASCIIPushTab(viewer));
395: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
396: PetscCall(MatView(map, viewer));
397: PetscCall(PetscViewerPopFormat(viewer));
398: PetscCall(PetscViewerASCIIPopTab(viewer));
399: }
400: PetscFunctionReturn(PETSC_SUCCESS);
401: }
403: static PetscErrorCode TaoTermView_Sum(TaoTerm term, PetscViewer viewer)
404: {
405: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
406: PetscBool iascii;
408: PetscFunctionBegin;
409: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
410: if (iascii) {
411: PetscViewerFormat format;
413: PetscCall(PetscViewerGetFormat(viewer, &format));
414: if (sum->n_terms <= 3 && format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
415: PetscCall(TaoTermView_Sum_ASCII_INFO(term, viewer));
416: PetscFunctionReturn(PETSC_SUCCESS);
417: }
418: PetscCall(PetscViewerASCIIPrintf(viewer, "Sum of %" PetscInt_FMT " terms:\n", sum->n_terms));
419: PetscCall(PetscViewerASCIIPushTab(viewer));
420: for (PetscInt i = 0; i < sum->n_terms; i++) {
421: PetscReal scale;
422: const char *subprefix;
423: TaoTerm subterm;
424: Mat map;
425: TaoTermMask mask;
427: PetscCall(TaoTermSumGetTerm(term, i, &subprefix, &scale, &subterm, &map));
429: PetscCall(PetscViewerASCIIPrintf(viewer, "Summand %" PetscInt_FMT ":\n", i));
430: PetscCall(PetscViewerASCIIPushTab(viewer));
432: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscCall(PetscViewerASCIIPrintf(viewer, "Scale (tao_term_sum_%sscale): %g\n", subprefix, (double)scale));
433: else PetscCall(PetscViewerASCIIPrintf(viewer, "Scale: %g\n", (double)scale));
434: PetscCall(PetscViewerASCIIPrintf(viewer, "Term:\n"));
435: PetscCall(PetscViewerASCIIPushTab(viewer));
436: PetscCall(TaoTermView(subterm, viewer));
437: PetscCall(PetscViewerASCIIPopTab(viewer));
438: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL && map == NULL) PetscCall(PetscViewerASCIIPrintf(viewer, "Map: unmapped\n"));
439: else if (map != NULL) {
440: PetscCall(PetscViewerASCIIPrintf(viewer, "Map:\n"));
441: PetscCall(PetscViewerASCIIPushTab(viewer));
442: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
443: PetscCall(MatView(map, viewer));
444: PetscCall(PetscViewerPopFormat(viewer));
445: PetscCall(PetscViewerASCIIPopTab(viewer));
446: }
447: PetscCall(TaoTermSumGetTermMask(term, i, &mask));
448: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL && mask != TAOTERM_MASK_NONE) {
449: PetscBool preceding = PETSC_FALSE;
451: PetscCall(PetscViewerASCIIPrintf(viewer, "Mask (tao_term_sum_%smask): ", subprefix));
452: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
453: if (TaoTermObjectiveMasked(mask)) {
454: PetscCall(PetscViewerASCIIPrintf(viewer, "objective"));
455: preceding = PETSC_TRUE;
456: }
457: if (TaoTermGradientMasked(mask)) {
458: PetscCall(PetscViewerASCIIPrintf(viewer, "%sgradient", preceding ? ", " : ""));
459: preceding = PETSC_TRUE;
460: }
461: if (TaoTermHessianMasked(mask)) PetscCall(PetscViewerASCIIPrintf(viewer, "%shessian", preceding ? ", " : ""));
462: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
463: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
464: }
466: PetscCall(PetscViewerASCIIPopTab(viewer));
467: }
468: PetscCall(PetscViewerASCIIPopTab(viewer));
469: }
470: PetscFunctionReturn(PETSC_SUCCESS);
471: }
473: /*@
474: TaoTermSumSetNumberTerms - Set the number of terms in the sum
476: Collective
478: Input Parameters:
479: + term - a `TaoTerm` of type `TAOTERMSUM`
480: - n_terms - the number of terms that will be in the sum
482: Level: developer
484: Note:
485: If `n_terms` is smaller than the current number of terms, the trailing terms will be dropped.
487: .seealso: [](sec_tao_term),
488: `TaoTerm`,
489: `TAOTERMSUM`,
490: `TaoTermSumGetNumberTerms()`
491: @*/
492: PetscErrorCode TaoTermSumSetNumberTerms(TaoTerm term, PetscInt n_terms)
493: {
494: PetscFunctionBegin;
497: PetscTryMethod(term, "TaoTermSumSetNumberTerms_C", (TaoTerm, PetscInt), (term, n_terms));
498: PetscFunctionReturn(PETSC_SUCCESS);
499: }
501: static PetscErrorCode TaoTermSumSetNumberTerms_Sum(TaoTerm term, PetscInt n_terms)
502: {
503: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
504: PetscInt n_terms_old = sum->n_terms;
505: PetscReal *new_values;
506: TaoTermMapping *new_summands;
508: PetscFunctionBegin;
509: if (n_terms == n_terms_old) PetscFunctionReturn(PETSC_SUCCESS);
510: for (PetscInt i = n_terms; i < n_terms_old; i++) PetscCall(TaoTermMappingReset(&sum->terms[i]));
511: PetscCall(PetscMalloc1(n_terms, &new_summands));
512: PetscCall(PetscCalloc1(n_terms, &new_values));
513: PetscCall(PetscArraycpy(new_summands, sum->terms, PetscMin(n_terms, n_terms_old)));
514: PetscCall(PetscArrayzero(&new_summands[n_terms_old], PetscMax(0, n_terms - n_terms_old)));
515: PetscCall(PetscFree(sum->terms));
516: PetscCall(PetscFree(sum->subterm_values));
517: sum->terms = new_summands;
518: sum->subterm_values = new_values;
519: sum->n_terms = n_terms;
520: for (PetscInt i = n_terms_old; i < n_terms; i++) PetscCall(TaoTermSumSetTerm(term, i, NULL, 1.0, NULL, NULL));
521: PetscFunctionReturn(PETSC_SUCCESS);
522: }
524: /*@
525: TaoTermSumGetNumberTerms - Get the number of terms in the sum
527: Not collective
529: Input Parameter:
530: . term - a `TaoTerm` of type `TAOTERMSUM`
532: Output Parameter:
533: . n_terms - the number of terms that will be in the sum
535: Level: developer
537: .seealso: [](sec_tao_term),
538: `TaoTerm`,
539: `TAOTERMSUM`,
540: `TaoTermSumSetNumberTerms()`
541: @*/
542: PetscErrorCode TaoTermSumGetNumberTerms(TaoTerm term, PetscInt *n_terms)
543: {
544: PetscFunctionBegin;
546: PetscAssertPointer(n_terms, 2);
547: PetscUseMethod(term, "TaoTermSumGetNumberTerms_C", (TaoTerm, PetscInt *), (term, n_terms));
548: PetscFunctionReturn(PETSC_SUCCESS);
549: }
551: static PetscErrorCode TaoTermSumGetNumberTerms_Sum(TaoTerm term, PetscInt *n_terms)
552: {
553: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
555: PetscFunctionBegin;
556: *n_terms = sum->n_terms;
557: PetscFunctionReturn(PETSC_SUCCESS);
558: }
560: /*@
561: TaoTermSumGetTerm - Get the data for a term in a `TAOTERMSUM`
563: Not collective
565: Input Parameters:
566: + sumterm - a `TaoTerm` of type `TAOTERMSUM`
567: - index - a number $0 \leq i < n$, where $n$ is the number of terms in `TaoTermSumGetNumberTerms()`
569: Output Parameters:
570: + prefix - (optional) the prefix used for configuring the term
571: . scale - (optional) the coefficient scaling the term in the sum
572: . term - the `TaoTerm` at given index of `TAOTERMSUM`
573: - map - (optional) a map from the `TAOTERMSUM` solution space to the `term` solution space; if `NULL` the map is assumed to be the identity
575: Level: developer
577: .seealso: [](sec_tao_term),
578: `TaoTerm`,
579: `TAOTERMSUM`,
580: `TaoTermSumSetTerm()`,
581: `TaoTermSumAddTerm()`
582: @*/
583: PetscErrorCode TaoTermSumGetTerm(TaoTerm sumterm, PetscInt index, const char *prefix[], PetscReal *scale, TaoTerm *term, Mat *map)
584: {
585: PetscFunctionBegin;
588: if (prefix) PetscAssertPointer(prefix, 3);
589: if (term) PetscAssertPointer(term, 5);
590: if (scale) PetscAssertPointer(scale, 4);
591: if (map) PetscAssertPointer(map, 6);
592: PetscUseMethod(sumterm, "TaoTermSumGetTerm_C", (TaoTerm, PetscInt, const char **, PetscReal *, TaoTerm *, Mat *), (sumterm, index, prefix, scale, term, map));
593: PetscFunctionReturn(PETSC_SUCCESS);
594: }
596: static PetscErrorCode TaoTermSumGetTerm_Sum(TaoTerm term, PetscInt index, const char **prefix, PetscReal *scale, TaoTerm *subterm, Mat *map)
597: {
598: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
599: TaoTermMapping *summand;
601: PetscFunctionBegin;
602: PetscCheck(index >= 0 && index < sum->n_terms, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Index %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", index, sum->n_terms);
603: summand = &sum->terms[index];
604: PetscCall(TaoTermMappingGetData(summand, prefix, scale, subterm, map));
605: PetscFunctionReturn(PETSC_SUCCESS);
606: }
608: /*@
609: TaoTermSumSetTerm - Set a term in a sum of terms
611: Collective
613: Input Parameters:
614: + sumterm - a `TaoTerm` of type `TAOTERMSUM`
615: . index - a number $0 \leq i < n$, where $n$ is the number of terms in `TaoTermSumSetNumberTerms()`
616: . prefix - (optional) the prefix used for configuring the term (if `NULL`, `term_x_` will be the prefix, e.g. "term_0_", "term_1_", etc.)
617: . scale - the coefficient scaling the term in the sum
618: . term - the `TaoTerm` to be set in `TAOTERMSUM`
619: - map - (optional) a map from the `TAOTERMSUM` solution space to the `term` solution space; if `NULL` the map is assumed to be the identity
621: Level: developer
623: .seealso: [](sec_tao_term),
624: `TaoTerm`,
625: `TAOTERMSUM`,
626: `TaoTermSumGetTerm()`,
627: `TaoTermSumAddTerm()`
628: @*/
629: PetscErrorCode TaoTermSumSetTerm(TaoTerm sumterm, PetscInt index, const char prefix[], PetscReal scale, TaoTerm term, Mat map)
630: {
631: PetscFunctionBegin;
634: if (prefix) PetscAssertPointer(prefix, 3);
638: PetscTryMethod(sumterm, "TaoTermSumSetTerm_C", (TaoTerm, PetscInt, const char[], PetscReal, TaoTerm, Mat), (sumterm, index, prefix, scale, term, map));
639: PetscFunctionReturn(PETSC_SUCCESS);
640: }
642: static PetscErrorCode TaoTermSumSetTerm_Sum(TaoTerm term, PetscInt index, const char prefix[], PetscReal scale, TaoTerm subterm, Mat map)
643: {
644: char subterm_x_[256];
645: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
646: TaoTermMapping *summand;
648: PetscFunctionBegin;
649: PetscCheck(index >= 0 && index < sum->n_terms, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Index %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", index, sum->n_terms);
650: summand = &sum->terms[index];
651: if (prefix == NULL) {
652: PetscCall(PetscSNPrintf(subterm_x_, 256, "term_%" PetscInt_FMT "_", index));
653: prefix = subterm_x_;
654: }
655: PetscCall(TaoTermMappingSetData(summand, prefix, scale, subterm, map));
656: PetscFunctionReturn(PETSC_SUCCESS);
657: }
659: /*@
660: TaoTermSumSetTermHessianMatrices - Set Hessian matrices that can be used internally by a `TAOTERMSUM`
662: Logically collective
664: Input Parameters:
665: + term - a `TaoTerm` of type `TAOTERMSUM`
666: . index - the index for the term from `TaoTermSumSetTerm()` or `TaoTermSumAddTerm()`
667: . unmapped_H - (optional) unmapped Hessian matrix
668: . unmapped_Hpre - (optional) unmapped matrix for constructing the preconditioner of `unmapped_H`
669: . mapped_H - (optional) Hessian matrix
670: - mapped_Hpre - (optional) matrix for constructing the preconditioner of `mapped_H`
672: Level: developer
674: Notes:
675: If the inner term has the form $g(x) = \alpha f(Ax; p)$, the "mapped" Hessians should be able to hold the Hessian
676: $\nabla^2 g$ and the unmapped Hessians should be able to hold the Hessian $\nabla_x^2 f$. If the term is not mapped,
677: just pass the unmapped Hessians (e.g. `TaoTermSumSetTermHessianMatrices(term, 0, H, Hpre, NULL, NULL)`).
679: .seealso: [](sec_tao_term),
680: `TaoTerm`,
681: `TAOTERMSUM`,
682: `TaoTermComputeHessian()`,
683: `TaoTermSumGetTermHessianMatrices()`
684: @*/
685: PetscErrorCode TaoTermSumSetTermHessianMatrices(TaoTerm term, PetscInt index, Mat unmapped_H, Mat unmapped_Hpre, Mat mapped_H, Mat mapped_Hpre)
686: {
687: PetscFunctionBegin;
694: PetscTryMethod(term, "TaoTermSumSetTermHessianMatrices_C", (TaoTerm, PetscInt, Mat, Mat, Mat, Mat), (term, index, unmapped_H, unmapped_Hpre, mapped_H, mapped_Hpre));
695: PetscFunctionReturn(PETSC_SUCCESS);
696: }
698: static PetscErrorCode TaoTermSumSetTermHessianMatrices_Sum(TaoTerm term, PetscInt index, Mat unmapped_H, Mat unmapped_Hpre, Mat mapped_H, Mat mapped_Hpre)
699: {
700: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
701: TaoTermMapping *summand;
702: PetscBool is_callback;
704: PetscFunctionBegin;
705: PetscCheck(index >= 0 && index < sum->n_terms, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Index %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", index, sum->n_terms);
706: summand = &sum->terms[index];
708: PetscCall(PetscObjectTypeCompare((PetscObject)summand->term, TAOTERMCALLBACKS, &is_callback));
709: if (is_callback) {
710: MatType H_type, Hpre_type;
711: PetscBool Hpre_is_H;
713: Hpre_is_H = (mapped_H == mapped_Hpre) ? PETSC_TRUE : PETSC_FALSE;
715: if (mapped_H) PetscCall(MatGetType(mapped_H, &H_type));
716: else H_type = NULL;
717: if (mapped_Hpre) PetscCall(MatGetType(mapped_Hpre, &Hpre_type));
718: else Hpre_type = NULL;
719: PetscCall(TaoTermSetCreateHessianMode(summand->term, Hpre_is_H, H_type, Hpre_type));
720: }
722: PetscCall(PetscObjectReference((PetscObject)unmapped_H));
723: PetscCall(MatDestroy(&summand->_unmapped_H));
724: summand->_unmapped_H = unmapped_H;
726: PetscCall(PetscObjectReference((PetscObject)unmapped_Hpre));
727: PetscCall(MatDestroy(&summand->_unmapped_Hpre));
728: summand->_unmapped_Hpre = unmapped_Hpre;
730: PetscCall(PetscObjectReference((PetscObject)mapped_H));
731: PetscCall(MatDestroy(&summand->_mapped_H));
732: summand->_mapped_H = mapped_H;
734: PetscCall(PetscObjectReference((PetscObject)mapped_Hpre));
735: PetscCall(MatDestroy(&summand->_mapped_Hpre));
736: summand->_mapped_Hpre = mapped_Hpre;
737: PetscFunctionReturn(PETSC_SUCCESS);
738: }
740: /*@
741: TaoTermSumGetTermHessianMatrices - Get Hessian matrices set with `TaoTermSumSetTermHessianMatrices()`.
743: Not collective
745: Input Parameters:
746: + term - a `TaoTerm` of type `TAOTERMSUM`
747: - index - the index for the term from `TaoTermSumSetTerm()` or `TaoTermSumAddTerm()`
749: Output Parameters:
750: + unmapped_H - (optional) unmapped Hessian matrix
751: . unmapped_Hpre - (optional) unmapped matrix for constructing the preconditioner for `unmapped_H`
752: . mapped_H - (optional) Hessian matrix
753: - mapped_Hpre - (optional) matrix for constructing the preconditioner for `mapped_H`
755: Level: developer
757: .seealso: [](sec_tao_term),
758: `TaoTerm`,
759: `TAOTERMSUM`,
760: `TaoTermComputeHessian()`,
761: `TaoTermSumSetTermHessianMatrices()`
762: @*/
763: PetscErrorCode TaoTermSumGetTermHessianMatrices(TaoTerm term, PetscInt index, Mat *unmapped_H, Mat *unmapped_Hpre, Mat *mapped_H, Mat *mapped_Hpre)
764: {
765: PetscFunctionBegin;
767: if (unmapped_H) PetscAssertPointer(unmapped_H, 3);
768: if (unmapped_Hpre) PetscAssertPointer(unmapped_Hpre, 4);
769: if (mapped_H) PetscAssertPointer(mapped_H, 5);
770: if (mapped_Hpre) PetscAssertPointer(mapped_Hpre, 6);
771: PetscTryMethod(term, "TaoTermSumGetTermHessianMatrices_C", (TaoTerm, PetscInt, Mat *, Mat *, Mat *, Mat *), (term, index, unmapped_H, unmapped_Hpre, mapped_H, mapped_Hpre));
772: PetscFunctionReturn(PETSC_SUCCESS);
773: }
775: static PetscErrorCode TaoTermSumGetTermHessianMatrices_Sum(TaoTerm term, PetscInt index, Mat *unmapped_H, Mat *unmapped_Hpre, Mat *mapped_H, Mat *mapped_Hpre)
776: {
777: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
778: TaoTermMapping *summand;
780: PetscFunctionBegin;
781: PetscCheck(index >= 0 && index < sum->n_terms, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Index %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", index, sum->n_terms);
782: summand = &sum->terms[index];
784: if (unmapped_H) *unmapped_H = summand->_unmapped_H;
785: if (unmapped_Hpre) *unmapped_Hpre = summand->_unmapped_Hpre;
786: if (mapped_H) *mapped_H = summand->_mapped_H;
787: if (mapped_Hpre) *mapped_Hpre = summand->_mapped_Hpre;
788: PetscFunctionReturn(PETSC_SUCCESS);
789: }
791: /*@
792: TaoTermSumGetTermMask - Get the `TaoTermMask` of a term in the sum
794: Not collective
796: Input Parameters:
797: + term - a `TaoTerm` of type `TAOTERMSUM`
798: - index - the index for the term from `TaoTermSumSetTerm()` or `TaoTermSumAddTerm()`
800: Output Parameter:
801: . mask - a bitmask of `TaoTermMask` evaluation methods to mask (e.g. just `TAOTERM_MASK_OBJECTIVE` or a bitwise-or like `TAOTERM_MASK_OBJECTIVE | TAOTERM_MASK_GRADIENT`)
803: Level: developer
805: .seealso: [](sec_tao_term),
806: `TaoTerm`,
807: `TAOTERMSUM`,
808: `TaoTermSumSetTermMask()`
809: @*/
810: PetscErrorCode TaoTermSumGetTermMask(TaoTerm term, PetscInt index, TaoTermMask *mask)
811: {
812: PetscFunctionBegin;
814: PetscAssertPointer(mask, 3);
815: PetscUseMethod(term, "TaoTermSumGetTermMask_C", (TaoTerm, PetscInt, TaoTermMask *), (term, index, mask));
816: PetscFunctionReturn(PETSC_SUCCESS);
817: }
819: static PetscErrorCode TaoTermSumGetTermMask_Sum(TaoTerm term, PetscInt index, TaoTermMask *mask)
820: {
821: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
822: TaoTermMapping *summand;
824: PetscFunctionBegin;
825: PetscCheck(index >= 0 && index < sum->n_terms, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Index %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", index, sum->n_terms);
826: summand = &sum->terms[index];
827: *mask = summand->mask;
828: PetscFunctionReturn(PETSC_SUCCESS);
829: }
831: /*@
832: TaoTermSumSetTermMask - Set a `TaoTermMask` on a term in the sum
834: Logically collective
836: Input Parameters:
837: + term - a `TaoTerm` of type `TAOTERMSUM`
838: . index - the index for the term from `TaoTermSumSetTerm()` or `TaoTermSumAddTerm()`
839: - mask - a bitmask of `TaoTermMask` evaluation methods to mask (e.g. just `TAOTERM_MASK_OBJECTIVE` or a bitwise-or like `TAOTERM_MASK_OBJECTIVE | TAOTERM_MASK_GRADIENT`)
841: Options Database Keys:
842: . -tao_term_sum_<prefix_>mask - a list containing any of `none`, `objective`, `gradient`, and `hessian` to indicate which evaluations to mask for a term with a given prefix (see `TaoTermSumSetTerm()`)
844: Level: developer
846: Note:
847: Some optimization methods may add a damping term to the Hessian of an
848: objective function without affecting the objective or gradient. If, e.g.,
849: the regularizer has index `1`, then this can be accomplished with
850: `TaoTermSumSetTermMask(term, 1, TAOTERM_MASK_OBJECTIVE | TAOTERM_MASK_GRADIENT)`.
852: .seealso: [](sec_tao_term),
853: `TaoTerm`,
854: `TAOTERMSUM`,
855: `TaoTermSumGetTermMask()`
856: @*/
857: PetscErrorCode TaoTermSumSetTermMask(TaoTerm term, PetscInt index, TaoTermMask mask)
858: {
859: PetscFunctionBegin;
862: PetscTryMethod(term, "TaoTermSumSetTermMask_C", (TaoTerm, PetscInt, TaoTermMask), (term, index, mask));
863: PetscFunctionReturn(PETSC_SUCCESS);
864: }
866: static PetscErrorCode TaoTermSumSetTermMask_Sum(TaoTerm term, PetscInt index, TaoTermMask mask)
867: {
868: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
869: TaoTermMapping *summand;
871: PetscFunctionBegin;
872: PetscCheck(index >= 0 && index < sum->n_terms, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Index %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", index, sum->n_terms);
873: summand = &sum->terms[index];
874: summand->mask = mask;
875: PetscFunctionReturn(PETSC_SUCCESS);
876: }
878: /*@
879: TaoTermSumAddTerm - Append a term to the terms being summed
881: Collective
883: Input Parameters:
884: + sumterm - a `TaoTerm` of type `TAOTERMSUM`
885: . prefix - (optional) the prefix used for configuring the term (if `NULL`, the index of the term will be used as a prefix, e.g. `term_0_`, `term_1_`, etc.)
886: . scale - the coefficient scaling the term in the sum
887: . term - the `TaoTerm` to add
888: - map - (optional) a map from the `TAOTERMSUM` solution space to the `term` solution space; if `NULL` the map is assumed to be the identity
890: Output Parameter:
891: . index - (optional) the index of the newly added term
893: Level: developer
895: .seealso: [](sec_tao_term), `TaoTerm`, `TAOTERMSUM`
896: @*/
897: PetscErrorCode TaoTermSumAddTerm(TaoTerm sumterm, const char prefix[], PetscReal scale, TaoTerm term, Mat map, PetscInt *index)
898: {
899: PetscFunctionBegin;
901: if (prefix) PetscAssertPointer(prefix, 2);
905: if (index) PetscAssertPointer(index, 6);
906: PetscTryMethod(sumterm, "TaoTermSumAddTerm_C", (TaoTerm, const char[], PetscReal, TaoTerm, Mat, PetscInt *), (sumterm, prefix, scale, term, map, index));
907: PetscFunctionReturn(PETSC_SUCCESS);
908: }
910: static PetscErrorCode TaoTermSumAddTerm_Sum(TaoTerm term, const char prefix[], PetscReal scale, TaoTerm subterm, Mat map, PetscInt *index)
911: {
912: PetscInt n_terms_old;
914: PetscFunctionBegin;
915: PetscCall(TaoTermSumGetNumberTerms(term, &n_terms_old));
916: PetscCall(TaoTermSumSetNumberTerms(term, n_terms_old + 1));
917: PetscCall(TaoTermSumSetTerm(term, n_terms_old, prefix, scale, subterm, map));
918: if (index) *index = n_terms_old;
919: PetscFunctionReturn(PETSC_SUCCESS);
920: }
922: static PetscErrorCode TaoTermSetFromOptions_Sum(TaoTerm term, PetscOptionItems PetscOptionsObject)
923: {
924: PetscInt n_terms;
925: const char *prefix;
927: PetscFunctionBegin;
928: PetscCall(TaoTermSumGetNumberTerms(term, &n_terms));
929: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)term, &prefix));
930: PetscOptionsHeadBegin(PetscOptionsObject, "TaoTerm sum options");
931: PetscCall(PetscOptionsBoundedInt("-tao_term_sum_number_terms", "The number of terms in the sum", "TaoTermSumSetNumberTerms", n_terms, &n_terms, NULL, 0));
932: PetscCall(TaoTermSumSetNumberTerms(term, n_terms));
933: for (PetscInt i = 0; i < n_terms; i++) {
934: const char *subprefix;
935: Mat map;
936: PetscReal scale;
937: TaoTerm subterm;
938: char arg[256];
939: PetscBool flg;
940: PetscEnum masks[4] = {ENUM_DUMMY, ENUM_DUMMY, ENUM_DUMMY, ENUM_DUMMY};
941: PetscInt n_masks = PETSC_STATIC_ARRAY_LENGTH(masks);
943: PetscCall(TaoTermSumGetTerm(term, i, &subprefix, &scale, &subterm, &map));
944: if (subterm == NULL) {
945: PetscCall(TaoTermDuplicate(term, TAOTERM_DUPLICATE_SIZEONLY, &subterm));
946: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)subterm, prefix));
947: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)subterm, subprefix));
948: } else PetscCall(PetscObjectReference((PetscObject)subterm));
949: PetscCall(TaoTermSetFromOptions(subterm));
951: PetscCall(PetscSNPrintf(arg, 256, "-tao_term_sum_%sscale", subprefix));
952: PetscCall(PetscOptionsReal(arg, "The scale of the term in the TaoTermSum", "TaoTermSumSetTerm", scale, &scale, NULL));
954: PetscCall(PetscSNPrintf(arg, 256, "-tao_term_sum_%smask", subprefix));
955: PetscCall(PetscOptionsEnumArray(arg, "The mask of the term in the TaoTermSum", "TaoTermSumSetTermMask", TaoTermMasks, masks, &n_masks, &flg));
956: if (flg) {
957: PetscEnum mask = (PetscEnum)TAOTERM_MASK_NONE;
959: for (PetscInt j = 0; j < n_masks; j++) {
960: PetscEnum this_mask = masks[j] ? (PetscEnum)(1 << (masks[j] - 1)) : (PetscEnum)TAOTERM_MASK_NONE;
962: mask = (PetscEnum)(mask | this_mask);
963: }
965: PetscCall(TaoTermSumSetTermMask(term, i, (TaoTermMask)mask));
966: }
968: if (map) PetscCall(MatSetFromOptions(map));
969: PetscCall(TaoTermSumSetTerm(term, i, subprefix, scale, subterm, map));
970: PetscCall(TaoTermDestroy(&subterm));
971: }
972: PetscOptionsHeadEnd();
973: PetscFunctionReturn(PETSC_SUCCESS);
974: }
976: /*@
977: TaoTermSumGetLastTermObjectives - Get the contributions from each term to the
978: last evaluation of `TaoTermComputeObjective()` or `TaoTermComputeObjectiveAndGradient()`
980: Not collective
982: Input Parameter:
983: . term - a `TaoTerm` of type `TAOTERMSUM`
985: Output Parameter:
986: . values - an array of the contributions to the last computed objective value
988: Level: developer
990: .seealso: [](sec_tao_term),
991: `TaoTerm`,
992: `TAOTERMSUM`
993: @*/
994: PetscErrorCode TaoTermSumGetLastTermObjectives(TaoTerm term, const PetscReal *values[])
995: {
996: PetscFunctionBegin;
998: PetscAssertPointer(values, 2);
999: PetscUseMethod(term, "TaoTermSumGetLastTermObjectives_C", (TaoTerm, const PetscReal *[]), (term, values));
1000: PetscFunctionReturn(PETSC_SUCCESS);
1001: }
1003: static PetscErrorCode TaoTermSumGetLastTermObjectives_Sum(TaoTerm term, const PetscReal *values[])
1004: {
1005: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
1007: PetscFunctionBegin;
1008: *values = sum->subterm_values;
1009: PetscFunctionReturn(PETSC_SUCCESS);
1010: }
1012: static PetscErrorCode TaoTermComputeObjective_Sum(TaoTerm term, Vec x, Vec params, PetscReal *value)
1013: {
1014: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
1015: Vec *sub_params = NULL;
1016: PetscBool *is_dummy = NULL;
1017: PetscReal value_;
1018: PetscReal *values = sum->subterm_values;
1020: PetscFunctionBegin;
1021: if (params) PetscCall(TaoTermSumVecNestGetSubVecsRead(params, NULL, &sub_params, &is_dummy));
1022: value_ = 0.0;
1023: for (PetscInt i = 0; i < sum->n_terms; i++) {
1024: TaoTermMapping *summand = &sum->terms[i];
1025: Vec sub_param = TaoTermSumGetSubVec(params, sub_params, is_dummy, i);
1027: PetscCall(TaoTermMappingComputeObjective(summand, x, sub_param, INSERT_VALUES, &values[i]));
1028: value_ += values[i];
1029: }
1030: if (params) PetscCall(TaoTermSumVecNestRestoreSubVecsRead(params, NULL, &sub_params, &is_dummy));
1031: *value = value_;
1032: PetscFunctionReturn(PETSC_SUCCESS);
1033: }
1035: static PetscErrorCode TaoTermComputeGradient_Sum(TaoTerm term, Vec x, Vec params, Vec g)
1036: {
1037: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
1038: Vec *sub_params = NULL;
1039: PetscBool *is_dummy = NULL;
1041: PetscFunctionBegin;
1042: if (params) PetscCall(TaoTermSumVecNestGetSubVecsRead(params, NULL, &sub_params, &is_dummy));
1043: for (PetscInt i = 0; i < sum->n_terms; i++) {
1044: TaoTermMapping *summand = &sum->terms[i];
1045: Vec sub_param = TaoTermSumGetSubVec(params, sub_params, is_dummy, i);
1047: PetscCall(TaoTermMappingComputeGradient(summand, x, sub_param, i == 0 ? INSERT_VALUES : ADD_VALUES, g));
1048: }
1049: if (params) PetscCall(TaoTermSumVecNestRestoreSubVecsRead(params, NULL, &sub_params, &is_dummy));
1050: PetscFunctionReturn(PETSC_SUCCESS);
1051: }
1053: static PetscErrorCode TaoTermComputeObjectiveAndGradient_Sum(TaoTerm term, Vec x, Vec params, PetscReal *value, Vec g)
1054: {
1055: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
1056: Vec *sub_params = NULL;
1057: PetscBool *is_dummy = NULL;
1058: PetscReal *values = sum->subterm_values;
1059: PetscReal value_;
1061: PetscFunctionBegin;
1062: if (params) PetscCall(TaoTermSumVecNestGetSubVecsRead(params, NULL, &sub_params, &is_dummy));
1063: value_ = 0.0;
1064: for (PetscInt i = 0; i < sum->n_terms; i++) {
1065: TaoTermMapping *summand = &sum->terms[i];
1066: Vec sub_param = TaoTermSumGetSubVec(params, sub_params, is_dummy, i);
1068: values[i] = 0.0;
1069: PetscCall(TaoTermMappingComputeObjectiveAndGradient(summand, x, sub_param, i == 0 ? INSERT_VALUES : ADD_VALUES, &values[i], g));
1070: value_ += values[i];
1071: }
1072: if (params) PetscCall(TaoTermSumVecNestRestoreSubVecsRead(params, NULL, &sub_params, &is_dummy));
1073: *value = value_;
1074: PetscFunctionReturn(PETSC_SUCCESS);
1075: }
1077: static PetscErrorCode TaoTermComputeHessian_Sum(TaoTerm term, Vec x, Vec params, Mat H, Mat Hpre)
1078: {
1079: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
1080: Vec *sub_params = NULL;
1081: PetscBool *is_dummy = NULL;
1083: PetscFunctionBegin;
1084: if (H == NULL && Hpre == NULL) PetscFunctionReturn(PETSC_SUCCESS);
1085: if (params) PetscCall(TaoTermSumVecNestGetSubVecsRead(params, NULL, &sub_params, &is_dummy));
1086: // If mattype dense, then after zero entries, H->assembled = true.
1087: // But for aij, H->assembled is still false.
1088: if (H) {
1089: PetscCall(MatZeroEntries(H));
1090: PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
1091: PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
1092: }
1093: if (Hpre && (Hpre != H)) {
1094: PetscCall(MatZeroEntries(Hpre));
1095: PetscCall(MatAssemblyBegin(Hpre, MAT_FINAL_ASSEMBLY));
1096: PetscCall(MatAssemblyEnd(Hpre, MAT_FINAL_ASSEMBLY));
1097: }
1098: for (PetscInt i = 0; i < sum->n_terms; i++) {
1099: TaoTermMapping *summand = &sum->terms[i];
1100: Vec sub_param = TaoTermSumGetSubVec(params, sub_params, is_dummy, i);
1102: PetscCall(TaoTermMappingComputeHessian(summand, x, sub_param, ADD_VALUES, H, Hpre == H ? NULL : Hpre));
1103: }
1104: if (params) PetscCall(TaoTermSumVecNestRestoreSubVecsRead(params, NULL, &sub_params, &is_dummy));
1105: PetscFunctionReturn(PETSC_SUCCESS);
1106: }
1108: static PetscErrorCode TaoTermSetUp_Sum(TaoTerm term)
1109: {
1110: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
1111: PetscBool all_none = PETSC_TRUE;
1112: PetscBool any_required = PETSC_FALSE;
1113: PetscInt k = 0, K = 0;
1114: Mat *mats, new_parameters_factory;
1115: PetscLayout layout = NULL, clayout;
1117: PetscFunctionBegin;
1118: PetscCall(PetscCalloc1(sum->n_terms, &mats));
1119: PetscCall(MatGetLayouts(term->solution_factory, &layout, &clayout));
1120: if (layout->setupcalled == PETSC_FALSE) layout = NULL;
1121: for (PetscInt i = 0; i < sum->n_terms; i++) {
1122: TaoTermMapping *summand = &sum->terms[i];
1123: TaoTermParametersMode submode;
1124: PetscLayout sub_layout;
1125: PetscBool congruent;
1127: PetscCall(TaoTermSetUp(summand->term));
1128: if (summand->map) {
1129: PetscCall(MatSetUp(summand->map));
1130: PetscCall(MatGetLayouts(summand->map, NULL, &sub_layout));
1131: } else PetscCall(TaoTermGetSolutionLayout(summand->term, &sub_layout));
1132: if (i == 0 && layout == NULL) layout = sub_layout;
1133: PetscCall(PetscLayoutCompare(layout, sub_layout, &congruent));
1134: if (congruent == PETSC_FALSE) {
1135: PetscInt N, sub_N;
1136: MPI_Comm comm = PetscObjectComm((PetscObject)term);
1138: PetscCall(PetscLayoutGetSize(layout, &N));
1139: PetscCall(PetscLayoutGetSize(sub_layout, &sub_N));
1141: SETERRQ(comm, PETSC_ERR_ARG_SIZ, "%sterm %" PetscInt_FMT " has solution layout (input size %" PetscInt_FMT ") that is incompatible with %s solution layout (size %" PetscInt_FMT ")", summand->map ? "mapped " : "", i, sub_N, i == 0 ? "the sum's" : "previous terms'", N);
1142: }
1143: PetscCall(TaoTermGetParametersMode(summand->term, &submode));
1144: if (submode == TAOTERM_PARAMETERS_REQUIRED) any_required = PETSC_TRUE;
1145: if (submode != TAOTERM_PARAMETERS_NONE) {
1146: PetscInt subk, subK;
1148: all_none = PETSC_FALSE;
1149: PetscCall(TaoTermGetParametersSizes(summand->term, &subk, &subK, NULL));
1150: k += subk;
1151: K += subK;
1152: }
1153: if (summand->term->parameters_mode != TAOTERM_PARAMETERS_NONE) {
1154: PetscCall(PetscObjectReference((PetscObject)summand->term->parameters_factory));
1155: mats[i] = summand->term->parameters_factory;
1156: } else {
1157: PetscCall(MatCreate(PetscObjectComm((PetscObject)term), &mats[i]));
1158: PetscCall(MatSetType(mats[i], MATDUMMY));
1159: PetscCall(MatSetSizes(mats[i], 0, 0, 0, 0));
1160: PetscCall(MatSetUp(mats[i]));
1161: }
1162: }
1163: PetscCall(MatSetLayouts(term->solution_factory, layout, clayout));
1164: if (all_none) {
1165: term->parameters_mode = TAOTERM_PARAMETERS_NONE;
1166: } else if (any_required) {
1167: term->parameters_mode = TAOTERM_PARAMETERS_REQUIRED;
1168: } else {
1169: term->parameters_mode = TAOTERM_PARAMETERS_OPTIONAL;
1170: }
1171: PetscCall(TaoTermSetParametersSizes(term, k, K, 1));
1172: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)term), sum->n_terms, NULL, 1, NULL, mats, &new_parameters_factory));
1173: PetscCall(MatSetVecType(new_parameters_factory, VECNEST));
1174: PetscCall(MatDestroy(&term->parameters_factory));
1175: term->parameters_factory = new_parameters_factory;
1176: for (PetscInt i = 0; i < sum->n_terms; i++) PetscCall(MatDestroy(&mats[i]));
1177: PetscCall(PetscFree(mats));
1178: PetscFunctionReturn(PETSC_SUCCESS);
1179: }
1181: static PetscErrorCode TaoTermCreateSolutionVec_Sum(TaoTerm term, Vec *solution_vec)
1182: {
1183: PetscFunctionBegin;
1184: if (solution_vec) PetscCall(MatCreateVecs(term->solution_factory, NULL, solution_vec));
1185: PetscFunctionReturn(PETSC_SUCCESS);
1186: }
1188: static PetscErrorCode TaoTermCreateParametersVec_Sum(TaoTerm term, Vec *parameters_vec)
1189: {
1190: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
1192: PetscFunctionBegin;
1193: if (parameters_vec) {
1194: Vec *vecs;
1196: PetscCall(PetscCalloc1(sum->n_terms, &vecs));
1197: for (PetscInt i = 0; i < sum->n_terms; i++) {
1198: TaoTermMapping *summand = &sum->terms[i];
1199: TaoTermParametersMode submode;
1201: PetscCall(TaoTermGetParametersMode(summand->term, &submode));
1202: if (submode != TAOTERM_PARAMETERS_NONE) PetscCall(TaoTermCreateParametersVec(summand->term, &vecs[i]));
1203: }
1204: PetscCall(TaoTermSumParametersPack(term, vecs, parameters_vec));
1205: for (PetscInt i = 0; i < sum->n_terms; i++) PetscCall(VecDestroy(&vecs[i]));
1206: PetscCall(PetscFree(vecs));
1207: }
1208: PetscFunctionReturn(PETSC_SUCCESS);
1209: }
1211: static PetscErrorCode TaoTermCreateHessianMatrices_Sum(TaoTerm term, Mat *H, Mat *Hpre)
1212: {
1213: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
1214: PetscBool Hpre_is_H, sub_Hpre_is_H;
1216: PetscFunctionBegin;
1217: Hpre_is_H = term->Hpre_is_H;
1218: // Need to create subterms' mapped Hessians and PtAP routines, if needed
1219: for (PetscInt i = 0; i < sum->n_terms; i++) {
1220: TaoTermMapping *summand = &sum->terms[i];
1221: PetscBool is_callback;
1223: PetscCall(PetscObjectTypeCompare((PetscObject)summand->term, TAOTERMCALLBACKS, &is_callback));
1224: if (is_callback) {
1225: Mat c_H;
1227: PetscCall(TaoTermSumGetTermHessianMatrices(term, i, NULL, NULL, &c_H, NULL));
1228: PetscCheck(c_H, PetscObjectComm((PetscObject)summand->term), PETSC_ERR_USER, "TAOTERMCALLBACKS does not have Hessian routines set. Call TaoSetHessian()");
1229: }
1230: PetscCall(TaoTermMappingCreateHessianMatrices(summand, &summand->_mapped_H, &summand->_mapped_Hpre));
1232: sub_Hpre_is_H = (summand->_mapped_H == summand->_mapped_Hpre) ? PETSC_TRUE : PETSC_FALSE;
1233: Hpre_is_H = (Hpre_is_H && sub_Hpre_is_H) ? PETSC_TRUE : PETSC_FALSE;
1234: }
1236: term->Hpre_is_H = Hpre_is_H;
1237: PetscCall(TaoTermCreateHessianMatricesDefault(term, H, Hpre));
1238: PetscFunctionReturn(PETSC_SUCCESS);
1239: }
1241: /*MC
1242: TAOTERMSUM - A `TaoTerm` that is a sum of multiple `TaoTerms`.
1244: Level: developer
1246: Note:
1247: The default Hessian creation mode (see `TaoTermGetCreateHessianMode()`) is `H == Hpre` and `TaoTermCreateHessianMatrices()`
1248: will create a `MATAIJ`.
1250: .seealso: [](sec_tao_term),
1251: `TaoTerm`,
1252: `TaoTermType`,
1253: `TaoTermSumGetNumberTerms()`,
1254: `TaoTermSumSetNumberTerms()`,
1255: `TaoTermSumGetTerm()`,
1256: `TaoTermSumSetTerm()`,
1257: `TaoTermSumAddTerm()`,
1258: `TaoTermSumGetTermHessianMatrices()`,
1259: `TaoTermSumSetTermHessianMatrices()`,
1260: `TaoTermSumGetTermMask()`,
1261: `TaoTermSumSetTermMask()`
1262: M*/
1263: PETSC_INTERN PetscErrorCode TaoTermCreate_Sum(TaoTerm term)
1264: {
1265: TaoTerm_Sum *sum;
1267: PetscFunctionBegin;
1268: PetscCall(PetscNew(&sum));
1269: term->data = (void *)sum;
1270: term->parameters_mode = TAOTERM_PARAMETERS_OPTIONAL;
1272: PetscCall(PetscFree(term->H_mattype));
1273: PetscCall(PetscFree(term->Hpre_mattype));
1275: PetscCall(PetscStrallocpy(MATAIJ, (char **)&term->H_mattype));
1276: PetscCall(PetscStrallocpy(MATAIJ, (char **)&term->Hpre_mattype));
1277: term->Hpre_is_H = PETSC_TRUE;
1279: term->ops->destroy = TaoTermDestroy_Sum;
1280: term->ops->view = TaoTermView_Sum;
1281: term->ops->setfromoptions = TaoTermSetFromOptions_Sum;
1282: term->ops->objective = TaoTermComputeObjective_Sum;
1283: term->ops->gradient = TaoTermComputeGradient_Sum;
1284: term->ops->objectiveandgradient = TaoTermComputeObjectiveAndGradient_Sum;
1285: term->ops->hessian = TaoTermComputeHessian_Sum;
1286: term->ops->setup = TaoTermSetUp_Sum;
1287: term->ops->createsolutionvec = TaoTermCreateSolutionVec_Sum;
1288: term->ops->createparametersvec = TaoTermCreateParametersVec_Sum;
1289: term->ops->createhessianmatrices = TaoTermCreateHessianMatrices_Sum;
1290: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumGetNumberTerms_C", TaoTermSumGetNumberTerms_Sum));
1291: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumSetNumberTerms_C", TaoTermSumSetNumberTerms_Sum));
1292: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumGetTerm_C", TaoTermSumGetTerm_Sum));
1293: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumSetTerm_C", TaoTermSumSetTerm_Sum));
1294: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumAddTerm_C", TaoTermSumAddTerm_Sum));
1295: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumGetTermHessianMatrices_C", TaoTermSumGetTermHessianMatrices_Sum));
1296: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumSetTermHessianMatrices_C", TaoTermSumSetTermHessianMatrices_Sum));
1297: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumGetTermMask_C", TaoTermSumGetTermMask_Sum));
1298: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumSetTermMask_C", TaoTermSumSetTermMask_Sum));
1299: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumGetLastTermObjectives_C", TaoTermSumGetLastTermObjectives_Sum));
1300: PetscFunctionReturn(PETSC_SUCCESS);
1301: }