1:  #include <../src/tao/bound/impls/bqnk/bqnk.h>

  3: static const char *BNK_AS[64] = {"none", "bertsekas"};

  5: static PetscErrorCode TaoBQNLSComputeHessian(Tao tao)
  6: {
  7:   TAO_BNK        *bnk = (TAO_BNK *)tao->data;
  8:   TAO_BQNK       *bqnk = (TAO_BQNK*)bnk->ctx;
 10:   PetscReal      gnorm2, delta;
 11: 
 13:   /* Compute the initial scaling and update the approximation */
 14:   gnorm2 = bnk->gnorm*bnk->gnorm;
 15:   if (gnorm2 == 0.0) gnorm2 = PETSC_MACHINE_EPSILON;
 16:   if (bnk->f == 0.0) {
 17:     delta = 2.0 / gnorm2;
 18:   } else {
 19:     delta = 2.0 * PetscAbsScalar(bnk->f) / gnorm2;
 20:   }
 21:   MatSymBrdnSetDelta(bqnk->B, delta);
 22:   MatLMVMUpdate(bqnk->B, tao->solution, bnk->unprojected_gradient);
 23:   return(0);
 24: }

 26: static PetscErrorCode TaoBQNLSComputeStep(Tao tao, PetscBool shift, KSPConvergedReason *ksp_reason, PetscInt *step_type)
 27: {
 28:   TAO_BNK        *bnk = (TAO_BNK *)tao->data;
 29:   TAO_BQNK       *bqnk = (TAO_BQNK*)bnk->ctx;
 31:   PetscInt       nupdates;

 34:   MatSolve(bqnk->B, tao->gradient, tao->stepdirection);
 35:   VecScale(tao->stepdirection, -1.0);
 36:   TaoBNKBoundStep(tao, bnk->as_type, tao->stepdirection);
 37:   *ksp_reason = KSP_CONVERGED_ATOL;
 38:   MatLMVMGetUpdateCount(bqnk->B, &nupdates);
 39:   if (nupdates == 0) {
 40:     *step_type = BNK_SCALED_GRADIENT;
 41:   } else {
 42:     *step_type = BNK_BFGS;
 43:   }
 44:   return(0);
 45: }

 47: static PetscErrorCode TaoSetFromOptions_BQNLS(PetscOptionItems *PetscOptionsObject,Tao tao)
 48: {
 49:   TAO_BNK        *bnk = (TAO_BNK *)tao->data;
 50:   TAO_BQNK       *bqnk = (TAO_BQNK*)bnk->ctx;
 52:   KSPType        ksp_type;
 53:   PetscBool      is_spd;

 56:   PetscOptionsHead(PetscOptionsObject,"Quasi-Newton-Krylov method for bound constrained optimization");
 57:   PetscOptionsEList("-tao_bqnls_as_type", "active set estimation method", "", BNK_AS, BNK_AS_TYPES, BNK_AS[bnk->as_type], &bnk->as_type, 0);
 58:   PetscOptionsReal("-tao_bqnls_epsilon", "(developer) tolerance used when computing actual and predicted reduction", "", bnk->epsilon, &bnk->epsilon,NULL);
 59:   PetscOptionsReal("-tao_bqnls_as_tol", "(developer) initial tolerance used when estimating actively bounded variables", "", bnk->as_tol, &bnk->as_tol,NULL);
 60:   PetscOptionsReal("-tao_bqnls_as_step", "(developer) step length used when estimating actively bounded variables", "", bnk->as_step, &bnk->as_step,NULL);
 61:   PetscOptionsInt("-tao_bqnls_max_cg_its", "number of BNCG iterations to take for each Newton step", "", bnk->max_cg_its, &bnk->max_cg_its,NULL);
 62:   PetscOptionsTail();
 63:   TaoSetFromOptions(bnk->bncg);
 64:   TaoLineSearchSetFromOptions(tao->linesearch);
 65:   KSPSetFromOptions(tao->ksp);
 66:   KSPGetType(tao->ksp,&ksp_type);
 67:   bnk->is_nash = bnk->is_gltr = bnk->is_stcg = PETSC_FALSE;
 68:   MatSetFromOptions(bqnk->B);
 69:   MatGetOption(bqnk->B, MAT_SPD, &is_spd);
 70:   if (!is_spd) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix must be symmetric positive-definite");
 71:   return(0);
 72: }

 74: /*MC
 75:   TAOBQNLS - Bounded Quasi-Newton Line Search method for nonlinear minimization with bound 
 76:              constraints. This method approximates the action of the inverse-Hessian with a 
 77:              limited memory quasi-Newton formula. The quasi-Newton matrix and its options are 
 78:              accessible via the prefix `-tao_bqnls_`

 80:   Options Database Keys:
 81: + -tao_bqnls_max_cg_its - maximum number of bounded conjugate-gradient iterations taken in each Newton loop
 82: - -tao_bqnls_as_type - active-set estimation method ("none", "bertsekas")

 84:   Level: beginner
 85: M*/
 86: PETSC_EXTERN PetscErrorCode TaoCreate_BQNLS(Tao tao)
 87: {
 88:   TAO_BNK        *bnk;
 89:   TAO_BQNK       *bqnk;
 91: 
 93:   TaoCreate_BQNK(tao);
 94:   KSPSetOptionsPrefix(tao->ksp, "unused");
 95:   tao->ops->solve = TaoSolve_BNLS;
 96:   tao->ops->setfromoptions = TaoSetFromOptions_BQNLS;
 97: 
 98:   bnk = (TAO_BNK*)tao->data;
 99:   bnk->update_type = BNK_UPDATE_STEP;
100:   bnk->computehessian = TaoBQNLSComputeHessian;
101:   bnk->computestep = TaoBQNLSComputeStep;
102: 
103:   bqnk = (TAO_BQNK*)bnk->ctx;
104:   MatSetOptionsPrefix(bqnk->B, "tao_bqnls_");
105:   MatSetType(bqnk->B, MATLMVMBFGS);
106:   return(0);
107: }