aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Eigen/LU7
-rw-r--r--Eigen/src/Core/DummyPacketMath.h2
-rw-r--r--Eigen/src/Core/MatrixBase.h1
-rw-r--r--Eigen/src/Core/Swap.h4
-rw-r--r--Eigen/src/Core/util/ForwardDeclarations.h2
-rw-r--r--Eigen/src/Core/util/Macros.h1
-rw-r--r--Eigen/src/LU/Determinant.h4
-rw-r--r--Eigen/src/LU/LU.h188
-rw-r--r--Eigen/src/QR/QR.h2
-rw-r--r--test/CMakeLists.txt4
-rw-r--r--test/determinant.cpp72
11 files changed, 237 insertions, 50 deletions
diff --git a/Eigen/LU b/Eigen/LU
index f452d8307..c91150035 100644
--- a/Eigen/LU
+++ b/Eigen/LU
@@ -1,5 +1,5 @@
-#ifndef EIGEN_LU_H
-#define EIGEN_LU_H
+#ifndef EIGEN_LU_MODULE_H
+#define EIGEN_LU_MODULE_H
#include "Core"
@@ -16,9 +16,10 @@ namespace Eigen {
* \endcode
*/
+#include "src/LU/LU.h"
#include "src/LU/Determinant.h"
#include "src/LU/Inverse.h"
} // namespace Eigen
-#endif // EIGEN_LU_H
+#endif // EIGEN_LU_MODULE_H
diff --git a/Eigen/src/Core/DummyPacketMath.h b/Eigen/src/Core/DummyPacketMath.h
index 4abd6e997..313001f2f 100644
--- a/Eigen/src/Core/DummyPacketMath.h
+++ b/Eigen/src/Core/DummyPacketMath.h
@@ -139,7 +139,7 @@ template<int Offset,typename PacketType>
struct ei_palign_impl
{
// by default data are aligned, so there is nothing to be done :)
- inline static void run(PacketType& first, const PacketType& second) {}
+ inline static void run(PacketType&, const PacketType&) {}
};
/** \internal update \a first using the concatenation of the \a Offset last elements
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index f15e560cb..c2b12a0a3 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -526,6 +526,7 @@ template<typename Derived> class MatrixBase
/////////// LU module ///////////
+ const LU<EvalType> lu() const;
const EvalType inverse() const;
void computeInverse(EvalType *result) const;
Scalar determinant() const;
diff --git a/Eigen/src/Core/Swap.h b/Eigen/src/Core/Swap.h
index 35590b56f..3b864789e 100644
--- a/Eigen/src/Core/Swap.h
+++ b/Eigen/src/Core/Swap.h
@@ -42,7 +42,9 @@ template<typename OtherDerived>
void MatrixBase<Derived>::swap(const MatrixBase<OtherDerived>& other)
{
MatrixBase<OtherDerived> *_other = const_cast<MatrixBase<OtherDerived>*>(&other);
- if(SizeAtCompileTime == Dynamic)
+
+ // disable that path: it makes LU decomposition fail ! I can't see the bug though.
+ if(false /*SizeAtCompileTime == Dynamic*/)
{
ei_swap_selector<Derived,OtherDerived>::run(derived(),other.const_cast_derived());
}
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h
index 76a2f60cf..067ccd0b0 100644
--- a/Eigen/src/Core/util/ForwardDeclarations.h
+++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -95,7 +95,7 @@ void ei_cache_friendly_product(
bool _rhsRowMajor, const Scalar* _rhs, int _rhsStride,
bool resRowMajor, Scalar* res, int resStride);
-template<typename ExpressionType, bool CheckExistence = true> class Inverse;
+template<typename MatrixType> class LU;
template<typename MatrixType> class QR;
template<typename MatrixType> class Cholesky;
template<typename MatrixType> class CholeskyWithoutSquareRoot;
diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h
index eaa5df27f..40db9a0a7 100644
--- a/Eigen/src/Core/util/Macros.h
+++ b/Eigen/src/Core/util/Macros.h
@@ -135,7 +135,6 @@ typedef typename Eigen::NumTraits<Scalar>::Real RealScalar; \
typedef typename Base::PacketScalar PacketScalar; \
typedef typename Eigen::ei_nested<Derived>::type Nested; \
typedef typename Eigen::ei_eval<Derived>::type Eval; \
-typedef typename Eigen::Inverse<Eval> InverseType; \
enum { RowsAtCompileTime = Eigen::ei_traits<Derived>::RowsAtCompileTime, \
ColsAtCompileTime = Eigen::ei_traits<Derived>::ColsAtCompileTime, \
MaxRowsAtCompileTime = Eigen::ei_traits<Derived>::MaxRowsAtCompileTime, \
diff --git a/Eigen/src/LU/Determinant.h b/Eigen/src/LU/Determinant.h
index fdb1673d8..96e400564 100644
--- a/Eigen/src/LU/Determinant.h
+++ b/Eigen/src/LU/Determinant.h
@@ -63,7 +63,7 @@ const typename Derived::Scalar ei_bruteforce_det(const MatrixBase<Derived>& m)
- ei_bruteforce_det4_helper(m,1,3,0,2)
+ ei_bruteforce_det4_helper(m,2,3,0,1);
default:
- assert(false);
+ ei_internal_assert(false);
}
}
@@ -85,7 +85,7 @@ typename ei_traits<Derived>::Scalar MatrixBase<Derived>::determinant() const
return derived().diagonal().redux(ei_scalar_product_op<Scalar>());
}
else if(rows() <= 4) return ei_bruteforce_det(derived());
- else assert(false); // unimplemented for now
+ else return lu().determinant();
}
#endif // EIGEN_DETERMINANT_H
diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h
new file mode 100644
index 000000000..aa701f2ce
--- /dev/null
+++ b/Eigen/src/LU/LU.h
@@ -0,0 +1,188 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2006-2008 Benoit Jacob <jacob@math.jussieu.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_LU_H
+#define EIGEN_LU_H
+
+/** \ingroup LU_Module
+ *
+ * \class LU
+ *
+ * \brief LU decomposition of a matrix with complete pivoting, and associated features
+ *
+ * \param MatrixType the type of the matrix of which we are computing the LU decomposition
+ *
+ * This class performs a LU decomposition of any matrix, with complete pivoting: the matrix A
+ * is decomposed as A = PLUQ where L is unit-lower-triangular, U is upper-triangular, and P and Q
+ * are permutation matrices.
+ *
+ * This decomposition provides the generic approach to solving systems of linear equations, computing
+ * the rank, invertibility, inverse, and determinant. However for the case when invertibility is
+ * assumed, we have a specialized variant (see MatrixBase::inverse()) achieving better performance.
+ *
+ * \sa MatrixBase::lu(), MatrixBase::determinant(), MatrixBase::rank(), MatrixBase::kernelDim(),
+ * MatrixBase::kernelBasis(), MatrixBase::solve(), MatrixBase::isInvertible(),
+ * MatrixBase::inverse(), MatrixBase::computeInverse()
+ */
+template<typename MatrixType> class LU
+{
+ public:
+
+ typedef typename MatrixType::Scalar Scalar;
+ typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
+ typedef Matrix<int, MatrixType::ColsAtCompileTime, 1> IntRowVectorType;
+ typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
+
+ LU(const MatrixType& matrix);
+
+ const MatrixType& matrixLU() const
+ {
+ return m_lu;
+ }
+
+ const Part<MatrixType, UnitLower> matrixL() const
+ {
+ return m_lu;
+ }
+
+ const Part<MatrixType, Upper> matrixU() const
+ {
+ return m_lu;
+ }
+
+ const IntColVectorType& permutationP() const
+ {
+ return m_p;
+ }
+
+ const IntRowVectorType& permutationQ() const
+ {
+ return m_q;
+ }
+
+ template<typename OtherDerived>
+ typename ProductReturnType<Transpose<MatrixType>, OtherDerived>::Type::Eval
+ solve(const MatrixBase<MatrixType> &b) const;
+
+ /**
+ * This method returns the determinant of the matrix of which
+ * *this is the LU decomposition. It has only linear complexity
+ * (that is, O(n) where n is the dimension of the square matrix)
+ * as the LU decomposition has already been computed.
+ *
+ * Warning: a determinant can be very big or small, so for matrices
+ * of large dimension (like a 50-by-50 matrix) there can be a risk of
+ * overflow/underflow.
+ */
+ typename ei_traits<MatrixType>::Scalar determinant() const;
+
+ protected:
+ MatrixType m_lu;
+ IntColVectorType m_p;
+ IntRowVectorType m_q;
+ int m_det_pq;
+ Scalar m_biggest_eigenvalue_of_u;
+ int m_dimension_of_kernel;
+};
+
+template<typename MatrixType>
+LU<MatrixType>::LU(const MatrixType& matrix)
+ : m_lu(matrix),
+ m_p(matrix.rows()),
+ m_q(matrix.cols())
+{
+ const int size = matrix.diagonal().size();
+ const int rows = matrix.rows();
+ const int cols = matrix.cols();
+
+ IntColVectorType rows_transpositions(matrix.rows());
+ IntRowVectorType cols_transpositions(matrix.cols());
+ int number_of_transpositions = 0;
+
+ for(int k = 0; k < size; k++)
+ {
+ int row_of_biggest, col_of_biggest;
+ const Scalar biggest = m_lu.corner(Eigen::BottomRight, rows-k, cols-k)
+ .cwise().abs()
+ .maxCoeff(&row_of_biggest, &col_of_biggest);
+ row_of_biggest += k;
+ col_of_biggest += k;
+ rows_transpositions.coeffRef(k) = row_of_biggest;
+ cols_transpositions.coeffRef(k) = col_of_biggest;
+ if(k != row_of_biggest) {
+ m_lu.row(k).swap(m_lu.row(row_of_biggest));
+ number_of_transpositions++;
+ }
+ if(k != col_of_biggest) {
+ m_lu.col(k).swap(m_lu.col(col_of_biggest));
+ number_of_transpositions++;
+ }
+ const Scalar lu_k_k = m_lu.coeff(k,k);
+ if(ei_isMuchSmallerThan(lu_k_k, biggest)) continue;
+ if(k<rows-1)
+ m_lu.col(k).end(rows-k-1) /= lu_k_k;
+ if(k<size-1)
+ m_lu.corner(BottomRight, rows-k-1, cols-k-1)
+ -= m_lu.col(k).end(rows-k-1) * m_lu.row(k).end(cols-k-1);
+ }
+
+ for(int k = 0; k < matrix.rows(); k++) m_p.coeffRef(k) = k;
+ for(int k = size-1; k >= 0; k--)
+ std::swap(m_p.coeffRef(k), m_p.coeffRef(rows_transpositions.coeff(k)));
+
+ for(int k = 0; k < matrix.cols(); k++) m_q.coeffRef(k) = k;
+ for(int k = 0; k < size; k++)
+ std::swap(m_q.coeffRef(k), m_q.coeffRef(cols_transpositions.coeff(k)));
+
+ m_det_pq = (number_of_transpositions%2) ? -1 : 1;
+
+ int index_of_biggest;
+ m_lu.diagonal().cwise().abs().maxCoeff(&index_of_biggest);
+ m_biggest_eigenvalue_of_u = m_lu.diagonal().coeff(index_of_biggest);
+
+ m_dimension_of_kernel = 0;
+ for(int k = 0; k < size; k++)
+ m_dimension_of_kernel += ei_isMuchSmallerThan(m_lu.diagonal().coeff(k), m_biggest_eigenvalue_of_u);
+}
+
+template<typename MatrixType>
+typename ei_traits<MatrixType>::Scalar LU<MatrixType>::determinant() const
+{
+ Scalar res = m_det_pq;
+ for(int k = 0; k < m_lu.diagonal().size(); k++) res *= m_lu.diagonal().coeff(k);
+ return res;
+}
+
+/** \return the LU decomposition of \c *this.
+ *
+ * \sa class LU
+ */
+template<typename Derived>
+const LU<typename MatrixBase<Derived>::EvalType>
+MatrixBase<Derived>::lu() const
+{
+ return eval();
+}
+
+#endif // EIGEN_LU_H
diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h
index 8ca787c2a..a6cd31f12 100644
--- a/Eigen/src/QR/QR.h
+++ b/Eigen/src/QR/QR.h
@@ -171,7 +171,7 @@ template<typename Derived>
const QR<typename MatrixBase<Derived>::EvalType>
MatrixBase<Derived>::qr() const
{
- return QR<typename ei_eval<Derived>::type>(derived());
+ return eval();
}
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 97c7f4937..f3f3a625c 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -2,7 +2,7 @@ IF(BUILD_TESTS)
IF(CMAKE_COMPILER_IS_GNUCXX)
IF(CMAKE_SYSTEM_NAME MATCHES Linux)
- SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O1 -g1")
+ SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g2")
SET(CMAKE_CXX_FLAGS_RELWITHDEBINFO "${CMAKE_CXX_FLAGS_RELWITHDEBINFO} -O2 -g2")
SET(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -fno-inline-functions")
SET(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -O0 -g2")
@@ -95,7 +95,7 @@ EI_ADD_TEST(map)
EI_ADD_TEST(array)
EI_ADD_TEST(triangular)
EI_ADD_TEST(cholesky)
-# EI_ADD_TEST(determinant)
+EI_ADD_TEST(determinant)
EI_ADD_TEST(inverse)
EI_ADD_TEST(qr)
EI_ADD_TEST(eigensolver)
diff --git a/test/determinant.cpp b/test/determinant.cpp
index e19aee918..dc9f5eb54 100644
--- a/test/determinant.cpp
+++ b/test/determinant.cpp
@@ -25,54 +25,50 @@
#include "main.h"
#include <Eigen/LU>
-template<typename MatrixType> void nullDeterminant(const MatrixType& m)
+template<typename MatrixType> void determinant(const MatrixType& m)
{
/* this test covers the following files:
Determinant.h
*/
- int rows = m.rows();
- int cols = m.cols();
+ int size = m.rows();
+ MatrixType m1(size, size), m2(size, size);
+ m1.setRandom();
+ m2.setRandom();
typedef typename MatrixType::Scalar Scalar;
- typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> SquareMatrixType;
- typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
-
- MatrixType dinv(rows, cols), dnotinv(rows, cols);
-
- dinv.col(0).setOnes();
- dinv.block(0,1, rows, cols-2).setRandom();
-
- dnotinv.col(0).setOnes();
- dnotinv.block(0,1, rows, cols-2).setRandom();
- dnotinv.col(cols-1).setOnes();
-
- for (int i=0 ; i<rows ; ++i)
- {
- dnotinv.row(i).block(0,1,1,cols-2) = ei_random<Scalar>(99.999999,100.00000001)*dnotinv.row(i).block(0,1,1,cols-2).normalized();
- dnotinv(i,cols-1) = dnotinv.row(i).block(0,1,1,cols-2).norm2();
- dinv(i,cols-1) = dinv.row(i).block(0,1,1,cols-2).norm2();
- }
-
- SquareMatrixType invertibleCovarianceMatrix = dinv.transpose() * dinv;
- SquareMatrixType notInvertibleCovarianceMatrix = dnotinv.transpose() * dnotinv;
-
- std::cout << notInvertibleCovarianceMatrix << "\n" << notInvertibleCovarianceMatrix.determinant() << "\n";
-
- VERIFY_IS_MUCH_SMALLER_THAN(notInvertibleCovarianceMatrix.determinant(),
- notInvertibleCovarianceMatrix.cwise().abs().maxCoeff());
-
- VERIFY(invertibleCovarianceMatrix.inverse().exists());
-
- VERIFY(!notInvertibleCovarianceMatrix.inverse().exists());
+ Scalar x = ei_random<Scalar>();
+ VERIFY(ei_isApprox(MatrixType::Identity(size, size).determinant(), Scalar(1)));
+ VERIFY(ei_isApprox((m1*m2).determinant(), m1.determinant() * m2.determinant()));
+ if(size==1) return;
+ int i = ei_random<int>(0, size-1);
+ int j;
+ do {
+ j = ei_random<int>(0, size-1);
+ } while(j==i);
+ m2 = m1;
+ m2.row(i).swap(m2.row(j));
+ VERIFY(ei_isApprox(m2.determinant(), -m1.determinant()));
+ m2 = m1;
+ m2.col(i).swap(m2.col(j));
+ VERIFY(ei_isApprox(m2.determinant(), -m1.determinant()));
+ VERIFY(ei_isApprox(m2.determinant(), m2.transpose().determinant()));
+ VERIFY(ei_isApprox(ei_conj(m2.determinant()), m2.adjoint().determinant()));
+ m2 = m1;
+ m2.row(i) += x*m2.row(j);
+ VERIFY(ei_isApprox(m2.determinant(), m1.determinant()));
+ m2 = m1;
+ m2.row(i) *= x;
+ VERIFY(ei_isApprox(m2.determinant(), m1.determinant() * x));
}
void test_determinant()
{
for(int i = 0; i < g_repeat; i++) {
- CALL_SUBTEST( nullDeterminant(Matrix<float, 30, 3>()) );
- CALL_SUBTEST( nullDeterminant(Matrix<double, 30, 3>()) );
- CALL_SUBTEST( nullDeterminant(Matrix<float, 20, 4>()) );
- CALL_SUBTEST( nullDeterminant(Matrix<double, 20, 4>()) );
-// CALL_SUBTEST( nullDeterminant(MatrixXd(20,4));
+ CALL_SUBTEST( determinant(Matrix<float, 1, 1>()) );
+ CALL_SUBTEST( determinant(Matrix<double, 2, 2>()) );
+ CALL_SUBTEST( determinant(Matrix<double, 3, 3>()) );
+ CALL_SUBTEST( determinant(Matrix<double, 4, 4>()) );
+ CALL_SUBTEST( determinant(Matrix<std::complex<double>, 10, 10>()) );
+ CALL_SUBTEST( determinant(MatrixXd(20, 20)) );
}
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-05 00:20:44 +0000