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, MATMPISBAIJ, ""));
    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 - matrix used to construct the preconditioner, usually the same as `H`

    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 - matrix used to compute the preconditioner, often the same as `J`

    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 - matrix used to compute the preconditioner, often the same as `J`

    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:   TaoGetJacobianEqualityRoutine - Gets the function used to compute equality constraint Jacobian.

    792:   Not Collective

    794:   Input Parameter:
    795: . tao - the `Tao` context

    797:   Output Parameters:
    798: + J    - the matrix to internally hold the constraint computation
    799: . Jpre - the matrix used to construct the preconditioner
    800: . func - Jacobian evaluation routine
    801: - ctx  - the (optional) user-defined context

    803:   Calling sequence of `func`:
    804: + tao  - the `Tao` context
    805: . x    - input vector
    806: . J    - Jacobian matrix
    807: . Jpre - matrix used to construct the preconditioner, usually the same as `J`
    808: - ctx  - [optional] user-defined Jacobian context

    810:   Level: intermediate

    812: .seealso: [](ch_tao), `Tao`, `TaoComputeJacobianEquality()`, `TaoSetJacobianEqualityRoutine()`
    813: @*/
    814: PetscErrorCode TaoGetJacobianEqualityRoutine(Tao tao, Mat *J, Mat *Jpre, PetscErrorCode (**func)(Tao tao, Vec x, Mat J, Mat Jpre, void *ctx), void **ctx)
    815: {
    816:   PetscFunctionBegin;
    818:   if (J) *J = tao->jacobian_equality;
    819:   if (Jpre) *Jpre = tao->jacobian_equality_pre;
    820:   if (func) *func = tao->ops->computejacobianequality;
    821:   if (ctx) *ctx = tao->user_jac_equalityP;
    822:   PetscFunctionReturn(PETSC_SUCCESS);
    823: }

    825: /*@C
    826:   TaoSetJacobianInequalityRoutine - Sets the function to compute the Jacobian
    827:   (and its inverse) of the constraint function with respect to the inequality variables.
    828:   Used only for PDE-constrained optimization.

    830:   Logically Collective

    832:   Input Parameters:
    833: + tao  - the `Tao` context
    834: . J    - Matrix used for the Jacobian
    835: . Jpre - Matrix that will be used to construct the preconditioner, can be same as `J`.
    836: . func - Jacobian evaluation routine
    837: - ctx  - [optional] user-defined context for private data for the
    838:           Jacobian evaluation routine (may be `NULL`)

    840:   Calling sequence of `func`:
    841: + tao  - the `Tao` context
    842: . x    - input vector
    843: . J    - Jacobian matrix
    844: . Jpre - matrix used to construct the preconditioner, usually the same as `J`
    845: - ctx  - [optional] user-defined Jacobian context

    847:   Level: intermediate

    849: .seealso: [](ch_tao), `Tao`, `TaoComputeJacobianInequality()`, `TaoSetJacobianDesignRoutine()`, `TaoSetInequalityDesignIS()`
    850: @*/
    851: PetscErrorCode TaoSetJacobianInequalityRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao tao, Vec x, Mat J, Mat Jpre, void *ctx), void *ctx)
    852: {
    853:   PetscFunctionBegin;
    855:   if (J) {
    857:     PetscCheckSameComm(tao, 1, J, 2);
    858:   }
    859:   if (Jpre) {
    861:     PetscCheckSameComm(tao, 1, Jpre, 3);
    862:   }
    863:   if (ctx) tao->user_jac_inequalityP = ctx;
    864:   if (func) tao->ops->computejacobianinequality = func;
    865:   if (J) {
    866:     PetscCall(PetscObjectReference((PetscObject)J));
    867:     PetscCall(MatDestroy(&tao->jacobian_inequality));
    868:     tao->jacobian_inequality = J;
    869:   }
    870:   if (Jpre) {
    871:     PetscCall(PetscObjectReference((PetscObject)Jpre));
    872:     PetscCall(MatDestroy(&tao->jacobian_inequality_pre));
    873:     tao->jacobian_inequality_pre = Jpre;
    874:   }
    875:   PetscFunctionReturn(PETSC_SUCCESS);
    876: }

    878: /*@C
    879:   TaoGetJacobianInequalityRoutine - Gets the function used to compute inequality constraint Jacobian.

    881:   Not Collective

    883:   Input Parameter:
    884: . tao - the `Tao` context

    886:   Output Parameters:
    887: + J    - the matrix to internally hold the constraint computation
    888: . Jpre - the matrix used to construct the preconditioner
    889: . func - Jacobian evaluation routine
    890: - ctx  - the (optional) user-defined context

    892:   Calling sequence of `func`:
    893: + tao  - the `Tao` context
    894: . x    - input vector
    895: . J    - Jacobian matrix
    896: . Jpre - matrix used to construct the preconditioner, usually the same as `J`
    897: - ctx  - [optional] user-defined Jacobian context

    899:   Level: intermediate

    901: .seealso: [](ch_tao), `Tao`, `TaoComputeJacobianInequality()`, `TaoSetJacobianInequalityRoutine()`
    902: @*/
    903: PetscErrorCode TaoGetJacobianInequalityRoutine(Tao tao, Mat *J, Mat *Jpre, PetscErrorCode (**func)(Tao tao, Vec x, Mat J, Mat Jpre, void *ctx), void **ctx)
    904: {
    905:   PetscFunctionBegin;
    907:   if (J) *J = tao->jacobian_inequality;
    908:   if (Jpre) *Jpre = tao->jacobian_inequality_pre;
    909:   if (func) *func = tao->ops->computejacobianinequality;
    910:   if (ctx) *ctx = tao->user_jac_inequalityP;
    911:   PetscFunctionReturn(PETSC_SUCCESS);
    912: }