Compare commits

..

105 Commits

Author SHA1 Message Date
Benoit Jacob
907fba9ea9 bump 2010-07-16 22:25:08 -04:00
Gael Guennebaud
7e0b7b1f25 uncomment tests (sorry) 2010-07-16 11:50:40 +02:00
Gael Guennebaud
5ac00624ed adapted pactch from Nick Lewycky fixing a couple of issue with CLang 2010-07-15 22:48:44 +02:00
Gael Guennebaud
c0a0d2d181 add unit tests for 0 matrix 2010-07-15 20:00:50 +02:00
Gael Guennebaud
ec39a39cb6 fix LLT for zero matrix 2010-07-15 20:00:34 +02:00
Gael Guennebaud
ccc6731f86 fix use of rank in QR 2010-07-15 19:59:21 +02:00
Gael Guennebaud
b09bb50aeb fix QR solving with m>n 2010-07-15 19:52:11 +02:00
Gael Guennebaud
c7b8de77c0 fix aligned_delete for null pointers 2010-07-15 09:28:29 +02:00
Gael Guennebaud
c44bbabdcc fix LU and QR solve when rank==0, fix LLT when the matrix is purely 0 2010-07-12 16:42:38 +02:00
Gael Guennebaud
468863f3a0 fix bad Map in row-vector by matrix products 2010-07-09 20:27:05 +02:00
Gael Guennebaud
81f36b0e21 disable MSVC optimization when the underlying compiler is ICC 2010-07-09 19:36:43 +02:00
Benoit Jacob
c196a49f67 Added tag 2.0.14 for changeset e18e51d891 2010-06-22 22:23:39 -04:00
Benoit Jacob
e18e51d891 bump 2010-06-22 22:22:11 -04:00
Stuart Glaser
58fc972ed3 LU on limited-size matrices no longer allocates for temporaries. 2010-06-21 21:29:45 -07:00
Benoit Jacob
e7b6a4bcba fix bug introduced yesterday preventing vectorization of vectors when the storage order is not "the right one".
expand a little the vectorization_logic test and backport EIGEN_DEBUG_ASSIGN.
2010-06-22 09:24:07 -04:00
Benoit Jacob
eaa81c135a fix brain dead computation of the aligned bit.
When using a max-size that is fixed and not a multiple of 16 bit, we're not aligned.
2010-06-21 21:07:53 -04:00
Benoit Jacob
ad8b6c2342 fix #127: remove static keywords that had no effect anyway since the forward-declaration wasn't static,
and that would have been bad if they had taken effect.
2010-06-16 08:28:34 -04:00
Benoit Jacob
37fe67372b Added tag 2.0.13 for changeset ee499a855c 2010-06-10 08:06:38 -04:00
Benoit Jacob
ee499a855c bump version 2010-06-10 08:06:06 -04:00
Benoit Jacob
bc0ce37395 mention that that issue has been solved in GCC 4.5 2010-06-10 08:02:20 -04:00
Benoit Jacob
65c0b2a04d install the Eigen and Dense public headers !!!!! 2010-06-10 08:02:02 -04:00
Benoit Jacob
93d8d0e1b5 add unit test for bug #132 2010-06-09 20:50:18 -04:00
Benoit Jacob
501063d9e9 Fix bug #132
In the matrix-vector products, we were calling coeffRef on the vector xpr without checking it has DirectAccess.

Will add unit test (since it's in 2.0, just import the test case provided in the bug report).

Confirming that this can't happen in the devel branch, and that if we tried to call coeffRef on an xpr
without DirectAccess, that would not compile (since the DenseCoeffsBase class was introduced).
2010-06-09 19:39:05 -04:00
Gael Guennebaud
c76c8d6917 fix issue #125 2010-05-31 22:53:02 +02:00
Thomas Capricelli
0bb8688d70 sync .hgignore with tip 2010-05-26 02:37:52 +02:00
Thomas Capricelli
ea87337647 fix readcost for complex types (backport from tip) 2010-05-26 02:34:55 +02:00
Piotr Trojanek
12557fb2a2 more std:: namespace fixups for 2.0 branch 2010-05-01 09:25:36 +09:00
Benoit Jacob
8a95876825 with QCC, don't try passing --version 2010-04-29 07:58:47 -04:00
Benoit Jacob
82d7c4e1d0 forgot hg add 2010-04-26 07:54:19 -04:00
Benoit Jacob
8d4b0aae04 Fix bug #93: add platform check for how to link to the standard math library.
This allows to support QNX.
2010-04-19 11:30:46 -04:00
Gael Guennebaud
4785e27d6a fix btl compilation 2010-04-01 12:34:55 +02:00
Benoit Jacob
20b544b444 disable all alignment with QCC/QNX in eigen 2.0 2010-03-06 00:31:55 -05:00
Piotr Trojanek
135013c608 EIGEN_ALIGN std:: fixup 2010-03-04 16:38:37 +01:00
Piotr Trojanek
1625a5e3f8 fixups for clean QCC compilation (add std:: in front of size_t, memcpy, etc.) 2010-03-05 11:21:11 +01:00
Piotr Trojanek
47a61bbd80 posix_memalign check for QNX 2010-03-04 17:12:30 +01:00
Benoit Jacob
12bcfae0c5 clarify EIGEN_DONT_ALIGN docs wrt versions 2010-02-28 15:56:50 -05:00
Benoit Jacob
5f42104e0a really fix the LDLt, at the expense of letting isPositiveDefinite() always return true (it was fundamentally broken anyway, especially as in 2.0 we don't even pivot at all).
also, fix compilation
2010-02-23 18:10:31 -05:00
Benoit Jacob
6b3f81b414 precision improvement. Wtf were we using sqrt(precision) for the cutoff? Replaced by precision*biggest. 2010-02-23 16:26:39 -05:00
Gael Guennebaud
1f4b8e6a36 compilation fix in ldlt() for non matrix types
(transplanted from 608959aa6f
)
2010-02-21 10:29:19 +01:00
Benoit Jacob
e1f61b40c8 oops, this had to be done here, not at the end. 2010-02-12 09:02:25 -05:00
Benoit Jacob
f369dc873e Piotr's patch was missing many occurences of size_t. So,
using std::size_t;
This is the only way that we can ensure QCC support in the long term without having to think about it everytime.
2010-02-12 08:57:04 -05:00
Benoit Jacob
c03bca21c4 Added tag 2.0.12 for changeset ed6eb5a625 2010-02-11 21:39:46 -05:00
Benoit Jacob
ed6eb5a625 bump 2010-02-11 21:39:41 -05:00
Benoit Jacob
9488a12125 work around brain dead ICC 2010-02-11 19:32:56 -05:00
Piotr Trojanek
7b44957c4b std:: namespace fixup for more restricive compilers such as QNX's QCC 2010-02-10 22:27:35 +01:00
Hauke Heibel
743ad75595 BenchTimer backport (clock_gettime & QueryPerformanceCounter). 2010-02-03 21:55:01 +01:00
Benoit Jacob
a9eabed421 Patch by 'Wolf' from the issue tracker:
Fix bug #90, missing type cast in LU, allow to use LU with MPFR.
2010-02-02 07:06:15 -05:00
Benoit Jacob
cd34a1d351 backport bug fix by Jitse. 2010-01-28 14:00:09 -05:00
Benoit Jacob
3e963ee69d EIGEN_ENUM_MIN ---> EIGEN_SIZE_MIN 2010-01-26 20:37:57 -05:00
Benoit Jacob
6cc9dc17f2 In LU / Inverse, decouple the output type from the input type.
This has long been done in the default branch
2010-01-26 18:45:23 -05:00
Gael Guennebaud
7852a48a2f fix matrix product with EIGEN_DEFAULT_TO_ROW_MAJOR 2010-01-25 21:56:01 +01:00
Benoit Jacob
d209120180 * Introduce EIGEN_DEFAULT_TO_ROW_MAJOR tests option
---> Now only product_large fails with EIGEN_DEFAULT_TO_ROW_MAJOR.
* Fix EIGEN_NO_ASSERTION_CHECKING tests option
* Fix a crash in Tridiagonalization on row-major matrices + SSE
* Fix inverse test (numeric stability noise)
* Extend map test (see previous fixes in MapBase)
* Fix vectorization_logic test for row-major
* Disable sparse tests with EIGEN_DEFAULT_TO_ROW_MAJOR
2010-01-25 14:00:02 -05:00
Thomas Capricelli
55c0707b1d fix the script again (definitely?) + cleaning 2010-01-22 19:28:33 +01:00
Benoit Jacob
72044ca925 fix a super nasty bug: on row-major expressions that are NOT vectors but that
do have LinearAccess, the MapBase::coeff(int) and MapBase::coeffRef(int)
methods were broken.
2010-01-21 23:33:20 -05:00
Benoit Jacob
c2b8ca7493 if EIGEN_DONT_ALIGN then don't try to vectorize (was giving a #error later on). 2010-01-21 22:32:16 -05:00
Gael Guennebaud
018cb8975a fix plugin doc 2010-01-17 19:55:08 +01:00
Benoit Jacob
3ab280ce4e add missing semicolon in the example 2010-01-17 12:40:19 -05:00
Benoit Jacob
b40030753b Added tag 2.0.11 for changeset 5f73a8df20 2010-01-10 11:30:40 -05:00
Benoit Jacob
5f73a8df20 bump 2010-01-10 11:30:10 -05:00
Thomas Capricelli
8a6d5f10dc backport from tip : actually stop on compile failure 2010-01-06 17:17:40 +01:00
Benoit Jacob
ba6ed5fa5f Fix CoeffReadCost in Part: it must account for the cost of the
conditional jump. This makes Part considered an "expensive" xpr
 that must be evaluated in operations such as Product.
This fixes bug #80.
2010-01-02 13:04:04 -05:00
Benoit Jacob
e4c88c14ec clarify docs as requested on the forum 2010-01-02 12:54:55 -05:00
Benoit Jacob
74207a31fa backport the fix to bug #79, and the unit test 2010-01-02 12:45:49 -05:00
Benoit Jacob
6fd9248c09 add Intel copyright info 2009-12-15 08:43:31 -05:00
Benoit Jacob
4262117f84 backport 4x4 inverse changes:
- use cofactors
 - use Intel's SSE code in the float case
2009-12-15 08:16:48 -05:00
Gael Guennebaud
b581cb870c fix #74: sparse triangular solver for lower/row-major matrices 2009-12-14 10:20:35 +01:00
Gael Guennebaud
72fc81dd9d backport quaternion slerp precision fix 2009-12-05 18:28:17 +01:00
Gael Guennebaud
f36650b00a fix MSVC10 compilation issues 2009-12-02 19:34:37 +01:00
Benoit Jacob
8d31f58ea1 fix bug #70
Was trying to apply stupid invertibility check to top-left 2x2 corner.
2009-11-26 15:33:07 -05:00
Benoit Jacob
a161a70696 Added tag 2.0.10 for changeset 8f1ce52e76 2009-11-25 08:54:17 -05:00
Benoit Jacob
8f1ce52e76 bump 2009-11-25 08:46:42 -05:00
Benoit Jacob
268df314f1 make the inverse_4x4 test pass more consistently 2009-11-25 08:43:20 -05:00
Benoit Jacob
522022ebfc wow, restore Gael's changeset 5455d6fbe8
that I had accidentally undone in my changeset c64ca6870e
.
2009-11-25 08:31:25 -05:00
Thomas Capricelli
d048d7e712 fix bitbucket url after recent change 2009-11-24 23:08:16 +01:00
Benoit Jacob
cd3c8a9404 typo 2009-11-23 11:23:30 -05:00
Benoit Jacob
ec8f37ac75 improve precision test 2009-11-23 11:20:48 -05:00
Benoit Jacob
fc7f39980c backport improvement in 4x4 inverse precisio, and rigorous precision test 2009-11-23 10:27:10 -05:00
Gael Guennebaud
70af59c455 an attempt to fix compilation with recent MSVC 2009-11-23 10:29:40 +01:00
Benoit Jacob
f4dd399499 fix warnings 2009-11-16 14:15:47 -05:00
Benoit Jacob
153557527e backport: init-by-zero option: resize with same size must be a NOP 2009-11-16 13:47:02 -05:00
Benoit Jacob
6aad7f80ff avoid infinite loop, optimization not important, so a plain for loop is the safe way 2009-11-12 14:09:53 -05:00
Benoit Jacob
e3f6c3115a backport the initialize-by-0 option 2009-11-12 12:53:24 -05:00
Benoit Jacob
a2c838ff8f fix PowerPC platform detection 2009-11-11 10:52:00 -05:00
Thomas Capricelli
1e2f56c35a backport of b53c2fcc99
: fix install dir for *.pc

Ingmar Vanhassel <ingmar@exherbo.org>:
Packages that don't install architecture-specific files should install
their pkg-config file to datadir, not libdir.
2009-11-11 15:35:12 +01:00
Hauke Heibel
808c4e9581 Fixed the packport of 62 - Packet4f/d/i does not exist in 2.0. 2009-11-05 10:49:49 +01:00
Hauke Heibel
65331c3884 backporting3979f6d8aad001174160774b49b747430a7686b5
: fixed bug #62
2009-11-04 17:49:34 +01:00
Benoit Jacob
e158cdd61d fix Matrix::Map/MapAligned documentation, and rephrase the tutorial on Map 2009-10-31 14:45:50 -04:00
Benoit Jacob
c64ca6870e this explicit keyword can't hurt 2009-10-31 11:49:20 -04:00
Benoit Jacob
6a90f6c5f0 * default MatrixBase ctor: make it protected, make it a static assert, only do the check when debugging eigen to avoid slowing down compilation for everybody (this check is paranoiac, it's very seldom useful)
* add private MatrixBase ctors to catch cases when the user tries to construct MatrixBase objects directly
2009-10-31 11:48:33 -04:00
Gael Guennebaud
22dd13fdb9 backporting fix of #65 2009-10-29 14:26:38 +01:00
Gael Guennebaud
5455d6fbe8 backporting fix of #65 2009-10-29 14:26:00 +01:00
Benoit Jacob
de693cf34a remove extra ; 2009-10-28 10:04:13 -04:00
Benoit Jacob
21c4e0802d fix potential warning 2009-10-28 09:45:09 -04:00
Benoit Jacob
241b9d34a7 Hey, I was insomniac too ;)
This restores much of the performance benefit of Euler's method, without compromising accuracy (tested on 1e+7 matrices). Namely, my benchmark now runs in 1.5 s instead of 2.2 s. The same in the default branch runs in 1.08 s instead of 1.9 s, so the default branch benefits even more!
2009-10-28 03:50:29 -04:00
Benoit Jacob
9e15a6da2e Fix 4x4 matrix inversion. Applying Euler's trick is more tricky than what "high performance matrix inversion" websites would have you believe. Our 4x4 matrix inversion wasn't numerically stable, because in applying the Euler trick we didn't take the 2x2 block of biggest determinant. As a result, with float, we got relative errors above 1% every 1000 random matrices, and we got completely wrong results every 10000 matrices.
Note that this decreases the performance, but we're still significantly faster than the brutal cofactors approach.
2009-10-27 07:35:25 -04:00
Benoit Jacob
3d365a75cd Added tag 2.0.9 for changeset 38bc82a6f7 2009-10-24 19:38:58 -04:00
Benoit Jacob
38bc82a6f7 bump 2009-10-24 16:35:46 -04:00
Benoit Jacob
6173eb67ff really fix pkgconfig support and make install.
* mistake: was using the install dir instead of binary dir
* was also using INCLUDE_INSTALL_DIR before it was set, so on a first cmake run, the pkgconfig file was bad
2009-10-24 16:16:48 -04:00
Benoit Jacob
9f89431cea install NewStdVector 2009-10-23 19:58:37 -04:00
Benoit Jacob
79e392472a Added tag 2.0.8 for changeset e1c96f3fe0 2009-10-23 18:47:33 -04:00
Benoit Jacob
e1c96f3fe0 bump 2009-10-23 18:37:05 -04:00
Benoit Jacob
46f0fe3b4b SVD:
* implement default ctor, users were relying on the compiler generater one and reported issues on the forum
* turns out that we crash on 1x1 matrices but work on Nx1. No interest in fixing that, so just guard by assert and unit test.
2009-10-23 18:05:55 -04:00
Benoit Jacob
e17e4f3654 just this is enough 2009-10-23 17:51:32 -04:00
Benoit Jacob
2b006ae430 fix "make install"
aaargh
my shiny birthday release, it's broken already!
2009-10-23 17:47:12 -04:00
Benoit Jacob
df6923fd2b Added tag 2.0.7 for changeset 21d081c6da 2009-10-22 16:15:15 -04:00
72 changed files with 1119 additions and 411 deletions

View File

@@ -5,6 +5,7 @@ qrc_*cxx
*.diff
diff
*.save
save
*.old
*.gmo
*.qm
@@ -12,7 +13,7 @@ core
core.*
*.bak
*~
build
build*
*.moc.*
*.moc
ui_*

View File

@@ -1,7 +1,13 @@
project(Eigen)
cmake_minimum_required(VERSION 2.6.2)
set(EIGEN_VERSION_NUMBER "2.0.7")
set(INCLUDE_INSTALL_DIR
"${CMAKE_INSTALL_PREFIX}/include/eigen2"
CACHE PATH
"The directory where we install the header files"
FORCE)
set(EIGEN_VERSION_NUMBER "2.0.15")
set(EIGEN_VERSION "${EIGEN_VERSION_NUMBER}")
set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)
@@ -20,6 +26,38 @@ if(EIGEN_BUILD_LIB)
option(EIGEN_TEST_LIB "Build the unit tests using the library (disable -pedantic)" OFF)
endif(EIGEN_BUILD_LIB)
#############################################################################
# find how to link to the standard libraries #
#############################################################################
find_package(StandardMathLibrary)
set(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO "")
if(NOT STANDARD_MATH_LIBRARY_FOUND)
message(FATAL_ERROR
"Can't link to the standard math library. Please report to the Eigen developers, telling them about your platform.")
else()
if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO)
set(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO "${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO} ${STANDARD_MATH_LIBRARY}")
else()
set(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO "${STANDARD_MATH_LIBRARY}")
endif()
endif()
if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO)
message(STATUS "Standard libraries to link to explicitly: ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO}")
else()
message(STATUS "Standard libraries to link to explicitly: none")
endif()
set(CMAKE_INCLUDE_CURRENT_DIR ON)
if(CMAKE_COMPILER_IS_GNUCXX)
@@ -80,9 +118,9 @@ endif(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION)
include_directories(${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR})
if(EIGEN_BUILD_PKGCONFIG)
configure_file(eigen2.pc.in eigen2.pc)
install(FILES eigen2.pc
DESTINATION lib/pkgconfig
configure_file(eigen2.pc.in eigen2.pc) # uses INCLUDE_INSTALL_DIR
install(FILES ${CMAKE_CURRENT_BINARY_DIR}/eigen2.pc
DESTINATION share/pkgconfig
)
endif(EIGEN_BUILD_PKGCONFIG)

View File

@@ -1,4 +1,7 @@
set(Eigen_HEADERS Core LU Cholesky QR Geometry Sparse Array SVD LeastSquares QtAlignedMalloc StdVector)
set(Eigen_HEADERS Core LU Cholesky QR Geometry
Sparse Array SVD LeastSquares
QtAlignedMalloc StdVector NewStdVector
Eigen Dense)
if(EIGEN_BUILD_LIB)
set(Eigen_SRCS
@@ -20,12 +23,6 @@ if(CMAKE_COMPILER_IS_GNUCXX)
set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "${CMAKE_CXX_FLAGS_RELWITHDEBINFO} -g1 -O2")
endif(CMAKE_COMPILER_IS_GNUCXX)
set(INCLUDE_INSTALL_DIR
"${CMAKE_INSTALL_PREFIX}/include/eigen2"
CACHE PATH
"The directory where we install the header files"
FORCE)
install(FILES
${Eigen_HEADERS}
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen

View File

@@ -15,7 +15,9 @@
#endif
#endif
#ifdef __GNUC__
// FIXME: this check should not be against __QNXNTO__, which is also defined
// while compiling with GCC for QNX target. Better solution is welcome!
#if defined(__GNUC__) && !defined(__QNXNTO__)
#define EIGEN_GNUC_AT_LEAST(x,y) ((__GNUC__>=x && __GNUC_MINOR__>=y) || __GNUC__>x)
#else
#define EIGEN_GNUC_AT_LEAST(x,y) 0
@@ -26,7 +28,7 @@
#define EIGEN_SSE2_BUT_NOT_OLD_GCC
#endif
#ifndef EIGEN_DONT_VECTORIZE
#if !defined(EIGEN_DONT_VECTORIZE) && !defined(EIGEN_DONT_ALIGN)
#if defined (EIGEN_SSE2_BUT_NOT_OLD_GCC) || defined(EIGEN_SSE2_ON_MSVC_2008_OR_LATER)
#define EIGEN_VECTORIZE
#define EIGEN_VECTORIZE_SSE
@@ -77,6 +79,10 @@
namespace Eigen {
// we use size_t frequently and we'll never remember to prepend it with std:: everytime just to
// ensure QNX/QCC support
using std::size_t;
/** \defgroup Core_Module Core module
* This is the main module of Eigen providing dense matrix and vector support
* (both fixed and dynamic size) with all the features corresponding to a BLAS library

View File

@@ -36,8 +36,8 @@ template <class T>
class aligned_allocator_indirection : public aligned_allocator<T>
{
public:
typedef size_t size_type;
typedef ptrdiff_t difference_type;
typedef std::size_t size_type;
typedef std::ptrdiff_t difference_type;
typedef T* pointer;
typedef const T* const_pointer;
typedef T& reference;

View File

@@ -139,7 +139,7 @@ inline bool MatrixBase<Derived>::any() const
template<typename Derived>
inline int MatrixBase<Derived>::count() const
{
return this->cast<bool>().cast<int>().sum();
return this->cast<bool>().template cast<int>().sum();
}
#endif // EIGEN_ALLANDANY_H

View File

@@ -160,13 +160,15 @@ template<typename ExpressionType, int Direction> class PartialRedux
public:
typedef typename ei_traits<ExpressionType>::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef typename ei_meta_if<ei_must_nest_by_value<ExpressionType>::ret,
ExpressionType, const ExpressionType&>::ret ExpressionTypeNested;
template<template<typename _Scalar> class Functor> struct ReturnType
template<template<typename _Scalar> class Functor,
typename Scalar = typename ei_traits<ExpressionType>::Scalar> struct ReturnType
{
typedef PartialReduxExpr<ExpressionType,
Functor<typename ei_traits<ExpressionType>::Scalar>,
Functor<Scalar>,
Direction
> Type;
};
@@ -217,7 +219,7 @@ template<typename ExpressionType, int Direction> class PartialRedux
* Output: \verbinclude PartialRedux_squaredNorm.out
*
* \sa MatrixBase::squaredNorm() */
const typename ReturnType<ei_member_squaredNorm>::Type squaredNorm() const
const typename ReturnType<ei_member_squaredNorm,RealScalar>::Type squaredNorm() const
{ return _expression(); }
/** \returns a row (or column) vector expression of the norm
@@ -227,7 +229,7 @@ template<typename ExpressionType, int Direction> class PartialRedux
* Output: \verbinclude PartialRedux_norm.out
*
* \sa MatrixBase::norm() */
const typename ReturnType<ei_member_norm>::Type norm() const
const typename ReturnType<ei_member_norm,RealScalar>::Type norm() const
{ return _expression(); }
/** \returns a row (or column) vector expression of the sum

View File

@@ -96,8 +96,7 @@ void LDLT<MatrixType>::compute(const MatrixType& a)
assert(a.rows()==a.cols());
const int size = a.rows();
m_matrix.resize(size, size);
m_isPositiveDefinite = true;
const RealScalar eps = ei_sqrt(precision<Scalar>());
m_isPositiveDefinite = true; // always true. This decomposition is not rank-revealing anyway.
if (size<=1)
{
@@ -121,12 +120,6 @@ void LDLT<MatrixType>::compute(const MatrixType& a)
RealScalar tmp = ei_real(a.coeff(j,j) - (m_matrix.row(j).start(j) * m_matrix.col(j).start(j).conjugate()).coeff(0,0));
m_matrix.coeffRef(j,j) = tmp;
if (tmp < eps)
{
m_isPositiveDefinite = false;
return;
}
int endSize = size-j-1;
if (endSize>0)
{
@@ -136,7 +129,8 @@ void LDLT<MatrixType>::compute(const MatrixType& a)
m_matrix.row(j).end(endSize) = a.row(j).end(endSize).conjugate()
- _temporary.end(endSize).transpose();
m_matrix.col(j).end(endSize) = m_matrix.row(j).end(endSize) / tmp;
if(tmp != RealScalar(0))
m_matrix.col(j).end(endSize) = m_matrix.row(j).end(endSize) / tmp;
}
}
}
@@ -192,7 +186,7 @@ template<typename Derived>
inline const LDLT<typename MatrixBase<Derived>::PlainMatrixType>
MatrixBase<Derived>::ldlt() const
{
return derived();
return LDLT<PlainMatrixType>(derived());
}
#endif // EIGEN_LDLT_H

View File

@@ -132,11 +132,12 @@ void LLT<MatrixType>::compute(const MatrixType& a)
m_isInitialized = true;
return;
}
m_matrix.col(0).end(size-1) = a.row(0).end(size-1).adjoint() / ei_real(m_matrix.coeff(0,0));
if(ei_real(m_matrix.coeff(0,0))>0)
m_matrix.col(0).end(size-1) = a.row(0).end(size-1).adjoint() / ei_real(m_matrix.coeff(0,0));
for (int j = 1; j < size; ++j)
{
x = ei_real(a.coeff(j,j)) - m_matrix.row(j).start(j).squaredNorm();
if (x < cutoff)
if (x <= cutoff)
{
m_isPositiveDefinite = false;
continue;

View File

@@ -90,6 +90,28 @@ public:
? ( int(MayUnrollCompletely) && int(DstIsAligned) ? int(CompleteUnrolling) : int(NoUnrolling) )
: int(NoUnrolling)
};
#ifdef EIGEN_DEBUG_ASSIGN
#define EIGEN_DEBUG_VAR(x) std::cerr << #x << " = " << x << std::endl;
static void debug()
{
EIGEN_DEBUG_VAR(DstIsAligned)
EIGEN_DEBUG_VAR(SrcIsAligned)
EIGEN_DEBUG_VAR(SrcAlignment)
EIGEN_DEBUG_VAR(InnerSize)
EIGEN_DEBUG_VAR(InnerMaxSize)
EIGEN_DEBUG_VAR(PacketSize)
EIGEN_DEBUG_VAR(MightVectorize)
EIGEN_DEBUG_VAR(MayInnerVectorize)
EIGEN_DEBUG_VAR(MayLinearVectorize)
EIGEN_DEBUG_VAR(MaySliceVectorize)
EIGEN_DEBUG_VAR(Vectorization)
EIGEN_DEBUG_VAR(UnrollingLimit)
EIGEN_DEBUG_VAR(MayUnrollCompletely)
EIGEN_DEBUG_VAR(MayUnrollInner)
EIGEN_DEBUG_VAR(Unrolling)
}
#endif
};
/***************************************************************************
@@ -400,6 +422,9 @@ template<typename OtherDerived>
EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>
::lazyAssign(const MatrixBase<OtherDerived>& other)
{
#ifdef EIGEN_DEBUG_ASSIGN
ei_assign_traits<Derived, OtherDerived>::debug();
#endif
EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Derived,OtherDerived)
EIGEN_STATIC_ASSERT((ei_is_same_type<typename Derived::Scalar, typename OtherDerived::Scalar>::ret),
YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)

View File

@@ -222,7 +222,7 @@ class Block<MatrixType,BlockRows,BlockCols,PacketAccess,HasDirectAccess>
class InnerIterator;
typedef typename ei_traits<Block>::AlignedDerivedType AlignedDerivedType;
friend class Block<MatrixType,BlockRows,BlockCols,PacketAccess==AsRequested?ForceAligned:AsRequested,HasDirectAccess>;
friend class Block<MatrixType,BlockRows,BlockCols,PacketAccess==int(AsRequested)?ForceAligned:AsRequested,HasDirectAccess>;
EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Block)

View File

@@ -33,7 +33,7 @@ struct ei_L2_block_traits {
#ifndef EIGEN_EXTERN_INSTANTIATIONS
template<typename Scalar>
static void ei_cache_friendly_product(
void ei_cache_friendly_product(
int _rows, int _cols, int depth,
bool _lhsRowMajor, const Scalar* _lhs, int _lhsStride,
bool _rhsRowMajor, const Scalar* _rhs, int _rhsStride,
@@ -84,7 +84,7 @@ static void ei_cache_friendly_product(
MaxL2BlockSize = ei_L2_block_traits<EIGEN_TUNE_FOR_CPU_CACHE_SIZE,Scalar>::width
};
const bool resIsAligned = (PacketSize==1) || (((resStride%PacketSize) == 0) && (size_t(res)%16==0));
const bool resIsAligned = (PacketSize==1) || (((resStride%PacketSize) == 0) && (std::size_t(res)%16==0));
const int remainingSize = depth % PacketSize;
const int size = depth - remainingSize; // third dimension of the product clamped to packet boundaries
@@ -92,7 +92,7 @@ static void ei_cache_friendly_product(
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));
const bool needRhsCopy = (PacketSize>1) && ((rhsStride%PacketSize!=0) || (std::size_t(rhs)%16!=0));
Scalar* EIGEN_RESTRICT block = 0;
const int allocBlockSize = l2BlockRows*size;
block = ei_aligned_stack_new(Scalar, allocBlockSize);
@@ -172,7 +172,7 @@ static void ei_cache_friendly_product(
for(int l1j=l2j; l1j<l2blockColEnd; l1j+=1)
{
ei_internal_assert(l2BlockSizeAligned*(l1j-l2j)+(l2blockSizeEnd-l2k) < l2BlockSizeAligned*l2BlockSizeAligned);
memcpy(rhsCopy+l2BlockSizeAligned*(l1j-l2j),&(rhs[l1j*rhsStride+l2k]),(l2blockSizeEnd-l2k)*sizeof(Scalar));
std::memcpy(rhsCopy+l2BlockSizeAligned*(l1j-l2j),&(rhs[l1j*rhsStride+l2k]),(l2blockSizeEnd-l2k)*sizeof(Scalar));
}
// for each bw x 1 result's block
@@ -352,7 +352,7 @@ static void ei_cache_friendly_product(
* TODO: since rhs gets evaluated only once, no need to evaluate it
*/
template<typename Scalar, typename RhsType>
static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
int size,
const Scalar* lhs, int lhsStride,
const RhsType& rhs,
@@ -397,7 +397,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
int skipColumns = 0;
if (PacketSize>1)
{
ei_internal_assert(size_t(lhs+lhsAlignmentOffset)%sizeof(Packet)==0 || size<PacketSize);
ei_internal_assert(std::size_t(lhs+lhsAlignmentOffset)%sizeof(Packet)==0 || size<PacketSize);
while (skipColumns<PacketSize &&
alignedStart != ((lhsAlignmentOffset + alignmentStep*skipColumns)%PacketSize))
@@ -414,7 +414,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
// note that the skiped columns are processed later.
}
ei_internal_assert((alignmentPattern==NoneAligned) || (size_t(lhs+alignedStart+lhsStride*skipColumns)%sizeof(Packet))==0);
ei_internal_assert((alignmentPattern==NoneAligned) || (std::size_t(lhs+alignedStart+lhsStride*skipColumns)%sizeof(Packet))==0);
}
int offset1 = (FirstAligned && alignmentStep==1?3:1);
@@ -516,7 +516,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
res[j] += ei_pfirst(ptmp0) * lhs0[j];
// process aligned result's coeffs
if ((size_t(lhs0+alignedStart)%sizeof(Packet))==0)
if ((std::size_t(lhs0+alignedStart)%sizeof(Packet))==0)
for (int j = alignedStart;j<alignedSize;j+=PacketSize)
ei_pstore(&res[j], ei_pmadd(ptmp0,ei_pload(&lhs0[j]),ei_pload(&res[j])));
else
@@ -542,7 +542,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
// TODO add peeling to mask unaligned load/stores
template<typename Scalar, typename ResType>
static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
const Scalar* lhs, int lhsStride,
const Scalar* rhs, int rhsSize,
ResType& res)
@@ -586,7 +586,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
int skipRows = 0;
if (PacketSize>1)
{
ei_internal_assert(size_t(lhs+lhsAlignmentOffset)%sizeof(Packet)==0 || size<PacketSize);
ei_internal_assert(std::size_t(lhs+lhsAlignmentOffset)%sizeof(Packet)==0 || size<PacketSize);
while (skipRows<PacketSize &&
alignedStart != ((lhsAlignmentOffset + alignmentStep*skipRows)%PacketSize))
@@ -603,7 +603,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
// note that the skiped columns are processed later.
}
ei_internal_assert((alignmentPattern==NoneAligned) || PacketSize==1
|| (size_t(lhs+alignedStart+lhsStride*skipRows)%sizeof(Packet))==0);
|| (std::size_t(lhs+alignedStart+lhsStride*skipRows)%sizeof(Packet))==0);
}
int offset1 = (FirstAligned && alignmentStep==1?3:1);
@@ -722,7 +722,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
if (alignedSize>alignedStart)
{
// process aligned rhs coeffs
if ((size_t(lhs0+alignedStart)%sizeof(Packet))==0)
if ((std::size_t(lhs0+alignedStart)%sizeof(Packet))==0)
for (int j = alignedStart;j<alignedSize;j+=PacketSize)
ptmp0 = ei_pmadd(ei_pload(&rhs[j]), ei_pload(&lhs0[j]), ptmp0);
else

View File

@@ -32,7 +32,7 @@ namespace Eigen
{
#define EIGEN_INSTANTIATE_PRODUCT(TYPE) \
template static void ei_cache_friendly_product<TYPE>( \
template void ei_cache_friendly_product<TYPE>( \
int _rows, int _cols, int depth, \
bool _lhsRowMajor, const TYPE* _lhs, int _lhsStride, \
bool _rhsRowMajor, const TYPE* _rhs, int _rhsStride, \

View File

@@ -47,11 +47,11 @@ struct ei_traits<DiagonalCoeffs<MatrixType> >
typedef typename ei_unref<MatrixTypeNested>::type _MatrixTypeNested;
enum {
RowsAtCompileTime = int(MatrixType::SizeAtCompileTime) == Dynamic ? Dynamic
: EIGEN_ENUM_MIN(MatrixType::RowsAtCompileTime,
: EIGEN_SIZE_MIN(MatrixType::RowsAtCompileTime,
MatrixType::ColsAtCompileTime),
ColsAtCompileTime = 1,
MaxRowsAtCompileTime = int(MatrixType::MaxSizeAtCompileTime) == Dynamic ? Dynamic
: EIGEN_ENUM_MIN(MatrixType::MaxRowsAtCompileTime,
: EIGEN_SIZE_MIN(MatrixType::MaxRowsAtCompileTime,
MatrixType::MaxColsAtCompileTime),
MaxColsAtCompileTime = 1,
Flags = (unsigned int)_MatrixTypeNested::Flags & (HereditaryBits | LinearAccessBit),

View File

@@ -99,7 +99,7 @@ template<typename Derived> class MapBase
inline const Scalar coeff(int index) const
{
ei_assert(Derived::IsVectorAtCompileTime || (ei_traits<Derived>::Flags & LinearAccessBit));
if ( ((RowsAtCompileTime == 1) == IsRowMajor) )
if ( ((RowsAtCompileTime == 1) == IsRowMajor) || !int(Derived::IsVectorAtCompileTime) )
return m_data[index];
else
return m_data[index*stride()];
@@ -108,7 +108,7 @@ template<typename Derived> class MapBase
inline Scalar& coeffRef(int index)
{
ei_assert(Derived::IsVectorAtCompileTime || (ei_traits<Derived>::Flags & LinearAccessBit));
if ( ((RowsAtCompileTime == 1) == IsRowMajor) )
if ( ((RowsAtCompileTime == 1) == IsRowMajor) || !int(Derived::IsVectorAtCompileTime) )
return const_cast<Scalar*>(m_data)[index];
else
return const_cast<Scalar*>(m_data)[index*stride()];

View File

@@ -35,16 +35,6 @@ template<typename T> inline T ei_random_amplitude()
else return static_cast<T>(10);
}
template<typename T> inline T ei_hypot(T x, T y)
{
T _x = ei_abs(x);
T _y = ei_abs(y);
T p = std::max(_x, _y);
T q = std::min(_x, _y);
T qp = q/p;
return p * ei_sqrt(T(1) + qp*qp);
}
/**************
*** int ***
**************/
@@ -54,7 +44,7 @@ template<> inline int machine_epsilon<int>() { return 0; }
inline int ei_real(int x) { return x; }
inline int ei_imag(int) { return 0; }
inline int ei_conj(int x) { return x; }
inline int ei_abs(int x) { return abs(x); }
inline int ei_abs(int x) { return std::abs(x); }
inline int ei_abs2(int x) { return x*x; }
inline int ei_sqrt(int) { ei_assert(false); return 0; }
inline int ei_exp(int) { ei_assert(false); return 0; }
@@ -67,7 +57,7 @@ inline int ei_pow(int x, int y) { return int(std::pow(double(x), y)); }
template<> inline int ei_random(int a, int b)
{
// We can't just do rand()%n as only the high-order bits are really random
return a + static_cast<int>((b-a+1) * (rand() / (RAND_MAX + 1.0)));
return a + static_cast<int>((b-a+1) * (std::rand() / (RAND_MAX + 1.0)));
}
template<> inline int ei_random()
{
@@ -292,4 +282,14 @@ inline bool ei_isApproxOrLessThan(long double a, long double b, long double prec
return a <= b || ei_isApprox(a, b, prec);
}
template<typename T> inline T ei_hypot(T x, T y)
{
T _x = ei_abs(x);
T _y = ei_abs(y);
T p = std::max(_x, _y);
T q = std::min(_x, _y);
T qp = q/p;
return p * ei_sqrt(T(1) + qp*qp);
}
#endif // EIGEN_MATHFUNCTIONS_H

View File

@@ -25,6 +25,11 @@
#ifndef EIGEN_MATRIX_H
#define EIGEN_MATRIX_H
#ifdef EIGEN_INITIALIZE_MATRICES_BY_ZERO
# define EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED for(int i=0;i<base().size();++i) coeffRef(i)=Scalar(0);
#else
# define EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED
#endif
/** \class Matrix
*
@@ -233,7 +238,14 @@ class Matrix
&& (RowsAtCompileTime == Dynamic || RowsAtCompileTime == rows)
&& (MaxColsAtCompileTime == Dynamic || MaxColsAtCompileTime >= cols)
&& (ColsAtCompileTime == Dynamic || ColsAtCompileTime == cols));
m_storage.resize(rows * cols, rows, cols);
#ifdef EIGEN_INITIALIZE_MATRICES_BY_ZERO
int size = rows*cols;
bool size_changed = size != this->size();
m_storage.resize(size, rows, cols);
if(size_changed) EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED
#else
m_storage.resize(rows*cols, rows, cols);
#endif
}
/** Resizes \c *this to a vector of length \a size
@@ -243,10 +255,17 @@ class Matrix
inline void resize(int size)
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Matrix)
ei_assert(SizeAtCompileTime == Dynamic || SizeAtCompileTime == size);
#ifdef EIGEN_INITIALIZE_MATRICES_BY_ZERO
bool size_changed = size != this->size();
#endif
if(RowsAtCompileTime == 1)
m_storage.resize(size, 1, size);
else
m_storage.resize(size, size, 1);
#ifdef EIGEN_INITIALIZE_MATRICES_BY_ZERO
if(size_changed) EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED
#endif
}
/** Copies the value of the expression \a other into \c *this with automatic resizing.
@@ -290,13 +309,14 @@ class Matrix
EIGEN_STRONG_INLINE explicit Matrix() : m_storage()
{
_check_template_params();
EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED
}
#ifndef EIGEN_PARSED_BY_DOXYGEN
/** \internal */
Matrix(ei_constructor_without_unaligned_array_assert)
: m_storage(ei_constructor_without_unaligned_array_assert())
{}
{EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED}
#endif
/** Constructs a vector or row-vector with given dimension. \only_for_vectors
@@ -312,6 +332,7 @@ class Matrix
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Matrix)
ei_assert(dim > 0);
ei_assert(SizeAtCompileTime == Dynamic || SizeAtCompileTime == dim);
EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED
}
/** This constructor has two very different behaviors, depending on the type of *this.
@@ -337,6 +358,7 @@ class Matrix
{
ei_assert(x > 0 && (RowsAtCompileTime == Dynamic || RowsAtCompileTime == x)
&& y > 0 && (ColsAtCompileTime == Dynamic || ColsAtCompileTime == y));
EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED
}
}
/** constructs an initialized 2D vector with given coefficients */
@@ -402,9 +424,10 @@ class Matrix
void swap(const MatrixBase<OtherDerived>& other);
/** \name Map
* These are convenience functions returning Map objects. The Map() static functions return unaligned Map objects,
* while the AlignedMap() functions return aligned Map objects and thus should be called only with 16-byte-aligned
* \a data pointers.
* These are convenience functions returning Map objects.
*
* \warning Do not use MapAligned in the Eigen 2.0. Mapping aligned arrays will be fully
* supported in Eigen 3.0 (already implemented in the development branch)
*
* \see class Map
*/
@@ -569,7 +592,8 @@ template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int
template<typename OtherDerived>
inline void Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::swap(const MatrixBase<OtherDerived>& other)
{
ei_matrix_swap_impl<Matrix, OtherDerived>::run(*this, *const_cast<MatrixBase<OtherDerived>*>(&other));
// the Eigen:: here is to work around a stupid ICC 11.1 bug.
Eigen::ei_matrix_swap_impl<Matrix, OtherDerived>::run(*this, *const_cast<MatrixBase<OtherDerived>*>(&other));
}

View File

@@ -136,12 +136,6 @@ template<typename Derived> class MatrixBase
*/
};
/** Default constructor. Just checks at compile-time for self-consistency of the flags. */
MatrixBase()
{
ei_assert(ei_are_flags_consistent<Flags>::ret);
}
#ifndef EIGEN_PARSED_BY_DOXYGEN
/** This is the "real scalar" type; if the \a Scalar type is already real numbers
* (e.g. int, float or double) then \a RealScalar is just the same as \a Scalar. If
@@ -165,7 +159,7 @@ template<typename Derived> class MatrixBase
inline int size() const { return rows() * cols(); }
/** \returns the number of nonzero coefficients which is in practice the number
* of stored coefficients. */
inline int nonZeros() const { return derived.nonZeros(); }
inline int nonZeros() const { return size(); }
/** \returns true if either the number of rows or the number of columns is equal to 1.
* In other words, this function returns
* \code rows()==1 || cols()==1 \endcode
@@ -589,7 +583,8 @@ template<typename Derived> class MatrixBase
const LU<PlainMatrixType> lu() const;
const PlainMatrixType inverse() const;
void computeInverse(PlainMatrixType *result) const;
template<typename ResultType>
void computeInverse(MatrixBase<ResultType> *result) const;
Scalar determinant() const;
/////////// Cholesky module ///////////
@@ -627,6 +622,24 @@ template<typename Derived> class MatrixBase
#ifdef EIGEN_MATRIXBASE_PLUGIN
#include EIGEN_MATRIXBASE_PLUGIN
#endif
protected:
/** Default constructor. Do nothing. */
MatrixBase()
{
/* Just checks for self-consistency of the flags.
* Only do it when debugging Eigen, as this borders on paranoiac and could slow compilation down
*/
#ifdef EIGEN_INTERNAL_DEBUGGING
EIGEN_STATIC_ASSERT(ei_are_flags_consistent<Flags>::ret,
INVALID_MATRIXBASE_TEMPLATE_PARAMETERS)
#endif
}
private:
explicit MatrixBase(int);
MatrixBase(int,int);
template<typename OtherDerived> explicit MatrixBase(const MatrixBase<OtherDerived>&);
};
#endif // EIGEN_MATRIXBASE_H

View File

@@ -94,7 +94,7 @@ template<typename _Real> struct NumTraits<std::complex<_Real> >
enum {
IsComplex = 1,
HasFloatingPoint = NumTraits<Real>::HasFloatingPoint,
ReadCost = 2,
ReadCost = 2 * NumTraits<_Real>::ReadCost,
AddCost = 2 * NumTraits<Real>::AddCost,
MulCost = 4 * NumTraits<Real>::MulCost + 2 * NumTraits<Real>::AddCost
};

View File

@@ -50,7 +50,7 @@ struct ei_traits<Part<MatrixType, Mode> > : ei_traits<MatrixType>
typedef typename ei_unref<MatrixTypeNested>::type _MatrixTypeNested;
enum {
Flags = (_MatrixTypeNested::Flags & (HereditaryBits) & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit))) | Mode,
CoeffReadCost = _MatrixTypeNested::CoeffReadCost
CoeffReadCost = _MatrixTypeNested::CoeffReadCost + ConditionalJumpCost
};
};

View File

@@ -66,11 +66,8 @@ struct ProductReturnType
template<typename Lhs, typename Rhs>
struct ProductReturnType<Lhs,Rhs,CacheFriendlyProduct>
{
typedef typename ei_nested<Lhs,Rhs::ColsAtCompileTime>::type LhsNested;
typedef typename ei_nested<Rhs,Lhs::RowsAtCompileTime,
typename ei_plain_matrix_type_column_major<Rhs>::type
>::type RhsNested;
typedef const Lhs& LhsNested;
typedef const Rhs& RhsNested;
typedef Product<LhsNested, RhsNested, CacheFriendlyProduct> Type;
};
@@ -128,7 +125,7 @@ struct ei_traits<Product<LhsNested, RhsNested, ProductMode> >
RowsAtCompileTime = _LhsNested::RowsAtCompileTime,
ColsAtCompileTime = _RhsNested::ColsAtCompileTime,
InnerSize = EIGEN_ENUM_MIN(_LhsNested::ColsAtCompileTime, _RhsNested::RowsAtCompileTime),
InnerSize = EIGEN_SIZE_MIN(_LhsNested::ColsAtCompileTime, _RhsNested::RowsAtCompileTime),
MaxRowsAtCompileTime = _LhsNested::MaxRowsAtCompileTime,
MaxColsAtCompileTime = _RhsNested::MaxColsAtCompileTime,
@@ -144,7 +141,7 @@ struct ei_traits<Product<LhsNested, RhsNested, ProductMode> >
EvalToRowMajor = RhsRowMajor && (ProductMode==(int)CacheFriendlyProduct ? LhsRowMajor : (!CanVectorizeLhs)),
RemovedBits = ~(EvalToRowMajor ? 0 : RowMajorBit),
RemovedBits = ~((EvalToRowMajor ? 0 : RowMajorBit)|DirectAccessBit),
Flags = ((unsigned int)(LhsFlags | RhsFlags) & HereditaryBits & RemovedBits)
| EvalBeforeAssigningBit
@@ -571,7 +568,7 @@ struct ei_cache_friendly_product_selector<ProductType,LhsRows,ColMajor,HasDirect
else
{
_res = ei_aligned_stack_new(Scalar,res.size());
Map<Matrix<Scalar,DestDerived::RowsAtCompileTime,1> >(_res, res.size()) = res;
Map<Matrix<Scalar,DestDerived::RowsAtCompileTime,1,ColMajor> >(_res, res.size()) = res;
}
ei_cache_friendly_product_colmajor_times_vector(res.size(),
&product.lhs().const_cast_derived().coeffRef(0,0), product.lhs().stride(),
@@ -579,7 +576,7 @@ struct ei_cache_friendly_product_selector<ProductType,LhsRows,ColMajor,HasDirect
if (!EvalToRes)
{
res = Map<Matrix<Scalar,DestDerived::SizeAtCompileTime,1> >(_res, res.size());
res = Map<Matrix<Scalar,DestDerived::SizeAtCompileTime,1,ColMajor> >(_res, res.size());
ei_aligned_stack_delete(Scalar, _res, res.size());
}
}
@@ -617,7 +614,7 @@ struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCo
else
{
_res = ei_aligned_stack_new(Scalar, res.size());
Map<Matrix<Scalar,DestDerived::SizeAtCompileTime,1> >(_res, res.size()) = res;
Map<Matrix<Scalar,1,DestDerived::SizeAtCompileTime,ColMajor> >(_res, res.size()) = res;
}
ei_cache_friendly_product_colmajor_times_vector(res.size(),
&product.rhs().const_cast_derived().coeffRef(0,0), product.rhs().stride(),
@@ -625,7 +622,7 @@ struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCo
if (!EvalToRes)
{
res = Map<Matrix<Scalar,DestDerived::SizeAtCompileTime,1> >(_res, res.size());
res = Map<Matrix<Scalar,1,DestDerived::SizeAtCompileTime,ColMajor> >(_res, res.size());
ei_aligned_stack_delete(Scalar, _res, res.size());
}
}
@@ -639,6 +636,7 @@ struct ei_cache_friendly_product_selector<ProductType,LhsRows,RowMajor,HasDirect
typedef typename ei_traits<ProductType>::_RhsNested Rhs;
enum {
UseRhsDirectly = ((ei_packet_traits<Scalar>::size==1) || (Rhs::Flags&ActualPacketAccessBit))
&& (Rhs::Flags&DirectAccessBit)
&& (!(Rhs::Flags & RowMajorBit)) };
template<typename DestDerived>
@@ -650,7 +648,7 @@ struct ei_cache_friendly_product_selector<ProductType,LhsRows,RowMajor,HasDirect
else
{
_rhs = ei_aligned_stack_new(Scalar, product.rhs().size());
Map<Matrix<Scalar,Rhs::SizeAtCompileTime,1> >(_rhs, product.rhs().size()) = product.rhs();
Map<Matrix<Scalar,Rhs::SizeAtCompileTime,1,ColMajor> >(_rhs, product.rhs().size()) = product.rhs();
}
ei_cache_friendly_product_rowmajor_times_vector(&product.lhs().const_cast_derived().coeffRef(0,0), product.lhs().stride(),
_rhs, product.rhs().size(), res);
@@ -667,6 +665,7 @@ struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCo
typedef typename ei_traits<ProductType>::_LhsNested Lhs;
enum {
UseLhsDirectly = ((ei_packet_traits<Scalar>::size==1) || (Lhs::Flags&ActualPacketAccessBit))
&& (Lhs::Flags&DirectAccessBit)
&& (Lhs::Flags & RowMajorBit) };
template<typename DestDerived>
@@ -678,7 +677,7 @@ struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCo
else
{
_lhs = ei_aligned_stack_new(Scalar, product.lhs().size());
Map<Matrix<Scalar,Lhs::SizeAtCompileTime,1> >(_lhs, product.lhs().size()) = product.lhs();
Map<Matrix<Scalar,1,Lhs::SizeAtCompileTime,ColMajor> >(_lhs, product.lhs().size()) = product.lhs();
}
ei_cache_friendly_product_rowmajor_times_vector(&product.rhs().const_cast_derived().coeffRef(0,0), product.rhs().stride(),
_lhs, product.lhs().size(), res);
@@ -709,7 +708,17 @@ MatrixBase<Derived>::operator+=(const Flagged<Product<Lhs,Rhs,CacheFriendlyProdu
if (other._expression()._useCacheFriendlyProduct())
ei_cache_friendly_product_selector<Product<Lhs,Rhs,CacheFriendlyProduct> >::run(const_cast_derived(), other._expression());
else
lazyAssign(derived() + other._expression());
{
typedef typename ei_cleantype<Lhs>::type _Lhs;
typedef typename ei_cleantype<Rhs>::type _Rhs;
typedef typename ei_nested<_Lhs,_Rhs::ColsAtCompileTime>::type LhsNested;
typedef typename ei_nested<_Rhs,_Lhs::RowsAtCompileTime>::type RhsNested;
Product<LhsNested,RhsNested,NormalProduct> prod(other._expression().lhs(),other._expression().rhs());
lazyAssign(derived() + prod);
}
return derived();
}
@@ -724,12 +733,21 @@ inline Derived& MatrixBase<Derived>::lazyAssign(const Product<Lhs,Rhs,CacheFrien
}
else
{
lazyAssign<Product<Lhs,Rhs,CacheFriendlyProduct> >(product);
typedef typename ei_cleantype<Lhs>::type _Lhs;
typedef typename ei_cleantype<Rhs>::type _Rhs;
typedef typename ei_nested<_Lhs,_Rhs::ColsAtCompileTime>::type LhsNested;
typedef typename ei_nested<_Rhs,_Lhs::RowsAtCompileTime>::type RhsNested;
typedef Product<LhsNested,RhsNested,NormalProduct> NormalProduct;
NormalProduct normal_prod(product.lhs(),product.rhs());
lazyAssign<NormalProduct>(normal_prod);
}
return derived();
}
template<typename T> struct ei_product_copy_rhs
template<typename T,int StorageOrder> struct ei_product_copy_rhs
{
typedef typename ei_meta_if<
(ei_traits<T>::Flags & RowMajorBit)
@@ -739,11 +757,30 @@ template<typename T> struct ei_product_copy_rhs
>::ret type;
};
template<typename T> struct ei_product_copy_lhs
template<typename T> struct ei_product_copy_rhs<T,RowMajorBit>
{
typedef typename ei_meta_if<
(!(ei_traits<T>::Flags & DirectAccessBit)),
typename ei_plain_matrix_type<T>::type,
const T&
>::ret type;
};
template<typename T,int StorageOrder> struct ei_product_copy_lhs
{
typedef typename ei_meta_if<
(!(int(ei_traits<T>::Flags) & DirectAccessBit)),
typename ei_plain_matrix_type<T>::type,
typename ei_plain_matrix_type_row_major<T>::type,
const T&
>::ret type;
};
template<typename T> struct ei_product_copy_lhs<T,RowMajorBit>
{
typedef typename ei_meta_if<
((ei_traits<T>::Flags & RowMajorBit)==0)
|| (!(int(ei_traits<T>::Flags) & DirectAccessBit)),
typename ei_plain_matrix_type_row_major<T>::type,
const T&
>::ret type;
};
@@ -752,9 +789,9 @@ template<typename Lhs, typename Rhs, int ProductMode>
template<typename DestDerived>
inline void Product<Lhs,Rhs,ProductMode>::_cacheFriendlyEvalAndAdd(DestDerived& res) const
{
typedef typename ei_product_copy_lhs<_LhsNested>::type LhsCopy;
typedef typename ei_product_copy_lhs<_LhsNested,DestDerived::Flags&RowMajorBit>::type LhsCopy;
typedef typename ei_unref<LhsCopy>::type _LhsCopy;
typedef typename ei_product_copy_rhs<_RhsNested>::type RhsCopy;
typedef typename ei_product_copy_rhs<_RhsNested,DestDerived::Flags&RowMajorBit>::type RhsCopy;
typedef typename ei_unref<RhsCopy>::type _RhsCopy;
LhsCopy lhs(m_lhs);
RhsCopy rhs(m_rhs);
@@ -764,6 +801,7 @@ inline void Product<Lhs,Rhs,ProductMode>::_cacheFriendlyEvalAndAdd(DestDerived&
_RhsCopy::Flags&RowMajorBit, (const Scalar*)&(rhs.const_cast_derived().coeffRef(0,0)), rhs.stride(),
DestDerived::Flags&RowMajorBit, (Scalar*)&(res.coeffRef(0,0)), res.stride()
);
}
#endif // EIGEN_PRODUCT_H

View File

@@ -35,7 +35,7 @@ template<typename Lhs, typename Rhs,
? UpperTriangular
: -1,
int StorageOrder = ei_is_part<Lhs>::value ? -1 // this is to solve ambiguous specializations
: int(Lhs::Flags) & (RowMajorBit|SparseBit)
: int(Lhs::Flags) & int(RowMajorBit|SparseBit)
>
struct ei_solve_triangular_selector;

View File

@@ -69,7 +69,6 @@ template<typename MatrixType> class Transpose
inline int rows() const { return m_matrix.cols(); }
inline int cols() const { return m_matrix.rows(); }
inline int nonZeros() const { return m_matrix.nonZeros(); }
inline int stride(void) const { return m_matrix.stride(); }
inline Scalar& coeffRef(int row, int col)

View File

@@ -114,9 +114,15 @@ template<> EIGEN_STRONG_INLINE void ei_pstoreu<float>(float* to, const __m128&
template<> EIGEN_STRONG_INLINE void ei_pstoreu<double>(double* to, const __m128d& from) { _mm_storeu_pd(to, from); }
template<> EIGEN_STRONG_INLINE void ei_pstoreu<int>(int* to, const __m128i& from) { _mm_storeu_si128(reinterpret_cast<__m128i*>(to), from); }
#ifdef _MSC_VER
// this fix internal compilation error
template<> EIGEN_STRONG_INLINE float ei_pfirst<__m128>(const __m128& a) { float x = _mm_cvtss_f32(a); return x; }
#if defined(_MSC_VER) && (_MSC_VER <= 1500) && defined(_WIN64) && !defined(__INTEL_COMPILER)
// The temporary variable fixes an internal compilation error.
// Direct of the struct members fixed bug #62.
template<> EIGEN_STRONG_INLINE float ei_pfirst<__m128>(const __m128& a) { return a.m128_f32[0]; }
template<> EIGEN_STRONG_INLINE double ei_pfirst<__m128d>(const __m128d& a) { return a.m128d_f64[0]; }
template<> EIGEN_STRONG_INLINE int ei_pfirst<__m128i>(const __m128i& a) { int x = _mm_cvtsi128_si32(a); return x; }
#elif defined(_MSC_VER) && (_MSC_VER <= 1500) && !defined(__INTEL_COMPILER)
// The temporary variable fixes an internal compilation error.
template<> EIGEN_STRONG_INLINE float ei_pfirst<__m128>(const __m128& a) { float x = _mm_cvtss_f32(a); return x; }
template<> EIGEN_STRONG_INLINE double ei_pfirst<__m128d>(const __m128d& a) { double x = _mm_cvtsd_f64(a); return x; }
template<> EIGEN_STRONG_INLINE int ei_pfirst<__m128i>(const __m128i& a) { int x = _mm_cvtsi128_si32(a); return x; }
#else

View File

@@ -30,7 +30,7 @@
#define EIGEN_WORLD_VERSION 2
#define EIGEN_MAJOR_VERSION 0
#define EIGEN_MINOR_VERSION 7
#define EIGEN_MINOR_VERSION 15
#define EIGEN_VERSION_AT_LEAST(x,y,z) (EIGEN_WORLD_VERSION>x || (EIGEN_WORLD_VERSION>=x && \
(EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \
@@ -39,7 +39,7 @@
// 16 byte alignment is only useful for vectorization. Since it affects the ABI, we need to enable 16 byte alignment on all
// platforms where vectorization might be enabled. In theory we could always enable alignment, but it can be a cause of problems
// on some platforms, so we just disable it in certain common platform (compiler+architecture combinations) to avoid these problems.
#if defined(__GNUC__) && !(defined(__i386__) || defined(__x86_64__) || defined(__powerpc__) || defined(__ia64__))
#if defined(__GNUC__) && !(defined(__i386__) || defined(__x86_64__) || defined(__powerpc__) || defined(__ppc__) || defined(__ia64__))
#define EIGEN_GCC_AND_ARCH_DOESNT_WANT_ALIGNMENT 1
#else
#define EIGEN_GCC_AND_ARCH_DOESNT_WANT_ALIGNMENT 0
@@ -52,7 +52,7 @@
#endif
// FIXME vectorization + alignment is completely disabled with sun studio
#if !EIGEN_GCC_AND_ARCH_DOESNT_WANT_ALIGNMENT && !EIGEN_GCC3_OR_OLDER && !defined(__SUNPRO_CC)
#if !EIGEN_GCC_AND_ARCH_DOESNT_WANT_ALIGNMENT && !EIGEN_GCC3_OR_OLDER && !defined(__SUNPRO_CC) && !defined(__QNXNTO__)
#define EIGEN_ARCH_WANTS_ALIGNMENT 1
#else
#define EIGEN_ARCH_WANTS_ALIGNMENT 0
@@ -257,6 +257,9 @@ enum { RowsAtCompileTime = Eigen::ei_traits<Derived>::RowsAtCompileTime, \
_EIGEN_GENERIC_PUBLIC_INTERFACE(Derived, Eigen::MatrixBase<Derived>)
#define EIGEN_ENUM_MIN(a,b) (((int)a <= (int)b) ? (int)a : (int)b)
#define EIGEN_SIZE_MIN(a,b) (((int)a == 1 || (int)b == 1) ? 1 \
: ((int)a == Dynamic || (int)b == Dynamic) ? Dynamic \
: ((int)a <= (int)b) ? (int)a : (int)b)
#define EIGEN_ENUM_MAX(a,b) (((int)a >= (int)b) ? (int)a : (int)b)
// just an empty macro !

View File

@@ -43,7 +43,7 @@
#define EIGEN_MALLOC_ALREADY_ALIGNED 0
#endif
#if ((defined _GNU_SOURCE) || ((defined _XOPEN_SOURCE) && (_XOPEN_SOURCE >= 600))) && (defined _POSIX_ADVISORY_INFO) && (_POSIX_ADVISORY_INFO > 0)
#if (defined __QNXNTO__) || (((defined _GNU_SOURCE) || ((defined _XOPEN_SOURCE) && (_XOPEN_SOURCE >= 600))) && (defined _POSIX_ADVISORY_INFO) && (_POSIX_ADVISORY_INFO > 0))
#define EIGEN_HAS_POSIX_MEMALIGN 1
#else
#define EIGEN_HAS_POSIX_MEMALIGN 0
@@ -59,10 +59,10 @@
* Fast, but wastes 16 additional bytes of memory.
* Does not throw any exception.
*/
inline void* ei_handmade_aligned_malloc(size_t size)
inline void* ei_handmade_aligned_malloc(std::size_t size)
{
void *original = malloc(size+16);
void *aligned = reinterpret_cast<void*>((reinterpret_cast<size_t>(original) & ~(size_t(15))) + 16);
void *original = std::malloc(size+16);
void *aligned = reinterpret_cast<void*>((reinterpret_cast<std::size_t>(original) & ~(std::size_t(15))) + 16);
*(reinterpret_cast<void**>(aligned) - 1) = original;
return aligned;
}
@@ -71,13 +71,13 @@ inline void* ei_handmade_aligned_malloc(size_t size)
inline void ei_handmade_aligned_free(void *ptr)
{
if(ptr)
free(*(reinterpret_cast<void**>(ptr) - 1));
std::free(*(reinterpret_cast<void**>(ptr) - 1));
}
/** \internal allocates \a size bytes. The returned pointer is guaranteed to have 16 bytes alignment.
* On allocation error, the returned pointer is null, and if exceptions are enabled then a std::bad_alloc is thrown.
*/
inline void* ei_aligned_malloc(size_t size)
inline void* ei_aligned_malloc(std::size_t size)
{
#ifdef EIGEN_NO_MALLOC
ei_assert(false && "heap allocation is forbidden (EIGEN_NO_MALLOC is defined)");
@@ -85,9 +85,9 @@ inline void* ei_aligned_malloc(size_t size)
void *result;
#if !EIGEN_ALIGN
result = malloc(size);
result = std::malloc(size);
#elif EIGEN_MALLOC_ALREADY_ALIGNED
result = malloc(size);
result = std::malloc(size);
#elif EIGEN_HAS_POSIX_MEMALIGN
if(posix_memalign(&result, 16, size)) result = 0;
#elif EIGEN_HAS_MM_MALLOC
@@ -108,18 +108,18 @@ inline void* ei_aligned_malloc(size_t size)
/** allocates \a size bytes. If Align is true, then the returned ptr is 16-byte-aligned.
* On allocation error, the returned pointer is null, and if exceptions are enabled then a std::bad_alloc is thrown.
*/
template<bool Align> inline void* ei_conditional_aligned_malloc(size_t size)
template<bool Align> inline void* ei_conditional_aligned_malloc(std::size_t size)
{
return ei_aligned_malloc(size);
}
template<> inline void* ei_conditional_aligned_malloc<false>(size_t size)
template<> inline void* ei_conditional_aligned_malloc<false>(std::size_t size)
{
#ifdef EIGEN_NO_MALLOC
ei_assert(false && "heap allocation is forbidden (EIGEN_NO_MALLOC is defined)");
#endif
void *result = malloc(size);
void *result = std::malloc(size);
#ifdef EIGEN_EXCEPTIONS
if(!result) throw std::bad_alloc();
#endif
@@ -129,9 +129,9 @@ template<> inline void* ei_conditional_aligned_malloc<false>(size_t size)
/** \internal construct the elements of an array.
* The \a size parameter tells on how many objects to call the constructor of T.
*/
template<typename T> inline T* ei_construct_elements_of_array(T *ptr, size_t size)
template<typename T> inline T* ei_construct_elements_of_array(T *ptr, std::size_t size)
{
for (size_t i=0; i < size; ++i) ::new (ptr + i) T;
for (std::size_t i=0; i < size; ++i) ::new (ptr + i) T;
return ptr;
}
@@ -139,13 +139,13 @@ template<typename T> inline T* ei_construct_elements_of_array(T *ptr, size_t siz
* On allocation error, the returned pointer is undefined, but if exceptions are enabled then a std::bad_alloc is thrown.
* The default constructor of T is called.
*/
template<typename T> inline T* ei_aligned_new(size_t size)
template<typename T> inline T* ei_aligned_new(std::size_t size)
{
T *result = reinterpret_cast<T*>(ei_aligned_malloc(sizeof(T)*size));
return ei_construct_elements_of_array(result, size);
}
template<typename T, bool Align> inline T* ei_conditional_aligned_new(size_t size)
template<typename T, bool Align> inline T* ei_conditional_aligned_new(std::size_t size)
{
T *result = reinterpret_cast<T*>(ei_conditional_aligned_malloc<Align>(sizeof(T)*size));
return ei_construct_elements_of_array(result, size);
@@ -156,11 +156,11 @@ template<typename T, bool Align> inline T* ei_conditional_aligned_new(size_t siz
inline void ei_aligned_free(void *ptr)
{
#if !EIGEN_ALIGN
free(ptr);
std::free(ptr);
#elif EIGEN_MALLOC_ALREADY_ALIGNED
free(ptr);
std::free(ptr);
#elif EIGEN_HAS_POSIX_MEMALIGN
free(ptr);
std::free(ptr);
#elif EIGEN_HAS_MM_MALLOC
_mm_free(ptr);
#elif defined(_MSC_VER)
@@ -179,22 +179,23 @@ template<bool Align> inline void ei_conditional_aligned_free(void *ptr)
template<> inline void ei_conditional_aligned_free<false>(void *ptr)
{
free(ptr);
std::free(ptr);
}
/** \internal destruct the elements of an array.
* The \a size parameters tells on how many objects to call the destructor of T.
*/
template<typename T> inline void ei_destruct_elements_of_array(T *ptr, size_t size)
template<typename T> inline void ei_destruct_elements_of_array(T *ptr, std::size_t size)
{
// always destruct an array starting from the end.
while(size) ptr[--size].~T();
if(ptr)
while(size) ptr[--size].~T();
}
/** \internal delete objects constructed with ei_aligned_new
* The \a size parameters tells on how many objects to call the destructor of T.
*/
template<typename T> inline void ei_aligned_delete(T *ptr, size_t size)
template<typename T> inline void ei_aligned_delete(T *ptr, std::size_t size)
{
ei_destruct_elements_of_array<T>(ptr, size);
ei_aligned_free(ptr);
@@ -203,24 +204,53 @@ template<typename T> inline void ei_aligned_delete(T *ptr, size_t size)
/** \internal delete objects constructed with ei_conditional_aligned_new
* The \a size parameters tells on how many objects to call the destructor of T.
*/
template<typename T, bool Align> inline void ei_conditional_aligned_delete(T *ptr, size_t size)
template<typename T, bool Align> inline void ei_conditional_aligned_delete(T *ptr, std::size_t size)
{
ei_destruct_elements_of_array<T>(ptr, size);
ei_conditional_aligned_free<Align>(ptr);
}
/** \internal \returns the number of elements which have to be skipped such that data are 16 bytes aligned */
template<typename Scalar>
inline static int ei_alignmentOffset(const Scalar* ptr, int maxOffset)
/** \internal \returns the index of the first element of the array that is well aligned for vectorization.
*
* \param array the address of the start of the array
* \param size the size of the array
*
* \note If no element of the array is well aligned, the size of the array is returned. Typically,
* for example with SSE, "well aligned" means 16-byte-aligned. If vectorization is disabled or if the
* packet size for the given scalar type is 1, then everything is considered well-aligned.
*
* \note If the scalar type is vectorizable, we rely on the following assumptions: sizeof(Scalar) is a
* power of 2, the packet size in bytes is also a power of 2, and is a multiple of sizeof(Scalar). On the
* other hand, we do not assume that the array address is a multiple of sizeof(Scalar), as that fails for
* example with Scalar=double on certain 32-bit platforms, see bug #79.
*
* There is also the variant ei_first_aligned(const MatrixBase&, Integer) defined in Coeffs.h.
*/
template<typename Scalar, typename Integer>
inline static Integer ei_alignmentOffset(const Scalar* array, Integer size)
{
typedef typename ei_packet_traits<Scalar>::type Packet;
const int PacketSize = ei_packet_traits<Scalar>::size;
const int PacketAlignedMask = PacketSize-1;
const bool Vectorized = PacketSize>1;
return Vectorized
? std::min<int>( (PacketSize - (int((size_t(ptr)/sizeof(Scalar))) & PacketAlignedMask))
& PacketAlignedMask, maxOffset)
: 0;
enum { PacketSize = ei_packet_traits<Scalar>::size,
PacketAlignedMask = PacketSize-1
};
if(PacketSize==1)
{
// Either there is no vectorization, or a packet consists of exactly 1 scalar so that all elements
// of the array have the same aligment.
return 0;
}
else if(std::size_t(array) & (sizeof(Scalar)-1))
{
// There is vectorization for this scalar type, but the array is not aligned to the size of a single scalar.
// Consequently, no element of the array is well aligned.
return size;
}
else
{
return std::min<Integer>( (PacketSize - (Integer((std::size_t(array)/sizeof(Scalar))) & PacketAlignedMask))
& PacketAlignedMask, size);
}
}
/** \internal
@@ -252,23 +282,23 @@ inline static int ei_alignmentOffset(const Scalar* ptr, int maxOffset)
#if EIGEN_ALIGN
#ifdef EIGEN_EXCEPTIONS
#define EIGEN_MAKE_ALIGNED_OPERATOR_NEW_NOTHROW(NeedsToAlign) \
void* operator new(size_t size, const std::nothrow_t&) throw() { \
void* operator new(std::size_t size, const std::nothrow_t&) throw() { \
try { return Eigen::ei_conditional_aligned_malloc<NeedsToAlign>(size); } \
catch (...) { return 0; } \
return 0; \
}
#else
#define EIGEN_MAKE_ALIGNED_OPERATOR_NEW_NOTHROW(NeedsToAlign) \
void* operator new(size_t size, const std::nothrow_t&) throw() { \
void* operator new(std::size_t size, const std::nothrow_t&) throw() { \
return Eigen::ei_conditional_aligned_malloc<NeedsToAlign>(size); \
}
#endif
#define EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(NeedsToAlign) \
void *operator new(size_t size) { \
void *operator new(std::size_t size) { \
return Eigen::ei_conditional_aligned_malloc<NeedsToAlign>(size); \
} \
void *operator new[](size_t size) { \
void *operator new[](std::size_t size) { \
return Eigen::ei_conditional_aligned_malloc<NeedsToAlign>(size); \
} \
void operator delete(void * ptr) throw() { Eigen::ei_conditional_aligned_free<NeedsToAlign>(ptr); } \
@@ -276,7 +306,7 @@ inline static int ei_alignmentOffset(const Scalar* ptr, int maxOffset)
/* in-place new and delete. since (at least afaik) there is no actual */ \
/* memory allocated we can safely let the default implementation handle */ \
/* this particular case. */ \
static void *operator new(size_t size, void *ptr) { return ::operator new(size,ptr); } \
static void *operator new(std::size_t size, void *ptr) { return ::operator new(size,ptr); } \
void operator delete(void * memory, void *ptr) throw() { return ::operator delete(memory,ptr); } \
/* nothrow-new (returns zero instead of std::bad_alloc) */ \
EIGEN_MAKE_ALIGNED_OPERATOR_NEW_NOTHROW(NeedsToAlign) \
@@ -310,8 +340,8 @@ template<class T>
class aligned_allocator
{
public:
typedef size_t size_type;
typedef ptrdiff_t difference_type;
typedef std::size_t size_type;
typedef std::ptrdiff_t difference_type;
typedef T* pointer;
typedef const T* const_pointer;
typedef T& reference;

View File

@@ -74,6 +74,7 @@
THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES,
THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES,
INVALID_MATRIX_TEMPLATE_PARAMETERS,
INVALID_MATRIXBASE_TEMPLATE_PARAMETERS,
BOTH_MATRICES_MUST_HAVE_THE_SAME_STORAGE_ORDER,
THIS_METHOD_IS_ONLY_FOR_DIAGONAL_MATRIX
};

View File

@@ -94,10 +94,16 @@ class ei_compute_matrix_flags
{
enum {
row_major_bit = Options&RowMajor ? RowMajorBit : 0,
inner_max_size = row_major_bit ? MaxCols : MaxRows,
inner_max_size = int(MaxRows==1) ? int(MaxCols)
: int(MaxCols==1) ? int(MaxRows)
: int(row_major_bit) ? int(MaxCols) : int(MaxRows),
is_big = inner_max_size == Dynamic,
is_packet_size_multiple = (Cols*Rows) % ei_packet_traits<Scalar>::size == 0,
aligned_bit = ((Options&AutoAlign) && (is_big || is_packet_size_multiple)) ? AlignedBit : 0,
storage_has_fixed_size = MaxRows != Dynamic && MaxCols != Dynamic,
storage_has_aligned_fixed_size = storage_has_fixed_size
&& ( (MaxCols*MaxRows) % ei_packet_traits<Scalar>::size == 0 ),
aligned_bit = ( (Options&AutoAlign)
&& (is_big || storage_has_aligned_fixed_size)
) ? AlignedBit : 0,
packet_access_bit = ei_packet_traits<Scalar>::size > 1 && aligned_bit ? PacketAccessBit : 0
};
@@ -161,6 +167,19 @@ template<typename T> struct ei_plain_matrix_type_column_major
> type;
};
/* ei_plain_matrix_type_row_major : same as ei_plain_matrix_type but guaranteed to be row-major
*/
template<typename T> struct ei_plain_matrix_type_row_major
{
typedef Matrix<typename ei_traits<T>::Scalar,
ei_traits<T>::RowsAtCompileTime,
ei_traits<T>::ColsAtCompileTime,
AutoAlign | RowMajor,
ei_traits<T>::MaxRowsAtCompileTime,
ei_traits<T>::MaxColsAtCompileTime
> type;
};
template<typename T> struct ei_must_nest_by_value { enum { ret = false }; };
template<typename T> struct ei_must_nest_by_value<NestByValue<T> > { enum { ret = true }; };

View File

@@ -52,9 +52,9 @@ public:
typedef _Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef Matrix<Scalar,AmbientDimAtCompileTime,1> VectorType;
typedef Matrix<Scalar,AmbientDimAtCompileTime==Dynamic
typedef Matrix<Scalar,int(AmbientDimAtCompileTime)==Dynamic
? Dynamic
: AmbientDimAtCompileTime+1,1> Coefficients;
: int(AmbientDimAtCompileTime)+1,1> Coefficients;
typedef Block<Coefficients,AmbientDimAtCompileTime,1> NormalReturnType;
/** Default constructor without initialization */

View File

@@ -450,22 +450,31 @@ inline Scalar Quaternion<Scalar>::angularDistance(const Quaternion& other) const
template <typename Scalar>
Quaternion<Scalar> Quaternion<Scalar>::slerp(Scalar t, const Quaternion& other) const
{
static const Scalar one = Scalar(1) - precision<Scalar>();
static const Scalar one = Scalar(1) - machine_epsilon<Scalar>();
Scalar d = this->dot(other);
Scalar absD = ei_abs(d);
Scalar scale0;
Scalar scale1;
if (absD>=one)
return *this;
{
scale0 = Scalar(1) - t;
scale1 = t;
}
else
{
// theta is the angle between the 2 quaternions
Scalar theta = std::acos(absD);
Scalar sinTheta = ei_sin(theta);
// theta is the angle between the 2 quaternions
Scalar theta = std::acos(absD);
Scalar sinTheta = ei_sin(theta);
scale0 = ei_sin( ( Scalar(1) - t ) * theta) / sinTheta;
scale1 = ei_sin( ( t * theta) ) / sinTheta;
if (d<0)
scale1 = -scale1;
}
Scalar scale0 = ei_sin( ( Scalar(1) - t ) * theta) / sinTheta;
Scalar scale1 = ei_sin( ( t * theta) ) / sinTheta;
if (d<0)
scale1 = -scale1;
return Quaternion(scale0 * m_coeffs + scale1 * other.m_coeffs);
return Quaternion<Scalar>(scale0 * coeffs() + scale1 * other.coeffs());
}
// set from a rotation matrix

View File

@@ -29,8 +29,8 @@
*** Part 1 : optimized implementations for fixed-size 2,3,4 cases ***
********************************************************************/
template<typename MatrixType>
void ei_compute_inverse_in_size2_case(const MatrixType& matrix, MatrixType* result)
template<typename XprType, typename MatrixType>
void ei_compute_inverse_in_size2_case(const XprType& matrix, MatrixType* result)
{
typedef typename MatrixType::Scalar Scalar;
const Scalar invdet = Scalar(1) / matrix.determinant();
@@ -54,10 +54,10 @@ bool ei_compute_inverse_in_size2_case_with_check(const XprType& matrix, MatrixTy
return true;
}
template<typename MatrixType>
void ei_compute_inverse_in_size3_case(const MatrixType& matrix, MatrixType* result)
template<typename Derived, typename OtherDerived>
void ei_compute_inverse_in_size3_case(const Derived& matrix, OtherDerived* result)
{
typedef typename MatrixType::Scalar Scalar;
typedef typename Derived::Scalar Scalar;
const Scalar det_minor00 = matrix.minor(0,0).determinant();
const Scalar det_minor10 = matrix.minor(1,0).determinant();
const Scalar det_minor20 = matrix.minor(2,0).determinant();
@@ -75,140 +75,204 @@ void ei_compute_inverse_in_size3_case(const MatrixType& matrix, MatrixType* resu
result->coeffRef(2, 2) = matrix.minor(2,2).determinant() * invdet;
}
template<typename MatrixType>
bool ei_compute_inverse_in_size4_case_helper(const MatrixType& matrix, MatrixType* result)
template<typename Derived, typename OtherDerived, typename Scalar = typename Derived::Scalar>
struct ei_compute_inverse_in_size4_case
{
/* Let's split M into four 2x2 blocks:
* (P Q)
* (R S)
* If P is invertible, with inverse denoted by P_inverse, and if
* (S - R*P_inverse*Q) is also invertible, then the inverse of M is
* (P' Q')
* (R' S')
* where
* S' = (S - R*P_inverse*Q)^(-1)
* P' = P1 + (P1*Q) * S' *(R*P_inverse)
* Q' = -(P_inverse*Q) * S'
* R' = -S' * (R*P_inverse)
*/
typedef Block<MatrixType,2,2> XprBlock22;
typedef typename MatrixBase<XprBlock22>::PlainMatrixType Block22;
Block22 P_inverse;
if(ei_compute_inverse_in_size2_case_with_check(matrix.template block<2,2>(0,0), &P_inverse))
static void run(const Derived& matrix, OtherDerived& result)
{
const Block22 Q = matrix.template block<2,2>(0,2);
const Block22 P_inverse_times_Q = P_inverse * Q;
const XprBlock22 R = matrix.template block<2,2>(2,0);
const Block22 R_times_P_inverse = R * P_inverse;
const Block22 R_times_P_inverse_times_Q = R_times_P_inverse * Q;
const XprBlock22 S = matrix.template block<2,2>(2,2);
const Block22 X = S - R_times_P_inverse_times_Q;
Block22 Y;
ei_compute_inverse_in_size2_case(X, &Y);
result->template block<2,2>(2,2) = Y;
result->template block<2,2>(2,0) = - Y * R_times_P_inverse;
const Block22 Z = P_inverse_times_Q * Y;
result->template block<2,2>(0,2) = - Z;
result->template block<2,2>(0,0) = P_inverse + Z * R_times_P_inverse;
return true;
result.coeffRef(0,0) = matrix.minor(0,0).determinant();
result.coeffRef(1,0) = -matrix.minor(0,1).determinant();
result.coeffRef(2,0) = matrix.minor(0,2).determinant();
result.coeffRef(3,0) = -matrix.minor(0,3).determinant();
result.coeffRef(0,2) = matrix.minor(2,0).determinant();
result.coeffRef(1,2) = -matrix.minor(2,1).determinant();
result.coeffRef(2,2) = matrix.minor(2,2).determinant();
result.coeffRef(3,2) = -matrix.minor(2,3).determinant();
result.coeffRef(0,1) = -matrix.minor(1,0).determinant();
result.coeffRef(1,1) = matrix.minor(1,1).determinant();
result.coeffRef(2,1) = -matrix.minor(1,2).determinant();
result.coeffRef(3,1) = matrix.minor(1,3).determinant();
result.coeffRef(0,3) = -matrix.minor(3,0).determinant();
result.coeffRef(1,3) = matrix.minor(3,1).determinant();
result.coeffRef(2,3) = -matrix.minor(3,2).determinant();
result.coeffRef(3,3) = matrix.minor(3,3).determinant();
result /= (matrix.col(0).cwise()*result.row(0).transpose()).sum();
}
else
{
return false;
}
}
};
template<typename MatrixType>
void ei_compute_inverse_in_size4_case(const MatrixType& matrix, MatrixType* result)
#ifdef EIGEN_VECTORIZE_SSE
// The SSE code for the 4x4 float matrix inverse in this file comes from the file
// ftp://download.intel.com/design/PentiumIII/sml/24504301.pdf
// its copyright information is:
// Copyright (C) 1999 Intel Corporation
// See page ii of that document for legal stuff. Not being lawyers, we just assume
// here that if Intel makes this document publically available, with source code
// and detailed explanations, it's because they want their CPUs to be fed with
// good code, and therefore they presumably don't mind us using it in Eigen.
template<typename Derived, typename OtherDerived>
struct ei_compute_inverse_in_size4_case<Derived, OtherDerived, float>
{
if(ei_compute_inverse_in_size4_case_helper(matrix, result))
static void run(const Derived& matrix, OtherDerived& result)
{
// good ! The topleft 2x2 block was invertible, so the 2x2 blocks approach is successful.
return;
// Variables (Streaming SIMD Extensions registers) which will contain cofactors and, later, the
// lines of the inverted matrix.
__m128 minor0, minor1, minor2, minor3;
// Variables which will contain the lines of the reference matrix and, later (after the transposition),
// the columns of the original matrix.
__m128 row0, row1, row2, row3;
// Temporary variables and the variable that will contain the matrix determinant.
__m128 det, tmp1;
// Matrix transposition
const float *src = matrix.data();
tmp1 = _mm_loadh_pi(_mm_castpd_ps(_mm_load_sd((double*)src)), (__m64*)(src+ 4));
row1 = _mm_loadh_pi(_mm_castpd_ps(_mm_load_sd((double*)(src+8))), (__m64*)(src+12));
row0 = _mm_shuffle_ps(tmp1, row1, 0x88);
row1 = _mm_shuffle_ps(row1, tmp1, 0xDD);
tmp1 = _mm_loadh_pi(_mm_castpd_ps(_mm_load_sd((double*)(src+ 2))), (__m64*)(src+ 6));
row3 = _mm_loadh_pi(_mm_castpd_ps(_mm_load_sd((double*)(src+10))), (__m64*)(src+14));
row2 = _mm_shuffle_ps(tmp1, row3, 0x88);
row3 = _mm_shuffle_ps(row3, tmp1, 0xDD);
// Cofactors calculation. Because in the process of cofactor computation some pairs in three-
// element products are repeated, it is not reasonable to load these pairs anew every time. The
// values in the registers with these pairs are formed using shuffle instruction. Cofactors are
// calculated row by row (4 elements are placed in 1 SP FP SIMD floating point register).
tmp1 = _mm_mul_ps(row2, row3);
tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1);
minor0 = _mm_mul_ps(row1, tmp1);
minor1 = _mm_mul_ps(row0, tmp1);
tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E);
minor0 = _mm_sub_ps(_mm_mul_ps(row1, tmp1), minor0);
minor1 = _mm_sub_ps(_mm_mul_ps(row0, tmp1), minor1);
minor1 = _mm_shuffle_ps(minor1, minor1, 0x4E);
// -----------------------------------------------
tmp1 = _mm_mul_ps(row1, row2);
tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1);
minor0 = _mm_add_ps(_mm_mul_ps(row3, tmp1), minor0);
minor3 = _mm_mul_ps(row0, tmp1);
tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E);
minor0 = _mm_sub_ps(minor0, _mm_mul_ps(row3, tmp1));
minor3 = _mm_sub_ps(_mm_mul_ps(row0, tmp1), minor3);
minor3 = _mm_shuffle_ps(minor3, minor3, 0x4E);
// -----------------------------------------------
tmp1 = _mm_mul_ps(_mm_shuffle_ps(row1, row1, 0x4E), row3);
tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1);
row2 = _mm_shuffle_ps(row2, row2, 0x4E);
minor0 = _mm_add_ps(_mm_mul_ps(row2, tmp1), minor0);
minor2 = _mm_mul_ps(row0, tmp1);
tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E);
minor0 = _mm_sub_ps(minor0, _mm_mul_ps(row2, tmp1));
minor2 = _mm_sub_ps(_mm_mul_ps(row0, tmp1), minor2);
minor2 = _mm_shuffle_ps(minor2, minor2, 0x4E);
// -----------------------------------------------
tmp1 = _mm_mul_ps(row0, row1);
tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1);
minor2 = _mm_add_ps(_mm_mul_ps(row3, tmp1), minor2);
minor3 = _mm_sub_ps(_mm_mul_ps(row2, tmp1), minor3);
tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E);
minor2 = _mm_sub_ps(_mm_mul_ps(row3, tmp1), minor2);
minor3 = _mm_sub_ps(minor3, _mm_mul_ps(row2, tmp1));
// -----------------------------------------------
tmp1 = _mm_mul_ps(row0, row3);
tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1);
minor1 = _mm_sub_ps(minor1, _mm_mul_ps(row2, tmp1));
minor2 = _mm_add_ps(_mm_mul_ps(row1, tmp1), minor2);
tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E);
minor1 = _mm_add_ps(_mm_mul_ps(row2, tmp1), minor1);
minor2 = _mm_sub_ps(minor2, _mm_mul_ps(row1, tmp1));
// -----------------------------------------------
tmp1 = _mm_mul_ps(row0, row2);
tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1);
minor1 = _mm_add_ps(_mm_mul_ps(row3, tmp1), minor1);
minor3 = _mm_sub_ps(minor3, _mm_mul_ps(row1, tmp1));
tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E);
minor1 = _mm_sub_ps(minor1, _mm_mul_ps(row3, tmp1));
minor3 = _mm_add_ps(_mm_mul_ps(row1, tmp1), minor3);
// Evaluation of determinant and its reciprocal value. In the original Intel document,
// 1/det was evaluated using a fast rcpps command with subsequent approximation using
// the Newton-Raphson algorithm. Here, we go for a IEEE-compliant division instead,
// so as to not compromise precision at all.
det = _mm_mul_ps(row0, minor0);
det = _mm_add_ps(_mm_shuffle_ps(det, det, 0x4E), det);
det = _mm_add_ss(_mm_shuffle_ps(det, det, 0xB1), det);
// tmp1= _mm_rcp_ss(det);
// det= _mm_sub_ss(_mm_add_ss(tmp1, tmp1), _mm_mul_ss(det, _mm_mul_ss(tmp1, tmp1)));
det = _mm_div_ss(_mm_set_ss(1.0f), det); // <--- yay, one original line not copied from Intel
det = _mm_shuffle_ps(det, det, 0x00);
// warning, Intel's variable naming is very confusing: now 'det' is 1/det !
// Multiplication of cofactors by 1/det. Storing the inverse matrix to the address in pointer src.
minor0 = _mm_mul_ps(det, minor0);
float *dst = result.data();
_mm_storel_pi((__m64*)(dst), minor0);
_mm_storeh_pi((__m64*)(dst+2), minor0);
minor1 = _mm_mul_ps(det, minor1);
_mm_storel_pi((__m64*)(dst+4), minor1);
_mm_storeh_pi((__m64*)(dst+6), minor1);
minor2 = _mm_mul_ps(det, minor2);
_mm_storel_pi((__m64*)(dst+ 8), minor2);
_mm_storeh_pi((__m64*)(dst+10), minor2);
minor3 = _mm_mul_ps(det, minor3);
_mm_storel_pi((__m64*)(dst+12), minor3);
_mm_storeh_pi((__m64*)(dst+14), minor3);
}
else
{
// rare case: the topleft 2x2 block is not invertible (but the matrix itself is assumed to be).
// since this is a rare case, we don't need to optimize it. We just want to handle it with little
// additional code.
MatrixType m(matrix);
m.row(0).swap(m.row(2));
m.row(1).swap(m.row(3));
if(ei_compute_inverse_in_size4_case_helper(m, result))
{
// good, the topleft 2x2 block of m is invertible. Since m is different from matrix in that some
// rows were permuted, the actual inverse of matrix is derived from the inverse of m by permuting
// the corresponding columns.
result->col(0).swap(result->col(2));
result->col(1).swap(result->col(3));
}
else
{
// last possible case. Since matrix is assumed to be invertible, this last case has to work.
// first, undo the swaps previously made
m.row(0).swap(m.row(2));
m.row(1).swap(m.row(3));
// swap row 0 with the the row among 0 and 1 that has the biggest 2 first coeffs
int swap0with = ei_abs(m.coeff(0,0))+ei_abs(m.coeff(0,1))>ei_abs(m.coeff(1,0))+ei_abs(m.coeff(1,1)) ? 0 : 1;
m.row(0).swap(m.row(swap0with));
// swap row 1 with the the row among 2 and 3 that has the biggest 2 first coeffs
int swap1with = ei_abs(m.coeff(2,0))+ei_abs(m.coeff(2,1))>ei_abs(m.coeff(3,0))+ei_abs(m.coeff(3,1)) ? 2 : 3;
m.row(1).swap(m.row(swap1with));
ei_compute_inverse_in_size4_case_helper(m, result);
result->col(1).swap(result->col(swap1with));
result->col(0).swap(result->col(swap0with));
}
}
}
};
#endif
/***********************************************
*** Part 2 : selector and MatrixBase methods ***
***********************************************/
template<typename MatrixType, int Size = MatrixType::RowsAtCompileTime>
template<typename Derived, typename OtherDerived, int Size = Derived::RowsAtCompileTime>
struct ei_compute_inverse
{
static inline void run(const MatrixType& matrix, MatrixType* result)
static inline void run(const Derived& matrix, OtherDerived* result)
{
LU<MatrixType> lu(matrix);
LU<Derived> lu(matrix);
lu.computeInverse(result);
}
};
template<typename MatrixType>
struct ei_compute_inverse<MatrixType, 1>
template<typename Derived, typename OtherDerived>
struct ei_compute_inverse<Derived, OtherDerived, 1>
{
static inline void run(const MatrixType& matrix, MatrixType* result)
static inline void run(const Derived& matrix, OtherDerived* result)
{
typedef typename MatrixType::Scalar Scalar;
typedef typename Derived::Scalar Scalar;
result->coeffRef(0,0) = Scalar(1) / matrix.coeff(0,0);
}
};
template<typename MatrixType>
struct ei_compute_inverse<MatrixType, 2>
template<typename Derived, typename OtherDerived>
struct ei_compute_inverse<Derived, OtherDerived, 2>
{
static inline void run(const MatrixType& matrix, MatrixType* result)
static inline void run(const Derived& matrix, OtherDerived* result)
{
ei_compute_inverse_in_size2_case(matrix, result);
}
};
template<typename MatrixType>
struct ei_compute_inverse<MatrixType, 3>
template<typename Derived, typename OtherDerived>
struct ei_compute_inverse<Derived, OtherDerived, 3>
{
static inline void run(const MatrixType& matrix, MatrixType* result)
static inline void run(const Derived& matrix, OtherDerived* result)
{
ei_compute_inverse_in_size3_case(matrix, result);
}
};
template<typename MatrixType>
struct ei_compute_inverse<MatrixType, 4>
template<typename Derived, typename OtherDerived>
struct ei_compute_inverse<Derived, OtherDerived, 4>
{
static inline void run(const MatrixType& matrix, MatrixType* result)
static inline void run(const Derived& matrix, OtherDerived* result)
{
ei_compute_inverse_in_size4_case(matrix, result);
ei_compute_inverse_in_size4_case<Derived, OtherDerived>::run(matrix, *result);
}
};
@@ -226,11 +290,12 @@ struct ei_compute_inverse<MatrixType, 4>
* \sa inverse()
*/
template<typename Derived>
inline void MatrixBase<Derived>::computeInverse(PlainMatrixType *result) const
template<typename OtherDerived>
inline void MatrixBase<Derived>::computeInverse(MatrixBase<OtherDerived> *result) const
{
ei_assert(rows() == cols());
EIGEN_STATIC_ASSERT(NumTraits<Scalar>::HasFloatingPoint,NUMERIC_TYPE_MUST_BE_FLOATING_POINT)
ei_compute_inverse<PlainMatrixType>::run(eval(), result);
ei_compute_inverse<PlainMatrixType, OtherDerived>::run(eval(), static_cast<OtherDerived*>(result));
}
/** \lu_module

View File

@@ -63,12 +63,12 @@ template<typename MatrixType> class LU
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType;
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVectorType;
typedef Matrix<int, 1, MatrixType::ColsAtCompileTime, MatrixType::Options, 1, MatrixType::MaxColsAtCompileTime> IntRowVectorType;
typedef Matrix<int, MatrixType::RowsAtCompileTime, 1, MatrixType::Options, MatrixType::MaxRowsAtCompileTime, 1> IntColVectorType;
typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime, MatrixType::Options, 1, MatrixType::MaxColsAtCompileTime> RowVectorType;
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1, MatrixType::Options, MatrixType::MaxRowsAtCompileTime, 1> ColVectorType;
enum { MaxSmallDimAtCompileTime = EIGEN_ENUM_MIN(
enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN(
MatrixType::MaxColsAtCompileTime,
MatrixType::MaxRowsAtCompileTime)
};
@@ -297,7 +297,8 @@ template<typename MatrixType> class LU
*
* \sa MatrixBase::computeInverse(), inverse()
*/
inline void computeInverse(MatrixType *result) const
template<typename ResultType>
inline void computeInverse(ResultType *result) const
{
solve(MatrixType::Identity(m_lu.rows(), m_lu.cols()), result);
}
@@ -508,15 +509,16 @@ bool LU<MatrixType>::solve(
if(!isSurjective())
{
// is c is in the image of U ?
RealScalar biggest_in_c = m_rank>0 ? c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff() : 0;
RealScalar biggest_in_c = m_rank>0 ? c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff() : RealScalar(0);
for(int col = 0; col < c.cols(); ++col)
for(int row = m_rank; row < c.rows(); ++row)
if(!ei_isMuchSmallerThan(c.coeff(row,col), biggest_in_c, m_precision))
return false;
}
m_lu.corner(TopLeft, m_rank, m_rank)
.template marked<UpperTriangular>()
.solveTriangularInPlace(c.corner(TopLeft, m_rank, c.cols()));
if(m_rank>0)
m_lu.corner(TopLeft, m_rank, m_rank)
.template marked<UpperTriangular>()
.solveTriangularInPlace(c.corner(TopLeft, m_rank, c.cols()));
// Step 4
result->resize(m_lu.cols(), b.cols());

View File

@@ -270,12 +270,18 @@ bool QR<MatrixType>::solve(
ei_assert(m_isInitialized && "QR is not initialized.");
const int rows = m_qr.rows();
ei_assert(b.rows() == rows);
result->resize(rows, b.cols());
// enforce the computation of the rank
rank();
result->resize(m_qr.cols(), b.cols());
// TODO(keir): There is almost certainly a faster way to multiply by
// Q^T without explicitly forming matrixQ(). Investigate.
*result = matrixQ().transpose()*b;
if(m_rank==0)
return result->isZero();
if(!isSurjective())
{
// is result is in the image of R ?

View File

@@ -293,7 +293,7 @@ void Tridiagonalization<MatrixType>::_compute(MatrixType& matA, CoeffVectorType&
{
int starti = i+1;
int alignedEnd = starti;
if (PacketSize>1)
if (PacketSize>1 && (int(MatrixType::Flags)&RowMajor) == 0)
{
int alignedStart = (starti) + ei_alignmentOffset(&matA.coeffRef(starti,j1), n-starti);
alignedEnd = alignedStart + ((n-alignedStart)/PacketSize)*PacketSize;

View File

@@ -49,7 +49,7 @@ template<typename MatrixType> class SVD
enum {
PacketSize = ei_packet_traits<Scalar>::size,
AlignmentMask = int(PacketSize)-1,
MinSize = EIGEN_ENUM_MIN(MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime)
MinSize = EIGEN_SIZE_MIN(MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime)
};
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVector;
@@ -61,6 +61,8 @@ template<typename MatrixType> class SVD
public:
SVD() {} // a user who relied on compiler-generated default compiler reported problems with MSVC in 2.0.7
SVD(const MatrixType& matrix)
: m_matU(matrix.rows(), std::min(matrix.rows(), matrix.cols())),
m_matV(matrix.cols(),matrix.cols()),
@@ -108,6 +110,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
const int n = matrix.cols();
const int nu = std::min(m,n);
ei_assert(m>=n && "In Eigen 2.0, SVD only works for MxN matrices with M>=N. Sorry!");
ei_assert(m>1 && "In Eigen 2.0, SVD doesn't work on 1x1 matrices");
m_matU.resize(m, nu);
m_matU.setZero();

View File

@@ -98,7 +98,7 @@ template<typename _Scalar> class AmbiVector
int allocSize = m_allocatedElements * sizeof(ListEl);
allocSize = allocSize/sizeof(Scalar) + (allocSize%sizeof(Scalar)>0?1:0);
Scalar* newBuffer = new Scalar[allocSize];
memcpy(newBuffer, m_buffer, copyElements * sizeof(ListEl));
std::memcpy(newBuffer, m_buffer, copyElements * sizeof(ListEl));
delete[] m_buffer;
m_buffer = newBuffer;
}

View File

@@ -37,7 +37,7 @@ class CompressedStorage
: m_values(0), m_indices(0), m_size(0), m_allocatedSize(0)
{}
CompressedStorage(size_t size)
CompressedStorage(std::size_t size)
: m_values(0), m_indices(0), m_size(0), m_allocatedSize(0)
{
resize(size);
@@ -52,8 +52,8 @@ class CompressedStorage
CompressedStorage& operator=(const CompressedStorage& other)
{
resize(other.size());
memcpy(m_values, other.m_values, m_size * sizeof(Scalar));
memcpy(m_indices, other.m_indices, m_size * sizeof(int));
std::memcpy(m_values, other.m_values, m_size * sizeof(Scalar));
std::memcpy(m_indices, other.m_indices, m_size * sizeof(int));
return *this;
}
@@ -71,9 +71,9 @@ class CompressedStorage
delete[] m_indices;
}
void reserve(size_t size)
void reserve(std::size_t size)
{
size_t newAllocatedSize = m_size + size;
std::size_t newAllocatedSize = m_size + size;
if (newAllocatedSize > m_allocatedSize)
reallocate(newAllocatedSize);
}
@@ -84,10 +84,10 @@ class CompressedStorage
reallocate(m_size);
}
void resize(size_t size, float reserveSizeFactor = 0)
void resize(std::size_t size, float reserveSizeFactor = 0)
{
if (m_allocatedSize<size)
reallocate(size + size_t(reserveSizeFactor*size));
reallocate(size + std::size_t(reserveSizeFactor*size));
m_size = size;
}
@@ -99,17 +99,17 @@ class CompressedStorage
m_indices[id] = i;
}
inline size_t size() const { return m_size; }
inline size_t allocatedSize() const { return m_allocatedSize; }
inline std::size_t size() const { return m_size; }
inline std::size_t allocatedSize() const { return m_allocatedSize; }
inline void clear() { m_size = 0; }
inline Scalar& value(size_t i) { return m_values[i]; }
inline const Scalar& value(size_t i) const { return m_values[i]; }
inline Scalar& value(std::size_t i) { return m_values[i]; }
inline const Scalar& value(std::size_t i) const { return m_values[i]; }
inline int& index(size_t i) { return m_indices[i]; }
inline const int& index(size_t i) const { return m_indices[i]; }
inline int& index(std::size_t i) { return m_indices[i]; }
inline const int& index(std::size_t i) const { return m_indices[i]; }
static CompressedStorage Map(int* indices, Scalar* values, size_t size)
static CompressedStorage Map(int* indices, Scalar* values, std::size_t size)
{
CompressedStorage res;
res.m_indices = indices;
@@ -125,11 +125,11 @@ class CompressedStorage
}
/** \returns the largest \c k in [start,end) such that for all \c j in [start,k) index[\c j]\<\a key */
inline int searchLowerIndex(size_t start, size_t end, int key) const
inline int searchLowerIndex(std::size_t start, std::size_t end, int key) const
{
while(end>start)
{
size_t mid = (end+start)>>1;
std::size_t mid = (end+start)>>1;
if (m_indices[mid]<key)
start = mid+1;
else
@@ -148,12 +148,12 @@ class CompressedStorage
return m_values[m_size-1];
// ^^ optimization: let's first check if it is the last coefficient
// (very common in high level algorithms)
const size_t id = searchLowerIndex(0,m_size-1,key);
const std::size_t id = searchLowerIndex(0,m_size-1,key);
return ((id<m_size) && (m_indices[id]==key)) ? m_values[id] : defaultValue;
}
/** Like at(), but the search is performed in the range [start,end) */
inline Scalar atInRange(size_t start, size_t end, int key, Scalar defaultValue = Scalar(0)) const
inline Scalar atInRange(std::size_t start, std::size_t end, int key, Scalar defaultValue = Scalar(0)) const
{
if (start==end)
return Scalar(0);
@@ -161,7 +161,7 @@ class CompressedStorage
return m_values[end-1];
// ^^ optimization: let's first check if it is the last coefficient
// (very common in high level algorithms)
const size_t id = searchLowerIndex(start,end-1,key);
const std::size_t id = searchLowerIndex(start,end-1,key);
return ((id<end) && (m_indices[id]==key)) ? m_values[id] : defaultValue;
}
@@ -170,11 +170,11 @@ class CompressedStorage
* such that the keys are sorted. */
inline Scalar& atWithInsertion(int key, Scalar defaultValue = Scalar(0))
{
size_t id = searchLowerIndex(0,m_size,key);
std::size_t id = searchLowerIndex(0,m_size,key);
if (id>=m_size || m_indices[id]!=key)
{
resize(m_size+1,1);
for (size_t j=m_size-1; j>id; --j)
for (std::size_t j=m_size-1; j>id; --j)
{
m_indices[j] = m_indices[j-1];
m_values[j] = m_values[j-1];
@@ -187,9 +187,9 @@ class CompressedStorage
void prune(Scalar reference, RealScalar epsilon = precision<RealScalar>())
{
size_t k = 0;
size_t n = size();
for (size_t i=0; i<n; ++i)
std::size_t k = 0;
std::size_t n = size();
for (std::size_t i=0; i<n; ++i)
{
if (!ei_isMuchSmallerThan(value(i), reference, epsilon))
{
@@ -203,14 +203,14 @@ class CompressedStorage
protected:
inline void reallocate(size_t size)
inline void reallocate(std::size_t size)
{
Scalar* newValues = new Scalar[size];
int* newIndices = new int[size];
size_t copySize = std::min(size, m_size);
std::size_t copySize = std::min(size, m_size);
// copy
memcpy(newValues, m_values, copySize * sizeof(Scalar));
memcpy(newIndices, m_indices, copySize * sizeof(int));
std::memcpy(newValues, m_values, copySize * sizeof(Scalar));
std::memcpy(newIndices, m_indices, copySize * sizeof(int));
// delete old stuff
delete[] m_values;
delete[] m_indices;
@@ -222,8 +222,8 @@ class CompressedStorage
protected:
Scalar* m_values;
int* m_indices;
size_t m_size;
size_t m_allocatedSize;
std::size_t m_size;
std::size_t m_allocatedSize;
};

View File

@@ -212,7 +212,7 @@ class DynamicSparseMatrix
// remove all coefficients with innerCoord>=innerSize
// TODO
std::cerr << "not implemented yet\n";
exit(2);
std::exit(2);
}
if (m_data.size() != outerSize)
{

View File

@@ -122,7 +122,7 @@ class SparseMatrix
{
m_data.clear();
//if (m_outerSize)
memset(m_outerIndex, 0, (m_outerSize+1)*sizeof(int));
std::memset(m_outerIndex, 0, (m_outerSize+1)*sizeof(int));
// for (int i=0; i<m_outerSize; ++i)
// m_outerIndex[i] = 0;
// if (m_outerSize)
@@ -164,7 +164,7 @@ class SparseMatrix
{
ei_assert(m_data.index(m_data.size()-1)<inner && "wrong sorted insertion");
}
assert(size_t(m_outerIndex[outer+1]) == m_data.size());
assert(std::size_t(m_outerIndex[outer+1]) == m_data.size());
int id = m_outerIndex[outer+1];
++m_outerIndex[outer+1];
@@ -190,10 +190,10 @@ class SparseMatrix
}
m_outerIndex[outer+1] = m_outerIndex[outer];
}
assert(size_t(m_outerIndex[outer+1]) == m_data.size() && "invalid outer index");
size_t startId = m_outerIndex[outer];
// FIXME let's make sure sizeof(long int) == sizeof(size_t)
size_t id = m_outerIndex[outer+1];
assert(std::size_t(m_outerIndex[outer+1]) == m_data.size() && "invalid outer index");
std::size_t startId = m_outerIndex[outer];
// FIXME let's make sure sizeof(long int) == sizeof(std::size_t)
std::size_t id = m_outerIndex[outer+1];
++m_outerIndex[outer+1];
float reallocRatio = 1;
@@ -273,7 +273,7 @@ class SparseMatrix
m_outerIndex = new int [outerSize+1];
m_outerSize = outerSize;
}
memset(m_outerIndex, 0, (m_outerSize+1)*sizeof(int));
std::memset(m_outerIndex, 0, (m_outerSize+1)*sizeof(int));
}
void resizeNonZeros(int size)
{
@@ -324,7 +324,7 @@ class SparseMatrix
else
{
resize(other.rows(), other.cols());
memcpy(m_outerIndex, other.m_outerIndex, (m_outerSize+1)*sizeof(int));
std::memcpy(m_outerIndex, other.m_outerIndex, (m_outerSize+1)*sizeof(int));
m_data = other.m_data;
}
return *this;

View File

@@ -97,7 +97,7 @@ struct ei_traits<SparseProduct<LhsNested, RhsNested, ProductMode> >
RowsAtCompileTime = _LhsNested::RowsAtCompileTime,
ColsAtCompileTime = _RhsNested::ColsAtCompileTime,
InnerSize = EIGEN_ENUM_MIN(_LhsNested::ColsAtCompileTime, _RhsNested::RowsAtCompileTime),
InnerSize = EIGEN_SIZE_MIN(_LhsNested::ColsAtCompileTime, _RhsNested::RowsAtCompileTime),
MaxRowsAtCompileTime = _LhsNested::MaxRowsAtCompileTime,
MaxColsAtCompileTime = _RhsNested::MaxColsAtCompileTime,

View File

@@ -43,8 +43,11 @@ struct ei_solve_triangular_selector<Lhs,Rhs,LowerTriangular,RowMajor|IsSparse>
{
lastVal = it.value();
lastIndex = it.index();
if(lastIndex == i)
break;
tmp -= lastVal * other.coeff(lastIndex,col);
}
if (Lhs::Flags & UnitDiagBit)
other.coeffRef(i,col) = tmp;
else

View File

@@ -1,8 +1,8 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra. Eigen itself is part of the KDE project.
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@@ -26,8 +26,15 @@
#ifndef EIGEN_BENCH_TIMER_H
#define EIGEN_BENCH_TIMER_H
#include <sys/time.h>
#if defined(_WIN32) || defined(__CYGWIN__)
#define NOMINMAX
#define WIN32_LEAN_AND_MEAN
#include <windows.h>
#else
#include <time.h>
#include <unistd.h>
#endif
#include <cstdlib>
#include <numeric>
@@ -35,12 +42,25 @@ namespace Eigen
{
/** Elapsed time timer keeping the best try.
*
* On POSIX platforms we use clock_gettime with CLOCK_PROCESS_CPUTIME_ID.
* On Windows we use QueryPerformanceCounter
*
* Important: on linux, you must link with -lrt
*/
class BenchTimer
{
public:
BenchTimer() { reset(); }
BenchTimer()
{
#if defined(_WIN32) || defined(__CYGWIN__)
LARGE_INTEGER freq;
QueryPerformanceFrequency(&freq);
m_frequency = (double)freq.QuadPart;
#endif
reset();
}
~BenchTimer() {}
@@ -51,23 +71,34 @@ public:
m_best = std::min(m_best, getTime() - m_start);
}
/** Return the best elapsed time.
/** Return the best elapsed time in seconds.
*/
inline double value(void)
{
return m_best;
return m_best;
}
#if defined(_WIN32) || defined(__CYGWIN__)
inline double getTime(void)
#else
static inline double getTime(void)
#endif
{
struct timeval tv;
struct timezone tz;
gettimeofday(&tv, &tz);
return (double)tv.tv_sec + 1.e-6 * (double)tv.tv_usec;
#ifdef WIN32
LARGE_INTEGER query_ticks;
QueryPerformanceCounter(&query_ticks);
return query_ticks.QuadPart/m_frequency;
#else
timespec ts;
clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &ts);
return double(ts.tv_sec) + 1e-9 * double(ts.tv_nsec);
#endif
}
protected:
#if defined(_WIN32) || defined(__CYGWIN__)
double m_frequency;
#endif
double m_best, m_start;
};

View File

@@ -142,7 +142,7 @@ public :
}
static inline void cholesky(const gene_matrix & X, gene_matrix & C, int N){
C = X.cholesky().matrixL();
C = X.llt().matrixL();
// C = X;
// Cholesky<gene_matrix>::computeInPlace(C);
// Cholesky<gene_matrix>::computeInPlaceBlock(C);

View File

@@ -0,0 +1,64 @@
# - Try to find how to link to the standard math library, if anything at all is needed to do.
# On most platforms this is automatic, but for example it's not automatic on QNX.
#
# Once done this will define
#
# STANDARD_MATH_LIBRARY_FOUND - we found how to successfully link to the standard math library
# STANDARD_MATH_LIBRARY - the name of the standard library that one has to link to.
# -- this will be left empty if it's automatic (most platforms).
# -- this will be set to "m" on platforms where one must explicitly
# pass the "-lm" linker flag.
#
# Copyright (c) 2010 Benoit Jacob <jacob.benoit.1@gmail.com>
# Redistribution and use is allowed according to the terms of the 2-clause BSD license.
include(CheckCXXSourceCompiles)
# a little test program for c++ math functions.
# notice the std:: is required on some platforms such as QNX
set(find_standard_math_library_test_program
"#include<cmath>
int main() { std::sin(0.0); std::log(0.0f); }")
# first try compiling/linking the test program without any linker flags
set(CMAKE_REQUIRED_FLAGS "")
set(CMAKE_REQUIRED_LIBRARIES "")
CHECK_CXX_SOURCE_COMPILES(
"${find_standard_math_library_test_program}"
standard_math_library_linked_to_automatically
)
if(standard_math_library_linked_to_automatically)
# the test program linked successfully without any linker flag.
set(STANDARD_MATH_LIBRARY "")
set(STANDARD_MATH_LIBRARY_FOUND TRUE)
else()
# the test program did not link successfully without any linker flag.
# This is a very uncommon case that so far we only saw on QNX. The next try is the
# standard name 'm' for the standard math library.
set(CMAKE_REQUIRED_LIBRARIES "m")
CHECK_CXX_SOURCE_COMPILES(
"${find_standard_math_library_test_program}"
standard_math_library_linked_to_as_m)
if(standard_math_library_linked_to_as_m)
# the test program linked successfully when linking to the 'm' library
set(STANDARD_MATH_LIBRARY "m")
set(STANDARD_MATH_LIBRARY_FOUND TRUE)
else()
# the test program still doesn't link successfully
set(STANDARD_MATH_LIBRARY_FOUND FALSE)
endif()
endif()

View File

@@ -57,10 +57,10 @@ void makeFloor(const MatrixBase<OtherDerived>& other) { derived() = derived().cw
template<typename OtherDerived>
void makeCeil(const MatrixBase<OtherDerived>& other) { derived() = derived().cwise().max(other.derived()); }
const typename Cwise<Derived>::ScalarAddReturnType
operator+(const Scalar& scalar) const { return cwise() + scalar }
const const CwiseUnaryOp<ei_scalar_add_op<Scalar>, Derived>
operator+(const Scalar& scalar) const { return cwise() + scalar; }
friend const typename Cwise<Derived>::ScalarAddReturnType
friend const CwiseUnaryOp<ei_scalar_add_op<Scalar>, Derived>
operator+(const Scalar& scalar, const MatrixBase<Derived>& mat) { return mat + scalar; }
\endcode

View File

@@ -4,6 +4,10 @@ namespace Eigen {
This is an issue that, so far, we met only with GCC on Windows: for instance, MinGW and TDM-GCC.
\htmlonly
<b><big><big>It seems that this GCC bug has been <a href="http://gcc.gnu.org/gcc-4.5/changes.html">fixed in GCC 4.5</a>. We encourage all GCC/Windows users to upgrade to GCC 4.5.</big></big></b>
\endhtmlonly
By default, in a function like this,
\code

View File

@@ -229,17 +229,24 @@ Of course, fixed-size matrices can't be resized.
\subsection TutorialMap Map
Any memory buffer can be mapped as an Eigen expression:
<table class="tutorial_code"><tr><td>
Any memory buffer can be mapped as an Eigen expression using the Map() static method:
\code
std::vector<float> stlarray(10);
Map<VectorXf>(&stlarray[0], stlarray.size()).setOnes();
int data[4] = 1, 2, 3, 4;
Matrix2i mat2x2(data);
MatrixXi mat2x2 = Map<Matrix2i>(data);
MatrixXi mat2x2 = Map<MatrixXi>(data,2,2);
VectorXf::Map(&stlarray[0], stlarray.size()).squaredNorm();
\endcode
Here VectorXf::Map returns an object of class Map<VectorXf>, which behaves like a VectorXf except that it uses the existing array. You can write to this object, that will write to the existing array. You can also construct a named obtect to reuse it:
\code
float array[rows*cols];
Map<MatrixXf> m(array,rows,cols);
m = othermatrix1 * othermatrix2;
m.eigenvalues();
\endcode
In the fixed-size case, no need to pass sizes:
\code
float array[9];
Map<Matrix3d> m(array);
Matrix3d::Map(array).setIdentity();
\endcode
</td></tr></table>

View File

@@ -11,7 +11,7 @@ namespace Eigen {
\section summary Executive summary
Using STL containers on \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types" requires taking the following two steps:
Using STL containers on \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types", or classes having members of such types, requires taking the following two steps:
\li A 16-byte-aligned allocator must be used. Eigen does provide one ready for use: aligned_allocator.
\li If you want to use the std::vector container, you need to \#include <Eigen/StdVector>.

View File

@@ -55,16 +55,17 @@ Note that here, Eigen::Vector2d is only used as an example, more generally the i
\section c2 Cause 2: STL Containers
If you use STL Containers such as std::vector, std::map, ..., with Eigen objects, like this,
If you use STL Containers such as std::vector, std::map, ..., with Eigen objects, or with classes containing Eigen objects, like this,
\code
std::vector<Eigen::Matrix2f> my_vector;
std::map<int, Eigen::Matrix2f> my_map;
struct my_class { ... Eigen::Matrix2f m; ... };
std::map<int, my_class> my_map;
\endcode
then you need to read this separate page: \ref StlContainers "Using STL Containers with Eigen".
Note that here, Eigen::Matrix2f is only used as an example, more generally the issue arises for all \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types".
Note that here, Eigen::Matrix2f is only used as an example, more generally the issue arises for all \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types" and \ref StructHavingEigenMembers "structures having such Eigen objects as member".
\section c3 Cause 3: Passing Eigen objects by value
@@ -106,7 +107,7 @@ However there are a few corner cases where these alignment settings get overridd
Two possibilities:
<ul>
<li>Define EIGEN_DONT_ALIGN. That disables all 128-bit alignment code, and in particular everything vectorization-related. But do note that this in particular breaks ABI compatibility with vectorized code.</li>
<li>Define EIGEN_DONT_ALIGN. That disables all 128-bit alignment code, and in particular everything vectorization-related. But do note that this in particular breaks ABI compatibility with vectorized code. This requires Eigen 2.0.6 or later, and with Eigen 2.0.12 or later it automatically implies EIGEN_DONT_VECTORIZE (before 2.0.12 you had to define both).</li>
<li>Or define both EIGEN_DONT_VECTORIZE and EIGEN_DISABLE_UNALIGNED_ARRAY_ASSERT. This keeps the 128-bit alignment code and thus preserves ABI compatibility.</li>
</ul>

View File

@@ -5,6 +5,9 @@ ADD_CUSTOM_TARGET(all_examples)
FOREACH(example_src ${examples_SRCS})
GET_FILENAME_COMPONENT(example ${example_src} NAME_WE)
ADD_EXECUTABLE(${example} ${example_src})
if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO)
target_link_libraries(${example} ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO})
endif()
GET_TARGET_PROPERTY(example_executable
${example} LOCATION)
ADD_CUSTOM_COMMAND(

View File

@@ -11,6 +11,9 @@ CONFIGURE_FILE(${CMAKE_CURRENT_SOURCE_DIR}/compile_snippet.cpp.in
${CMAKE_CURRENT_BINARY_DIR}/${compile_snippet_src})
ADD_EXECUTABLE(${compile_snippet_target}
${CMAKE_CURRENT_BINARY_DIR}/${compile_snippet_src})
if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO)
target_link_libraries(${compile_snippet_target} ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO})
endif()
GET_TARGET_PROPERTY(compile_snippet_executable
${compile_snippet_target} LOCATION)
ADD_CUSTOM_COMMAND(

View File

@@ -1,8 +1,5 @@
#!/bin/sh
# TODO : actually exit on exit, currently it only exit from the ()
# TODO : display error msg on stderr instead of stdout
# configuration
# You should call this script with USER set as you want, else some default
# will be used
@@ -11,16 +8,12 @@ USER=${USER:-'orzel'}
# step 1 : build
# todo if 'build is not there, create one:
#mkdir build
(cd build && cmake .. && make -j3 doc) || (echo "make failed"; exit 1)
(cd build && cmake .. && make -j3 doc) || { echo "make failed"; exit 1; }
#todo: n+1 where n = number of cpus
#step 2 : upload
BRANCH=`hg branch`
if [ $BRANCH == "default" ]
then
BRANCH='devel'
fi
# (the '/' at the end of path are very important, see rsync documentation)
rsync -az build/doc/html/ $USER@ssh.tuxfamily.org:eigen/eigen.tuxfamily.org-web/htdocs/dox-$BRANCH/ || (echo "upload failed"; exit 1)
rsync -az build/doc/html/ $USER@ssh.tuxfamily.org:eigen/eigen.tuxfamily.org-web/htdocs/dox-2.0/ || { echo "upload failed"; exit 1; }
echo "Uploaded successfully"

View File

@@ -1,3 +1,7 @@
option(EIGEN_DEFAULT_TO_ROW_MAJOR "Use row-major as default matrix storage order" OFF)
if(EIGEN_DEFAULT_TO_ROW_MAJOR)
add_definitions("-DEIGEN_DEFAULT_TO_ROW_MAJOR")
endif()
find_package(GSL)
if(GSL_FOUND AND GSL_VERSION_MINOR LESS 9)
@@ -93,12 +97,20 @@ else(CMAKE_COMPILER_IS_GNUCXX)
endif(CMAKE_COMPILER_IS_GNUCXX)
option(EIGEN_NO_ASSERTION_CHECKING "Disable checking of assertions" OFF)
if(EIGEN_NO_ASSERTION_CHECKING)
add_definitions("-DEIGEN_NO_ASSERTION_CHECKING=1")
endif()
# similar to set_target_properties but append the property instead of overwriting it
macro(ei_add_target_property target prop value)
get_target_property(previous ${target} ${prop})
set_target_properties(${target} PROPERTIES ${prop} "${previous} ${value}")
if(previous MATCHES "NOTFOUND")
set_target_properties(${target} PROPERTIES ${prop} "${value}")
else()
set_target_properties(${target} PROPERTIES ${prop} "${previous} ${value}")
endif()
endmacro(ei_add_target_property)
@@ -134,13 +146,9 @@ macro(ei_add_test testname)
option(EIGEN_DEBUG_ASSERTS "Enable debuging of assertions" OFF)
if(EIGEN_DEBUG_ASSERTS)
set_target_properties(${targetname} PROPERTIES COMPILE_DEFINITIONS "-DEIGEN_DEBUG_ASSERTS=1")
set_target_properties(${targetname} PROPERTIES COMPILE_DEFINITIONS "EIGEN_DEBUG_ASSERTS=1")
endif(EIGEN_DEBUG_ASSERTS)
else(NOT EIGEN_NO_ASSERTION_CHECKING)
set_target_properties(${targetname} PROPERTIES COMPILE_DEFINITIONS "-DEIGEN_NO_ASSERTION_CHECKING=1")
endif(NOT EIGEN_NO_ASSERTION_CHECKING)
if(${ARGC} GREATER 1)
@@ -153,6 +161,10 @@ macro(ei_add_test testname)
target_link_libraries(${targetname} Eigen2)
endif(TEST_LIB)
if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO)
target_link_libraries(${targetname} ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO})
endif()
target_link_libraries(${targetname} ${EXTERNAL_LIBS})
if(${ARGC} GREATER 2)
string(STRIP "${ARGV2}" ARGV2_stripped)
@@ -180,6 +192,7 @@ ei_add_test(meta)
ei_add_test(sizeof)
ei_add_test(dynalloc)
ei_add_test(nomalloc)
ei_add_test(first_aligned)
ei_add_test(mixingtypes)
ei_add_test(packetmath)
ei_add_test(unalignedassert)
@@ -200,7 +213,7 @@ ei_add_test(array)
ei_add_test(triangular)
ei_add_test(cholesky " " "${GSL_LIBRARIES}")
ei_add_test(lu ${EI_OFLAG})
ei_add_test(determinant)
ei_add_test(determinant ${EI_OFLAG})
ei_add_test(inverse)
ei_add_test(qr)
ei_add_test(eigensolver " " "${GSL_LIBRARIES}")
@@ -215,12 +228,17 @@ ei_add_test(newstdvector)
if(QT4_FOUND)
ei_add_test(qtvector " " "${QT_QTCORE_LIBRARY}")
endif(QT4_FOUND)
ei_add_test(sparse_vector)
ei_add_test(sparse_basic)
ei_add_test(sparse_solvers " " "${SPARSE_LIBS}")
ei_add_test(sparse_product)
if(NOT EIGEN_DEFAULT_TO_ROW_MAJOR)
ei_add_test(sparse_vector)
ei_add_test(sparse_basic)
ei_add_test(sparse_solvers " " "${SPARSE_LIBS}")
ei_add_test(sparse_product)
endif()
ei_add_test(swap)
ei_add_test(visitor)
ei_add_test(bug_132)
ei_add_test(prec_inverse_4x4 ${EI_OFLAG})
# print a summary of the different options
message("************************************************************")
@@ -265,11 +283,22 @@ else(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION)
message("Explicit vec: AUTO")
endif(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION)
if(EIGEN_DEFAULT_TO_ROW_MAJOR)
message("Default order: Row-major")
else()
message("Default order: Column-major")
endif()
string(TOLOWER "${CMAKE_CXX_COMPILER}" cmake_cxx_compiler_tolower)
if(cmake_cxx_compiler_tolower MATCHES "qcc")
set(CXX_IS_QCC "ON")
endif()
message("CXX: ${CMAKE_CXX_COMPILER}")
if(CMAKE_COMPILER_IS_GNUCXX)
if(CMAKE_COMPILER_IS_GNUCXX AND NOT CXX_IS_QCC)
execute_process(COMMAND ${CMAKE_CXX_COMPILER} --version COMMAND head -n 1 OUTPUT_VARIABLE EIGEN_CXX_VERSION_STRING OUTPUT_STRIP_TRAILING_WHITESPACE)
message("CXX_VERSION: ${EIGEN_CXX_VERSION_STRING}")
endif(CMAKE_COMPILER_IS_GNUCXX)
endif()
message("CXX_FLAGS: ${CMAKE_CXX_FLAGS}")
message("Sparse lib flags: ${SPARSE_LIBS}")

41
test/bug_132.cpp Normal file
View File

@@ -0,0 +1,41 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra. Eigen itself is part of the KDE project.
//
// Copyright (C) 2010 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// 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 <http://www.gnu.org/licenses/>.
#include "main.h"
void test_bug_132() {
enum { size = 100 };
MatrixXd A(size, size);
VectorXd b(size), c(size);
{
VectorXd y = A.transpose() * (b-c); // bug 132: infinite recursion in coeffRef
VectorXd z = (b-c).transpose() * A; // bug 132: infinite recursion in coeffRef
}
// the following ones weren't failing, but let's include them for completeness:
{
VectorXd y = A * (b-c);
VectorXd z = (b-c).transpose() * A.transpose();
}
}

View File

@@ -100,6 +100,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
VERIFY_IS_APPROX(symm * matX, matB);
}
#if 0 // cholesky is not rank-revealing anyway
// test isPositiveDefinite on non definite matrix
if (rows>4)
{
@@ -109,6 +110,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
LDLT<SquareMatrixType> cholnosqrt(symm);
VERIFY(!cholnosqrt.isPositiveDefinite());
}
#endif
}
void test_cholesky()
@@ -122,4 +124,8 @@ void test_cholesky()
CALL_SUBTEST( cholesky(MatrixXf(17,17)) );
CALL_SUBTEST( cholesky(MatrixXd(33,33)) );
}
MatrixXf m = MatrixXf::Zero(10,10);
VectorXf b = VectorXf::Zero(10);
VERIFY(!m.llt().isPositiveDefinite());
}

View File

@@ -35,7 +35,7 @@ void check_handmade_aligned_malloc()
for(int i = 1; i < 1000; i++)
{
char *p = (char*)ei_handmade_aligned_malloc(i);
VERIFY(size_t(p)%ALIGNMENT==0);
VERIFY(std::size_t(p)%ALIGNMENT==0);
// if the buffer is wrongly allocated this will give a bad write --> check with valgrind
for(int j = 0; j < i; j++) p[j]=0;
ei_handmade_aligned_free(p);
@@ -47,7 +47,7 @@ void check_aligned_malloc()
for(int i = 1; i < 1000; i++)
{
char *p = (char*)ei_aligned_malloc(i);
VERIFY(size_t(p)%ALIGNMENT==0);
VERIFY(std::size_t(p)%ALIGNMENT==0);
// if the buffer is wrongly allocated this will give a bad write --> check with valgrind
for(int j = 0; j < i; j++) p[j]=0;
ei_aligned_free(p);
@@ -59,7 +59,7 @@ void check_aligned_new()
for(int i = 1; i < 1000; i++)
{
float *p = ei_aligned_new<float>(i);
VERIFY(size_t(p)%ALIGNMENT==0);
VERIFY(std::size_t(p)%ALIGNMENT==0);
// if the buffer is wrongly allocated this will give a bad write --> check with valgrind
for(int j = 0; j < i; j++) p[j]=0;
ei_aligned_delete(p,i);
@@ -71,7 +71,7 @@ void check_aligned_stack_alloc()
for(int i = 1; i < 1000; i++)
{
float *p = ei_aligned_stack_new(float,i);
VERIFY(size_t(p)%ALIGNMENT==0);
VERIFY(std::size_t(p)%ALIGNMENT==0);
// if the buffer is wrongly allocated this will give a bad write --> check with valgrind
for(int j = 0; j < i; j++) p[j]=0;
ei_aligned_stack_delete(float,p,i);
@@ -98,7 +98,7 @@ class MyClassA
template<typename T> void check_dynaligned()
{
T* obj = new T;
VERIFY(size_t(obj)%ALIGNMENT==0);
VERIFY(std::size_t(obj)%ALIGNMENT==0);
delete obj;
}
@@ -121,15 +121,15 @@ void test_dynalloc()
// check static allocation, who knows ?
{
MyStruct foo0; VERIFY(size_t(foo0.avec.data())%ALIGNMENT==0);
MyClassA fooA; VERIFY(size_t(fooA.avec.data())%ALIGNMENT==0);
MyStruct foo0; VERIFY(std::size_t(foo0.avec.data())%ALIGNMENT==0);
MyClassA fooA; VERIFY(std::size_t(fooA.avec.data())%ALIGNMENT==0);
}
// dynamic allocation, single object
for (int i=0; i<g_repeat*100; ++i)
{
MyStruct *foo0 = new MyStruct(); VERIFY(size_t(foo0->avec.data())%ALIGNMENT==0);
MyClassA *fooA = new MyClassA(); VERIFY(size_t(fooA->avec.data())%ALIGNMENT==0);
MyStruct *foo0 = new MyStruct(); VERIFY(std::size_t(foo0->avec.data())%ALIGNMENT==0);
MyClassA *fooA = new MyClassA(); VERIFY(std::size_t(fooA->avec.data())%ALIGNMENT==0);
delete foo0;
delete fooA;
}
@@ -138,8 +138,8 @@ void test_dynalloc()
const int N = 10;
for (int i=0; i<g_repeat*100; ++i)
{
MyStruct *foo0 = new MyStruct[N]; VERIFY(size_t(foo0->avec.data())%ALIGNMENT==0);
MyClassA *fooA = new MyClassA[N]; VERIFY(size_t(fooA->avec.data())%ALIGNMENT==0);
MyStruct *foo0 = new MyStruct[N]; VERIFY(std::size_t(foo0->avec.data())%ALIGNMENT==0);
MyClassA *fooA = new MyClassA[N]; VERIFY(std::size_t(fooA->avec.data())%ALIGNMENT==0);
delete[] foo0;
delete[] fooA;
}

64
test/first_aligned.cpp Normal file
View File

@@ -0,0 +1,64 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// 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 <http://www.gnu.org/licenses/>.
#include "main.h"
template<typename Scalar>
void test_first_aligned_helper(Scalar *array, int size)
{
const int packet_size = sizeof(Scalar) * ei_packet_traits<Scalar>::size;
VERIFY(((std::size_t(array) + sizeof(Scalar) * ei_alignmentOffset(array, size)) % packet_size) == 0);
}
template<typename Scalar>
void test_none_aligned_helper(Scalar *array, int size)
{
VERIFY(ei_packet_traits<Scalar>::size == 1 || ei_alignmentOffset(array, size) == size);
}
struct some_non_vectorizable_type { float x; };
void test_first_aligned()
{
EIGEN_ALIGN_128 float array_float[100];
test_first_aligned_helper(array_float, 50);
test_first_aligned_helper(array_float+1, 50);
test_first_aligned_helper(array_float+2, 50);
test_first_aligned_helper(array_float+3, 50);
test_first_aligned_helper(array_float+4, 50);
test_first_aligned_helper(array_float+5, 50);
EIGEN_ALIGN_128 double array_double[100];
test_first_aligned_helper(array_double, 50);
test_first_aligned_helper(array_double+1, 50);
test_first_aligned_helper(array_double+2, 50);
double *array_double_plus_4_bytes = (double*)(std::size_t(array_double)+4);
test_none_aligned_helper(array_double_plus_4_bytes, 50);
test_none_aligned_helper(array_double_plus_4_bytes+1, 50);
some_non_vectorizable_type array_nonvec[100];
test_first_aligned_helper(array_nonvec, 100);
test_none_aligned_helper(array_nonvec, 100);
}

View File

@@ -43,11 +43,9 @@ template<typename MatrixType> void inverse(const MatrixType& m)
mzero = MatrixType::Zero(rows, cols),
identity = MatrixType::Identity(rows, rows);
if (ei_is_same_type<RealScalar,float>::ret)
while(ei_abs(m1.determinant()) < RealScalar(0.1) && rows <= 8)
{
// let's build a more stable to inverse matrix
MatrixType a = MatrixType::Random(rows,cols);
m1 += m1 * m1.adjoint() + a * a.adjoint();
m1 = MatrixType::Random(rows, cols);
}
m2 = m1.inverse();

View File

@@ -132,4 +132,10 @@ void test_lu()
CALL_SUBTEST( lu_invertible<MatrixXcf>() );
CALL_SUBTEST( lu_invertible<MatrixXcd>() );
}
MatrixXf m = MatrixXf::Zero(10,10);
VectorXf b = VectorXf::Zero(10);
VectorXf x = VectorXf::Random(10);
VERIFY(m.lu().solve(b,&x));
VERIFY(x.isZero());
}

View File

@@ -142,7 +142,7 @@ namespace Eigen
#define VERIFY(a) do { if (!(a)) { \
std::cerr << "Test " << g_test_stack.back() << " failed in "EI_PP_MAKE_STRING(__FILE__) << " (" << EI_PP_MAKE_STRING(__LINE__) << ")" \
<< std::endl << " " << EI_PP_MAKE_STRING(a) << std::endl << std::endl; \
exit(2); \
std::exit(2); \
} } while (0)
#define VERIFY_IS_APPROX(a, b) VERIFY(test_ei_isApprox(a, b))
@@ -257,7 +257,7 @@ int main(int argc, char *argv[])
std::cout << "Argument " << argv[i] << " conflicting with a former argument" << std::endl;
return 1;
}
repeat = atoi(argv[i]+1);
repeat = std::atoi(argv[i]+1);
has_set_repeat = true;
if(repeat <= 0)
{
@@ -272,7 +272,7 @@ int main(int argc, char *argv[])
std::cout << "Argument " << argv[i] << " conflicting with a former argument" << std::endl;
return 1;
}
seed = int(strtoul(argv[i]+1, 0, 10));
seed = int(std::strtoul(argv[i]+1, 0, 10));
has_set_seed = true;
bool ok = seed!=0;
if(!ok)
@@ -295,11 +295,11 @@ int main(int argc, char *argv[])
return 1;
}
if(!has_set_seed) seed = (unsigned int) time(NULL);
if(!has_set_seed) seed = (unsigned int) std::time(NULL);
if(!has_set_repeat) repeat = DEFAULT_REPEAT;
std::cout << "Initializing random number generator with seed " << seed << std::endl;
srand(seed);
std::srand(seed);
std::cout << "Repeating each test " << repeat << " times" << std::endl;
Eigen::g_repeat = repeat;

View File

@@ -1,7 +1,7 @@
// 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 <jacob.benoit.1@gmail.com>
// Copyright (C) 2007-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@@ -24,7 +24,7 @@
#include "main.h"
template<typename VectorType> void map_class(const VectorType& m)
template<typename VectorType> void map_class_vector(const VectorType& m)
{
typedef typename VectorType::Scalar Scalar;
@@ -34,7 +34,7 @@ template<typename VectorType> void map_class(const VectorType& m)
Scalar* array1 = ei_aligned_new<Scalar>(size);
Scalar* array2 = ei_aligned_new<Scalar>(size);
Scalar* array3 = new Scalar[size+1];
Scalar* array3unaligned = size_t(array3)%16 == 0 ? array3+1 : array3;
Scalar* array3unaligned = std::size_t(array3)%16 == 0 ? array3+1 : array3;
Map<VectorType, Aligned>(array1, size) = VectorType::Random(size);
Map<VectorType>(array2, size) = Map<VectorType>(array1, size);
@@ -50,6 +50,34 @@ template<typename VectorType> void map_class(const VectorType& m)
delete[] array3;
}
template<typename MatrixType> void map_class_matrix(const MatrixType& m)
{
typedef typename MatrixType::Scalar Scalar;
int rows = m.rows(), cols = m.cols(), size = rows*cols;
// test Map.h
Scalar* array1 = ei_aligned_new<Scalar>(size);
for(int i = 0; i < size; i++) array1[i] = Scalar(1);
Scalar* array2 = ei_aligned_new<Scalar>(size);
for(int i = 0; i < size; i++) array2[i] = Scalar(1);
Scalar* array3 = new Scalar[size+1];
for(int i = 0; i < size+1; i++) array3[i] = Scalar(1);
Scalar* array3unaligned = std::size_t(array3)%16 == 0 ? array3+1 : array3;
Map<MatrixType, Aligned>(array1, rows, cols) = MatrixType::Ones(rows,cols);
Map<MatrixType>(array2, rows, cols) = Map<MatrixType>(array1, rows, cols);
Map<MatrixType>(array3unaligned, rows, cols) = Map<MatrixType>(array1, rows, cols);
MatrixType ma1 = Map<MatrixType>(array1, rows, cols);
MatrixType ma2 = Map<MatrixType, Aligned>(array2, rows, cols);
VERIFY_IS_APPROX(ma1, ma2);
MatrixType ma3 = Map<MatrixType>(array3unaligned, rows, cols);
VERIFY_IS_APPROX(ma1, ma3);
ei_aligned_delete(array1, size);
ei_aligned_delete(array2, size);
delete[] array3;
}
template<typename VectorType> void map_static_methods(const VectorType& m)
{
typedef typename VectorType::Scalar Scalar;
@@ -60,7 +88,7 @@ template<typename VectorType> void map_static_methods(const VectorType& m)
Scalar* array1 = ei_aligned_new<Scalar>(size);
Scalar* array2 = ei_aligned_new<Scalar>(size);
Scalar* array3 = new Scalar[size+1];
Scalar* array3unaligned = size_t(array3)%16 == 0 ? array3+1 : array3;
Scalar* array3unaligned = std::size_t(array3)%16 == 0 ? array3+1 : array3;
VectorType::MapAligned(array1, size) = VectorType::Random(size);
VectorType::Map(array2, size) = VectorType::Map(array1, size);
@@ -80,11 +108,17 @@ template<typename VectorType> void map_static_methods(const VectorType& m)
void test_map()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( map_class(Matrix<float, 1, 1>()) );
CALL_SUBTEST( map_class(Vector4d()) );
CALL_SUBTEST( map_class(RowVector4f()) );
CALL_SUBTEST( map_class(VectorXcf(8)) );
CALL_SUBTEST( map_class(VectorXi(12)) );
CALL_SUBTEST( map_class_vector(Matrix<float, 1, 1>()) );
CALL_SUBTEST( map_class_vector(Vector4d()) );
CALL_SUBTEST( map_class_vector(RowVector4f()) );
CALL_SUBTEST( map_class_vector(VectorXcf(8)) );
CALL_SUBTEST( map_class_vector(VectorXi(12)) );
CALL_SUBTEST( map_class_matrix(Matrix<float, 1, 1>()) );
CALL_SUBTEST( map_class_matrix(Matrix4d()) );
CALL_SUBTEST( map_class_matrix(Matrix<float,3,5>()) );
CALL_SUBTEST( map_class_matrix(MatrixXcf(ei_random<int>(1,10),ei_random<int>(1,10))) );
CALL_SUBTEST( map_class_matrix(MatrixXi(ei_random<int>(1,10),ei_random<int>(1,10))) );
CALL_SUBTEST( map_static_methods(Matrix<double, 1, 1>()) );
CALL_SUBTEST( map_static_methods(Vector3f()) );

View File

@@ -50,7 +50,7 @@ void check_stdvector_matrix(const MatrixType& m)
VERIFY_IS_APPROX(v[21], y);
v.push_back(x);
VERIFY_IS_APPROX(v[22], x);
VERIFY((size_t)&(v[22]) == (size_t)&(v[21]) + sizeof(MatrixType));
VERIFY((std::size_t)&(v[22]) == (std::size_t)&(v[21]) + sizeof(MatrixType));
// do a lot of push_back such that the vector gets internally resized
// (with memory reallocation)
@@ -85,7 +85,7 @@ void check_stdvector_transform(const TransformType&)
VERIFY_IS_APPROX(v[21], y);
v.push_back(x);
VERIFY_IS_APPROX(v[22], x);
VERIFY((size_t)&(v[22]) == (size_t)&(v[21]) + sizeof(TransformType));
VERIFY((std::size_t)&(v[22]) == (std::size_t)&(v[21]) + sizeof(TransformType));
// do a lot of push_back such that the vector gets internally resized
// (with memory reallocation)
@@ -120,7 +120,7 @@ void check_stdvector_quaternion(const QuaternionType&)
VERIFY_IS_APPROX(v[21], y);
v.push_back(x);
VERIFY_IS_APPROX(v[22], x);
VERIFY((size_t)&(v[22]) == (size_t)&(v[21]) + sizeof(QuaternionType));
VERIFY((std::size_t)&(v[22]) == (std::size_t)&(v[21]) + sizeof(QuaternionType));
// do a lot of push_back such that the vector gets internally resized
// (with memory reallocation)

99
test/prec_inverse_4x4.cpp Normal file
View File

@@ -0,0 +1,99 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// 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 <http://www.gnu.org/licenses/>.
#include "main.h"
#include <Eigen/LU>
#include <algorithm>
template<typename T> std::string type_name() { return "other"; }
template<> std::string type_name<float>() { return "float"; }
template<> std::string type_name<double>() { return "double"; }
template<> std::string type_name<int>() { return "int"; }
template<> std::string type_name<std::complex<float> >() { return "complex<float>"; }
template<> std::string type_name<std::complex<double> >() { return "complex<double>"; }
template<> std::string type_name<std::complex<int> >() { return "complex<int>"; }
#define EIGEN_DEBUG_VAR(x) std::cerr << #x << " = " << x << std::endl;
template<typename T> inline typename NumTraits<T>::Real epsilon()
{
return std::numeric_limits<typename NumTraits<T>::Real>::epsilon();
}
template<typename MatrixType> void inverse_permutation_4x4()
{
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
Vector4i indices(0,1,2,3);
for(int i = 0; i < 24; ++i)
{
MatrixType m = MatrixType::Zero();
m(indices(0),0) = 1;
m(indices(1),1) = 1;
m(indices(2),2) = 1;
m(indices(3),3) = 1;
MatrixType inv = m.inverse();
double error = double( (m*inv-MatrixType::Identity()).norm() / epsilon<Scalar>() );
VERIFY(error == 0.0);
std::next_permutation(indices.data(),indices.data()+4);
}
}
template<typename MatrixType> void inverse_general_4x4(int repeat)
{
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
double error_sum = 0., error_max = 0.;
for(int i = 0; i < repeat; ++i)
{
MatrixType m;
RealScalar absdet;
do {
m = MatrixType::Random();
absdet = ei_abs(m.determinant());
} while(absdet < 10 * epsilon<Scalar>());
MatrixType inv = m.inverse();
double error = double( (m*inv-MatrixType::Identity()).norm() * absdet / epsilon<Scalar>() );
error_sum += error;
error_max = std::max(error_max, error);
}
std::cerr << "inverse_general_4x4, Scalar = " << type_name<Scalar>() << std::endl;
double error_avg = error_sum / repeat;
EIGEN_DEBUG_VAR(error_avg);
EIGEN_DEBUG_VAR(error_max);
VERIFY(error_avg < (NumTraits<Scalar>::IsComplex ? 8.0 : 1.0));
VERIFY(error_max < (NumTraits<Scalar>::IsComplex ? 64.0 : 20.0));
}
void test_prec_inverse_4x4()
{
CALL_SUBTEST((inverse_permutation_4x4<Matrix4f>()));
CALL_SUBTEST(( inverse_general_4x4<Matrix4f>(200000 * g_repeat) ));
CALL_SUBTEST((inverse_permutation_4x4<Matrix<double,4,4,RowMajor> >()));
CALL_SUBTEST(( inverse_general_4x4<Matrix<double,4,4,RowMajor> >(200000 * g_repeat) ));
CALL_SUBTEST((inverse_permutation_4x4<Matrix4cf>()));
CALL_SUBTEST((inverse_general_4x4<Matrix4cf>(50000 * g_repeat)));
}

View File

@@ -46,7 +46,7 @@ template<typename MatrixType> void product(const MatrixType& m)
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> RowSquareMatrixType;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> ColSquareMatrixType;
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime,
MatrixType::Flags&RowMajorBit> OtherMajorMatrixType;
MatrixType::Options^RowMajor> OtherMajorMatrixType;
int rows = m.rows();
int cols = m.cols();
@@ -77,6 +77,7 @@ template<typename MatrixType> void product(const MatrixType& m)
// begin testing Product.h: only associativity for now
// (we use Transpose.h but this doesn't count as a test for it)
VERIFY_IS_APPROX((m1*m1.transpose())*m2, m1*(m1.transpose()*m2));
m3 = m1;
m3 *= m1.transpose() * m2;
@@ -137,6 +138,7 @@ template<typename MatrixType> void product(const MatrixType& m)
res2 = square2;
res2 += (m1.transpose() * m2).lazy();
VERIFY_IS_APPROX(res2, square2 + m1.transpose() * m2);
if (NumTraits<Scalar>::HasFloatingPoint && std::min(rows,cols)>1)
{
VERIFY(areNotApprox(res2,square2 + m2.transpose() * m1));

View File

@@ -48,4 +48,11 @@ void test_product_large()
MatrixXf a = MatrixXf::Random(10,4), b = MatrixXf::Random(4,10), c = a;
VERIFY_IS_APPROX((a = a * b), (c * b).eval());
}
{
MatrixXf mat1(10,10); mat1.setRandom();
MatrixXf mat2(32,10); mat2.setRandom();
MatrixXf result = mat1.row(2)*mat2.transpose();
VERIFY_IS_APPROX(result, (mat1.row(2)*mat2.transpose()).eval());
}
}

View File

@@ -75,4 +75,11 @@ void test_qr()
mat << 1, 1, 1, 2, 2, 2, 1, 2, 3;
VERIFY(!mat.qr().isFullRank());
}
{
MatrixXf m = MatrixXf::Zero(10,10);
VectorXf b = VectorXf::Zero(10);
VectorXf x = VectorXf::Random(10);
VERIFY(m.qr().solve(b,&x));
VERIFY(x.isZero());
}
}

View File

@@ -68,12 +68,20 @@ template<typename Scalar> void sparse_solvers(int rows, int cols)
VERIFY_IS_APPROX(refMat2.template marked<LowerTriangular>().solveTriangular(vec2),
m2.template marked<LowerTriangular>().solveTriangular(vec3));
// lower - transpose
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
VERIFY_IS_APPROX(refMat2.template marked<LowerTriangular>().transpose().solveTriangular(vec2),
m2.template marked<LowerTriangular>().transpose().solveTriangular(vec3));
// upper
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
VERIFY_IS_APPROX(refMat2.template marked<UpperTriangular>().solveTriangular(vec2),
m2.template marked<UpperTriangular>().solveTriangular(vec3));
// TODO test row major
// upper - transpose
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
VERIFY_IS_APPROX(refMat2.template marked<UpperTriangular>().transpose().solveTriangular(vec2),
m2.template marked<UpperTriangular>().transpose().solveTriangular(vec3));
}
// test LLT

View File

@@ -49,7 +49,7 @@ void check_stdvector_matrix(const MatrixType& m)
VERIFY_IS_APPROX(v[21], y);
v.push_back(x);
VERIFY_IS_APPROX(v[22], x);
VERIFY((size_t)&(v[22]) == (size_t)&(v[21]) + sizeof(MatrixType));
VERIFY((std::size_t)&(v[22]) == (std::size_t)&(v[21]) + sizeof(MatrixType));
// do a lot of push_back such that the vector gets internally resized
// (with memory reallocation)
@@ -84,7 +84,7 @@ void check_stdvector_transform(const TransformType&)
VERIFY_IS_APPROX(v[21], y);
v.push_back(x);
VERIFY_IS_APPROX(v[22], x);
VERIFY((size_t)&(v[22]) == (size_t)&(v[21]) + sizeof(TransformType));
VERIFY((std::size_t)&(v[22]) == (std::size_t)&(v[21]) + sizeof(TransformType));
// do a lot of push_back such that the vector gets internally resized
// (with memory reallocation)
@@ -119,7 +119,7 @@ void check_stdvector_quaternion(const QuaternionType&)
VERIFY_IS_APPROX(v[21], y);
v.push_back(x);
VERIFY_IS_APPROX(v[22], x);
VERIFY((size_t)&(v[22]) == (size_t)&(v[21]) + sizeof(QuaternionType));
VERIFY((std::size_t)&(v[22]) == (std::size_t)&(v[21]) + sizeof(QuaternionType));
// do a lot of push_back such that the vector gets internally resized
// (with memory reallocation)

View File

@@ -95,5 +95,8 @@ void test_svd()
// complex are not implemented yet
// CALL_SUBTEST( svd(MatrixXcd(6,6)) );
// CALL_SUBTEST( svd(MatrixXcf(3,3)) );
SVD<MatrixXf> s;
MatrixXf m = MatrixXf::Random(10,1);
s.compute(m);
}
}

View File

@@ -39,7 +39,7 @@
# VERSION=opensuse-11.1
# WORK_DIR=/home/gael/Coding/eigen2/cdash
# # get the last version of the script
# wget http://bitbucket.org/eigen/eigen2/raw/tip/test/testsuite.cmake -o $WORK_DIR/testsuite.cmake
# wget http://bitbucket.org/eigen/eigen/raw/tip/test/testsuite.cmake -o $WORK_DIR/testsuite.cmake
# COMMON="ctest -S $WORK_DIR/testsuite.cmake,EIGEN_WORK_DIR=$WORK_DIR,EIGEN_SITE=$SITE,EIGEN_MODE=$1,EIGEN_BUILD_STRING=$OS_VERSION-$ARCH"
# $COMMON-gcc-3.4.6,EIGEN_CXX=g++-3.4
# $COMMON-gcc-4.0.1,EIGEN_CXX=g++-4.0.1
@@ -133,7 +133,7 @@ endif(NOT EIGEN_MODE)
## mandatory variables (the default should be ok in most cases):
SET (CTEST_CVS_COMMAND "hg")
SET (CTEST_CVS_CHECKOUT "${CTEST_CVS_COMMAND} clone -r 2.0 http://bitbucket.org/eigen/eigen2 \"${CTEST_SOURCE_DIRECTORY}\"")
SET (CTEST_CVS_CHECKOUT "${CTEST_CVS_COMMAND} clone -r 2.0 http://bitbucket.org/eigen/eigen \"${CTEST_SOURCE_DIRECTORY}\"")
# which ctest command to use for running the dashboard
SET (CTEST_COMMAND "${EIGEN_CMAKE_DIR}ctest -D ${EIGEN_MODE}")

View File

@@ -98,7 +98,7 @@ void check_unalignedassert_bad()
{
float buf[sizeof(T)+16];
float *unaligned = buf;
while((reinterpret_cast<size_t>(unaligned)&0xf)==0) ++unaligned; // make sure unaligned is really unaligned
while((reinterpret_cast<std::size_t>(unaligned)&0xf)==0) ++unaligned; // make sure unaligned is really unaligned
T *x = ::new(static_cast<void*>(unaligned)) T;
x->~T();
}

View File

@@ -44,12 +44,21 @@ void test_vectorization_logic()
#ifdef EIGEN_VECTORIZE
#ifdef EIGEN_DEFAULT_TO_ROW_MAJOR
VERIFY(test_assign(Vector4f(),Vector4f(),
LinearVectorization,CompleteUnrolling));
VERIFY(test_assign(Vector4f(),Vector4f()+Vector4f(),
LinearVectorization,CompleteUnrolling));
VERIFY(test_assign(Vector4f(),Vector4f().cwise() * Vector4f(),
LinearVectorization,CompleteUnrolling));
#else
VERIFY(test_assign(Vector4f(),Vector4f(),
InnerVectorization,CompleteUnrolling));
VERIFY(test_assign(Vector4f(),Vector4f()+Vector4f(),
InnerVectorization,CompleteUnrolling));
VERIFY(test_assign(Vector4f(),Vector4f().cwise() * Vector4f(),
InnerVectorization,CompleteUnrolling));
#endif
VERIFY(test_assign(Matrix4f(),Matrix4f(),
InnerVectorization,CompleteUnrolling));
@@ -76,6 +85,8 @@ void test_vectorization_logic()
VERIFY(test_assign(MatrixXf(10,10),MatrixXf(20,20).block(10,10,2,3),
SliceVectorization,NoUnrolling));
VERIFY(test_assign(VectorXf(10),VectorXf(10)+VectorXf(10),
LinearVectorization,NoUnrolling));
VERIFY(test_sum(VectorXf(10),
LinearVectorization,NoUnrolling));
@@ -92,8 +103,10 @@ void test_vectorization_logic()
VERIFY(test_sum(Matrix<float,16,16>().block<4,4>(1,2),
NoVectorization,CompleteUnrolling));
#ifndef EIGEN_DEFAULT_TO_ROW_MAJOR
VERIFY(test_sum(Matrix<float,16,16>().block<8,1>(1,2),
LinearVectorization,CompleteUnrolling));
#endif
VERIFY(test_sum(Matrix<double,7,3>(),
NoVectorization,CompleteUnrolling));