diff options
-rw-r--r-- | src/Core/Conjugate.h | 2 | ||||
-rw-r--r-- | src/Core/Dot.h | 12 | ||||
-rw-r--r-- | src/Core/Fuzzy.h | 4 | ||||
-rw-r--r-- | src/Core/MatrixBase.h | 12 | ||||
-rw-r--r-- | src/Core/NumTraits.h | 252 | ||||
-rw-r--r-- | src/Core/Random.h | 2 | ||||
-rw-r--r-- | test/adjoint.cpp | 32 | ||||
-rw-r--r-- | test/basicstuff.cpp | 13 | ||||
-rw-r--r-- | test/main.cpp | 7 | ||||
-rw-r--r-- | test/main.h | 6 |
10 files changed, 193 insertions, 149 deletions
diff --git a/src/Core/Conjugate.h b/src/Core/Conjugate.h index 1d763ed1b..dc5eb4653 100644 --- a/src/Core/Conjugate.h +++ b/src/Core/Conjugate.h @@ -52,7 +52,7 @@ template<typename MatrixType> class Conjugate Scalar _read(int row, int col) const { - return NumTraits<Scalar>::conj(m_matrix.read(row, col)); + return conj(m_matrix.read(row, col)); } protected: diff --git a/src/Core/Dot.h b/src/Core/Dot.h index 4667d8129..4a0629e56 100644 --- a/src/Core/Dot.h +++ b/src/Core/Dot.h @@ -32,7 +32,7 @@ struct DotUnroller static void run(const Derived1 &v1, const Derived2& v2, typename Derived1::Scalar &dot) { DotUnroller<Index-1, Size, Derived1, Derived2>::run(v1, v2, dot); - dot += v1[Index] * NumTraits<typename Derived1::Scalar>::conj(v2[Index]); + dot += v1[Index] * conj(v2[Index]); } }; @@ -41,7 +41,7 @@ struct DotUnroller<0, Size, Derived1, Derived2> { static void run(const Derived1 &v1, const Derived2& v2, typename Derived1::Scalar &dot) { - dot = v1[0] * NumTraits<typename Derived1::Scalar>::conj(v2[0]); + dot = v1[0] * conj(v2[0]); } }; @@ -67,9 +67,9 @@ Scalar MatrixBase<Scalar, Derived>::dot(const OtherDerived& other) const ::run(*static_cast<const Derived*>(this), other, res); else { - res = (*this)[0] * NumTraits<Scalar>::conj(other[0]); + res = (*this)[0] * conj(other[0]); for(int i = 1; i < size(); i++) - res += (*this)[i]* NumTraits<Scalar>::conj(other[i]); + res += (*this)[i]* conj(other[i]); } return res; } @@ -77,13 +77,13 @@ Scalar MatrixBase<Scalar, Derived>::dot(const OtherDerived& other) const template<typename Scalar, typename Derived> typename NumTraits<Scalar>::Real MatrixBase<Scalar, Derived>::norm2() const { - return NumTraits<Scalar>::real(dot(*this)); + return real(dot(*this)); } template<typename Scalar, typename Derived> typename NumTraits<Scalar>::Real MatrixBase<Scalar, Derived>::norm() const { - return std::sqrt(norm2()); + return sqrt(norm2()); } template<typename Scalar, typename Derived> diff --git a/src/Core/Fuzzy.h b/src/Core/Fuzzy.h index 1939217bf..8aa6dd8b2 100644 --- a/src/Core/Fuzzy.h +++ b/src/Core/Fuzzy.h @@ -56,12 +56,12 @@ bool MatrixBase<Scalar, Derived>::isMuchSmallerThan( { if(IsVector) { - return(norm2() <= NumTraits<Scalar>::abs2(other) * prec * prec); + return(norm2() <= abs2(other) * prec * prec); } else { for(int i = 0; i < cols(); i++) - if(col(i).norm2() > NumTraits<Scalar>::abs2(other) * prec * prec) + if(col(i).norm2() > abs2(other) * prec * prec) return false; return true; } diff --git a/src/Core/MatrixBase.h b/src/Core/MatrixBase.h index 7f1a34b44..f9a668b06 100644 --- a/src/Core/MatrixBase.h +++ b/src/Core/MatrixBase.h @@ -37,16 +37,16 @@ template<typename Scalar, typename Derived> class MatrixBase template<typename OtherDerived> bool _isApprox_helper( const OtherDerived& other, - const typename NumTraits<Scalar>::Real& prec = NumTraits<Scalar>::precision() + const typename NumTraits<Scalar>::Real& prec = precision<Scalar>() ) const; bool _isMuchSmallerThan_helper( const Scalar& other, - const typename NumTraits<Scalar>::Real& prec = NumTraits<Scalar>::precision() + const typename NumTraits<Scalar>::Real& prec = precision<Scalar>() ) const; template<typename OtherDerived> bool _isMuchSmallerThan_helper( const MatrixBase<Scalar, OtherDerived>& other, - const typename NumTraits<Scalar>::Real& prec = NumTraits<Scalar>::precision() + const typename NumTraits<Scalar>::Real& prec = precision<Scalar>() ) const; public: @@ -121,16 +121,16 @@ template<typename Scalar, typename Derived> class MatrixBase template<typename OtherDerived> bool isApprox( const OtherDerived& other, - const typename NumTraits<Scalar>::Real& prec = NumTraits<Scalar>::precision() + const typename NumTraits<Scalar>::Real& prec = precision<Scalar>() ) const; bool isMuchSmallerThan( const Scalar& other, - const typename NumTraits<Scalar>::Real& prec = NumTraits<Scalar>::precision() + const typename NumTraits<Scalar>::Real& prec = precision<Scalar>() ) const; template<typename OtherDerived> bool isMuchSmallerThan( const MatrixBase<Scalar, OtherDerived>& other, - const typename NumTraits<Scalar>::Real& prec = NumTraits<Scalar>::precision() + const typename NumTraits<Scalar>::Real& prec = precision<Scalar>() ) const; template<typename OtherDerived> diff --git a/src/Core/NumTraits.h b/src/Core/NumTraits.h index 132c9205e..d4563c8e0 100644 --- a/src/Core/NumTraits.h +++ b/src/Core/NumTraits.h @@ -23,144 +23,176 @@ // License. This exception does not invalidate any other reasons why a work // based on this file might be covered by the GNU General Public License. -#ifndef EIGEN_NUMERIC_H -#define EIGEN_NUMERIC_H +#ifndef EIGEN_NUMTRAITS_H +#define EIGEN_NUMTRAITS_H template<typename T> struct NumTraits; template<> struct NumTraits<int> { typedef int Real; - typedef double FloatingPoint; - typedef double RealFloatingPoint; - static const bool IsComplex = false; static const bool HasFloatingPoint = false; - - static int precision() { return 0; } - static int real(const int& x) { return x; } - static int imag(const int& x) { EIGEN_UNUSED(x); return 0; } - static int conj(const int& x) { return x; } - static int abs2(const int& x) { return x*x; } - static int random() - { - // "random() % n" is bad, they say, because the low-order bits are not random enough. - // However here, 21 is odd, so random() % 21 uses the high-order bits - // as well, so there's no problem. - return (std::rand() % 21) - 10; - } - static bool isMuchSmallerThan(const int& a, const int& b, const int& prec = precision()) - { - EIGEN_UNUSED(b); - EIGEN_UNUSED(prec); - return a == 0; - } - static bool isApprox(const int& a, const int& b, const int& prec = precision()) - { - EIGEN_UNUSED(prec); - return a == b; - } - static bool isApproxOrLessThan(const int& a, const int& b, const int& prec = precision()) - { - EIGEN_UNUSED(prec); - return a <= b; - } }; template<> struct NumTraits<float> { typedef float Real; - typedef float FloatingPoint; - typedef float RealFloatingPoint; - static const bool IsComplex = false; static const bool HasFloatingPoint = true; - - static float precision() { return 1e-5f; } - static float real(const float& x) { return x; } - static float imag(const float& x) { EIGEN_UNUSED(x); return 0; } - static float conj(const float& x) { return x; } - static float abs2(const float& x) { return x*x; } - static float random() - { - return std::rand() / (RAND_MAX/20.0f) - 10.0f; - } - static bool isMuchSmallerThan(const float& a, const float& b, const float& prec = precision()) - { - return std::abs(a) <= std::abs(b) * prec; - } - static bool isApprox(const float& a, const float& b, const float& prec = precision()) - { - return std::abs(a - b) <= std::min(std::abs(a), std::abs(b)) * prec; - } - static bool isApproxOrLessThan(const float& a, const float& b, const float& prec = precision()) - { - return a <= b || isApprox(a, b, prec); - } }; template<> struct NumTraits<double> { typedef double Real; - typedef double FloatingPoint; - typedef double RealFloatingPoint; - static const bool IsComplex = false; static const bool HasFloatingPoint = true; - - static double precision() { return 1e-11; } - static double real(const double& x) { return x; } - static double imag(const double& x) { EIGEN_UNUSED(x); return 0; } - static double conj(const double& x) { return x; } - static double abs2(const double& x) { return x*x; } - static double random() - { - return std::rand() / (RAND_MAX/20.0) - 10.0; - } - static bool isMuchSmallerThan(const double& a, const double& b, const double& prec = precision()) - { - return std::abs(a) <= std::abs(b) * prec; - } - static bool isApprox(const double& a, const double& b, const double& prec = precision()) - { - return std::abs(a - b) <= std::min(std::abs(a), std::abs(b)) * prec; - } - static bool isApproxOrLessThan(const double& a, const double& b, const double& prec = precision()) - { - return a <= b || isApprox(a, b, prec); - } }; template<typename _Real> struct NumTraits<std::complex<_Real> > { typedef _Real Real; - typedef std::complex<Real> Complex; - typedef std::complex<double> FloatingPoint; - typedef typename NumTraits<Real>::FloatingPoint RealFloatingPoint; - static const bool IsComplex = true; static const bool HasFloatingPoint = NumTraits<Real>::HasFloatingPoint; - - static Real precision() { return NumTraits<Real>::precision(); } - static Real real(const Complex& x) { return std::real(x); } - static Real imag(const Complex& x) { return std::imag(x); } - static Complex conj(const Complex& x) { return std::conj(x); } - static Real abs2(const Complex& x) - { return std::norm(x); } - static Complex random() - { - return Complex(NumTraits<Real>::random(), NumTraits<Real>::random()); - } - static bool isMuchSmallerThan(const Complex& a, const Complex& b, const Real& prec = precision()) - { - return abs2(a) <= abs2(b) * prec * prec; - } - static bool isApprox(const Complex& a, const Complex& b, const Real& prec = precision()) - { - return NumTraits<Real>::isApprox(std::real(a), std::real(b), prec) - && NumTraits<Real>::isApprox(std::imag(a), std::imag(b), prec); - } - // isApproxOrLessThan wouldn't make sense for complex numbers }; -#endif // EIGEN_NUMERIC_H +template<typename T> inline typename NumTraits<T>::Real precision(); +template<typename T> inline T random(); + +template<> inline int precision<int>() { return 0; } +inline int real(const int& x) { return x; } +inline int imag(const int& x) { EIGEN_UNUSED(x); return 0; } +inline int conj(const int& x) { return x; } +inline int abs(const int& x) { return std::abs(x); } +inline int abs2(const int& x) { return x*x; } +inline int sqrt(const int& x) +{ + EIGEN_UNUSED(x); + // Taking the square root of integers is not allowed + // (the square root does not always exist within the integers). + // Please cast to a floating-point type. + assert(false); +} +template<> inline int random() +{ + // "rand() % n" is bad, they say, because the low-order bits are not random enough. + // However here, 21 is odd, so random() % 21 uses the high-order bits + // as well, so there's no problem. + return (std::rand() % 21) - 10; +} +inline bool isMuchSmallerThan(const int& a, const int& b, const int& prec = precision<int>()) +{ + EIGEN_UNUSED(b); + EIGEN_UNUSED(prec); + return a == 0; +} +inline bool isApprox(const int& a, const int& b, const int& prec = precision<int>()) +{ + EIGEN_UNUSED(prec); + return a == b; +} +inline bool isApproxOrLessThan(const int& a, const int& b, const int& prec = precision<int>()) +{ + EIGEN_UNUSED(prec); + return a <= b; +} + +template<> inline float precision<float>() { return 1e-5f; } +inline float real(const float& x) { return x; } +inline float imag(const float& x) { EIGEN_UNUSED(x); return 0.f; } +inline float conj(const float& x) { return x; } +inline float abs(const float& x) { return std::abs(x); } +inline float abs2(const float& x) { return x*x; } +inline float sqrt(const float& x) { return std::sqrt(x); } +template<> inline float random() +{ + return std::rand() / (RAND_MAX/20.0f) - 10.0f; +} +inline bool isMuchSmallerThan(const float& a, const float& b, const float& prec = precision<float>()) +{ + return std::abs(a) <= std::abs(b) * prec; +} +inline bool isApprox(const float& a, const float& b, const float& prec = precision<float>()) +{ + return std::abs(a - b) <= std::min(std::abs(a), std::abs(b)) * prec; +} +inline bool isApproxOrLessThan(const float& a, const float& b, const float& prec = precision<float>()) +{ + return a <= b || isApprox(a, b, prec); +} + +template<> inline double precision<double>() { return 1e-11; } +inline double real(const double& x) { return x; } +inline double imag(const double& x) { EIGEN_UNUSED(x); return 0.; } +inline double conj(const double& x) { return x; } +inline double abs(const double& x) { return std::abs(x); } +inline double abs2(const double& x) { return x*x; } +inline double sqrt(const double& x) { return std::sqrt(x); } +template<> inline double random() +{ + return std::rand() / (RAND_MAX/20.0) - 10.0; +} +inline bool isMuchSmallerThan(const double& a, const double& b, const double& prec = precision<double>()) +{ + return std::abs(a) <= std::abs(b) * prec; +} +inline bool isApprox(const double& a, const double& b, const double& prec = precision<double>()) +{ + return std::abs(a - b) <= std::min(std::abs(a), std::abs(b)) * prec; +} +inline bool isApproxOrLessThan(const double& a, const double& b, const double& prec = precision<double>()) +{ + return a <= b || isApprox(a, b, prec); +} + +template<> inline float precision<std::complex<float> >() { return precision<float>(); } +inline float real(const std::complex<float>& x) { return std::real(x); } +inline float imag(const std::complex<float>& x) { return std::imag(x); } +inline std::complex<float> conj(const std::complex<float>& x) { return std::conj(x); } +inline float abs(const std::complex<float>& x) { return std::abs(x); } +inline float abs2(const std::complex<float>& x) { return std::norm(x); } +inline std::complex<float> sqrt(const std::complex<float>& x) +{ + EIGEN_UNUSED(x); + // Taking the square roots of complex numbers is not allowed, + // as this is ambiguous (there are two square roots). + // What were you trying to do? + assert(false); +} +template<> inline std::complex<float> random() +{ + return std::complex<float>(random<float>(), random<float>()); +} +inline bool isMuchSmallerThan(const std::complex<float>& a, const std::complex<float>& b, const float& prec = precision<float>()) +{ + return abs2(a) <= abs2(b) * prec * prec; +} +inline bool isApprox(const std::complex<float>& a, const std::complex<float>& b, const float& prec = precision<float>()) +{ + return isApprox(std::real(a), std::real(b), prec) + && isApprox(std::imag(a), std::imag(b), prec); +} +// isApproxOrLessThan wouldn't make sense for complex numbers + +template<> inline double precision<std::complex<double> >() { return precision<double>(); } +inline double real(const std::complex<double>& x) { return std::real(x); } +inline double imag(const std::complex<double>& x) { return std::imag(x); } +inline std::complex<double> conj(const std::complex<double>& x) { return std::conj(x); } +inline double abs(const std::complex<double>& x) { return std::abs(x); } +inline double abs2(const std::complex<double>& x) { return std::norm(x); } +template<> inline std::complex<double> random() +{ + return std::complex<double>(random<double>(), random<double>()); +} +inline bool isMuchSmallerThan(const std::complex<double>& a, const std::complex<double>& b, const double& prec = precision<double>()) +{ + return abs2(a) <= abs2(b) * prec * prec; +} +inline bool isApprox(const std::complex<double>& a, const std::complex<double>& b, const double& prec = precision<double>()) +{ + return isApprox(std::real(a), std::real(b), prec) + && isApprox(std::imag(a), std::imag(b), prec); +} +// isApproxOrLessThan wouldn't make sense for complex numbers + +#endif // EIGEN_NUMTRAITS_H diff --git a/src/Core/Random.h b/src/Core/Random.h index ba61306ab..a7e1d0f80 100644 --- a/src/Core/Random.h +++ b/src/Core/Random.h @@ -53,7 +53,7 @@ template<typename MatrixType> class Random { EIGEN_UNUSED(row); EIGEN_UNUSED(col); - return NumTraits<Scalar>::random(); + return random<Scalar>(); } protected: diff --git a/test/adjoint.cpp b/test/adjoint.cpp index 046bbbc5b..2d6aac535 100644 --- a/test/adjoint.cpp +++ b/test/adjoint.cpp @@ -25,6 +25,8 @@ #include "main.h" +namespace Eigen { + template<typename MatrixType> void adjoint(const MatrixType& m) { /* this test covers the following files: @@ -49,8 +51,8 @@ template<typename MatrixType> void adjoint(const MatrixType& m) v3 = VectorType::random(rows), vzero = VectorType::zero(rows); - Scalar s1 = NumTraits<Scalar>::random(), - s2 = NumTraits<Scalar>::random(); + Scalar s1 = random<Scalar>(), + s2 = random<Scalar>(); // check involutivity of adjoint, transpose, conjugate QVERIFY(m1.transpose().transpose().isApprox(m1)); @@ -67,28 +69,32 @@ template<typename MatrixType> void adjoint(const MatrixType& m) QVERIFY((m1.adjoint() * m2).adjoint().isApprox(m2.adjoint() * m1)); QVERIFY((m1.transpose() * m2).conjugate().isApprox(m1.adjoint() * m2.conjugate())); QVERIFY((s1 * m1).transpose().isApprox(s1 * m1.transpose())); - QVERIFY((s1 * m1).conjugate().isApprox(NumTraits<Scalar>::conj(s1) * m1.conjugate())); - QVERIFY((s1 * m1).adjoint().isApprox(NumTraits<Scalar>::conj(s1) * m1.adjoint())); + QVERIFY((s1 * m1).conjugate().isApprox(conj(s1) * m1.conjugate())); + QVERIFY((s1 * m1).adjoint().isApprox(conj(s1) * m1.adjoint())); // check basic properties of dot, norm, norm2 typedef typename NumTraits<Scalar>::Real RealScalar; - QVERIFY(NumTraits<Scalar>::isApprox((s1 * v1 + s2 * v2).dot(v3), s1 * v1.dot(v3) + s2 * v2.dot(v3))); - QVERIFY(NumTraits<Scalar>::isApprox(v3.dot(s1 * v1 + s2 * v2), NumTraits<Scalar>::conj(s1) * v3.dot(v1) + NumTraits<Scalar>::conj(s2) * v3.dot(v2))); - QVERIFY(NumTraits<Scalar>::isApprox(NumTraits<Scalar>::conj(v1.dot(v2)), v2.dot(v1))); - QVERIFY(NumTraits<RealScalar>::isApprox(abs(v1.dot(v1)), v1.norm2())); - if(NumTraits<Scalar>::HasFloatingPoint) QVERIFY(NumTraits<RealScalar>::isApprox(v1.norm2(), v1.norm() * v1.norm())); - QVERIFY(NumTraits<RealScalar>::isMuchSmallerThan(abs(vzero.dot(v1)), 1)); - QVERIFY(NumTraits<RealScalar>::isMuchSmallerThan(vzero.norm(), 1)); + QVERIFY(isApprox((s1 * v1 + s2 * v2).dot(v3), s1 * v1.dot(v3) + s2 * v2.dot(v3))); + QVERIFY(isApprox(v3.dot(s1 * v1 + s2 * v2), conj(s1) * v3.dot(v1) + conj(s2) * v3.dot(v2))); + QVERIFY(isApprox(conj(v1.dot(v2)), v2.dot(v1))); + QVERIFY(isApprox(abs(v1.dot(v1)), v1.norm2())); + if(NumTraits<Scalar>::HasFloatingPoint) + QVERIFY(isApprox(v1.norm2(), v1.norm() * v1.norm())); + QVERIFY(isMuchSmallerThan(abs(vzero.dot(v1)), static_cast<RealScalar>(1))); + if(NumTraits<Scalar>::HasFloatingPoint) + QVERIFY(isMuchSmallerThan(vzero.norm(), static_cast<RealScalar>(1))); // check compatibility of dot and adjoint - QVERIFY(NumTraits<Scalar>::isApprox(v1.dot(square * v2), (square.adjoint() * v1).dot(v2))); + QVERIFY(isApprox(v1.dot(square * v2), (square.adjoint() * v1).dot(v2))); } void EigenTest::testAdjoint() { adjoint(Matrix<float, 1, 1>()); - adjoint(Matrix<complex<double>, 4, 4>()); + adjoint(Matrix4cd()); adjoint(MatrixXcf(3, 3)); adjoint(MatrixXi(8, 12)); adjoint(MatrixXd(20, 20)); } + +} // namespace Eigen diff --git a/test/basicstuff.cpp b/test/basicstuff.cpp index d929b0883..babddb670 100644 --- a/test/basicstuff.cpp +++ b/test/basicstuff.cpp @@ -25,6 +25,8 @@ #include "main.h" +namespace Eigen { + template<typename MatrixType> void basicStuff(const MatrixType& m) { /* this test covers the following files: @@ -56,14 +58,15 @@ template<typename MatrixType> void basicStuff(const MatrixType& m) v2 = VectorType::random(rows), vzero = VectorType::zero(rows); - Scalar s1 = NumTraits<Scalar>::random(), - s2 = NumTraits<Scalar>::random(); + Scalar s1 = random<Scalar>(), + s2 = random<Scalar>(); // test Fuzzy.h and Zero.h. QVERIFY(v1.isApprox(v1)); QVERIFY(!v1.isApprox(2*v1)); QVERIFY(vzero.isMuchSmallerThan(v1)); - QVERIFY(vzero.isMuchSmallerThan(v1.norm())); + if(NumTraits<Scalar>::HasFloatingPoint) + QVERIFY(vzero.isMuchSmallerThan(v1.norm())); QVERIFY(!v1.isMuchSmallerThan(v1)); QVERIFY(vzero.isApprox(v1-v1)); QVERIFY(m1.isApprox(m1)); @@ -134,8 +137,10 @@ template<typename MatrixType> void basicStuff(const MatrixType& m) void EigenTest::testBasicStuff() { basicStuff(Matrix<float, 1, 1>()); - basicStuff(Matrix<complex<double>, 4, 4>()); + basicStuff(Matrix4cd()); basicStuff(MatrixXcf(3, 3)); basicStuff(MatrixXi(8, 12)); basicStuff(MatrixXd(20, 20)); } + +} // namespace Eigen diff --git a/test/main.cpp b/test/main.cpp index f5c758961..e1272087a 100644 --- a/test/main.cpp +++ b/test/main.cpp @@ -25,13 +25,14 @@ #include "main.h" -EigenTest::EigenTest() +Eigen::EigenTest::EigenTest() { - unsigned int t = (unsigned int) time( NULL ); + unsigned int t = (unsigned int) time(NULL); qDebug() << "Initializing random number generator with seed" << t; srand(t); } -QTEST_APPLESS_MAIN( EigenTest ) +QTEST_APPLESS_MAIN(Eigen::EigenTest) + #include "main.moc" diff --git a/test/main.h b/test/main.h index c752b6ec6..3d90a2f48 100644 --- a/test/main.h +++ b/test/main.h @@ -29,12 +29,10 @@ #include <QtTest/QtTest> #include "../src/Core.h" -using namespace Eigen; - #include <cstdlib> #include <ctime> -using namespace std; +namespace Eigen { class EigenTest : public QObject { @@ -48,4 +46,6 @@ class EigenTest : public QObject void testAdjoint(); }; +} // end namespace Eigen + #endif // EIGEN_TEST_MAIN_H |