Actual source code: taosolver_hj.c
1: #include <petsc/private/taoimpl.h>
3: /*@C
4: TaoSetHessian - Sets the function to compute the Hessian as well as the location to store the matrix.
6: Logically Collective
8: Input Parameters:
9: + tao - the `Tao` context
10: . H - Matrix used for the hessian
11: . Hpre - Matrix that will be used to construct the preconditioner, can be same as `H`
12: . func - Hessian evaluation routine
13: - ctx - [optional] user-defined context for private data for the
14: Hessian evaluation routine (may be `NULL`)
16: Calling sequence of `func`:
17: + tao - the `Tao` context
18: . x - input vector
19: . H - Hessian matrix
20: . Hpre - matrix used to construct the preconditioner, usually the same as `H`
21: - ctx - [optional] user-defined Hessian context
23: Level: beginner
25: .seealso: [](ch_tao), `Tao`, `TaoType`, `TaoSetObjective()`, `TaoSetGradient()`, `TaoSetObjectiveAndGradient()`, `TaoGetHessian()`
26: @*/
27: PetscErrorCode TaoSetHessian(Tao tao, Mat H, Mat Hpre, PetscErrorCode (*func)(Tao tao, Vec x, Mat H, Mat Hpre, PetscCtx ctx), PetscCtx ctx)
28: {
29: PetscFunctionBegin;
31: if (H) {
33: PetscCheckSameComm(tao, 1, H, 2);
34: }
35: if (Hpre) {
37: PetscCheckSameComm(tao, 1, Hpre, 3);
38: }
39: PetscCall(TaoTermCallbacksSetHessian(tao->callbacks, func, ctx));
40: if (H) {
41: PetscCall(PetscObjectReference((PetscObject)H));
42: PetscCall(MatDestroy(&tao->hessian));
43: tao->hessian = H;
44: }
45: if (Hpre) {
46: PetscCall(PetscObjectReference((PetscObject)Hpre));
47: PetscCall(MatDestroy(&tao->hessian_pre));
48: tao->hessian_pre = Hpre;
49: }
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: /*@C
54: TaoGetHessian - Gets the function to compute the Hessian as well as the location to store the matrix.
56: Not Collective
58: Input Parameter:
59: . tao - the `Tao` context
61: Output Parameters:
62: + H - Matrix used for the hessian
63: . Hpre - Matrix that will be used to construct the preconditioner, can be the same as `H`
64: . func - Hessian evaluation routine
65: - ctx - user-defined context for private data for the Hessian evaluation routine
67: Calling sequence of `func`:
68: + tao - the `Tao` context
69: . x - input vector
70: . H - Hessian matrix
71: . Hpre - matrix used to construct the preconditioner, usually the same as `H`
72: - ctx - [optional] user-defined Hessian context
74: Level: beginner
76: Notes:
77: In addition to specifying an objective function using callbacks such as
78: `TaoSetObjectiveAndGradient()` and `TaoSetHessian()`, users can specify
79: objective functions with `TaoAddTerm()`.
81: `TaoGetHessian()` will always return the callback specified with
82: `TaoSetHessian()`, even if the objective function has been changed by
83: calling `TaoAddTerm()`.
85: .seealso: [](ch_tao), `Tao`, `TaoType`, `TaoGetObjective()`, `TaoGetGradient()`, `TaoGetObjectiveAndGradient()`, `TaoSetHessian()`, `TaoGetHessianMatrices()`
86: @*/
87: PetscErrorCode TaoGetHessian(Tao tao, Mat *H, Mat *Hpre, PetscErrorCode (**func)(Tao tao, Vec x, Mat H, Mat Hpre, PetscCtx ctx), PetscCtxRt ctx)
88: {
89: PetscFunctionBegin;
91: PetscCall(TaoGetHessianMatrices(tao, H, Hpre));
92: if (func || ctx) PetscCall(TaoTermCallbacksGetHessian(tao->callbacks, func, ctx));
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
96: /*@
97: TaoGetHessianMatrices - Get the matrices that store the Hessian matrix and its (optional) approximation that is used to construct the preconditioner
99: Not collective
101: Input Parameter:
102: . tao - the `Tao` context
104: Output Parameters:
105: + H - the Hessian matrix
106: - Hpre - approximation to the Hessian matrix used to construct the preconditioner (often `H`)
108: Level: intermediate
110: .seealso: [](ch_tao), `Tao`, `TaoType`, `TaoGetObjective()`, `TaoGetGradient()`, `TaoGetObjectiveAndGradient()`, `TaoSetHessian()`, `TaoGetHessian()`
111: @*/
112: PetscErrorCode TaoGetHessianMatrices(Tao tao, Mat *H, Mat *Hpre)
113: {
114: PetscFunctionBegin;
116: if (H) *H = tao->hessian;
117: if (Hpre) *Hpre = tao->hessian_pre;
118: PetscFunctionReturn(PETSC_SUCCESS);
119: }
121: /*@
122: TaoTestHessian - Compare the user-supplied Hessian with a finite-difference approximation, when requested via
123: the options database, and print the difference.
125: Collective
127: Input Parameter:
128: . tao - the `Tao` context
130: Options Database Keys:
131: + -tao_test_hessian threshold - enable the comparison, optionally overriding the reporting threshold (default `1e-5`)
132: - -tao_test_hessian_view - display the user-supplied Hessian, the finite-difference Hessian, and their difference
134: Level: intermediate
136: Note:
137: If `-tao_test_hessian` is not set, this routine returns immediately without performing any work.
139: .seealso: [](ch_tao), `Tao`, `TaoTestGradient()`, `TaoComputeHessian()`
140: @*/
141: PetscErrorCode TaoTestHessian(Tao tao)
142: {
143: Mat A, B, C, D, hessian;
144: Vec x = tao->solution;
145: PetscReal nrm, gnorm;
146: PetscReal threshold = 1.e-5;
147: PetscBool complete_print = PETSC_FALSE, test = PETSC_FALSE, flg;
148: PetscViewer viewer, mviewer;
149: MPI_Comm comm;
150: PetscInt tabs;
151: static PetscBool directionsprinted = PETSC_FALSE;
152: PetscViewerFormat format;
154: PetscFunctionBegin;
155: PetscObjectOptionsBegin((PetscObject)tao);
156: PetscCall(PetscOptionsName("-tao_test_hessian", "Compare hand-coded and finite difference Hessians", "None", &test));
157: PetscCall(PetscOptionsReal("-tao_test_hessian", "Threshold for element difference between hand-coded and finite difference being meaningful", "None", threshold, &threshold, NULL));
158: PetscCall(PetscOptionsViewer("-tao_test_hessian_view", "View difference between hand-coded and finite difference Hessians element entries", "None", &mviewer, &format, &complete_print));
159: PetscOptionsEnd();
160: if (!test) PetscFunctionReturn(PETSC_SUCCESS);
162: PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
163: PetscCall(PetscViewerASCIIGetStdout(comm, &viewer));
164: PetscCall(PetscViewerASCIIGetTab(viewer, &tabs));
165: PetscCall(PetscViewerASCIISetTab(viewer, ((PetscObject)tao)->tablevel));
166: PetscCall(PetscViewerASCIIPrintf(viewer, " ---------- Testing Hessian -------------\n"));
167: if (!complete_print && !directionsprinted) {
168: PetscCall(PetscViewerASCIIPrintf(viewer, " Run with -tao_test_hessian_view and optionally -tao_test_hessian <threshold> to show difference\n"));
169: PetscCall(PetscViewerASCIIPrintf(viewer, " of hand-coded and finite difference Hessian entries greater than <threshold>.\n"));
170: }
171: if (!directionsprinted) {
172: PetscCall(PetscViewerASCIIPrintf(viewer, " Testing hand-coded Hessian, if (for double precision runs) ||J - Jfd||_F/||J||_F is\n"));
173: PetscCall(PetscViewerASCIIPrintf(viewer, " O(1.e-8), the hand-coded Hessian is probably correct.\n"));
174: directionsprinted = PETSC_TRUE;
175: }
176: if (complete_print) PetscCall(PetscViewerPushFormat(mviewer, format));
178: PetscCall(PetscObjectTypeCompare((PetscObject)tao->hessian, MATMFFD, &flg));
179: if (!flg) hessian = tao->hessian;
180: else hessian = tao->hessian_pre;
182: while (hessian) {
183: PetscLayout rmap, cmap;
184: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)hessian, &flg, MATSEQAIJ, MATMPIAIJ, MATSEQDENSE, MATMPIDENSE, MATSEQBAIJ, MATMPIBAIJ, MATSEQSBAIJ, MATMPISBAIJ, ""));
185: if (flg) {
186: A = hessian;
187: PetscCall(PetscObjectReference((PetscObject)A));
188: } else {
189: PetscCall(MatComputeOperator(hessian, MATAIJ, &A));
190: }
192: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
193: PetscCall(MatGetLayouts(A, &rmap, &cmap));
194: PetscCall(MatSetLayouts(B, rmap, cmap));
195: PetscCall(MatSetType(B, ((PetscObject)A)->type_name));
196: PetscCall(MatSetUp(B));
197: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
199: PetscCall(TaoDefaultComputeHessian(tao, x, B, B, NULL));
201: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &D));
202: PetscCall(MatAYPX(D, -1.0, A, DIFFERENT_NONZERO_PATTERN));
203: PetscCall(MatNorm(D, NORM_FROBENIUS, &nrm));
204: PetscCall(MatNorm(A, NORM_FROBENIUS, &gnorm));
205: PetscCall(MatDestroy(&D));
206: if (!gnorm) gnorm = 1; /* just in case */
207: PetscCall(PetscViewerASCIIPrintf(viewer, " ||H - Hfd||_F/||H||_F = %g, ||H - Hfd||_F = %g\n", (double)(nrm / gnorm), (double)nrm));
209: if (complete_print) {
210: PetscCall(PetscViewerASCIIPrintf(viewer, " Hand-coded Hessian ----------\n"));
211: PetscCall(MatView(A, mviewer));
212: PetscCall(PetscViewerASCIIPrintf(viewer, " Finite difference Hessian ----------\n"));
213: PetscCall(MatView(B, mviewer));
214: }
216: if (complete_print) {
217: PetscInt Istart, Iend, *ccols, bncols, cncols, j, row;
218: PetscScalar *cvals;
219: const PetscInt *bcols;
220: const PetscScalar *bvals;
222: PetscCall(MatAYPX(B, -1.0, A, DIFFERENT_NONZERO_PATTERN));
223: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
224: PetscCall(MatSetLayouts(C, rmap, cmap));
225: PetscCall(MatSetType(C, ((PetscObject)A)->type_name));
226: PetscCall(MatSetUp(C));
227: PetscCall(MatSetOption(C, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
228: PetscCall(MatGetOwnershipRange(B, &Istart, &Iend));
230: for (row = Istart; row < Iend; row++) {
231: PetscCall(MatGetRow(B, row, &bncols, &bcols, &bvals));
232: PetscCall(PetscMalloc2(bncols, &ccols, bncols, &cvals));
233: for (j = 0, cncols = 0; j < bncols; j++) {
234: if (PetscAbsScalar(bvals[j]) > threshold) {
235: ccols[cncols] = bcols[j];
236: cvals[cncols] = bvals[j];
237: cncols += 1;
238: }
239: }
240: if (cncols) PetscCall(MatSetValues(C, 1, &row, cncols, ccols, cvals, INSERT_VALUES));
241: PetscCall(MatRestoreRow(B, row, &bncols, &bcols, &bvals));
242: PetscCall(PetscFree2(ccols, cvals));
243: }
244: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
245: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
246: PetscCall(PetscViewerASCIIPrintf(viewer, " Finite-difference minus hand-coded Hessian with tolerance %g ----------\n", (double)threshold));
247: PetscCall(MatView(C, mviewer));
248: PetscCall(MatDestroy(&C));
249: }
250: PetscCall(MatDestroy(&A));
251: PetscCall(MatDestroy(&B));
253: if (hessian != tao->hessian_pre) {
254: hessian = tao->hessian_pre;
255: PetscCall(PetscViewerASCIIPrintf(viewer, " ---------- Testing Hessian for preconditioner -------------\n"));
256: } else hessian = NULL;
257: }
258: if (complete_print) {
259: PetscCall(PetscViewerPopFormat(mviewer));
260: PetscCall(PetscViewerDestroy(&mviewer));
261: }
262: PetscCall(PetscViewerASCIISetTab(viewer, tabs));
263: PetscFunctionReturn(PETSC_SUCCESS);
264: }
266: /*@
267: TaoComputeHessian - Computes the Hessian matrix that has been
268: set with `TaoSetHessian()`.
270: Collective
272: Input Parameters:
273: + tao - the Tao solver context
274: - X - input vector
276: Output Parameters:
277: + H - Hessian matrix
278: - Hpre - matrix used to construct the preconditioner, usually the same as `H`
280: Options Database Keys:
281: + -tao_test_hessian - compare the user provided Hessian with one compute via finite differences to check for errors
282: . -tao_test_hessian numerical value - display entries in the difference between the user provided Hessian and finite difference Hessian that are greater than a certain value to help users detect errors
283: - -tao_test_hessian_view - display the user provided Hessian, the finite difference Hessian and the difference between them to help users detect the location of errors in the user provided Hessian
285: Level: developer
287: Notes:
288: Most users should not need to explicitly call this routine, as it
289: is used internally within the minimization solvers.
291: `TaoComputeHessian()` is typically used within optimization algorithms,
292: so most users would not generally call this routine
293: themselves.
295: Developer Notes:
296: The Hessian test mechanism follows `SNESTestJacobian()`.
298: If there is no separate preconditioning matrix, `TaoComputeHessian(tao, X, H, NULL)` and `TaoComputeHessian(tao, X, H, H)` are equivalent.
300: .seealso: [](ch_tao), `Tao`, `TaoComputeObjective()`, `TaoComputeObjectiveAndGradient()`, `TaoSetHessian()`
301: @*/
302: PetscErrorCode TaoComputeHessian(Tao tao, Vec X, Mat H, Mat Hpre)
303: {
304: PetscFunctionBegin;
307: PetscCheckSameComm(tao, 1, X, 2);
308: PetscCall(TaoTermMappingComputeHessian(&tao->objective_term, X, tao->objective_parameters, INSERT_VALUES, H, Hpre));
309: PetscCall(TaoTestHessian(tao));
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: /*@
314: TaoComputeJacobian - Computes the Jacobian matrix that has been
315: set with TaoSetJacobianRoutine().
317: Collective
319: Input Parameters:
320: + tao - the Tao solver context
321: - X - input vector
323: Output Parameters:
324: + J - Jacobian matrix
325: - Jpre - matrix used to compute the preconditioner, often the same as `J`
327: Level: developer
329: Notes:
330: Most users should not need to explicitly call this routine, as it
331: is used internally within the minimization solvers.
333: `TaoComputeJacobian()` is typically used within minimization
334: implementations, so most users would not generally call this routine
335: themselves.
337: .seealso: [](ch_tao), `TaoComputeObjective()`, `TaoComputeObjectiveAndGradient()`, `TaoSetJacobianRoutine()`
338: @*/
339: PetscErrorCode TaoComputeJacobian(Tao tao, Vec X, Mat J, Mat Jpre)
340: {
341: PetscFunctionBegin;
344: PetscCheckSameComm(tao, 1, X, 2);
345: ++tao->njac;
346: PetscCall(VecLockReadPush(X));
347: PetscCall(PetscLogEventBegin(TAO_JacobianEval, tao, X, J, Jpre));
348: PetscCallBack("Tao callback Jacobian", (*tao->ops->computejacobian)(tao, X, J, Jpre, tao->user_jacP));
349: PetscCall(PetscLogEventEnd(TAO_JacobianEval, tao, X, J, Jpre));
350: PetscCall(VecLockReadPop(X));
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }
354: /*@
355: TaoComputeResidualJacobian - Computes the least-squares residual Jacobian matrix that has been
356: set with `TaoSetJacobianResidual()`.
358: Collective
360: Input Parameters:
361: + tao - the Tao solver context
362: - X - input vector
364: Output Parameters:
365: + J - Jacobian matrix
366: - Jpre - matrix used to compute the preconditioner, often the same as `J`
368: Level: developer
370: Notes:
371: Most users should not need to explicitly call this routine, as it
372: is used internally within the minimization solvers.
374: `TaoComputeResidualJacobian()` is typically used within least-squares
375: implementations, so most users would not generally call this routine
376: themselves.
378: .seealso: [](ch_tao), `Tao`, `TaoComputeResidual()`, `TaoSetJacobianResidual()`
379: @*/
380: PetscErrorCode TaoComputeResidualJacobian(Tao tao, Vec X, Mat J, Mat Jpre)
381: {
382: PetscFunctionBegin;
385: PetscCheckSameComm(tao, 1, X, 2);
386: ++tao->njac;
387: PetscCall(VecLockReadPush(X));
388: PetscCall(PetscLogEventBegin(TAO_JacobianEval, tao, X, J, Jpre));
389: PetscCallBack("Tao callback least-squares residual Jacobian", (*tao->ops->computeresidualjacobian)(tao, X, J, Jpre, tao->user_lsjacP));
390: PetscCall(PetscLogEventEnd(TAO_JacobianEval, tao, X, J, Jpre));
391: PetscCall(VecLockReadPop(X));
392: PetscFunctionReturn(PETSC_SUCCESS);
393: }
395: /*@
396: TaoComputeJacobianState - Computes the Jacobian matrix that has been
397: set with `TaoSetJacobianStateRoutine()`.
399: Collective
401: Input Parameters:
402: + tao - the `Tao` solver context
403: - X - input vector
405: Output Parameters:
406: + J - Jacobian matrix
407: . Jpre - matrix used to construct the preconditioner, often the same as `J`
408: - Jinv - unknown
410: Level: developer
412: Note:
413: Most users should not need to explicitly call this routine, as it
414: is used internally within the optimization algorithms.
416: .seealso: [](ch_tao), `Tao`, `TaoComputeObjective()`, `TaoComputeObjectiveAndGradient()`, `TaoSetJacobianStateRoutine()`, `TaoComputeJacobianDesign()`, `TaoSetStateDesignIS()`
417: @*/
418: PetscErrorCode TaoComputeJacobianState(Tao tao, Vec X, Mat J, Mat Jpre, Mat Jinv)
419: {
420: PetscFunctionBegin;
423: PetscCheckSameComm(tao, 1, X, 2);
424: ++tao->njac_state;
425: PetscCall(VecLockReadPush(X));
426: PetscCall(PetscLogEventBegin(TAO_JacobianEval, tao, X, J, Jpre));
427: PetscCallBack("Tao callback Jacobian(state)", (*tao->ops->computejacobianstate)(tao, X, J, Jpre, Jinv, tao->user_jac_stateP));
428: PetscCall(PetscLogEventEnd(TAO_JacobianEval, tao, X, J, Jpre));
429: PetscCall(VecLockReadPop(X));
430: PetscFunctionReturn(PETSC_SUCCESS);
431: }
433: /*@
434: TaoComputeJacobianDesign - Computes the Jacobian matrix that has been
435: set with `TaoSetJacobianDesignRoutine()`.
437: Collective
439: Input Parameters:
440: + tao - the Tao solver context
441: - X - input vector
443: Output Parameter:
444: . J - Jacobian matrix
446: Level: developer
448: Note:
449: Most users should not need to explicitly call this routine, as it
450: is used internally within the optimization algorithms.
452: .seealso: [](ch_tao), `Tao`, `TaoComputeObjective()`, `TaoComputeObjectiveAndGradient()`, `TaoSetJacobianDesignRoutine()`, `TaoSetStateDesignIS()`
453: @*/
454: PetscErrorCode TaoComputeJacobianDesign(Tao tao, Vec X, Mat J)
455: {
456: PetscFunctionBegin;
459: PetscCheckSameComm(tao, 1, X, 2);
460: ++tao->njac_design;
461: PetscCall(VecLockReadPush(X));
462: PetscCall(PetscLogEventBegin(TAO_JacobianEval, tao, X, J, NULL));
463: PetscCallBack("Tao callback Jacobian(design)", (*tao->ops->computejacobiandesign)(tao, X, J, tao->user_jac_designP));
464: PetscCall(PetscLogEventEnd(TAO_JacobianEval, tao, X, J, NULL));
465: PetscCall(VecLockReadPop(X));
466: PetscFunctionReturn(PETSC_SUCCESS);
467: }
469: /*@C
470: TaoSetJacobianRoutine - Sets the function to compute the Jacobian as well as the location to store the matrix.
472: Logically Collective
474: Input Parameters:
475: + tao - the `Tao` context
476: . J - Matrix used for the Jacobian
477: . Jpre - Matrix that will be used to construct the preconditioner, can be same as `J`
478: . func - Jacobian evaluation routine
479: - ctx - [optional] user-defined context for private data for the
480: Jacobian evaluation routine (may be `NULL`)
482: Calling sequence of `func`:
483: + tao - the `Tao` context
484: . x - input vector
485: . J - Jacobian matrix
486: . Jpre - matrix used to construct the preconditioner, usually the same as `J`
487: - ctx - [optional] user-defined Jacobian context
489: Level: intermediate
491: .seealso: [](ch_tao), `Tao`, `TaoSetGradient()`, `TaoSetObjective()`
492: @*/
493: PetscErrorCode TaoSetJacobianRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao tao, Vec x, Mat J, Mat Jpre, PetscCtx ctx), PetscCtx ctx)
494: {
495: PetscFunctionBegin;
497: if (J) {
499: PetscCheckSameComm(tao, 1, J, 2);
500: }
501: if (Jpre) {
503: PetscCheckSameComm(tao, 1, Jpre, 3);
504: }
505: if (ctx) tao->user_jacP = ctx;
506: if (func) tao->ops->computejacobian = func;
507: if (J) {
508: PetscCall(PetscObjectReference((PetscObject)J));
509: PetscCall(MatDestroy(&tao->jacobian));
510: tao->jacobian = J;
511: }
512: if (Jpre) {
513: PetscCall(PetscObjectReference((PetscObject)Jpre));
514: PetscCall(MatDestroy(&tao->jacobian_pre));
515: tao->jacobian_pre = Jpre;
516: }
517: PetscFunctionReturn(PETSC_SUCCESS);
518: }
520: /*@C
521: TaoSetJacobianResidualRoutine - Sets the function to compute the least-squares residual Jacobian as well as the
522: location to store the matrix.
524: Logically Collective
526: Input Parameters:
527: + tao - the `Tao` context
528: . J - Matrix used for the jacobian
529: . Jpre - Matrix that will be used to construct the preconditioner, can be same as `J`
530: . func - Jacobian evaluation routine
531: - ctx - [optional] user-defined context for private data for the
532: Jacobian evaluation routine (may be `NULL`)
534: Calling sequence of `func`:
535: + tao - the `Tao` context
536: . x - input vector
537: . J - Jacobian matrix
538: . Jpre - matrix used to construct the preconditioner, usually the same as `J`
539: - ctx - [optional] user-defined Jacobian context
541: Level: intermediate
543: .seealso: [](ch_tao), `Tao`, `TaoSetGradient()`, `TaoSetObjective()`
544: @*/
545: PetscErrorCode TaoSetJacobianResidualRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao tao, Vec x, Mat J, Mat Jpre, PetscCtx ctx), PetscCtx ctx)
546: {
547: PetscFunctionBegin;
549: if (J) {
551: PetscCheckSameComm(tao, 1, J, 2);
552: }
553: if (Jpre) {
555: PetscCheckSameComm(tao, 1, Jpre, 3);
556: }
557: if (ctx) tao->user_lsjacP = ctx;
558: if (func) tao->ops->computeresidualjacobian = func;
559: if (J) {
560: PetscCall(PetscObjectReference((PetscObject)J));
561: PetscCall(MatDestroy(&tao->ls_jac));
562: tao->ls_jac = J;
563: }
564: if (Jpre) {
565: PetscCall(PetscObjectReference((PetscObject)Jpre));
566: PetscCall(MatDestroy(&tao->ls_jac_pre));
567: tao->ls_jac_pre = Jpre;
568: }
569: PetscFunctionReturn(PETSC_SUCCESS);
570: }
572: /*@C
573: TaoSetJacobianStateRoutine - Sets the function to compute the Jacobian
574: (and its inverse) of the constraint function with respect to the state variables.
575: Used only for PDE-constrained optimization.
577: Logically Collective
579: Input Parameters:
580: + tao - the `Tao` context
581: . J - Matrix used for the Jacobian
582: . Jpre - Matrix that will be used to construct the preconditioner, can be same as `J`. Only used if `Jinv` is `NULL`
583: . Jinv - [optional] Matrix used to apply the inverse of the state Jacobian. Use `NULL` to default to PETSc `KSP` solvers to apply the inverse.
584: . func - Jacobian evaluation routine
585: - ctx - [optional] user-defined context for private data for the
586: Jacobian evaluation routine (may be `NULL`)
588: Calling sequence of `func`:
589: + tao - the `Tao` context
590: . x - input vector
591: . J - Jacobian matrix
592: . Jpre - matrix used to construct the preconditioner, usually the same as `J`
593: . Jinv - inverse of `J`
594: - ctx - [optional] user-defined Jacobian context
596: Level: intermediate
598: .seealso: [](ch_tao), `Tao`, `TaoComputeJacobianState()`, `TaoSetJacobianDesignRoutine()`, `TaoSetStateDesignIS()`
599: @*/
600: PetscErrorCode TaoSetJacobianStateRoutine(Tao tao, Mat J, Mat Jpre, Mat Jinv, PetscErrorCode (*func)(Tao tao, Vec x, Mat J, Mat Jpre, Mat Jinv, PetscCtx ctx), PetscCtx ctx)
601: {
602: PetscFunctionBegin;
604: if (J) {
606: PetscCheckSameComm(tao, 1, J, 2);
607: }
608: if (Jpre) {
610: PetscCheckSameComm(tao, 1, Jpre, 3);
611: }
612: if (Jinv) {
614: PetscCheckSameComm(tao, 1, Jinv, 4);
615: }
616: if (ctx) tao->user_jac_stateP = ctx;
617: if (func) tao->ops->computejacobianstate = func;
618: if (J) {
619: PetscCall(PetscObjectReference((PetscObject)J));
620: PetscCall(MatDestroy(&tao->jacobian_state));
621: tao->jacobian_state = J;
622: }
623: if (Jpre) {
624: PetscCall(PetscObjectReference((PetscObject)Jpre));
625: PetscCall(MatDestroy(&tao->jacobian_state_pre));
626: tao->jacobian_state_pre = Jpre;
627: }
628: if (Jinv) {
629: PetscCall(PetscObjectReference((PetscObject)Jinv));
630: PetscCall(MatDestroy(&tao->jacobian_state_inv));
631: tao->jacobian_state_inv = Jinv;
632: }
633: PetscFunctionReturn(PETSC_SUCCESS);
634: }
636: /*@C
637: TaoSetJacobianDesignRoutine - Sets the function to compute the Jacobian of
638: the constraint function with respect to the design variables. Used only for
639: PDE-constrained optimization.
641: Logically Collective
643: Input Parameters:
644: + tao - the `Tao` context
645: . J - Matrix used for the Jacobian
646: . func - Jacobian evaluation routine
647: - ctx - [optional] user-defined context for private data for the
648: Jacobian evaluation routine (may be `NULL`)
650: Calling sequence of `func`:
651: + tao - the `Tao` context
652: . x - input vector
653: . J - Jacobian matrix
654: - ctx - [optional] user-defined Jacobian context
656: Level: intermediate
658: .seealso: [](ch_tao), `Tao`, `TaoComputeJacobianDesign()`, `TaoSetJacobianStateRoutine()`, `TaoSetStateDesignIS()`
659: @*/
660: PetscErrorCode TaoSetJacobianDesignRoutine(Tao tao, Mat J, PetscErrorCode (*func)(Tao tao, Vec x, Mat J, PetscCtx ctx), PetscCtx ctx)
661: {
662: PetscFunctionBegin;
664: if (J) {
666: PetscCheckSameComm(tao, 1, J, 2);
667: }
668: if (ctx) tao->user_jac_designP = ctx;
669: if (func) tao->ops->computejacobiandesign = func;
670: if (J) {
671: PetscCall(PetscObjectReference((PetscObject)J));
672: PetscCall(MatDestroy(&tao->jacobian_design));
673: tao->jacobian_design = J;
674: }
675: PetscFunctionReturn(PETSC_SUCCESS);
676: }
678: /*@
679: TaoSetStateDesignIS - Indicate to the `Tao` object which variables in the
680: solution vector are state variables and which are design. Only applies to
681: PDE-constrained optimization.
683: Logically Collective
685: Input Parameters:
686: + tao - The `Tao` context
687: . s_is - the index set corresponding to the state variables
688: - d_is - the index set corresponding to the design variables
690: Level: intermediate
692: .seealso: [](ch_tao), `Tao`, `TaoSetJacobianStateRoutine()`, `TaoSetJacobianDesignRoutine()`
693: @*/
694: PetscErrorCode TaoSetStateDesignIS(Tao tao, IS s_is, IS d_is)
695: {
696: PetscFunctionBegin;
697: PetscCall(PetscObjectReference((PetscObject)s_is));
698: PetscCall(ISDestroy(&tao->state_is));
699: tao->state_is = s_is;
700: PetscCall(PetscObjectReference((PetscObject)d_is));
701: PetscCall(ISDestroy(&tao->design_is));
702: tao->design_is = d_is;
703: PetscFunctionReturn(PETSC_SUCCESS);
704: }
706: /*@
707: TaoComputeJacobianEquality - Computes the Jacobian matrix that has been
708: set with `TaoSetJacobianEqualityRoutine()`.
710: Collective
712: Input Parameters:
713: + tao - the `Tao` solver context
714: - X - input vector
716: Output Parameters:
717: + J - Jacobian matrix
718: - Jpre - matrix used to construct the preconditioner, often the same as `J`
720: Level: developer
722: Notes:
723: Most users should not need to explicitly call this routine, as it
724: is used internally within the optimization algorithms.
726: .seealso: [](ch_tao), `TaoComputeObjective()`, `TaoComputeObjectiveAndGradient()`, `TaoSetJacobianStateRoutine()`, `TaoComputeJacobianDesign()`, `TaoSetStateDesignIS()`
727: @*/
728: PetscErrorCode TaoComputeJacobianEquality(Tao tao, Vec X, Mat J, Mat Jpre)
729: {
730: PetscFunctionBegin;
733: PetscCheckSameComm(tao, 1, X, 2);
734: ++tao->njac_equality;
735: PetscCall(VecLockReadPush(X));
736: PetscCall(PetscLogEventBegin(TAO_JacobianEval, tao, X, J, Jpre));
737: PetscCallBack("Tao callback Jacobian(equality)", (*tao->ops->computejacobianequality)(tao, X, J, Jpre, tao->user_jac_equalityP));
738: PetscCall(PetscLogEventEnd(TAO_JacobianEval, tao, X, J, Jpre));
739: PetscCall(VecLockReadPop(X));
740: PetscFunctionReturn(PETSC_SUCCESS);
741: }
743: /*@
744: TaoComputeJacobianInequality - Computes the Jacobian matrix that has been
745: set with `TaoSetJacobianInequalityRoutine()`.
747: Collective
749: Input Parameters:
750: + tao - the `Tao` solver context
751: - X - input vector
753: Output Parameters:
754: + J - Jacobian matrix
755: - Jpre - matrix used to construct the preconditioner
757: Level: developer
759: Note:
760: Most users should not need to explicitly call this routine, as it
761: is used internally within the minimization solvers.
763: .seealso: [](ch_tao), `Tao`, `TaoComputeObjective()`, `TaoComputeObjectiveAndGradient()`, `TaoSetJacobianStateRoutine()`, `TaoComputeJacobianDesign()`, `TaoSetStateDesignIS()`
764: @*/
765: PetscErrorCode TaoComputeJacobianInequality(Tao tao, Vec X, Mat J, Mat Jpre)
766: {
767: PetscFunctionBegin;
770: PetscCheckSameComm(tao, 1, X, 2);
771: ++tao->njac_inequality;
772: PetscCall(VecLockReadPush(X));
773: PetscCall(PetscLogEventBegin(TAO_JacobianEval, tao, X, J, Jpre));
774: PetscCallBack("Tao callback Jacobian (inequality)", (*tao->ops->computejacobianinequality)(tao, X, J, Jpre, tao->user_jac_inequalityP));
775: PetscCall(PetscLogEventEnd(TAO_JacobianEval, tao, X, J, Jpre));
776: PetscCall(VecLockReadPop(X));
777: PetscFunctionReturn(PETSC_SUCCESS);
778: }
780: /*@C
781: TaoSetJacobianEqualityRoutine - Sets the function to compute the Jacobian
782: (and its inverse) of the constraint function with respect to the equality variables.
783: Used only for PDE-constrained optimization.
785: Logically Collective
787: Input Parameters:
788: + tao - the `Tao` context
789: . J - Matrix used for the Jacobian
790: . Jpre - Matrix that will be used to construct the preconditioner, can be same as `J`.
791: . func - Jacobian evaluation routine
792: - ctx - [optional] user-defined context for private data for the
793: Jacobian evaluation routine (may be `NULL`)
795: Calling sequence of `func`:
796: + tao - the `Tao` context
797: . x - input vector
798: . J - Jacobian matrix
799: . Jpre - matrix used to construct the preconditioner, usually the same as `J`
800: - ctx - [optional] user-defined Jacobian context
802: Level: intermediate
804: .seealso: [](ch_tao), `Tao`, `TaoComputeJacobianEquality()`, `TaoSetJacobianDesignRoutine()`, `TaoSetEqualityDesignIS()`
805: @*/
806: PetscErrorCode TaoSetJacobianEqualityRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao tao, Vec x, Mat J, Mat Jpre, PetscCtx ctx), PetscCtx ctx)
807: {
808: PetscFunctionBegin;
810: if (J) {
812: PetscCheckSameComm(tao, 1, J, 2);
813: }
814: if (Jpre) {
816: PetscCheckSameComm(tao, 1, Jpre, 3);
817: }
818: if (ctx) tao->user_jac_equalityP = ctx;
819: if (func) tao->ops->computejacobianequality = func;
820: if (J) {
821: PetscCall(PetscObjectReference((PetscObject)J));
822: PetscCall(MatDestroy(&tao->jacobian_equality));
823: tao->jacobian_equality = J;
824: }
825: if (Jpre) {
826: PetscCall(PetscObjectReference((PetscObject)Jpre));
827: PetscCall(MatDestroy(&tao->jacobian_equality_pre));
828: tao->jacobian_equality_pre = Jpre;
829: }
830: PetscFunctionReturn(PETSC_SUCCESS);
831: }
833: /*@C
834: TaoGetJacobianEqualityRoutine - Gets the function used to compute equality constraint Jacobian.
836: Not Collective
838: Input Parameter:
839: . tao - the `Tao` context
841: Output Parameters:
842: + J - the matrix to internally hold the constraint computation
843: . Jpre - the matrix used to construct the preconditioner
844: . func - Jacobian evaluation routine
845: - ctx - the (optional) user-defined context
847: Calling sequence of `func`:
848: + tao - the `Tao` context
849: . x - input vector
850: . J - Jacobian matrix
851: . Jpre - matrix used to construct the preconditioner, usually the same as `J`
852: - ctx - [optional] user-defined Jacobian context
854: Level: intermediate
856: .seealso: [](ch_tao), `Tao`, `TaoComputeJacobianEquality()`, `TaoSetJacobianEqualityRoutine()`
857: @*/
858: PetscErrorCode TaoGetJacobianEqualityRoutine(Tao tao, Mat *J, Mat *Jpre, PetscErrorCode (**func)(Tao tao, Vec x, Mat J, Mat Jpre, PetscCtx ctx), PetscCtxRt ctx)
859: {
860: PetscFunctionBegin;
862: if (J) *J = tao->jacobian_equality;
863: if (Jpre) *Jpre = tao->jacobian_equality_pre;
864: if (func) *func = tao->ops->computejacobianequality;
865: if (ctx) *(void **)ctx = tao->user_jac_equalityP;
866: PetscFunctionReturn(PETSC_SUCCESS);
867: }
869: /*@C
870: TaoSetJacobianInequalityRoutine - Sets the function to compute the Jacobian
871: (and its inverse) of the constraint function with respect to the inequality variables.
872: Used only for PDE-constrained optimization.
874: Logically Collective
876: Input Parameters:
877: + tao - the `Tao` context
878: . J - Matrix used for the Jacobian
879: . Jpre - Matrix that will be used to construct the preconditioner, can be same as `J`.
880: . func - Jacobian evaluation routine
881: - ctx - [optional] user-defined context for private data for the
882: Jacobian evaluation routine (may be `NULL`)
884: Calling sequence of `func`:
885: + tao - the `Tao` context
886: . x - input vector
887: . J - Jacobian matrix
888: . Jpre - matrix used to construct the preconditioner, usually the same as `J`
889: - ctx - [optional] user-defined Jacobian context
891: Level: intermediate
893: .seealso: [](ch_tao), `Tao`, `TaoComputeJacobianInequality()`, `TaoSetJacobianDesignRoutine()`, `TaoSetInequalityDesignIS()`
894: @*/
895: PetscErrorCode TaoSetJacobianInequalityRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao tao, Vec x, Mat J, Mat Jpre, PetscCtx ctx), PetscCtx ctx)
896: {
897: PetscFunctionBegin;
899: if (J) {
901: PetscCheckSameComm(tao, 1, J, 2);
902: }
903: if (Jpre) {
905: PetscCheckSameComm(tao, 1, Jpre, 3);
906: }
907: if (ctx) tao->user_jac_inequalityP = ctx;
908: if (func) tao->ops->computejacobianinequality = func;
909: if (J) {
910: PetscCall(PetscObjectReference((PetscObject)J));
911: PetscCall(MatDestroy(&tao->jacobian_inequality));
912: tao->jacobian_inequality = J;
913: }
914: if (Jpre) {
915: PetscCall(PetscObjectReference((PetscObject)Jpre));
916: PetscCall(MatDestroy(&tao->jacobian_inequality_pre));
917: tao->jacobian_inequality_pre = Jpre;
918: }
919: PetscFunctionReturn(PETSC_SUCCESS);
920: }
922: /*@C
923: TaoGetJacobianInequalityRoutine - Gets the function used to compute inequality constraint Jacobian.
925: Not Collective
927: Input Parameter:
928: . tao - the `Tao` context
930: Output Parameters:
931: + J - the matrix to internally hold the constraint computation
932: . Jpre - the matrix used to construct the preconditioner
933: . func - Jacobian evaluation routine
934: - ctx - the (optional) user-defined context
936: Calling sequence of `func`:
937: + tao - the `Tao` context
938: . x - input vector
939: . J - Jacobian matrix
940: . Jpre - matrix used to construct the preconditioner, usually the same as `J`
941: - ctx - [optional] user-defined Jacobian context
943: Level: intermediate
945: .seealso: [](ch_tao), `Tao`, `TaoComputeJacobianInequality()`, `TaoSetJacobianInequalityRoutine()`
946: @*/
947: PetscErrorCode TaoGetJacobianInequalityRoutine(Tao tao, Mat *J, Mat *Jpre, PetscErrorCode (**func)(Tao tao, Vec x, Mat J, Mat Jpre, PetscCtx ctx), PetscCtxRt ctx)
948: {
949: PetscFunctionBegin;
951: if (J) *J = tao->jacobian_inequality;
952: if (Jpre) *Jpre = tao->jacobian_inequality_pre;
953: if (func) *func = tao->ops->computejacobianinequality;
954: if (ctx) *(void **)ctx = tao->user_jac_inequalityP;
955: PetscFunctionReturn(PETSC_SUCCESS);
956: }