Actual source code: taotermquadratic.c
1: #include <petsc/private/taoimpl.h>
3: typedef struct _n_TaoTerm_Quadratic TaoTerm_Quadratic;
5: struct _n_TaoTerm_Quadratic {
6: Mat A;
7: Vec _diff;
8: Vec Adiff;
9: };
11: static PetscErrorCode TaoTermDestroy_Quadratic(TaoTerm term)
12: {
13: TaoTerm_Quadratic *quad = (TaoTerm_Quadratic *)term->data;
15: PetscFunctionBegin;
16: PetscCall(MatDestroy(&quad->A));
17: PetscCall(VecDestroy(&quad->_diff));
18: PetscCall(VecDestroy(&quad->Adiff));
19: PetscCall(PetscFree(quad));
20: term->data = NULL;
21: PetscCall(TaoTermDestroy_ElementwiseDivergence_Internal(term));
22: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermQuadraticSetMat_C", NULL));
23: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermQuadraticGetMat_C", NULL));
24: PetscFunctionReturn(PETSC_SUCCESS);
25: }
27: static PetscErrorCode TaoTermView_Quadratic(TaoTerm term, PetscViewer viewer)
28: {
29: TaoTerm_Quadratic *quad = (TaoTerm_Quadratic *)term->data;
31: PetscFunctionBegin;
32: if (quad->A) {
33: PetscBool iascii;
35: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
37: if (iascii) {
38: PetscViewerFormat format;
39: PetscBool pop = PETSC_FALSE;
41: PetscCall(PetscViewerGetFormat(viewer, &format));
42: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
43: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
44: pop = PETSC_TRUE;
45: }
46: PetscCall(MatView(quad->A, viewer));
47: if (pop) PetscCall(PetscViewerPopFormat(viewer));
48: }
49: }
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: static PetscErrorCode TaoTermQuadraticDiff(TaoTerm term, Vec x, Vec params, Vec *diff)
54: {
55: PetscFunctionBegin;
56: *diff = x;
57: if (params) {
58: TaoTerm_Quadratic *quad = (TaoTerm_Quadratic *)term->data;
59: if (quad->_diff == NULL) PetscCall(VecDuplicate(x, &quad->_diff));
60: PetscCall(VecCopy(x, quad->_diff));
61: PetscCall(VecAXPY(quad->_diff, -1.0, params));
62: *diff = quad->_diff;
63: }
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
67: static PetscErrorCode TaoTermComputeObjective_Quadratic(TaoTerm term, Vec x, Vec params, PetscReal *value)
68: {
69: TaoTerm_Quadratic *quad = (TaoTerm_Quadratic *)term->data;
70: Vec diff;
71: PetscScalar sval;
73: PetscFunctionBegin;
74: PetscCheck(quad->A, PetscObjectComm((PetscObject)term), PETSC_ERR_ORDER, "Quadratic matrix not set, call TaoTermQuadraticSetMat() first");
75: PetscCall(TaoTermQuadraticDiff(term, x, params, &diff));
76: if (quad->Adiff == NULL) PetscCall(VecDuplicate(diff, &quad->Adiff));
77: PetscCall(MatMult(quad->A, diff, quad->Adiff));
78: PetscCall(VecDot(diff, quad->Adiff, &sval));
79: *value = 0.5 * PetscRealPart(sval);
80: PetscFunctionReturn(PETSC_SUCCESS);
81: }
83: static PetscErrorCode TaoTermComputeGradient_Quadratic(TaoTerm term, Vec x, Vec params, Vec g)
84: {
85: TaoTerm_Quadratic *quad = (TaoTerm_Quadratic *)term->data;
86: Vec diff;
88: PetscFunctionBegin;
89: PetscCheck(quad->A, PetscObjectComm((PetscObject)term), PETSC_ERR_ORDER, "Quadratic matrix not set, call TaoTermQuadraticSetMat() first");
90: PetscCall(TaoTermQuadraticDiff(term, x, params, &diff));
91: PetscCall(MatMult(quad->A, diff, g));
92: PetscFunctionReturn(PETSC_SUCCESS);
93: }
95: static PetscErrorCode TaoTermComputeObjectiveAndGradient_Quadratic(TaoTerm term, Vec x, Vec params, PetscReal *value, Vec g)
96: {
97: TaoTerm_Quadratic *quad = (TaoTerm_Quadratic *)term->data;
98: Vec diff;
99: PetscScalar sval;
101: PetscFunctionBegin;
102: PetscCheck(quad->A, PetscObjectComm((PetscObject)term), PETSC_ERR_ORDER, "Quadratic matrix not set, call TaoTermQuadraticSetMat() first");
103: PetscCall(TaoTermQuadraticDiff(term, x, params, &diff));
104: PetscCall(MatMult(quad->A, diff, g));
105: PetscCall(VecDot(diff, g, &sval));
106: *value = 0.5 * PetscRealPart(sval);
107: PetscFunctionReturn(PETSC_SUCCESS);
108: }
110: static PetscErrorCode TaoTermComputeHessian_Quadratic(TaoTerm term, Vec x, Vec params, Mat H, Mat Hpre)
111: {
112: TaoTerm_Quadratic *quad = (TaoTerm_Quadratic *)term->data;
114: PetscFunctionBegin;
115: PetscCheck(quad->A, PetscObjectComm((PetscObject)term), PETSC_ERR_ORDER, "Quadratic matrix not set, call TaoTermQuadraticSetMat() first");
116: // TODO caching to avoid unnecessary computation
117: if (H) PetscCall(MatCopy(quad->A, H, UNKNOWN_NONZERO_PATTERN));
118: if (Hpre && Hpre != H) PetscCall(MatCopy(quad->A, Hpre, UNKNOWN_NONZERO_PATTERN));
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: static PetscErrorCode TaoTermCreateHessianMatrices_Quadratic(TaoTerm term, Mat *H, Mat *Hpre)
123: {
124: TaoTerm_Quadratic *quad = (TaoTerm_Quadratic *)term->data;
126: PetscFunctionBegin;
127: PetscCheck(quad->A, PetscObjectComm((PetscObject)term), PETSC_ERR_ORDER, "Quadratic matrix not set, call TaoTermQuadraticSetMat() first");
128: PetscCall(PetscInfo(term, "%s: Creating TAOTERMQUADRATIC Hessian Matrices by duplicating quadratic matrix set by TaoTermQuadraticSetMat, overriding custom MatType options.\n", ((PetscObject)term)->prefix));
129: if (H) PetscCall(MatDuplicate(quad->A, MAT_DO_NOT_COPY_VALUES, H));
130: if (Hpre) {
131: if (term->Hpre_is_H && H) {
132: PetscCall(PetscObjectReference((PetscObject)*H));
133: *Hpre = *H;
134: } else PetscCall(MatDuplicate(quad->A, MAT_DO_NOT_COPY_VALUES, Hpre));
135: }
136: PetscFunctionReturn(PETSC_SUCCESS);
137: }
139: /*@
140: TaoTermQuadraticGetMat - Get the matrix defining a `TaoTerm` of type `TAOTERMQUADRATIC`
142: Not collective
144: Input Parameter:
145: . term - a `TaoTerm` of type `TAOTERMQUADRATIC`
147: Output Parameter:
148: . A - the matrix
150: Level: intermediate
152: Note:
153: This function will return `NULL` if the term is not a `TAOTERMQUADRATIC`.
155: .seealso: [](sec_tao_term),
156: `TaoTerm`,
157: `TAOTERMQUADRATIC`,
158: `TaoTermQuadraticSetMat()`
159: @*/
160: PetscErrorCode TaoTermQuadraticGetMat(TaoTerm term, Mat *A)
161: {
162: PetscFunctionBegin;
164: PetscAssertPointer(A, 2);
165: *A = NULL;
166: PetscTryMethod(term, "TaoTermQuadraticGetMat_C", (TaoTerm, Mat *), (term, A));
167: PetscFunctionReturn(PETSC_SUCCESS);
168: }
170: static PetscErrorCode TaoTermQuadraticGetMat_Quadratic(TaoTerm term, Mat *A)
171: {
172: TaoTerm_Quadratic *quad = (TaoTerm_Quadratic *)term->data;
174: PetscFunctionBegin;
175: *A = quad->A;
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: /*@
180: TaoTermQuadraticSetMat - Set the matrix defining a `TaoTerm` of type `TAOTERMQUADRATIC`
182: Collective
184: Input Parameters:
185: + term - a `TaoTerm` of type `TAOTERMQUADRATIC`
186: - A - the matrix
188: Level: intermediate
190: .seealso: [](sec_tao_term),
191: `TaoTerm`,
192: `TAOTERMQUADRATIC`,
193: `TaoTermQuadraticGetMat()`
194: @*/
195: PetscErrorCode TaoTermQuadraticSetMat(TaoTerm term, Mat A)
196: {
197: PetscFunctionBegin;
200: PetscCheckSameComm(term, 1, A, 2);
201: PetscTryMethod(term, "TaoTermQuadraticSetMat_C", (TaoTerm, Mat), (term, A));
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
205: static PetscErrorCode TaoTermQuadraticSetMat_Quadratic(TaoTerm term, Mat A)
206: {
207: TaoTerm_Quadratic *quad = (TaoTerm_Quadratic *)term->data;
208: PetscLayout rmap, cmap;
209: VecType vec_type;
210: PetscBool is_square;
212: PetscFunctionBegin;
213: if (A != quad->A) {
214: MatType mat_type;
215: PetscCall(MatGetLayouts(A, &rmap, &cmap));
216: PetscCall(PetscLayoutCompare(rmap, cmap, &is_square));
217: PetscCheck(is_square, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_SIZ, "Matrix for quadratic term must be square");
218: PetscCall(PetscObjectReference((PetscObject)A));
219: PetscCall(MatDestroy(&quad->A));
220: PetscCall(VecDestroy(&quad->_diff));
221: PetscCall(VecDestroy(&quad->Adiff));
222: quad->A = A;
223: PetscCall(MatGetVecType(A, &vec_type));
224: PetscCall(TaoTermSetSolutionVecType(term, vec_type));
225: PetscCall(TaoTermSetParametersVecType(term, vec_type));
226: PetscCall(TaoTermSetSolutionLayout(term, rmap));
227: PetscCall(TaoTermSetParametersLayout(term, rmap));
228: PetscCall(PetscFree(term->H_mattype));
229: PetscCall(PetscFree(term->Hpre_mattype));
230: PetscCall(MatGetType(A, &mat_type));
231: PetscCall(PetscStrallocpy(mat_type, (char **)&term->H_mattype));
232: PetscCall(PetscStrallocpy(mat_type, (char **)&term->Hpre_mattype));
233: }
234: PetscFunctionReturn(PETSC_SUCCESS);
235: }
237: static PetscErrorCode TaoTermIsComputeHessianFDPossible_Quadratic(TaoTerm term, PetscBool3 *ispossible)
238: {
239: PetscFunctionBegin;
240: *ispossible = PETSC_BOOL3_FALSE;
241: PetscFunctionReturn(PETSC_SUCCESS);
242: }
244: /*MC
245: TAOTERMQUADRATIC - A `TaoTerm` that computes $\tfrac{1}{2}(x - p)^T A (x - p)$, for a fixed matrix $A$, solution $x$ and parameters $p$.
247: Level: intermediate
249: Notes:
250: This term is `TAOTERM_PARAMETERS_OPTIONAL`. If the parameters argument is `NULL` for
251: evaluation routines the term computes $\tfrac{1}{2}x^T A x$.
253: The matrix $A$ must be symmetric.
255: The default Hessian creation mode (see `TaoTermGetCreateHessianMode()`) is `H == Hpre` and `TaoTermCreateHessianMatrices()`
256: will create a matrix with the same type as $A$.
258: .seealso: [](sec_tao_term),
259: `TaoTerm`,
260: `TaoTermType`,
261: `TaoTermCreateQuadratic()`,
262: `TAOTERMHALFL2SQUARED`,
263: `TAOTERML1`, `TaoTermQuadraticSetMat()`
264: M*/
265: PETSC_INTERN PetscErrorCode TaoTermCreate_Quadratic(TaoTerm term)
266: {
267: TaoTerm_Quadratic *quad;
269: PetscFunctionBegin;
270: PetscCall(TaoTermCreate_ElementwiseDivergence_Internal(term));
271: PetscCall(PetscNew(&quad));
272: term->data = (void *)quad;
274: PetscCall(PetscFree(term->H_mattype));
275: PetscCall(PetscFree(term->Hpre_mattype));
277: term->ops->destroy = TaoTermDestroy_Quadratic;
278: term->ops->view = TaoTermView_Quadratic;
279: term->ops->objective = TaoTermComputeObjective_Quadratic;
280: term->ops->gradient = TaoTermComputeGradient_Quadratic;
281: term->ops->objectiveandgradient = TaoTermComputeObjectiveAndGradient_Quadratic;
282: term->ops->hessian = TaoTermComputeHessian_Quadratic;
283: term->ops->createhessianmatrices = TaoTermCreateHessianMatrices_Quadratic;
284: term->ops->iscomputehessianfdpossible = TaoTermIsComputeHessianFDPossible_Quadratic;
286: term->Hpre_is_H = PETSC_TRUE;
288: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermQuadraticGetMat_C", TaoTermQuadraticGetMat_Quadratic));
289: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermQuadraticSetMat_C", TaoTermQuadraticSetMat_Quadratic));
290: PetscFunctionReturn(PETSC_SUCCESS);
291: }
293: /*@
294: TaoTermCreateQuadratic - Create a `TAOTERMQUADRATIC` for a given matrix
296: Collective
298: Input Parameter:
299: . A - a square matrix
301: Output Parameter:
302: . term - a `TaoTerm` that implements $\tfrac{1}{2}(x - p)^T A (x - p)$
304: Level: beginner
306: .seealso: [](sec_tao_term),
307: `TaoTerm`,
308: `TaoTermCreate()`,
309: `TAOTERMQUADRATIC`,
310: `TaoTermCreateHalfL2Squared()`,
311: `TaoTermCreateL1()`, `TaoTermQuadraticSetMat()`
312: @*/
313: PetscErrorCode TaoTermCreateQuadratic(Mat A, TaoTerm *term)
314: {
315: TaoTerm _term;
317: PetscFunctionBegin;
319: PetscAssertPointer(term, 2);
320: PetscCall(TaoTermCreate(PetscObjectComm((PetscObject)A), &_term));
321: PetscCall(TaoTermSetType(_term, TAOTERMQUADRATIC));
322: PetscCall(TaoTermQuadraticSetMat(_term, A));
323: *term = _term;
324: PetscFunctionReturn(PETSC_SUCCESS);
325: }