Actual source code: taosolver_hj.c

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