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: }