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 operated on by 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: $    func(Tao tao,Vec x,Mat H,Mat Hpre,void *ctx);

 19: +  tao  - the Tao  context
 20: .  x    - input vector
 21: .  H    - Hessian matrix
 22: .  Hpre - preconditioner matrix, usually the same as H
 23: -  ctx  - [optional] user-defined Hessian context

 25:    Level: beginner

 27: .seealso: `Tao`, `TaoTypes`, `TaoSetObjective()`, `TaoSetGradient()`, `TaoSetObjectiveAndGradient()`, `TaoGetHessian()`
 28: @*/
 29: PetscErrorCode TaoSetHessian(Tao tao, Mat H, Mat Hpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, void *), void *ctx)
 30: {
 32:   if (H) {
 35:   }
 36:   if (Hpre) {
 39:   }
 40:   if (ctx) tao->user_hessP = ctx;
 41:   if (func) tao->ops->computehessian = func;
 42:   if (H) {
 43:     PetscObjectReference((PetscObject)H);
 44:     MatDestroy(&tao->hessian);
 45:     tao->hessian = H;
 46:   }
 47:   if (Hpre) {
 48:     PetscObjectReference((PetscObject)Hpre);
 49:     MatDestroy(&tao->hessian_pre);
 50:     tao->hessian_pre = Hpre;
 51:   }
 52:   return 0;
 53: }

 55: /*@C
 56:    TaoGetHessian - Gets the function to compute the Hessian as well as the location to store the matrix.

 58:    Not collective

 60:    Input Parameter:
 61: .  tao  - the Tao context

 63:    OutputParameters:
 64: +  H    - Matrix used for the hessian
 65: .  Hpre - Matrix that will be used operated on by preconditioner, can be the same as H
 66: .  func - Hessian evaluation routine
 67: -  ctx  - user-defined context for private data for the Hessian evaluation routine

 69:    Calling sequence of func:
 70: $    func(Tao tao,Vec x,Mat H,Mat Hpre,void *ctx);

 72: +  tao  - the Tao  context
 73: .  x    - input vector
 74: .  H    - Hessian matrix
 75: .  Hpre - preconditioner matrix, usually the same as H
 76: -  ctx  - [optional] user-defined Hessian context

 78:    Level: beginner

 80: .seealso: `Tao`, TaoType`, `TaoGetObjective()`, `TaoGetGradient()`, `TaoGetObjectiveAndGradient()`, `TaoSetHessian()`
 81: @*/
 82: PetscErrorCode TaoGetHessian(Tao tao, Mat *H, Mat *Hpre, PetscErrorCode (**func)(Tao, Vec, Mat, Mat, void *), void **ctx)
 83: {
 85:   if (H) *H = tao->hessian;
 86:   if (Hpre) *Hpre = tao->hessian_pre;
 87:   if (ctx) *ctx = tao->user_hessP;
 88:   if (func) *func = tao->ops->computehessian;
 89:   return 0;
 90: }

 92: PetscErrorCode TaoTestHessian(Tao tao)
 93: {
 94:   Mat               A, B, C, D, hessian;
 95:   Vec               x = tao->solution;
 96:   PetscReal         nrm, gnorm;
 97:   PetscReal         threshold = 1.e-5;
 98:   PetscInt          m, n, M, N;
 99:   PetscBool         complete_print = PETSC_FALSE, test = PETSC_FALSE, flg;
100:   PetscViewer       viewer, mviewer;
101:   MPI_Comm          comm;
102:   PetscInt          tabs;
103:   static PetscBool  directionsprinted = PETSC_FALSE;
104:   PetscViewerFormat format;

106:   PetscObjectOptionsBegin((PetscObject)tao);
107:   PetscOptionsName("-tao_test_hessian", "Compare hand-coded and finite difference Hessians", "None", &test);
108:   PetscOptionsReal("-tao_test_hessian", "Threshold for element difference between hand-coded and finite difference being meaningful", "None", threshold, &threshold, NULL);
109:   PetscOptionsViewer("-tao_test_hessian_view", "View difference between hand-coded and finite difference Hessians element entries", "None", &mviewer, &format, &complete_print);
110:   PetscOptionsEnd();
111:   if (!test) return 0;

113:   PetscObjectGetComm((PetscObject)tao, &comm);
114:   PetscViewerASCIIGetStdout(comm, &viewer);
115:   PetscViewerASCIIGetTab(viewer, &tabs);
116:   PetscViewerASCIISetTab(viewer, ((PetscObject)tao)->tablevel);
117:   PetscViewerASCIIPrintf(viewer, "  ---------- Testing Hessian -------------\n");
118:   if (!complete_print && !directionsprinted) {
119:     PetscViewerASCIIPrintf(viewer, "  Run with -tao_test_hessian_view and optionally -tao_test_hessian <threshold> to show difference\n");
120:     PetscViewerASCIIPrintf(viewer, "    of hand-coded and finite difference Hessian entries greater than <threshold>.\n");
121:   }
122:   if (!directionsprinted) {
123:     PetscViewerASCIIPrintf(viewer, "  Testing hand-coded Hessian, if (for double precision runs) ||J - Jfd||_F/||J||_F is\n");
124:     PetscViewerASCIIPrintf(viewer, "    O(1.e-8), the hand-coded Hessian is probably correct.\n");
125:     directionsprinted = PETSC_TRUE;
126:   }
127:   if (complete_print) PetscViewerPushFormat(mviewer, format);

129:   PetscObjectTypeCompare((PetscObject)tao->hessian, MATMFFD, &flg);
130:   if (!flg) hessian = tao->hessian;
131:   else hessian = tao->hessian_pre;

133:   while (hessian) {
134:     PetscObjectBaseTypeCompareAny((PetscObject)hessian, &flg, MATSEQAIJ, MATMPIAIJ, MATSEQDENSE, MATMPIDENSE, MATSEQBAIJ, MATMPIBAIJ, MATSEQSBAIJ, MATMPIBAIJ, "");
135:     if (flg) {
136:       A = hessian;
137:       PetscObjectReference((PetscObject)A);
138:     } else {
139:       MatComputeOperator(hessian, MATAIJ, &A);
140:     }

142:     MatCreate(PetscObjectComm((PetscObject)A), &B);
143:     MatGetSize(A, &M, &N);
144:     MatGetLocalSize(A, &m, &n);
145:     MatSetSizes(B, m, n, M, N);
146:     MatSetType(B, ((PetscObject)A)->type_name);
147:     MatSetUp(B);
148:     MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

150:     TaoDefaultComputeHessian(tao, x, B, B, NULL);

152:     MatDuplicate(B, MAT_COPY_VALUES, &D);
153:     MatAYPX(D, -1.0, A, DIFFERENT_NONZERO_PATTERN);
154:     MatNorm(D, NORM_FROBENIUS, &nrm);
155:     MatNorm(A, NORM_FROBENIUS, &gnorm);
156:     MatDestroy(&D);
157:     if (!gnorm) gnorm = 1; /* just in case */
158:     PetscViewerASCIIPrintf(viewer, "  ||H - Hfd||_F/||H||_F = %g, ||H - Hfd||_F = %g\n", (double)(nrm / gnorm), (double)nrm);

160:     if (complete_print) {
161:       PetscViewerASCIIPrintf(viewer, "  Hand-coded Hessian ----------\n");
162:       MatView(A, mviewer);
163:       PetscViewerASCIIPrintf(viewer, "  Finite difference Hessian ----------\n");
164:       MatView(B, mviewer);
165:     }

167:     if (complete_print) {
168:       PetscInt           Istart, Iend, *ccols, bncols, cncols, j, row;
169:       PetscScalar       *cvals;
170:       const PetscInt    *bcols;
171:       const PetscScalar *bvals;

173:       MatAYPX(B, -1.0, A, DIFFERENT_NONZERO_PATTERN);
174:       MatCreate(PetscObjectComm((PetscObject)A), &C);
175:       MatSetSizes(C, m, n, M, N);
176:       MatSetType(C, ((PetscObject)A)->type_name);
177:       MatSetUp(C);
178:       MatSetOption(C, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
179:       MatGetOwnershipRange(B, &Istart, &Iend);

181:       for (row = Istart; row < Iend; row++) {
182:         MatGetRow(B, row, &bncols, &bcols, &bvals);
183:         PetscMalloc2(bncols, &ccols, bncols, &cvals);
184:         for (j = 0, cncols = 0; j < bncols; j++) {
185:           if (PetscAbsScalar(bvals[j]) > threshold) {
186:             ccols[cncols] = bcols[j];
187:             cvals[cncols] = bvals[j];
188:             cncols += 1;
189:           }
190:         }
191:         if (cncols) MatSetValues(C, 1, &row, cncols, ccols, cvals, INSERT_VALUES);
192:         MatRestoreRow(B, row, &bncols, &bcols, &bvals);
193:         PetscFree2(ccols, cvals);
194:       }
195:       MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
196:       MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
197:       PetscViewerASCIIPrintf(viewer, "  Finite-difference minus hand-coded Hessian with tolerance %g ----------\n", (double)threshold);
198:       MatView(C, mviewer);
199:       MatDestroy(&C);
200:     }
201:     MatDestroy(&A);
202:     MatDestroy(&B);

204:     if (hessian != tao->hessian_pre) {
205:       hessian = tao->hessian_pre;
206:       PetscViewerASCIIPrintf(viewer, "  ---------- Testing Hessian for preconditioner -------------\n");
207:     } else hessian = NULL;
208:   }
209:   if (complete_print) {
210:     PetscViewerPopFormat(mviewer);
211:     PetscViewerDestroy(&mviewer);
212:   }
213:   PetscViewerASCIISetTab(viewer, tabs);
214:   return 0;
215: }

217: /*@C
218:    TaoComputeHessian - Computes the Hessian matrix that has been
219:    set with `TaoSetHessian()`.

221:    Collective

223:    Input Parameters:
224: +  tao - the Tao solver context
225: -  X   - input vector

227:    Output Parameters:
228: +  H    - Hessian matrix
229: -  Hpre - Preconditioning matrix

231:    Options Database Keys:
232: +     -tao_test_hessian - compare the user provided Hessian with one compute via finite differences to check for errors
233: .     -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
234: -     -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

236:    Notes:
237:    Most users should not need to explicitly call this routine, as it
238:    is used internally within the minimization solvers.

240:    `TaoComputeHessian()` is typically used within optimization algorithms,
241:    so most users would not generally call this routine
242:    themselves.

244:    Developer Note:
245:    The Hessian test mechanism follows `SNESTestJacobian()`.

247:    Level: developer

249: .seealso: `Tao`, `TaoComputeObjective()`, `TaoComputeObjectiveAndGradient()`, `TaoSetHessian()`
250: @*/
251: PetscErrorCode TaoComputeHessian(Tao tao, Vec X, Mat H, Mat Hpre)
252: {

258:   ++tao->nhess;
259:   VecLockReadPush(X);
260:   PetscLogEventBegin(TAO_HessianEval, tao, X, H, Hpre);
261:   PetscCallBack("Tao callback Hessian", (*tao->ops->computehessian)(tao, X, H, Hpre, tao->user_hessP));
262:   PetscLogEventEnd(TAO_HessianEval, tao, X, H, Hpre);
263:   VecLockReadPop(X);

265:   TaoTestHessian(tao);
266:   return 0;
267: }

269: /*@C
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:    Notes:
284:    Most users should not need to explicitly call this routine, as it
285:    is used internally within the minimization solvers.

287:    `TaoComputeJacobian()` is typically used within minimization
288:    implementations, so most users would not generally call this routine
289:    themselves.

291:    Level: developer

293: .seealso: `TaoComputeObjective()`, `TaoComputeObjectiveAndGradient()`, `TaoSetJacobianRoutine()`
294: @*/
295: PetscErrorCode TaoComputeJacobian(Tao tao, Vec X, Mat J, Mat Jpre)
296: {
300:   ++tao->njac;
301:   VecLockReadPush(X);
302:   PetscLogEventBegin(TAO_JacobianEval, tao, X, J, Jpre);
303:   PetscCallBack("Tao callback Jacobian", (*tao->ops->computejacobian)(tao, X, J, Jpre, tao->user_jacP));
304:   PetscLogEventEnd(TAO_JacobianEval, tao, X, J, Jpre);
305:   VecLockReadPop(X);
306:   return 0;
307: }

309: /*@C
310:    TaoComputeResidualJacobian - Computes the least-squares residual Jacobian matrix that has been
311:    set with `TaoSetJacobianResidual()`.

313:    Collective

315:    Input Parameters:
316: +  tao - the Tao solver context
317: -  X   - input vector

319:    Output Parameters:
320: +  J    - Jacobian matrix
321: -  Jpre - Preconditioning matrix

323:    Notes:
324:    Most users should not need to explicitly call this routine, as it
325:    is used internally within the minimization solvers.

327:    `TaoComputeResidualJacobian()` is typically used within least-squares
328:    implementations, so most users would not generally call this routine
329:    themselves.

331:    Level: developer

333: .seealso: `Tao`, `TaoComputeResidual()`, `TaoSetJacobianResidual()`
334: @*/
335: PetscErrorCode TaoComputeResidualJacobian(Tao tao, Vec X, Mat J, Mat Jpre)
336: {
340:   ++tao->njac;
341:   VecLockReadPush(X);
342:   PetscLogEventBegin(TAO_JacobianEval, tao, X, J, Jpre);
343:   PetscCallBack("Tao callback least-squares residual Jacobian", (*tao->ops->computeresidualjacobian)(tao, X, J, Jpre, tao->user_lsjacP));
344:   PetscLogEventEnd(TAO_JacobianEval, tao, X, J, Jpre);
345:   VecLockReadPop(X);
346:   return 0;
347: }

349: /*@C
350:    TaoComputeJacobianState - Computes the Jacobian matrix that has been
351:    set with `TaoSetJacobianStateRoutine()`.

353:    Collective

355:    Input Parameters:
356: +  tao - the Tao solver context
357: -  X   - input vector

359:    Output Parameters:
360: +  J    - Jacobian matrix
361: .  Jpre - Preconditioning matrix
362: -  Jinv - unknown

364:    Notes:
365:    Most users should not need to explicitly call this routine, as it
366:    is used internally within the optimization algorithms.

368:    Level: developer

370: .seealso: `Tao`, `TaoComputeObjective()`, `TaoComputeObjectiveAndGradient()`, `TaoSetJacobianStateRoutine()`, `TaoComputeJacobianDesign()`, `TaoSetStateDesignIS()`
371: @*/
372: PetscErrorCode TaoComputeJacobianState(Tao tao, Vec X, Mat J, Mat Jpre, Mat Jinv)
373: {
377:   ++tao->njac_state;
378:   VecLockReadPush(X);
379:   PetscLogEventBegin(TAO_JacobianEval, tao, X, J, Jpre);
380:   PetscCallBack("Tao callback Jacobian(state)", (*tao->ops->computejacobianstate)(tao, X, J, Jpre, Jinv, tao->user_jac_stateP));
381:   PetscLogEventEnd(TAO_JacobianEval, tao, X, J, Jpre);
382:   VecLockReadPop(X);
383:   return 0;
384: }

386: /*@C
387:    TaoComputeJacobianDesign - Computes the Jacobian matrix that has been
388:    set with `TaoSetJacobianDesignRoutine()`.

390:    Collective

392:    Input Parameters:
393: +  tao - the Tao solver context
394: -  X   - input vector

396:    Output Parameters:
397: .  J - Jacobian matrix

399:    Notes:
400:    Most users should not need to explicitly call this routine, as it
401:    is used internally within the optimization algorithms.

403:    Level: developer

405: .seealso: `Tao`, `TaoComputeObjective()`, `TaoComputeObjectiveAndGradient()`, `TaoSetJacobianDesignRoutine()`, `TaoComputeJacobianDesign()`, `TaoSetStateDesignIS()`
406: @*/
407: PetscErrorCode TaoComputeJacobianDesign(Tao tao, Vec X, Mat J)
408: {
412:   ++tao->njac_design;
413:   VecLockReadPush(X);
414:   PetscLogEventBegin(TAO_JacobianEval, tao, X, J, NULL);
415:   PetscCallBack("Tao callback Jacobian(design)", (*tao->ops->computejacobiandesign)(tao, X, J, tao->user_jac_designP));
416:   PetscLogEventEnd(TAO_JacobianEval, tao, X, J, NULL);
417:   VecLockReadPop(X);
418:   return 0;
419: }

421: /*@C
422:    TaoSetJacobianRoutine - Sets the function to compute the Jacobian as well as the location to store the matrix.

424:    Logically collective

426:    Input Parameters:
427: +  tao  - the Tao context
428: .  J    - Matrix used for the jacobian
429: .  Jpre - Matrix that will be used operated on by preconditioner, can be same as J
430: .  func - Jacobian evaluation routine
431: -  ctx  - [optional] user-defined context for private data for the
432:           Jacobian evaluation routine (may be NULL)

434:    Calling sequence of func:
435: $    func(Tao tao,Vec x,Mat J,Mat Jpre,void *ctx);

437: +  tao  - the Tao  context
438: .  x    - input vector
439: .  J    - Jacobian matrix
440: .  Jpre - preconditioning matrix, usually the same as J
441: -  ctx  - [optional] user-defined Jacobian context

443:    Level: intermediate

445: .seealso: `Tao`, `TaoSetGradient()`, `TaoSetObjective()`
446: @*/
447: PetscErrorCode TaoSetJacobianRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, void *), void *ctx)
448: {
450:   if (J) {
453:   }
454:   if (Jpre) {
457:   }
458:   if (ctx) tao->user_jacP = ctx;
459:   if (func) tao->ops->computejacobian = func;
460:   if (J) {
461:     PetscObjectReference((PetscObject)J);
462:     MatDestroy(&tao->jacobian);
463:     tao->jacobian = J;
464:   }
465:   if (Jpre) {
466:     PetscObjectReference((PetscObject)Jpre);
467:     MatDestroy(&tao->jacobian_pre);
468:     tao->jacobian_pre = Jpre;
469:   }
470:   return 0;
471: }

473: /*@C
474:    TaoSetJacobianResidualRoutine - Sets the function to compute the least-squares residual Jacobian as well as the
475:    location to store the matrix.

477:    Logically collective

479:    Input Parameters:
480: +  tao  - the Tao context
481: .  J    - Matrix used for the jacobian
482: .  Jpre - Matrix that will be used operated on by preconditioner, can be same as J
483: .  func - Jacobian evaluation routine
484: -  ctx  - [optional] user-defined context for private data for the
485:           Jacobian evaluation routine (may be NULL)

487:    Calling sequence of func:
488: $    func(Tao tao,Vec x,Mat J,Mat Jpre,void *ctx);

490: +  tao  - the Tao  context
491: .  x    - input vector
492: .  J    - Jacobian matrix
493: .  Jpre - preconditioning matrix, usually the same as J
494: -  ctx  - [optional] user-defined Jacobian context

496:    Level: intermediate

498: .seealso: `Tao`, `TaoSetGradient()`, `TaoSetObjective()`
499: @*/
500: PetscErrorCode TaoSetJacobianResidualRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, void *), void *ctx)
501: {
503:   if (J) {
506:   }
507:   if (Jpre) {
510:   }
511:   if (ctx) tao->user_lsjacP = ctx;
512:   if (func) tao->ops->computeresidualjacobian = func;
513:   if (J) {
514:     PetscObjectReference((PetscObject)J);
515:     MatDestroy(&tao->ls_jac);
516:     tao->ls_jac = J;
517:   }
518:   if (Jpre) {
519:     PetscObjectReference((PetscObject)Jpre);
520:     MatDestroy(&tao->ls_jac_pre);
521:     tao->ls_jac_pre = Jpre;
522:   }
523:   return 0;
524: }

526: /*@C
527:    TaoSetJacobianStateRoutine - Sets the function to compute the Jacobian
528:    (and its inverse) of the constraint function with respect to the state variables.
529:    Used only for PDE-constrained optimization.

531:    Logically collective

533:    Input Parameters:
534: +  tao  - the Tao context
535: .  J    - Matrix used for the jacobian
536: .  Jpre - Matrix that will be used operated on by PETSc preconditioner, can be same as J.  Only used if Jinv is NULL
537: .  Jinv - [optional] Matrix used to apply the inverse of the state jacobian. Use NULL to default to PETSc KSP solvers to apply the inverse.
538: .  func - Jacobian evaluation routine
539: -  ctx  - [optional] user-defined context for private data for the
540:           Jacobian evaluation routine (may be NULL)

542:    Calling sequence of func:
543: $    func(Tao tao,Vec x,Mat J,Mat Jpre,Mat Jinv,void *ctx);

545: +  tao  - the Tao  context
546: .  x    - input vector
547: .  J    - Jacobian matrix
548: .  Jpre - preconditioner matrix, usually the same as J
549: .  Jinv - inverse of J
550: -  ctx  - [optional] user-defined Jacobian context

552:    Level: intermediate

554: .seealso: `Tao`, `TaoComputeJacobianState()`, `TaoSetJacobianDesignRoutine()`, `TaoSetStateDesignIS()`
555: @*/
556: PetscErrorCode TaoSetJacobianStateRoutine(Tao tao, Mat J, Mat Jpre, Mat Jinv, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, Mat, void *), void *ctx)
557: {
559:   if (J) {
562:   }
563:   if (Jpre) {
566:   }
567:   if (Jinv) {
570:   }
571:   if (ctx) tao->user_jac_stateP = ctx;
572:   if (func) tao->ops->computejacobianstate = func;
573:   if (J) {
574:     PetscObjectReference((PetscObject)J);
575:     MatDestroy(&tao->jacobian_state);
576:     tao->jacobian_state = J;
577:   }
578:   if (Jpre) {
579:     PetscObjectReference((PetscObject)Jpre);
580:     MatDestroy(&tao->jacobian_state_pre);
581:     tao->jacobian_state_pre = Jpre;
582:   }
583:   if (Jinv) {
584:     PetscObjectReference((PetscObject)Jinv);
585:     MatDestroy(&tao->jacobian_state_inv);
586:     tao->jacobian_state_inv = Jinv;
587:   }
588:   return 0;
589: }

591: /*@C
592:    TaoSetJacobianDesignRoutine - Sets the function to compute the Jacobian of
593:    the constraint function with respect to the design variables.  Used only for
594:    PDE-constrained optimization.

596:    Logically collective

598:    Input Parameters:
599: +  tao  - the Tao context
600: .  J    - Matrix used for the jacobian
601: .  func - Jacobian evaluation routine
602: -  ctx  - [optional] user-defined context for private data for the
603:           Jacobian evaluation routine (may be NULL)

605:    Calling sequence of func:
606: $    func(Tao tao,Vec x,Mat J,void *ctx);

608: +  tao - the Tao  context
609: .  x   - input vector
610: .  J   - Jacobian matrix
611: -  ctx - [optional] user-defined Jacobian context

613:    Level: intermediate

615: .seealso: `Tao`, `TaoComputeJacobianDesign()`, `TaoSetJacobianStateRoutine()`, `TaoSetStateDesignIS()`
616: @*/
617: PetscErrorCode TaoSetJacobianDesignRoutine(Tao tao, Mat J, PetscErrorCode (*func)(Tao, Vec, Mat, void *), void *ctx)
618: {
620:   if (J) {
623:   }
624:   if (ctx) tao->user_jac_designP = ctx;
625:   if (func) tao->ops->computejacobiandesign = func;
626:   if (J) {
627:     PetscObjectReference((PetscObject)J);
628:     MatDestroy(&tao->jacobian_design);
629:     tao->jacobian_design = J;
630:   }
631:   return 0;
632: }

634: /*@
635:    TaoSetStateDesignIS - Indicate to the Tao 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: `Tao`, `TaoSetJacobianStateRoutine()`, `TaoSetJacobianDesignRoutine()`
649: @*/
650: PetscErrorCode TaoSetStateDesignIS(Tao tao, IS s_is, IS d_is)
651: {
652:   PetscObjectReference((PetscObject)s_is);
653:   ISDestroy(&tao->state_is);
654:   tao->state_is = s_is;
655:   PetscObjectReference((PetscObject)(d_is));
656:   ISDestroy(&tao->design_is);
657:   tao->design_is = d_is;
658:   return 0;
659: }

661: /*@C
662:    TaoComputeJacobianEquality - Computes the Jacobian matrix that has been
663:    set with `TaoSetJacobianEqualityRoutine()`.

665:    Collective

667:    Input Parameters:
668: +  tao - the Tao solver context
669: -  X   - input vector

671:    Output Parameters:
672: +  J    - Jacobian matrix
673: -  Jpre - Preconditioning matrix

675:    Notes:
676:    Most users should not need to explicitly call this routine, as it
677:    is used internally within the optimization algorithms.

679:    Level: developer

681: .seealso: `TaoComputeObjective()`, `TaoComputeObjectiveAndGradient()`, `TaoSetJacobianStateRoutine()`, `TaoComputeJacobianDesign()`, `TaoSetStateDesignIS()`
682: @*/
683: PetscErrorCode TaoComputeJacobianEquality(Tao tao, Vec X, Mat J, Mat Jpre)
684: {
688:   ++tao->njac_equality;
689:   VecLockReadPush(X);
690:   PetscLogEventBegin(TAO_JacobianEval, tao, X, J, Jpre);
691:   PetscCallBack("Tao callback Jacobian(equality)", (*tao->ops->computejacobianequality)(tao, X, J, Jpre, tao->user_jac_equalityP));
692:   PetscLogEventEnd(TAO_JacobianEval, tao, X, J, Jpre);
693:   VecLockReadPop(X);
694:   return 0;
695: }

697: /*@C
698:    TaoComputeJacobianInequality - Computes the Jacobian matrix that has been
699:    set with `TaoSetJacobianInequalityRoutine()`.

701:    Collective

703:    Input Parameters:
704: +  tao - the Tao solver context
705: -  X   - input vector

707:    Output Parameters:
708: +  J    - Jacobian matrix
709: -  Jpre - Preconditioning matrix

711:    Notes:
712:    Most users should not need to explicitly call this routine, as it
713:    is used internally within the minimization solvers.

715:    Level: developer

717: .seealso: `Tao`, `TaoComputeObjective()`, `TaoComputeObjectiveAndGradient()`, `TaoSetJacobianStateRoutine()`, `TaoComputeJacobianDesign()`, `TaoSetStateDesignIS()`
718: @*/
719: PetscErrorCode TaoComputeJacobianInequality(Tao tao, Vec X, Mat J, Mat Jpre)
720: {
724:   ++tao->njac_inequality;
725:   VecLockReadPush(X);
726:   PetscLogEventBegin(TAO_JacobianEval, tao, X, J, Jpre);
727:   PetscCallBack("Tao callback Jacobian (inequality)", (*tao->ops->computejacobianinequality)(tao, X, J, Jpre, tao->user_jac_inequalityP));
728:   PetscLogEventEnd(TAO_JacobianEval, tao, X, J, Jpre);
729:   VecLockReadPop(X);
730:   return 0;
731: }

733: /*@C
734:    TaoSetJacobianEqualityRoutine - Sets the function to compute the Jacobian
735:    (and its inverse) of the constraint function with respect to the equality variables.
736:    Used only for PDE-constrained optimization.

738:    Logically collective

740:    Input Parameters:
741: +  tao  - the Tao context
742: .  J    - Matrix used for the jacobian
743: .  Jpre - Matrix that will be used operated on by PETSc preconditioner, can be same as J.
744: .  func - Jacobian evaluation routine
745: -  ctx  - [optional] user-defined context for private data for the
746:           Jacobian evaluation routine (may be NULL)

748:    Calling sequence of func:
749: $    func(Tao tao,Vec x,Mat J,Mat Jpre,void *ctx);

751: +  tao  - the Tao  context
752: .  x    - input vector
753: .  J    - Jacobian matrix
754: .  Jpre - preconditioner matrix, usually the same as J
755: -  ctx  - [optional] user-defined Jacobian context

757:    Level: intermediate

759: .seealso: `Tao`, `TaoComputeJacobianEquality()`, `TaoSetJacobianDesignRoutine()`, `TaoSetEqualityDesignIS()`
760: @*/
761: PetscErrorCode TaoSetJacobianEqualityRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, void *), void *ctx)
762: {
764:   if (J) {
767:   }
768:   if (Jpre) {
771:   }
772:   if (ctx) tao->user_jac_equalityP = ctx;
773:   if (func) tao->ops->computejacobianequality = func;
774:   if (J) {
775:     PetscObjectReference((PetscObject)J);
776:     MatDestroy(&tao->jacobian_equality);
777:     tao->jacobian_equality = J;
778:   }
779:   if (Jpre) {
780:     PetscObjectReference((PetscObject)Jpre);
781:     MatDestroy(&tao->jacobian_equality_pre);
782:     tao->jacobian_equality_pre = Jpre;
783:   }
784:   return 0;
785: }

787: /*@C
788:    TaoSetJacobianInequalityRoutine - Sets the function to compute the Jacobian
789:    (and its inverse) of the constraint function with respect to the inequality variables.
790:    Used only for PDE-constrained optimization.

792:    Logically collective

794:    Input Parameters:
795: +  tao  - the Tao context
796: .  J    - Matrix used for the jacobian
797: .  Jpre - Matrix that will be used operated on by PETSc preconditioner, can be same as J.
798: .  func - Jacobian evaluation routine
799: -  ctx  - [optional] user-defined context for private data for the
800:           Jacobian evaluation routine (may be NULL)

802:    Calling sequence of func:
803: $    func(Tao tao,Vec x,Mat J,Mat Jpre,void *ctx);

805: +  tao  - the Tao  context
806: .  x    - input vector
807: .  J    - Jacobian matrix
808: .  Jpre - preconditioner matrix, usually the same as J
809: -  ctx  - [optional] user-defined Jacobian context

811:    Level: intermediate

813: .seealso: `Tao`, `TaoComputeJacobianInequality()`, `TaoSetJacobianDesignRoutine()`, `TaoSetInequalityDesignIS()`
814: @*/
815: PetscErrorCode TaoSetJacobianInequalityRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, void *), void *ctx)
816: {
818:   if (J) {
821:   }
822:   if (Jpre) {
825:   }
826:   if (ctx) tao->user_jac_inequalityP = ctx;
827:   if (func) tao->ops->computejacobianinequality = func;
828:   if (J) {
829:     PetscObjectReference((PetscObject)J);
830:     MatDestroy(&tao->jacobian_inequality);
831:     tao->jacobian_inequality = J;
832:   }
833:   if (Jpre) {
834:     PetscObjectReference((PetscObject)Jpre);
835:     MatDestroy(&tao->jacobian_inequality_pre);
836:     tao->jacobian_inequality_pre = Jpre;
837:   }
838:   return 0;
839: }