GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-12 19:59:05 +08:00
Code Issues Packages Projects Releases Wiki Activity

fix MKL interface with LLT::rankUpdate

Browse Source
This commit is contained in:
Gael Guennebaud 2012-02-28 16:19:40 +01:00
parent 309b27b545
commit fc85f91df0
2 changed files with 80 additions and 97 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

142
Eigen/src/Cholesky/LLT.h
View File

@ -184,7 +184,7 @@ template<typename _MatrixType, int _UpLo> class LLT
inline Index cols() const { return m_matrix.cols(); } inline Index cols() const { return m_matrix.cols(); }
template<typename VectorType> template<typename VectorType>
LLT& rankUpdate(const VectorType& vec, const RealScalar& sigma = 1); LLT rankUpdate(const VectorType& vec, const RealScalar& sigma = 1);
protected: protected:
/** \internal /** \internal
@ -200,6 +200,76 @@ namespace internal {
template<typename Scalar, int UpLo> struct llt_inplace; template<typename Scalar, int UpLo> struct llt_inplace;
template<typename MatrixType, typename VectorType>
static typename MatrixType::Index llt_rank_update_lower(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma)
{
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef typename MatrixType::Index Index;
typedef typename MatrixType::ColXpr ColXpr;
typedef typename internal::remove_all<ColXpr>::type ColXprCleaned;
typedef typename ColXprCleaned::SegmentReturnType ColXprSegment;
typedef Matrix<Scalar,Dynamic,1> TempVectorType;
typedef typename TempVectorType::SegmentReturnType TempVecSegment;
int n = mat.cols();
eigen_assert(mat.rows()==n && vec.size()==n);
TempVectorType temp;
if(sigma>0)
{
// This version is based on Givens rotations.
// It is faster than the other one below, but only works for updates,
// i.e., for sigma > 0
temp = sqrt(sigma) * vec;
for(int i=0; i<n; ++i)
{
JacobiRotation<Scalar> g;
g.makeGivens(mat(i,i), -temp(i), &mat(i,i));
int rs = n-i-1;
if(rs>0)
{
ColXprSegment x(mat.col(i).tail(rs));
TempVecSegment y(temp.tail(rs));
apply_rotation_in_the_plane(x, y, g);
}
}
}
else
{
temp = vec;
RealScalar beta = 1;
for(int j=0; j<n; ++j)
{
RealScalar Ljj = real(mat.coeff(j,j));
RealScalar dj = abs2(Ljj);
Scalar wj = temp.coeff(j);
RealScalar swj2 = sigma*abs2(wj);
RealScalar gamma = dj*beta + swj2;
RealScalar x = dj + swj2/beta;
if (x<=RealScalar(0))
return j;
RealScalar nLjj = sqrt(x);
mat.coeffRef(j,j) = nLjj;
beta += swj2/dj;
// Update the terms of L
Index rs = n-j-1;
if(rs)
{
temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
if(gamma != 0)
mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*conj(wj)/gamma)*temp.tail(rs);
}
}
}
return -1;
}
template<typename Scalar> struct llt_inplace<Scalar, Lower> template<typename Scalar> struct llt_inplace<Scalar, Lower>
{ {
typedef typename NumTraits<Scalar>::Real RealScalar; typedef typename NumTraits<Scalar>::Real RealScalar;
@ -265,72 +335,10 @@ template<typename Scalar> struct llt_inplace<Scalar, Lower>
template<typename MatrixType, typename VectorType> template<typename MatrixType, typename VectorType>
static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma) static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
{ {
typedef typename MatrixType::Index Index; return llt_rank_update_lower(mat, vec, sigma);
typedef typename MatrixType::ColXpr ColXpr;
typedef typename internal::remove_all<ColXpr>::type ColXprCleaned;
typedef typename ColXprCleaned::SegmentReturnType ColXprSegment;
typedef Matrix<Scalar,Dynamic,1> TempVectorType;
typedef typename TempVectorType::SegmentReturnType TempVecSegment;
int n = mat.cols();
eigen_assert(mat.rows()==n && vec.size()==n);
TempVectorType temp;
if(sigma>0)
{
// This version is based on Givens rotations.
// It is faster than the other one below, but only works for updates,
// i.e., for sigma > 0
temp = sqrt(sigma) * vec;
for(int i=0; i<n; ++i)
{
JacobiRotation<Scalar> g;
g.makeGivens(mat(i,i), -temp(i), &mat(i,i));
int rs = n-i-1;
if(rs>0)
{
ColXprSegment x(mat.col(i).tail(rs));
TempVecSegment y(temp.tail(rs));
apply_rotation_in_the_plane(x, y, g);
}
}
}
else
{
temp = vec;
RealScalar beta = 1;
for(int j=0; j<n; ++j)
{
RealScalar Ljj = real(mat.coeff(j,j));
RealScalar dj = abs2(Ljj);
Scalar wj = temp.coeff(j);
RealScalar swj2 = sigma*abs2(wj);
RealScalar gamma = dj*beta + swj2;
RealScalar x = dj + swj2/beta;
if (x<=RealScalar(0))
return j;
RealScalar nLjj = sqrt(x);
mat.coeffRef(j,j) = nLjj;
beta += swj2/dj;
// Update the terms of L
Index rs = n-j-1;
if(rs)
{
temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
if(gamma != 0)
mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*conj(wj)/gamma)*temp.tail(rs);
}
}
}
return -1;
} }
}; };
template<typename Scalar> struct llt_inplace<Scalar, Upper> template<typename Scalar> struct llt_inplace<Scalar, Upper>
{ {
typedef typename NumTraits<Scalar>::Real RealScalar; typedef typename NumTraits<Scalar>::Real RealScalar;
@ -404,9 +412,9 @@ LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
* then after it we have LL^* = A + sigma * v v^* where \a v must be a vector * then after it we have LL^* = A + sigma * v v^* where \a v must be a vector
* of same dimension. * of same dimension.
*/ */
template<typename MatrixType, int _UpLo> template<typename _MatrixType, int _UpLo>
template<typename VectorType> template<typename VectorType>
LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma) LLT<_MatrixType,_UpLo> LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma)
{ {
EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType); EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType);
eigen_assert(v.size()==m_matrix.cols()); eigen_assert(v.size()==m_matrix.cols());

35
Eigen/src/Cholesky/LLT_MKL.h
View File

@ -50,7 +50,7 @@ template<> struct mkl_llt<EIGTYPE> \
lapack_int size, lda, info, StorageOrder; \ lapack_int size, lda, info, StorageOrder; \
EIGTYPE* a; \ EIGTYPE* a; \
eigen_assert(m.rows()==m.cols()); \ eigen_assert(m.rows()==m.cols()); \
/* Set up parameters for ?potrf */ \ /* Set up parameters for ?potrf */ \
size = m.rows(); \ size = m.rows(); \
StorageOrder = MatrixType::Flags&RowMajorBit?RowMajor:ColMajor; \ StorageOrder = MatrixType::Flags&RowMajorBit?RowMajor:ColMajor; \
matrix_order = StorageOrder==RowMajor ? LAPACK_ROW_MAJOR : LAPACK_COL_MAJOR; \ matrix_order = StorageOrder==RowMajor ? LAPACK_ROW_MAJOR : LAPACK_COL_MAJOR; \
@ -70,33 +70,8 @@ template<> struct llt_inplace<EIGTYPE, Lower> \
return mkl_llt<EIGTYPE>::potrf(m, 'L'); \ return mkl_llt<EIGTYPE>::potrf(m, 'L'); \
} \ } \
template<typename MatrixType, typename VectorType> \ template<typename MatrixType, typename VectorType> \
static void rankUpdate(MatrixType& mat, const VectorType& vec) \ static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma) \
{ \ { return llt_rank_update_lower(mat, vec, sigma); } \
typedef typename MatrixType::ColXpr ColXpr; \
typedef typename internal::remove_all<ColXpr>::type ColXprCleaned; \
typedef typename ColXprCleaned::SegmentReturnType ColXprSegment; \
typedef typename MatrixType::Scalar Scalar; \
typedef Matrix<Scalar,Dynamic,1> TempVectorType; \
typedef typename TempVectorType::SegmentReturnType TempVecSegment; \
\
int n = mat.cols(); \
eigen_assert(mat.rows()==n && vec.size()==n); \
TempVectorType temp(vec); \
\
for(int i=0; i<n; ++i) \
{ \
JacobiRotation<Scalar> g; \
g.makeGivens(mat(i,i), -temp(i), &mat(i,i)); \
\
int rs = n-i-1; \
if(rs>0) \
{ \
ColXprSegment x(mat.col(i).tail(rs)); \
TempVecSegment y(temp.tail(rs)); \
apply_rotation_in_the_plane(x, y, g); \
} \
} \
} \
}; \ }; \
template<> struct llt_inplace<EIGTYPE, Upper> \ template<> struct llt_inplace<EIGTYPE, Upper> \
{ \ { \
@ -106,10 +81,10 @@ template<> struct llt_inplace<EIGTYPE, Upper> \
return mkl_llt<EIGTYPE>::potrf(m, 'U'); \ return mkl_llt<EIGTYPE>::potrf(m, 'U'); \
} \ } \
template<typename MatrixType, typename VectorType> \ template<typename MatrixType, typename VectorType> \
static void rankUpdate(MatrixType& mat, const VectorType& vec) \ static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma) \
{ \ { \
Transpose<MatrixType> matt(mat); \ Transpose<MatrixType> matt(mat); \
return llt_inplace<EIGTYPE, Lower>::rankUpdate(matt, vec.conjugate()); \ return llt_inplace<EIGTYPE, Lower>::rankUpdate(matt, vec.conjugate(), sigma); \
} \ } \
}; };