C++ DMatrixSlice::col方法代码示例

bool DenseLinAlgPack::assert_print_nan_inf( const DMatrixSlice& m
  , const std::string & name, bool throw_excpt, std::ostream* out )
{

  bool has_nan_or_inf = false;
  bool printed_header = false;

  for( size_type j = 1; j <= m.cols(); ++j ) {
    const DVectorSlice& v = m.col(j);
    for( DVectorSlice::const_iterator v_itr = v.begin(); v_itr != v.end(); ++v_itr ) {
      if( RTOp_is_nan_inf(*v_itr) ) {
        if(out) {
          if(!printed_header) {
            *out
              << "The matrix \"" << name
              << "\" has the following NaN or Inf entries\n";
            printed_header = true;
          }
          *out
            << name << "(" << v_itr - v.begin() + 1 << "," << j << ") = "
            << *v_itr << std::endl;
        }
        has_nan_or_inf = true;
      }
    }
  }

  if( has_nan_or_inf && throw_excpt ) {
    if(out)
      out->flush();
    std::ostringstream omsg;
    omsg
      << "assert_print_nan_inf(...) : Error, the matrix named "
      << name << " has at least one element which is NaN or Inf";
    throw NaNInfException( omsg.str() );
  }

  return has_nan_or_inf;
}

void DenseLinAlgPack::delete_row_col( size_type kd, DMatrixSliceTriEle* tri_M )
{
  // Validate input
  TEUCHOS_TEST_FOR_EXCEPT( !(  tri_M  ) );
  TEUCHOS_TEST_FOR_EXCEPT( !(  tri_M->rows()  ) );
  TEUCHOS_TEST_FOR_EXCEPT( !(  1 <= kd && kd <= tri_M->rows()  ) );

  DMatrixSlice   M = tri_M->gms();
  const size_type  n = M.rows();

  if( tri_M->uplo() == BLAS_Cpp::lower ) {
    // Move M31 up one row at a time
    if( 1 < kd && kd < n ) {
      Range1D rng(1,kd-1);
      for( size_type i = kd; i < n; ++i )
        M.row(i)(rng) = M.row(i+1)(rng);
    }
    // Move M33 up and to the left one column at a time
    if( kd < n ) {
      for( size_type i = kd; i < n; ++i )
        M.col(i)(i,n-1) = M.col(i+1)(i+1,n);
    }
  }
  else if(  tri_M->uplo() == BLAS_Cpp::upper ) {
    // Move M13 left one column at a time.
    if( 1 < kd && kd < n ) {
      Range1D rng(1,kd-1);
      for( size_type j = kd; j < n; ++j )
        M.col(j)(rng) = M.col(j+1)(rng);
    }
    // Move the updated U33 up and left one column at a time.
    if( kd < n ) {
      for( size_type j = kd; j < n; ++j )
        M.col(j)(kd,j) = M.col(j+1)(kd+1,j+1);
    }
  }
  else {
    TEUCHOS_TEST_FOR_EXCEPT(true); // Invalid input
  }
}

bool DenseLinAlgPack::comp(const DMatrixSlice& gms, value_type alpha)
{
  for(size_type i = 1; i < gms.cols(); ++i)
    if( !comp( gms.col(i) , alpha ) ) return false;
  return true;
}

inline 
/** \brief . */
const DVectorSlice col(const DMatrixSlice& gms, BLAS_Cpp::Transp trans, size_type j) {
  return (trans ==  BLAS_Cpp::no_trans) ? gms.col(j) : gms.row(j);
} 

inline 
/** \brief . */
DVectorSlice row(DMatrixSlice& gms, BLAS_Cpp::Transp trans, size_type i) {
  return (trans ==  BLAS_Cpp::no_trans) ? gms.row(i) : gms.col(i);
} 

//.........这里部分代码省略.........
            rHL_RR->init_identity( n_pz_R, rHL_scale );
            if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
              out << "\nrHL_RR after intialized to rHL_RR = rHL_scale*I but before performing current BFGS update:\nrHL_RR =\n"
                << dynamic_cast<MatrixOp&>(*rHL_RR); // I know this is okay!
            }
            // Initialize rHL_XX = rHL_scale*I
#ifdef _WINDOWS
            MatrixSymInitDiag
              &rHL_XX = dynamic_cast<MatrixSymInitDiag&>(
                const_cast<MatrixSymOp&>(*rHL_super.B_XX_ptr()));
#else
            MatrixSymInitDiag
              &rHL_XX = dyn_cast<MatrixSymInitDiag>(
                const_cast<MatrixSymOp&>(*rHL_super.B_XX_ptr()));
#endif
            rHL_XX.init_identity( n_pz_X, rHL_scale );
            // Reinitialize rHL
            rHL_super.initialize(
              n_pz, n_pz_R, &i_x_free[0], &i_x_fixed[0], &bnd_fixed[0]
              ,Teuchos::rcp_const_cast<const MatrixSymWithOpFactorized>(
                Teuchos::rcp_dynamic_cast<MatrixSymWithOpFactorized>(rHL_RR))
              ,NULL,BLAS_Cpp::no_trans,rHL_super.B_XX_ptr()
              );
            //
            // Perform the current BFGS update first
            //
            MatrixSymOp
              &rHL_RR_op = dynamic_cast<MatrixSymOp&>(*rHL_RR);
            const GenPermMatrixSlice
              &Q_R = rHL_super.Q_R(),
              &Q_X = rHL_super.Q_X();
            // Get projected BFGS update vectors y_bfgs_R, s_bfgs_R
            Workspace<value_type>
              y_bfgs_R_ws(wss,Q_R.cols()),
              s_bfgs_R_ws(wss,Q_R.cols()),
              y_bfgs_X_ws(wss,Q_X.cols()),
              s_bfgs_X_ws(wss,Q_X.cols());
            DVectorSlice y_bfgs_R(&y_bfgs_R_ws[0],y_bfgs_R_ws.size());
            DVectorSlice s_bfgs_R(&s_bfgs_R_ws[0],s_bfgs_R_ws.size());
            DVectorSlice y_bfgs_X(&y_bfgs_X_ws[0],y_bfgs_X_ws.size());
            DVectorSlice s_bfgs_X(&s_bfgs_X_ws[0],s_bfgs_X_ws.size());
            V_MtV( &y_bfgs_R, Q_R, BLAS_Cpp::trans, *y_bfgs );  // y_bfgs_R = Q_R'*y_bfgs
            V_MtV( &s_bfgs_R, Q_R, BLAS_Cpp::trans, *s_bfgs );  // s_bfgs_R = Q_R'*s_bfgs
            V_MtV( &y_bfgs_X, Q_X, BLAS_Cpp::trans, *y_bfgs );  // y_bfgs_X = Q_X'*y_bfgs
            V_MtV( &s_bfgs_X, Q_X, BLAS_Cpp::trans, *s_bfgs );  // s_bfgs_X = Q_X'*s_bfgs
            // Update rHL_RR
            if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
              out << "\nPerform current BFGS update on " << n_pz_R << " x " << n_pz_R << " projected "
                << "reduced Hessian for the superbasic variables where B = rHL_RR...\n";
            }
            QuasiNewtonStats quasi_newton_stats;
            proj_bfgs_updater().bfgs_update().perform_update(
              &s_bfgs_R(),&y_bfgs_R(),false,out,olevel,algo->algo_cntr().check_results()
              ,&rHL_RR_op, &quasi_newton_stats );
            // Perform previous updates if possible
            if( lbfgs_rHL_RR->m_bar() ) {
              const size_type num_add_updates = std::_MIN(num_add_recent_updates(),lbfgs_rHL_RR->m_bar());
              if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
                out	<< "\nAdd the min(num_previous_updates,num_add_recent_updates) = min(" << lbfgs_rHL_RR->m_bar()
                  << "," << num_add_recent_updates() << ") = " << num_add_updates << " most recent BFGS updates if possible ...\n";
              }
              // Now add previous update vectors
              const value_type
                project_error_tol = proj_bfgs_updater().project_error_tol();
              const DMatrixSlice
                S = lbfgs_rHL_RR->S(),