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)) {
462: PetscCall(PetscViewerASCIIPrintf(viewer, "%shessian", preceding ? ", " : ""));
463: preceding = PETSC_TRUE;
464: }
465: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
466: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
467: }
469: PetscCall(PetscViewerASCIIPopTab(viewer));
470: }
471: PetscCall(PetscViewerASCIIPopTab(viewer));
472: }
473: PetscFunctionReturn(PETSC_SUCCESS);
474: }
476: /*@
477: TaoTermSumSetNumberTerms - Set the number of terms in the sum
479: Collective
481: Input Parameters:
482: + term - a `TaoTerm` of type `TAOTERMSUM`
483: - n_terms - the number of terms that will be in the sum
485: Level: developer
487: Note:
488: If `n_terms` is smaller than the current number of terms, the trailing terms will be dropped.
490: .seealso: [](sec_tao_term),
491: `TaoTerm`,
492: `TAOTERMSUM`,
493: `TaoTermSumGetNumberTerms()`
494: @*/
495: PetscErrorCode TaoTermSumSetNumberTerms(TaoTerm term, PetscInt n_terms)
496: {
497: PetscFunctionBegin;
500: PetscTryMethod(term, "TaoTermSumSetNumberTerms_C", (TaoTerm, PetscInt), (term, n_terms));
501: PetscFunctionReturn(PETSC_SUCCESS);
502: }
504: static PetscErrorCode TaoTermSumSetNumberTerms_Sum(TaoTerm term, PetscInt n_terms)
505: {
506: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
507: PetscInt n_terms_old = sum->n_terms;
508: PetscReal *new_values;
509: TaoTermMapping *new_summands;
511: PetscFunctionBegin;
512: if (n_terms == n_terms_old) PetscFunctionReturn(PETSC_SUCCESS);
513: for (PetscInt i = n_terms; i < n_terms_old; i++) PetscCall(TaoTermMappingReset(&sum->terms[i]));
514: PetscCall(PetscMalloc1(n_terms, &new_summands));
515: PetscCall(PetscCalloc1(n_terms, &new_values));
516: PetscCall(PetscArraycpy(new_summands, sum->terms, PetscMin(n_terms, n_terms_old)));
517: PetscCall(PetscArrayzero(&new_summands[n_terms_old], PetscMax(0, n_terms - n_terms_old)));
518: PetscCall(PetscFree(sum->terms));
519: PetscCall(PetscFree(sum->subterm_values));
520: sum->terms = new_summands;
521: sum->subterm_values = new_values;
522: sum->n_terms = n_terms;
523: for (PetscInt i = n_terms_old; i < n_terms; i++) PetscCall(TaoTermSumSetTerm(term, i, NULL, 1.0, NULL, NULL));
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }
527: /*@
528: TaoTermSumGetNumberTerms - Get the number of terms in the sum
530: Not collective
532: Input Parameter:
533: . term - a `TaoTerm` of type `TAOTERMSUM`
535: Output Parameter:
536: . n_terms - the number of terms that will be in the sum
538: Level: developer
540: .seealso: [](sec_tao_term),
541: `TaoTerm`,
542: `TAOTERMSUM`,
543: `TaoTermSumSetNumberTerms()`
544: @*/
545: PetscErrorCode TaoTermSumGetNumberTerms(TaoTerm term, PetscInt *n_terms)
546: {
547: PetscFunctionBegin;
549: PetscAssertPointer(n_terms, 2);
550: PetscUseMethod(term, "TaoTermSumGetNumberTerms_C", (TaoTerm, PetscInt *), (term, n_terms));
551: PetscFunctionReturn(PETSC_SUCCESS);
552: }
554: static PetscErrorCode TaoTermSumGetNumberTerms_Sum(TaoTerm term, PetscInt *n_terms)
555: {
556: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
558: PetscFunctionBegin;
559: *n_terms = sum->n_terms;
560: PetscFunctionReturn(PETSC_SUCCESS);
561: }
563: /*@
564: TaoTermSumGetTerm - Get the data for a term in a `TAOTERMSUM`
566: Not collective
568: Input Parameters:
569: + sumterm - a `TaoTerm` of type `TAOTERMSUM`
570: - index - a number $0 \leq i < n$, where $n$ is the number of terms in `TaoTermSumGetNumberTerms()`
572: Output Parameters:
573: + prefix - (optional) the prefix used for configuring the term
574: . scale - (optional) the coefficient scaling the term in the sum
575: . term - the `TaoTerm` at given index of `TAOTERMSUM`
576: - map - (optional) a map from the `TAOTERMSUM` solution space to the `term` solution space; if `NULL` the map is assumed to be the identity
578: Level: developer
580: .seealso: [](sec_tao_term),
581: `TaoTerm`,
582: `TAOTERMSUM`,
583: `TaoTermSumSetTerm()`,
584: `TaoTermSumAddTerm()`
585: @*/
586: PetscErrorCode TaoTermSumGetTerm(TaoTerm sumterm, PetscInt index, const char **prefix, PetscReal *scale, TaoTerm *term, Mat *map)
587: {
588: PetscFunctionBegin;
591: if (prefix) PetscAssertPointer(prefix, 3);
592: if (term) PetscAssertPointer(term, 5);
593: if (scale) PetscAssertPointer(scale, 4);
594: if (map) PetscAssertPointer(map, 6);
595: PetscUseMethod(sumterm, "TaoTermSumGetTerm_C", (TaoTerm, PetscInt, const char **, PetscReal *, TaoTerm *, Mat *), (sumterm, index, prefix, scale, term, map));
596: PetscFunctionReturn(PETSC_SUCCESS);
597: }
599: static PetscErrorCode TaoTermSumGetTerm_Sum(TaoTerm term, PetscInt index, const char **prefix, PetscReal *scale, TaoTerm *subterm, Mat *map)
600: {
601: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
602: TaoTermMapping *summand;
604: PetscFunctionBegin;
605: 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);
606: summand = &sum->terms[index];
607: PetscCall(TaoTermMappingGetData(summand, prefix, scale, subterm, map));
608: PetscFunctionReturn(PETSC_SUCCESS);
609: }
611: /*@
612: TaoTermSumSetTerm - Set a term in a sum of terms
614: Collective
616: Input Parameters:
617: + sumterm - a `TaoTerm` of type `TAOTERMSUM`
618: . index - a number $0 \leq i < n$, where $n$ is the number of terms in `TaoTermSumSetNumberTerms()`
619: . prefix - (optional) the prefix used for configuring the term (if `NULL`, `term_x_` will be the prefix, e.g. "term_0_", "term_1_", etc.)
620: . scale - the coefficient scaling the term in the sum
621: . term - the `TaoTerm` to be set in `TAOTERMSUM`
622: - map - (optional) a map from the `TAOTERMSUM` solution space to the `term` solution space; if `NULL` the map is assumed to be the identity
624: Level: developer
626: .seealso: [](sec_tao_term),
627: `TaoTerm`,
628: `TAOTERMSUM`,
629: `TaoTermSumGetTerm()`,
630: `TaoTermSumAddTerm()`
631: @*/
632: PetscErrorCode TaoTermSumSetTerm(TaoTerm sumterm, PetscInt index, const char prefix[], PetscReal scale, TaoTerm term, Mat map)
633: {
634: PetscFunctionBegin;
637: if (prefix) PetscAssertPointer(prefix, 3);
641: PetscTryMethod(sumterm, "TaoTermSumSetTerm_C", (TaoTerm, PetscInt, const char[], PetscReal, TaoTerm, Mat), (sumterm, index, prefix, scale, term, map));
642: PetscFunctionReturn(PETSC_SUCCESS);
643: }
645: static PetscErrorCode TaoTermSumSetTerm_Sum(TaoTerm term, PetscInt index, const char prefix[], PetscReal scale, TaoTerm subterm, Mat map)
646: {
647: char subterm_x_[256];
648: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
649: TaoTermMapping *summand;
651: PetscFunctionBegin;
652: 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);
653: summand = &sum->terms[index];
654: if (prefix == NULL) {
655: PetscCall(PetscSNPrintf(subterm_x_, 256, "term_%" PetscInt_FMT "_", index));
656: prefix = subterm_x_;
657: }
658: PetscCall(TaoTermMappingSetData(summand, prefix, scale, subterm, map));
659: PetscFunctionReturn(PETSC_SUCCESS);
660: }
662: /*@
663: TaoTermSumSetTermHessianMatrices - Set Hessian matrices that can be used internally by a `TAOTERMSUM`
665: Logically collective
667: Input Parameters:
668: + term - a `TaoTerm` of type `TAOTERMSUM`
669: . index - the index for the term from `TaoTermSumSetTerm()` or `TaoTermSumAddTerm()`
670: . unmapped_H - (optional) unmapped Hessian matrix
671: . unmapped_Hpre - (optional) unmapped matrix for constructing the preconditioner of `unmapped_H`
672: . mapped_H - (optional) Hessian matrix
673: - mapped_Hpre - (optional) matrix for constructing the preconditioner of `mapped_H`
675: Level: developer
677: Notes:
678: If the inner term has the form $g(x) = \alpha f(Ax; p)$, the "mapped" Hessians should be able to hold the Hessian
679: $\nabla^2 g$ and the unmapped Hessians should be able to hold the Hessian $\nabla_x^2 f$. If the term is not mapped,
680: just pass the unmapped Hessians (e.g. `TaoTermSumSetTermHessianMatrices(term, 0, H, Hpre, NULL, NULL)`).
682: .seealso: [](sec_tao_term),
683: `TaoTerm`,
684: `TAOTERMSUM`,
685: `TaoTermComputeHessian()`,
686: `TaoTermSumGetTermHessianMatrices()`
687: @*/
688: PetscErrorCode TaoTermSumSetTermHessianMatrices(TaoTerm term, PetscInt index, Mat unmapped_H, Mat unmapped_Hpre, Mat mapped_H, Mat mapped_Hpre)
689: {
690: PetscFunctionBegin;
697: PetscTryMethod(term, "TaoTermSumSetTermHessianMatrices_C", (TaoTerm, PetscInt, Mat, Mat, Mat, Mat), (term, index, unmapped_H, unmapped_Hpre, mapped_H, mapped_Hpre));
698: PetscFunctionReturn(PETSC_SUCCESS);
699: }
701: static PetscErrorCode TaoTermSumSetTermHessianMatrices_Sum(TaoTerm term, PetscInt index, Mat unmapped_H, Mat unmapped_Hpre, Mat mapped_H, Mat mapped_Hpre)
702: {
703: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
704: TaoTermMapping *summand;
705: PetscBool is_callback;
707: PetscFunctionBegin;
708: 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);
709: summand = &sum->terms[index];
711: PetscCall(PetscObjectTypeCompare((PetscObject)summand->term, TAOTERMCALLBACKS, &is_callback));
712: if (is_callback) {
713: MatType H_type, Hpre_type;
714: PetscBool Hpre_is_H;
716: Hpre_is_H = (mapped_H == mapped_Hpre) ? PETSC_TRUE : PETSC_FALSE;
718: if (mapped_H) PetscCall(MatGetType(mapped_H, &H_type));
719: else H_type = NULL;
720: if (mapped_Hpre) PetscCall(MatGetType(mapped_Hpre, &Hpre_type));
721: else Hpre_type = NULL;
722: PetscCall(TaoTermSetCreateHessianMode(summand->term, Hpre_is_H, H_type, Hpre_type));
723: }
725: PetscCall(PetscObjectReference((PetscObject)unmapped_H));
726: PetscCall(MatDestroy(&summand->_unmapped_H));
727: summand->_unmapped_H = unmapped_H;
729: PetscCall(PetscObjectReference((PetscObject)unmapped_Hpre));
730: PetscCall(MatDestroy(&summand->_unmapped_Hpre));
731: summand->_unmapped_Hpre = unmapped_Hpre;
733: PetscCall(PetscObjectReference((PetscObject)mapped_H));
734: PetscCall(MatDestroy(&summand->_mapped_H));
735: summand->_mapped_H = mapped_H;
737: PetscCall(PetscObjectReference((PetscObject)mapped_Hpre));
738: PetscCall(MatDestroy(&summand->_mapped_Hpre));
739: summand->_mapped_Hpre = mapped_Hpre;
740: PetscFunctionReturn(PETSC_SUCCESS);
741: }
743: /*@
744: TaoTermSumGetTermHessianMatrices - Get Hessian matrices set with `TaoTermSumSetTermHessianMatrices()`.
746: Not collective
748: Input Parameters:
749: + term - a `TaoTerm` of type `TAOTERMSUM`
750: - index - the index for the term from `TaoTermSumSetTerm()` or `TaoTermSumAddTerm()`
752: Output Parameters:
753: + unmapped_H - (optional) unmapped Hessian matrix
754: . unmapped_Hpre - (optional) unmapped matrix for constructing the preconditioner for `unmapped_H`
755: . mapped_H - (optional) Hessian matrix
756: - mapped_Hpre - (optional) matrix for constructing the preconditioner for `mapped_H`
758: Level: developer
760: .seealso: [](sec_tao_term),
761: `TaoTerm`,
762: `TAOTERMSUM`,
763: `TaoTermComputeHessian()`,
764: `TaoTermSumSetTermHessianMatrices()`
765: @*/
766: PetscErrorCode TaoTermSumGetTermHessianMatrices(TaoTerm term, PetscInt index, Mat *unmapped_H, Mat *unmapped_Hpre, Mat *mapped_H, Mat *mapped_Hpre)
767: {
768: PetscFunctionBegin;
770: if (unmapped_H) PetscAssertPointer(unmapped_H, 3);
771: if (unmapped_Hpre) PetscAssertPointer(unmapped_Hpre, 4);
772: if (mapped_H) PetscAssertPointer(mapped_H, 5);
773: if (mapped_Hpre) PetscAssertPointer(mapped_Hpre, 6);
774: PetscTryMethod(term, "TaoTermSumGetTermHessianMatrices_C", (TaoTerm, PetscInt, Mat *, Mat *, Mat *, Mat *), (term, index, unmapped_H, unmapped_Hpre, mapped_H, mapped_Hpre));
775: PetscFunctionReturn(PETSC_SUCCESS);
776: }
778: static PetscErrorCode TaoTermSumGetTermHessianMatrices_Sum(TaoTerm term, PetscInt index, Mat *unmapped_H, Mat *unmapped_Hpre, Mat *mapped_H, Mat *mapped_Hpre)
779: {
780: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
781: TaoTermMapping *summand;
783: PetscFunctionBegin;
784: 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);
785: summand = &sum->terms[index];
787: if (unmapped_H) *unmapped_H = summand->_unmapped_H;
788: if (unmapped_Hpre) *unmapped_Hpre = summand->_unmapped_Hpre;
789: if (mapped_H) *mapped_H = summand->_mapped_H;
790: if (mapped_Hpre) *mapped_Hpre = summand->_mapped_Hpre;
791: PetscFunctionReturn(PETSC_SUCCESS);
792: }
794: /*@
795: TaoTermSumGetTermMask - Get the `TaoTermMask` of a term in the sum
797: Not collective
799: Input Parameters:
800: + term - a `TaoTerm` of type `TAOTERMSUM`
801: - index - the index for the term from `TaoTermSumSetTerm()` or `TaoTermSumAddTerm()`
803: Output Parameter:
804: . 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`)
806: Level: developer
808: .seealso: [](sec_tao_term),
809: `TaoTerm`,
810: `TAOTERMSUM`,
811: `TaoTermSumSetTermMask()`
812: @*/
813: PetscErrorCode TaoTermSumGetTermMask(TaoTerm term, PetscInt index, TaoTermMask *mask)
814: {
815: PetscFunctionBegin;
817: PetscAssertPointer(mask, 3);
818: PetscUseMethod(term, "TaoTermSumGetTermMask_C", (TaoTerm, PetscInt, TaoTermMask *), (term, index, mask));
819: PetscFunctionReturn(PETSC_SUCCESS);
820: }
822: static PetscErrorCode TaoTermSumGetTermMask_Sum(TaoTerm term, PetscInt index, TaoTermMask *mask)
823: {
824: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
825: TaoTermMapping *summand;
827: PetscFunctionBegin;
828: 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);
829: summand = &sum->terms[index];
830: *mask = summand->mask;
831: PetscFunctionReturn(PETSC_SUCCESS);
832: }
834: /*@
835: TaoTermSumSetTermMask - Set a `TaoTermMask` on a term in the sum
837: Logically collective
839: Input Parameters:
840: + term - a `TaoTerm` of type `TAOTERMSUM`
841: . index - the index for the term from `TaoTermSumSetTerm()` or `TaoTermSumAddTerm()`
842: - 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`)
844: Options Database Keys:
845: . -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()`)
847: Level: developer
849: Note:
850: Some optimization methods may add a damping term to the Hessian of an
851: objective function without affecting the objective or gradient. If, e.g.,
852: the regularizer has index `1`, then this can be accomplished with
853: `TaoTermSumSetTermMask(term, 1, TAOTERM_MASK_OBJECTIVE | TAOTERM_MASK_GRADIENT)`.
855: .seealso: [](sec_tao_term),
856: `TaoTerm`,
857: `TAOTERMSUM`,
858: `TaoTermSumGetTermMask()`
859: @*/
860: PetscErrorCode TaoTermSumSetTermMask(TaoTerm term, PetscInt index, TaoTermMask mask)
861: {
862: PetscFunctionBegin;
865: PetscTryMethod(term, "TaoTermSumSetTermMask_C", (TaoTerm, PetscInt, TaoTermMask), (term, index, mask));
866: PetscFunctionReturn(PETSC_SUCCESS);
867: }
869: static PetscErrorCode TaoTermSumSetTermMask_Sum(TaoTerm term, PetscInt index, TaoTermMask mask)
870: {
871: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
872: TaoTermMapping *summand;
874: PetscFunctionBegin;
875: 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);
876: summand = &sum->terms[index];
877: summand->mask = mask;
878: PetscFunctionReturn(PETSC_SUCCESS);
879: }
881: /*@
882: TaoTermSumAddTerm - Append a term to the terms being summed
884: Collective
886: Input Parameters:
887: + sumterm - a `TaoTerm` of type `TAOTERMSUM`
888: . 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.)
889: . scale - the coefficient scaling the term in the sum
890: . term - the `TaoTerm` to add
891: - map - (optional) a map from the `TAOTERMSUM` solution space to the `term` solution space; if `NULL` the map is assumed to be the identity
893: Output Parameter:
894: . index - (optional) the index of the newly added term
896: Level: developer
898: .seealso: [](sec_tao_term), `TaoTerm`, `TAOTERMSUM`
899: @*/
900: PetscErrorCode TaoTermSumAddTerm(TaoTerm sumterm, const char prefix[], PetscReal scale, TaoTerm term, Mat map, PetscInt *index)
901: {
902: PetscFunctionBegin;
904: if (prefix) PetscAssertPointer(prefix, 2);
908: if (index) PetscAssertPointer(index, 6);
909: PetscTryMethod(sumterm, "TaoTermSumAddTerm_C", (TaoTerm, const char[], PetscReal, TaoTerm, Mat, PetscInt *), (sumterm, prefix, scale, term, map, index));
910: PetscFunctionReturn(PETSC_SUCCESS);
911: }
913: static PetscErrorCode TaoTermSumAddTerm_Sum(TaoTerm term, const char prefix[], PetscReal scale, TaoTerm subterm, Mat map, PetscInt *index)
914: {
915: PetscInt n_terms_old;
917: PetscFunctionBegin;
918: PetscCall(TaoTermSumGetNumberTerms(term, &n_terms_old));
919: PetscCall(TaoTermSumSetNumberTerms(term, n_terms_old + 1));
920: PetscCall(TaoTermSumSetTerm(term, n_terms_old, prefix, scale, subterm, map));
921: if (index) *index = n_terms_old;
922: PetscFunctionReturn(PETSC_SUCCESS);
923: }
925: static PetscErrorCode TaoTermSetFromOptions_Sum(TaoTerm term, PetscOptionItems PetscOptionsObject)
926: {
927: PetscInt n_terms;
928: const char *prefix;
930: PetscFunctionBegin;
931: PetscCall(TaoTermSumGetNumberTerms(term, &n_terms));
932: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)term, &prefix));
933: PetscOptionsHeadBegin(PetscOptionsObject, "TaoTerm sum options");
934: PetscCall(PetscOptionsBoundedInt("-tao_term_sum_number_terms", "The number of terms in the sum", "TaoTermSumSetNumberTerms", n_terms, &n_terms, NULL, 0));
935: PetscCall(TaoTermSumSetNumberTerms(term, n_terms));
936: for (PetscInt i = 0; i < n_terms; i++) {
937: const char *subprefix;
938: Mat map;
939: PetscReal scale;
940: TaoTerm subterm;
941: char arg[256];
942: PetscBool flg;
943: PetscEnum masks[4] = {ENUM_DUMMY, ENUM_DUMMY, ENUM_DUMMY, ENUM_DUMMY};
944: PetscInt n_masks = PETSC_STATIC_ARRAY_LENGTH(masks);
946: PetscCall(TaoTermSumGetTerm(term, i, &subprefix, &scale, &subterm, &map));
947: if (subterm == NULL) {
948: PetscCall(TaoTermDuplicate(term, TAOTERM_DUPLICATE_SIZEONLY, &subterm));
949: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)subterm, prefix));
950: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)subterm, subprefix));
951: } else PetscCall(PetscObjectReference((PetscObject)subterm));
952: PetscCall(TaoTermSetFromOptions(subterm));
954: PetscCall(PetscSNPrintf(arg, 256, "-tao_term_sum_%sscale", subprefix));
955: PetscCall(PetscOptionsReal(arg, "The scale of the term in the TaoTermSum", "TaoTermSumSetTerm", scale, &scale, NULL));
957: PetscCall(PetscSNPrintf(arg, 256, "-tao_term_sum_%smask", subprefix));
958: PetscCall(PetscOptionsEnumArray(arg, "The mask of the term in the TaoTermSum", "TaoTermSumSetTermMask", TaoTermMasks, masks, &n_masks, &flg));
959: if (flg) {
960: PetscEnum mask = (PetscEnum)TAOTERM_MASK_NONE;
962: for (PetscInt j = 0; j < n_masks; j++) {
963: PetscEnum this_mask = masks[j] ? (PetscEnum)(1 << (masks[j] - 1)) : (PetscEnum)TAOTERM_MASK_NONE;
965: mask = (PetscEnum)(mask | this_mask);
966: }
968: PetscCall(TaoTermSumSetTermMask(term, i, (TaoTermMask)mask));
969: }
971: if (map) PetscCall(MatSetFromOptions(map));
972: PetscCall(TaoTermSumSetTerm(term, i, subprefix, scale, subterm, map));
973: PetscCall(TaoTermDestroy(&subterm));
974: }
975: PetscOptionsHeadEnd();
976: PetscFunctionReturn(PETSC_SUCCESS);
977: }
979: /*@
980: TaoTermSumGetLastTermObjectives - Get the contributions from each term to the
981: last evaluation of `TaoTermComputeObjective()` or `TaoTermComputeObjectiveAndGradient()`
983: Not collective
985: Input Parameter:
986: . term - a `TaoTerm` of type `TAOTERMSUM`
988: Output Parameter:
989: . values - an array of the contributions to the last computed objective value
991: Level: developer
993: .seealso: [](sec_tao_term),
994: `TaoTerm`,
995: `TAOTERMSUM`
996: @*/
997: PetscErrorCode TaoTermSumGetLastTermObjectives(TaoTerm term, const PetscReal *values[])
998: {
999: PetscFunctionBegin;
1001: PetscAssertPointer(values, 2);
1002: PetscUseMethod(term, "TaoTermSumGetLastTermObjectives_C", (TaoTerm, const PetscReal *[]), (term, values));
1003: PetscFunctionReturn(PETSC_SUCCESS);
1004: }
1006: static PetscErrorCode TaoTermSumGetLastTermObjectives_Sum(TaoTerm term, const PetscReal *values[])
1007: {
1008: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
1010: PetscFunctionBegin;
1011: *values = sum->subterm_values;
1012: PetscFunctionReturn(PETSC_SUCCESS);
1013: }
1015: static PetscErrorCode TaoTermComputeObjective_Sum(TaoTerm term, Vec x, Vec params, PetscReal *value)
1016: {
1017: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
1018: Vec *sub_params = NULL;
1019: PetscBool *is_dummy = NULL;
1020: PetscReal value_;
1021: PetscReal *values = sum->subterm_values;
1023: PetscFunctionBegin;
1024: if (params) PetscCall(TaoTermSumVecNestGetSubVecsRead(params, NULL, &sub_params, &is_dummy));
1025: value_ = 0.0;
1026: for (PetscInt i = 0; i < sum->n_terms; i++) {
1027: TaoTermMapping *summand = &sum->terms[i];
1028: Vec sub_param = TaoTermSumGetSubVec(params, sub_params, is_dummy, i);
1030: PetscCall(TaoTermMappingComputeObjective(summand, x, sub_param, INSERT_VALUES, &values[i]));
1031: value_ += values[i];
1032: }
1033: if (params) PetscCall(TaoTermSumVecNestRestoreSubVecsRead(params, NULL, &sub_params, &is_dummy));
1034: *value = value_;
1035: PetscFunctionReturn(PETSC_SUCCESS);
1036: }
1038: static PetscErrorCode TaoTermComputeGradient_Sum(TaoTerm term, Vec x, Vec params, Vec g)
1039: {
1040: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
1041: Vec *sub_params = NULL;
1042: PetscBool *is_dummy = NULL;
1044: PetscFunctionBegin;
1045: if (params) PetscCall(TaoTermSumVecNestGetSubVecsRead(params, NULL, &sub_params, &is_dummy));
1046: for (PetscInt i = 0; i < sum->n_terms; i++) {
1047: TaoTermMapping *summand = &sum->terms[i];
1048: Vec sub_param = TaoTermSumGetSubVec(params, sub_params, is_dummy, i);
1050: PetscCall(TaoTermMappingComputeGradient(summand, x, sub_param, i == 0 ? INSERT_VALUES : ADD_VALUES, g));
1051: }
1052: if (params) PetscCall(TaoTermSumVecNestRestoreSubVecsRead(params, NULL, &sub_params, &is_dummy));
1053: PetscFunctionReturn(PETSC_SUCCESS);
1054: }
1056: static PetscErrorCode TaoTermComputeObjectiveAndGradient_Sum(TaoTerm term, Vec x, Vec params, PetscReal *value, Vec g)
1057: {
1058: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
1059: Vec *sub_params = NULL;
1060: PetscBool *is_dummy = NULL;
1061: PetscReal *values = sum->subterm_values;
1062: PetscReal value_;
1064: PetscFunctionBegin;
1065: if (params) PetscCall(TaoTermSumVecNestGetSubVecsRead(params, NULL, &sub_params, &is_dummy));
1066: value_ = 0.0;
1067: for (PetscInt i = 0; i < sum->n_terms; i++) {
1068: TaoTermMapping *summand = &sum->terms[i];
1069: Vec sub_param = TaoTermSumGetSubVec(params, sub_params, is_dummy, i);
1071: values[i] = 0.0;
1072: PetscCall(TaoTermMappingComputeObjectiveAndGradient(summand, x, sub_param, i == 0 ? INSERT_VALUES : ADD_VALUES, &values[i], g));
1073: value_ += values[i];
1074: }
1075: if (params) PetscCall(TaoTermSumVecNestRestoreSubVecsRead(params, NULL, &sub_params, &is_dummy));
1076: *value = value_;
1077: PetscFunctionReturn(PETSC_SUCCESS);
1078: }
1080: static PetscErrorCode TaoTermComputeHessian_Sum(TaoTerm term, Vec x, Vec params, Mat H, Mat Hpre)
1081: {
1082: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
1083: Vec *sub_params = NULL;
1084: PetscBool *is_dummy = NULL;
1086: PetscFunctionBegin;
1087: if (H == NULL && Hpre == NULL) PetscFunctionReturn(PETSC_SUCCESS);
1088: if (params) PetscCall(TaoTermSumVecNestGetSubVecsRead(params, NULL, &sub_params, &is_dummy));
1089: // If mattype dense, then after zero entries, H->assembled = true.
1090: // But for aij, H->assembled is still false.
1091: if (H) {
1092: PetscCall(MatZeroEntries(H));
1093: PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
1094: PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
1095: }
1096: if (Hpre && (Hpre != H)) {
1097: PetscCall(MatZeroEntries(Hpre));
1098: PetscCall(MatAssemblyBegin(Hpre, MAT_FINAL_ASSEMBLY));
1099: PetscCall(MatAssemblyEnd(Hpre, MAT_FINAL_ASSEMBLY));
1100: }
1101: for (PetscInt i = 0; i < sum->n_terms; i++) {
1102: TaoTermMapping *summand = &sum->terms[i];
1103: Vec sub_param = TaoTermSumGetSubVec(params, sub_params, is_dummy, i);
1105: PetscCall(TaoTermMappingComputeHessian(summand, x, sub_param, ADD_VALUES, H, Hpre == H ? NULL : Hpre));
1106: }
1107: if (params) PetscCall(TaoTermSumVecNestRestoreSubVecsRead(params, NULL, &sub_params, &is_dummy));
1108: PetscFunctionReturn(PETSC_SUCCESS);
1109: }
1111: static PetscErrorCode TaoTermSetUp_Sum(TaoTerm term)
1112: {
1113: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
1114: PetscBool all_none = PETSC_TRUE;
1115: PetscBool any_required = PETSC_FALSE;
1116: PetscInt k = 0, K = 0;
1117: Mat *mats, new_parameters_factory;
1118: PetscLayout layout = NULL, clayout;
1120: PetscFunctionBegin;
1121: PetscCall(PetscCalloc1(sum->n_terms, &mats));
1122: PetscCall(MatGetLayouts(term->solution_factory, &layout, &clayout));
1123: if (layout->setupcalled == PETSC_FALSE) layout = NULL;
1124: for (PetscInt i = 0; i < sum->n_terms; i++) {
1125: TaoTermMapping *summand = &sum->terms[i];
1126: TaoTermParametersMode submode;
1127: PetscLayout sub_layout;
1128: PetscBool congruent;
1130: PetscCall(TaoTermSetUp(summand->term));
1131: if (summand->map) {
1132: PetscCall(MatSetUp(summand->map));
1133: PetscCall(MatGetLayouts(summand->map, NULL, &sub_layout));
1134: } else PetscCall(TaoTermGetSolutionLayout(summand->term, &sub_layout));
1135: if (i == 0 && layout == NULL) layout = sub_layout;
1136: PetscCall(PetscLayoutCompare(layout, sub_layout, &congruent));
1137: if (congruent == PETSC_FALSE) {
1138: PetscInt N, sub_N;
1139: MPI_Comm comm = PetscObjectComm((PetscObject)term);
1141: PetscCall(PetscLayoutGetSize(layout, &N));
1142: PetscCall(PetscLayoutGetSize(sub_layout, &sub_N));
1144: 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);
1145: }
1146: PetscCall(TaoTermGetParametersMode(summand->term, &submode));
1147: if (submode == TAOTERM_PARAMETERS_REQUIRED) any_required = PETSC_TRUE;
1148: if (submode != TAOTERM_PARAMETERS_NONE) {
1149: PetscInt subk, subK;
1151: all_none = PETSC_FALSE;
1152: PetscCall(TaoTermGetParametersSizes(summand->term, &subk, &subK, NULL));
1153: k += subk;
1154: K += subK;
1155: }
1156: if (summand->term->parameters_mode != TAOTERM_PARAMETERS_NONE) {
1157: PetscCall(PetscObjectReference((PetscObject)summand->term->parameters_factory));
1158: mats[i] = summand->term->parameters_factory;
1159: } else {
1160: PetscCall(MatCreate(PetscObjectComm((PetscObject)term), &mats[i]));
1161: PetscCall(MatSetType(mats[i], MATDUMMY));
1162: PetscCall(MatSetSizes(mats[i], 0, 0, 0, 0));
1163: PetscCall(MatSetUp(mats[i]));
1164: }
1165: }
1166: PetscCall(MatSetLayouts(term->solution_factory, layout, clayout));
1167: if (all_none) {
1168: term->parameters_mode = TAOTERM_PARAMETERS_NONE;
1169: } else if (any_required) {
1170: term->parameters_mode = TAOTERM_PARAMETERS_REQUIRED;
1171: } else {
1172: term->parameters_mode = TAOTERM_PARAMETERS_OPTIONAL;
1173: }
1174: PetscCall(TaoTermSetParametersSizes(term, k, K, 1));
1175: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)term), sum->n_terms, NULL, 1, NULL, mats, &new_parameters_factory));
1176: PetscCall(MatSetVecType(new_parameters_factory, VECNEST));
1177: PetscCall(MatDestroy(&term->parameters_factory));
1178: term->parameters_factory = new_parameters_factory;
1179: for (PetscInt i = 0; i < sum->n_terms; i++) PetscCall(MatDestroy(&mats[i]));
1180: PetscCall(PetscFree(mats));
1181: PetscFunctionReturn(PETSC_SUCCESS);
1182: }
1184: static PetscErrorCode TaoTermCreateSolutionVec_Sum(TaoTerm term, Vec *solution_vec)
1185: {
1186: PetscFunctionBegin;
1187: if (solution_vec) PetscCall(MatCreateVecs(term->solution_factory, NULL, solution_vec));
1188: PetscFunctionReturn(PETSC_SUCCESS);
1189: }
1191: static PetscErrorCode TaoTermCreateParametersVec_Sum(TaoTerm term, Vec *parameters_vec)
1192: {
1193: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
1195: PetscFunctionBegin;
1196: if (parameters_vec) {
1197: Vec *vecs;
1199: PetscCall(PetscCalloc1(sum->n_terms, &vecs));
1200: for (PetscInt i = 0; i < sum->n_terms; i++) {
1201: TaoTermMapping *summand = &sum->terms[i];
1202: TaoTermParametersMode submode;
1204: PetscCall(TaoTermGetParametersMode(summand->term, &submode));
1205: if (submode != TAOTERM_PARAMETERS_NONE) PetscCall(TaoTermCreateParametersVec(summand->term, &vecs[i]));
1206: }
1207: PetscCall(TaoTermSumParametersPack(term, vecs, parameters_vec));
1208: for (PetscInt i = 0; i < sum->n_terms; i++) PetscCall(VecDestroy(&vecs[i]));
1209: PetscCall(PetscFree(vecs));
1210: }
1211: PetscFunctionReturn(PETSC_SUCCESS);
1212: }
1214: static PetscErrorCode TaoTermCreateHessianMatrices_Sum(TaoTerm term, Mat *H, Mat *Hpre)
1215: {
1216: TaoTerm_Sum *sum = (TaoTerm_Sum *)term->data;
1217: PetscBool Hpre_is_H, sub_Hpre_is_H;
1219: PetscFunctionBegin;
1220: Hpre_is_H = term->Hpre_is_H;
1221: // Need to create subterms' mapped Hessians and PtAP routines, if needed
1222: for (PetscInt i = 0; i < sum->n_terms; i++) {
1223: TaoTermMapping *summand = &sum->terms[i];
1224: PetscBool is_callback;
1226: PetscCall(PetscObjectTypeCompare((PetscObject)summand->term, TAOTERMCALLBACKS, &is_callback));
1227: if (is_callback) {
1228: Mat c_H;
1230: PetscCall(TaoTermSumGetTermHessianMatrices(term, i, NULL, NULL, &c_H, NULL));
1231: PetscCheck(c_H, PetscObjectComm((PetscObject)summand->term), PETSC_ERR_USER, "TAOTERMCALLBACKS does not have Hessian routines set. Call TaoSetHessian()");
1232: }
1233: PetscCall(TaoTermMappingCreateHessianMatrices(summand, &summand->_mapped_H, &summand->_mapped_Hpre));
1235: sub_Hpre_is_H = (summand->_mapped_H == summand->_mapped_Hpre) ? PETSC_TRUE : PETSC_FALSE;
1236: Hpre_is_H = (Hpre_is_H && sub_Hpre_is_H) ? PETSC_TRUE : PETSC_FALSE;
1237: }
1239: term->Hpre_is_H = Hpre_is_H;
1240: PetscCall(TaoTermCreateHessianMatricesDefault(term, H, Hpre));
1241: PetscFunctionReturn(PETSC_SUCCESS);
1242: }
1244: /*MC
1245: TAOTERMSUM - A `TaoTerm` that is a sum of multiple `TaoTerms`.
1247: Level: developer
1249: Note:
1250: The default Hessian creation mode (see `TaoTermGetCreateHessianMode()`) is `H == Hpre` and `TaoTermCreateHessianMatrices()`
1251: will create a `MATAIJ`.
1253: .seealso: [](sec_tao_term),
1254: `TaoTerm`,
1255: `TaoTermType`,
1256: `TaoTermSumGetNumberTerms()`,
1257: `TaoTermSumSetNumberTerms()`,
1258: `TaoTermSumGetTerm()`,
1259: `TaoTermSumSetTerm()`,
1260: `TaoTermSumAddTerm()`,
1261: `TaoTermSumGetTermHessianMatrices()`,
1262: `TaoTermSumSetTermHessianMatrices()`,
1263: `TaoTermSumGetTermMask()`,
1264: `TaoTermSumSetTermMask()`
1265: M*/
1266: PETSC_INTERN PetscErrorCode TaoTermCreate_Sum(TaoTerm term)
1267: {
1268: TaoTerm_Sum *sum;
1270: PetscFunctionBegin;
1271: PetscCall(PetscNew(&sum));
1272: term->data = (void *)sum;
1273: term->parameters_mode = TAOTERM_PARAMETERS_OPTIONAL;
1275: PetscCall(PetscFree(term->H_mattype));
1276: PetscCall(PetscFree(term->Hpre_mattype));
1278: PetscCall(PetscStrallocpy(MATAIJ, (char **)&term->H_mattype));
1279: PetscCall(PetscStrallocpy(MATAIJ, (char **)&term->Hpre_mattype));
1280: term->Hpre_is_H = PETSC_TRUE;
1282: term->ops->destroy = TaoTermDestroy_Sum;
1283: term->ops->view = TaoTermView_Sum;
1284: term->ops->setfromoptions = TaoTermSetFromOptions_Sum;
1285: term->ops->objective = TaoTermComputeObjective_Sum;
1286: term->ops->gradient = TaoTermComputeGradient_Sum;
1287: term->ops->objectiveandgradient = TaoTermComputeObjectiveAndGradient_Sum;
1288: term->ops->hessian = TaoTermComputeHessian_Sum;
1289: term->ops->setup = TaoTermSetUp_Sum;
1290: term->ops->createsolutionvec = TaoTermCreateSolutionVec_Sum;
1291: term->ops->createparametersvec = TaoTermCreateParametersVec_Sum;
1292: term->ops->createhessianmatrices = TaoTermCreateHessianMatrices_Sum;
1293: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumGetNumberTerms_C", TaoTermSumGetNumberTerms_Sum));
1294: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumSetNumberTerms_C", TaoTermSumSetNumberTerms_Sum));
1295: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumGetTerm_C", TaoTermSumGetTerm_Sum));
1296: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumSetTerm_C", TaoTermSumSetTerm_Sum));
1297: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumAddTerm_C", TaoTermSumAddTerm_Sum));
1298: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumGetTermHessianMatrices_C", TaoTermSumGetTermHessianMatrices_Sum));
1299: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumSetTermHessianMatrices_C", TaoTermSumSetTermHessianMatrices_Sum));
1300: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumGetTermMask_C", TaoTermSumGetTermMask_Sum));
1301: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumSetTermMask_C", TaoTermSumSetTermMask_Sum));
1302: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermSumGetLastTermObjectives_C", TaoTermSumGetLastTermObjectives_Sum));
1303: PetscFunctionReturn(PETSC_SUCCESS);
1304: }