Actual source code: taosolver_hj.c
1: #include "include/private/taosolver_impl.h" /*I "taosolver.h" I*/
5: /*@C
6: TaoSetHessianRoutine - Sets the function to compute the Hessian as well as the location to store the matrix.
8: Logically collective on TaoSolver
10: Input Parameters:
11: + tao - the TaoSolver context
12: . H - Matrix used for the hessian
13: . Hpre - Matrix that will be used operated on by preconditioner, can be same as H
14: . hess - Hessian evaluation routine
15: - ctx - [optional] user-defined context for private data for the
16: Hessian evaluation routine (may be PETSC_NULL)
18: Calling sequence of hess:
19: $ hess (TaoSolver tao,Vec x,Mat *H,Mat *Hpre,MatStructure *flag,void *ctx);
21: + tao - the TaoSolver context
22: . x - input vector
23: . H - Hessian matrix
24: . Hpre - preconditioner matrix, usually the same as H
25: . flag - flag indicating information about the preconditioner matrix
26: structure (see below)
27: - ctx - [optional] user-defined Hessian context
30: Notes:
32: The function hess() takes Mat * as the matrix arguments rather than Mat.
33: This allows the Hessian evaluation routine to replace A and/or B with a
34: completely new new matrix structure (not just different matrix elements)
35: when appropriate, for instance, if the nonzero structure is changing
36: throughout the global iterations.
38: The flag can be used to eliminate unnecessary work in the preconditioner
39: during the repeated solution of linear systems of the same size. The
40: available options are
41: $ SAME_PRECONDITIONER -
42: $ Hpre is identical during successive linear solves.
43: $ This option is intended for folks who are using
44: $ different Amat and Pmat matrices and want to reuse the
45: $ same preconditioner matrix. For example, this option
46: $ saves work by not recomputing incomplete factorization
47: $ for ILU/ICC preconditioners.
48: $ SAME_NONZERO_PATTERN -
49: $ Hpre has the same nonzero structure during
50: $ successive linear solves.
51: $ DIFFERENT_NONZERO_PATTERN -
52: $ Hpre does not have the same nonzero structure.
54: Caution:
55: If you specify SAME_NONZERO_PATTERN, the software believes your assertion
56: and does not check the structure of the matrix. If you erroneously
57: claim that the structure is the same when it actually is not, the new
58: preconditioner will not function correctly. Thus, use this optimization
59: feature carefully!
61: If in doubt about whether your preconditioner matrix has changed
62: structure or not, use the flag DIFFERENT_NONZERO_PATTERN.
64: Level: beginner
66: @*/
67: PetscErrorCode TaoSetHessianRoutine(TaoSolver tao, Mat H, Mat Hpre, PetscErrorCode (*func)(TaoSolver, Vec, Mat*, Mat *, MatStructure *, void*), void *ctx)
68: {
72: if (H) {
75: }
76: if (Hpre) {
79: }
80: if (ctx) {
81: tao->user_hessP = ctx;
82: }
83: if (func) {
84: tao->ops->computehessian = func;
85: }
87:
88: if (H) {
89: PetscObjectReference((PetscObject)H);
90: if (tao->hessian) { MatDestroy(&tao->hessian); }
91: tao->hessian = H;
92: }
93: if (Hpre) {
94: PetscObjectReference((PetscObject)Hpre);
95: if (tao->hessian_pre) { MatDestroy(&tao->hessian_pre); }
96: tao->hessian_pre=Hpre;
97: }
98: return(0);
99: }
103: /*@C
104: TaoComputeHessian - Computes the Hessian matrix that has been
105: set with TaoSetHessianRoutine().
107: Collective on TaoSolver
109: Input Parameters:
110: + solver - the TaoSolver solver context
111: - xx - input vector
113: Output Parameters:
114: + H - Hessian matrix
115: . Hpre - Preconditioning matrix
116: - flag - flag indicating matrix structure (SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN, or SAME_PRECONDITIONER)
118: Notes:
119: Most users should not need to explicitly call this routine, as it
120: is used internally within the minimization solvers.
122: TaoComputeHessian() is typically used within minimization
123: implementations, so most users would not generally call this routine
124: themselves.
126: Level: developer
128: .seealso: TaoComputeObjective(), TaoComputeObjectiveAndGradient(), TaoSetHessian()
130: @*/
131: PetscErrorCode TaoComputeHessian(TaoSolver tao, Vec X, Mat *H, Mat *Hpre, MatStructure *flg)
132: {
139:
140: if (!tao->ops->computehessian) {
141: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TaoSetHessian() first");
142: }
143: *flg = DIFFERENT_NONZERO_PATTERN;
144: ++tao->nhess;
145: PetscLogEventBegin(TaoSolver_HessianEval,tao,X,*H,*Hpre);
146: PetscStackPush("TaoSolver user Hessian function");
147: CHKMEMQ;
148: (*tao->ops->computehessian)(tao,X,H,Hpre,flg,tao->user_hessP);
149: CHKMEMQ;
150: PetscStackPop;
151: PetscLogEventEnd(TaoSolver_HessianEval,tao,X,*H,*Hpre);
152:
153: return(0);
154: }
158: /*@C
159: TaoComputeJacobian - Computes the Jacobian matrix that has been
160: set with TaoSetJacobianRoutine().
162: Collective on TaoSolver
164: Input Parameters:
165: + solver - the TaoSolver solver context
166: - xx - input vector
168: Output Parameters:
169: + H - Jacobian matrix
170: . Hpre - Preconditioning matrix
171: - flag - flag indicating matrix structure (SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN, or SAME_PRECONDITIONER)
173: Notes:
174: Most users should not need to explicitly call this routine, as it
175: is used internally within the minimization solvers.
177: TaoComputeJacobian() is typically used within minimization
178: implementations, so most users would not generally call this routine
179: themselves.
181: Level: developer
183: .seealso: TaoComputeObjective(), TaoComputeObjectiveAndGradient(), TaoSetJacobian()
185: @*/
186: PetscErrorCode TaoComputeJacobian(TaoSolver tao, Vec X, Mat *J, Mat *Jpre, MatStructure *flg)
187: {
194:
195: if (!tao->ops->computejacobian) {
196: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TaoSetJacobian() first");
197: }
198: *flg = DIFFERENT_NONZERO_PATTERN;
199: ++tao->njac;
200: PetscLogEventBegin(TaoSolver_JacobianEval,tao,X,*J,*Jpre);
201: PetscStackPush("TaoSolver user Jacobian function");
202: CHKMEMQ;
203: (*tao->ops->computejacobian)(tao,X,J,Jpre,flg,tao->user_jacP);
204: CHKMEMQ;
205: PetscStackPop;
206: PetscLogEventEnd(TaoSolver_JacobianEval,tao,X,*J,*Jpre);
207:
208: return(0);
209: }
213: /*@C
214: TaoComputeJacobianState - Computes the Jacobian matrix that has been
215: set with TaoSetJacobianStateRoutine().
217: Collective on TaoSolver
219: Input Parameters:
220: + solver - the TaoSolver solver context
221: - xx - input vector
223: Output Parameters:
224: + H - Jacobian matrix
225: . Hpre - Preconditioning matrix
226: - flag - flag indicating matrix structure (SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN, or SAME_PRECONDITIONER)
228: Notes:
229: Most users should not need to explicitly call this routine, as it
230: is used internally within the minimization solvers.
232: TaoComputeJacobianState() is typically used within minimization
233: implementations, so most users would not generally call this routine
234: themselves.
236: Level: developer
238: .seealso: TaoComputeObjective(), TaoComputeObjectiveAndGradient(), TaoSetJacobianStateRoutine(), TaoComputeJacobianDesign(), TaoSetStateDesignIS()
240: @*/
241: PetscErrorCode TaoComputeJacobianState(TaoSolver tao, Vec X, Mat *J, Mat *Jpre, Mat *Jinv, MatStructure *flg)
242: {
249:
250: if (!tao->ops->computejacobianstate) {
251: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TaoSetJacobianState() first");
252: }
253: *flg = DIFFERENT_NONZERO_PATTERN;
254: ++tao->njac_state;
255: PetscLogEventBegin(TaoSolver_JacobianEval,tao,X,*J,*Jpre);
256: PetscStackPush("TaoSolver user Jacobian(state) function");
257: CHKMEMQ;
258: (*tao->ops->computejacobianstate)(tao,X,J,Jpre,Jinv,flg,tao->user_jac_stateP);
259: CHKMEMQ;
260: PetscStackPop;
261: PetscLogEventEnd(TaoSolver_JacobianEval,tao,X,*J,*Jpre);
262:
263: return(0);
264: }
268: /*@C
269: TaoComputeJacobianDesign - Computes the Jacobian matrix that has been
270: set with TaoSetJacobianDesignRoutine().
272: Collective on TaoSolver
274: Input Parameters:
275: + solver - the TaoSolver solver context
276: - xx - input vector
278: Output Parameters:
279: . H - Jacobian matrix
281: Notes:
282: Most users should not need to explicitly call this routine, as it
283: is used internally within the minimization solvers.
285: TaoComputeJacobianDesign() is typically used within minimization
286: implementations, so most users would not generally call this routine
287: themselves.
289: Level: developer
291: .seealso: TaoComputeObjective(), TaoComputeObjectiveAndGradient(), TaoSetJacobianDesignRoutine(), TaoComputeJacobianDesign(), TaoSetStateDesignIS()
293: @*/
294: PetscErrorCode TaoComputeJacobianDesign(TaoSolver tao, Vec X, Mat *J)
295: {
301:
302: if (!tao->ops->computejacobiandesign) {
303: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TaoSetJacobianDesign() first");
304: }
305: ++tao->njac_design;
306: PetscLogEventBegin(TaoSolver_JacobianEval,tao,X,*J,PETSC_NULL);
307: PetscStackPush("TaoSolver user Jacobian(design) function");
308: CHKMEMQ;
309: (*tao->ops->computejacobiandesign)(tao,X,J,tao->user_jac_designP);
310: CHKMEMQ;
311: PetscStackPop;
312: PetscLogEventEnd(TaoSolver_JacobianEval,tao,X,*J,PETSC_NULL);
313:
314: return(0);
315: }
321: /*@C
322: TaoSetJacobianRoutine - Sets the function to compute the Jacobian as well as the location to store the matrix.
324: Logically collective on TaoSolver
326: Input Parameters:
327: + tao - the TaoSolver context
328: . J - Matrix used for the jacobian
329: . Jpre - Matrix that will be used operated on by preconditioner, can be same as J
330: . jac - Jacobian evaluation routine
331: - ctx - [optional] user-defined context for private data for the
332: Jacobian evaluation routine (may be PETSC_NULL)
334: Calling sequence of jac:
335: $ jac (TaoSolver tao,Vec x,Mat *J,Mat *Jpre,MatStructure *flag,void *ctx);
337: + tao - the TaoSolver context
338: . x - input vector
339: . J - Jacobian matrix
340: . Jpre - preconditioner matrix, usually the same as J
341: . flag - flag indicating information about the preconditioner matrix
342: structure (see below)
343: - ctx - [optional] user-defined Jacobian context
346: Notes:
348: The function jac() takes Mat * as the matrix arguments rather than Mat.
349: This allows the Jacobian evaluation routine to replace A and/or B with a
350: completely new new matrix structure (not just different matrix elements)
351: when appropriate, for instance, if the nonzero structure is changing
352: throughout the global iterations.
354: The flag can be used to eliminate unnecessary work in the preconditioner
355: during the repeated solution of linear systems of the same size. The
356: available options are
357: $ SAME_PRECONDITIONER -
358: $ Jpre is identical during successive linear solves.
359: $ This option is intended for folks who are using
360: $ different Amat and Pmat matrices and want to reuse the
361: $ same preconditioner matrix. For example, this option
362: $ saves work by not recomputing incomplete factorization
363: $ for ILU/ICC preconditioners.
364: $ SAME_NONZERO_PATTERN -
365: $ Jpre has the same nonzero structure during
366: $ successive linear solves.
367: $ DIFFERENT_NONZERO_PATTERN -
368: $ Jpre does not have the same nonzero structure.
370: Caution:
371: If you specify SAME_NONZERO_PATTERN, the software believes your assertion
372: and does not check the structure of the matrix. If you erroneously
373: claim that the structure is the same when it actually is not, the new
374: preconditioner will not function correctly. Thus, use this optimization
375: feature carefully!
377: If in doubt about whether your preconditioner matrix has changed
378: structure or not, use the flag DIFFERENT_NONZERO_PATTERN.
380: Level: intermediate
382: @*/
383: PetscErrorCode TaoSetJacobianRoutine(TaoSolver tao, Mat J, Mat Jpre, PetscErrorCode (*func)(TaoSolver, Vec, Mat*, Mat *, MatStructure *, void*), void *ctx)
384: {
388: if (J) {
391: }
392: if (Jpre) {
395: }
396: if (ctx) {
397: tao->user_jacP = ctx;
398: }
399: if (func) {
400: tao->ops->computejacobian = func;
401: }
403:
404: if (J) {
405: PetscObjectReference((PetscObject)J);
406: if (tao->jacobian) { MatDestroy(&tao->jacobian); }
407: tao->jacobian = J;
408: }
409: if (Jpre) {
410: PetscObjectReference((PetscObject)Jpre);
411: if (tao->jacobian_pre) { MatDestroy(&tao->jacobian_pre); }
412: tao->jacobian_pre=Jpre;
413: }
414: return(0);
415: }
420: /*@C
421: TaoSetJacobianStateRoutine - Sets the function to compute the Jacobian
422: (and its inverse) of the constraint function with respect to the state variables.
423: Used only for pde-constrained optimization.
425: Logically collective on TaoSolver
427: Input Parameters:
428: + tao - the TaoSolver context
429: . J - Matrix used for the jacobian
430: . Jpre - Matrix that will be used operated on by PETSc preconditioner, can be same as J. Only used if Jinv is PETSC_NULL
431: . Jinv - [optional] Matrix used to apply the inverse of the state jacobian. Use PETSC_NULL to default to PETSc KSP solvers to apply the inverse.
432: . jac - Jacobian evaluation routine
433: - ctx - [optional] user-defined context for private data for the
434: Jacobian evaluation routine (may be PETSC_NULL)
436: Calling sequence of jac:
437: $ jac (TaoSolver tao,Vec x,Mat *J,Mat *Jpre,MatStructure *flag,void *ctx);
439: + tao - the TaoSolver context
440: . x - input vector
441: . J - Jacobian matrix
442: . Jpre - preconditioner matrix, usually the same as J
443: . Jinv - inverse of J
444: . flag - flag indicating information about the preconditioner matrix
445: structure (see below)
446: - ctx - [optional] user-defined Jacobian context
449: Notes:
450: Because of the structure of the jacobian matrix,
451: It may be more efficient for a pde-constrained application to provide
452: its own Jinv matrix.
454: The function jac() takes Mat * as the matrix arguments rather than Mat.
455: This allows the Jacobian evaluation routine to replace A and/or B with a
456: completely new new maitrix structure (not just different matrix elements)
457: when appropriate, for instance, if the nonzero structure is changing
458: throughout the global iterations.
460: The flag can be used to eliminate unnecessary work in the preconditioner
461: during the repeated solution of linear systems of the same size. The
462: available options are
463: $ SAME_PRECONDITIONER -
464: $ Jpre is identical during successive linear solves.
465: $ This option is intended for folks who are using
466: $ different Amat and Pmat matrices and want to reuse the
467: $ same preconditioner matrix. For example, this option
468: $ saves work by not recomputing incomplete factorization
469: $ for ILU/ICC preconditioners.
470: $ SAME_NONZERO_PATTERN -
471: $ Jpre has the same nonzero structure during
472: $ successive linear solves.
473: $ DIFFERENT_NONZERO_PATTERN -
474: $ Jpre does not have the same nonzero structure.
476: Caution:
477: If you specify SAME_NONZERO_PATTERN, the software believes your assertion
478: and does not check the structure of the matrix. If you erroneously
479: claim that the structure is the same when it actually is not, the new
480: preconditioner will not function correctly. Thus, use this optimization
481: feature carefully!
483: If in doubt about whether your preconditioner matrix has changed
484: structure or not, use the flag DIFFERENT_NONZERO_PATTERN.
486: Level: intermediate
487: .seealse: TaoComputeJacobianState(), TaoSetJacobianDesignRoutine(), TaoSetStateDesignIS()
488: @*/
489: PetscErrorCode TaoSetJacobianStateRoutine(TaoSolver tao, Mat J, Mat Jpre, Mat Jinv, PetscErrorCode (*func)(TaoSolver, Vec, Mat*, Mat *, Mat *, MatStructure *, void*), void *ctx)
490: {
494: if (J) {
497: }
498: if (Jpre) {
501: }
502: if (Jinv) {
505: }
506: if (ctx) {
507: tao->user_jac_stateP = ctx;
508: }
509: if (func) {
510: tao->ops->computejacobianstate = func;
511: }
513:
514: if (J) {
515: PetscObjectReference((PetscObject)J);
516: if (tao->jacobian_state) { MatDestroy(&tao->jacobian_state); }
517: tao->jacobian_state = J;
518: }
519: if (Jpre) {
520: PetscObjectReference((PetscObject)Jpre);
521: if (tao->jacobian_state_pre) { MatDestroy(&tao->jacobian_state_pre); }
522: tao->jacobian_state_pre=Jpre;
523: }
524: if (Jinv) {
525: PetscObjectReference((PetscObject)Jinv);
526: if (tao->jacobian_state_inv) {MatDestroy(&tao->jacobian_state_inv); }
527: tao->jacobian_state_inv=Jinv;
528: }
529: return(0);
530: }
534: /*@C
535: TaoSetJacobianDesignRoutine - Sets the function to compute the Jacobian of
536: the constraint function with respect to the design variables. Used only for
537: pde-constrained optimization.
539: Logically collective on TaoSolver
541: Input Parameters:
542: + tao - the TaoSolver context
543: . J - Matrix used for the jacobian
544: . jac - Jacobian evaluation routine
545: - ctx - [optional] user-defined context for private data for the
546: Jacobian evaluation routine (may be PETSC_NULL)
548: Calling sequence of jac:
549: $ jac (TaoSolver tao,Vec x,Mat *J,void *ctx);
551: + tao - the TaoSolver context
552: . x - input vector
553: . J - Jacobian matrix
554: - ctx - [optional] user-defined Jacobian context
557: Notes:
559: The function jac() takes Mat * as the matrix arguments rather than Mat.
560: This allows the Jacobian evaluation routine to replace A and/or B with a
561: completely new new matrix structure (not just different matrix elements)
562: when appropriate, for instance, if the nonzero structure is changing
563: throughout the global iterations.
565: Level: intermediate
566: .seealso: TaoComputeJacobianDesign(), TaoSetJacobianStateRoutine(), TaoSetStateDesignIS()
567: @*/
568: PetscErrorCode TaoSetJacobianDesignRoutine(TaoSolver tao, Mat J, PetscErrorCode (*func)(TaoSolver, Vec, Mat*, void*), void *ctx)
569: {
573: if (J) {
576: }
577: if (ctx) {
578: tao->user_jac_designP = ctx;
579: }
580: if (func) {
581: tao->ops->computejacobiandesign = func;
582: }
584:
585: if (J) {
586: PetscObjectReference((PetscObject)J);
587: if (tao->jacobian_design) { MatDestroy(&tao->jacobian_design); }
588: tao->jacobian_design = J;
589: }
590: return(0);
591: }
595: /*@
596: TaoSetStateDesignIS - Indicate to the TaoSolver which variables in the
597: solution vector are state variables and which are design. Only applies to
598: pde-constrained optimization.
600: Logically Collective on TaoSolver
602: Input Parameters:
603: + tao - The TaoSolver context
604: . s_is - the index set corresponding to the state variables
605: - d_is - the index set corresponding to the design variables
606:
607: Level: intermediate
609: .seealso: TaoSetJacobianStateRoutine(), TaoSetJacobianDesignRoutine()
610: @*/
611: PetscErrorCode TaoSetStateDesignIS(TaoSolver tao, IS s_is, IS d_is)
612: {
614: if (tao->state_is) {
615: PetscObjectDereference((PetscObject)(tao->state_is));
616: }
617: if (tao->design_is) {
618: PetscObjectDereference((PetscObject)(tao->design_is));
619: }
620: tao->state_is = s_is;
621: tao->design_is = d_is;
622: if (s_is) {
623: PetscObjectReference((PetscObject)(tao->state_is));
624: }
625: if (d_is) {
626: PetscObjectReference((PetscObject)(tao->design_is));
627: }
628: return(0);
629: }