From 7a9519a9be524a879b44a1912ae694cfe4ee4ef7 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 14 Jul 2009 23:00:53 +0200 Subject: [PATCH 1/2] fix typo in blue norm --- Eigen/src/Core/Dot.h | 2 +- test/adjoint.cpp | 32 +++++++++----------------------- 2 files changed, 10 insertions(+), 24 deletions(-) diff --git a/Eigen/src/Core/Dot.h b/Eigen/src/Core/Dot.h index 4f185ea5b..c5f2e8505 100644 --- a/Eigen/src/Core/Dot.h +++ b/Eigen/src/Core/Dot.h @@ -394,7 +394,7 @@ MatrixBase::blueNorm() const { ax = ei_abs(coeff(j)); if(ax > b2) abig += ei_abs2(ax*s2m); - else if(ax < b2) asml += ei_abs2(ax*s1m); + else if(ax < b1) asml += ei_abs2(ax*s1m); else amed += ei_abs2(ax); } if(abig > Scalar(0)) diff --git a/test/adjoint.cpp b/test/adjoint.cpp index 14ee44a0f..403f16a45 100644 --- a/test/adjoint.cpp +++ b/test/adjoint.cpp @@ -114,29 +114,15 @@ template void adjoint(const MatrixType& m) void test_adjoint() { -// for(int i = 0; i < g_repeat; i++) { -// CALL_SUBTEST( adjoint(Matrix()) ); -// CALL_SUBTEST( adjoint(Matrix3d()) ); -// CALL_SUBTEST( adjoint(Matrix4f()) ); -// CALL_SUBTEST( adjoint(MatrixXcf(4, 4)) ); -// CALL_SUBTEST( adjoint(MatrixXi(8, 12)) ); -// CALL_SUBTEST( adjoint(MatrixXf(21, 21)) ); -// } - // test a large matrix only once -// CALL_SUBTEST( adjoint(Matrix()) ); - for(int i = 0; i < g_repeat; i++) - { - std::cerr.precision(20); - int s = 1000000; - double y = 1.131242353467546478463457843445677435233e23 * ei_abs(ei_random()); - VectorXf v = VectorXf::Ones(s) * y; -// Vector4f x(v.segment(0,s/4).blueNorm(), v.segment(s/4+1,s/4).blueNorm(), -// v.segment((s/2)+1,s/4).blueNorm(), v.segment(3*s/4+1,s - 3*s/4-1).blueNorm()); -// std::cerr << v.norm() << " == " << v.stableNorm() << " == " << v.blueNorm() << " == " << x.norm() << "\n"; - std::cerr << v.norm() << "\n" << v.stableNorm() << "\n" << v.blueNorm() << "\n" << ei_sqrt(double(s)) * y << "\n\n\n"; - -// VectorXd d = VectorXd::Ones(s) * y;//v.cast(); -// std::cerr << d.norm() << "\n" << d.stableNorm() << "\n" << d.blueNorm() << "\n" << ei_sqrt(double(s)) * y << "\n\n\n"; + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST( adjoint(Matrix()) ); + CALL_SUBTEST( adjoint(Matrix3d()) ); + CALL_SUBTEST( adjoint(Matrix4f()) ); + CALL_SUBTEST( adjoint(MatrixXcf(4, 4)) ); + CALL_SUBTEST( adjoint(MatrixXi(8, 12)) ); + CALL_SUBTEST( adjoint(MatrixXf(21, 21)) ); } + // test a large matrix only once + CALL_SUBTEST( adjoint(Matrix()) ); } From 1578421ed14c23fa5c7ab3c818a069f1c1cefb8a Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 15 Jul 2009 14:20:45 +0200 Subject: [PATCH 2/2] fix issue #25 : the problem was that we assumed Dynamic was a multiple of a packet size (also disable the test of blueNorm) --- Eigen/src/Core/Product.h | 4 ++-- Eigen/src/Core/util/XprHelper.h | 6 ++++-- Eigen/src/QR/Tridiagonalization.h | 1 + test/adjoint.cpp | 5 +++-- 4 files changed, 10 insertions(+), 6 deletions(-) diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index c2317e425..05a8221ee 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h @@ -138,10 +138,10 @@ struct ei_traits > RhsRowMajor = RhsFlags & RowMajorBit, CanVectorizeRhs = RhsRowMajor && (RhsFlags & PacketAccessBit) - && (ColsAtCompileTime % ei_packet_traits::size == 0), + && (ColsAtCompileTime == Dynamic || (ColsAtCompileTime % ei_packet_traits::size) == 0), CanVectorizeLhs = (!LhsRowMajor) && (LhsFlags & PacketAccessBit) - && (RowsAtCompileTime % ei_packet_traits::size == 0), + && (RowsAtCompileTime == Dynamic || (RowsAtCompileTime % ei_packet_traits::size) == 0), EvalToRowMajor = RhsRowMajor && (ProductMode==(int)CacheFriendlyProduct ? LhsRowMajor : (!CanVectorizeLhs)), diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h index d8bbe11d5..4ebf33e8a 100644 --- a/Eigen/src/Core/util/XprHelper.h +++ b/Eigen/src/Core/util/XprHelper.h @@ -86,9 +86,11 @@ class ei_compute_matrix_flags { enum { row_major_bit = Options&RowMajor ? RowMajorBit : 0, - inner_max_size = row_major_bit ? MaxCols : MaxRows, + inner_max_size = MaxCols==1 ? MaxRows + : MaxRows==1 ? MaxCols + : row_major_bit ? MaxCols : MaxRows, is_big = inner_max_size == Dynamic, - is_packet_size_multiple = (Cols*Rows) % ei_packet_traits::size == 0, + is_packet_size_multiple = Rows==Dynamic || Cols==Dynamic || ((Cols*Rows) % ei_packet_traits::size) == 0, aligned_bit = (((Options&DontAlign)==0) && (is_big || is_packet_size_multiple)) ? AlignedBit : 0, packet_access_bit = ei_packet_traits::size > 1 && aligned_bit ? PacketAccessBit : 0 }; diff --git a/Eigen/src/QR/Tridiagonalization.h b/Eigen/src/QR/Tridiagonalization.h index bd8ff4fe3..59bffc12a 100644 --- a/Eigen/src/QR/Tridiagonalization.h +++ b/Eigen/src/QR/Tridiagonalization.h @@ -279,6 +279,7 @@ Tridiagonalization::matrixQ(void) const Scalar tmp = m_matrix.coeff(i+1,i); m_matrix.const_cast_derived().coeffRef(i+1,i) = 1; + // TODO this product could be optimized by processing the submatrix per panel of at least 4 columns matQ.corner(BottomRight,n-i-1,n-i-1) -= ((m_hCoeffs.coeff(i) * m_matrix.col(i).end(n-i-1)) * (m_matrix.col(i).end(n-i-1).adjoint() * matQ.corner(BottomRight,n-i-1,n-i-1)).lazy()).lazy(); diff --git a/test/adjoint.cpp b/test/adjoint.cpp index 403f16a45..1f4aa7427 100644 --- a/test/adjoint.cpp +++ b/test/adjoint.cpp @@ -76,7 +76,8 @@ template void adjoint(const MatrixType& m) { VERIFY_IS_MUCH_SMALLER_THAN(vzero.norm(), static_cast(1)); VERIFY_IS_APPROX(v1.norm(), v1.stableNorm()); - VERIFY_IS_APPROX(v1.blueNorm(), v1.stableNorm()); + // NOTE disabled because it currently compiles for float and double only + // VERIFY_IS_APPROX(v1.blueNorm(), v1.stableNorm()); } // check compatibility of dot and adjoint @@ -109,7 +110,7 @@ template void adjoint(const MatrixType& m) m3.transposeInPlace(); VERIFY_IS_APPROX(m3,m1.conjugate()); - + } void test_adjoint()