Actual source code: taoimpl.h

      1: #pragma once

      3: #include <petsctao.h>
      4: #include <petsctaolinesearch.h>
      5: #include <petsc/private/petscimpl.h>

      7: PETSC_EXTERN PetscBool      TaoRegisterAllCalled;
      8: PETSC_EXTERN PetscErrorCode TaoRegisterAll(void);

     10: typedef struct _TaoOps *TaoOps;

     12: struct _TaoOps {
     13:   /* Methods set by application */
     14:   PetscErrorCode (*computeresidual)(Tao, Vec, Vec, void *);
     15:   PetscErrorCode (*computeresidualjacobian)(Tao, Vec, Mat, Mat, void *);
     16:   PetscErrorCode (*computeconstraints)(Tao, Vec, Vec, void *);
     17:   PetscErrorCode (*computeinequalityconstraints)(Tao, Vec, Vec, void *);
     18:   PetscErrorCode (*computeequalityconstraints)(Tao, Vec, Vec, void *);
     19:   PetscErrorCode (*computejacobian)(Tao, Vec, Mat, Mat, void *);
     20:   PetscErrorCode (*computejacobianstate)(Tao, Vec, Mat, Mat, Mat, void *);
     21:   PetscErrorCode (*computejacobiandesign)(Tao, Vec, Mat, void *);
     22:   PetscErrorCode (*computejacobianinequality)(Tao, Vec, Mat, Mat, void *);
     23:   PetscErrorCode (*computejacobianequality)(Tao, Vec, Mat, Mat, void *);
     24:   PetscErrorCode (*computebounds)(Tao, Vec, Vec, void *);
     25:   PetscErrorCode (*update)(Tao, PetscInt, void *);
     26:   PetscErrorCode (*convergencetest)(Tao, void *);
     27:   PetscErrorCode (*convergencedestroy)(void *);

     29:   /* Methods set by solver */
     30:   PetscErrorCode (*computedual)(Tao, Vec, Vec);
     31:   PetscErrorCode (*setup)(Tao);
     32:   PetscErrorCode (*solve)(Tao);
     33:   PetscErrorCode (*view)(Tao, PetscViewer);
     34:   PetscErrorCode (*setfromoptions)(Tao, PetscOptionItems);
     35:   PetscErrorCode (*destroy)(Tao);
     36: };

     38: #define MAXTAOMONITORS 10

     40: typedef struct _n_TaoTermMapping TaoTermMapping;

     42: /*S
     43:    TaoTermMapping - Object held by either `Tao` or `TaoTerm` with `TAOTERMSUM` type
     44:    that contain necessary information regarding mapping matrix

     46:    Level: developer

     48:    Notes:
     49:    This object is internally handled in `TaoAddTerm()`, which couples `TaoTerm` with
     50:    appropriately mapped gradients, mapped Hessians, and other necessary information.

     52:    This object is also used in `TAOTERMSUM` to handle mapping of summands of `TAOTERMSUM`.

     54:    Users would not need to work directly with the `TaoTermMapping`, rather, they work with
     55:    `Tao` with added terms via `TaoAddTerm()`, or with `TaoTerm` with type `TAOTERMSUM`.

     57:    Developer Note:
     58:    Currently, as `MatPtAP` does not support diagonal matrices, internal work matrices was added
     59:    for a workaround.

     61: .seealso: [](ch_tao), `Tao`, `TaoAddTerm()`, `TAOTERMSUM`,
     62: S*/
     63: struct _n_TaoTermMapping {
     64:   char       *prefix;
     65:   TaoTerm     term;
     66:   PetscReal   scale;
     67:   Mat         map;
     68:   Vec         _map_output;
     69:   Vec         _unmapped_gradient;
     70:   Vec         _mapped_gradient;
     71:   Mat         _unmapped_H;
     72:   Mat         _unmapped_Hpre;
     73:   Mat         _mapped_H;
     74:   Mat         _mapped_Hpre;
     75:   Mat         _mapped_H_work; /* Temporary work matrices for PtAP for diagonal A */
     76:   Mat         _mapped_Hpre_work;
     77:   TaoTermMask mask;
     78: };

     80: #define TaoTermObjectiveMasked(a) ((a) & TAOTERM_MASK_OBJECTIVE)
     81: #define TaoTermGradientMasked(a)  ((a) & TAOTERM_MASK_GRADIENT)
     82: #define TaoTermHessianMasked(a)   ((a) & TAOTERM_MASK_HESSIAN)

     84: struct _p_Tao {
     85:   PETSCHEADER(struct _TaoOps);
     86:   PetscCtx ctx; /* user provided context */
     87:   void    *user_lsresP;
     88:   void    *user_lsjacP;
     89:   void    *user_conP;
     90:   void    *user_con_equalityP;
     91:   void    *user_con_inequalityP;
     92:   void    *user_jacP;
     93:   void    *user_jac_equalityP;
     94:   void    *user_jac_inequalityP;
     95:   void    *user_jac_stateP;
     96:   void    *user_jac_designP;
     97:   void    *user_boundsP;
     98:   void    *user_update;

    100:   PetscErrorCode (*monitor[MAXTAOMONITORS])(Tao, void *);
    101:   PetscCtxDestroyFn *monitordestroy[MAXTAOMONITORS];
    102:   void              *monitorcontext[MAXTAOMONITORS];
    103:   PetscInt           numbermonitors;
    104:   void              *cnvP;
    105:   TaoConvergedReason reason;

    107:   PetscBool setupcalled;
    108:   void     *data;

    110:   Vec        solution;
    111:   Vec        gradient;
    112:   Vec        stepdirection;
    113:   Vec        XL;
    114:   Vec        XU;
    115:   Vec        IL;
    116:   Vec        IU;
    117:   Vec        DI;
    118:   Vec        DE;
    119:   Mat        hessian;
    120:   Mat        hessian_pre;
    121:   Mat        gradient_norm;
    122:   Vec        gradient_norm_tmp;
    123:   Vec        ls_res;
    124:   Mat        ls_jac;
    125:   Mat        ls_jac_pre;
    126:   Vec        res_weights_v;
    127:   PetscInt   res_weights_n;
    128:   PetscInt  *res_weights_rows;
    129:   PetscInt  *res_weights_cols;
    130:   PetscReal *res_weights_w;
    131:   Vec        constraints;
    132:   Vec        constraints_equality;
    133:   Vec        constraints_inequality;
    134:   Mat        jacobian;
    135:   Mat        jacobian_pre;
    136:   Mat        jacobian_inequality;
    137:   Mat        jacobian_inequality_pre;
    138:   Mat        jacobian_equality;
    139:   Mat        jacobian_equality_pre;
    140:   Mat        jacobian_state;
    141:   Mat        jacobian_state_inv;
    142:   Mat        jacobian_design;
    143:   Mat        jacobian_state_pre;
    144:   Mat        jacobian_design_pre;
    145:   IS         state_is;
    146:   IS         design_is;
    147:   PetscReal  step;
    148:   PetscReal  residual;
    149:   PetscReal  gnorm0;
    150:   PetscReal  cnorm;
    151:   PetscReal  cnorm0;
    152:   PetscReal  fc;

    154:   PetscInt max_constraints;
    155:   PetscInt niter;
    156:   PetscInt ntotalits;
    157:   PetscInt nconstraints;
    158:   PetscInt niconstraints;
    159:   PetscInt neconstraints;
    160:   PetscInt nres;
    161:   PetscInt njac;
    162:   PetscInt njac_equality;
    163:   PetscInt njac_inequality;
    164:   PetscInt njac_state;
    165:   PetscInt njac_design;

    167:   PetscInt ksp_its;     /* KSP iterations for this solver iteration */
    168:   PetscInt ksp_tot_its; /* Total (cumulative) KSP iterations */

    170:   TaoLineSearch linesearch;
    171:   PetscBool     lsflag; /* goes up when line search fails */
    172:   KSP           ksp;
    173:   PetscReal     trust; /* Current trust region */

    175:   /* EW type forcing term */
    176:   PetscBool ksp_ewconv;
    177:   SNES      snes_ewdummy;

    179:   PetscObjectParameterDeclare(PetscReal, gatol);
    180:   PetscObjectParameterDeclare(PetscReal, grtol);
    181:   PetscObjectParameterDeclare(PetscReal, gttol);
    182:   PetscObjectParameterDeclare(PetscReal, catol);
    183:   PetscObjectParameterDeclare(PetscReal, crtol);
    184:   PetscObjectParameterDeclare(PetscReal, steptol);
    185:   PetscObjectParameterDeclare(PetscReal, fmin);
    186:   PetscObjectParameterDeclare(PetscInt, max_it);
    187:   PetscObjectParameterDeclare(PetscInt, max_funcs);
    188:   PetscObjectParameterDeclare(PetscReal, trust0); /* initial trust region radius */

    190:   PetscBool printreason;
    191:   PetscBool viewsolution;
    192:   PetscBool viewgradient;
    193:   PetscBool viewconstraints;
    194:   PetscBool viewhessian;
    195:   PetscBool viewjacobian;
    196:   PetscBool bounded;
    197:   PetscBool constrained;
    198:   PetscBool eq_constrained;
    199:   PetscBool ineq_constrained;
    200:   PetscBool ineq_doublesided;
    201:   PetscBool header_printed;
    202:   PetscBool recycle;

    204:   TaoSubsetType subset_type;
    205:   PetscInt      hist_max;   /* Number of iteration histories to keep */
    206:   PetscReal    *hist_obj;   /* obj value at each iteration */
    207:   PetscReal    *hist_resid; /* residual at each iteration */
    208:   PetscReal    *hist_cnorm; /* constraint norm at each iteration */
    209:   PetscInt     *hist_lits;  /* number of ksp its at each TAO iteration */
    210:   PetscInt      hist_len;
    211:   PetscBool     hist_reset;
    212:   PetscBool     hist_malloc;

    214:   TaoTermMapping objective_term; /* TaoTerm in use */
    215:   Vec            objective_parameters;
    216:   PetscInt       num_terms;
    217:   PetscBool      term_set;

    219:   TaoTerm   callbacks; /* TAOTERMCALLBACKS for the original callbacks */
    220:   PetscBool uses_hessian_matrices;
    221:   PetscBool uses_gradient;
    222: };

    224: PETSC_EXTERN PetscLogEvent TAO_Solve;
    225: PETSC_EXTERN PetscLogEvent TAO_ConstraintsEval;
    226: PETSC_EXTERN PetscLogEvent TAO_JacobianEval;
    227: PETSC_INTERN PetscLogEvent TAO_ResidualEval;

    229: PETSC_INTERN PetscLogEvent TAOTERM_ObjectiveEval;
    230: PETSC_INTERN PetscLogEvent TAOTERM_GradientEval;
    231: PETSC_INTERN PetscLogEvent TAOTERM_ObjGradEval;
    232: PETSC_INTERN PetscLogEvent TAOTERM_HessianEval;

    234: static inline PetscErrorCode TaoLogConvergenceHistory(Tao tao, PetscReal obj, PetscReal resid, PetscReal cnorm, PetscInt totits)
    235: {
    236:   PetscFunctionBegin;
    237:   if (tao->hist_max > tao->hist_len) {
    238:     if (tao->hist_obj) tao->hist_obj[tao->hist_len] = obj;
    239:     if (tao->hist_resid) tao->hist_resid[tao->hist_len] = resid;
    240:     if (tao->hist_cnorm) tao->hist_cnorm[tao->hist_len] = cnorm;
    241:     if (tao->hist_lits) {
    242:       PetscInt sits = totits;
    243:       PetscCheck(tao->hist_len >= 0, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "History length cannot be negative");
    244:       for (PetscInt i = 0; i < tao->hist_len; i++) sits -= tao->hist_lits[i];
    245:       tao->hist_lits[tao->hist_len] = sits;
    246:     }
    247:     tao->hist_len++;
    248:   }
    249:   PetscFunctionReturn(PETSC_SUCCESS);
    250: }

    252: PETSC_INTERN PetscErrorCode TaoTestGradient_Internal(Tao, Vec, Vec, PetscViewer, PetscViewer);

    254: typedef struct _TaoTermOps *TaoTermOps;

    256: struct _TaoTermOps {
    257:   PetscErrorCode (*setfromoptions)(TaoTerm, PetscOptionItems);
    258:   PetscErrorCode (*setup)(TaoTerm);
    259:   PetscErrorCode (*view)(TaoTerm, PetscViewer);
    260:   PetscErrorCode (*destroy)(TaoTerm);

    262:   TaoTermObjectiveFn            *objective;
    263:   TaoTermObjectiveAndGradientFn *objectiveandgradient;
    264:   TaoTermGradientFn             *gradient;
    265:   TaoTermHessianFn              *hessian;

    267:   PetscErrorCode (*isobjectivedefined)(TaoTerm, PetscBool *);
    268:   PetscErrorCode (*isgradientdefined)(TaoTerm, PetscBool *);
    269:   PetscErrorCode (*isobjectiveandgradientdefined)(TaoTerm, PetscBool *);
    270:   PetscErrorCode (*ishessiandefined)(TaoTerm, PetscBool *);
    271:   PetscErrorCode (*iscreatehessianmatricesdefined)(TaoTerm, PetscBool *);
    272:   PetscErrorCode (*iscomputehessianfdpossible)(TaoTerm, PetscBool3 *);

    274:   PetscErrorCode (*createsolutionvec)(TaoTerm, Vec *);
    275:   PetscErrorCode (*createparametersvec)(TaoTerm, Vec *);
    276:   PetscErrorCode (*createhessianmatrices)(TaoTerm, Mat *, Mat *);
    277: };

    279: struct _p_TaoTerm {
    280:   PETSCHEADER(struct _TaoTermOps);
    281:   void                 *data;
    282:   PetscBool             setup_called;
    283:   Mat                   solution_factory; // dummies used to create vectors
    284:   Mat                   parameters_factory;
    285:   Mat                   parameters_factory_orig; // copy so that parameters_factory can be made a reference of solution_factory if parameter space == vector space
    286:   TaoTermParametersMode parameters_mode;
    287:   PetscBool             Hpre_is_H; // Hessian mode data
    288:   MatType               H_mattype;
    289:   MatType               Hpre_mattype;

    291:   PetscInt ngrad_mffd;
    292:   PetscInt nobj;     // actual objective callback invocations
    293:   PetscInt ngrad;    // actual gradient callback invocations
    294:   PetscInt nobjgrad; // actual objective+gradient callback invocations
    295:   PetscInt nhess;    // actual Hessian callback invocations

    297:   PetscReal fd_delta;    // for TaoTermComputeGradientFD()
    298:   PetscBool fd_gradient; // use finite differences for the gradient
    299:   PetscBool fd_hessian;  // use finite differences for the Hessian
    300: };

    302: PETSC_INTERN PetscErrorCode TaoTermRegisterAll(void);

    304: PETSC_INTERN PetscErrorCode TaoTermCreateCallbacks(Tao, TaoTerm *);

    306: PETSC_INTERN PetscErrorCode TaoTermCreate_ElementwiseDivergence_Internal(TaoTerm);
    307: PETSC_INTERN PetscErrorCode TaoTermDestroy_ElementwiseDivergence_Internal(TaoTerm);

    309: PETSC_INTERN PetscErrorCode TaoTermCallbacksSetObjective(TaoTerm, PetscErrorCode (*)(Tao, Vec, PetscReal *, PetscCtx), PetscCtx);
    310: PETSC_INTERN PetscErrorCode TaoTermCallbacksSetGradient(TaoTerm, PetscErrorCode (*)(Tao, Vec, Vec, PetscCtx), PetscCtx);
    311: PETSC_INTERN PetscErrorCode TaoTermCallbacksSetObjectiveAndGradient(TaoTerm, PetscErrorCode (*)(Tao, Vec, PetscReal *, Vec, PetscCtx), PetscCtx);
    312: PETSC_INTERN PetscErrorCode TaoTermCallbacksSetHessian(TaoTerm, PetscErrorCode (*)(Tao, Vec, Mat, Mat, PetscCtx), PetscCtx);

    314: PETSC_INTERN PetscErrorCode TaoTermCallbacksGetObjective(TaoTerm, PetscErrorCode (**)(Tao, Vec, PetscReal *, PetscCtx), PetscCtxRt);
    315: PETSC_INTERN PetscErrorCode TaoTermCallbacksGetGradient(TaoTerm, PetscErrorCode (**)(Tao, Vec, Vec, PetscCtx), PetscCtxRt);
    316: PETSC_INTERN PetscErrorCode TaoTermCallbacksGetObjectiveAndGradient(TaoTerm, PetscErrorCode (**)(Tao, Vec, PetscReal *, Vec, PetscCtx), PetscCtxRt);
    317: PETSC_INTERN PetscErrorCode TaoTermCallbacksGetHessian(TaoTerm, PetscErrorCode (**)(Tao, Vec, Mat, Mat, PetscCtx), PetscCtxRt);

    319: PETSC_INTERN PetscErrorCode TaoTermMappingSetData(TaoTermMapping *, const char *, PetscReal, TaoTerm, Mat);
    320: PETSC_INTERN PetscErrorCode TaoTermMappingGetData(TaoTermMapping *, const char **, PetscReal *, TaoTerm *, Mat *);
    321: PETSC_INTERN PetscErrorCode TaoTermMappingReset(TaoTermMapping *);
    322: PETSC_INTERN PetscErrorCode TaoTermMappingComputeObjective(TaoTermMapping *, Vec, Vec, InsertMode, PetscReal *);
    323: PETSC_INTERN PetscErrorCode TaoTermMappingComputeGradient(TaoTermMapping *, Vec, Vec, InsertMode, Vec);
    324: PETSC_INTERN PetscErrorCode TaoTermMappingComputeObjectiveAndGradient(TaoTermMapping *, Vec, Vec, InsertMode, PetscReal *, Vec);
    325: PETSC_INTERN PetscErrorCode TaoTermMappingComputeHessian(TaoTermMapping *, Vec, Vec, InsertMode, Mat, Mat);
    326: PETSC_INTERN PetscErrorCode TaoTermMappingSetUp(TaoTermMapping *);
    327: PETSC_INTERN PetscErrorCode TaoTermMappingCreateSolutionVec(TaoTermMapping *, Vec *);
    328: PETSC_INTERN PetscErrorCode TaoTermMappingCreateParametersVec(TaoTermMapping *, Vec *);
    329: PETSC_INTERN PetscErrorCode TaoTermMappingCreateHessianMatrices(TaoTermMapping *, Mat *, Mat *);

    331: PETSC_INTERN PetscErrorCode VecIfNotCongruentGetSameLayoutVec(Vec, Vec *);

    333: PETSC_INTERN PetscErrorCode TaoTermCreateHessianMatricesDefault_H_Internal(TaoTerm, Mat *, Mat *, PetscBool, MatType);
    334: PETSC_INTERN PetscErrorCode TaoTermCreateHessianMatricesDefault_Hpre_Internal(TaoTerm, Mat *, Mat *, PetscBool, MatType);