Actual source code: taotermhalfl2squared.c
1: #include <petsc/private/taoimpl.h>
3: typedef struct _n_TaoTerm_HalfL2Squared TaoTerm_HalfL2Squared;
5: struct _n_TaoTerm_HalfL2Squared {
6: Vec pdiff_work;
7: };
9: static PetscErrorCode TaoTermDestroy_Halfl2squared(TaoTerm term)
10: {
11: TaoTerm_HalfL2Squared *l2 = (TaoTerm_HalfL2Squared *)term->data;
13: PetscFunctionBegin;
14: PetscCall(VecDestroy(&l2->pdiff_work));
15: PetscCall(PetscFree(l2));
16: term->data = NULL;
17: PetscCall(TaoTermDestroy_ElementwiseDivergence_Internal(term));
18: PetscFunctionReturn(PETSC_SUCCESS);
19: }
21: static PetscErrorCode TaoTermComputeObjective_Halfl2squared(TaoTerm term, Vec x, Vec params, PetscReal *value)
22: {
23: PetscScalar v;
24: TaoTerm_HalfL2Squared *l2 = (TaoTerm_HalfL2Squared *)term->data;
26: PetscFunctionBegin;
27: if (params) {
28: if (l2->pdiff_work == NULL) PetscCall(VecDuplicate(x, &l2->pdiff_work));
30: PetscCall(VecCopy(x, l2->pdiff_work));
31: PetscCall(VecAXPY(l2->pdiff_work, -1.0, params));
32: PetscCall(VecDot(l2->pdiff_work, l2->pdiff_work, &v));
33: } else PetscCall(VecDot(x, x, &v));
34: *value = 0.5 * PetscRealPart(v);
35: PetscFunctionReturn(PETSC_SUCCESS);
36: }
38: static PetscErrorCode TaoTermComputeObjectiveAndGradient_Halfl2squared(TaoTerm term, Vec x, Vec params, PetscReal *value, Vec g)
39: {
40: PetscScalar v;
42: PetscFunctionBegin;
43: if (params) PetscCall(VecWAXPY(g, -1.0, params, x));
44: else PetscCall(VecCopy(x, g));
45: PetscCall(VecDot(g, g, &v));
46: *value = 0.5 * PetscRealPart(v);
47: PetscFunctionReturn(PETSC_SUCCESS);
48: }
50: static PetscErrorCode TaoTermComputeGradient_Halfl2squared(TaoTerm term, Vec x, Vec params, Vec g)
51: {
52: PetscFunctionBegin;
53: if (params) PetscCall(VecWAXPY(g, -1.0, params, x));
54: else PetscCall(VecCopy(x, g));
55: PetscFunctionReturn(PETSC_SUCCESS);
56: }
58: static PetscErrorCode TaoTermComputeHessian_Halfl2squared(TaoTerm term, Vec x, Vec params, Mat H, Mat Hpre)
59: {
60: PetscFunctionBegin;
61: if (H) {
62: PetscCall(MatZeroEntries(H));
63: PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
64: PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
65: PetscCall(MatShift(H, 1.0));
66: }
67: if (Hpre && Hpre != H) {
68: PetscCall(MatZeroEntries(Hpre));
69: PetscCall(MatAssemblyBegin(Hpre, MAT_FINAL_ASSEMBLY));
70: PetscCall(MatAssemblyEnd(Hpre, MAT_FINAL_ASSEMBLY));
71: PetscCall(MatShift(Hpre, 1.0));
72: }
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: static PetscErrorCode TaoTermCreateHessianMatrices_Halfl2squared(TaoTerm term, Mat *H, Mat *Hpre)
77: {
78: PetscBool is_hdiag, is_hprediag;
80: PetscFunctionBegin;
81: PetscCall(PetscInfo(term, "%s: Creating TAOTERMHALFL2SQUARED Hessian Matrices. TAOTERMHALFL2SQUARED only accepts MATDIAGONAL for MatType, overriding any user-set MatType.\n", ((PetscObject)term)->prefix));
82: PetscCall(PetscStrcmp(term->H_mattype, MATDIAGONAL, &is_hdiag));
83: PetscCall(PetscStrcmp(term->Hpre_mattype, MATDIAGONAL, &is_hprediag));
84: if (!is_hdiag) {
85: PetscCall(PetscFree(term->H_mattype));
86: PetscCall(PetscStrallocpy(MATDIAGONAL, (char **)&term->H_mattype));
87: }
88: if (!is_hprediag) {
89: PetscCall(PetscFree(term->Hpre_mattype));
90: PetscCall(PetscStrallocpy(MATDIAGONAL, (char **)&term->Hpre_mattype));
91: }
92: PetscCall(TaoTermCreateHessianMatricesDefault(term, H, Hpre));
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
96: static PetscErrorCode TaoTermIsComputeHessianFDPossible_Halfl2squared(TaoTerm term, PetscBool3 *ispossible)
97: {
98: PetscFunctionBegin;
99: *ispossible = PETSC_BOOL3_FALSE;
100: PetscFunctionReturn(PETSC_SUCCESS);
101: }
103: /*MC
104: TAOTERMHALFL2SQUARED - A `TaoTerm` that computes $\tfrac{1}{2}\|x - p\|_2^2$, for solution $x$ and parameters $p$.
106: Level: intermediate
108: Notes:
109: By default this term is `TAOTERM_PARAMETERS_OPTIONAL`. If the parameters
110: argument is `NULL` in the evaluation routines (`TaoTermComputeObjective()`,
111: `TaoTermComputeGradient()`, etc.), then it is assumed $p = 0$ and the term computes
112: $\tfrac{1}{2}\|x\|_2^2$.
114: The default Hessian creation mode (see `TaoTermGetCreateHessianMode()`) is `H == Hpre` and `TaoTermCreateHessianMatrices()`
115: will create a `MATDIAGONAL` for the Hessian.
117: .seealso: [](sec_tao_term),
118: `TaoTerm`,
119: `TaoTermType`,
120: `TaoTermCreateHalfL2Squared()`,
121: `TAOTERML1`,
122: `TAOTERMQUADRATIC`
123: M*/
124: PETSC_INTERN PetscErrorCode TaoTermCreate_Halfl2squared(TaoTerm term)
125: {
126: TaoTerm_HalfL2Squared *l2;
128: PetscFunctionBegin;
129: PetscCall(TaoTermCreate_ElementwiseDivergence_Internal(term));
130: PetscCall(PetscNew(&l2));
131: term->data = (void *)l2;
133: PetscCall(PetscFree(term->H_mattype));
134: PetscCall(PetscFree(term->Hpre_mattype));
136: PetscCall(PetscStrallocpy(MATDIAGONAL, (char **)&term->H_mattype));
137: PetscCall(PetscStrallocpy(MATDIAGONAL, (char **)&term->Hpre_mattype));
139: term->ops->destroy = TaoTermDestroy_Halfl2squared;
140: term->ops->objective = TaoTermComputeObjective_Halfl2squared;
141: term->ops->gradient = TaoTermComputeGradient_Halfl2squared;
142: term->ops->objectiveandgradient = TaoTermComputeObjectiveAndGradient_Halfl2squared;
143: term->ops->hessian = TaoTermComputeHessian_Halfl2squared;
144: term->ops->createhessianmatrices = TaoTermCreateHessianMatrices_Halfl2squared;
145: term->ops->iscomputehessianfdpossible = TaoTermIsComputeHessianFDPossible_Halfl2squared;
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: /*@
150: TaoTermCreateHalfL2Squared - Create a `TaoTerm` for the objective term $\tfrac{1}{2}\|x - p\|_2^2$, for solution $x$ and parameters $p$.
152: Collective
154: Input Parameters:
155: + comm - the MPI communicator where the `TaoTerm` will be computed
156: . n - the local size of the $x$ and $p$ vectors (or `PETSC_DECIDE`)
157: - N - the global size of the $x$ and $p$ vectors (or `PETSC_DECIDE`)
159: Output Parameter:
160: . term - the `TaoTerm`
162: Level: beginner
164: Note:
165: If you would like to add a Tikhonov regularization term $\alpha \tfrac{1}{2}\|x\|_2^2$ to the objective function of a `Tao`, do the following\:
166: .vb
167: VecGetSizes(x, &n, &N);
168: TaoTermCreateHalfL2Squared(PetscObjectComm((PetscObject)x), n, N, &term);
169: TaoAddTerm(tao, "reg_", alpha, term, NULL, NULL);
170: TaoTermDestroy(&term);
171: .ve
172: If you would like to add a biased regularization term $\alpha \tfrac{1}{2}\|x - p \|_2^2$, do the same but pass `p` as the parameters of the term\:
173: .vb
174: TaoAddTerm(tao, "reg_", alpha, term, p, NULL);
175: .ve
177: .seealso: [](sec_tao_term),
178: `TaoTerm`,
179: `TAOTERMHALFL2SQUARED`,
180: `TaoTermCreateL1()`,
181: `TaoTermCreateQuadratic()`
182: @*/
183: PetscErrorCode TaoTermCreateHalfL2Squared(MPI_Comm comm, PetscInt n, PetscInt N, TaoTerm *term)
184: {
185: TaoTerm _term;
187: PetscFunctionBegin;
188: PetscAssertPointer(term, 4);
189: PetscCall(TaoTermCreate(comm, &_term));
190: PetscCall(TaoTermSetType(_term, TAOTERMHALFL2SQUARED));
191: PetscCall(TaoTermSetSolutionSizes(_term, n, N, 1));
192: *term = _term;
193: PetscFunctionReturn(PETSC_SUCCESS);
194: }