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

bool TangentialStepIP_Step::do_step(
  Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
  ,poss_type assoc_step_poss
  )
  {
  using BLAS_Cpp::no_trans;
  using Teuchos::dyn_cast;
  using AbstractLinAlgPack::assert_print_nan_inf;
  using LinAlgOpPack::Vt_S;
  using LinAlgOpPack::Vp_StV;
  using LinAlgOpPack::V_StV;
  using LinAlgOpPack::V_MtV;
  using LinAlgOpPack::V_InvMtV;
   using LinAlgOpPack::M_StM;
  using LinAlgOpPack::Mp_StM;
  using LinAlgOpPack::assign;

  NLPAlgo	&algo	= rsqp_algo(_algo);
  IpState	    &s      = dyn_cast<IpState>(_algo.state());

  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 );
  }

  // Compute qp_grad which is an approximation to rGf + Z'*(mu*(invXu*e-invXl*e) + no_cross_term
  // minimize round off error by calc'ing Z'*(Gf + mu*(invXu*e-invXl*e))

  // qp_grad_k = Z'*(Gf + mu*(invXu*e-invXl*e))
  const MatrixSymDiagStd  &invXu = s.invXu().get_k(0);
  const MatrixSymDiagStd  &invXl = s.invXl().get_k(0);
  const value_type            &mu    = s.barrier_parameter().get_k(0);
  const MatrixOp          &Z_k   = s.Z().get_k(0);

  Teuchos::RCP<VectorMutable> rhs = s.Gf().get_k(0).clone();
  Vp_StV( rhs.get(), mu,      invXu.diag() );
  Vp_StV( rhs.get(), -1.0*mu, invXl.diag() );
  
  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) 
    {
    out << "\n||Gf_k + mu_k*(invXu_k-invXl_k)||inf = " << rhs->norm_inf() << std::endl;
    }
  if( (int)olevel >= (int)PRINT_VECTORS)
    {
    out << "\nGf_k + mu_k*(invXu_k-invXl_k) =\n" << *rhs;
    }

  VectorMutable &qp_grad_k = s.qp_grad().set_k(0);
  V_MtV(&qp_grad_k, Z_k, BLAS_Cpp::trans, *rhs);
  
  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) 
    {
    out << "\n||Z_k'*(Gf_k + mu_k*(invXu_k-invXl_k))||inf = " << qp_grad_k.norm_inf() << std::endl;
    }
  if( (int)olevel >= (int)PRINT_VECTORS )
    {
    out << "\nZ_k'*(Gf_k + mu_k*(invXu_k-invXl_k)) =\n" << qp_grad_k;
    }

  // error check for cross term
  value_type         &zeta    = s.zeta().set_k(0);
  const Vector &w_sigma = s.w_sigma().get_k(0);
  
  // need code to calculate damping parameter
  zeta = 1.0;

  Vp_StV(&qp_grad_k, zeta, w_sigma);

  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) 
    {
    out << "\n||qp_grad_k||inf = " << qp_grad_k.norm_inf() << std::endl;
    }
  if( (int)olevel >= (int)PRINT_VECTORS ) 
    {
    out << "\nqp_grad_k =\n" << qp_grad_k;
    }

  // build the "Hessian" term B = rHL + rHB
  // should this be MatrixSymOpNonsing
  const MatrixSymOp      &rHL_k = s.rHL().get_k(0);
  const MatrixSymOp      &rHB_k = s.rHB().get_k(0);
  MatrixSymOpNonsing &B_k   = dyn_cast<MatrixSymOpNonsing>(s.B().set_k(0));
  if (B_k.cols() != Z_k.cols())
    {
    // Initialize space in rHB
    dyn_cast<MatrixSymInitDiag>(B_k).init_identity(Z_k.space_rows(), 0.0);
    }

  //	M_StM(&B_k, 1.0, rHL_k, no_trans);
  assign(&B_k, rHL_k, BLAS_Cpp::no_trans);
  if( (int)olevel >= (int)PRINT_VECTORS ) 
    {
    out << "\nB_k = rHL_k =\n" << B_k;
    }
  Mp_StM(&B_k, 1.0, rHB_k, BLAS_Cpp::no_trans);
  if( (int)olevel >= (int)PRINT_VECTORS ) 
//.........这里部分代码省略.........

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

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****"
//.........这里部分代码省略.........

bool CheckConvergenceStd_Strategy::Converged(
  Algorithm& _algo
  )
  {
  using AbstractLinAlgPack::assert_print_nan_inf;
  using AbstractLinAlgPack::combined_nu_comp_err;
  
  NLPAlgo        &algo = rsqp_algo(_algo);
  NLPAlgoState   &s    = algo.rsqp_state();
  NLP            &nlp  = algo.nlp();

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

  const size_type
    n  = nlp.n(),
    m  = nlp.m(),
    nb = nlp.num_bounded_x();

  // Get the iteration quantities
  IterQuantityAccess<value_type>
    &opt_kkt_err_iq  = s.opt_kkt_err(),
    &feas_kkt_err_iq = s.feas_kkt_err(),
      &comp_kkt_err_iq = s.comp_kkt_err();
  
  IterQuantityAccess<VectorMutable>
    &x_iq       = s.x(),
    &d_iq       = s.d(),
    &Gf_iq      = s.Gf(),
    *c_iq       = m     ? &s.c()      : NULL,
    *rGL_iq     = n > m ? &s.rGL()    : NULL,
    *GL_iq      = n > m ? &s.GL()     : NULL,
    *nu_iq      = n > m ? &s.nu()     : NULL;

  // opt_err = (||rGL||inf or ||GL||) / (||Gf|| + scale_kkt_factor)
  value_type 
    norm_inf_Gf_k = 0.0,
    norm_inf_GLrGL_k = 0.0;

  if( n > m && scale_opt_error_by_Gf() && Gf_iq.updated_k(0) ) {
    assert_print_nan_inf(
      norm_inf_Gf_k = Gf_iq.get_k(0).norm_inf(),
      "||Gf_k||inf",true,&out
      );
  }

  // NOTE:
  // The strategy object CheckConvergenceIP_Strategy assumes
  // that this will always be the gradient of the lagrangian
  // of the original problem, not the gradient of the lagrangian
  // for psi. (don't use augmented nlp info here)
  if( n > m ) {
    if( opt_error_check() == OPT_ERROR_REDUCED_GRADIENT_LAGR ) {
      assert_print_nan_inf( norm_inf_GLrGL_k = rGL_iq->get_k(0).norm_inf(),
                  "||rGL_k||inf",true,&out);
    }
    else {
      assert_print_nan_inf( norm_inf_GLrGL_k = GL_iq->get_k(0).norm_inf(),
                  "||GL_k||inf",true,&out);
    }
  }

  const value_type
    opt_scale_factor = 1.0 + norm_inf_Gf_k,
    opt_err = norm_inf_GLrGL_k / opt_scale_factor;
  
  // feas_err
  const value_type feas_err = ( ( m ? c_iq->get_k(0).norm_inf() : 0.0 ) );

  // comp_err
  value_type comp_err = 0.0;
  if ( n > m ) {
    if (nb > 0) {
      comp_err = combined_nu_comp_err(nu_iq->get_k(0), x_iq.get_k(0), nlp.xl(), nlp.xu());
    }
    if(m) {
      assert_print_nan_inf( feas_err,"||c_k||inf",true,&out);
    }
  }

  // scaling factors
  const value_type 
    scale_opt_factor = CalculateScalingFactor(s, scale_opt_error_by()),
    scale_feas_factor = CalculateScalingFactor(s, scale_feas_error_by()),
    scale_comp_factor = CalculateScalingFactor(s, scale_comp_error_by());
  
  // kkt_err
  const value_type
    opt_kkt_err_k  = opt_err/scale_opt_factor,
     feas_kkt_err_k = feas_err/scale_feas_factor,
    comp_kkt_err_k = comp_err/scale_comp_factor;

  // update the iteration quantities
  if(n > m) opt_kkt_err_iq.set_k(0) = opt_kkt_err_k;
  feas_kkt_err_iq.set_k(0) = feas_kkt_err_k;
  comp_kkt_err_iq.set_k(0) = comp_kkt_err_k;

  // step_err
  value_type step_err = 0.0;
  if( d_iq.updated_k(0) ) {
//.........这里部分代码省略.........

bool PostEvalNewPointBarrier_Step::do_step(
  Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
  ,poss_type assoc_step_poss
  )
  {
  using Teuchos::dyn_cast;
  using IterationPack::print_algorithm_step;
  using AbstractLinAlgPack::inv_of_difference;
  using AbstractLinAlgPack::correct_upper_bound_multipliers;
  using AbstractLinAlgPack::correct_lower_bound_multipliers;
  using LinAlgOpPack::Vp_StV;

  NLPAlgo            &algo   = dyn_cast<NLPAlgo>(_algo);
  IpState             &s      = dyn_cast<IpState>(_algo.state());
  NLP                 &nlp    = algo.nlp();
  
  EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
  std::ostream& out = algo.track().journal_out();
  
  if(!nlp.is_initialized())
    nlp.initialize(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 );
    }

  IterQuantityAccess<VectorMutable> &x_iq = s.x();
  IterQuantityAccess<MatrixSymDiagStd> &Vl_iq = s.Vl();
  IterQuantityAccess<MatrixSymDiagStd> &Vu_iq = s.Vu();

  ///***********************************************************
  // Calculate invXl = diag(1/(x-xl)) 
  //  and invXu = diag(1/(xu-x)) matrices
  ///***********************************************************

  // get references to x, invXl, and invXu
  VectorMutable& x = x_iq.get_k(0);
  MatrixSymDiagStd& invXu = s.invXu().set_k(0);
  MatrixSymDiagStd& invXl = s.invXl().set_k(0);
  
  //std::cout << "xu=\n";
  //nlp.xu().output(std::cout);

  inv_of_difference(1.0, nlp.xu(), x, &invXu.diag());
  inv_of_difference(1.0, x, nlp.xl(), &invXl.diag());

  //std::cout << "invXu=\v";
  //invXu.output(std::cout);

  //std::cout << "\ninvXl=\v";
  //invXl.output(std::cout);
  
  // Check for divide by zeros - 
    using AbstractLinAlgPack::assert_print_nan_inf;
    assert_print_nan_inf(invXu.diag(), "invXu", true, &out); 
    assert_print_nan_inf(invXl.diag(), "invXl", true, &out); 
  // These should never go negative either - could be a useful check

  // Initialize Vu and Vl if necessary
  if ( /*!Vu_iq.updated_k(0) */ Vu_iq.last_updated() == IterQuantity::NONE_UPDATED )
    {
    VectorMutable& vu = Vu_iq.set_k(0).diag();		
    vu = 0;
    Vp_StV(&vu, s.barrier_parameter().get_k(-1), invXu.diag());

    if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) 
      {
      out << "\nInitialize Vu with barrier_parameter * invXu ...\n";
      }
    }

if ( /*!Vl_iq.updated_k(0) */ Vl_iq.last_updated() == IterQuantity::NONE_UPDATED  )
    {
    VectorMutable& vl = Vl_iq.set_k(0).diag();
    vl = 0;
    Vp_StV(&vl, s.barrier_parameter().get_k(-1), invXl.diag());

    if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) 
      {
      out << "\nInitialize Vl with barrier_parameter * invXl ...\n";
      }
    }

  if (s.num_basis().updated_k(0))
    {
    // Basis changed, reorder Vl and Vu
    if (Vu_iq.updated_k(-1))
      { Vu_iq.set_k(0,-1); }
    if (Vl_iq.updated_k(-1))
      { Vl_iq.set_k(0,-1); }
      
    VectorMutable& vu = Vu_iq.set_k(0).diag();
    VectorMutable& vl = Vl_iq.set_k(0).diag();

    s.P_var_last().permute( BLAS_Cpp::trans, &vu ); // Permute back to original order
    s.P_var_last().permute( BLAS_Cpp::trans, &vl ); // Permute back to original order

//.........这里部分代码省略.........

bool EvalNewPointStd_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 AbstractLinAlgPack::assert_print_nan_inf;
  using IterationPack::print_algorithm_step;
  using NLPInterfacePack::NLPFirstOrder;

  NLPAlgo         &algo   = rsqp_algo(_algo);
  NLPAlgoState    &s      = algo.rsqp_state();
  NLPFirstOrder   &nlp    = dyn_cast<NLPFirstOrder>(algo.nlp());

  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();

  // 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 );
  }

  if(!nlp.is_initialized())
    nlp.initialize(algo.algo_cntr().check_results());

  Teuchos::VerboseObjectTempState<NLP>
    nlpOutputTempState(rcp(&nlp,false),Teuchos::getFancyOStream(rcp(&out,false)),convertToVerbLevel(olevel));

  const size_type
    n  = nlp.n(),
    nb = nlp.num_bounded_x(),
    m  = nlp.m();
  size_type
    r  = 0;

  // Get the iteration quantity container objects
  IterQuantityAccess<value_type>
    &f_iq   = s.f();
  IterQuantityAccess<VectorMutable>
    &x_iq   = s.x(),
    *c_iq   = m > 0  ? &s.c() : NULL,
    &Gf_iq  = s.Gf();
  IterQuantityAccess<MatrixOp>
    *Gc_iq  = m  > 0 ? &s.Gc() : NULL,
    *Z_iq   = NULL,
    *Y_iq   = NULL,
    *Uz_iq  = NULL,
    *Uy_iq  = NULL;
  IterQuantityAccess<MatrixOpNonsing>
    *R_iq   = NULL;

  MatrixOp::EMatNormType mat_nrm_inf = MatrixOp::MAT_NORM_INF;
  const bool calc_matrix_norms = algo.algo_cntr().calc_matrix_norms();
  const bool calc_matrix_info_null_space_only = algo.algo_cntr().calc_matrix_info_null_space_only();
  
  if( x_iq.last_updated() == IterQuantity::NONE_UPDATED ) {
    if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
      out << "\nx is not updated for any k so set x_k = nlp.xinit() ...\n";
    }
    x_iq.set_k(0) = nlp.xinit();
  }
  
  // Validate x
  if( nb && algo.algo_cntr().check_results() ) {
    assert_print_nan_inf(
      x_iq.get_k(0), "x_k", true
      , int(olevel) >= int(PRINT_ALGORITHM_STEPS) ? &out : NULL
      );
    if( nlp.num_bounded_x() > 0 ) {
      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_iq.get_k(0),   "x_k"
           ))
      {
        TEUCHOS_TEST_FOR_EXCEPTION(
          true, TestFailed
          ,"EvalNewPointStd_Step::do_step(...) : Error, "
          "the variables bounds xl <= x_k <= xu where violated!" );
      }
    }
  }

  Vector &x = x_iq.get_k(0);

  Range1D  var_dep(Range1D::INVALID), var_indep(Range1D::INVALID);
  if( s.get_decomp_sys().get() ) {
    const ConstrainedOptPack::DecompositionSystemVarReduct
      *decomp_sys_vr = dynamic_cast<ConstrainedOptPack::DecompositionSystemVarReduct*>(&s.decomp_sys());
    if(decomp_sys_vr) {
      var_dep   = decomp_sys_vr->var_dep();
      var_indep = decomp_sys_vr->var_indep();
    }
    s.var_dep(var_dep);
//.........这里部分代码省略.........

bool PreEvalNewPointBarrier_Step::do_step(
  Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
  ,poss_type assoc_step_poss
  )
  {
  using Teuchos::dyn_cast;
  using IterationPack::print_algorithm_step;
  using AbstractLinAlgPack::force_in_bounds_buffer;

  NLPAlgo            &algo   = dyn_cast<NLPAlgo>(_algo);
  IpState             &s      = dyn_cast<IpState>(_algo.state());
  NLP                 &nlp    = algo.nlp();
  NLPFirstOrder   *nlp_foi = dynamic_cast<NLPFirstOrder*>(&nlp);
  
  EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
  std::ostream& out = algo.track().journal_out();
  
  if(!nlp.is_initialized())
    nlp.initialize(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 );
    }

  IterQuantityAccess<value_type>     &barrier_parameter_iq = s.barrier_parameter();
  IterQuantityAccess<VectorMutable>  &x_iq  = s.x();

  if( x_iq.last_updated() == IterQuantity::NONE_UPDATED ) 
    {
    if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) 
      {
      out << "\nInitialize x with x_k = nlp.xinit() ...\n"
        << " and push x_k within bounds.\n";
      }
    VectorMutable& x_k = x_iq.set_k(0) = nlp.xinit();
  
    // apply transformation operator to push x sufficiently within bounds
    force_in_bounds_buffer(relative_bound_push_, 
                 absolute_bound_push_,
                 nlp.xl(),
                 nlp.xu(),
                 &x_k);

    // evaluate the func and constraints
    IterQuantityAccess<value_type>
      &f_iq    = s.f();
    IterQuantityAccess<VectorMutable>
      &Gf_iq   = s.Gf(),
      *c_iq    = nlp.m() > 0 ? &s.c() : NULL;
    IterQuantityAccess<MatrixOp>
      *Gc_iq   = nlp_foi ? &s.Gc() : NULL;

    using AbstractLinAlgPack::assert_print_nan_inf;
    assert_print_nan_inf(x_k, "x", true, NULL); // With throw exception if Inf or NaN!

    // Wipe out storage for computed iteration quantities (just in case?) : RAB: 7/29/2002
    if(f_iq.updated_k(0))
      f_iq.set_not_updated_k(0);
    if(Gf_iq.updated_k(0))
      Gf_iq.set_not_updated_k(0);
    if (c_iq)
      {
      if(c_iq->updated_k(0))
        c_iq->set_not_updated_k(0);
      }
    if (nlp_foi)
      {
      if(Gc_iq->updated_k(0))
        Gc_iq->set_not_updated_k(0);
      }
    }

  if (barrier_parameter_iq.last_updated() == IterQuantity::NONE_UPDATED)
    {
    barrier_parameter_iq.set_k(-1) = 0.1; // RAB: 7/29/2002: This should be parameterized! (allow user to set this!)
    }

  // Print vector information
  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) 
    {
    out << "x_k =\n" << x_iq.get_k(0);
    }

  return true;
  }

bool MatrixOpNonsingTester::test_matrix(
  const MatrixOpNonsing   &M
  ,const char                     M_name[]
  ,std::ostream                   *out
  )
{
  namespace rcp = MemMngPack;
  using BLAS_Cpp::no_trans;
  using BLAS_Cpp::trans;
  using BLAS_Cpp::left;
  using BLAS_Cpp::right;
  using AbstractLinAlgPack::sum;
  using AbstractLinAlgPack::dot;
  using AbstractLinAlgPack::Vp_StV;
  using AbstractLinAlgPack::assert_print_nan_inf;
  using AbstractLinAlgPack::random_vector;
  using LinAlgOpPack::V_StMtV;
  using LinAlgOpPack::V_MtV;
  using LinAlgOpPack::V_StV;
  using LinAlgOpPack::V_VpV;
  using LinAlgOpPack::Vp_V;
  
  // ToDo: Check the preconditions
  
  bool success = true, result, lresult;
  const value_type
    rand_y_l  = -1.0,
    rand_y_u  = 1.0,
    small_num = ::pow(std::numeric_limits<value_type>::epsilon(),0.25),
    alpha     = 2.0;
  
  //
  // Perform the tests
  //

  if(out && print_tests() >= PRINT_BASIC)
    *out
      << "\nCheck: alpha*op(op(inv("<<M_name<<"))*op("<<M_name<<"))*v == alpha*v ...";
  if(out && print_tests() > PRINT_BASIC)
    *out << std::endl;

  VectorSpace::vec_mut_ptr_t
    v_c1 = M.space_cols().create_member(),
    v_c2 = M.space_cols().create_member(),
    v_r1 = M.space_rows().create_member(),
    v_r2 = M.space_rows().create_member();

  // Side of the matrix inverse	
  const BLAS_Cpp::Side    a_side[2]  = { BLAS_Cpp::left,     BLAS_Cpp::right };
  // If the matrices are transposed or not
  const BLAS_Cpp::Transp  a_trans[2] = { BLAS_Cpp::no_trans, BLAS_Cpp::trans };

  for( int side_i = 0; side_i < 2; ++side_i ) {
    for( int trans_i = 0; trans_i < 2; ++trans_i ) {
      const BLAS_Cpp::Side    t_side  = a_side[side_i];
      const BLAS_Cpp::Transp  t_trans = a_trans[trans_i];
      if(out && print_tests() >= PRINT_MORE)
        *out
          << "\n" << side_i+1 << "." << trans_i+1 << ") "
          << "Check: (t2 = "<<(t_side==left?"inv(":"alpha * ")<< M_name<<(t_trans==trans?"\'":"")<<(t_side==left?")":"")
          << " * (t1 = "<<(t_side==right?"inv(":"alpha * ")<<M_name<<(t_trans==trans?"\'":"")<<(t_side==right?")":"")
          << " * v)) == alpha * v ...";
      if(out && print_tests() > PRINT_MORE)
        *out << std::endl;
      result = true;
      VectorMutable
        *v  = NULL,
        *t1 = NULL,
        *t2 = NULL;
      if( (t_side == left && t_trans == no_trans) || (t_side == right && t_trans == trans) ) {
        // (inv(R)*R*v or R'*inv(R')*v
        v  = v_r1.get();
        t1 = v_c1.get();
        t2 = v_r2.get();
      }
      else {
        // (inv(R')*R'*v or R*inv(R)*v
        v  = v_c1.get();
        t1 = v_r1.get();
        t2 = v_c2.get();
      }
      for( int k = 1; k <= num_random_tests(); ++k ) {
        lresult = true;
        random_vector( rand_y_l, rand_y_u, v );
          if(out && print_tests() >= PRINT_ALL) {
          *out
            << "\n"<<side_i+1<<"."<<trans_i+1<<"."<<k<<") random vector " << k
            << " ( ||v||_1 / n = " << (v->norm_1() / v->dim()) << " )\n";
          if(dump_all() && print_tests() >= PRINT_ALL)
            *out << "\nv =\n" << *v;
        }
        // t1
        if( t_side == right ) {
          // t1 = inv(op(M))*v
          V_InvMtV( t1, M, t_trans, *v );
        }
        else {
          // t1 = alpha*op(M)*v
          V_StMtV( t1, alpha, M, t_trans, *v );
        }
//.........这里部分代码省略.........