diff --git a/Eigen/src/Cholesky/Cholesky.h b/Eigen/src/Cholesky/Cholesky.h index db7c6e1ed..aafcd156c 100644 --- a/Eigen/src/Cholesky/Cholesky.h +++ b/Eigen/src/Cholesky/Cholesky.h @@ -54,11 +54,6 @@ template class Cholesky compute(matrix); } - Triangular > > matrixU(void) const - { - return m_matrix.transpose().temporary().upper(); - } - Triangular matrixL(void) const { return m_matrix.lower(); @@ -88,24 +83,9 @@ void Cholesky::compute(const MatrixType& matrix) { assert(matrix.rows()==matrix.cols()); const int size = matrix.rows(); - m_matrix = matrix; + m_matrix = matrix.conjugate(); - #if 1 - // this version looks faster for large matrices - m_isPositiveDefinite = m_matrix(0,0) > Scalar(0); - m_matrix(0,0) = ei_sqrt(m_matrix(0,0)); - m_matrix.col(0).end(size-1) = m_matrix.row(0).end(size-1) / m_matrix(0,0); - for (int j = 1; j < size; ++j) - { - Scalar tmp = m_matrix(j,j) - m_matrix.row(j).start(j).norm2(); - m_isPositiveDefinite = m_isPositiveDefinite && tmp > Scalar(0); - m_matrix(j,j) = ei_sqrt(tmp::compute(const MatrixType& matrix) m_matrix.col(j).end(size-j) -= m_matrix(j,i) * m_matrix.col(i).end(size-j); } } + #else + // this version looks faster for large matrices +// m_isPositiveDefinite = m_matrix(0,0) > Scalar(0); + m_matrix(0,0) = ei_sqrt(m_matrix(0,0)); + m_matrix.col(0).end(size-1) = m_matrix.row(0).end(size-1) / m_matrix(0,0); + for (int j = 1; j < size; ++j) + { +// Scalar tmp = m_matrix(j,j) - m_matrix.row(j).start(j).norm2(); + Scalar tmp = m_matrix(j,j) - (m_matrix.row(j).start(j) * m_matrix.row(j).start(j).adjoint())(0,0); +// m_isPositiveDefinite = m_isPositiveDefinite && tmp > Scalar(0); +// m_matrix(j,j) = ei_sqrt(tmp::solve(MatrixBase &ve // FIXME .inverseProduct creates a temporary that is not nice since it is called twice // add a .inverseProductInPlace ?? - return m_matrix.transpose().upper() - .inverseProduct(m_matrix.lower().inverseProduct(vecB)); + return m_matrix.adjoint().upper().inverseProduct(m_matrix.lower().inverseProduct(vecB)); } diff --git a/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h b/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h index 669728fc3..a8125b957 100644 --- a/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h +++ b/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h @@ -55,12 +55,7 @@ template class CholeskyWithoutSquareRoot compute(matrix); } - Triangular > > matrixU(void) const - { - return m_matrix.transpose().temporary().upperWithUnitDiag(); - } - - Triangular matrixL(void) const + Triangular matrixL(void) const { return m_matrix.lowerWithUnitDiag(); } @@ -78,7 +73,7 @@ template class CholeskyWithoutSquareRoot template typename DerivedVec::Eval solve(MatrixBase &vecB); - /** Compute the Cholesky decomposition A = U'DU = LDL' of \a matrix + /** Compute / recompute the Cholesky decomposition A = U'DU = LDL' of \a matrix */ void compute(const MatrixType& matrix); @@ -97,7 +92,7 @@ void CholeskyWithoutSquareRoot::compute(const MatrixType& matrix) { assert(matrix.rows()==matrix.cols()); const int size = matrix.rows(); - m_matrix = matrix; + m_matrix = matrix.conjugate(); #if 0 for (int i = 0; i < size; ++i) { @@ -113,12 +108,12 @@ void CholeskyWithoutSquareRoot::compute(const MatrixType& matrix) m_matrix.col(0).end(size-1) = m_matrix.row(0).end(size-1) / m_matrix(0,0); for (int j = 1; j < size; ++j) { - Scalar tmp = m_matrix(j,j) - (m_matrix.row(j).start(j) * m_matrix.col(j).start(j))(0,0); + Scalar tmp = m_matrix(j,j) - (m_matrix.row(j).start(j) * m_matrix.col(j).start(j).conjugate())(0,0); m_matrix(j,j) = tmp; tmp = Scalar(1) / tmp; for (int i = j+1; i < size; ++i) { - m_matrix(j,i) = (m_matrix(j,i) - (m_matrix.row(i).start(j) * m_matrix.col(j).start(j))(0,0) ); + m_matrix(j,i) = (m_matrix(j,i) - (m_matrix.row(i).start(j) * m_matrix.col(j).start(j).conjugate())(0,0) ); m_matrix(i,j) = tmp * m_matrix(j,i); } } @@ -136,12 +131,16 @@ typename DerivedVec::Eval CholeskyWithoutSquareRoot::solve(MatrixBas // FIXME .inverseProduct creates a temporary that is not nice since it is called twice // maybe add a .inverseProductInPlace() ?? - return m_matrix.transpose().upperWithUnitDiag() + return m_matrix.adjoint().upperWithUnitDiag() .inverseProduct( (m_matrix.lowerWithUnitDiag() .inverseProduct(vecB)) .cwiseQuotient(m_matrix.diagonal()) ); + +// return m_matrix.adjoint().upperWithUnitDiag() +// .inverseProduct( +// (m_matrix.lowerWithUnitDiag() * (m_matrix.diagonal().asDiagonal())).lower().inverseProduct(vecB)); } diff --git a/Eigen/src/Core/Functors.h b/Eigen/src/Core/Functors.h index bb4b8bb54..5af460dac 100644 --- a/Eigen/src/Core/Functors.h +++ b/Eigen/src/Core/Functors.h @@ -158,7 +158,8 @@ struct ei_functor_traits > * \sa class CwiseUnaryOp, MatrixBase::cwiseAbs */ template struct ei_scalar_abs_op EIGEN_EMPTY_STRUCT { - const Scalar operator() (const Scalar& a) const { return ei_abs(a); } + typedef typename NumTraits::Real result_type; + const result_type operator() (const Scalar& a) const { return ei_abs(a); } }; template struct ei_functor_traits > @@ -170,8 +171,8 @@ struct ei_functor_traits > * \sa class CwiseUnaryOp, MatrixBase::cwiseAbs2 */ template struct ei_scalar_abs2_op EIGEN_EMPTY_STRUCT { - const Scalar operator() (const Scalar& a) const { return ei_abs2(a); } - enum { Cost = NumTraits::MulCost }; + typedef typename NumTraits::Real result_type; + const result_type operator() (const Scalar& a) const { return ei_abs2(a); } }; template struct ei_functor_traits > diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index 72cf916f8..f646d64e4 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h @@ -181,6 +181,43 @@ inline std::complex ei_exp(std::complex x) { return std::exp(x) inline std::complex ei_sin(std::complex x) { return std::sin(x); } inline std::complex ei_cos(std::complex x) { return std::cos(x); } +template +inline std::complex ei_sqrt(const std::complex& x) +{ + if (std::real(x) == 0.0 && std::imag(x) == 0.0) + return std::complex(0); + else + { + T a = ei_abs(std::real(x)); + T b = ei_abs(std::imag(x)); + T c; + + if (a >= b) + { + T t = b / a; + c = ei_sqrt(a) * ei_sqrt(0.5 * (1.0 + ei_sqrt(1.0 + t * t))); + } + else + { + T t = a / b; + c = ei_sqrt(b) * ei_sqrt(0.5 * (t + ei_sqrt (1.0 + t * t))); + } + + T d = std::imag(x) / (2.0 * c); + if (std::real(x) >= 0.0) + { + return std::complex(c, d); + } + else + { + std::complex res(d, c); + if (std::imag(x)<0.0) + res = -res; + return res; + } + } +} + template<> inline std::complex ei_random() { return std::complex(ei_random(), ei_random()); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 4ff34f351..b675ed947 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -9,17 +9,17 @@ INCLUDE_DIRECTORIES( ${QT_INCLUDE_DIR} ) SET(test_SRCS cholesky.cpp main.cpp -# basicstuff.cpp -# linearstructure.cpp -# product.cpp -# adjoint.cpp -# submatrices.cpp -# miscmatrices.cpp -# smallvectors.cpp -# map.cpp -# cwiseop.cpp -# determinant.cpp -# triangular.cpp + basicstuff.cpp + linearstructure.cpp + product.cpp + adjoint.cpp + submatrices.cpp + miscmatrices.cpp + smallvectors.cpp + map.cpp + cwiseop.cpp + determinant.cpp + triangular.cpp ) QT4_AUTOMOC(${test_SRCS}) diff --git a/test/cholesky.cpp b/test/cholesky.cpp index 3d8351d8c..f4367ba6e 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -42,14 +42,15 @@ template void cholesky(const MatrixType& m) MatrixType a = MatrixType::random(rows,cols).transpose(); VectorType b = VectorType::random(cols); - SquareMatrixType covMat = a.transpose() * a; + SquareMatrixType covMat = a.adjoint() * a; CholeskyWithoutSquareRoot cholnosqrt(covMat); - VERIFY_IS_APPROX(covMat, cholnosqrt.matrixU().transpose() * cholnosqrt.vectorD().asDiagonal() * cholnosqrt.matrixU()); + VERIFY_IS_APPROX(covMat, cholnosqrt.matrixL() * cholnosqrt.vectorD().asDiagonal() * cholnosqrt.matrixL().adjoint()); VERIFY_IS_APPROX(covMat * cholnosqrt.solve(b), b); + Cholesky chol(covMat); - VERIFY_IS_APPROX(covMat, chol.matrixU().transpose() * chol.matrixU()); + VERIFY_IS_APPROX(covMat, chol.matrixL() * chol.matrixL().adjoint()); VERIFY_IS_APPROX(covMat * chol.solve(b), b); } @@ -58,7 +59,7 @@ void EigenTest::testCholesky() for(int i = 0; i < 1; i++) { cholesky(Matrix3f()); cholesky(Matrix4d()); - cholesky(MatrixXd(17,17)); + cholesky(MatrixXcd(7,7)); } } diff --git a/test/main.h b/test/main.h index 5d19c3b6b..17e7a7740 100644 --- a/test/main.h +++ b/test/main.h @@ -204,17 +204,17 @@ class EigenTest : public QObject EigenTest(int repeat) : m_repeat(repeat) {} private slots: -// void testBasicStuff(); -// void testLinearStructure(); -// void testProduct(); -// void testAdjoint(); -// void testSubmatrices(); -// void testMiscMatrices(); -// void testSmallVectors(); -// void testMap(); -// void testCwiseops(); -// void testDeterminant(); -// void testTriangular(); + void testBasicStuff(); + void testLinearStructure(); + void testProduct(); + void testAdjoint(); + void testSubmatrices(); + void testMiscMatrices(); + void testSmallVectors(); + void testMap(); + void testCwiseops(); + void testDeterminant(); + void testTriangular(); void testCholesky(); protected: int m_repeat;