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: PetscErrorCode TaoTestHessian(Tao tao)
122: {
123: Mat A, B, C, D, hessian;
124: Vec x = tao->solution;
125: PetscReal nrm, gnorm;
126: PetscReal threshold = 1.e-5;
127: PetscBool complete_print = PETSC_FALSE, test = PETSC_FALSE, flg;
128: PetscViewer viewer, mviewer;
129: MPI_Comm comm;
130: PetscInt tabs;
131: static PetscBool directionsprinted = PETSC_FALSE;
132: PetscViewerFormat format;
134: PetscFunctionBegin;
135: PetscObjectOptionsBegin((PetscObject)tao);
136: PetscCall(PetscOptionsName("-tao_test_hessian", "Compare hand-coded and finite difference Hessians", "None", &test));
137: PetscCall(PetscOptionsReal("-tao_test_hessian", "Threshold for element difference between hand-coded and finite difference being meaningful", "None", threshold, &threshold, NULL));
138: PetscCall(PetscOptionsViewer("-tao_test_hessian_view", "View difference between hand-coded and finite difference Hessians element entries", "None", &mviewer, &format, &complete_print));
139: PetscOptionsEnd();
140: if (!test) PetscFunctionReturn(PETSC_SUCCESS);
142: PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
143: PetscCall(PetscViewerASCIIGetStdout(comm, &viewer));
144: PetscCall(PetscViewerASCIIGetTab(viewer, &tabs));
145: PetscCall(PetscViewerASCIISetTab(viewer, ((PetscObject)tao)->tablevel));
146: PetscCall(PetscViewerASCIIPrintf(viewer, " ---------- Testing Hessian -------------\n"));
147: if (!complete_print && !directionsprinted) {
148: PetscCall(PetscViewerASCIIPrintf(viewer, " Run with -tao_test_hessian_view and optionally -tao_test_hessian <threshold> to show difference\n"));
149: PetscCall(PetscViewerASCIIPrintf(viewer, " of hand-coded and finite difference Hessian entries greater than <threshold>.\n"));
150: }
151: if (!directionsprinted) {
152: PetscCall(PetscViewerASCIIPrintf(viewer, " Testing hand-coded Hessian, if (for double precision runs) ||J - Jfd||_F/||J||_F is\n"));
153: PetscCall(PetscViewerASCIIPrintf(viewer, " O(1.e-8), the hand-coded Hessian is probably correct.\n"));
154: directionsprinted = PETSC_TRUE;
155: }
156: if (complete_print) PetscCall(PetscViewerPushFormat(mviewer, format));
158: PetscCall(PetscObjectTypeCompare((PetscObject)tao->hessian, MATMFFD, &flg));
159: if (!flg) hessian = tao->hessian;
160: else hessian = tao->hessian_pre;
162: while (hessian) {
163: PetscLayout rmap, cmap;
164: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)hessian, &flg, MATSEQAIJ, MATMPIAIJ, MATSEQDENSE, MATMPIDENSE, MATSEQBAIJ, MATMPIBAIJ, MATSEQSBAIJ, MATMPISBAIJ, ""));
165: if (flg) {
166: A = hessian;
167: PetscCall(PetscObjectReference((PetscObject)A));
168: } else {
169: PetscCall(MatComputeOperator(hessian, MATAIJ, &A));
170: }
172: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
173: PetscCall(MatGetLayouts(A, &rmap, &cmap));
174: PetscCall(MatSetLayouts(B, rmap, cmap));
175: PetscCall(MatSetType(B, ((PetscObject)A)->type_name));
176: PetscCall(MatSetUp(B));
177: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
179: PetscCall(TaoDefaultComputeHessian(tao, x, B, B, NULL));
181: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &D));
182: PetscCall(MatAYPX(D, -1.0, A, DIFFERENT_NONZERO_PATTERN));
183: PetscCall(MatNorm(D, NORM_FROBENIUS, &nrm));
184: PetscCall(MatNorm(A, NORM_FROBENIUS, &gnorm));
185: PetscCall(MatDestroy(&D));
186: if (!gnorm) gnorm = 1; /* just in case */
187: PetscCall(PetscViewerASCIIPrintf(viewer, " ||H - Hfd||_F/||H||_F = %g, ||H - Hfd||_F = %g\n", (double)(nrm / gnorm), (double)nrm));
189: if (complete_print) {
190: PetscCall(PetscViewerASCIIPrintf(viewer, " Hand-coded Hessian ----------\n"));
191: PetscCall(MatView(A, mviewer));
192: PetscCall(PetscViewerASCIIPrintf(viewer, " Finite difference Hessian ----------\n"));
193: PetscCall(MatView(B, mviewer));
194: }
196: if (complete_print) {
197: PetscInt Istart, Iend, *ccols, bncols, cncols, j, row;
198: PetscScalar *cvals;
199: const PetscInt *bcols;
200: const PetscScalar *bvals;
202: PetscCall(MatAYPX(B, -1.0, A, DIFFERENT_NONZERO_PATTERN));
203: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
204: PetscCall(MatSetLayouts(C, rmap, cmap));
205: PetscCall(MatSetType(C, ((PetscObject)A)->type_name));
206: PetscCall(MatSetUp(C));
207: PetscCall(MatSetOption(C, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
208: PetscCall(MatGetOwnershipRange(B, &Istart, &Iend));
210: for (row = Istart; row < Iend; row++) {
211: PetscCall(MatGetRow(B, row, &bncols, &bcols, &bvals));
212: PetscCall(PetscMalloc2(bncols, &ccols, bncols, &cvals));
213: for (j = 0, cncols = 0; j < bncols; j++) {
214: if (PetscAbsScalar(bvals[j]) > threshold) {
215: ccols[cncols] = bcols[j];
216: cvals[cncols] = bvals[j];
217: cncols += 1;
218: }
219: }
220: if (cncols) PetscCall(MatSetValues(C, 1, &row, cncols, ccols, cvals, INSERT_VALUES));
221: PetscCall(MatRestoreRow(B, row, &bncols, &bcols, &bvals));
222: PetscCall(PetscFree2(ccols, cvals));
223: }
224: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
225: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
226: PetscCall(PetscViewerASCIIPrintf(viewer, " Finite-difference minus hand-coded Hessian with tolerance %g ----------\n", (double)threshold));
227: PetscCall(MatView(C, mviewer));
228: PetscCall(MatDestroy(&C));
229: }
230: PetscCall(MatDestroy(&A));
231: PetscCall(MatDestroy(&B));
233: if (hessian != tao->hessian_pre) {
234: hessian = tao->hessian_pre;
235: PetscCall(PetscViewerASCIIPrintf(viewer, " ---------- Testing Hessian for preconditioner -------------\n"));
236: } else hessian = NULL;
237: }
238: if (complete_print) {
239: PetscCall(PetscViewerPopFormat(mviewer));
240: PetscCall(PetscViewerDestroy(&mviewer));
241: }
242: PetscCall(PetscViewerASCIISetTab(viewer, tabs));
243: PetscFunctionReturn(PETSC_SUCCESS);
244: }
246: /*@
247: TaoComputeHessian - Computes the Hessian matrix that has been
248: set with `TaoSetHessian()`.
250: Collective
252: Input Parameters:
253: + tao - the Tao solver context
254: - X - input vector
256: Output Parameters:
257: + H - Hessian matrix
258: - Hpre - matrix used to construct the preconditioner, usually the same as `H`
260: Options Database Keys:
261: + -tao_test_hessian - compare the user provided Hessian with one compute via finite differences to check for errors
262: . -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
263: - -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
265: Level: developer
267: Notes:
268: Most users should not need to explicitly call this routine, as it
269: is used internally within the minimization solvers.
271: `TaoComputeHessian()` is typically used within optimization algorithms,
272: so most users would not generally call this routine
273: themselves.
275: Developer Notes:
276: The Hessian test mechanism follows `SNESTestJacobian()`.
278: If there is no separate preconditioning matrix, `TaoComputeHessian(tao, X, H, NULL)` and `TaoComputeHessian(tao, X, H, H)` are equivalent.
280: .seealso: [](ch_tao), `Tao`, `TaoComputeObjective()`, `TaoComputeObjectiveAndGradient()`, `TaoSetHessian()`
281: @*/
282: PetscErrorCode TaoComputeHessian(Tao tao, Vec X, Mat H, Mat Hpre)
283: {
284: PetscFunctionBegin;
287: PetscCheckSameComm(tao, 1, X, 2);
288: PetscCall(TaoTermMappingComputeHessian(&tao->objective_term, X, tao->objective_parameters, INSERT_VALUES, H, Hpre));
289: PetscCall(TaoTestHessian(tao));
290: PetscFunctionReturn(PETSC_SUCCESS);
291: }
293: /*@
294: TaoComputeJacobian - Computes the Jacobian matrix that has been
295: set with TaoSetJacobianRoutine().
297: Collective
299: Input Parameters:
300: + tao - the Tao solver context
301: - X - input vector
303: Output Parameters:
304: + J - Jacobian matrix
305: - Jpre - matrix used to compute the preconditioner, often the same as `J`
307: Level: developer
309: Notes:
310: Most users should not need to explicitly call this routine, as it
311: is used internally within the minimization solvers.
313: `TaoComputeJacobian()` is typically used within minimization
314: implementations, so most users would not generally call this routine
315: themselves.
317: .seealso: [](ch_tao), `TaoComputeObjective()`, `TaoComputeObjectiveAndGradient()`, `TaoSetJacobianRoutine()`
318: @*/
319: PetscErrorCode TaoComputeJacobian(Tao tao, Vec X, Mat J, Mat Jpre)
320: {
321: PetscFunctionBegin;
324: PetscCheckSameComm(tao, 1, X, 2);
325: ++tao->njac;
326: PetscCall(VecLockReadPush(X));
327: PetscCall(PetscLogEventBegin(TAO_JacobianEval, tao, X, J, Jpre));
328: PetscCallBack("Tao callback Jacobian", (*tao->ops->computejacobian)(tao, X, J, Jpre, tao->user_jacP));
329: PetscCall(PetscLogEventEnd(TAO_JacobianEval, tao, X, J, Jpre));
330: PetscCall(VecLockReadPop(X));
331: PetscFunctionReturn(PETSC_SUCCESS);
332: }
334: /*@
335: TaoComputeResidualJacobian - Computes the least-squares residual Jacobian matrix that has been
336: set with `TaoSetJacobianResidual()`.
338: Collective
340: Input Parameters:
341: + tao - the Tao solver context
342: - X - input vector
344: Output Parameters:
345: + J - Jacobian matrix
346: - Jpre - matrix used to compute the preconditioner, often the same as `J`
348: Level: developer
350: Notes:
351: Most users should not need to explicitly call this routine, as it
352: is used internally within the minimization solvers.
354: `TaoComputeResidualJacobian()` is typically used within least-squares
355: implementations, so most users would not generally call this routine
356: themselves.
358: .seealso: [](ch_tao), `Tao`, `TaoComputeResidual()`, `TaoSetJacobianResidual()`
359: @*/
360: PetscErrorCode TaoComputeResidualJacobian(Tao tao, Vec X, Mat J, Mat Jpre)
361: {
362: PetscFunctionBegin;
365: PetscCheckSameComm(tao, 1, X, 2);
366: ++tao->njac;
367: PetscCall(VecLockReadPush(X));
368: PetscCall(PetscLogEventBegin(TAO_JacobianEval, tao, X, J, Jpre));
369: PetscCallBack("Tao callback least-squares residual Jacobian", (*tao->ops->computeresidualjacobian)(tao, X, J, Jpre, tao->user_lsjacP));
370: PetscCall(PetscLogEventEnd(TAO_JacobianEval, tao, X, J, Jpre));
371: PetscCall(VecLockReadPop(X));
372: PetscFunctionReturn(PETSC_SUCCESS);
373: }
375: /*@
376: TaoComputeJacobianState - Computes the Jacobian matrix that has been
377: set with `TaoSetJacobianStateRoutine()`.
379: Collective
381: Input Parameters:
382: + tao - the `Tao` solver context
383: - X - input vector
385: Output Parameters:
386: + J - Jacobian matrix
387: . Jpre - matrix used to construct the preconditioner, often the same as `J`
388: - Jinv - unknown
390: Level: developer
392: Note:
393: Most users should not need to explicitly call this routine, as it
394: is used internally within the optimization algorithms.
396: .seealso: [](ch_tao), `Tao`, `TaoComputeObjective()`, `TaoComputeObjectiveAndGradient()`, `TaoSetJacobianStateRoutine()`, `TaoComputeJacobianDesign()`, `TaoSetStateDesignIS()`
397: @*/
398: PetscErrorCode TaoComputeJacobianState(Tao tao, Vec X, Mat J, Mat Jpre, Mat Jinv)
399: {
400: PetscFunctionBegin;
403: PetscCheckSameComm(tao, 1, X, 2);
404: ++tao->njac_state;
405: PetscCall(VecLockReadPush(X));
406: PetscCall(PetscLogEventBegin(TAO_JacobianEval, tao, X, J, Jpre));
407: PetscCallBack("Tao callback Jacobian(state)", (*tao->ops->computejacobianstate)(tao, X, J, Jpre, Jinv, tao->user_jac_stateP));
408: PetscCall(PetscLogEventEnd(TAO_JacobianEval, tao, X, J, Jpre));
409: PetscCall(VecLockReadPop(X));
410: PetscFunctionReturn(PETSC_SUCCESS);
411: }
413: /*@
414: TaoComputeJacobianDesign - Computes the Jacobian matrix that has been
415: set with `TaoSetJacobianDesignRoutine()`.
417: Collective
419: Input Parameters:
420: + tao - the Tao solver context
421: - X - input vector
423: Output Parameter:
424: . J - Jacobian matrix
426: Level: developer
428: Note:
429: Most users should not need to explicitly call this routine, as it
430: is used internally within the optimization algorithms.
432: .seealso: [](ch_tao), `Tao`, `TaoComputeObjective()`, `TaoComputeObjectiveAndGradient()`, `TaoSetJacobianDesignRoutine()`, `TaoSetStateDesignIS()`
433: @*/
434: PetscErrorCode TaoComputeJacobianDesign(Tao tao, Vec X, Mat J)
435: {
436: PetscFunctionBegin;
439: PetscCheckSameComm(tao, 1, X, 2);
440: ++tao->njac_design;
441: PetscCall(VecLockReadPush(X));
442: PetscCall(PetscLogEventBegin(TAO_JacobianEval, tao, X, J, NULL));
443: PetscCallBack("Tao callback Jacobian(design)", (*tao->ops->computejacobiandesign)(tao, X, J, tao->user_jac_designP));
444: PetscCall(PetscLogEventEnd(TAO_JacobianEval, tao, X, J, NULL));
445: PetscCall(VecLockReadPop(X));
446: PetscFunctionReturn(PETSC_SUCCESS);
447: }
449: /*@C
450: TaoSetJacobianRoutine - Sets the function to compute the Jacobian as well as the location to store the matrix.
452: Logically Collective
454: Input Parameters:
455: + tao - the `Tao` context
456: . J - Matrix used for the Jacobian
457: . Jpre - Matrix that will be used to construct the preconditioner, can be same as `J`
458: . func - Jacobian evaluation routine
459: - ctx - [optional] user-defined context for private data for the
460: Jacobian evaluation routine (may be `NULL`)
462: Calling sequence of `func`:
463: + tao - the `Tao` context
464: . x - input vector
465: . J - Jacobian matrix
466: . Jpre - matrix used to construct the preconditioner, usually the same as `J`
467: - ctx - [optional] user-defined Jacobian context
469: Level: intermediate
471: .seealso: [](ch_tao), `Tao`, `TaoSetGradient()`, `TaoSetObjective()`
472: @*/
473: PetscErrorCode TaoSetJacobianRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao tao, Vec x, Mat J, Mat Jpre, PetscCtx ctx), PetscCtx ctx)
474: {
475: PetscFunctionBegin;
477: if (J) {
479: PetscCheckSameComm(tao, 1, J, 2);
480: }
481: if (Jpre) {
483: PetscCheckSameComm(tao, 1, Jpre, 3);
484: }
485: if (ctx) tao->user_jacP = ctx;
486: if (func) tao->ops->computejacobian = func;
487: if (J) {
488: PetscCall(PetscObjectReference((PetscObject)J));
489: PetscCall(MatDestroy(&tao->jacobian));
490: tao->jacobian = J;
491: }
492: if (Jpre) {
493: PetscCall(PetscObjectReference((PetscObject)Jpre));
494: PetscCall(MatDestroy(&tao->jacobian_pre));
495: tao->jacobian_pre = Jpre;
496: }
497: PetscFunctionReturn(PETSC_SUCCESS);
498: }
500: /*@C
501: TaoSetJacobianResidualRoutine - Sets the function to compute the least-squares residual Jacobian as well as the
502: location to store the matrix.
504: Logically Collective
506: Input Parameters:
507: + tao - the `Tao` context
508: . J - Matrix used for the jacobian
509: . Jpre - Matrix that will be used to construct the preconditioner, can be same as `J`
510: . func - Jacobian evaluation routine
511: - ctx - [optional] user-defined context for private data for the
512: Jacobian evaluation routine (may be `NULL`)
514: Calling sequence of `func`:
515: + tao - the `Tao` context
516: . x - input vector
517: . J - Jacobian matrix
518: . Jpre - matrix used to construct the preconditioner, usually the same as `J`
519: - ctx - [optional] user-defined Jacobian context
521: Level: intermediate
523: .seealso: [](ch_tao), `Tao`, `TaoSetGradient()`, `TaoSetObjective()`
524: @*/
525: PetscErrorCode TaoSetJacobianResidualRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao tao, Vec x, Mat J, Mat Jpre, PetscCtx ctx), PetscCtx ctx)
526: {
527: PetscFunctionBegin;
529: if (J) {
531: PetscCheckSameComm(tao, 1, J, 2);
532: }
533: if (Jpre) {
535: PetscCheckSameComm(tao, 1, Jpre, 3);
536: }
537: if (ctx) tao->user_lsjacP = ctx;
538: if (func) tao->ops->computeresidualjacobian = func;
539: if (J) {
540: PetscCall(PetscObjectReference((PetscObject)J));
541: PetscCall(MatDestroy(&tao->ls_jac));
542: tao->ls_jac = J;
543: }
544: if (Jpre) {
545: PetscCall(PetscObjectReference((PetscObject)Jpre));
546: PetscCall(MatDestroy(&tao->ls_jac_pre));
547: tao->ls_jac_pre = Jpre;
548: }
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }
552: /*@C
553: TaoSetJacobianStateRoutine - Sets the function to compute the Jacobian
554: (and its inverse) of the constraint function with respect to the state variables.
555: Used only for PDE-constrained optimization.
557: Logically Collective
559: Input Parameters:
560: + tao - the `Tao` context
561: . J - Matrix used for the Jacobian
562: . Jpre - Matrix that will be used to construct the preconditioner, can be same as `J`. Only used if `Jinv` is `NULL`
563: . Jinv - [optional] Matrix used to apply the inverse of the state Jacobian. Use `NULL` to default to PETSc `KSP` solvers to apply the inverse.
564: . func - Jacobian evaluation routine
565: - ctx - [optional] user-defined context for private data for the
566: Jacobian evaluation routine (may be `NULL`)
568: Calling sequence of `func`:
569: + tao - the `Tao` context
570: . x - input vector
571: . J - Jacobian matrix
572: . Jpre - matrix used to construct the preconditioner, usually the same as `J`
573: . Jinv - inverse of `J`
574: - ctx - [optional] user-defined Jacobian context
576: Level: intermediate
578: .seealso: [](ch_tao), `Tao`, `TaoComputeJacobianState()`, `TaoSetJacobianDesignRoutine()`, `TaoSetStateDesignIS()`
579: @*/
580: 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)
581: {
582: PetscFunctionBegin;
584: if (J) {
586: PetscCheckSameComm(tao, 1, J, 2);
587: }
588: if (Jpre) {
590: PetscCheckSameComm(tao, 1, Jpre, 3);
591: }
592: if (Jinv) {
594: PetscCheckSameComm(tao, 1, Jinv, 4);
595: }
596: if (ctx) tao->user_jac_stateP = ctx;
597: if (func) tao->ops->computejacobianstate = func;
598: if (J) {
599: PetscCall(PetscObjectReference((PetscObject)J));
600: PetscCall(MatDestroy(&tao->jacobian_state));
601: tao->jacobian_state = J;
602: }
603: if (Jpre) {
604: PetscCall(PetscObjectReference((PetscObject)Jpre));
605: PetscCall(MatDestroy(&tao->jacobian_state_pre));
606: tao->jacobian_state_pre = Jpre;
607: }
608: if (Jinv) {
609: PetscCall(PetscObjectReference((PetscObject)Jinv));
610: PetscCall(MatDestroy(&tao->jacobian_state_inv));
611: tao->jacobian_state_inv = Jinv;
612: }
613: PetscFunctionReturn(PETSC_SUCCESS);
614: }
616: /*@C
617: TaoSetJacobianDesignRoutine - Sets the function to compute the Jacobian of
618: the constraint function with respect to the design variables. Used only for
619: PDE-constrained optimization.
621: Logically Collective
623: Input Parameters:
624: + tao - the `Tao` context
625: . J - Matrix used for the Jacobian
626: . func - Jacobian evaluation routine
627: - ctx - [optional] user-defined context for private data for the
628: Jacobian evaluation routine (may be `NULL`)
630: Calling sequence of `func`:
631: + tao - the `Tao` context
632: . x - input vector
633: . J - Jacobian matrix
634: - ctx - [optional] user-defined Jacobian context
636: Level: intermediate
638: .seealso: [](ch_tao), `Tao`, `TaoComputeJacobianDesign()`, `TaoSetJacobianStateRoutine()`, `TaoSetStateDesignIS()`
639: @*/
640: PetscErrorCode TaoSetJacobianDesignRoutine(Tao tao, Mat J, PetscErrorCode (*func)(Tao tao, Vec x, Mat J, PetscCtx ctx), PetscCtx ctx)
641: {
642: PetscFunctionBegin;
644: if (J) {
646: PetscCheckSameComm(tao, 1, J, 2);
647: }
648: if (ctx) tao->user_jac_designP = ctx;
649: if (func) tao->ops->computejacobiandesign = func;
650: if (J) {
651: PetscCall(PetscObjectReference((PetscObject)J));
652: PetscCall(MatDestroy(&tao->jacobian_design));
653: tao->jacobian_design = J;
654: }
655: PetscFunctionReturn(PETSC_SUCCESS);
656: }
658: /*@
659: TaoSetStateDesignIS - Indicate to the `Tao` object which variables in the
660: solution vector are state variables and which are design. Only applies to
661: PDE-constrained optimization.
663: Logically Collective
665: Input Parameters:
666: + tao - The `Tao` context
667: . s_is - the index set corresponding to the state variables
668: - d_is - the index set corresponding to the design variables
670: Level: intermediate
672: .seealso: [](ch_tao), `Tao`, `TaoSetJacobianStateRoutine()`, `TaoSetJacobianDesignRoutine()`
673: @*/
674: PetscErrorCode TaoSetStateDesignIS(Tao tao, IS s_is, IS d_is)
675: {
676: PetscFunctionBegin;
677: PetscCall(PetscObjectReference((PetscObject)s_is));
678: PetscCall(ISDestroy(&tao->state_is));
679: tao->state_is = s_is;
680: PetscCall(PetscObjectReference((PetscObject)d_is));
681: PetscCall(ISDestroy(&tao->design_is));
682: tao->design_is = d_is;
683: PetscFunctionReturn(PETSC_SUCCESS);
684: }
686: /*@
687: TaoComputeJacobianEquality - Computes the Jacobian matrix that has been
688: set with `TaoSetJacobianEqualityRoutine()`.
690: Collective
692: Input Parameters:
693: + tao - the `Tao` solver context
694: - X - input vector
696: Output Parameters:
697: + J - Jacobian matrix
698: - Jpre - matrix used to construct the preconditioner, often the same as `J`
700: Level: developer
702: Notes:
703: Most users should not need to explicitly call this routine, as it
704: is used internally within the optimization algorithms.
706: .seealso: [](ch_tao), `TaoComputeObjective()`, `TaoComputeObjectiveAndGradient()`, `TaoSetJacobianStateRoutine()`, `TaoComputeJacobianDesign()`, `TaoSetStateDesignIS()`
707: @*/
708: PetscErrorCode TaoComputeJacobianEquality(Tao tao, Vec X, Mat J, Mat Jpre)
709: {
710: PetscFunctionBegin;
713: PetscCheckSameComm(tao, 1, X, 2);
714: ++tao->njac_equality;
715: PetscCall(VecLockReadPush(X));
716: PetscCall(PetscLogEventBegin(TAO_JacobianEval, tao, X, J, Jpre));
717: PetscCallBack("Tao callback Jacobian(equality)", (*tao->ops->computejacobianequality)(tao, X, J, Jpre, tao->user_jac_equalityP));
718: PetscCall(PetscLogEventEnd(TAO_JacobianEval, tao, X, J, Jpre));
719: PetscCall(VecLockReadPop(X));
720: PetscFunctionReturn(PETSC_SUCCESS);
721: }
723: /*@
724: TaoComputeJacobianInequality - Computes the Jacobian matrix that has been
725: set with `TaoSetJacobianInequalityRoutine()`.
727: Collective
729: Input Parameters:
730: + tao - the `Tao` solver context
731: - X - input vector
733: Output Parameters:
734: + J - Jacobian matrix
735: - Jpre - matrix used to construct the preconditioner
737: Level: developer
739: Note:
740: Most users should not need to explicitly call this routine, as it
741: is used internally within the minimization solvers.
743: .seealso: [](ch_tao), `Tao`, `TaoComputeObjective()`, `TaoComputeObjectiveAndGradient()`, `TaoSetJacobianStateRoutine()`, `TaoComputeJacobianDesign()`, `TaoSetStateDesignIS()`
744: @*/
745: PetscErrorCode TaoComputeJacobianInequality(Tao tao, Vec X, Mat J, Mat Jpre)
746: {
747: PetscFunctionBegin;
750: PetscCheckSameComm(tao, 1, X, 2);
751: ++tao->njac_inequality;
752: PetscCall(VecLockReadPush(X));
753: PetscCall(PetscLogEventBegin(TAO_JacobianEval, tao, X, J, Jpre));
754: PetscCallBack("Tao callback Jacobian (inequality)", (*tao->ops->computejacobianinequality)(tao, X, J, Jpre, tao->user_jac_inequalityP));
755: PetscCall(PetscLogEventEnd(TAO_JacobianEval, tao, X, J, Jpre));
756: PetscCall(VecLockReadPop(X));
757: PetscFunctionReturn(PETSC_SUCCESS);
758: }
760: /*@C
761: TaoSetJacobianEqualityRoutine - Sets the function to compute the Jacobian
762: (and its inverse) of the constraint function with respect to the equality variables.
763: Used only for PDE-constrained optimization.
765: Logically Collective
767: Input Parameters:
768: + tao - the `Tao` context
769: . J - Matrix used for the Jacobian
770: . Jpre - Matrix that will be used to construct the preconditioner, can be same as `J`.
771: . func - Jacobian evaluation routine
772: - ctx - [optional] user-defined context for private data for the
773: Jacobian evaluation routine (may be `NULL`)
775: Calling sequence of `func`:
776: + tao - the `Tao` context
777: . x - input vector
778: . J - Jacobian matrix
779: . Jpre - matrix used to construct the preconditioner, usually the same as `J`
780: - ctx - [optional] user-defined Jacobian context
782: Level: intermediate
784: .seealso: [](ch_tao), `Tao`, `TaoComputeJacobianEquality()`, `TaoSetJacobianDesignRoutine()`, `TaoSetEqualityDesignIS()`
785: @*/
786: PetscErrorCode TaoSetJacobianEqualityRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao tao, Vec x, Mat J, Mat Jpre, PetscCtx ctx), PetscCtx ctx)
787: {
788: PetscFunctionBegin;
790: if (J) {
792: PetscCheckSameComm(tao, 1, J, 2);
793: }
794: if (Jpre) {
796: PetscCheckSameComm(tao, 1, Jpre, 3);
797: }
798: if (ctx) tao->user_jac_equalityP = ctx;
799: if (func) tao->ops->computejacobianequality = func;
800: if (J) {
801: PetscCall(PetscObjectReference((PetscObject)J));
802: PetscCall(MatDestroy(&tao->jacobian_equality));
803: tao->jacobian_equality = J;
804: }
805: if (Jpre) {
806: PetscCall(PetscObjectReference((PetscObject)Jpre));
807: PetscCall(MatDestroy(&tao->jacobian_equality_pre));
808: tao->jacobian_equality_pre = Jpre;
809: }
810: PetscFunctionReturn(PETSC_SUCCESS);
811: }
813: /*@C
814: TaoGetJacobianEqualityRoutine - Gets the function used to compute equality constraint Jacobian.
816: Not Collective
818: Input Parameter:
819: . tao - the `Tao` context
821: Output Parameters:
822: + J - the matrix to internally hold the constraint computation
823: . Jpre - the matrix used to construct the preconditioner
824: . func - Jacobian evaluation routine
825: - ctx - the (optional) user-defined context
827: Calling sequence of `func`:
828: + tao - the `Tao` context
829: . x - input vector
830: . J - Jacobian matrix
831: . Jpre - matrix used to construct the preconditioner, usually the same as `J`
832: - ctx - [optional] user-defined Jacobian context
834: Level: intermediate
836: .seealso: [](ch_tao), `Tao`, `TaoComputeJacobianEquality()`, `TaoSetJacobianEqualityRoutine()`
837: @*/
838: PetscErrorCode TaoGetJacobianEqualityRoutine(Tao tao, Mat *J, Mat *Jpre, PetscErrorCode (**func)(Tao tao, Vec x, Mat J, Mat Jpre, PetscCtx ctx), PetscCtxRt ctx)
839: {
840: PetscFunctionBegin;
842: if (J) *J = tao->jacobian_equality;
843: if (Jpre) *Jpre = tao->jacobian_equality_pre;
844: if (func) *func = tao->ops->computejacobianequality;
845: if (ctx) *(void **)ctx = tao->user_jac_equalityP;
846: PetscFunctionReturn(PETSC_SUCCESS);
847: }
849: /*@C
850: TaoSetJacobianInequalityRoutine - Sets the function to compute the Jacobian
851: (and its inverse) of the constraint function with respect to the inequality variables.
852: Used only for PDE-constrained optimization.
854: Logically Collective
856: Input Parameters:
857: + tao - the `Tao` context
858: . J - Matrix used for the Jacobian
859: . Jpre - Matrix that will be used to construct the preconditioner, can be same as `J`.
860: . func - Jacobian evaluation routine
861: - ctx - [optional] user-defined context for private data for the
862: Jacobian evaluation routine (may be `NULL`)
864: Calling sequence of `func`:
865: + tao - the `Tao` context
866: . x - input vector
867: . J - Jacobian matrix
868: . Jpre - matrix used to construct the preconditioner, usually the same as `J`
869: - ctx - [optional] user-defined Jacobian context
871: Level: intermediate
873: .seealso: [](ch_tao), `Tao`, `TaoComputeJacobianInequality()`, `TaoSetJacobianDesignRoutine()`, `TaoSetInequalityDesignIS()`
874: @*/
875: PetscErrorCode TaoSetJacobianInequalityRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao tao, Vec x, Mat J, Mat Jpre, PetscCtx ctx), PetscCtx ctx)
876: {
877: PetscFunctionBegin;
879: if (J) {
881: PetscCheckSameComm(tao, 1, J, 2);
882: }
883: if (Jpre) {
885: PetscCheckSameComm(tao, 1, Jpre, 3);
886: }
887: if (ctx) tao->user_jac_inequalityP = ctx;
888: if (func) tao->ops->computejacobianinequality = func;
889: if (J) {
890: PetscCall(PetscObjectReference((PetscObject)J));
891: PetscCall(MatDestroy(&tao->jacobian_inequality));
892: tao->jacobian_inequality = J;
893: }
894: if (Jpre) {
895: PetscCall(PetscObjectReference((PetscObject)Jpre));
896: PetscCall(MatDestroy(&tao->jacobian_inequality_pre));
897: tao->jacobian_inequality_pre = Jpre;
898: }
899: PetscFunctionReturn(PETSC_SUCCESS);
900: }
902: /*@C
903: TaoGetJacobianInequalityRoutine - Gets the function used to compute inequality constraint Jacobian.
905: Not Collective
907: Input Parameter:
908: . tao - the `Tao` context
910: Output Parameters:
911: + J - the matrix to internally hold the constraint computation
912: . Jpre - the matrix used to construct the preconditioner
913: . func - Jacobian evaluation routine
914: - ctx - the (optional) user-defined context
916: Calling sequence of `func`:
917: + tao - the `Tao` context
918: . x - input vector
919: . J - Jacobian matrix
920: . Jpre - matrix used to construct the preconditioner, usually the same as `J`
921: - ctx - [optional] user-defined Jacobian context
923: Level: intermediate
925: .seealso: [](ch_tao), `Tao`, `TaoComputeJacobianInequality()`, `TaoSetJacobianInequalityRoutine()`
926: @*/
927: PetscErrorCode TaoGetJacobianInequalityRoutine(Tao tao, Mat *J, Mat *Jpre, PetscErrorCode (**func)(Tao tao, Vec x, Mat J, Mat Jpre, PetscCtx ctx), PetscCtxRt ctx)
928: {
929: PetscFunctionBegin;
931: if (J) *J = tao->jacobian_inequality;
932: if (Jpre) *Jpre = tao->jacobian_inequality_pre;
933: if (func) *func = tao->ops->computejacobianinequality;
934: if (ctx) *(void **)ctx = tao->user_jac_inequalityP;
935: PetscFunctionReturn(PETSC_SUCCESS);
936: }