Actual source code: taosolver_hj.c

petsc-master 2020-07-14
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:     PetscViewerDestroy(&mviewer);
187:   }
188:   PetscViewerASCIISetTab(viewer,tabs);
189:   return(0);
190: }

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

196:    Collective on Tao

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

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

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

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

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

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

222:    Level: developer

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

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

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

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

252:    Collective on Tao

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

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

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

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

270:    Level: developer

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

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

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

298:    Collective on Tao

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

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

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

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

316:    Level: developer

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

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

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

344:    Collective on Tao

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

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

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

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

362:    Level: developer

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

374:   if (!tao->ops->computejacobianstate) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TaoSetJacobianState() first");
375:   ++tao->njac_state;
376:   VecLockReadPush(X);
377:   PetscLogEventBegin(TAO_JacobianEval,tao,X,J,Jpre);
378:   PetscStackPush("Tao user Jacobian(state) function");
379:   (*tao->ops->computejacobianstate)(tao,X,J,Jpre,Jinv,tao->user_jac_stateP);
380:   PetscStackPop;
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 on Tao

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 minimization solvers.

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

407:    Level: developer

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

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

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

434:    Logically collective on Tao

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

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

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

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

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

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

492:    Logically collective on Tao

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

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

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

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

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

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

551:    Logically collective on Tao

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

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

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

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

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

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

622:    Logically collective on Tao

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

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

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

639:    Level: intermediate

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

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

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

672:    Logically Collective on Tao

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

679:    Level: intermediate

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

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

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

701:    Collective on Tao

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:  TaoComputeObjective(), TaoComputeObjectiveAndGradient(), TaoSetJacobianStateRoutine(), TaoComputeJacobianDesign(), TaoSetStateDesignIS()
718: @*/
719: PetscErrorCode TaoComputeJacobianEquality(Tao tao, Vec X, Mat J, Mat Jpre)
720: {

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

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

743:    Collective on Tao

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

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

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

757:    Level: developer

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

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

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

786:    Logically collective on Tao

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

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

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

805:    Level: intermediate

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

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

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

847:    Logically collective on Tao

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

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

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

866:    Level: intermediate

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

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