diff --git a/Eigen/src/Core/CacheFriendlyProduct.h b/Eigen/src/Core/CacheFriendlyProduct.h index 3ce72f5ba..a9ef66cdd 100644 --- a/Eigen/src/Core/CacheFriendlyProduct.h +++ b/Eigen/src/Core/CacheFriendlyProduct.h @@ -86,13 +86,15 @@ static void ei_cache_friendly_product( const int l2BlockRows = MaxL2BlockSize > rows ? rows : MaxL2BlockSize; const int l2BlockCols = MaxL2BlockSize > cols ? cols : MaxL2BlockSize; const int l2BlockSize = MaxL2BlockSize > size ? size : MaxL2BlockSize; + const int l2BlockSizeAligned = (1 + std::max(l2BlockSize,l2BlockCols)/PacketSize)*PacketSize; + const bool needRhsCopy = (PacketSize>1) && ((rhsStride%PacketSize!=0) || (size_t(rhs)%16!=0)); Scalar* __restrict__ block = 0; const int allocBlockSize = sizeof(Scalar)*l2BlockRows*size; if (allocBlockSize>16000000) block = (Scalar*)malloc(allocBlockSize); else block = (Scalar*)alloca(allocBlockSize); - Scalar* __restrict__ rhsCopy = (Scalar*)alloca(sizeof(Scalar)*l2BlockSize); + Scalar* __restrict__ rhsCopy = (Scalar*)alloca(sizeof(Scalar)*l2BlockSizeAligned*l2BlockSizeAligned); // loops on each L2 cache friendly blocks of the result for(int l2i=0; l2i1 && size_t(rhsColumn)%16) - { - int count = 0; - // FIXME this loop get vectorized by the compiler (ICC) - // I'm not sure thats good or not - for (int k = l2k; k1 && resIsAligned) { - ei_pstore(&(localRes[0]), ei_padd(ei_pload(&(localRes[0])), ei_preduxp(dst))); + // the result is aligned: let's do packet reduction + ei_pstore(&(localRes[0]), ei_padd(ei_pload(&(localRes[0])), ei_preduxp(&dst[0]))); if (PacketSize==2) ei_pstore(&(localRes[2]), ei_padd(ei_pload(&(localRes[2])), ei_preduxp(&(dst[2])))); if (MaxBlockRows==8) @@ -239,6 +238,7 @@ static void ei_cache_friendly_product( } else { + // not aligned => per coeff packet reduction localRes[0] += ei_predux(dst[0]); localRes[1] += ei_predux(dst[1]); localRes[2] += ei_predux(dst[2]); @@ -262,32 +262,16 @@ static void ei_cache_friendly_product( asm("#eigen begin dynkernel"); for(int l1j=l2j; l1j1 && size_t(rhsColumn)%16) - { - int count = 0; - // FIXME this loop get vectorized by the compiler ! - for (int k = l2k; k4) - { - dst[4] = dst[0]; - dst[5] = dst[0]; - dst[6] = dst[0]; - dst[7] = dst[0]; - } + dst[3] = dst[2] = dst[1] = dst[0] = ei_pset1(Scalar(0.)); + if (MaxBlockRows==8) + dst[7] = dst[6] = dst[5] = dst[4] = dst[0]; // let's declare a few other temporary registers PacketType tmp; @@ -300,7 +284,7 @@ static void ei_cache_friendly_product( if (l2blockRemainingRows>=2) dst[1] = ei_pmadd(tmp, ei_pload(&(localB[k*l2blockRemainingRows+ PacketSize])), dst[1]); if (l2blockRemainingRows>=3) dst[2] = ei_pmadd(tmp, ei_pload(&(localB[k*l2blockRemainingRows+2*PacketSize])), dst[2]); if (l2blockRemainingRows>=4) dst[3] = ei_pmadd(tmp, ei_pload(&(localB[k*l2blockRemainingRows+3*PacketSize])), dst[3]); - if (MaxBlockRows>4) + if (MaxBlockRows==8) { if (l2blockRemainingRows>=5) dst[4] = ei_pmadd(tmp, ei_pload(&(localB[k*l2blockRemainingRows+4*PacketSize])), dst[4]); if (l2blockRemainingRows>=6) dst[5] = ei_pmadd(tmp, ei_pload(&(localB[k*l2blockRemainingRows+5*PacketSize])), dst[5]); @@ -316,7 +300,7 @@ static void ei_cache_friendly_product( if (l2blockRemainingRows>=2) localRes[1] += ei_predux(dst[1]); if (l2blockRemainingRows>=3) localRes[2] += ei_predux(dst[2]); if (l2blockRemainingRows>=4) localRes[3] += ei_predux(dst[3]); - if (MaxBlockRows>4) + if (MaxBlockRows==8) { if (l2blockRemainingRows>=5) localRes[4] += ei_predux(dst[4]); if (l2blockRemainingRows>=6) localRes[5] += ei_predux(dst[5]); @@ -573,16 +557,16 @@ EIGEN_DONT_INLINE static void ei_cache_friendly_product_rowmajor_times_vector( enum { AllAligned, EvenAligned, FirstAligned, NoneAligned }; const int rowsAtOnce = 4; -// const int peels = 2; + const int peels = 2; const int PacketAlignedMask = PacketSize-1; -// const int PeelAlignedMask = PacketSize*peels-1; + const int PeelAlignedMask = PacketSize*peels-1; const int size = rhsSize; // How many coeffs of the result do we have to skip to be aligned. // Here we assume data are at least aligned on the base scalar type that is mandatory anyway. const int alignedStart = ei_alignmentOffset(rhs, size); const int alignedSize = PacketSize>1 ? alignedStart + ((size-alignedStart) & ~PacketAlignedMask) : 0; - //const int peeledSize = peels>1 ? alignedStart + ((alignedSize-alignedStart) & ~PeelAlignedMask) : 0; + const int peeledSize = peels>1 ? alignedStart + ((alignedSize-alignedStart) & ~PeelAlignedMask) : alignedStart; const int alignmentStep = PacketSize>1 ? (PacketSize - lhsStride % PacketSize) & PacketAlignedMask : 0; int alignmentPattern = alignmentStep==0 ? AllAligned @@ -650,7 +634,32 @@ EIGEN_DONT_INLINE static void ei_cache_friendly_product_rowmajor_times_vector( _EIGEN_ACCUMULATE_PACKETS(,u,,); break; case FirstAligned: - for (int j = alignedStart; j1) + { + Packet A01, A02, A03, b; + for (int j = alignedStart; j > ColsAtCompileTime = MatrixType::ColsAtCompileTime, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, - Flags = MatrixType::Flags - & ( (HereditaryBits | LinearAccessBit | DirectAccessBit) - | (Alignment == Aligned ? PacketAccessBit : 0) ), + Flags = MatrixType::Flags, CoeffReadCost = NumTraits::ReadCost }; }; diff --git a/Eigen/src/Core/arch/SSE/PacketMath.h b/Eigen/src/Core/arch/SSE/PacketMath.h index 9deb19d32..28825ceee 100644 --- a/Eigen/src/Core/arch/SSE/PacketMath.h +++ b/Eigen/src/Core/arch/SSE/PacketMath.h @@ -192,12 +192,12 @@ inline __m128i ei_preduxp(const __m128i* vecs) } #if (defined __GNUC__) -template <> inline __m128 ei_pmadd(const __m128& a, const __m128& b, const __m128& c) -{ - __m128 res = b; - asm("mulps %[a], %[b] \n\taddps %[c], %[b]" : [b] "+x" (res) : [a] "x" (a), [c] "x" (c)); - return res; -} +// template <> inline __m128 ei_pmadd(const __m128& a, const __m128& b, const __m128& c) +// { +// __m128 res = b; +// asm("mulps %[a], %[b] \n\taddps %[c], %[b]" : [b] "+x" (res) : [a] "x" (a), [c] "x" (c)); +// return res; +// } #endif #endif // EIGEN_PACKET_MATH_SSE_H diff --git a/Eigen/src/QR/Tridiagonalization.h b/Eigen/src/QR/Tridiagonalization.h index 7834c1aca..765a87130 100755 --- a/Eigen/src/QR/Tridiagonalization.h +++ b/Eigen/src/QR/Tridiagonalization.h @@ -34,7 +34,7 @@ * \param MatrixType the type of the matrix of which we are performing the tridiagonalization * * This class performs a tridiagonal decomposition of a selfadjoint matrix \f$ A \f$ such that: - * \f$ A = Q T Q^* \f$ where \f$ Q \f$ is unitatry and \f$ T \f$ a real symmetric tridiagonal matrix + * \f$ A = Q T Q^* \f$ where \f$ Q \f$ is unitary and \f$ T \f$ a real symmetric tridiagonal matrix. * * \sa MatrixBase::tridiagonalize() */ @@ -45,12 +45,15 @@ template class Tridiagonalization typedef _MatrixType MatrixType; typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; + typedef typename ei_packet_traits::type Packet; enum { Size = MatrixType::RowsAtCompileTime, SizeMinusOne = MatrixType::RowsAtCompileTime==Dynamic - ? Dynamic - : MatrixType::RowsAtCompileTime-1}; + ? Dynamic + : MatrixType::RowsAtCompileTime-1, + PacketSize = ei_packet_traits::size + }; typedef Matrix CoeffVectorType; typedef Matrix DiagonalType; @@ -59,8 +62,7 @@ template class Tridiagonalization typedef typename NestByValue >::RealReturnType DiagonalReturnType; typedef typename NestByValue > > >::RealReturnType SubDiagonalReturnType; + NestByValue > > >::RealReturnType SubDiagonalReturnType; /** This constructor initializes a Tridiagonalization object for * further use with Tridiagonalization::compute() @@ -103,7 +105,7 @@ template class Tridiagonalization * Householder coefficients returned by householderCoefficients(), * allows to reconstruct the matrix Q as follow: * Q = H_{N-1} ... H_1 H_0 - * where the matrices H are the Householder transformation: + * where the matrices H are the Householder transformations: * H_i = (I - h_i * v_i * v_i') * where h_i == householderCoefficients()[i] and v_i is a Householder vector: * v_i = [ 0, ..., 0, 1, M(i+2,i), ..., M(N-1,i) ] @@ -157,8 +159,8 @@ template typename Tridiagonalization::MatrixType Tridiagonalization::matrixT(void) const { - // FIXME should this function (and other similar) rather take a matrix as argument - // and fill it (avoids temporaries) + // FIXME should this function (and other similar ones) rather take a matrix as argument + // and fill it ? (to avoid temporaries) int n = m_matrix.rows(); MatrixType matT = m_matrix; matT.corner(TopRight,n-1, n-1).diagonal() = subDiagonal().conjugate(); @@ -189,6 +191,7 @@ void Tridiagonalization::_compute(MatrixType& matA, CoeffVectorType& { assert(matA.rows()==matA.cols()); int n = matA.rows(); +// std::cerr << matA << "\n\n"; for (int i = 0; i::_compute(MatrixType& matA, CoeffVectorType& // i.e., A = H' A H where H = I - h v v' and v = matA.col(i).end(n-i-1) matA.col(i).coeffRef(i+1) = 1; - // let's use the end of hCoeffs to store temporary values - hCoeffs.end(n-i-1) = h * (matA.corner(BottomRight,n-i-1,n-i-1).template part() - * matA.col(i).end(n-i-1)); + + /* This is the initial algorithm which minimize operation counts and maximize + * the use of Eigen's expression. Unfortunately, the first matrix-vector product + * using Part is very very slow */ + #ifdef EIGEN_NEVER_DEFINED + // matrix - vector product + hCoeffs.end(n-i-1) = (matA.corner(BottomRight,n-i-1,n-i-1).template part() + * (h * matA.col(i).end(n-i-1))).lazy(); + // simple axpy + hCoeffs.end(n-i-1) += (h * Scalar(-0.5) * matA.col(i).end(n-i-1).dot(hCoeffs.end(n-i-1))) + * matA.col(i).end(n-i-1); + // rank-2 update + //Block B(matA,i+1,i,n-i-1,1); + matA.corner(BottomRight,n-i-1,n-i-1).template part() -= + (matA.col(i).end(n-i-1) * hCoeffs.end(n-i-1).adjoint()).lazy() + + (hCoeffs.end(n-i-1) * matA.col(i).end(n-i-1).adjoint()).lazy(); + #endif + /* end initial algorithm */ + /* If we still want to minimize operation count (i.e., perform operation on the lower part only) + * then we could provide the following algorithm for selfadjoint - vector product. However, a full + * matrix-vector product is still faster (at least for dynamic size, and not too small, did not check + * small matrices). The algo performs block matrix-vector and transposed matrix vector products. */ + #ifdef EIGEN_NEVER_DEFINED + int n4 = (std::max(0,n-4)/4)*4; + hCoeffs.end(n-i-1).setZero(); + for (int b=i+1; b(matA,b+4,b,n-b-4,4) * matA.template block<4,1>(b,i); + // the respective transposed part: + Block(hCoeffs, b, 0, 4,1) += + Block(matA,b+4,b,n-b-4,4).adjoint() * Block(matA,b+4,i,n-b-4,1); + // the 4x4 block diagonal: + Block(hCoeffs, b, 0, 4,1) += + (Block(matA,b,b,4,4).template part() + * (h * Block(matA,b,i,4,1))).lazy(); + } + #endif + // todo: handle the remaining part + /* end optimized selfadjoint - vector product */ + + /* Another interesting note: the above rank-2 update is much slower than the following hand written loop. + * After an analyse of the ASM, it seems GCC (4.2) generate poor code because of the Block. Moreover, + * if we remove the specialization of Block for Matrix then it is even worse, much worse ! */ + #ifdef EIGEN_NEVER_DEFINED + for (int j1=i+1; j1() = - matA.corner(BottomRight,n-i-1,n-i-1) - ( - (matA.col(i).end(n-i-1) * hCoeffs.end(n-i-1).adjoint()).lazy() - + (hCoeffs.end(n-i-1) * matA.col(i).end(n-i-1).adjoint()).lazy() ); - // FIXME check that the above expression does follow the lazy path (no temporary and - // only lower products are evaluated) - // FIXME can we avoid to evaluate twice the diagonal products ? - // (in a simple way otherwise it's overkill) + + const Scalar* __restrict__ pb = &matA.coeffRef(0,i); + const Scalar* __restrict__ pa = (&hCoeffs.coeffRef(0)) - 1; + for (int j1=i+1; j11) + { + int alignedStart = (starti) + ei_alignmentOffset(&matA.coeffRef(starti,j1), n-starti); + alignedEnd = alignedStart + ((n-alignedStart)/PacketSize)*PacketSize; + + for (int i1=starti; i1::decomposeInPlace(MatrixType& mat, DiagonalT ei_assert(mat.cols()==n && diag.size()==n && subdiag.size()==n-1); if (n==3 && (!NumTraits::IsComplex) ) { - Tridiagonalization tridiag(mat); _decomposeInPlace3x3(mat, diag, subdiag, extractQ); } else @@ -301,7 +381,7 @@ void Tridiagonalization::decomposeInPlace(MatrixType& mat, DiagonalT /** \internal * Optimized path for 3x3 matrices. - * Especially usefull for plane fit. + * Especially useful for plane fitting. */ template void Tridiagonalization::_decomposeInPlace3x3(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 14ed29a3d..97c7f4937 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -38,7 +38,12 @@ MACRO(EI_ADD_TEST testname) SET(targetname test_${testname}) - ADD_EXECUTABLE(${targetname} ${testname}.cpp) + IF(${ARGC} EQUAL 2) + SET(filename ${ARGV1}) + ELSE(${ARGC} EQUAL 2) + SET(filename ${testname}.cpp) + ENDIF(${ARGC} EQUAL 2) + ADD_EXECUTABLE(${targetname} ${filename}) IF(NOT EIGEN_NO_ASSERTION_CHECKING) @@ -80,7 +85,8 @@ EI_ADD_TEST(nomalloc) EI_ADD_TEST(basicstuff) EI_ADD_TEST(linearstructure) EI_ADD_TEST(cwiseop) -EI_ADD_TEST(product) +EI_ADD_TEST(product_small) +EI_ADD_TEST(product_large) EI_ADD_TEST(adjoint) EI_ADD_TEST(submatrices) EI_ADD_TEST(miscmatrices) diff --git a/test/eigensolver.cpp b/test/eigensolver.cpp index 9837162f6..a1ab4a685 100644 --- a/test/eigensolver.cpp +++ b/test/eigensolver.cpp @@ -63,8 +63,8 @@ void test_eigensolver() // very important to test a 3x3 matrix since we provide a special path for it CALL_SUBTEST( eigensolver(Matrix3f()) ); CALL_SUBTEST( eigensolver(Matrix4d()) ); - CALL_SUBTEST( eigensolver(MatrixXd(7,7)) ); + CALL_SUBTEST( eigensolver(MatrixXf(7,7)) ); CALL_SUBTEST( eigensolver(MatrixXcd(6,6)) ); - CALL_SUBTEST( eigensolver(MatrixXcd(3,3)) ); + CALL_SUBTEST( eigensolver(MatrixXcf(3,3)) ); } } diff --git a/test/product.cpp b/test/product.h similarity index 85% rename from test/product.cpp rename to test/product.h index 2f6677ff1..374994576 100644 --- a/test/product.cpp +++ b/test/product.h @@ -144,25 +144,3 @@ template void product(const MatrixType& m) } } -void test_product() -{ - for(int i = 0; i < g_repeat; i++) { - CALL_SUBTEST( product(Matrix3i()) ); - CALL_SUBTEST( product(Matrix()) ); - CALL_SUBTEST( product(Matrix4d()) ); - CALL_SUBTEST( product(Matrix4f()) ); - CALL_SUBTEST( product(MatrixXf(3,5)) ); - CALL_SUBTEST( product(MatrixXi(28,39)) ); - } - for(int i = 0; i < g_repeat; i++) { - CALL_SUBTEST( product(MatrixXf(ei_random(1,320), ei_random(1,320))) ); - CALL_SUBTEST( product(MatrixXd(ei_random(1,320), ei_random(1,320))) ); - CALL_SUBTEST( product(MatrixXi(ei_random(1,256), ei_random(1,256))) ); - CALL_SUBTEST( product(MatrixXcf(ei_random(1,50), ei_random(1,50))) ); - #ifndef EIGEN_DEFAULT_TO_ROW_MAJOR - CALL_SUBTEST( product(Matrix(ei_random(1,320), ei_random(1,320))) ); - #else - CALL_SUBTEST( product(Matrix(ei_random(1,320), ei_random(1,320))) ); - #endif - } -} diff --git a/test/product_large.cpp b/test/product_large.cpp new file mode 100644 index 000000000..904cf5a0b --- /dev/null +++ b/test/product_large.cpp @@ -0,0 +1,40 @@ +// 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 +// +// 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 . + +#include "product.h" + +void test_product_large() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST( product(MatrixXf(ei_random(1,320), ei_random(1,320))) ); + CALL_SUBTEST( product(MatrixXd(ei_random(1,320), ei_random(1,320))) ); + CALL_SUBTEST( product(MatrixXi(ei_random(1,320), ei_random(1,320))) ); + CALL_SUBTEST( product(MatrixXcf(ei_random(1,50), ei_random(1,50))) ); + #ifndef EIGEN_DEFAULT_TO_ROW_MAJOR + CALL_SUBTEST( product(Matrix(ei_random(1,320), ei_random(1,320))) ); + #else + CALL_SUBTEST( product(Matrix(ei_random(1,320), ei_random(1,320))) ); + #endif + } +} diff --git a/test/product_small.cpp b/test/product_small.cpp new file mode 100644 index 000000000..ef44b0826 --- /dev/null +++ b/test/product_small.cpp @@ -0,0 +1,35 @@ +// 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 +// +// 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 . + +#include "product.h" + +void test_product_small() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST( product(Matrix()) ); + CALL_SUBTEST( product(Matrix()) ); + CALL_SUBTEST( product(Matrix4d()) ); + CALL_SUBTEST( product(Matrix4f()) ); + } +}