diff options
author | Gael Guennebaud <g.gael@free.fr> | 2012-02-28 16:19:40 +0100 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2012-02-28 16:19:40 +0100 |
commit | fc85f91df06699cbc643409df4c3aacaabc6f484 (patch) | |
tree | bd830004cf04fde5732a8ac0528ed5720b922737 /Eigen/src/Cholesky | |
parent | 309b27b545a79a9fefd05a2d7927553e950418b3 (diff) |
fix MKL interface with LLT::rankUpdate
Diffstat (limited to 'Eigen/src/Cholesky')
-rw-r--r-- | Eigen/src/Cholesky/LLT.h | 142 | ||||
-rw-r--r-- | Eigen/src/Cholesky/LLT_MKL.h | 35 |
2 files changed, 80 insertions, 97 deletions
diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index 0ce496c35..dd705c57e 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -184,7 +184,7 @@ template<typename _MatrixType, int _UpLo> class LLT inline Index cols() const { return m_matrix.cols(); } template<typename VectorType> - LLT& rankUpdate(const VectorType& vec, const RealScalar& sigma = 1); + LLT rankUpdate(const VectorType& vec, const RealScalar& sigma = 1); protected: /** \internal @@ -200,6 +200,76 @@ namespace internal { 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> { typedef typename NumTraits<Scalar>::Real RealScalar; @@ -265,72 +335,10 @@ template<typename Scalar> struct llt_inplace<Scalar, Lower> template<typename MatrixType, typename VectorType> static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma) { - 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; + return llt_rank_update_lower(mat, vec, sigma); } }; - + template<typename Scalar> struct llt_inplace<Scalar, Upper> { 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 * of same dimension. */ -template<typename MatrixType, int _UpLo> +template<typename _MatrixType, int _UpLo> 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_assert(v.size()==m_matrix.cols()); diff --git a/Eigen/src/Cholesky/LLT_MKL.h b/Eigen/src/Cholesky/LLT_MKL.h index c5ed77d33..10aa746dc 100644 --- a/Eigen/src/Cholesky/LLT_MKL.h +++ b/Eigen/src/Cholesky/LLT_MKL.h @@ -50,7 +50,7 @@ template<> struct mkl_llt<EIGTYPE> \ lapack_int size, lda, info, StorageOrder; \ EIGTYPE* a; \ eigen_assert(m.rows()==m.cols()); \ -/* Set up parameters for ?potrf */ \ + /* Set up parameters for ?potrf */ \ size = m.rows(); \ StorageOrder = MatrixType::Flags&RowMajorBit?RowMajor:ColMajor; \ 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'); \ } \ template<typename MatrixType, typename VectorType> \ - static void rankUpdate(MatrixType& mat, const VectorType& vec) \ - { \ - 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); \ - } \ - } \ - } \ + static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma) \ + { return llt_rank_update_lower(mat, vec, sigma); } \ }; \ template<> struct llt_inplace<EIGTYPE, Upper> \ { \ @@ -106,10 +81,10 @@ template<> struct llt_inplace<EIGTYPE, Upper> \ return mkl_llt<EIGTYPE>::potrf(m, 'U'); \ } \ 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); \ - return llt_inplace<EIGTYPE, Lower>::rankUpdate(matt, vec.conjugate()); \ + return llt_inplace<EIGTYPE, Lower>::rankUpdate(matt, vec.conjugate(), sigma); \ } \ }; |