C++ AbstractLinAlgPack类代码示例

void LinAlgOpPack::Mp_StM(
  DMatrixSlice* C, value_type a
  ,const MatrixOp& B, BLAS_Cpp::Transp B_trans
  )
{
  using AbstractLinAlgPack::VectorSpace;
  using AbstractLinAlgPack::VectorDenseEncap;
  using AbstractLinAlgPack::MatrixOpGetGMS;
  using AbstractLinAlgPack::MatrixDenseEncap;
  const MatrixOpGetGMS
    *B_get_gms = dynamic_cast<const MatrixOpGetGMS*>(&B);
  if(B_get_gms) {
    DenseLinAlgPack::Mp_StM( C, a, MatrixDenseEncap(*B_get_gms)(), B_trans );		
  }
  else {
    const size_type num_cols = C->cols();
    VectorSpace::multi_vec_mut_ptr_t
      B_mv = ( B_trans == BLAS_Cpp::no_trans 
          ? B.space_cols() : B.space_rows()
          ).create_members(num_cols);
    assign(B_mv.get(),B,B_trans);
    for( size_type j = 1; j <= num_cols; ++j ) {
      DenseLinAlgPack::Vp_StV(&C->col(j),a,VectorDenseEncap(*B_mv->col(j))());
    }
  }
}

bool VariableBoundsTester::check_in_bounds(
  std::ostream* out, bool print_all_warnings, bool print_vectors
  ,const Vector& xL, const char xL_name[]
  ,const Vector& xU, const char xU_name[]
  ,const Vector& x,  const char x_name[]
  )
{
  using AbstractLinAlgPack::max_near_feas_step;

  if(out)
    *out
      << "\n*** Checking that variables are in bounds\n";

  VectorSpace::vec_mut_ptr_t zero = x.space().create_member(0.0);
  std::pair<value_type,value_type>
    u = max_near_feas_step( x, *zero, xL, xU, warning_tol() );
  if(u.first < 0.0) {
    if(out)
      *out << "\nWarning! the variables " << xL_name << " <= " << x_name << " <= " << xU_name
        << " are out of bounds by more than warning_tol = "	<< warning_tol() << "\n";
    u = max_near_feas_step( x, *zero, xL, xU, error_tol() );
    if(u.first < 0.0) {
      if(out)
        *out << "\nError! the variables " << xL_name << " <= " << x_name << " <= " << xU_name
          << " are out of bounds by more than error_tol = " << error_tol() << "\n";
      return false;
    }
  }
  return true;
}

void LinAlgOpPack::Vp_StPtMtV(
  DVectorSlice* y, value_type a
  ,const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
  ,const MatrixOp& M, BLAS_Cpp::Transp M_trans
  ,const DVectorSlice& x, value_type b
  )
{
  namespace mmp = MemMngPack;
  using BLAS_Cpp::no_trans;
  using AbstractLinAlgPack::VectorMutableDense;
  using AbstractLinAlgPack::VectorDenseMutableEncap;
  using AbstractLinAlgPack::Vp_StPtMtV;
  VectorSpace::space_ptr_t
    ay_space = ( M_trans == no_trans ? M.space_cols() : M.space_rows() ).space(P,P_trans);
  VectorSpace::vec_mut_ptr_t
    ay =  ( ay_space.get()
        ? ay_space->create_member()
        : Teuchos::rcp_implicit_cast<VectorMutable>(
          Teuchos::rcp(new VectorMutableDense(BLAS_Cpp::rows(P.rows(),P.cols(),P_trans)))
          ) ),
    ax = ( M_trans == no_trans ? M.space_rows() : M.space_cols() ).create_member();
  (VectorDenseMutableEncap(*ay))() = *y;
  (VectorDenseMutableEncap(*ax))() = x;
  Vp_StPtMtV( ay.get(), a, P, P_trans, M, M_trans, *ax, b );
  *y = VectorDenseMutableEncap(*ay)();
}

void LinAlgOpPack::Vp_StMtV(
  DVectorSlice* y, value_type a, const MatrixOp& M
  ,BLAS_Cpp::Transp M_trans, const SpVectorSlice& x, value_type b
  )
{
  using BLAS_Cpp::no_trans;
  using AbstractLinAlgPack::VectorDenseMutableEncap;
  VectorSpace::vec_mut_ptr_t
    ay = ( M_trans == no_trans ? M.space_cols() : M.space_rows() ).create_member();
  (VectorDenseMutableEncap(*ay))() = *y;
  AbstractLinAlgPack::Vp_StMtV( ay.get(), a, M, M_trans, x, b );
  *y = VectorDenseMutableEncap(*ay)();
}

void LinAlgOpPack::V_InvMtV(
  DVector* y, const MatrixOpNonsing& M
  ,BLAS_Cpp::Transp M_trans, const DVectorSlice& x
  )
{
  using BLAS_Cpp::trans;
  using AbstractLinAlgPack::VectorDenseMutableEncap;
  VectorSpace::vec_mut_ptr_t
    ay = ( M_trans == trans ? M.space_cols() : M.space_rows() ).create_member(),
    ax = ( M_trans == trans ? M.space_rows() : M.space_cols() ).create_member();
  (VectorDenseMutableEncap(*ax))() = x;
  AbstractLinAlgPack::V_InvMtV( ay.get(), M, M_trans, *ax );
  *y = VectorDenseMutableEncap(*ay)();
}

value_type MeritFuncCalc1DQuadratic::operator()(value_type alpha) const
{
  using AbstractLinAlgPack::Vp_StV;
  *x_ = *d_[0];
  value_type alpha_i = alpha;
  for( size_type i = 1; i <= p_; ++i, alpha_i *= alpha ) {
    Vp_StV( x_, alpha_i, *d_[i] );
  }
  return phi_( *x_ );
}

void ExampleNLPObjGrad::imp_calc_f(const Vector& x, bool newx
  , const ZeroOrderInfo& zero_order_info) const
{
  using AbstractLinAlgPack::dot;
  assert_is_initialized();
  f(); // assert f is set
  TEUCHOS_TEST_FOR_EXCEPTION( n() != x.dim(), std::length_error, "ExampleNLPObjGrad::imp_calc_f(...)"  );
  // f(x) = (obj_scale/2) * sum( x(i)^2, for i = 1..n )
  *zero_order_info.f = obj_scale_ / 2.0 * dot(x,x);
}

void MatrixHessianRelaxed::Vp_StMtV(
    DVectorSlice* y, value_type a, BLAS_Cpp::Transp M_trans
  , const SpVectorSlice& x, value_type b ) const
{
  using BLAS_Cpp::no_trans;
  using BLAS_Cpp::trans;
  using AbstractLinAlgPack::Vp_StMtV;
  //
  // y = b*y + a * M * x
  // 
  //   = b*y + a * [ H  0    ] * [ x1 ]
  //               [ 0  bigM ]   [ x2 ]
  //               
  // =>              
  //               
  // y1 = b*y1 + a*H*x1
  // 
  // y2 = b*y2 + bigM * x2
  //
  LinAlgOpPack::Vp_MtV_assert_sizes(y->size(),rows(),cols(),M_trans,x.size());

  DVectorSlice
    y1 = (*y)(1,n_);
  value_type
    &y2 = (*y)(n_+1);
  const SpVectorSlice
    x1 = x(1,n_);
  const SpVectorSlice::element_type
    *x2_ele = x.lookup_element(n_+1);
  const value_type
    x2 = x2_ele ? x2_ele->value() : 0.0;

  // y1 = b*y1 + a*H*x1
  Vp_StMtV( &y1, a, *H_, no_trans, x1, b );

  // y2 = b*y2 + bigM * x2
  if( b == 0.0 )
    y2 = bigM_ * x2;
  else
    y2 = b*y2 + bigM_ * x2;
  
}

void MultiVectorMutableCols::Vp_StMtV(
  VectorMutable* y, value_type a, BLAS_Cpp::Transp M_trans
  ,const SpVectorSlice& x, value_type b
  ) const
{
  using AbstractLinAlgPack::dot;

  // y = b*y
  LinAlgOpPack::Vt_S(y,b);

  if( M_trans == BLAS_Cpp::no_trans ) {
    //
    // y += a*M*x
    //
    // =>
    //
    // y += a * x(1) * M(:,1) + a * x(2) * M(:,2) + ...
    //
    SpVectorSlice::difference_type o = x.offset();
    for( SpVectorSlice::const_iterator itr = x.begin(); itr != x.end(); ++itr ) {
      const size_type j = itr->index() + o;
      LinAlgOpPack::Vp_StV( y, a * itr->value(), *col_vecs_[j-1] );
    }
  }
  else {
    //
    // y += a*M'*x
    //
    // =>
    //
    // y(1) += a M(:,1)*x
    // y(2) += a M(:,2)*x
    // ...
    //
    for( size_type j = 1; j <= col_vecs_.size(); ++j )
      y->set_ele(
        j
        ,y->get_ele(j) + a * dot(*col_vecs_[j-1],x)
        );
  }
}

bool ReducedHessianExactStd_Step::do_step(
    Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
  , poss_type assoc_step_poss)
{
  using Teuchos::dyn_cast;
  using DenseLinAlgPack::nonconst_sym;
  using AbstractLinAlgPack::Mp_StMtMtM;
  typedef AbstractLinAlgPack::MatrixSymDenseInitialize	MatrixSymDenseInitialize;
  typedef AbstractLinAlgPack::MatrixSymOp			MatrixSymOp;
  using ConstrainedOptPack::NLPSecondOrder;

  NLPAlgo	&algo	= rsqp_algo(_algo);
  NLPAlgoState	&s		= algo.rsqp_state();
  NLPSecondOrder
#ifdef _WINDOWS
        &nlp	= dynamic_cast<NLPSecondOrder&>(algo.nlp());
#else
        &nlp	= dyn_cast<NLPSecondOrder>(algo.nlp());
#endif
  MatrixSymOp
    *HL_sym_op = dynamic_cast<MatrixSymOp*>(&s.HL().get_k(0));

  EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
  std::ostream& out = algo.track().journal_out();

  // print step header.
  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
    using IterationPack::print_algorithm_step;
    print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
  }

  // problem size
  size_type	n		= nlp.n(),
        r		= nlp.r(),
        nind	= n - r;

  // Compute HL first (You may want to move this into its own step later?)

  if( !s.lambda().updated_k(-1) ) {
    if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
      out << "Initializing lambda_km1 = nlp.get_lambda_init ... \n";
    }
    nlp.get_init_lagrange_mult( &s.lambda().set_k(-1).v(), NULL );
    if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
      out << "||lambda_km1||inf = " << s.lambda().get_k(-1).norm_inf() << std::endl;
    }
    if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
      out << "lambda_km1 = \n" << s.lambda().get_k(-1)();
    }
  }

  nlp.set_HL(	HL_sym_op );
  nlp.calc_HL( s.x().get_k(0)(), s.lambda().get_k(-1)(), false );

  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
    s.HL().get_k(0).output( out << "\nHL_k = \n" );
  }

  // If rHL has already been updated for this iteration then just leave it.
  if( !s.rHL().updated_k(0) ) {

    if( !HL_sym_op ) {
      std::ostringstream omsg;
      omsg
        << "ReducedHessianExactStd_Step::do_step(...) : Error, "
        << "The matrix HL with the concrete type "
        << typeName(s.HL().get_k(0)) << " does not support the "
        << "MatrixSymOp iterface";
      throw std::logic_error( omsg.str() );
    }		

    MatrixSymDenseInitialize
      *rHL_sym_init = dynamic_cast<MatrixSymDenseInitialize*>(&s.rHL().set_k(0));
    if( !rHL_sym_init ) {
      std::ostringstream omsg;
      omsg
        << "ReducedHessianExactStd_Step::do_step(...) : Error, "
        << "The matrix rHL with the concrete type "
        << typeName(s.rHL().get_k(0)) << " does not support the "
        << "MatrixSymDenseInitialize iterface";
      throw std::logic_error( omsg.str() );
    }		

    // Compute the dense reduced Hessian
    DMatrix rHL_sym_store(nind,nind);
    DMatrixSliceSym rHL_sym(rHL_sym_store(),BLAS_Cpp::lower);
    Mp_StMtMtM( &rHL_sym, 1.0, MatrixSymOp::DUMMY_ARG, *HL_sym_op
          , s.Z().get_k(0), BLAS_Cpp::no_trans, 0.0 );

    if( (int)olevel >= (int)PRINT_ITERATION_QUANTITIES ) {
      out << "\nLower triangular partion of dense reduced Hessian (ignore nonzeros above diagonal):\nrHL_dense = \n" << rHL_sym_store(); 
    }
  
    // Set the reduced Hessain
    rHL_sym_init->initialize( rHL_sym );

    if( (int)olevel >= (int)PRINT_ITERATION_QUANTITIES ) {
      s.rHL().get_k(0).output( out << "\nrHL_k = \n" );
    }
  }
//.........这里部分代码省略.........

bool TangentialStepWithInequStd_Step::do_step(
  Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
  ,poss_type assoc_step_poss
  )
{

  using Teuchos::RCP;
  using Teuchos::dyn_cast;
  using ::fabs;
  using LinAlgOpPack::Vt_S;
  using LinAlgOpPack::V_VpV;
  using LinAlgOpPack::V_VmV;
  using LinAlgOpPack::Vp_StV;
  using LinAlgOpPack::Vp_V;
  using LinAlgOpPack::V_StV;
  using LinAlgOpPack::V_MtV;
//	using ConstrainedOptPack::min_abs;
  using AbstractLinAlgPack::max_near_feas_step;
  typedef VectorMutable::vec_mut_ptr_t   vec_mut_ptr_t;

  NLPAlgo &algo = rsqp_algo(_algo);
  NLPAlgoState &s = algo.rsqp_state();
  EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
  EJournalOutputLevel ns_olevel = algo.algo_cntr().null_space_journal_output_level();
  std::ostream &out = algo.track().journal_out();
  //const bool check_results = algo.algo_cntr().check_results();

  // print step header.
  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
    using IterationPack::print_algorithm_step;
    print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
  }

  // problem dimensions
  const size_type
    //n  = algo.nlp().n(),
    m  = algo.nlp().m(),
    r  = s.equ_decomp().size();

  // Get the iteration quantity container objects
  IterQuantityAccess<value_type>
    &alpha_iq = s.alpha(),
    &zeta_iq  = s.zeta(),
    &eta_iq   = s.eta();
  IterQuantityAccess<VectorMutable>
    &dl_iq      = dl_iq_(s),
    &du_iq      = du_iq_(s),
    &nu_iq      = s.nu(),
    *c_iq       = m > 0  ? &s.c()       : NULL,
    *lambda_iq  = m > 0  ? &s.lambda()  : NULL,
    &rGf_iq     = s.rGf(),
    &w_iq       = s.w(),
    &qp_grad_iq = s.qp_grad(),
    &py_iq      = s.py(),
    &pz_iq      = s.pz(),
    &Ypy_iq     = s.Ypy(),
    &Zpz_iq     = s.Zpz();
  IterQuantityAccess<MatrixOp>
    &Z_iq   = s.Z(),
    //*Uz_iq  = (m > r)  ? &s.Uz() : NULL,
    *Uy_iq  = (m > r)  ? &s.Uy() : NULL;
  IterQuantityAccess<MatrixSymOp>
    &rHL_iq = s.rHL();
  IterQuantityAccess<ActSetStats>
    &act_set_stats_iq = act_set_stats_(s);
  
  // Accessed/modified/updated (just some)
  VectorMutable  *Ypy_k = (m ? &Ypy_iq.get_k(0) : NULL);
  const MatrixOp  &Z_k   = Z_iq.get_k(0);
  VectorMutable  &pz_k  = pz_iq.set_k(0);
  VectorMutable  &Zpz_k = Zpz_iq.set_k(0);

  // Comupte qp_grad which is an approximation to rGf + Z'*HL*Y*py

  // qp_grad = rGf
  VectorMutable
    &qp_grad_k = ( qp_grad_iq.set_k(0) = rGf_iq.get_k(0) );

  // qp_grad += zeta * w
  if( w_iq.updated_k(0) ) {
    if(zeta_iq.updated_k(0))
      Vp_StV( &qp_grad_k, zeta_iq.get_k(0), w_iq.get_k(0) );
    else
      Vp_V( &qp_grad_k, w_iq.get_k(0) );
  }

  //
  // Set the bounds for:
  //
  //   dl <= Z*pz + Y*py <= du  ->  dl - Ypy <= Z*pz <= du - Ypz

  vec_mut_ptr_t
    bl = s.space_x().create_member(),
    bu = s.space_x().create_member();

  if(m) {
    // bl = dl_k - Ypy_k
    V_VmV( bl.get(), dl_iq.get_k(0), *Ypy_k );
    // bu = du_k - Ypy_k
    V_VmV( bu.get(), du_iq.get_k(0), *Ypy_k );
//.........这里部分代码省略.........

bool NLPDirectTester::finite_diff_check(
  NLPDirect         *nlp
  ,const Vector     &xo
  ,const Vector     *xl
  ,const Vector     *xu
  ,const Vector     *c
  ,const Vector     *Gf
  ,const Vector     *py
  ,const Vector     *rGf
  ,const MatrixOp   *GcU
  ,const MatrixOp   *D
  ,const MatrixOp   *Uz
  ,bool             print_all_warnings
  ,std::ostream     *out
  ) const
{

  using std::setw;
  using std::endl;
  using std::right;

  using AbstractLinAlgPack::sum;
  using AbstractLinAlgPack::dot;
  using AbstractLinAlgPack::Vp_StV;
  using AbstractLinAlgPack::random_vector;
  using AbstractLinAlgPack::assert_print_nan_inf;
  using LinAlgOpPack::V_StV;
  using LinAlgOpPack::V_StMtV;
  using LinAlgOpPack::Vp_MtV;
  using LinAlgOpPack::M_StM;
  using LinAlgOpPack::M_StMtM;

  typedef VectorSpace::vec_mut_ptr_t  vec_mut_ptr_t;

//  using AbstractLinAlgPack::TestingPack::CompareDenseVectors;
//  using AbstractLinAlgPack::TestingPack::CompareDenseSparseMatrices;

  using TestingHelperPack::update_success;

  bool success = true, preformed_fd;
  if(out) {
    *out << std::boolalpha
       << std::endl
       << "*********************************************************\n"
       << "*** NLPDirectTester::finite_diff_check(...) ***\n"
       << "*********************************************************\n";
  }

  const Range1D
    var_dep      = nlp->var_dep(),
    var_indep    = nlp->var_indep(),
    con_decomp   = nlp->con_decomp(),
    con_undecomp = nlp->con_undecomp();
  NLP::vec_space_ptr_t
    space_x = nlp->space_x(),
    space_c = nlp->space_c();

  // //////////////////////////////////////////////
  // Validate the input

  TEST_FOR_EXCEPTION(
    py && !c, std::invalid_argument
    ,"NLPDirectTester::finite_diff_check(...) : "
    "Error, if py!=NULL then c!=NULL must also be true!" );

  const CalcFiniteDiffProd
    &fd_deriv_prod = this->calc_fd_prod();

  const value_type
    rand_y_l = -1.0, rand_y_u = 1.0,
    small_num = ::sqrt(std::numeric_limits<value_type>::epsilon());

  try {

  // ///////////////////////////////////////////////
  // (1) Check Gf

  if(Gf) {
    switch( Gf_testing_method() ) {
      case FD_COMPUTE_ALL: {
        // Compute FDGf outright
        TEST_FOR_EXCEPT(true); // ToDo: update above!
        break;
      }
      case FD_DIRECTIONAL: {
        // Compute FDGF'*y using random y's
        if(out)
          *out
            << "\nComparing products Gf'*y with finite difference values FDGf'*y for "
            << "random y's ...\n";
        vec_mut_ptr_t y = space_x->create_member();
        value_type max_warning_viol = 0.0;
        int num_warning_viol = 0;
        const int num_fd_directions_used = ( num_fd_directions() > 0 ? num_fd_directions() : 1 );
        for( int direc_i = 1; direc_i <= num_fd_directions_used; ++direc_i ) {
          if( num_fd_directions() > 0 ) {
            random_vector( rand_y_l, rand_y_u, y.get() );
            if(out)
              *out
                << "\n****"
//.........这里部分代码省略.........

void MatrixSymPosDefLBFGS::secant_update(
  VectorMutable* s, VectorMutable* y, VectorMutable* Bs
  )
{
  using AbstractLinAlgPack::BFGS_sTy_suff_p_d;
  using AbstractLinAlgPack::dot;
  using LinAlgOpPack::V_MtV;
  using Teuchos::Workspace;
  Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();

  assert_initialized();

  // Check skipping the BFGS update
  const value_type
    sTy	      = dot(*s,*y);
  std::ostringstream omsg;
  if( !BFGS_sTy_suff_p_d(*s,*y,&sTy,&omsg,"MatrixSymPosDefLBFGS::secant_update(...)" ) ) {
    throw UpdateSkippedException( omsg.str() );	
  }

  try {

  // Update counters
  if( m_bar_ == m_ ) {
    // We are at the end of the storage so remove the oldest stored update
    // and move updates to make room for the new update.  This has to be done for the
    // the matrix to behave properly
    {for( size_type k = 1; k <= m_-1; ++k ) {
      S_->col(k) = S_->col(k+1);                            // Shift S.col() to the left
      Y_->col(k) = Y_->col(k+1);                            // Shift Y.col() to the left
      STY_.col(k)(1,m_-1) = STY_.col(k+1)(2,m_);            // Move submatrix STY(2,m-1,2,m-1) up and left
      STSYTY_.col(k)(k+1,m_) = STSYTY_.col(k+1)(k+2,m_+1);  // Move triangular submatrix STS(2,m-1,2,m-1) up and left
      STSYTY_.col(k+1)(1,k) = STSYTY_.col(k+2)(2,k+1);      // Move triangular submatrix YTY(2,m-1,2,m-1) up and left
    }}
    // ToDo: Create an abstract interface, call it MultiVectorShiftVecs, to rearrange S and Y all at once.
    // This will be important for parallel performance.
  }
  else {
    m_bar_++;
  }
  // Set the update vectors
  *S_->col(m_bar_) = *s;
  *Y_->col(m_bar_) = *y;

  // /////////////////////////////////////////////////////////////////////////////////////
  // Update S'Y
  //
  // Update the row and column m_bar
  //
  //	S'Y = 
  //
  //	[	s(1)'*y(1)		...		s(1)'*y(m_bar)		...		s(1)'*y(m_bar)		]
  //	[	.						.							.					] row
  //	[	s(m_bar)'*y(1)	...		s(m_bar)'*y(m_bar)	...		s(m_bar)'*y(m_bar)	] m_bar
  //	[	.						.							.					]
  //	[	s(m_bar)'*y(1)	...		s(m_bar)'*y(m_bar)	...		s(m_bar)'*y(m_bar)	]
  //
  //								col m_bar
  //
  // Therefore we set:
  //	(S'Y)(:,m_bar) =  S'*y(m_bar)
  //	(S'Y)(m_bar,:) =  s(m_bar)'*Y

  const multi_vec_ptr_t
    S = this->S(),
    Y = this->Y();

  VectorSpace::vec_mut_ptr_t
    t = S->space_rows().create_member();  // temporary workspace

  //	(S'Y)(:,m_bar) =  S'*y(m_bar)
  V_MtV( t.get(), *S, BLAS_Cpp::trans, *y );
  STY_.col(m_bar_)(1,m_bar_) = VectorDenseEncap(*t)();

  //	(S'Y)(m_bar,:)' =  Y'*s(m_bar)
  V_MtV( t.get(), *Y, BLAS_Cpp::trans, *s );
  STY_.row(m_bar_)(1,m_bar_) = VectorDenseEncap(*t)();

  // /////////////////////////////////////////////////////////////////
  // Update S'S
  //
  //	S'S = 
  //
  //	[	s(1)'*s(1)		...		symmetric					symmetric			]
  //	[	.						.							.					] row
  //	[	s(m_bar)'*s(1)	...		s(m_bar)'*s(m_bar)	...		symmetric			] m_bar
  //	[	.						.							.					]
  //	[	s(m_bar)'*s(1)	...		s(m_bar)'*s(m_bar)	...		s(m_bar)'*s(m_bar)	]
  //
  //								col m_bar
  //
  // Here we will update the lower triangular part of S'S.  To do this we
  // only need to compute:
  //		t = S'*s(m_bar) = { s(m_bar)' * [ s(1),..,s(m_bar),..,s(m_bar) ]  }'
  // then set the appropriate rows and columns of S'S.

  Workspace<value_type>   work_ws(wss,m_bar_);
  DVectorSlice                  work(&work_ws[0],work_ws.size());

  // work = S'*s(m_bar)
//.........这里部分代码省略.........

bool LineSearchFullStep_Step::do_step(Algorithm& _algo
  , poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss)
{
  using AbstractLinAlgPack::Vp_StV;
  using AbstractLinAlgPack::assert_print_nan_inf;
  using LinAlgOpPack::V_VpV;

  NLPAlgo        &algo   = rsqp_algo(_algo);
  NLPAlgoState   &s      = algo.rsqp_state();
  NLP            &nlp    = algo.nlp();

  const size_type
    m = nlp.m();

  EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
  std::ostream& out = algo.track().journal_out();

  // print step header.
  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
    using IterationPack::print_algorithm_step;
    print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
  }
  
  // alpha_k = 1.0
  IterQuantityAccess<value_type>
    &alpha_iq   = s.alpha();
  if( !alpha_iq.updated_k(0) )
    alpha_iq.set_k(0) = 1.0;

  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
    out	<< "\nf_k        = " << s.f().get_k(0);
    if(m)
      out << "\n||c_k||inf = " << s.c().get_k(0).norm_inf();
    out << "\nalpha_k    = " << alpha_iq.get_k(0) << std::endl;
  }

  // x_kp1 = x_k + d_k
  IterQuantityAccess<VectorMutable>  &x_iq = s.x();
  const Vector                       &x_k   = x_iq.get_k(0);
  VectorMutable                      &x_kp1 = x_iq.set_k(+1);
  x_kp1 = x_k;
  Vp_StV( &x_kp1, alpha_iq.get_k(0), s.d().get_k(0) );

  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
    out	<< "\n||x_kp1||inf   = " << s.x().get_k(+1).norm_inf() << std::endl;
  }
  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
    out << "\nx_kp1 =\n" << s.x().get_k(+1);
  }

  if(algo.algo_cntr().check_results()) {
    assert_print_nan_inf(
      x_kp1, "x_kp1",true
      ,int(olevel) >= int(PRINT_ALGORITHM_STEPS) ? &out : NULL
      );
    if( nlp.num_bounded_x() ) {
      if(!bounds_tester().check_in_bounds(
          int(olevel) >= int(PRINT_ALGORITHM_STEPS) ? &out : NULL
        , int(olevel) >= int(PRINT_VECTORS)					// print_all_warnings
        , int(olevel) >= int(PRINT_ITERATION_QUANTITIES)	// print_vectors
        , nlp.xl(), "xl"
        , nlp.xu(), "xu"
        , x_kp1, "x_kp1"
        ))
      {
        TEST_FOR_EXCEPTION(
          true, TestFailed
          ,"LineSearchFullStep_Step::do_step(...) : Error, "
          "the variables bounds xl <= x_k(+1) <= xu where violated!" );
      }
    }
  }

  // Calcuate f and c at the new point.
  nlp.unset_quantities();
  nlp.set_f( &s.f().set_k(+1) );
  if(m) nlp.set_c( &s.c().set_k(+1) );
  nlp.calc_f(x_kp1);
  if(m) nlp.calc_c( x_kp1, false );
  nlp.unset_quantities();

  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
    out	<< "\nf_kp1        = " << s.f().get_k(+1);
    if(m)
      out << "\n||c_kp1||inf = " << s.c().get_k(+1).norm_inf();
    out << std::endl;
  }

  if( m && static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
    out << "\nc_kp1 =\n" << s.c().get_k(+1); 
  }

  if(algo.algo_cntr().check_results()) {
    assert_print_nan_inf( s.f().get_k(+1), "f(x_kp1)", true, &out );
    if(m)
      assert_print_nan_inf( s.c().get_k(+1), "c(x_kp1)", true, &out );
  }

  return true;
}

void ExampleNLPFirstOrder::imp_calc_Gc(
    const Vector& x, bool newx, const FirstOrderInfo& first_order_info) const
{
    namespace rcp = MemMngPack;
    using Teuchos::dyn_cast;
    using AbstractLinAlgPack::Vp_S; // Should not have to do this!

    const index_type
    n = this->n(),
    m = this->m();
    const Range1D
    var_dep   = this->var_dep(),
    var_indep = this->var_indep();

    // Get references to aggregate C and N matrices (if allocated)
    MatrixOpNonsing
    *C_aggr = NULL;
    MatrixOp
    *N_aggr = NULL;
    BasisSystemComposite::get_C_N(
        first_order_info.Gc, &C_aggr, &N_aggr ); // Will return NULLs if Gc is not initialized

    // Allocate C and N matrix objects if not done yet!
    Teuchos::RCP<MatrixOpNonsing>
    C_ptr = Teuchos::null;
    Teuchos::RCP<MatrixOp>
    N_ptr = Teuchos::null;
    if( C_aggr == NULL ) {
        const VectorSpace::space_ptr_t
        space_x  = this->space_x(),
        space_xD = space_x->sub_space(var_dep);
        C_ptr  = Teuchos::rcp(new MatrixSymDiagStd(space_xD->create_member()));
        N_ptr  = Teuchos::rcp(new MatrixSymDiagStd(space_xD->create_member()));
        C_aggr = C_ptr.get();
        N_aggr = N_ptr.get();
    }

    // Get references to concreate C and N matrices
    MatrixSymDiagStd
    &C = dyn_cast<MatrixSymDiagStd>(*C_aggr);
    MatrixSymDiagStd
    &N = dyn_cast<MatrixSymDiagStd>(*N_aggr);
    // Get x = [ x_D' x_I ]
    Vector::vec_ptr_t
    x_D = x.sub_view(var_dep),
    x_I = x.sub_view(var_indep);
    // Set the diagonals of C and N (this is the only computation going on here)
    C.diag() = *x_I;          // C.diag = x_I - 1.0
    Vp_S( &C.diag(),  -1.0 ); // ...
    N.diag() = *x_D;          // N.diag = x_D - 10.0
    Vp_S( &N.diag(), -10.0 ); // ...
    // Initialize the matrix object Gc if not done so yet
    if( C_ptr.get() != NULL ) {
        BasisSystemComposite::initialize_Gc(
            this->space_x(), var_dep, var_indep
            ,this->space_c()
            ,C_ptr, N_ptr
            ,first_order_info.Gc
        );
    }
}