From f6003f08737eee960a70541a750e1675a470cdcf Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Thu, 14 Apr 2016 09:47:26 -0700 Subject: Made the test msvc friendly --- unsupported/test/cxx11_float16.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/unsupported/test/cxx11_float16.cpp b/unsupported/test/cxx11_float16.cpp index 2dc0872d8..db5f2130b 100644 --- a/unsupported/test/cxx11_float16.cpp +++ b/unsupported/test/cxx11_float16.cpp @@ -122,6 +122,8 @@ void test_comparison() VERIFY(half(1.0f) != half(2.0f)); // Comparisons with NaNs and infinities. +#if !EIGEN_COMP_MSVC + // Visual Studio errors out on divisions by 0 VERIFY(!(half(0.0 / 0.0) == half(0.0 / 0.0))); VERIFY(half(0.0 / 0.0) != half(0.0 / 0.0)); @@ -132,6 +134,7 @@ void test_comparison() VERIFY(half(1.0) < half(1.0 / 0.0)); VERIFY(half(1.0) > half(-1.0 / 0.0)); +#endif } void test_functions() -- cgit v1.2.3 From c7167fee0eaf188fbe7ce7b877971b928f798a7a Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Thu, 14 Apr 2016 10:08:33 -0700 Subject: Added support for fp16 to the sigmoid function --- unsupported/Eigen/CXX11/src/Tensor/TensorFunctors.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/unsupported/Eigen/CXX11/src/Tensor/TensorFunctors.h b/unsupported/Eigen/CXX11/src/Tensor/TensorFunctors.h index b7c13f67f..4ef48c64c 100644 --- a/unsupported/Eigen/CXX11/src/Tensor/TensorFunctors.h +++ b/unsupported/Eigen/CXX11/src/Tensor/TensorFunctors.h @@ -64,7 +64,7 @@ struct scalar_sigmoid_op { EIGEN_EMPTY_STRUCT_CTOR(scalar_sigmoid_op) EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T operator()(const T& x) const { const T one = T(1); - return one / (one + std::exp(-x)); + return one / (one + numext::exp(-x)); } template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE @@ -799,7 +799,7 @@ class GaussianGenerator { T offset = coordinates[i] - m_means[i]; tmp += offset * offset / m_two_sigmas[i]; } - return std::exp(-tmp); + return numext::exp(-tmp); } private: -- cgit v1.2.3 From 5c13765ee333bb78b5c7baeb515eed97c59b6c1d Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Thu, 14 Apr 2016 10:24:52 -0700 Subject: Added ability to printf fp16 --- Eigen/src/Core/arch/CUDA/Half.h | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/Eigen/src/Core/arch/CUDA/Half.h b/Eigen/src/Core/arch/CUDA/Half.h index 482b654ac..8249ce2eb 100644 --- a/Eigen/src/Core/arch/CUDA/Half.h +++ b/Eigen/src/Core/arch/CUDA/Half.h @@ -510,6 +510,11 @@ static EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC int (isfinite)(const Eigen::half& a namespace std { +ostream& operator << (ostream& os, const Eigen::half& v) { + os << static_cast(v); + return os; +} + #if __cplusplus > 199711L template <> struct hash { -- cgit v1.2.3 From 7b3d7acebeadb443d8e3ac9756359d507324cc82 Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Thu, 14 Apr 2016 10:25:50 -0700 Subject: Added support for fp16 to test_isApprox, test_isMuchSmallerThan, and test_isApproxOrLessThan --- test/main.h | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/test/main.h b/test/main.h index bba5e7570..dbb496b89 100644 --- a/test/main.h +++ b/test/main.h @@ -316,9 +316,9 @@ inline bool test_isMuchSmallerThan(const float& a, const float& b) { return internal::isMuchSmallerThan(a, b, test_precision()); } inline bool test_isApproxOrLessThan(const float& a, const float& b) { return internal::isApproxOrLessThan(a, b, test_precision()); } + inline bool test_isApprox(const double& a, const double& b) { return internal::isApprox(a, b, test_precision()); } - inline bool test_isMuchSmallerThan(const double& a, const double& b) { return internal::isMuchSmallerThan(a, b, test_precision()); } inline bool test_isApproxOrLessThan(const double& a, const double& b) @@ -359,6 +359,12 @@ inline bool test_isApproxOrLessThan(const long double& a, const long double& b) { return internal::isApproxOrLessThan(a, b, test_precision()); } #endif // EIGEN_TEST_NO_LONGDOUBLE +inline bool test_isApprox(const half& a, const half& b) +{ return internal::isApprox(a, b, test_precision()); } +inline bool test_isMuchSmallerThan(const half& a, const half& b) +{ return internal::isMuchSmallerThan(a, b, test_precision()); } +inline bool test_isApproxOrLessThan(const half& a, const half& b) +{ return internal::isApproxOrLessThan(a, b, test_precision()); } // test_relative_error returns the relative difference between a and b as a real scalar as used in isApprox. template @@ -426,9 +432,7 @@ template typename NumTraits::Real test_relative_error(const T1 &a, const T2 &b, typename internal::enable_if::Real>::value, T1>::type* = 0) { typedef typename NumTraits::Real RealScalar; - using std::min; - using std::sqrt; - return sqrt(RealScalar(numext::abs2(a-b))/RealScalar((min)(numext::abs2(a),numext::abs2(b)))); + return numext::sqrt(RealScalar(numext::abs2(a-b))/RealScalar((numext::mini)(numext::abs2(a),numext::abs2(b)))); } template -- cgit v1.2.3 From 72510c80e1a7406af915f0851e6bfbe605d3f436 Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Thu, 14 Apr 2016 10:27:24 -0700 Subject: Added basic test for trigonometric functions on fp16 --- unsupported/test/cxx11_float16.cpp | 26 ++++++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/unsupported/test/cxx11_float16.cpp b/unsupported/test/cxx11_float16.cpp index db5f2130b..2d9376d29 100644 --- a/unsupported/test/cxx11_float16.cpp +++ b/unsupported/test/cxx11_float16.cpp @@ -137,7 +137,7 @@ void test_comparison() #endif } -void test_functions() +void test_basic_functions() { VERIFY_IS_EQUAL(float(numext::abs(half(3.5f))), 3.5f); VERIFY_IS_EQUAL(float(numext::abs(half(-3.5f))), 3.5f); @@ -149,10 +149,32 @@ void test_functions() VERIFY_IS_APPROX(float(numext::log(half(10.0f))), 2.30273f); } +void test_trigonometric_functions() +{ + VERIFY_IS_APPROX(numext::cos(half(0.0f)), half(cosf(0.0f))); + VERIFY_IS_APPROX(numext::cos(half(EIGEN_PI)), half(cosf(EIGEN_PI))); + //VERIFY_IS_APPROX(numext::cos(half(EIGEN_PI/2)), half(cosf(EIGEN_PI/2))); + //VERIFY_IS_APPROX(numext::cos(half(3*EIGEN_PI/2)), half(cosf(3*EIGEN_PI/2))); + VERIFY_IS_APPROX(numext::cos(half(3.5f)), half(cosf(3.5f))); + + VERIFY_IS_APPROX(numext::sin(half(0.0f)), half(sinf(0.0f))); + // VERIFY_IS_APPROX(numext::sin(half(EIGEN_PI)), half(sinf(EIGEN_PI))); + VERIFY_IS_APPROX(numext::sin(half(EIGEN_PI/2)), half(sinf(EIGEN_PI/2))); + VERIFY_IS_APPROX(numext::sin(half(3*EIGEN_PI/2)), half(sinf(3*EIGEN_PI/2))); + VERIFY_IS_APPROX(numext::sin(half(3.5f)), half(sinf(3.5f))); + + VERIFY_IS_APPROX(numext::tan(half(0.0f)), half(tanf(0.0f))); + // VERIFY_IS_APPROX(numext::tan(half(EIGEN_PI)), half(tanf(EIGEN_PI))); + // VERIFY_IS_APPROX(numext::tan(half(EIGEN_PI/2)), half(tanf(EIGEN_PI/2))); + //VERIFY_IS_APPROX(numext::tan(half(3*EIGEN_PI/2)), half(tanf(3*EIGEN_PI/2))); + VERIFY_IS_APPROX(numext::tan(half(3.5f)), half(tanf(3.5f))); +} + void test_cxx11_float16() { CALL_SUBTEST(test_conversion()); CALL_SUBTEST(test_arithmetic()); CALL_SUBTEST(test_comparison()); - CALL_SUBTEST(test_functions()); + CALL_SUBTEST(test_basic_functions()); + CALL_SUBTEST(test_trigonometric_functions()); } -- cgit v1.2.3 From 6f23e945f6fd75d8d7b48d83f01976f91da0c24f Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Thu, 14 Apr 2016 10:32:52 -0700 Subject: Added simple test for numext::sqrt and numext::pow on fp16 --- unsupported/test/cxx11_float16.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/unsupported/test/cxx11_float16.cpp b/unsupported/test/cxx11_float16.cpp index 2d9376d29..b437868d7 100644 --- a/unsupported/test/cxx11_float16.cpp +++ b/unsupported/test/cxx11_float16.cpp @@ -142,6 +142,12 @@ void test_basic_functions() VERIFY_IS_EQUAL(float(numext::abs(half(3.5f))), 3.5f); VERIFY_IS_EQUAL(float(numext::abs(half(-3.5f))), 3.5f); + VERIFY_IS_APPROX(float(numext::sqrt(half(0.0f))), 0.0f); + VERIFY_IS_APPROX(float(numext::sqrt(half(4.0f))), 2.0f); + + VERIFY_IS_APPROX(float(numext::pow(half(0.0f), half(1.0f))), 0.0f); + VERIFY_IS_APPROX(float(numext::pow(half(2.0f), half(2.0f))), 4.0f); + VERIFY_IS_EQUAL(float(numext::exp(half(0.0f))), 1.0f); VERIFY_IS_APPROX(float(numext::exp(half(EIGEN_PI))), float(20.0 + EIGEN_PI)); -- cgit v1.2.3 From 2b6e3de02f6f141c6bab523c54cda432d796eec7 Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Thu, 14 Apr 2016 11:39:18 -0700 Subject: Added tests to validate flooring and ceiling of fp16 --- unsupported/test/cxx11_float16.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/unsupported/test/cxx11_float16.cpp b/unsupported/test/cxx11_float16.cpp index b437868d7..273dcbc11 100644 --- a/unsupported/test/cxx11_float16.cpp +++ b/unsupported/test/cxx11_float16.cpp @@ -142,6 +142,12 @@ void test_basic_functions() VERIFY_IS_EQUAL(float(numext::abs(half(3.5f))), 3.5f); VERIFY_IS_EQUAL(float(numext::abs(half(-3.5f))), 3.5f); + VERIFY_IS_EQUAL(float(numext::floor(half(3.5f))), 3.0f); + VERIFY_IS_EQUAL(float(numext::floor(half(-3.5f))), -4.0f); + + VERIFY_IS_EQUAL(float(numext::ceil(half(3.5f))), 4.0f); + VERIFY_IS_EQUAL(float(numext::ceil(half(-3.5f))), -3.0f); + VERIFY_IS_APPROX(float(numext::sqrt(half(0.0f))), 0.0f); VERIFY_IS_APPROX(float(numext::sqrt(half(4.0f))), 2.0f); -- cgit v1.2.3 From 5912ad877c6fe0072c56e8d2f70b315a1f4da6ce Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Thu, 14 Apr 2016 11:40:14 -0700 Subject: Silenced a compilation warning --- unsupported/Eigen/CXX11/src/Tensor/TensorUInt128.h | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/unsupported/Eigen/CXX11/src/Tensor/TensorUInt128.h b/unsupported/Eigen/CXX11/src/Tensor/TensorUInt128.h index 3e56589c3..5950f38e2 100644 --- a/unsupported/Eigen/CXX11/src/Tensor/TensorUInt128.h +++ b/unsupported/Eigen/CXX11/src/Tensor/TensorUInt128.h @@ -53,9 +53,7 @@ struct TensorUInt128 template EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE explicit TensorUInt128(const T& x) : high(0), low(x) { - typedef typename conditional::type UnsignedT; - typedef typename conditional::type UnsignedLow; - eigen_assert(static_cast(x) <= static_cast(NumTraits::highest())); + eigen_assert((static_cast::type>(x) <= static_cast::type>(NumTraits::highest()))); eigen_assert(x >= 0); } -- cgit v1.2.3 From 5379d2b5944f2c26ff0ddce65fc6f99f5182f7b7 Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Thu, 14 Apr 2016 11:40:48 -0700 Subject: Inline the << operator on half floats --- Eigen/src/Core/arch/CUDA/Half.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Eigen/src/Core/arch/CUDA/Half.h b/Eigen/src/Core/arch/CUDA/Half.h index 8249ce2eb..bdf97dcd6 100644 --- a/Eigen/src/Core/arch/CUDA/Half.h +++ b/Eigen/src/Core/arch/CUDA/Half.h @@ -510,7 +510,7 @@ static EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC int (isfinite)(const Eigen::half& a namespace std { -ostream& operator << (ostream& os, const Eigen::half& v) { +EIGEN_STRONG_INLINE ostream& operator << (ostream& os, const Eigen::half& v) { os << static_cast(v); return os; } -- cgit v1.2.3 From 7718749fee835095f0671fa6ce5d257609f8e56b Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Thu, 14 Apr 2016 11:51:54 -0700 Subject: Force the inlining of the << operator on half floats --- Eigen/src/Core/arch/CUDA/Half.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Eigen/src/Core/arch/CUDA/Half.h b/Eigen/src/Core/arch/CUDA/Half.h index bdf97dcd6..9ecc4fd88 100644 --- a/Eigen/src/Core/arch/CUDA/Half.h +++ b/Eigen/src/Core/arch/CUDA/Half.h @@ -510,7 +510,7 @@ static EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC int (isfinite)(const Eigen::half& a namespace std { -EIGEN_STRONG_INLINE ostream& operator << (ostream& os, const Eigen::half& v) { +EIGEN_ALWAYS_INLINE ostream& operator << (ostream& os, const Eigen::half& v) { os << static_cast(v); return os; } -- cgit v1.2.3 From 20f387fafa5dbab90c240612e33e5c13d215ac5f Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Thu, 14 Apr 2016 22:46:55 +0200 Subject: Improve numerical robustness of JacoviSVD: - avoid noise amplification in complex to real conversion - compare off-diagonal entries to the current biggest diagonal entry: no need to bother about a 2x2 block containing ridiculously small entries compared to the rest of the matrix. --- Eigen/src/SVD/JacobiSVD.h | 42 ++++++++++++++++++++++++------------------ 1 file changed, 24 insertions(+), 18 deletions(-) diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index f08776bc6..1940c8294 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -350,7 +350,8 @@ template struct svd_precondition_2x2_block_to_be_real { typedef JacobiSVD SVD; - static bool run(typename SVD::WorkMatrixType&, SVD&, Index, Index, const typename MatrixType::RealScalar&) { return true; } + typedef typename MatrixType::RealScalar RealScalar; + static bool run(typename SVD::WorkMatrixType&, SVD&, Index, Index, RealScalar&) { return true; } }; template @@ -359,25 +360,30 @@ struct svd_precondition_2x2_block_to_be_real typedef JacobiSVD SVD; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; - static bool run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q, const typename MatrixType::RealScalar& precision) + static bool run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q, RealScalar& maxDiagEntry) { using std::sqrt; + using std::abs; Scalar z; JacobiRotation rot; RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p,p)) + numext::abs2(work_matrix.coeff(q,p))); - + + const RealScalar considerAsZero = (std::numeric_limits::min)(); + const RealScalar precision = NumTraits::epsilon(); + if(n==0) { // make sure first column is zero work_matrix.coeffRef(p,p) = work_matrix.coeffRef(q,p) = Scalar(0); - if(work_matrix.coeff(p,q)!=Scalar(0)) + + if(abs(numext::imag(work_matrix.coeff(p,q)))>considerAsZero) { // work_matrix.coeff(p,q) can be zero if work_matrix.coeff(q,p) is not zero but small enough to underflow when computing n z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); work_matrix.row(p) *= z; if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z); } - if(work_matrix.coeff(q,q)!=Scalar(0)) + if(abs(numext::imag(work_matrix.coeff(q,q)))>considerAsZero) { z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); work_matrix.row(q) *= z; @@ -391,13 +397,13 @@ struct svd_precondition_2x2_block_to_be_real rot.s() = work_matrix.coeff(q,p) / n; work_matrix.applyOnTheLeft(p,q,rot); if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint()); - if(work_matrix.coeff(p,q) != Scalar(0)) + if(abs(numext::imag(work_matrix.coeff(p,q)))>considerAsZero) { z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); work_matrix.col(q) *= z; if(svd.computeV()) svd.m_matrixV.col(q) *= z; } - if(work_matrix.coeff(q,q) != Scalar(0)) + if(abs(numext::imag(work_matrix.coeff(q,q)))>considerAsZero) { z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); work_matrix.row(q) *= z; @@ -405,11 +411,11 @@ struct svd_precondition_2x2_block_to_be_real } } - const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits::denorm_min(); - RealScalar threshold = numext::maxi(considerAsZero, - precision * numext::maxi(abs(work_matrix.coeff(p,p)), abs(work_matrix.coeff(q,q)))); - // return true if we still have some work to do - return numext::abs(work_matrix(p,q)) > threshold || numext::abs(work_matrix(q,p)) > threshold; + // update largest diagonal entry + maxDiagEntry = numext::maxi(maxDiagEntry,numext::maxi(abs(work_matrix.coeff(p,p)), abs(work_matrix.coeff(q,q)))); + // and check whether the 2x2 block is already diagonal + RealScalar threshold = numext::maxi(considerAsZero, precision * maxDiagEntry); + return abs(work_matrix.coeff(p,q))>threshold || abs(work_matrix.coeff(q,p)) > threshold; } }; @@ -426,7 +432,6 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, JacobiRotation rot1; RealScalar t = m.coeff(0,0) + m.coeff(1,1); RealScalar d = m.coeff(1,0) - m.coeff(0,1); - if(d == RealScalar(0)) { rot1.s() = RealScalar(0); @@ -719,6 +724,7 @@ JacobiSVD::compute(const MatrixType& matrix, unsig } /*** step 2. The main Jacobi SVD iteration. ***/ + RealScalar maxDiagEntry = m_workMatrix.cwiseAbs().diagonal().maxCoeff(); bool finished = false; while(!finished) @@ -734,16 +740,13 @@ JacobiSVD::compute(const MatrixType& matrix, unsig // if this 2x2 sub-matrix is not diagonal already... // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't // keep us iterating forever. Similarly, small denormal numbers are considered zero. - RealScalar threshold = numext::maxi(considerAsZero, - precision * numext::maxi(abs(m_workMatrix.coeff(p,p)), - abs(m_workMatrix.coeff(q,q)))); - // We compare both values to threshold instead of calling max to be robust to NaN (See bug 791) + RealScalar threshold = numext::maxi(considerAsZero, precision * maxDiagEntry); if(abs(m_workMatrix.coeff(p,q))>threshold || abs(m_workMatrix.coeff(q,p)) > threshold) { finished = false; // perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal // the complex to real operation returns true is the updated 2x2 block is not already diagonal - if(internal::svd_precondition_2x2_block_to_be_real::run(m_workMatrix, *this, p, q, precision)) + if(internal::svd_precondition_2x2_block_to_be_real::run(m_workMatrix, *this, p, q, maxDiagEntry)) { JacobiRotation j_left, j_right; internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right); @@ -754,6 +757,9 @@ JacobiSVD::compute(const MatrixType& matrix, unsig m_workMatrix.applyOnTheRight(p,q,j_right); if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right); + + // keep track of the largest diagonal coefficient + maxDiagEntry = numext::maxi(maxDiagEntry,numext::maxi(abs(m_workMatrix.coeff(p,p)), abs(m_workMatrix.coeff(q,q)))); } } } -- cgit v1.2.3 From 68897c52f3c8cd37824f87cd9582cf98a5c9eb32 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Thu, 14 Apr 2016 22:47:30 +0200 Subject: Add extreme values to the imaginary part for SVD unit tests. --- test/svd_fill.h | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/test/svd_fill.h b/test/svd_fill.h index 7e44b3d05..1bbe645ee 100644 --- a/test/svd_fill.h +++ b/test/svd_fill.h @@ -80,6 +80,8 @@ void svd_fill_random(MatrixType &m, int Option = 0) Index i = internal::random(0,m.rows()-1); Index j = internal::random(0,m.cols()-1); m(j,i) = m(i,j) = samples(internal::random(0,samples.size()-1)); + if(NumTraits::IsComplex) + *(&numext::real_ref(m(j,i))+1) = *(&numext::real_ref(m(i,j))+1) = samples.real()(internal::random(0,samples.size()-1)); } } } @@ -91,8 +93,14 @@ void svd_fill_random(MatrixType &m, int Option = 0) if(!(dup && unit_uv)) { Index n = internal::random(0,m.size()-1); - for(Index i=0; i(0,m.rows()-1), internal::random(0,m.cols()-1)) = samples(internal::random(0,samples.size()-1)); + for(Index k=0; k(0,m.rows()-1); + Index j = internal::random(0,m.cols()-1); + m(i,j) = samples(internal::random(0,samples.size()-1)); + if(NumTraits::IsComplex) + *(&numext::real_ref(m(i,j))+1) = samples.real()(internal::random(0,samples.size()-1)); + } } } } -- cgit v1.2.3