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: }