C++ AbstractLinAlgPack::transVtMtV方法代码示例

QPSolverStats::ESolutionType
QPSolverRelaxedQPKWIK::imp_solve_qp(
  std::ostream* out, EOutputLevel olevel, ERunTests test_what
  ,const Vector& g, const MatrixSymOp& G
  ,value_type etaL
  ,const Vector* dL, const Vector* dU
  ,const MatrixOp* E, BLAS_Cpp::Transp trans_E, const Vector* b
  ,const Vector* eL, const Vector* eU
  ,const MatrixOp* F, BLAS_Cpp::Transp trans_F, const Vector* f
  ,value_type* obj_d
  ,value_type* eta, VectorMutable* d
  ,VectorMutable* nu
  ,VectorMutable* mu, VectorMutable* Ed
  ,VectorMutable* lambda, VectorMutable* Fd
  )
{
  using Teuchos::dyn_cast;
  using DenseLinAlgPack::nonconst_tri_ele;
  using LinAlgOpPack::dot;
  using LinAlgOpPack::V_StV;
  using LinAlgOpPack::assign;
  using LinAlgOpPack::V_StV;
  using LinAlgOpPack::V_MtV;
  using AbstractLinAlgPack::EtaVector;
  using AbstractLinAlgPack::transVtMtV;
  using AbstractLinAlgPack::num_bounded;
  using ConstrainedOptPack::MatrixExtractInvCholFactor;

  // /////////////////////////
  // Map to QPKWIK input

  // Validate that rHL is of the proper type.
  const MatrixExtractInvCholFactor &cG
    = dyn_cast<const MatrixExtractInvCholFactor>(G);

  // Determine the number of sparse bounds on variables and inequalities.
  // By default set for the dense case
  const value_type inf = this->infinite_bound();
  const size_type
    nd              = d->dim(),
    m_in            = E  ? b->dim()                  : 0,
    m_eq            = F  ? f->dim()                  : 0,
    nvarbounds      = dL ? num_bounded(*dL,*dU,inf)  : 0,
    ninequbounds    = E  ? num_bounded(*eL,*eU,inf)  : 0,
    nequalities     = F  ? f->dim()                  : 0;

  // Determine if this is a QP with a structure different from the
  // one just solved.
  
  const bool same_qp_struct = (  N_ == nd && M1_ == nvarbounds && M2_ == ninequbounds && M3_ == nequalities );

  /////////////////////////////////////////////////////////////////
  // Set the input parameters to be sent to QPKWIKNEW

  // N
  N_ = nd;

  // M1
  M1_ = nvarbounds;

  // M2
  M2_ = ninequbounds;

  // M3
  M3_ = nequalities;

  // GRAD
  GRAD_ = VectorDenseEncap(g)();

  // UINV_AUG
  //
  // UINV_AUG = [ sqrt(bigM)  0  ]
  //            [ 0           L' ]
  //
  UINV_AUG_.resize(N_+1,N_+1);
  cG.extract_inv_chol( &nonconst_tri_ele( UINV_AUG_(2,N_+1,2,N_+1), BLAS_Cpp::upper ) );
  UINV_AUG_(1,1) = 1.0 / ::sqrt( NUMPARAM_[2] );
  UINV_AUG_.col(1)(2,N_+1) = 0.0;
  UINV_AUG_.row(1)(2,N_+1) = 0.0;

  // LDUINV_AUG
  LDUINV_AUG_ = UINV_AUG_.rows();

  // IBND, BL , BU, A, LDA, YPY

  IBND_INV_.resize( nd + m_in);
  std::fill( IBND_INV_.begin(), IBND_INV_.end(), 0 ); // Initialize the zero
  IBND_.resize( my_max( 1, M1_ + M2_ ) );
  BL_.resize( my_max( 1, M1_ + M2_ ) );
  BU_.resize( my_max( 1, M1_ + M2_ + M3_ ) );
  LDA_ = my_max( 1, M2_ + M3_ );
  A_.resize( LDA_, (  M2_ + M3_ > 0 ? N_ : 1 ) );
  YPY_.resize( my_max( 1, M1_ + M2_ ) );
  if(M1_)
    YPY_(1,M1_) = 0.0; // Must be for this QP interface

  // Initialize variable bound constraints
  if( dL ) {
    VectorDenseEncap dL_de(*dL);
    VectorDenseEncap dU_de(*dU);
//.........这里部分代码省略.........

bool ReducedHessianSecantUpdateLPBFGS_Strategy::perform_update(
  DVectorSlice* s_bfgs, DVectorSlice* y_bfgs, bool first_update
  ,std::ostream& out, EJournalOutputLevel olevel, NLPAlgo *algo, NLPAlgoState *s
  ,MatrixOp *rHL_k
  )
{
  using std::setw;
  using std::endl;
  using std::right;
  using Teuchos::dyn_cast;
  namespace rcp = MemMngPack;
  using Teuchos::RCP;
  using LinAlgOpPack::V_MtV;
  using DenseLinAlgPack::dot;
  using AbstractLinAlgPack::norm_inf;
  using AbstractLinAlgPack::transVtMtV;
  typedef ConstrainedOptPack::MatrixHessianSuperBasic MHSB_t;
  using Teuchos::Workspace;
  Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();

  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
    out << "\n*** (LPBFGS) Full limited memory BFGS to projected BFGS ...\n";
  }

#ifdef _WINDOWS
  MHSB_t &rHL_super = dynamic_cast<MHSB_t&>(*rHL_k);
#else
  MHSB_t &rHL_super = dyn_cast<MHSB_t>(*rHL_k);
#endif

  const size_type
    n    = algo->nlp().n(),
    r    = algo->nlp().r(),
    n_pz = n-r;

  bool do_projected_rHL_RR = false;

  // See if we still have a limited memory BFGS update matrix
  RCP<MatrixSymPosDefLBFGS> // We don't want this to be deleted until we are done with it
    lbfgs_rHL_RR = Teuchos::rcp_const_cast<MatrixSymPosDefLBFGS>(
      Teuchos::rcp_dynamic_cast<const MatrixSymPosDefLBFGS>(rHL_super.B_RR_ptr()) );

  if( lbfgs_rHL_RR.get() && rHL_super.Q_R().is_identity()  ) {
    //
    // We have a limited memory BFGS matrix and have not started projected BFGS updating
    // yet so lets determine if it is time to consider switching
    //
    // Check that the current update is sufficiently p.d. before we do anything
    const value_type
      sTy = dot(*s_bfgs,*y_bfgs),
      yTy = dot(*y_bfgs,*y_bfgs);
    if( !ConstrainedOptPack::BFGS_sTy_suff_p_d(
      *s_bfgs,*y_bfgs,&sTy
      ,  int(olevel) >= int(PRINT_ALGORITHM_STEPS) ? &out : NULL )
      && !proj_bfgs_updater().bfgs_update().use_dampening()
      )
    {
      if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
        out	<< "\nWarning!  use_damening == false so there is no way we can perform any"
            " kind of BFGS update (projected or not) so we will skip it!\n";
      }
      quasi_newton_stats_(*s).set_k(0).set_updated_stats(
        QuasiNewtonStats::INDEF_SKIPED );
      return true;
    }
    // Consider if we can even look at the active set yet.
    const bool consider_switch  = lbfgs_rHL_RR->num_secant_updates() >= min_num_updates_proj_start();
    if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
      out << "\nnum_previous_updates = " << lbfgs_rHL_RR->num_secant_updates()
        << ( consider_switch ? " >= " : " < " )
        << "min_num_updates_proj_start = " << min_num_updates_proj_start()
        << ( consider_switch
           ? "\nConsidering if we should switch to projected BFGS updating of superbasics ...\n"
           : "\nNot time to consider switching to projected BFGS updating of superbasics yet!" );
    }
    if( consider_switch ) {
      // 
      // Our job here is to determine if it is time to switch to projected projected
      // BFGS updating.
      //
      if( act_set_stats_(*s).updated_k(-1) ) {
        if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
          out	<< "\nDetermining if projected BFGS updating of superbasics should begin ...\n";
        }
        // Determine if the active set has calmed down enough
        const SpVector
          &nu_km1 = s->nu().get_k(-1);
        const SpVectorSlice
          nu_indep = nu_km1(s->var_indep());
        const bool 
          act_set_calmed_down
          = PBFGSPack::act_set_calmed_down(
            act_set_stats_(*s).get_k(-1)
            ,proj_bfgs_updater().act_set_frac_proj_start()
            ,olevel,out
            ),
          max_num_updates_exceeded
          = lbfgs_rHL_RR->m_bar() >= max_num_updates_proj_start();
        do_projected_rHL_RR = act_set_calmed_down || max_num_updates_exceeded;
        if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
//.........这里部分代码省略.........