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`, `TaoTypes`, `TaoSetObjective()`, `TaoSetGradient()`, `TaoSetObjectiveAndGradient()`, `TaoGetHessian()`
26: @*/
27: PetscErrorCode TaoSetHessian(Tao tao, Mat H, Mat Hpre, PetscErrorCode (*func)(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx), void *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: if (ctx) tao->user_hessP = ctx;
40: if (func) tao->ops->computehessian = func;
41: if (H) {
42: PetscCall(PetscObjectReference((PetscObject)H));
43: PetscCall(MatDestroy(&tao->hessian));
44: tao->hessian = H;
45: }
46: if (Hpre) {
47: PetscCall(PetscObjectReference((PetscObject)Hpre));
48: PetscCall(MatDestroy(&tao->hessian_pre));
49: tao->hessian_pre = Hpre;
50: }
51: PetscFunctionReturn(PETSC_SUCCESS);
52: }
54: /*@C
55: TaoGetHessian - Gets the function to compute the Hessian as well as the location to store the matrix.
57: Not Collective
59: Input Parameter:
60: . tao - the `Tao` context
62: Output Parameters:
63: + H - Matrix used for the hessian
64: . Hpre - Matrix that will be used to construct the preconditioner, can be the same as `H`
65: . func - Hessian evaluation routine
66: - ctx - user-defined context for private data for the Hessian evaluation routine
68: Calling sequence of `func`:
69: + tao - the `Tao` context
70: . x - input vector
71: . H - Hessian matrix
72: . Hpre - matrix used to construct the preconditioner, usually the same as `H`
73: - ctx - [optional] user-defined Hessian context
75: Level: beginner
77: .seealso: [](ch_tao), `Tao`, `TaoType`, `TaoGetObjective()`, `TaoGetGradient()`, `TaoGetObjectiveAndGradient()`, `TaoSetHessian()`
78: @*/
79: PetscErrorCode TaoGetHessian(Tao tao, Mat *H, Mat *Hpre, PetscErrorCode (**func)(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx), void **ctx)
80: {
81: PetscFunctionBegin;
83: if (H) *H = tao->hessian;
84: if (Hpre) *Hpre = tao->hessian_pre;
85: if (ctx) *ctx = tao->user_hessP;
86: if (func) *func = tao->ops->computehessian;
87: PetscFunctionReturn(PETSC_SUCCESS);
88: }
90: PetscErrorCode TaoTestHessian(Tao tao)
91: {
92: Mat A, B, C, D, hessian;
93: Vec x = tao->solution;
94: PetscReal nrm, gnorm;
95: PetscReal threshold = 1.e-5;
96: PetscInt m, n, M, N;
97: PetscBool complete_print = PETSC_FALSE, test = PETSC_FALSE, flg;
98: PetscViewer viewer, mviewer;
99: MPI_Comm comm;
100: PetscInt tabs;
101: static PetscBool directionsprinted = PETSC_FALSE;
102: PetscViewerFormat format;
104: PetscFunctionBegin;
105: PetscObjectOptionsBegin((PetscObject)tao);
106: PetscCall(PetscOptionsName("-tao_test_hessian", "Compare hand-coded and finite difference Hessians", "None", &test));
107: PetscCall(PetscOptionsReal("-tao_test_hessian", "Threshold for element difference between hand-coded and finite difference being meaningful", "None", threshold, &threshold, NULL));
108: PetscCall(PetscOptionsViewer("-tao_test_hessian_view", "View difference between hand-coded and finite difference Hessians element entries", "None", &mviewer, &format, &complete_print));
109: PetscOptionsEnd();
110: if (!test) PetscFunctionReturn(PETSC_SUCCESS);
112: PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
113: PetscCall(PetscViewerASCIIGetStdout(comm, &viewer));
114: PetscCall(PetscViewerASCIIGetTab(viewer, &tabs));
115: PetscCall(PetscViewerASCIISetTab(viewer, ((PetscObject)tao)->tablevel));
116: PetscCall(PetscViewerASCIIPrintf(viewer, " ---------- Testing Hessian -------------\n"));
117: if (!complete_print && !directionsprinted) {
118: PetscCall(PetscViewerASCIIPrintf(viewer, " Run with -tao_test_hessian_view and optionally -tao_test_hessian <threshold> to show difference\n"));
119: PetscCall(PetscViewerASCIIPrintf(viewer, " of hand-coded and finite difference Hessian entries greater than <threshold>.\n"));
120: }
121: if (!directionsprinted) {
122: PetscCall(PetscViewerASCIIPrintf(viewer, " Testing hand-coded Hessian, if (for double precision runs) ||J - Jfd||_F/||J||_F is\n"));
123: PetscCall(PetscViewerASCIIPrintf(viewer, " O(1.e-8), the hand-coded Hessian is probably correct.\n"));
124: directionsprinted = PETSC_TRUE;
125: }
126: if (complete_print) PetscCall(PetscViewerPushFormat(mviewer, format));
128: PetscCall(PetscObjectTypeCompare((PetscObject)tao->hessian, MATMFFD, &flg));
129: if (!flg) hessian = tao->hessian;
130: else hessian = tao->hessian_pre;
132: while (hessian) {
133: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)hessian, &flg, MATSEQAIJ, MATMPIAIJ, MATSEQDENSE, MATMPIDENSE, MATSEQBAIJ, MATMPIBAIJ, MATSEQSBAIJ, MATMPIBAIJ, ""));
134: if (flg) {
135: A = hessian;
136: PetscCall(PetscObjectReference((PetscObject)A));
137: } else {
138: PetscCall(MatComputeOperator(hessian, MATAIJ, &A));
139: }
141: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
142: PetscCall(MatGetSize(A, &M, &N));
143: PetscCall(MatGetLocalSize(A, &m, &n));
144: PetscCall(MatSetSizes(B, m, n, M, N));
145: PetscCall(MatSetType(B, ((PetscObject)A)->type_name));
146: PetscCall(MatSetUp(B));
147: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
149: PetscCall(TaoDefaultComputeHessian(tao, x, B, B, NULL));
151: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &D));
152: PetscCall(MatAYPX(D, -1.0, A, DIFFERENT_NONZERO_PATTERN));
153: PetscCall(MatNorm(D, NORM_FROBENIUS, &nrm));
154: PetscCall(MatNorm(A, NORM_FROBENIUS, &gnorm));
155: PetscCall(MatDestroy(&D));
156: if (!gnorm) gnorm = 1; /* just in case */
157: PetscCall(PetscViewerASCIIPrintf(viewer, " ||H - Hfd||_F/||H||_F = %g, ||H - Hfd||_F = %g\n", (double)(nrm / gnorm), (double)nrm));
159: if (complete_print) {
160: PetscCall(PetscViewerASCIIPrintf(viewer, " Hand-coded Hessian ----------\n"));
161: PetscCall(MatView(A, mviewer));
162: PetscCall(PetscViewerASCIIPrintf(viewer, " Finite difference Hessian ----------\n"));
163: PetscCall(MatView(B, mviewer));
164: }
166: if (complete_print) {
167: PetscInt Istart, Iend, *ccols, bncols, cncols, j, row;
168: PetscScalar *cvals;
169: const PetscInt *bcols;
170: const PetscScalar *bvals;
172: PetscCall(MatAYPX(B, -1.0, A, DIFFERENT_NONZERO_PATTERN));
173: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
174: PetscCall(MatSetSizes(C, m, n, M, N));
175: PetscCall(MatSetType(C, ((PetscObject)A)->type_name));
176: PetscCall(MatSetUp(C));
177: PetscCall(MatSetOption(C, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
178: PetscCall(MatGetOwnershipRange(B, &Istart, &Iend));
180: for (row = Istart; row < Iend; row++) {
181: PetscCall(MatGetRow(B, row, &bncols, &bcols, &bvals));
182: PetscCall(PetscMalloc2(bncols, &ccols, bncols, &cvals));
183: for (j = 0, cncols = 0; j < bncols; j++) {
184: if (PetscAbsScalar(bvals[j]) > threshold) {
185: ccols[cncols] = bcols[j];
186: cvals[cncols] = bvals[j];
187: cncols += 1;
188: }
189: }
190: if (cncols) PetscCall(MatSetValues(C, 1, &row, cncols, ccols, cvals, INSERT_VALUES));
191: PetscCall(MatRestoreRow(B, row, &bncols, &bcols, &bvals));
192: PetscCall(PetscFree2(ccols, cvals));
193: }
194: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
195: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
196: PetscCall(PetscViewerASCIIPrintf(viewer, " Finite-difference minus hand-coded Hessian with tolerance %g ----------\n", (double)threshold));
197: PetscCall(MatView(C, mviewer));
198: PetscCall(MatDestroy(&C));
199: }
200: PetscCall(MatDestroy(&A));
201: PetscCall(MatDestroy(&B));
203: if (hessian != tao->hessian_pre) {
204: hessian = tao->hessian_pre;
205: PetscCall(PetscViewerASCIIPrintf(viewer, " ---------- Testing Hessian for preconditioner -------------\n"));
206: } else hessian = NULL;
207: }
208: if (complete_print) {
209: PetscCall(PetscViewerPopFormat(mviewer));
210: PetscCall(PetscViewerDestroy(&mviewer));
211: }
212: PetscCall(PetscViewerASCIISetTab(viewer, tabs));
213: PetscFunctionReturn(PETSC_SUCCESS);
214: }
216: /*@
217: TaoComputeHessian - Computes the Hessian matrix that has been
218: set with `TaoSetHessian()`.
220: Collective
222: Input Parameters:
223: + tao - the Tao solver context
224: - X - input vector
226: Output Parameters:
227: + H - Hessian matrix
228: - Hpre - Preconditioning matrix
230: Options Database Keys:
231: + -tao_test_hessian - compare the user provided Hessian with one compute via finite differences to check for errors
232: . -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
233: - -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
235: Level: developer
237: Notes:
238: Most users should not need to explicitly call this routine, as it
239: is used internally within the minimization solvers.
241: `TaoComputeHessian()` is typically used within optimization algorithms,
242: so most users would not generally call this routine
243: themselves.
245: Developer Notes:
246: The Hessian test mechanism follows `SNESTestJacobian()`.
248: .seealso: [](ch_tao), `Tao`, `TaoComputeObjective()`, `TaoComputeObjectiveAndGradient()`, `TaoSetHessian()`
249: @*/
250: PetscErrorCode TaoComputeHessian(Tao tao, Vec X, Mat H, Mat Hpre)
251: {
252: PetscFunctionBegin;
255: PetscCheckSameComm(tao, 1, X, 2);
256: PetscCheck(tao->ops->computehessian, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "TaoSetHessian() not called");
258: ++tao->nhess;
259: PetscCall(VecLockReadPush(X));
260: PetscCall(PetscLogEventBegin(TAO_HessianEval, tao, X, H, Hpre));
261: PetscCallBack("Tao callback Hessian", (*tao->ops->computehessian)(tao, X, H, Hpre, tao->user_hessP));
262: PetscCall(PetscLogEventEnd(TAO_HessianEval, tao, X, H, Hpre));
263: PetscCall(VecLockReadPop(X));
265: PetscCall(TaoTestHessian(tao));
266: PetscFunctionReturn(PETSC_SUCCESS);
267: }
269: /*@
270: TaoComputeJacobian - Computes the Jacobian matrix that has been
271: set with TaoSetJacobianRoutine().
273: Collective
275: Input Parameters:
276: + tao - the Tao solver context
277: - X - input vector
279: Output Parameters:
280: + J - Jacobian matrix
281: - Jpre - Preconditioning matrix
283: Level: developer
285: Notes:
286: Most users should not need to explicitly call this routine, as it
287: is used internally within the minimization solvers.
289: `TaoComputeJacobian()` is typically used within minimization
290: implementations, so most users would not generally call this routine
291: themselves.
293: .seealso: [](ch_tao), `TaoComputeObjective()`, `TaoComputeObjectiveAndGradient()`, `TaoSetJacobianRoutine()`
294: @*/
295: PetscErrorCode TaoComputeJacobian(Tao tao, Vec X, Mat J, Mat Jpre)
296: {
297: PetscFunctionBegin;
300: PetscCheckSameComm(tao, 1, X, 2);
301: ++tao->njac;
302: PetscCall(VecLockReadPush(X));
303: PetscCall(PetscLogEventBegin(TAO_JacobianEval, tao, X, J, Jpre));
304: PetscCallBack("Tao callback Jacobian", (*tao->ops->computejacobian)(tao, X, J, Jpre, tao->user_jacP));
305: PetscCall(PetscLogEventEnd(TAO_JacobianEval, tao, X, J, Jpre));
306: PetscCall(VecLockReadPop(X));
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: /*@
311: TaoComputeResidualJacobian - Computes the least-squares residual Jacobian matrix that has been
312: set with `TaoSetJacobianResidual()`.
314: Collective
316: Input Parameters:
317: + tao - the Tao solver context
318: - X - input vector
320: Output Parameters:
321: + J - Jacobian matrix
322: - Jpre - Preconditioning matrix
324: Level: developer
326: Notes:
327: Most users should not need to explicitly call this routine, as it
328: is used internally within the minimization solvers.
330: `TaoComputeResidualJacobian()` is typically used within least-squares
331: implementations, so most users would not generally call this routine
332: themselves.
334: .seealso: [](ch_tao), `Tao`, `TaoComputeResidual()`, `TaoSetJacobianResidual()`
335: @*/
336: PetscErrorCode TaoComputeResidualJacobian(Tao tao, Vec X, Mat J, Mat Jpre)
337: {
338: PetscFunctionBegin;
341: PetscCheckSameComm(tao, 1, X, 2);
342: ++tao->njac;
343: PetscCall(VecLockReadPush(X));
344: PetscCall(PetscLogEventBegin(TAO_JacobianEval, tao, X, J, Jpre));
345: PetscCallBack("Tao callback least-squares residual Jacobian", (*tao->ops->computeresidualjacobian)(tao, X, J, Jpre, tao->user_lsjacP));
346: PetscCall(PetscLogEventEnd(TAO_JacobianEval, tao, X, J, Jpre));
347: PetscCall(VecLockReadPop(X));
348: PetscFunctionReturn(PETSC_SUCCESS);
349: }
351: /*@
352: TaoComputeJacobianState - Computes the Jacobian matrix that has been
353: set with `TaoSetJacobianStateRoutine()`.
355: Collective
357: Input Parameters:
358: + tao - the `Tao` solver context
359: - X - input vector
361: Output Parameters:
362: + J - Jacobian matrix
363: . Jpre - matrix used to construct the preconditioner, often the same as `J`
364: - Jinv - unknown
366: Level: developer
368: Note:
369: Most users should not need to explicitly call this routine, as it
370: is used internally within the optimization algorithms.
372: .seealso: [](ch_tao), `Tao`, `TaoComputeObjective()`, `TaoComputeObjectiveAndGradient()`, `TaoSetJacobianStateRoutine()`, `TaoComputeJacobianDesign()`, `TaoSetStateDesignIS()`
373: @*/
374: PetscErrorCode TaoComputeJacobianState(Tao tao, Vec X, Mat J, Mat Jpre, Mat Jinv)
375: {
376: PetscFunctionBegin;
379: PetscCheckSameComm(tao, 1, X, 2);
380: ++tao->njac_state;
381: PetscCall(VecLockReadPush(X));
382: PetscCall(PetscLogEventBegin(TAO_JacobianEval, tao, X, J, Jpre));
383: PetscCallBack("Tao callback Jacobian(state)", (*tao->ops->computejacobianstate)(tao, X, J, Jpre, Jinv, tao->user_jac_stateP));
384: PetscCall(PetscLogEventEnd(TAO_JacobianEval, tao, X, J, Jpre));
385: PetscCall(VecLockReadPop(X));
386: PetscFunctionReturn(PETSC_SUCCESS);
387: }
389: /*@
390: TaoComputeJacobianDesign - Computes the Jacobian matrix that has been
391: set with `TaoSetJacobianDesignRoutine()`.
393: Collective
395: Input Parameters:
396: + tao - the Tao solver context
397: - X - input vector
399: Output Parameter:
400: . J - Jacobian matrix
402: Level: developer
404: Note:
405: Most users should not need to explicitly call this routine, as it
406: is used internally within the optimization algorithms.
408: .seealso: [](ch_tao), `Tao`, `TaoComputeObjective()`, `TaoComputeObjectiveAndGradient()`, `TaoSetJacobianDesignRoutine()`, `TaoSetStateDesignIS()`
409: @*/
410: PetscErrorCode TaoComputeJacobianDesign(Tao tao, Vec X, Mat J)
411: {
412: PetscFunctionBegin;
415: PetscCheckSameComm(tao, 1, X, 2);
416: ++tao->njac_design;
417: PetscCall(VecLockReadPush(X));
418: PetscCall(PetscLogEventBegin(TAO_JacobianEval, tao, X, J, NULL));
419: PetscCallBack("Tao callback Jacobian(design)", (*tao->ops->computejacobiandesign)(tao, X, J, tao->user_jac_designP));
420: PetscCall(PetscLogEventEnd(TAO_JacobianEval, tao, X, J, NULL));
421: PetscCall(VecLockReadPop(X));
422: PetscFunctionReturn(PETSC_SUCCESS);
423: }
425: /*@C
426: TaoSetJacobianRoutine - Sets the function to compute the Jacobian as well as the location to store the matrix.
428: Logically Collective
430: Input Parameters:
431: + tao - the `Tao` context
432: . J - Matrix used for the Jacobian
433: . Jpre - Matrix that will be used to construct the preconditioner, can be same as `J`
434: . func - Jacobian evaluation routine
435: - ctx - [optional] user-defined context for private data for the
436: Jacobian evaluation routine (may be `NULL`)
438: Calling sequence of `func`:
439: + tao - the `Tao` context
440: . x - input vector
441: . J - Jacobian matrix
442: . Jpre - matrix used to construct the preconditioner, usually the same as `J`
443: - ctx - [optional] user-defined Jacobian context
445: Level: intermediate
447: .seealso: [](ch_tao), `Tao`, `TaoSetGradient()`, `TaoSetObjective()`
448: @*/
449: PetscErrorCode TaoSetJacobianRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao tao, Vec x, Mat J, Mat Jpre, void *ctx), void *ctx)
450: {
451: PetscFunctionBegin;
453: if (J) {
455: PetscCheckSameComm(tao, 1, J, 2);
456: }
457: if (Jpre) {
459: PetscCheckSameComm(tao, 1, Jpre, 3);
460: }
461: if (ctx) tao->user_jacP = ctx;
462: if (func) tao->ops->computejacobian = func;
463: if (J) {
464: PetscCall(PetscObjectReference((PetscObject)J));
465: PetscCall(MatDestroy(&tao->jacobian));
466: tao->jacobian = J;
467: }
468: if (Jpre) {
469: PetscCall(PetscObjectReference((PetscObject)Jpre));
470: PetscCall(MatDestroy(&tao->jacobian_pre));
471: tao->jacobian_pre = Jpre;
472: }
473: PetscFunctionReturn(PETSC_SUCCESS);
474: }
476: /*@C
477: TaoSetJacobianResidualRoutine - Sets the function to compute the least-squares residual Jacobian as well as the
478: location to store the matrix.
480: Logically Collective
482: Input Parameters:
483: + tao - the `Tao` context
484: . J - Matrix used for the jacobian
485: . Jpre - Matrix that will be used to construct the preconditioner, can be same as `J`
486: . func - Jacobian evaluation routine
487: - ctx - [optional] user-defined context for private data for the
488: Jacobian evaluation routine (may be `NULL`)
490: Calling sequence of `func`:
491: + tao - the `Tao` context
492: . x - input vector
493: . J - Jacobian matrix
494: . Jpre - matrix used to construct the preconditioner, usually the same as `J`
495: - ctx - [optional] user-defined Jacobian context
497: Level: intermediate
499: .seealso: [](ch_tao), `Tao`, `TaoSetGradient()`, `TaoSetObjective()`
500: @*/
501: PetscErrorCode TaoSetJacobianResidualRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao tao, Vec x, Mat J, Mat Jpre, void *ctx), void *ctx)
502: {
503: PetscFunctionBegin;
505: if (J) {
507: PetscCheckSameComm(tao, 1, J, 2);
508: }
509: if (Jpre) {
511: PetscCheckSameComm(tao, 1, Jpre, 3);
512: }
513: if (ctx) tao->user_lsjacP = ctx;
514: if (func) tao->ops->computeresidualjacobian = func;
515: if (J) {
516: PetscCall(PetscObjectReference((PetscObject)J));
517: PetscCall(MatDestroy(&tao->ls_jac));
518: tao->ls_jac = J;
519: }
520: if (Jpre) {
521: PetscCall(PetscObjectReference((PetscObject)Jpre));
522: PetscCall(MatDestroy(&tao->ls_jac_pre));
523: tao->ls_jac_pre = Jpre;
524: }
525: PetscFunctionReturn(PETSC_SUCCESS);
526: }
528: /*@C
529: TaoSetJacobianStateRoutine - Sets the function to compute the Jacobian
530: (and its inverse) of the constraint function with respect to the state variables.
531: Used only for PDE-constrained optimization.
533: Logically Collective
535: Input Parameters:
536: + tao - the `Tao` context
537: . J - Matrix used for the Jacobian
538: . Jpre - Matrix that will be used to construct the preconditioner, can be same as `J`. Only used if `Jinv` is `NULL`
539: . Jinv - [optional] Matrix used to apply the inverse of the state Jacobian. Use `NULL` to default to PETSc `KSP` solvers to apply the inverse.
540: . func - Jacobian evaluation routine
541: - ctx - [optional] user-defined context for private data for the
542: Jacobian evaluation routine (may be `NULL`)
544: Calling sequence of `func`:
545: + tao - the `Tao` context
546: . x - input vector
547: . J - Jacobian matrix
548: . Jpre - matrix used to construct the preconditioner, usually the same as `J`
549: . Jinv - inverse of `J`
550: - ctx - [optional] user-defined Jacobian context
552: Level: intermediate
554: .seealso: [](ch_tao), `Tao`, `TaoComputeJacobianState()`, `TaoSetJacobianDesignRoutine()`, `TaoSetStateDesignIS()`
555: @*/
556: PetscErrorCode TaoSetJacobianStateRoutine(Tao tao, Mat J, Mat Jpre, Mat Jinv, PetscErrorCode (*func)(Tao tao, Vec x, Mat J, Mat Jpre, Mat Jinv, void *ctx), void *ctx)
557: {
558: PetscFunctionBegin;
560: if (J) {
562: PetscCheckSameComm(tao, 1, J, 2);
563: }
564: if (Jpre) {
566: PetscCheckSameComm(tao, 1, Jpre, 3);
567: }
568: if (Jinv) {
570: PetscCheckSameComm(tao, 1, Jinv, 4);
571: }
572: if (ctx) tao->user_jac_stateP = ctx;
573: if (func) tao->ops->computejacobianstate = func;
574: if (J) {
575: PetscCall(PetscObjectReference((PetscObject)J));
576: PetscCall(MatDestroy(&tao->jacobian_state));
577: tao->jacobian_state = J;
578: }
579: if (Jpre) {
580: PetscCall(PetscObjectReference((PetscObject)Jpre));
581: PetscCall(MatDestroy(&tao->jacobian_state_pre));
582: tao->jacobian_state_pre = Jpre;
583: }
584: if (Jinv) {
585: PetscCall(PetscObjectReference((PetscObject)Jinv));
586: PetscCall(MatDestroy(&tao->jacobian_state_inv));
587: tao->jacobian_state_inv = Jinv;
588: }
589: PetscFunctionReturn(PETSC_SUCCESS);
590: }
592: /*@C
593: TaoSetJacobianDesignRoutine - Sets the function to compute the Jacobian of
594: the constraint function with respect to the design variables. Used only for
595: PDE-constrained optimization.
597: Logically Collective
599: Input Parameters:
600: + tao - the `Tao` context
601: . J - Matrix used for the Jacobian
602: . func - Jacobian evaluation routine
603: - ctx - [optional] user-defined context for private data for the
604: Jacobian evaluation routine (may be `NULL`)
606: Calling sequence of `func`:
607: + tao - the `Tao` context
608: . x - input vector
609: . J - Jacobian matrix
610: - ctx - [optional] user-defined Jacobian context
612: Level: intermediate
614: .seealso: [](ch_tao), `Tao`, `TaoComputeJacobianDesign()`, `TaoSetJacobianStateRoutine()`, `TaoSetStateDesignIS()`
615: @*/
616: PetscErrorCode TaoSetJacobianDesignRoutine(Tao tao, Mat J, PetscErrorCode (*func)(Tao tao, Vec x, Mat J, void *ctx), void *ctx)
617: {
618: PetscFunctionBegin;
620: if (J) {
622: PetscCheckSameComm(tao, 1, J, 2);
623: }
624: if (ctx) tao->user_jac_designP = ctx;
625: if (func) tao->ops->computejacobiandesign = func;
626: if (J) {
627: PetscCall(PetscObjectReference((PetscObject)J));
628: PetscCall(MatDestroy(&tao->jacobian_design));
629: tao->jacobian_design = J;
630: }
631: PetscFunctionReturn(PETSC_SUCCESS);
632: }
634: /*@
635: TaoSetStateDesignIS - Indicate to the `Tao` object which variables in the
636: solution vector are state variables and which are design. Only applies to
637: PDE-constrained optimization.
639: Logically Collective
641: Input Parameters:
642: + tao - The `Tao` context
643: . s_is - the index set corresponding to the state variables
644: - d_is - the index set corresponding to the design variables
646: Level: intermediate
648: .seealso: [](ch_tao), `Tao`, `TaoSetJacobianStateRoutine()`, `TaoSetJacobianDesignRoutine()`
649: @*/
650: PetscErrorCode TaoSetStateDesignIS(Tao tao, IS s_is, IS d_is)
651: {
652: PetscFunctionBegin;
653: PetscCall(PetscObjectReference((PetscObject)s_is));
654: PetscCall(ISDestroy(&tao->state_is));
655: tao->state_is = s_is;
656: PetscCall(PetscObjectReference((PetscObject)d_is));
657: PetscCall(ISDestroy(&tao->design_is));
658: tao->design_is = d_is;
659: PetscFunctionReturn(PETSC_SUCCESS);
660: }
662: /*@
663: TaoComputeJacobianEquality - Computes the Jacobian matrix that has been
664: set with `TaoSetJacobianEqualityRoutine()`.
666: Collective
668: Input Parameters:
669: + tao - the `Tao` solver context
670: - X - input vector
672: Output Parameters:
673: + J - Jacobian matrix
674: - Jpre - matrix used to construct the preconditioner, often the same as `J`
676: Level: developer
678: Notes:
679: Most users should not need to explicitly call this routine, as it
680: is used internally within the optimization algorithms.
682: .seealso: [](ch_tao), `TaoComputeObjective()`, `TaoComputeObjectiveAndGradient()`, `TaoSetJacobianStateRoutine()`, `TaoComputeJacobianDesign()`, `TaoSetStateDesignIS()`
683: @*/
684: PetscErrorCode TaoComputeJacobianEquality(Tao tao, Vec X, Mat J, Mat Jpre)
685: {
686: PetscFunctionBegin;
689: PetscCheckSameComm(tao, 1, X, 2);
690: ++tao->njac_equality;
691: PetscCall(VecLockReadPush(X));
692: PetscCall(PetscLogEventBegin(TAO_JacobianEval, tao, X, J, Jpre));
693: PetscCallBack("Tao callback Jacobian(equality)", (*tao->ops->computejacobianequality)(tao, X, J, Jpre, tao->user_jac_equalityP));
694: PetscCall(PetscLogEventEnd(TAO_JacobianEval, tao, X, J, Jpre));
695: PetscCall(VecLockReadPop(X));
696: PetscFunctionReturn(PETSC_SUCCESS);
697: }
699: /*@
700: TaoComputeJacobianInequality - Computes the Jacobian matrix that has been
701: set with `TaoSetJacobianInequalityRoutine()`.
703: Collective
705: Input Parameters:
706: + tao - the `Tao` solver context
707: - X - input vector
709: Output Parameters:
710: + J - Jacobian matrix
711: - Jpre - matrix used to construct the preconditioner
713: Level: developer
715: Note:
716: Most users should not need to explicitly call this routine, as it
717: is used internally within the minimization solvers.
719: .seealso: [](ch_tao), `Tao`, `TaoComputeObjective()`, `TaoComputeObjectiveAndGradient()`, `TaoSetJacobianStateRoutine()`, `TaoComputeJacobianDesign()`, `TaoSetStateDesignIS()`
720: @*/
721: PetscErrorCode TaoComputeJacobianInequality(Tao tao, Vec X, Mat J, Mat Jpre)
722: {
723: PetscFunctionBegin;
726: PetscCheckSameComm(tao, 1, X, 2);
727: ++tao->njac_inequality;
728: PetscCall(VecLockReadPush(X));
729: PetscCall(PetscLogEventBegin(TAO_JacobianEval, tao, X, J, Jpre));
730: PetscCallBack("Tao callback Jacobian (inequality)", (*tao->ops->computejacobianinequality)(tao, X, J, Jpre, tao->user_jac_inequalityP));
731: PetscCall(PetscLogEventEnd(TAO_JacobianEval, tao, X, J, Jpre));
732: PetscCall(VecLockReadPop(X));
733: PetscFunctionReturn(PETSC_SUCCESS);
734: }
736: /*@C
737: TaoSetJacobianEqualityRoutine - Sets the function to compute the Jacobian
738: (and its inverse) of the constraint function with respect to the equality variables.
739: Used only for PDE-constrained optimization.
741: Logically Collective
743: Input Parameters:
744: + tao - the `Tao` context
745: . J - Matrix used for the Jacobian
746: . Jpre - Matrix that will be used to construct the preconditioner, can be same as `J`.
747: . func - Jacobian evaluation routine
748: - ctx - [optional] user-defined context for private data for the
749: Jacobian evaluation routine (may be `NULL`)
751: Calling sequence of `func`:
752: + tao - the `Tao` context
753: . x - input vector
754: . J - Jacobian matrix
755: . Jpre - matrix used to construct the preconditioner, usually the same as `J`
756: - ctx - [optional] user-defined Jacobian context
758: Level: intermediate
760: .seealso: [](ch_tao), `Tao`, `TaoComputeJacobianEquality()`, `TaoSetJacobianDesignRoutine()`, `TaoSetEqualityDesignIS()`
761: @*/
762: PetscErrorCode TaoSetJacobianEqualityRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao tao, Vec x, Mat J, Mat Jpre, void *ctx), void *ctx)
763: {
764: PetscFunctionBegin;
766: if (J) {
768: PetscCheckSameComm(tao, 1, J, 2);
769: }
770: if (Jpre) {
772: PetscCheckSameComm(tao, 1, Jpre, 3);
773: }
774: if (ctx) tao->user_jac_equalityP = ctx;
775: if (func) tao->ops->computejacobianequality = func;
776: if (J) {
777: PetscCall(PetscObjectReference((PetscObject)J));
778: PetscCall(MatDestroy(&tao->jacobian_equality));
779: tao->jacobian_equality = J;
780: }
781: if (Jpre) {
782: PetscCall(PetscObjectReference((PetscObject)Jpre));
783: PetscCall(MatDestroy(&tao->jacobian_equality_pre));
784: tao->jacobian_equality_pre = Jpre;
785: }
786: PetscFunctionReturn(PETSC_SUCCESS);
787: }
789: /*@C
790: TaoSetJacobianInequalityRoutine - Sets the function to compute the Jacobian
791: (and its inverse) of the constraint function with respect to the inequality variables.
792: Used only for PDE-constrained optimization.
794: Logically Collective
796: Input Parameters:
797: + tao - the `Tao` context
798: . J - Matrix used for the Jacobian
799: . Jpre - Matrix that will be used to construct the preconditioner, can be same as `J`.
800: . func - Jacobian evaluation routine
801: - ctx - [optional] user-defined context for private data for the
802: Jacobian evaluation routine (may be `NULL`)
804: Calling sequence of `func`:
805: + tao - the `Tao` context
806: . x - input vector
807: . J - Jacobian matrix
808: . Jpre - matrix used to construct the preconditioner, usually the same as `J`
809: - ctx - [optional] user-defined Jacobian context
811: Level: intermediate
813: .seealso: [](ch_tao), `Tao`, `TaoComputeJacobianInequality()`, `TaoSetJacobianDesignRoutine()`, `TaoSetInequalityDesignIS()`
814: @*/
815: PetscErrorCode TaoSetJacobianInequalityRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao tao, Vec x, Mat J, Mat Jpre, void *ctx), void *ctx)
816: {
817: PetscFunctionBegin;
819: if (J) {
821: PetscCheckSameComm(tao, 1, J, 2);
822: }
823: if (Jpre) {
825: PetscCheckSameComm(tao, 1, Jpre, 3);
826: }
827: if (ctx) tao->user_jac_inequalityP = ctx;
828: if (func) tao->ops->computejacobianinequality = func;
829: if (J) {
830: PetscCall(PetscObjectReference((PetscObject)J));
831: PetscCall(MatDestroy(&tao->jacobian_inequality));
832: tao->jacobian_inequality = J;
833: }
834: if (Jpre) {
835: PetscCall(PetscObjectReference((PetscObject)Jpre));
836: PetscCall(MatDestroy(&tao->jacobian_inequality_pre));
837: tao->jacobian_inequality_pre = Jpre;
838: }
839: PetscFunctionReturn(PETSC_SUCCESS);
840: }