Actual source code: taosolver_hj.c

petsc-3.12.0 2019-09-29
Report Typos and Errors
  1:  #include <petsc/private/taoimpl.h>


  4: /*@C
  5:    TaoSetHessianRoutine - Sets the function to compute the Hessian as well as the location to store the matrix.

  7:    Logically collective on Tao

  9:    Input Parameters:
 10: +  tao  - the Tao context
 11: .  H    - Matrix used for the hessian
 12: .  Hpre - Matrix that will be used operated on by preconditioner, can be same as H
 13: .  func - Hessian evaluation routine
 14: -  ctx  - [optional] user-defined context for private data for the
 15:          Hessian evaluation routine (may be NULL)

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

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

 26:    Level: beginner
 27: @*/
 28: PetscErrorCode TaoSetHessianRoutine(Tao tao, Mat H, Mat Hpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, void*), void *ctx)
 29: {

 34:   if (H) {
 37:   }
 38:   if (Hpre) {
 41:   }
 42:   if (ctx) {
 43:     tao->user_hessP = ctx;
 44:   }
 45:   if (func) {
 46:     tao->ops->computehessian = func;
 47:   }
 48:   if (H) {
 49:     PetscObjectReference((PetscObject)H);
 50:     MatDestroy(&tao->hessian);
 51:     tao->hessian = H;
 52:   }
 53:   if (Hpre) {
 54:     PetscObjectReference((PetscObject)Hpre);
 55:     MatDestroy(&tao->hessian_pre);
 56:     tao->hessian_pre = Hpre;
 57:   }
 58:   return(0);
 59: }

 61: PetscErrorCode TaoTestHessian(Tao tao)
 62: {
 63:   Mat               A,B,C,D,hessian;
 64:   Vec               x = tao->solution;
 65:   PetscErrorCode    ierr;
 66:   PetscReal         nrm,gnorm;
 67:   PetscReal         threshold = 1.e-5;
 68:   PetscInt          m,n,M,N;
 69:   PetscBool         complete_print = PETSC_FALSE,test = PETSC_FALSE,flg;
 70:   PetscViewer       viewer,mviewer;
 71:   MPI_Comm          comm;
 72:   PetscInt          tabs;
 73:   static PetscBool  directionsprinted = PETSC_FALSE;
 74:   PetscViewerFormat format;

 77:   PetscObjectOptionsBegin((PetscObject)tao);
 78:   PetscOptionsName("-tao_test_hessian","Compare hand-coded and finite difference Hessians","None",&test);
 79:   PetscOptionsReal("-tao_test_hessian", "Threshold for element difference between hand-coded and finite difference being meaningful","None",threshold,&threshold,NULL);
 80:   PetscOptionsViewer("-tao_test_hessian_view","View difference between hand-coded and finite difference Hessians element entries","None",&mviewer,&format,&complete_print);
 81:   PetscOptionsEnd();
 82:   if (!test) return(0);

 84:   PetscObjectGetComm((PetscObject)tao,&comm);
 85:   PetscViewerASCIIGetStdout(comm,&viewer);
 86:   PetscViewerASCIIGetTab(viewer, &tabs);
 87:   PetscViewerASCIISetTab(viewer, ((PetscObject)tao)->tablevel);
 88:   PetscViewerASCIIPrintf(viewer,"  ---------- Testing Hessian -------------\n");
 89:   if (!complete_print && !directionsprinted) {
 90:     PetscViewerASCIIPrintf(viewer,"  Run with -tao_test_hessian_view and optionally -tao_test_hessian <threshold> to show difference\n");
 91:     PetscViewerASCIIPrintf(viewer,"    of hand-coded and finite difference Hessian entries greater than <threshold>.\n");
 92:   }
 93:   if (!directionsprinted) {
 94:     PetscViewerASCIIPrintf(viewer,"  Testing hand-coded Hessian, if (for double precision runs) ||J - Jfd||_F/||J||_F is\n");
 95:     PetscViewerASCIIPrintf(viewer,"    O(1.e-8), the hand-coded Hessian is probably correct.\n");
 96:     directionsprinted = PETSC_TRUE;
 97:   }
 98:   if (complete_print) {
 99:     PetscViewerPushFormat(mviewer,format);
100:   }

102:   PetscObjectTypeCompare((PetscObject)tao->hessian,MATMFFD,&flg);
103:   if (!flg) hessian = tao->hessian;
104:   else hessian = tao->hessian_pre;

106:   while (hessian) {
107:     PetscObjectBaseTypeCompareAny((PetscObject)hessian,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");
108:     if (flg) {
109:       A    = hessian;
110:       PetscObjectReference((PetscObject)A);
111:     } else {
112:       MatComputeOperator(hessian,MATAIJ,&A);
113:     }

115:     MatCreate(PetscObjectComm((PetscObject)A),&B);
116:     MatGetSize(A,&M,&N);
117:     MatGetLocalSize(A,&m,&n);
118:     MatSetSizes(B,m,n,M,N);
119:     MatSetType(B,((PetscObject)A)->type_name);
120:     MatSetUp(B);
121:     MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);

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

125:     MatDuplicate(B,MAT_COPY_VALUES,&D);
126:     MatAYPX(D,-1.0,A,DIFFERENT_NONZERO_PATTERN);
127:     MatNorm(D,NORM_FROBENIUS,&nrm);
128:     MatNorm(A,NORM_FROBENIUS,&gnorm);
129:     MatDestroy(&D);
130:     if (!gnorm) gnorm = 1; /* just in case */
131:     PetscViewerASCIIPrintf(viewer,"  ||H - Hfd||_F/||H||_F = %g, ||H - Hfd||_F = %g\n",(double)(nrm/gnorm),(double)nrm);

133:     if (complete_print) {
134:       PetscViewerASCIIPrintf(viewer,"  Hand-coded Hessian ----------\n");
135:       MatView(A,mviewer);
136:       PetscViewerASCIIPrintf(viewer,"  Finite difference Hessian ----------\n");
137:       MatView(B,mviewer);
138:     }

140:     if (complete_print) {
141:       PetscInt          Istart, Iend, *ccols, bncols, cncols, j, row;
142:       PetscScalar       *cvals;
143:       const PetscInt    *bcols;
144:       const PetscScalar *bvals;

146:       MatAYPX(B,-1.0,A,DIFFERENT_NONZERO_PATTERN);
147:       MatCreate(PetscObjectComm((PetscObject)A),&C);
148:       MatSetSizes(C,m,n,M,N);
149:       MatSetType(C,((PetscObject)A)->type_name);
150:       MatSetUp(C);
151:       MatSetOption(C,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
152:       MatGetOwnershipRange(B,&Istart,&Iend);

154:       for (row = Istart; row < Iend; row++) {
155:         MatGetRow(B,row,&bncols,&bcols,&bvals);
156:         PetscMalloc2(bncols,&ccols,bncols,&cvals);
157:         for (j = 0, cncols = 0; j < bncols; j++) {
158:           if (PetscAbsScalar(bvals[j]) > threshold) {
159:             ccols[cncols] = bcols[j];
160:             cvals[cncols] = bvals[j];
161:             cncols += 1;
162:           }
163:         }
164:         if (cncols) {
165:           MatSetValues(C,1,&row,cncols,ccols,cvals,INSERT_VALUES);
166:         }
167:         MatRestoreRow(B,row,&bncols,&bcols,&bvals);
168:         PetscFree2(ccols,cvals);
169:       }
170:       MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
171:       MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
172:       PetscViewerASCIIPrintf(viewer,"  Finite-difference minus hand-coded Hessian with tolerance %g ----------\n",(double)threshold);
173:       MatView(C,mviewer);
174:       MatDestroy(&C);
175:     }
176:     MatDestroy(&A);
177:     MatDestroy(&B);

179:     if (hessian != tao->hessian_pre) {
180:       hessian = tao->hessian_pre;
181:       PetscViewerASCIIPrintf(viewer,"  ---------- Testing Hessian for preconditioner -------------\n");
182:     } else hessian = NULL;
183:   }
184:   if (complete_print) {
185:     PetscViewerPopFormat(mviewer);
186:   }
187:   PetscViewerASCIISetTab(viewer,tabs);
188:   return(0);
189: }

191: /*@C
192:    TaoComputeHessian - Computes the Hessian matrix that has been
193:    set with TaoSetHessianRoutine().

195:    Collective on Tao

197:    Input Parameters:
198: +  tao - the Tao solver context
199: -  X   - input vector

201:    Output Parameters:
202: +  H    - Hessian matrix
203: -  Hpre - Preconditioning matrix

205:    Options Database Keys:
206: +     -tao_test_hessian - compare the user provided Hessian with one compute via finite differences to check for errors
207: .     -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
208: -     -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

210:    Notes:
211:    Most users should not need to explicitly call this routine, as it
212:    is used internally within the minimization solvers.

214:    TaoComputeHessian() is typically used within minimization
215:    implementations, so most users would not generally call this routine
216:    themselves.

218:    Developer Notes:
219:    The Hessian test mechanism follows SNESTestJacobian().

221:    Level: developer

223: .seealso: TaoComputeObjective(), TaoComputeObjectiveAndGradient(), TaoSetHessianRoutine()
224: @*/
225: PetscErrorCode TaoComputeHessian(Tao tao, Vec X, Mat H, Mat Hpre)
226: {

233:   if (!tao->ops->computehessian) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TaoSetHessianRoutine() first");
234:   ++tao->nhess;
235:   VecLockReadPush(X);
236:   PetscLogEventBegin(TAO_HessianEval,tao,X,H,Hpre);
237:   PetscStackPush("Tao user Hessian function");
238:   (*tao->ops->computehessian)(tao,X,H,Hpre,tao->user_hessP);
239:   PetscStackPop;
240:   PetscLogEventEnd(TAO_HessianEval,tao,X,H,Hpre);
241:   VecLockReadPop(X);

243:   TaoTestHessian(tao);
244:   return(0);
245: }

247: /*@C
248:    TaoComputeJacobian - Computes the Jacobian matrix that has been
249:    set with TaoSetJacobianRoutine().

251:    Collective on Tao

253:    Input Parameters:
254: +  tao - the Tao solver context
255: -  X   - input vector

257:    Output Parameters:
258: +  J    - Jacobian matrix
259: -  Jpre - Preconditioning matrix

261:    Notes:
262:    Most users should not need to explicitly call this routine, as it
263:    is used internally within the minimization solvers.

265:    TaoComputeJacobian() is typically used within minimization
266:    implementations, so most users would not generally call this routine
267:    themselves.

269:    Level: developer

271: .seealso: TaoComputeObjective(), TaoComputeObjectiveAndGradient(), TaoSetJacobianRoutine()
272: @*/
273: PetscErrorCode TaoComputeJacobian(Tao tao, Vec X, Mat J, Mat Jpre)
274: {

281:   if (!tao->ops->computejacobian) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TaoSetJacobian() first");
282:   ++tao->njac;
283:   VecLockReadPush(X);
284:   PetscLogEventBegin(TAO_JacobianEval,tao,X,J,Jpre);
285:   PetscStackPush("Tao user Jacobian function");
286:   (*tao->ops->computejacobian)(tao,X,J,Jpre,tao->user_jacP);
287:   PetscStackPop;
288:   PetscLogEventEnd(TAO_JacobianEval,tao,X,J,Jpre);
289:   VecLockReadPop(X);
290:   return(0);
291: }

293: /*@C
294:    TaoComputeResidualJacobian - Computes the least-squares residual Jacobian matrix that has been
295:    set with TaoSetJacobianResidual().

297:    Collective on Tao

299:    Input Parameters:
300: +  tao - the Tao solver context
301: -  X   - input vector

303:    Output Parameters:
304: +  J    - Jacobian matrix
305: -  Jpre - Preconditioning matrix

307:    Notes:
308:    Most users should not need to explicitly call this routine, as it
309:    is used internally within the minimization solvers.

311:    TaoComputeResidualJacobian() is typically used within least-squares
312:    implementations, so most users would not generally call this routine
313:    themselves.

315:    Level: developer

317: .seealso: TaoComputeResidual(), TaoSetJacobianResidual()
318: @*/
319: PetscErrorCode TaoComputeResidualJacobian(Tao tao, Vec X, Mat J, Mat Jpre)
320: {

327:   if (!tao->ops->computeresidualjacobian) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TaoSetResidualJacobian() first");
328:   ++tao->njac;
329:   VecLockReadPush(X);
330:   PetscLogEventBegin(TAO_JacobianEval,tao,X,J,Jpre);
331:   PetscStackPush("Tao user least-squares residual Jacobian function");
332:   (*tao->ops->computeresidualjacobian)(tao,X,J,Jpre,tao->user_lsjacP);
333:   PetscStackPop;
334:   PetscLogEventEnd(TAO_JacobianEval,tao,X,J,Jpre);
335:   VecLockReadPop(X);
336:   return(0);
337: }

339: /*@C
340:    TaoComputeJacobianState - Computes the Jacobian matrix that has been
341:    set with TaoSetJacobianStateRoutine().

343:    Collective on Tao

345:    Input Parameters:
346: +  tao - the Tao solver context
347: -  X   - input vector

349:    Output Parameters:
350: +  Jpre - Jacobian matrix
351: -  Jinv - Preconditioning matrix

353:    Notes:
354:    Most users should not need to explicitly call this routine, as it
355:    is used internally within the minimization solvers.

357:    TaoComputeJacobianState() is typically used within minimization
358:    implementations, so most users would not generally call this routine
359:    themselves.

361:    Level: developer

363: .seealso: TaoComputeObjective(), TaoComputeObjectiveAndGradient(), TaoSetJacobianStateRoutine(), TaoComputeJacobianDesign(), TaoSetStateDesignIS()
364: @*/
365: PetscErrorCode TaoComputeJacobianState(Tao tao, Vec X, Mat J, Mat Jpre, Mat Jinv)
366: {

373:   if (!tao->ops->computejacobianstate) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TaoSetJacobianState() first");
374:   ++tao->njac_state;
375:   VecLockReadPush(X);
376:   PetscLogEventBegin(TAO_JacobianEval,tao,X,J,Jpre);
377:   PetscStackPush("Tao user Jacobian(state) function");
378:   (*tao->ops->computejacobianstate)(tao,X,J,Jpre,Jinv,tao->user_jac_stateP);
379:   PetscStackPop;
380:   PetscLogEventEnd(TAO_JacobianEval,tao,X,J,Jpre);
381:   VecLockReadPop(X);
382:   return(0);
383: }

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

389:    Collective on Tao

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

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

398:    Notes:
399:    Most users should not need to explicitly call this routine, as it
400:    is used internally within the minimization solvers.

402:    TaoComputeJacobianDesign() is typically used within minimization
403:    implementations, so most users would not generally call this routine
404:    themselves.

406:    Level: developer

408: .seealso:  TaoComputeObjective(), TaoComputeObjectiveAndGradient(), TaoSetJacobianDesignRoutine(), TaoComputeJacobianDesign(), TaoSetStateDesignIS()
409: @*/
410: PetscErrorCode TaoComputeJacobianDesign(Tao tao, Vec X, Mat J)
411: {

418:   if (!tao->ops->computejacobiandesign) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TaoSetJacobianDesign() first");
419:   ++tao->njac_design;
420:   VecLockReadPush(X);
421:   PetscLogEventBegin(TAO_JacobianEval,tao,X,J,NULL);
422:   PetscStackPush("Tao user Jacobian(design) function");
423:   (*tao->ops->computejacobiandesign)(tao,X,J,tao->user_jac_designP);
424:   PetscStackPop;
425:   PetscLogEventEnd(TAO_JacobianEval,tao,X,J,NULL);
426:   VecLockReadPop(X);
427:   return(0);
428: }

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

433:    Logically collective on Tao

435:    Input Parameters:
436: +  tao  - the Tao context
437: .  J    - Matrix used for the jacobian
438: .  Jpre - Matrix that will be used operated on by preconditioner, can be same as J
439: .  func - Jacobian evaluation routine
440: -  ctx  - [optional] user-defined context for private data for the
441:           Jacobian evaluation routine (may be NULL)

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

446: +  tao  - the Tao  context
447: .  x    - input vector
448: .  J    - Jacobian matrix
449: .  Jpre - preconditioning matrix, usually the same as J
450: -  ctx  - [optional] user-defined Jacobian context

452:    Level: intermediate
453: @*/
454: PetscErrorCode TaoSetJacobianRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, void*), void *ctx)
455: {

460:   if (J) {
463:   }
464:   if (Jpre) {
467:   }
468:   if (ctx) {
469:     tao->user_jacP = ctx;
470:   }
471:   if (func) {
472:     tao->ops->computejacobian = func;
473:   }
474:   if (J) {
475:     PetscObjectReference((PetscObject)J);
476:     MatDestroy(&tao->jacobian);
477:     tao->jacobian = J;
478:   }
479:   if (Jpre) {
480:     PetscObjectReference((PetscObject)Jpre);
481:     MatDestroy(&tao->jacobian_pre);
482:     tao->jacobian_pre=Jpre;
483:   }
484:   return(0);
485: }

487: /*@C
488:    TaoSetJacobianResidualRoutine - Sets the function to compute the least-squares residual Jacobian as well as the 
489:    location to store the matrix.

491:    Logically collective on Tao

493:    Input Parameters:
494: +  tao  - the Tao context
495: .  J    - Matrix used for the jacobian
496: .  Jpre - Matrix that will be used operated on by preconditioner, can be same as J
497: .  func - Jacobian evaluation routine
498: -  ctx  - [optional] user-defined context for private data for the
499:           Jacobian evaluation routine (may be NULL)

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

504: +  tao  - the Tao  context
505: .  x    - input vector
506: .  J    - Jacobian matrix
507: .  Jpre - preconditioning matrix, usually the same as J
508: -  ctx  - [optional] user-defined Jacobian context

510:    Level: intermediate
511: @*/
512: PetscErrorCode TaoSetJacobianResidualRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, void*), void *ctx)
513: {

518:   if (J) {
521:   }
522:   if (Jpre) {
525:   }
526:   if (ctx) {
527:     tao->user_lsjacP = ctx;
528:   }
529:   if (func) {
530:     tao->ops->computeresidualjacobian = func;
531:   }
532:   if (J) {
533:     PetscObjectReference((PetscObject)J);
534:     MatDestroy(&tao->ls_jac);
535:     tao->ls_jac = J;
536:   }
537:   if (Jpre) {
538:     PetscObjectReference((PetscObject)Jpre);
539:     MatDestroy(&tao->ls_jac_pre);
540:     tao->ls_jac_pre=Jpre;
541:   }
542:   return(0);
543: }

545: /*@C
546:    TaoSetJacobianStateRoutine - Sets the function to compute the Jacobian
547:    (and its inverse) of the constraint function with respect to the state variables.
548:    Used only for pde-constrained optimization.

550:    Logically collective on Tao

552:    Input Parameters:
553: +  tao  - the Tao context
554: .  J    - Matrix used for the jacobian
555: .  Jpre - Matrix that will be used operated on by PETSc preconditioner, can be same as J.  Only used if Jinv is NULL
556: .  Jinv - [optional] Matrix used to apply the inverse of the state jacobian. Use NULL to default to PETSc KSP solvers to apply the inverse.
557: .  func - Jacobian evaluation routine
558: -  ctx  - [optional] user-defined context for private data for the
559:           Jacobian evaluation routine (may be NULL)

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

564: +  tao  - the Tao  context
565: .  x    - input vector
566: .  J    - Jacobian matrix
567: .  Jpre - preconditioner matrix, usually the same as J
568: .  Jinv - inverse of J
569: -  ctx  - [optional] user-defined Jacobian context

571:    Level: intermediate
572: .seealso: TaoComputeJacobianState(), TaoSetJacobianDesignRoutine(), TaoSetStateDesignIS()
573: @*/
574: PetscErrorCode TaoSetJacobianStateRoutine(Tao tao, Mat J, Mat Jpre, Mat Jinv, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, Mat, void*), void *ctx)
575: {

580:   if (J) {
583:   }
584:   if (Jpre) {
587:   }
588:   if (Jinv) {
591:   }
592:   if (ctx) {
593:     tao->user_jac_stateP = ctx;
594:   }
595:   if (func) {
596:     tao->ops->computejacobianstate = func;
597:   }
598:   if (J) {
599:     PetscObjectReference((PetscObject)J);
600:     MatDestroy(&tao->jacobian_state);
601:     tao->jacobian_state = J;
602:   }
603:   if (Jpre) {
604:     PetscObjectReference((PetscObject)Jpre);
605:     MatDestroy(&tao->jacobian_state_pre);
606:     tao->jacobian_state_pre=Jpre;
607:   }
608:   if (Jinv) {
609:     PetscObjectReference((PetscObject)Jinv);
610:     MatDestroy(&tao->jacobian_state_inv);
611:     tao->jacobian_state_inv=Jinv;
612:   }
613:   return(0);
614: }

616: /*@C
617:    TaoSetJacobianDesignRoutine - Sets the function to compute the Jacobian of
618:    the constraint function with respect to the design variables.  Used only for
619:    pde-constrained optimization.

621:    Logically collective on Tao

623:    Input Parameters:
624: +  tao  - the Tao context
625: .  J    - Matrix used for the jacobian
626: .  func - Jacobian evaluation routine
627: -  ctx  - [optional] user-defined context for private data for the
628:           Jacobian evaluation routine (may be NULL)

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

633: +  tao - the Tao  context
634: .  x   - input vector
635: .  J   - Jacobian matrix
636: -  ctx - [optional] user-defined Jacobian context

638:    Level: intermediate

640: .seealso: TaoComputeJacobianDesign(), TaoSetJacobianStateRoutine(), TaoSetStateDesignIS()
641: @*/
642: PetscErrorCode TaoSetJacobianDesignRoutine(Tao tao, Mat J, PetscErrorCode (*func)(Tao, Vec, Mat, void*), void *ctx)
643: {

648:   if (J) {
651:   }
652:   if (ctx) {
653:     tao->user_jac_designP = ctx;
654:   }
655:   if (func) {
656:     tao->ops->computejacobiandesign = func;
657:   }
658:   if (J) {
659:     PetscObjectReference((PetscObject)J);
660:     MatDestroy(&tao->jacobian_design);
661:     tao->jacobian_design = J;
662:   }
663:   return(0);
664: }

666: /*@
667:    TaoSetStateDesignIS - Indicate to the Tao which variables in the
668:    solution vector are state variables and which are design.  Only applies to
669:    pde-constrained optimization.

671:    Logically Collective on Tao

673:    Input Parameters:
674: +  tao  - The Tao context
675: .  s_is - the index set corresponding to the state variables
676: -  d_is - the index set corresponding to the design variables

678:    Level: intermediate

680: .seealso: TaoSetJacobianStateRoutine(), TaoSetJacobianDesignRoutine()
681: @*/
682: PetscErrorCode TaoSetStateDesignIS(Tao tao, IS s_is, IS d_is)
683: {

687:   PetscObjectReference((PetscObject)s_is);
688:   ISDestroy(&tao->state_is);
689:   tao->state_is = s_is;
690:   PetscObjectReference((PetscObject)(d_is));
691:   ISDestroy(&tao->design_is);
692:   tao->design_is = d_is;
693:   return(0);
694: }

696: /*@C
697:    TaoComputeJacobianEquality - Computes the Jacobian matrix that has been
698:    set with TaoSetJacobianEqualityRoutine().

700:    Collective on Tao

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

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

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

714:    Level: developer

716: .seealso:  TaoComputeObjective(), TaoComputeObjectiveAndGradient(), TaoSetJacobianStateRoutine(), TaoComputeJacobianDesign(), TaoSetStateDesignIS()
717: @*/
718: PetscErrorCode TaoComputeJacobianEquality(Tao tao, Vec X, Mat J, Mat Jpre)
719: {

726:   if (!tao->ops->computejacobianequality) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TaoSetJacobianEquality() first");
727:   ++tao->njac_equality;
728:   VecLockReadPush(X);
729:   PetscLogEventBegin(TAO_JacobianEval,tao,X,J,Jpre);
730:   PetscStackPush("Tao user Jacobian(equality) function");
731:   (*tao->ops->computejacobianequality)(tao,X,J,Jpre,tao->user_jac_equalityP);
732:   PetscStackPop;
733:   PetscLogEventEnd(TAO_JacobianEval,tao,X,J,Jpre);
734:   VecLockReadPop(X);
735:   return(0);
736: }

738: /*@C
739:    TaoComputeJacobianInequality - Computes the Jacobian matrix that has been
740:    set with TaoSetJacobianInequalityRoutine().

742:    Collective on Tao

744:    Input Parameters:
745: +  tao - the Tao solver context
746: -  X   - input vector

748:    Output Parameters:
749: +  J    - Jacobian matrix
750: -  Jpre - Preconditioning matrix

752:    Notes:
753:    Most users should not need to explicitly call this routine, as it
754:    is used internally within the minimization solvers.

756:    Level: developer

758: .seealso:  TaoComputeObjective(), TaoComputeObjectiveAndGradient(), TaoSetJacobianStateRoutine(), TaoComputeJacobianDesign(), TaoSetStateDesignIS()
759: @*/
760: PetscErrorCode TaoComputeJacobianInequality(Tao tao, Vec X, Mat J, Mat Jpre)
761: {

768:   if (!tao->ops->computejacobianinequality) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TaoSetJacobianInequality() first");
769:   ++tao->njac_inequality;
770:   VecLockReadPush(X);
771:   PetscLogEventBegin(TAO_JacobianEval,tao,X,J,Jpre);
772:   PetscStackPush("Tao user Jacobian(inequality) function");
773:   (*tao->ops->computejacobianinequality)(tao,X,J,Jpre,tao->user_jac_inequalityP);
774:   PetscStackPop;
775:   PetscLogEventEnd(TAO_JacobianEval,tao,X,J,Jpre);
776:   VecLockReadPop(X);
777:   return(0);
778: }

780: /*@C
781:    TaoSetJacobianEqualityRoutine - Sets the function to compute the Jacobian
782:    (and its inverse) of the constraint function with respect to the equality variables.
783:    Used only for pde-constrained optimization.

785:    Logically collective on Tao

787:    Input Parameters:
788: +  tao  - the Tao context
789: .  J    - Matrix used for the jacobian
790: .  Jpre - Matrix that will be used operated on by PETSc preconditioner, can be same as J.
791: .  func - Jacobian evaluation routine
792: -  ctx  - [optional] user-defined context for private data for the
793:           Jacobian evaluation routine (may be NULL)

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

798: +  tao  - the Tao  context
799: .  x    - input vector
800: .  J    - Jacobian matrix
801: .  Jpre - preconditioner matrix, usually the same as J
802: -  ctx  - [optional] user-defined Jacobian context

804:    Level: intermediate

806: .seealso: TaoComputeJacobianEquality(), TaoSetJacobianDesignRoutine(), TaoSetEqualityDesignIS()
807: @*/
808: PetscErrorCode TaoSetJacobianEqualityRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat,void*), void *ctx)
809: {

814:   if (J) {
817:   }
818:   if (Jpre) {
821:   }
822:   if (ctx) {
823:     tao->user_jac_equalityP = ctx;
824:   }
825:   if (func) {
826:     tao->ops->computejacobianequality = func;
827:   }
828:   if (J) {
829:     PetscObjectReference((PetscObject)J);
830:     MatDestroy(&tao->jacobian_equality);
831:     tao->jacobian_equality = J;
832:   }
833:   if (Jpre) {
834:     PetscObjectReference((PetscObject)Jpre);
835:     MatDestroy(&tao->jacobian_equality_pre);
836:     tao->jacobian_equality_pre=Jpre;
837:   }
838:   return(0);
839: }

841: /*@C
842:    TaoSetJacobianInequalityRoutine - Sets the function to compute the Jacobian
843:    (and its inverse) of the constraint function with respect to the inequality variables.
844:    Used only for pde-constrained optimization.

846:    Logically collective on Tao

848:    Input Parameters:
849: +  tao  - the Tao context
850: .  J    - Matrix used for the jacobian
851: .  Jpre - Matrix that will be used operated on by PETSc preconditioner, can be same as J.
852: .  func - Jacobian evaluation routine
853: -  ctx  - [optional] user-defined context for private data for the
854:           Jacobian evaluation routine (may be NULL)

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

859: +  tao  - the Tao  context
860: .  x    - input vector
861: .  J    - Jacobian matrix
862: .  Jpre - preconditioner matrix, usually the same as J
863: -  ctx  - [optional] user-defined Jacobian context

865:    Level: intermediate

867: .seealso: TaoComputeJacobianInequality(), TaoSetJacobianDesignRoutine(), TaoSetInequalityDesignIS()
868: @*/
869: PetscErrorCode TaoSetJacobianInequalityRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat,void*), void *ctx)
870: {

875:   if (J) {
878:   }
879:   if (Jpre) {
882:   }
883:   if (ctx) {
884:     tao->user_jac_inequalityP = ctx;
885:   }
886:   if (func) {
887:     tao->ops->computejacobianinequality = func;
888:   }
889:   if (J) {
890:     PetscObjectReference((PetscObject)J);
891:     MatDestroy(&tao->jacobian_inequality);
892:     tao->jacobian_inequality = J;
893:   }
894:   if (Jpre) {
895:     PetscObjectReference((PetscObject)Jpre);
896:     MatDestroy(&tao->jacobian_inequality_pre);
897:     tao->jacobian_inequality_pre=Jpre;
898:   }
899:   return(0);
900: }