C++ DMatrixSlice类代码示例

bool DenseLinAlgPack::comp(const DMatrixSlice& gms1, BLAS_Cpp::Transp trans1
  , const DMatrixSlice& gms2, BLAS_Cpp::Transp trans2)
{
  for(size_type i = 1; i < my_min(gms1.cols(),gms2.cols()); ++i)
    if( !comp( col(gms1,trans1,i) , col( gms2, trans2, i ) ) ) return false;
  return true;
}

void DenseLinAlgPack::M_diagVtM( DMatrixSlice* gms_lhs, const DVectorSlice& vs_rhs
                , const DMatrixSlice& gms_rhs, BLAS_Cpp::Transp trans_rhs )
{
  Mp_M_assert_sizes(gms_lhs->rows(), gms_lhs->cols(), BLAS_Cpp::no_trans
    , gms_rhs.rows(), gms_rhs.cols(), trans_rhs);
  for(DMatrixSlice::size_type j = 1; j <= gms_lhs->cols(); ++j)
    prod( &gms_lhs->col(j), vs_rhs, col(gms_rhs,trans_rhs,j) );
}

inline 
/// Assert a matrix is square and throws length_error if it is not (LINALGPACK_CHECK_SLICE_SETUP).
void assert_gms_square(const DMatrixSlice& gms) {
#ifdef LINALGPACK_CHECK_SLICE_SETUP
  if(gms.rows() != gms.cols())
    throw std::length_error("Matrix must be square");
#endif
} 

void DenseLinAlgPack::M_StMtInvM(DMatrixSlice* gms_lhs, value_type alpha, const DMatrixSlice& gms_rhs1
  , BLAS_Cpp::Transp trans_rhs1, const DMatrixSliceTri& tri_rhs2
  , BLAS_Cpp::Transp trans_rhs2)
{
  Mp_MtM_assert_sizes(	  gms_lhs->rows(), gms_lhs->cols(), BLAS_Cpp::no_trans
              , gms_rhs1.rows(), gms_rhs1.cols(), trans_rhs1
              , tri_rhs2.gms().rows(), tri_rhs2.gms().cols(), trans_rhs2 );
  i_trsm_alt(BLAS_Cpp::right,alpha,tri_rhs2,trans_rhs2,gms_rhs1,trans_rhs1,gms_lhs);
}

void DenseLinAlgPack::M_StInvMtM(DMatrix* gm_lhs, value_type alpha, const DMatrixSliceTri& tri_rhs1
  , BLAS_Cpp::Transp trans_rhs1, const DMatrixSlice& gms_rhs2
  , BLAS_Cpp::Transp trans_rhs2)
{
  MtM_assert_sizes(  tri_rhs1.gms().rows(), tri_rhs1.gms().cols(), trans_rhs1
           , gms_rhs2.rows(), gms_rhs2.cols(), trans_rhs2					);
  gm_lhs->resize(	  rows(tri_rhs1.gms().rows(), tri_rhs1.gms().cols(), trans_rhs1)
          , cols(gms_rhs2.rows(), gms_rhs2.cols(), trans_rhs2)				);
  i_trsm_alt(BLAS_Cpp::left,alpha,tri_rhs1,trans_rhs1,gms_rhs2,trans_rhs2,&(*gm_lhs)());
}

inline
DMatrixSlice::DMatrixSlice( DMatrixSlice& gms, const Range1D& I
    , const Range1D& J)
  : ptr_( gms.col_ptr(1) + (I.lbound() - 1) + (J.lbound() - 1) * gms.max_rows() )
  , max_rows_(gms.max_rows())
  , rows_(I.size())
  , cols_(J.size())
{	
  gms.validate_row_subscript(I.ubound());
  gms.validate_col_subscript(J.ubound());
}

void DenseLinAlgPack::syr2k(BLAS_Cpp::Transp trans,value_type alpha, const DMatrixSlice& gms_rhs1
  , const DMatrixSlice& gms_rhs2, value_type beta, DMatrixSliceSym* sym_lhs)
{
  Mp_MtM_assert_sizes(	  sym_lhs->gms().rows(), sym_lhs->gms().cols(), BLAS_Cpp::no_trans
              , gms_rhs1.rows(), gms_rhs1.cols(), trans
              , gms_rhs2.rows(), gms_rhs2.cols(), trans_not(trans) );
  BLAS_Cpp::syr2k(sym_lhs->uplo(),trans,sym_lhs->gms().rows()
    ,cols(gms_rhs1.rows(), gms_rhs1.cols(), trans),alpha
    ,gms_rhs1.col_ptr(1),gms_rhs1.max_rows()
    ,gms_rhs2.col_ptr(1),gms_rhs2.max_rows(),beta
    ,sym_lhs->gms().col_ptr(1),sym_lhs->gms().max_rows() );
}

void DenseLinAlgPack::Mp_StMtM(DMatrixSlice* gms_lhs, value_type alpha, const DMatrixSlice& gms_rhs1
  , BLAS_Cpp::Transp trans_rhs1, const DMatrixSlice& gms_rhs2
  , BLAS_Cpp::Transp trans_rhs2, value_type beta)
{
  Mp_MtM_assert_sizes(	  gms_lhs->rows(), gms_lhs->cols(), BLAS_Cpp::no_trans
              , gms_rhs1.rows(), gms_rhs1.cols(), trans_rhs1
              , gms_rhs2.rows(), gms_rhs2.cols(), trans_rhs2);
  BLAS_Cpp::gemm(trans_rhs1,trans_rhs2,gms_lhs->rows(),gms_lhs->cols()
    ,cols(gms_rhs1.rows(),gms_rhs1.cols(),trans_rhs1)
    ,alpha,gms_rhs1.col_ptr(1),gms_rhs1.max_rows()
    ,gms_rhs2.col_ptr(1),gms_rhs2.max_rows()
    ,beta,gms_lhs->col_ptr(1),gms_lhs->max_rows() );
}

// gm_lhs = op(gms_rhs)
void DenseLinAlgPack::assign(DMatrix* gm_lhs, const DMatrixSlice& gms_rhs, BLAS_Cpp::Transp trans_rhs)
{
  if(gm_lhs->overlap(gms_rhs) == SAME_MEM && trans_rhs == BLAS_Cpp::no_trans) return;	// assignment to self
  if(gm_lhs->overlap(gms_rhs) != NO_OVERLAP) {
    // some overlap so we must create a copy
    DMatrix tmp(gms_rhs);
    resize_gm_lhs(gm_lhs,gms_rhs.rows(),gms_rhs.cols(),trans_rhs);
    i_assign(&(*gm_lhs)(), tmp(), trans_rhs);
  }
  else {
    // no overlap so just assign
    resize_gm_lhs(gm_lhs,gms_rhs.rows(),gms_rhs.cols(),trans_rhs);
    i_assign(&(*gm_lhs)(), gms_rhs, trans_rhs);
  }
}

inline 
/** \brief . */
/* * Utility to check if a lhs matrix slice is the same size as a rhs matrix slice.
  *
  * A DMatrixSlice can not be resized since the rows_ property of the
  * DMatrix it came from will not be updated.  Allowing a DMatrixSlice
  * to resize from unsized would require that the DMatrixSlice carry
  * a reference to the DMatrix it was created from.  If this is needed
  * then it will be added.
  */
void assert_gms_lhs(const DMatrixSlice& gms_lhs, size_type rows, size_type cols
  , BLAS_Cpp::Transp trans_rhs = BLAS_Cpp::no_trans)
{
  if(trans_rhs == BLAS_Cpp::trans) std::swap(rows,cols);
  if(gms_lhs.rows() == rows && gms_lhs.cols() == cols) return; // same size
  // not the same size so is an error
  throw std::length_error("assert_gms_lhs(...):  lhs DMatrixSlice dim does not match rhs dim");
}

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::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::Vp_StMtV(DVectorSlice* vs_lhs, value_type alpha, const DMatrixSlice& gms_rhs1
  , BLAS_Cpp::Transp trans_rhs1, const DVectorSlice& vs_rhs2, value_type beta)
{
  Vp_MtV_assert_sizes(vs_lhs->dim(), gms_rhs1.rows()	, gms_rhs1.cols(), trans_rhs1
    , vs_rhs2.dim());
  BLAS_Cpp::gemv(trans_rhs1,gms_rhs1.rows(),gms_rhs1.cols(),alpha,gms_rhs1.col_ptr(1)
    ,gms_rhs1.max_rows(), vs_rhs2.raw_ptr(),vs_rhs2.stride(),beta,vs_lhs->raw_ptr()
    ,vs_lhs->stride());
}

void MatrixSymPosDefLBFGS::update_Q() const
{
  using DenseLinAlgPack::tri;
  using DenseLinAlgPack::tri_ele;
  using DenseLinAlgPack::Mp_StM;

  //
  // We need update the factorizations to solve for:
  //
  // x = inv(Q) * y
  //
  //	[ y1 ]	=	[ (1/gk)*S'S	 L	] * [ x1 ]
  //	[ y2 ]		[      L'		-D	]   [ x2 ]
  //
  // We will solve the above system using the schur complement:
  //
  // C = (1/gk)*S'S + L*inv(D)*L'
  //
  // According to the referenced paper, C is p.d. so:
  //
  // C = J*J'
  //
  // We then compute the solution as:
  //
  // x1 = inv(C) * ( y1 + L*inv(D)*y2 )
  // x2 = - inv(D) * ( y2 - L'*x1 )
  //
  // Therefore we will just update the factorization C = J*J'
  //

  // Form the upper triangular part of C which will become J
  // which we are using storage of QJ

  if( QJ_.rows() < m_ )
    QJ_.resize( m_, m_ );

  const size_type
    mb = m_bar_;

  DMatrixSlice
    C = QJ_(1,mb,1,mb);

  // C = L * inv(D) * L'
  //
  // Here L is a strictly lower triangular (zero diagonal) matrix where:
  //
  // L = [ 0  0 ]
  //     [ Lb 0 ]
  //
  // Lb is lower triangular (nonzero diagonal)
  //
  // Therefore we compute L*inv(D)*L' as:
  //
  // C = [ 0	0 ] * [ Db  0  ] * [ 0  Lb' ]
  //	   [ Lb 0 ]   [ 0   db ]   [ 0   0  ]
  //
  //   = [ 0  0  ] = [ 0      0     ]
  //     [ 0  Cb ]   [ 0  Lb*Db*Lb' ]
  //
  // We need to compute the upper triangular part of Cb.

  C.row(1) = 0.0;
  if( mb > 1 )
    comp_Cb( STY_(2,mb,1,mb-1), STY_.diag(0)(1,mb-1), &C(2,mb,2,mb) );

  // C += (1/gk)*S'S

  const DMatrixSliceSym &STS = this->STS();
  Mp_StM( &C, (1/gamma_k_), tri( STS.gms(), STS.uplo(), BLAS_Cpp::nonunit )
    , BLAS_Cpp::trans );

  // Now perform a cholesky factorization of C
  // After this factorization the upper triangular part of QJ
  // (through C) will contain the cholesky factor.

  DMatrixSliceTriEle C_upper = tri_ele( C, BLAS_Cpp::upper );
  try {
    DenseLinAlgLAPack::potrf( &C_upper );
  }
  catch( const DenseLinAlgLAPack::FactorizationException &fe ) {
    TEUCHOS_TEST_FOR_EXCEPTION(
      true, UpdateFailedException
      ,"Error, the factorization of Q which should be s.p.d. failed with"
      " the error message: {" << fe.what() << "}";
      );
  }

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