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

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

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

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

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