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 on Tao
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: 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: 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: PetscErrorCode ierr;
97: PetscReal nrm,gnorm;
98: PetscReal threshold = 1.e-5;
99: PetscInt m,n,M,N;
100: PetscBool complete_print = PETSC_FALSE,test = PETSC_FALSE,flg;
101: PetscViewer viewer,mviewer;
102: MPI_Comm comm;
103: PetscInt tabs;
104: static PetscBool directionsprinted = PETSC_FALSE;
105: PetscViewerFormat format;
107: PetscObjectOptionsBegin((PetscObject)tao);
108: PetscOptionsName("-tao_test_hessian","Compare hand-coded and finite difference Hessians","None",&test);
109: PetscOptionsReal("-tao_test_hessian", "Threshold for element difference between hand-coded and finite difference being meaningful","None",threshold,&threshold,NULL);
110: PetscOptionsViewer("-tao_test_hessian_view","View difference between hand-coded and finite difference Hessians element entries","None",&mviewer,&format,&complete_print);
111: PetscOptionsEnd();
112: if (!test) return 0;
114: PetscObjectGetComm((PetscObject)tao,&comm);
115: PetscViewerASCIIGetStdout(comm,&viewer);
116: PetscViewerASCIIGetTab(viewer, &tabs);
117: PetscViewerASCIISetTab(viewer, ((PetscObject)tao)->tablevel);
118: PetscViewerASCIIPrintf(viewer," ---------- Testing Hessian -------------\n");
119: if (!complete_print && !directionsprinted) {
120: PetscViewerASCIIPrintf(viewer," Run with -tao_test_hessian_view and optionally -tao_test_hessian <threshold> to show difference\n");
121: PetscViewerASCIIPrintf(viewer," of hand-coded and finite difference Hessian entries greater than <threshold>.\n");
122: }
123: if (!directionsprinted) {
124: PetscViewerASCIIPrintf(viewer," Testing hand-coded Hessian, if (for double precision runs) ||J - Jfd||_F/||J||_F is\n");
125: PetscViewerASCIIPrintf(viewer," O(1.e-8), the hand-coded Hessian is probably correct.\n");
126: directionsprinted = PETSC_TRUE;
127: }
128: if (complete_print) {
129: PetscViewerPushFormat(mviewer,format);
130: }
132: PetscObjectTypeCompare((PetscObject)tao->hessian,MATMFFD,&flg);
133: if (!flg) hessian = tao->hessian;
134: else hessian = tao->hessian_pre;
136: while (hessian) {
137: PetscObjectBaseTypeCompareAny((PetscObject)hessian,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");
138: if (flg) {
139: A = hessian;
140: PetscObjectReference((PetscObject)A);
141: } else {
142: MatComputeOperator(hessian,MATAIJ,&A);
143: }
145: MatCreate(PetscObjectComm((PetscObject)A),&B);
146: MatGetSize(A,&M,&N);
147: MatGetLocalSize(A,&m,&n);
148: MatSetSizes(B,m,n,M,N);
149: MatSetType(B,((PetscObject)A)->type_name);
150: MatSetUp(B);
151: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
153: TaoDefaultComputeHessian(tao,x,B,B,NULL);
155: MatDuplicate(B,MAT_COPY_VALUES,&D);
156: MatAYPX(D,-1.0,A,DIFFERENT_NONZERO_PATTERN);
157: MatNorm(D,NORM_FROBENIUS,&nrm);
158: MatNorm(A,NORM_FROBENIUS,&gnorm);
159: MatDestroy(&D);
160: if (!gnorm) gnorm = 1; /* just in case */
161: PetscViewerASCIIPrintf(viewer," ||H - Hfd||_F/||H||_F = %g, ||H - Hfd||_F = %g\n",(double)(nrm/gnorm),(double)nrm);
163: if (complete_print) {
164: PetscViewerASCIIPrintf(viewer," Hand-coded Hessian ----------\n");
165: MatView(A,mviewer);
166: PetscViewerASCIIPrintf(viewer," Finite difference Hessian ----------\n");
167: MatView(B,mviewer);
168: }
170: if (complete_print) {
171: PetscInt Istart, Iend, *ccols, bncols, cncols, j, row;
172: PetscScalar *cvals;
173: const PetscInt *bcols;
174: const PetscScalar *bvals;
176: MatAYPX(B,-1.0,A,DIFFERENT_NONZERO_PATTERN);
177: MatCreate(PetscObjectComm((PetscObject)A),&C);
178: MatSetSizes(C,m,n,M,N);
179: MatSetType(C,((PetscObject)A)->type_name);
180: MatSetUp(C);
181: MatSetOption(C,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
182: MatGetOwnershipRange(B,&Istart,&Iend);
184: for (row = Istart; row < Iend; row++) {
185: MatGetRow(B,row,&bncols,&bcols,&bvals);
186: PetscMalloc2(bncols,&ccols,bncols,&cvals);
187: for (j = 0, cncols = 0; j < bncols; j++) {
188: if (PetscAbsScalar(bvals[j]) > threshold) {
189: ccols[cncols] = bcols[j];
190: cvals[cncols] = bvals[j];
191: cncols += 1;
192: }
193: }
194: if (cncols) {
195: MatSetValues(C,1,&row,cncols,ccols,cvals,INSERT_VALUES);
196: }
197: MatRestoreRow(B,row,&bncols,&bcols,&bvals);
198: PetscFree2(ccols,cvals);
199: }
200: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
201: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
202: PetscViewerASCIIPrintf(viewer," Finite-difference minus hand-coded Hessian with tolerance %g ----------\n",(double)threshold);
203: MatView(C,mviewer);
204: MatDestroy(&C);
205: }
206: MatDestroy(&A);
207: MatDestroy(&B);
209: if (hessian != tao->hessian_pre) {
210: hessian = tao->hessian_pre;
211: PetscViewerASCIIPrintf(viewer," ---------- Testing Hessian for preconditioner -------------\n");
212: } else hessian = NULL;
213: }
214: if (complete_print) {
215: PetscViewerPopFormat(mviewer);
216: PetscViewerDestroy(&mviewer);
217: }
218: PetscViewerASCIISetTab(viewer,tabs);
219: return 0;
220: }
222: /*@C
223: TaoComputeHessian - Computes the Hessian matrix that has been
224: set with TaoSetHessian().
226: Collective on Tao
228: Input Parameters:
229: + tao - the Tao solver context
230: - X - input vector
232: Output Parameters:
233: + H - Hessian matrix
234: - Hpre - Preconditioning matrix
236: Options Database Keys:
237: + -tao_test_hessian - compare the user provided Hessian with one compute via finite differences to check for errors
238: . -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
239: - -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
241: Notes:
242: Most users should not need to explicitly call this routine, as it
243: is used internally within the minimization solvers.
245: TaoComputeHessian() is typically used within minimization
246: implementations, so most users would not generally call this routine
247: themselves.
249: Developer Notes:
250: The Hessian test mechanism follows SNESTestJacobian().
252: Level: developer
254: .seealso: TaoComputeObjective(), TaoComputeObjectiveAndGradient(), TaoSetHessian()
255: @*/
256: PetscErrorCode TaoComputeHessian(Tao tao, Vec X, Mat H, Mat Hpre)
257: {
262: ++tao->nhess;
263: VecLockReadPush(X);
264: PetscLogEventBegin(TAO_HessianEval,tao,X,H,Hpre);
265: PetscStackPush("Tao user Hessian function");
266: (*tao->ops->computehessian)(tao,X,H,Hpre,tao->user_hessP);
267: PetscStackPop;
268: PetscLogEventEnd(TAO_HessianEval,tao,X,H,Hpre);
269: VecLockReadPop(X);
271: TaoTestHessian(tao);
272: return 0;
273: }
275: /*@C
276: TaoComputeJacobian - Computes the Jacobian matrix that has been
277: set with TaoSetJacobianRoutine().
279: Collective on Tao
281: Input Parameters:
282: + tao - the Tao solver context
283: - X - input vector
285: Output Parameters:
286: + J - Jacobian matrix
287: - Jpre - Preconditioning matrix
289: Notes:
290: Most users should not need to explicitly call this routine, as it
291: is used internally within the minimization solvers.
293: TaoComputeJacobian() is typically used within minimization
294: implementations, so most users would not generally call this routine
295: themselves.
297: Level: developer
299: .seealso: TaoComputeObjective(), TaoComputeObjectiveAndGradient(), TaoSetJacobianRoutine()
300: @*/
301: PetscErrorCode TaoComputeJacobian(Tao tao, Vec X, Mat J, Mat Jpre)
302: {
307: ++tao->njac;
308: VecLockReadPush(X);
309: PetscLogEventBegin(TAO_JacobianEval,tao,X,J,Jpre);
310: PetscStackPush("Tao user Jacobian function");
311: (*tao->ops->computejacobian)(tao,X,J,Jpre,tao->user_jacP);
312: PetscStackPop;
313: PetscLogEventEnd(TAO_JacobianEval,tao,X,J,Jpre);
314: VecLockReadPop(X);
315: return 0;
316: }
318: /*@C
319: TaoComputeResidualJacobian - Computes the least-squares residual Jacobian matrix that has been
320: set with TaoSetJacobianResidual().
322: Collective on Tao
324: Input Parameters:
325: + tao - the Tao solver context
326: - X - input vector
328: Output Parameters:
329: + J - Jacobian matrix
330: - Jpre - Preconditioning matrix
332: Notes:
333: Most users should not need to explicitly call this routine, as it
334: is used internally within the minimization solvers.
336: TaoComputeResidualJacobian() is typically used within least-squares
337: implementations, so most users would not generally call this routine
338: themselves.
340: Level: developer
342: .seealso: TaoComputeResidual(), TaoSetJacobianResidual()
343: @*/
344: PetscErrorCode TaoComputeResidualJacobian(Tao tao, Vec X, Mat J, Mat Jpre)
345: {
350: ++tao->njac;
351: VecLockReadPush(X);
352: PetscLogEventBegin(TAO_JacobianEval,tao,X,J,Jpre);
353: PetscStackPush("Tao user least-squares residual Jacobian function");
354: (*tao->ops->computeresidualjacobian)(tao,X,J,Jpre,tao->user_lsjacP);
355: PetscStackPop;
356: PetscLogEventEnd(TAO_JacobianEval,tao,X,J,Jpre);
357: VecLockReadPop(X);
358: return 0;
359: }
361: /*@C
362: TaoComputeJacobianState - Computes the Jacobian matrix that has been
363: set with TaoSetJacobianStateRoutine().
365: Collective on Tao
367: Input Parameters:
368: + tao - the Tao solver context
369: - X - input vector
371: Output Parameters:
372: + J - Jacobian matrix
373: . Jpre - Preconditioning matrix
374: - Jinv -
376: Notes:
377: Most users should not need to explicitly call this routine, as it
378: is used internally within the minimization solvers.
380: TaoComputeJacobianState() is typically used within minimization
381: implementations, so most users would not generally call this routine
382: themselves.
384: Level: developer
386: .seealso: TaoComputeObjective(), TaoComputeObjectiveAndGradient(), TaoSetJacobianStateRoutine(), TaoComputeJacobianDesign(), TaoSetStateDesignIS()
387: @*/
388: PetscErrorCode TaoComputeJacobianState(Tao tao, Vec X, Mat J, Mat Jpre, Mat Jinv)
389: {
394: ++tao->njac_state;
395: VecLockReadPush(X);
396: PetscLogEventBegin(TAO_JacobianEval,tao,X,J,Jpre);
397: PetscStackPush("Tao user Jacobian(state) function");
398: (*tao->ops->computejacobianstate)(tao,X,J,Jpre,Jinv,tao->user_jac_stateP);
399: PetscStackPop;
400: PetscLogEventEnd(TAO_JacobianEval,tao,X,J,Jpre);
401: VecLockReadPop(X);
402: return 0;
403: }
405: /*@C
406: TaoComputeJacobianDesign - Computes the Jacobian matrix that has been
407: set with TaoSetJacobianDesignRoutine().
409: Collective on Tao
411: Input Parameters:
412: + tao - the Tao solver context
413: - X - input vector
415: Output Parameters:
416: . J - Jacobian matrix
418: Notes:
419: Most users should not need to explicitly call this routine, as it
420: is used internally within the minimization solvers.
422: TaoComputeJacobianDesign() is typically used within minimization
423: implementations, so most users would not generally call this routine
424: themselves.
426: Level: developer
428: .seealso: TaoComputeObjective(), TaoComputeObjectiveAndGradient(), TaoSetJacobianDesignRoutine(), TaoComputeJacobianDesign(), TaoSetStateDesignIS()
429: @*/
430: PetscErrorCode TaoComputeJacobianDesign(Tao tao, Vec X, Mat J)
431: {
436: ++tao->njac_design;
437: VecLockReadPush(X);
438: PetscLogEventBegin(TAO_JacobianEval,tao,X,J,NULL);
439: PetscStackPush("Tao user Jacobian(design) function");
440: (*tao->ops->computejacobiandesign)(tao,X,J,tao->user_jac_designP);
441: PetscStackPop;
442: PetscLogEventEnd(TAO_JacobianEval,tao,X,J,NULL);
443: VecLockReadPop(X);
444: return 0;
445: }
447: /*@C
448: TaoSetJacobianRoutine - Sets the function to compute the Jacobian as well as the location to store the matrix.
450: Logically collective on Tao
452: Input Parameters:
453: + tao - the Tao context
454: . J - Matrix used for the jacobian
455: . Jpre - Matrix that will be used operated on by preconditioner, can be same as J
456: . func - Jacobian evaluation routine
457: - ctx - [optional] user-defined context for private data for the
458: Jacobian evaluation routine (may be NULL)
460: Calling sequence of func:
461: $ func(Tao tao,Vec x,Mat J,Mat Jpre,void *ctx);
463: + tao - the Tao context
464: . x - input vector
465: . J - Jacobian matrix
466: . Jpre - preconditioning matrix, usually the same as J
467: - ctx - [optional] user-defined Jacobian context
469: Level: intermediate
470: @*/
471: PetscErrorCode TaoSetJacobianRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, void*), void *ctx)
472: {
474: if (J) {
477: }
478: if (Jpre) {
481: }
482: if (ctx) {
483: tao->user_jacP = ctx;
484: }
485: if (func) {
486: tao->ops->computejacobian = func;
487: }
488: if (J) {
489: PetscObjectReference((PetscObject)J);
490: MatDestroy(&tao->jacobian);
491: tao->jacobian = J;
492: }
493: if (Jpre) {
494: PetscObjectReference((PetscObject)Jpre);
495: MatDestroy(&tao->jacobian_pre);
496: tao->jacobian_pre=Jpre;
497: }
498: return 0;
499: }
501: /*@C
502: TaoSetJacobianResidualRoutine - Sets the function to compute the least-squares residual Jacobian as well as the
503: location to store the matrix.
505: Logically collective on Tao
507: Input Parameters:
508: + tao - the Tao context
509: . J - Matrix used for the jacobian
510: . Jpre - Matrix that will be used operated on by preconditioner, can be same as J
511: . func - Jacobian evaluation routine
512: - ctx - [optional] user-defined context for private data for the
513: Jacobian evaluation routine (may be NULL)
515: Calling sequence of func:
516: $ func(Tao tao,Vec x,Mat J,Mat Jpre,void *ctx);
518: + tao - the Tao context
519: . x - input vector
520: . J - Jacobian matrix
521: . Jpre - preconditioning matrix, usually the same as J
522: - ctx - [optional] user-defined Jacobian context
524: Level: intermediate
525: @*/
526: PetscErrorCode TaoSetJacobianResidualRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, void*), void *ctx)
527: {
529: if (J) {
532: }
533: if (Jpre) {
536: }
537: if (ctx) {
538: tao->user_lsjacP = ctx;
539: }
540: if (func) {
541: tao->ops->computeresidualjacobian = func;
542: }
543: if (J) {
544: PetscObjectReference((PetscObject)J);
545: MatDestroy(&tao->ls_jac);
546: tao->ls_jac = J;
547: }
548: if (Jpre) {
549: PetscObjectReference((PetscObject)Jpre);
550: MatDestroy(&tao->ls_jac_pre);
551: tao->ls_jac_pre=Jpre;
552: }
553: return 0;
554: }
556: /*@C
557: TaoSetJacobianStateRoutine - Sets the function to compute the Jacobian
558: (and its inverse) of the constraint function with respect to the state variables.
559: Used only for pde-constrained optimization.
561: Logically collective on Tao
563: Input Parameters:
564: + tao - the Tao context
565: . J - Matrix used for the jacobian
566: . Jpre - Matrix that will be used operated on by PETSc preconditioner, can be same as J. Only used if Jinv is NULL
567: . Jinv - [optional] Matrix used to apply the inverse of the state jacobian. Use NULL to default to PETSc KSP solvers to apply the inverse.
568: . func - Jacobian evaluation routine
569: - ctx - [optional] user-defined context for private data for the
570: Jacobian evaluation routine (may be NULL)
572: Calling sequence of func:
573: $ func(Tao tao,Vec x,Mat J,Mat Jpre,Mat Jinv,void *ctx);
575: + tao - the Tao context
576: . x - input vector
577: . J - Jacobian matrix
578: . Jpre - preconditioner matrix, usually the same as J
579: . Jinv - inverse of J
580: - ctx - [optional] user-defined Jacobian context
582: Level: intermediate
583: .seealso: TaoComputeJacobianState(), TaoSetJacobianDesignRoutine(), TaoSetStateDesignIS()
584: @*/
585: PetscErrorCode TaoSetJacobianStateRoutine(Tao tao, Mat J, Mat Jpre, Mat Jinv, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, Mat, void*), void *ctx)
586: {
588: if (J) {
591: }
592: if (Jpre) {
595: }
596: if (Jinv) {
599: }
600: if (ctx) {
601: tao->user_jac_stateP = ctx;
602: }
603: if (func) {
604: tao->ops->computejacobianstate = func;
605: }
606: if (J) {
607: PetscObjectReference((PetscObject)J);
608: MatDestroy(&tao->jacobian_state);
609: tao->jacobian_state = J;
610: }
611: if (Jpre) {
612: PetscObjectReference((PetscObject)Jpre);
613: MatDestroy(&tao->jacobian_state_pre);
614: tao->jacobian_state_pre=Jpre;
615: }
616: if (Jinv) {
617: PetscObjectReference((PetscObject)Jinv);
618: MatDestroy(&tao->jacobian_state_inv);
619: tao->jacobian_state_inv=Jinv;
620: }
621: return 0;
622: }
624: /*@C
625: TaoSetJacobianDesignRoutine - Sets the function to compute the Jacobian of
626: the constraint function with respect to the design variables. Used only for
627: pde-constrained optimization.
629: Logically collective on Tao
631: Input Parameters:
632: + tao - the Tao context
633: . J - Matrix used for the jacobian
634: . func - Jacobian evaluation routine
635: - ctx - [optional] user-defined context for private data for the
636: Jacobian evaluation routine (may be NULL)
638: Calling sequence of func:
639: $ func(Tao tao,Vec x,Mat J,void *ctx);
641: + tao - the Tao context
642: . x - input vector
643: . J - Jacobian matrix
644: - ctx - [optional] user-defined Jacobian context
646: Level: intermediate
648: .seealso: TaoComputeJacobianDesign(), TaoSetJacobianStateRoutine(), TaoSetStateDesignIS()
649: @*/
650: PetscErrorCode TaoSetJacobianDesignRoutine(Tao tao, Mat J, PetscErrorCode (*func)(Tao, Vec, Mat, void*), void *ctx)
651: {
653: if (J) {
656: }
657: if (ctx) {
658: tao->user_jac_designP = ctx;
659: }
660: if (func) {
661: tao->ops->computejacobiandesign = func;
662: }
663: if (J) {
664: PetscObjectReference((PetscObject)J);
665: MatDestroy(&tao->jacobian_design);
666: tao->jacobian_design = J;
667: }
668: return 0;
669: }
671: /*@
672: TaoSetStateDesignIS - Indicate to the Tao which variables in the
673: solution vector are state variables and which are design. Only applies to
674: pde-constrained optimization.
676: Logically Collective on Tao
678: Input Parameters:
679: + tao - The Tao context
680: . s_is - the index set corresponding to the state variables
681: - d_is - the index set corresponding to the design variables
683: Level: intermediate
685: .seealso: TaoSetJacobianStateRoutine(), TaoSetJacobianDesignRoutine()
686: @*/
687: PetscErrorCode TaoSetStateDesignIS(Tao tao, IS s_is, IS d_is)
688: {
689: PetscObjectReference((PetscObject)s_is);
690: ISDestroy(&tao->state_is);
691: tao->state_is = s_is;
692: PetscObjectReference((PetscObject)(d_is));
693: ISDestroy(&tao->design_is);
694: tao->design_is = d_is;
695: return 0;
696: }
698: /*@C
699: TaoComputeJacobianEquality - Computes the Jacobian matrix that has been
700: set with TaoSetJacobianEqualityRoutine().
702: Collective on Tao
704: Input Parameters:
705: + tao - the Tao solver context
706: - X - input vector
708: Output Parameters:
709: + J - Jacobian matrix
710: - Jpre - Preconditioning matrix
712: Notes:
713: Most users should not need to explicitly call this routine, as it
714: is used internally within the minimization solvers.
716: Level: developer
718: .seealso: TaoComputeObjective(), TaoComputeObjectiveAndGradient(), TaoSetJacobianStateRoutine(), TaoComputeJacobianDesign(), TaoSetStateDesignIS()
719: @*/
720: PetscErrorCode TaoComputeJacobianEquality(Tao tao, Vec X, Mat J, Mat Jpre)
721: {
726: ++tao->njac_equality;
727: VecLockReadPush(X);
728: PetscLogEventBegin(TAO_JacobianEval,tao,X,J,Jpre);
729: PetscStackPush("Tao user Jacobian(equality) function");
730: (*tao->ops->computejacobianequality)(tao,X,J,Jpre,tao->user_jac_equalityP);
731: PetscStackPop;
732: PetscLogEventEnd(TAO_JacobianEval,tao,X,J,Jpre);
733: VecLockReadPop(X);
734: return 0;
735: }
737: /*@C
738: TaoComputeJacobianInequality - Computes the Jacobian matrix that has been
739: set with TaoSetJacobianInequalityRoutine().
741: Collective on Tao
743: Input Parameters:
744: + tao - the Tao solver context
745: - X - input vector
747: Output Parameters:
748: + J - Jacobian matrix
749: - Jpre - Preconditioning matrix
751: Notes:
752: Most users should not need to explicitly call this routine, as it
753: is used internally within the minimization solvers.
755: Level: developer
757: .seealso: TaoComputeObjective(), TaoComputeObjectiveAndGradient(), TaoSetJacobianStateRoutine(), TaoComputeJacobianDesign(), TaoSetStateDesignIS()
758: @*/
759: PetscErrorCode TaoComputeJacobianInequality(Tao tao, Vec X, Mat J, Mat Jpre)
760: {
765: ++tao->njac_inequality;
766: VecLockReadPush(X);
767: PetscLogEventBegin(TAO_JacobianEval,tao,X,J,Jpre);
768: PetscStackPush("Tao user Jacobian(inequality) function");
769: (*tao->ops->computejacobianinequality)(tao,X,J,Jpre,tao->user_jac_inequalityP);
770: PetscStackPop;
771: PetscLogEventEnd(TAO_JacobianEval,tao,X,J,Jpre);
772: VecLockReadPop(X);
773: return 0;
774: }
776: /*@C
777: TaoSetJacobianEqualityRoutine - Sets the function to compute the Jacobian
778: (and its inverse) of the constraint function with respect to the equality variables.
779: Used only for pde-constrained optimization.
781: Logically collective on Tao
783: Input Parameters:
784: + tao - the Tao context
785: . J - Matrix used for the jacobian
786: . Jpre - Matrix that will be used operated on by PETSc preconditioner, can be same as J.
787: . func - Jacobian evaluation routine
788: - ctx - [optional] user-defined context for private data for the
789: Jacobian evaluation routine (may be NULL)
791: Calling sequence of func:
792: $ func(Tao tao,Vec x,Mat J,Mat Jpre,void *ctx);
794: + tao - the Tao context
795: . x - input vector
796: . J - Jacobian matrix
797: . Jpre - preconditioner matrix, usually the same as J
798: - ctx - [optional] user-defined Jacobian context
800: Level: intermediate
802: .seealso: TaoComputeJacobianEquality(), TaoSetJacobianDesignRoutine(), TaoSetEqualityDesignIS()
803: @*/
804: PetscErrorCode TaoSetJacobianEqualityRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, void*), void *ctx)
805: {
807: if (J) {
810: }
811: if (Jpre) {
814: }
815: if (ctx) {
816: tao->user_jac_equalityP = ctx;
817: }
818: if (func) {
819: tao->ops->computejacobianequality = func;
820: }
821: if (J) {
822: PetscObjectReference((PetscObject)J);
823: MatDestroy(&tao->jacobian_equality);
824: tao->jacobian_equality = J;
825: }
826: if (Jpre) {
827: PetscObjectReference((PetscObject)Jpre);
828: MatDestroy(&tao->jacobian_equality_pre);
829: tao->jacobian_equality_pre=Jpre;
830: }
831: return 0;
832: }
834: /*@C
835: TaoSetJacobianInequalityRoutine - Sets the function to compute the Jacobian
836: (and its inverse) of the constraint function with respect to the inequality variables.
837: Used only for pde-constrained optimization.
839: Logically collective on Tao
841: Input Parameters:
842: + tao - the Tao context
843: . J - Matrix used for the jacobian
844: . Jpre - Matrix that will be used operated on by PETSc preconditioner, can be same as J.
845: . func - Jacobian evaluation routine
846: - ctx - [optional] user-defined context for private data for the
847: Jacobian evaluation routine (may be NULL)
849: Calling sequence of func:
850: $ func(Tao tao,Vec x,Mat J,Mat Jpre,void *ctx);
852: + tao - the Tao context
853: . x - input vector
854: . J - Jacobian matrix
855: . Jpre - preconditioner matrix, usually the same as J
856: - ctx - [optional] user-defined Jacobian context
858: Level: intermediate
860: .seealso: TaoComputeJacobianInequality(), TaoSetJacobianDesignRoutine(), TaoSetInequalityDesignIS()
861: @*/
862: PetscErrorCode TaoSetJacobianInequalityRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat,void*), void *ctx)
863: {
865: if (J) {
868: }
869: if (Jpre) {
872: }
873: if (ctx) {
874: tao->user_jac_inequalityP = ctx;
875: }
876: if (func) {
877: tao->ops->computejacobianinequality = func;
878: }
879: if (J) {
880: PetscObjectReference((PetscObject)J);
881: MatDestroy(&tao->jacobian_inequality);
882: tao->jacobian_inequality = J;
883: }
884: if (Jpre) {
885: PetscObjectReference((PetscObject)Jpre);
886: MatDestroy(&tao->jacobian_inequality_pre);
887: tao->jacobian_inequality_pre=Jpre;
888: }
889: return 0;
890: }