diff --git a/Eigen/src/Core/Block.h b/Eigen/src/Core/Block.h index c8f4ac6b1..c026b174c 100644 --- a/Eigen/src/Core/Block.h +++ b/Eigen/src/Core/Block.h @@ -67,9 +67,11 @@ struct ei_traits > : (BlockRows==Dynamic ? MatrixType::MaxRowsAtCompileTime : BlockRows), MaxColsAtCompileTime = ColsAtCompileTime == 1 ? 1 : (BlockCols==Dynamic ? MatrixType::MaxColsAtCompileTime : BlockCols), - Flags = (RowsAtCompileTime == Dynamic || ColsAtCompileTime == Dynamic - ? (unsigned int)MatrixType::Flags - : (unsigned int)MatrixType::Flags &~ LargeBit) & ~VectorizableBit, + FlagsLargeBit = ((RowsAtCompileTime != Dynamic && MatrixType::RowsAtCompileTime == Dynamic) + || (ColsAtCompileTime != Dynamic && MatrixType::ColsAtCompileTime == Dynamic)) + ? ~LargeBit + : ~(unsigned int)0, + Flags = MatrixType::Flags & FlagsLargeBit & ~VectorizableBit, CoeffReadCost = MatrixType::CoeffReadCost }; }; @@ -236,8 +238,7 @@ const Block MatrixBase * \sa class Block, block(int,int) */ template -Block MatrixBase - ::start(int size) +Block MatrixBase::start(int size) { ei_assert(IsVectorAtCompileTime); return Block(derived(), 0, 0, @@ -247,8 +248,7 @@ Block MatrixBase /** This is the const version of start(int).*/ template -const Block MatrixBase - ::start(int size) const +const Block MatrixBase::start(int size) const { ei_assert(IsVectorAtCompileTime); return Block(derived(), 0, 0, @@ -272,8 +272,7 @@ const Block MatrixBase * \sa class Block, block(int,int) */ template -Block MatrixBase - ::end(int size) +Block MatrixBase::end(int size) { ei_assert(IsVectorAtCompileTime); return Block(derived(), @@ -285,8 +284,7 @@ Block MatrixBase /** This is the const version of end(int).*/ template -const Block MatrixBase - ::end(int size) const +const Block MatrixBase::end(int size) const { ei_assert(IsVectorAtCompileTime); return Block(derived(), @@ -296,6 +294,80 @@ const Block MatrixBase ColsAtCompileTime == 1 ? 1 : size); } +/** \returns a fixed-size expression of the first coefficients of *this. + * + * \only_for_vectors + * + * The template parameter \a Size is the number of coefficients in the block + * + * Example: \include MatrixBase_template_int_start.cpp + * Output: \verbinclude MatrixBase_template_int_start.out + * + * \sa class Block + */ +template +template +Block::RowsAtCompileTime == 1 ? 1 : Size, + ei_traits::ColsAtCompileTime == 1 ? 1 : Size> +MatrixBase::start() +{ + ei_assert(IsVectorAtCompileTime); + return Block(derived(), 0, 0); +} + +/** This is the const version of start().*/ +template +template +const Block::RowsAtCompileTime == 1 ? 1 : Size, + ei_traits::ColsAtCompileTime == 1 ? 1 : Size> +MatrixBase::start() const +{ + ei_assert(IsVectorAtCompileTime); + return Block(derived(), 0, 0); +} + +/** \returns a fixed-size expression of the last coefficients of *this. + * + * \only_for_vectors + * + * The template parameter \a Size is the number of coefficients in the block + * + * Example: \include MatrixBase_template_int_end.cpp + * Output: \verbinclude MatrixBase_template_int_end.out + * + * \sa class Block + */ +template +template +Block::RowsAtCompileTime == 1 ? 1 : Size, + ei_traits::ColsAtCompileTime == 1 ? 1 : Size> +MatrixBase::end() +{ + ei_assert(IsVectorAtCompileTime); + return Block + (derived(), + RowsAtCompileTime == 1 ? 0 : rows() - Size, + ColsAtCompileTime == 1 ? 0 : cols() - Size); +} + +/** This is the const version of end.*/ +template +template +const Block::RowsAtCompileTime == 1 ? 1 : Size, + ei_traits::ColsAtCompileTime == 1 ? 1 : Size> +MatrixBase::end() const +{ + ei_assert(IsVectorAtCompileTime); + return Block + (derived(), + RowsAtCompileTime == 1 ? 0 : rows() - Size, + ColsAtCompileTime == 1 ? 0 : cols() - Size); +} + /** \returns a dynamic-size expression of a corner of *this. * * \param type the type of corner. Can be \a Eigen::TopLeft, \a Eigen::TopRight, @@ -316,14 +388,23 @@ template Block MatrixBase ::corner(CornerType type, int cRows, int cCols) { - if(type == TopLeft) - return Block(derived(), 0, 0, cRows, cCols); - else if(type == TopRight) - return Block(derived(), 0, cols() - cCols, cRows, cCols); - else if(type == BottomLeft) - return Block(derived(), rows() - cRows, 0, cRows, cCols); - else - return Block(derived(), rows() - cRows, cols() - cCols, cRows, cCols); + switch(type) + { + case TopLeft: + return Block(derived(), 0, 0, cRows, cCols); + break; + case TopRight: + return Block(derived(), 0, cols() - cCols, cRows, cCols); + break; + case BottomLeft: + return Block(derived(), rows() - cRows, 0, cRows, cCols); + break; + case BottomRight: + return Block(derived(), rows() - cRows, cols() - cCols, cRows, cCols); + break; + default: + ei_assert(false && "Bad corner type."); + } } /** This is the const version of corner(CornerType, int, int).*/ @@ -331,14 +412,84 @@ template const Block MatrixBase ::corner(CornerType type, int cRows, int cCols) const { - if(type == TopLeft) - return Block(derived(), 0, 0, cRows, cCols); - else if(type == TopRight) - return Block(derived(), 0, cols() - cCols, cRows, cCols); - else if(type == BottomLeft) - return Block(derived(), rows() - cRows, 0, cRows, cCols); - else - return Block(derived(), rows() - cRows, cols() - cCols, cRows, cCols); + switch(type) + { + case TopLeft: + return Block(derived(), 0, 0, cRows, cCols); + break; + case TopRight: + return Block(derived(), 0, cols() - cCols, cRows, cCols); + break; + case BottomLeft: + return Block(derived(), rows() - cRows, 0, cRows, cCols); + break; + case BottomRight: + return Block(derived(), rows() - cRows, cols() - cCols, cRows, cCols); + break; + default: + ei_assert(false && "Bad corner type."); + } +} + +/** \returns a fixed-size expression of a corner of *this. + * + * \param type the type of corner. Can be \a Eigen::TopLeft, \a Eigen::TopRight, + * \a Eigen::BottomLeft, \a Eigen::BottomRight. + * + * The template parameters CRows and CCols arethe number of rows and columns in the corner. + * + * Example: \include MatrixBase_template_int_int_corner_enum.cpp + * Output: \verbinclude MatrixBase_template_int_int_corner_enum.out + * + * \sa class Block, block(int,int,int,int) + */ +template +template +Block MatrixBase + ::corner(CornerType type) +{ + switch(type) + { + case TopLeft: + return Block(derived(), 0, 0); + break; + case TopRight: + return Block(derived(), 0, cols() - CCols); + break; + case BottomLeft: + return Block(derived(), rows() - CRows, 0); + break; + case BottomRight: + return Block(derived(), rows() - CRows, cols() - CCols); + break; + default: + ei_assert(false && "Bad corner type."); + } +} + +/** This is the const version of corner(CornerType).*/ +template +template +const Block MatrixBase + ::corner(CornerType type) const +{ + switch(type) + { + case TopLeft: + return Block(derived(), 0, 0); + break; + case TopRight: + return Block(derived(), 0, cols() - CCols); + break; + case BottomLeft: + return Block(derived(), rows() - CRows, 0); + break; + case BottomRight: + return Block(derived(), rows() - CRows, cols() - CCols); + break; + default: + ei_assert(false && "Bad corner type."); + } } /** \returns a fixed-size expression of a block in *this. diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index f2d665872..72cf916f8 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h @@ -28,6 +28,11 @@ template inline typename NumTraits::Real precision(); template inline T ei_random(T a, T b); template inline T ei_random(); +template inline T ei_random_amplitude() +{ + if(NumTraits::HasFloatingPoint) return static_cast(1); + else return static_cast(10); +} template<> inline int precision() { return 0; } inline int ei_real(int x) { return x; } @@ -54,7 +59,7 @@ template<> inline int ei_random(int a, int b) } template<> inline int ei_random() { - return ei_random(-10, 10); + return ei_random(-ei_random_amplitude(), ei_random_amplitude()); } inline bool ei_isMuchSmallerThan(int a, int, int = precision()) { @@ -88,7 +93,7 @@ template<> inline float ei_random(float a, float b) } template<> inline float ei_random() { - return ei_random(-10.0f, 10.0f); + return ei_random(-ei_random_amplitude(), ei_random_amplitude()); } inline bool ei_isMuchSmallerThan(float a, float b, float prec = precision()) { @@ -122,7 +127,7 @@ template<> inline double ei_random(double a, double b) } template<> inline double ei_random() { - return ei_random(-10.0, 10.0); + return ei_random(-ei_random_amplitude(), ei_random_amplitude()); } inline bool ei_isMuchSmallerThan(double a, double b, double prec = precision()) { diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 5c2350d6c..151ee74d9 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -316,6 +316,23 @@ template class MatrixBase template const Block block(int startRow, int startCol) const; + template Block corner(CornerType type); + template const Block corner(CornerType type) const; + + template + Block::RowsAtCompileTime == 1 ? 1 : Size, + ei_traits::ColsAtCompileTime == 1 ? 1 : Size> start(); + template + const Block::RowsAtCompileTime == 1 ? 1 : Size, + ei_traits::ColsAtCompileTime == 1 ? 1 : Size> start() const; + + template + Block::RowsAtCompileTime == 1 ? 1 : Size, + ei_traits::ColsAtCompileTime == 1 ? 1 : Size> end(); + template + const Block::RowsAtCompileTime == 1 ? 1 : Size, + ei_traits::ColsAtCompileTime == 1 ? 1 : Size> end() const; + DiagonalCoeffs diagonal(); const DiagonalCoeffs diagonal() const; //@} diff --git a/Eigen/src/Core/PacketMath.h b/Eigen/src/Core/PacketMath.h index 3697f262e..cc925b50c 100644 --- a/Eigen/src/Core/PacketMath.h +++ b/Eigen/src/Core/PacketMath.h @@ -79,7 +79,7 @@ inline float ei_pfirst(const __m128& a) { return _mm_cvtss_f32(a); } inline double ei_pfirst(const __m128d& a) { return _mm_cvtsd_f64(a); } inline int ei_pfirst(const __m128i& a) { return _mm_cvtsi128_si32(a); } -#endif +#endif // EIGEN_VECTORIZE_SSE #endif // EIGEN_PACKET_MATH_H diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h index 2f39b48fd..43d595451 100644 --- a/Eigen/src/Core/util/Macros.h +++ b/Eigen/src/Core/util/Macros.h @@ -142,6 +142,7 @@ EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, /=) #define _EIGEN_GENERIC_PUBLIC_INTERFACE(Derived, BaseClass) \ typedef BaseClass Base; \ typedef typename Eigen::ei_traits::Scalar Scalar; \ +typedef typename Eigen::NumTraits::Real RealScalar; \ typedef typename Base::PacketScalar PacketScalar; \ typedef typename Eigen::ei_nested::type Nested; \ typedef typename Eigen::ei_eval::type Eval; \ diff --git a/bench/benchmark.cpp b/bench/benchmark.cpp index 39ed317aa..418bf72f2 100644 --- a/bench/benchmark.cpp +++ b/bench/benchmark.cpp @@ -28,7 +28,7 @@ int main(int argc, char *argv[]) asm("#begin"); for(int a = 0; a < REPEAT; a++) { - m = I + 0.00005 * (m + m*m); + m = Matrix::ones() + 0.00005 * (m + m*m); } asm("#end"); cout << m << endl; diff --git a/bench/benchmarkXcwise.cpp b/bench/benchmarkXcwise.cpp index dd29743cd..b2a7fc24c 100644 --- a/bench/benchmarkXcwise.cpp +++ b/bench/benchmarkXcwise.cpp @@ -10,24 +10,24 @@ USING_PART_OF_NAMESPACE_EIGEN #endif #ifndef MATSIZE -#define MATSIZE 400 +#define MATSIZE 1000000 #endif #ifndef REPEAT -#define REPEAT 10000 +#define REPEAT 1000 #endif int main(int argc, char *argv[]) { - MATTYPE I = MATTYPE::ones(MATSIZE,MATSIZE); - MATTYPE m(MATSIZE,MATSIZE); - for(int i = 0; i < MATSIZE; i++) for(int j = 0; j < MATSIZE; j++) + MATTYPE I = MATTYPE::ones(MATSIZE,1); + MATTYPE m(MATSIZE,1); + for(int i = 0; i < MATSIZE; i++) for(int j = 0; j < 1; j++) { - m(i,j) = 0.1 * (i+j+1)/(MATSIZE*MATSIZE); + m(i,j) = 0.1 * (i+j+1)/MATSIZE/MATSIZE; } for(int a = 0; a < REPEAT; a++) { - m = I + 0.00005 * (m + m/4); + m = MATTYPE::ones(MATSIZE,1) + 0.00005 * (m.cwiseProduct(m) + m/4); } cout << m(0,0) << endl; return 0; diff --git a/doc/echelon.cpp b/doc/echelon.cpp index e305eb238..04b1907cd 100644 --- a/doc/echelon.cpp +++ b/doc/echelon.cpp @@ -4,21 +4,72 @@ USING_PART_OF_NAMESPACE_EIGEN namespace Eigen { -template -void echelon(MatrixBase& m) +/* Echelon a matrix in-place: + * + * Meta-Unrolled version, for small fixed-size matrices + */ +template +struct unroll_echelon { - for(int k = 0; k < m.diagonal().size(); k++) + enum { k = Step - 1, + Rows = Derived::RowsAtCompileTime, + Cols = Derived::ColsAtCompileTime, + CornerRows = Rows - k, + CornerCols = Cols - k + }; + static void run(MatrixBase& m) { + unroll_echelon::run(m); int rowOfBiggest, colOfBiggest; - int cornerRows = m.rows()-k, cornerCols = m.cols()-k; - m.corner(BottomRight, cornerRows, cornerCols) + m.template corner(BottomRight) .cwiseAbs() .maxCoeff(&rowOfBiggest, &colOfBiggest); m.row(k).swap(m.row(k+rowOfBiggest)); m.col(k).swap(m.col(k+colOfBiggest)); - m.corner(BottomRight, cornerRows-1, cornerCols) - -= m.col(k).end(cornerRows-1) * (m.row(k).end(cornerCols) / m(k,k)); + m.template corner(BottomRight) + -= m.col(k).template end() + * (m.row(k).template end() / m(k,k)); } +}; + +template +struct unroll_echelon +{ + static void run(MatrixBase& m) {} +}; + +/* Echelon a matrix in-place: + * + * Non-unrolled version, for dynamic-size matrices. + * (this version works for all matrices, but in the fixed-size case the other + * version is faster). + */ +template +struct unroll_echelon +{ + static void run(MatrixBase& m) + { + for(int k = 0; k < m.diagonal().size(); k++) + { + int rowOfBiggest, colOfBiggest; + int cornerRows = m.rows()-k, cornerCols = m.cols()-k; + m.corner(BottomRight, cornerRows, cornerCols) + .cwiseAbs() + .maxCoeff(&rowOfBiggest, &colOfBiggest); + m.row(k).swap(m.row(k+rowOfBiggest)); + m.col(k).swap(m.col(k+colOfBiggest)); + m.corner(BottomRight, cornerRows-1, cornerCols) + -= m.col(k).end(cornerRows-1) * (m.row(k).end(cornerCols) / m(k,k)); + } + } +}; +using namespace std; +template +void echelon(MatrixBase& m) +{ + const int size = DiagonalCoeffs::SizeAtCompileTime; + const bool unroll = size <= 4; + unroll_echelon::run(m); } template @@ -49,7 +100,7 @@ using namespace std; int main(int, char **) { srand((unsigned int)time(0)); - const int Rows = 6, Cols = 4; + const int Rows = 6, Cols = 6; typedef Matrix Mat; const int N = Rows < Cols ? Rows : Cols;