diff --git a/CMakeLists.txt b/CMakeLists.txt
index eaee5d5e2..abae3b23c 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -120,7 +120,7 @@ endmacro(ei_add_cxx_compiler_flag)
if(NOT MSVC)
# We assume that other compilers are partly compatible with GNUCC
- # clang outputs some warnings for unknwon flags that are not caught by check_cxx_compiler_flag
+ # clang outputs some warnings for unknown flags that are not caught by check_cxx_compiler_flag
# adding -Werror turns such warnings into errors
check_cxx_compiler_flag("-Werror" COMPILER_SUPPORT_WERROR)
if(COMPILER_SUPPORT_WERROR)
@@ -141,8 +141,11 @@ if(NOT MSVC)
ei_add_cxx_compiler_flag("-Wwrite-strings")
ei_add_cxx_compiler_flag("-Wformat-security")
ei_add_cxx_compiler_flag("-Wshorten-64-to-32")
+ ei_add_cxx_compiler_flag("-Wlogical-op")
ei_add_cxx_compiler_flag("-Wenum-conversion")
ei_add_cxx_compiler_flag("-Wc++11-extensions")
+ ei_add_cxx_compiler_flag("-Wdouble-promotion")
+# ei_add_cxx_compiler_flag("-Wconversion")
# -Wshadow is insanely too strict with gcc, hopefully it will become usable with gcc 6
# if(NOT CMAKE_COMPILER_IS_GNUCXX OR (CMAKE_CXX_COMPILER_VERSION VERSION_GREATER "5.0.0"))
@@ -158,7 +161,7 @@ if(NOT MSVC)
ei_add_cxx_compiler_flag("-fno-common")
ei_add_cxx_compiler_flag("-fstrict-aliasing")
ei_add_cxx_compiler_flag("-wd981") # disable ICC's "operands are evaluated in unspecified order" remark
- ei_add_cxx_compiler_flag("-wd2304") # disbale ICC's "warning #2304: non-explicit constructor with single argument may cause implicit type conversion" produced by -Wnon-virtual-dtor
+ ei_add_cxx_compiler_flag("-wd2304") # disable ICC's "warning #2304: non-explicit constructor with single argument may cause implicit type conversion" produced by -Wnon-virtual-dtor
# The -ansi flag must be added last, otherwise it is also used as a linker flag by check_cxx_compiler_flag making it fails
@@ -221,6 +224,12 @@ if(NOT MSVC)
message(STATUS "Enabling FMA in tests/examples")
endif()
+ option(EIGEN_TEST_F16C "Enable/Disable F16C in tests/examples" OFF)
+ if(EIGEN_TEST_F16C)
+ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -mf16c")
+ message(STATUS "Enabling F16C in tests/examples")
+ endif()
+
option(EIGEN_TEST_ALTIVEC "Enable/Disable AltiVec in tests/examples" OFF)
if(EIGEN_TEST_ALTIVEC)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -maltivec -mabi=altivec")
@@ -240,7 +249,7 @@ if(NOT MSVC)
else()
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -mfpu=neon")
endif()
- set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -mfloat-abi=softfp")
+ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -mfloat-abi=hard")
message(STATUS "Enabling NEON in tests/examples")
endif()
@@ -250,7 +259,11 @@ if(NOT MSVC)
message(STATUS "Enabling NEON in tests/examples")
endif()
-
+ option(EIGEN_TEST_ZVECTOR "Enable/Disable S390X(zEC13) ZVECTOR in tests/examples" OFF)
+ if(EIGEN_TEST_ZVECTOR)
+ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=z13 -mzvector")
+ message(STATUS "Enabling S390X(zEC13) ZVECTOR in tests/examples")
+ endif()
check_cxx_compiler_flag("-fopenmp" COMPILER_SUPPORT_OPENMP)
if(COMPILER_SUPPORT_OPENMP)
@@ -336,6 +349,8 @@ endif()
option(EIGEN_TEST_CXX11 "Enable testing with C++11 and C++11 features (e.g. Tensor module)." OFF)
+set(EIGEN_CUDA_COMPUTE_ARCH 30 CACHE STRING "The CUDA compute architecture level to target when compiling CUDA code")
+
include_directories(${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR})
# Backward compatibility support for EIGEN_INCLUDE_INSTALL_DIR
diff --git a/Eigen/CholmodSupport b/Eigen/CholmodSupport
index 83e2c1da4..bed8924d3 100644
--- a/Eigen/CholmodSupport
+++ b/Eigen/CholmodSupport
@@ -19,7 +19,7 @@ extern "C" {
/** \ingroup Support_modules
* \defgroup CholmodSupport_Module CholmodSupport module
*
- * This module provides an interface to the Cholmod library which is part of the suitesparse package.
+ * This module provides an interface to the Cholmod library which is part of the suitesparse package.
* It provides the two following main factorization classes:
* - class CholmodSupernodalLLT: a supernodal LLT Cholesky factorization.
* - class CholmodDecomposiiton: a general L(D)LT Cholesky factorization with automatic or explicit runtime selection of the underlying factorization method (supernodal or simplicial).
diff --git a/Eigen/Core b/Eigen/Core
index 63602f4c3..6c79bcfae 100644
--- a/Eigen/Core
+++ b/Eigen/Core
@@ -33,16 +33,19 @@
#ifdef EIGEN_EXCEPTIONS
#undef EIGEN_EXCEPTIONS
#endif
-
+
// All functions callable from CUDA code must be qualified with __device__
#define EIGEN_DEVICE_FUNC __host__ __device__
-
+
#else
#define EIGEN_DEVICE_FUNC
-
+
#endif
-#if defined(__CUDA_ARCH__)
+// When compiling CUDA device code with NVCC, pull in math functions from the
+// global namespace. In host mode, and when device doee with clang, use the
+// std versions.
+#if defined(__CUDA_ARCH__) && defined(__NVCC__)
#define EIGEN_USING_STD_MATH(FUNC) using ::FUNC;
#else
#define EIGEN_USING_STD_MATH(FUNC) using std::FUNC;
@@ -153,6 +156,7 @@
#if EIGEN_COMP_ICC >= 1110
#include
#else
+ #include
#include
#include
#ifdef EIGEN_VECTORIZE_SSE3
@@ -194,12 +198,29 @@
#define EIGEN_VECTORIZE
#define EIGEN_VECTORIZE_NEON
#include
+ #elif (defined __s390x__ && defined __VEC__)
+ #define EIGEN_VECTORIZE
+ #define EIGEN_VECTORIZE_ZVECTOR
+ #include
#endif
#endif
+#if defined(__F16C__) && !defined(EIGEN_COMP_CLANG)
+ // We can use the optimized fp16 to float and float to fp16 conversion routines
+ #define EIGEN_HAS_FP16_C
+#endif
+
#if defined __CUDACC__
#define EIGEN_VECTORIZE_CUDA
#include
+ #if defined __CUDACC_VER__ && __CUDACC_VER__ >= 70500
+ #define EIGEN_HAS_CUDA_FP16
+ #endif
+#endif
+
+#if defined EIGEN_HAS_CUDA_FP16
+ #include
+ #include
#endif
#if (defined _OPENMP) && (!defined EIGEN_DONT_PARALLELIZE)
@@ -263,6 +284,8 @@ inline static const char *SimdInstructionSetsInUse(void) {
return "VSX";
#elif defined(EIGEN_VECTORIZE_NEON)
return "ARM NEON";
+#elif defined(EIGEN_VECTORIZE_ZVECTOR)
+ return "S390X ZVECTOR";
#else
return "None";
#endif
@@ -278,7 +301,7 @@ inline static const char *SimdInstructionSetsInUse(void) {
// 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;
-// gcc 4.6.0 wants std:: for ptrdiff_t
+// gcc 4.6.0 wants std:: for ptrdiff_t
using std::ptrdiff_t;
/** \defgroup Core_Module Core module
@@ -325,8 +348,17 @@ using std::ptrdiff_t;
#include "src/Core/arch/NEON/PacketMath.h"
#include "src/Core/arch/NEON/MathFunctions.h"
#include "src/Core/arch/NEON/Complex.h"
+#elif defined EIGEN_VECTORIZE_ZVECTOR
+ #include "src/Core/arch/ZVector/PacketMath.h"
+ #include "src/Core/arch/ZVector/MathFunctions.h"
+ #include "src/Core/arch/ZVector/Complex.h"
#endif
+// Half float support
+#include "src/Core/arch/CUDA/Half.h"
+#include "src/Core/arch/CUDA/PacketMathHalf.h"
+#include "src/Core/arch/CUDA/TypeCasting.h"
+
#if defined EIGEN_VECTORIZE_CUDA
#include "src/Core/arch/CUDA/PacketMath.h"
#include "src/Core/arch/CUDA/MathFunctions.h"
@@ -414,6 +446,7 @@ using std::ptrdiff_t;
#include "src/Core/products/TriangularSolverVector.h"
#include "src/Core/BandMatrix.h"
#include "src/Core/CoreIterators.h"
+#include "src/Core/ConditionEstimator.h"
#include "src/Core/BooleanRedux.h"
#include "src/Core/Select.h"
@@ -424,14 +457,14 @@ using std::ptrdiff_t;
#include "src/Core/ArrayWrapper.h"
#ifdef EIGEN_USE_BLAS
-#include "src/Core/products/GeneralMatrixMatrix_MKL.h"
-#include "src/Core/products/GeneralMatrixVector_MKL.h"
-#include "src/Core/products/GeneralMatrixMatrixTriangular_MKL.h"
-#include "src/Core/products/SelfadjointMatrixMatrix_MKL.h"
-#include "src/Core/products/SelfadjointMatrixVector_MKL.h"
-#include "src/Core/products/TriangularMatrixMatrix_MKL.h"
-#include "src/Core/products/TriangularMatrixVector_MKL.h"
-#include "src/Core/products/TriangularSolverMatrix_MKL.h"
+#include "src/Core/products/GeneralMatrixMatrix_BLAS.h"
+#include "src/Core/products/GeneralMatrixVector_BLAS.h"
+#include "src/Core/products/GeneralMatrixMatrixTriangular_BLAS.h"
+#include "src/Core/products/SelfadjointMatrixMatrix_BLAS.h"
+#include "src/Core/products/SelfadjointMatrixVector_BLAS.h"
+#include "src/Core/products/TriangularMatrixMatrix_BLAS.h"
+#include "src/Core/products/TriangularMatrixVector_BLAS.h"
+#include "src/Core/products/TriangularSolverMatrix_BLAS.h"
#endif // EIGEN_USE_BLAS
#ifdef EIGEN_USE_MKL_VML
diff --git a/Eigen/Geometry b/Eigen/Geometry
index 06b736e3f..716d52952 100644
--- a/Eigen/Geometry
+++ b/Eigen/Geometry
@@ -17,16 +17,16 @@
#include
/** \defgroup Geometry_Module Geometry module
- *
- *
*
* This module provides support for:
* - fixed-size homogeneous transformations
* - translation, scaling, 2D and 3D rotations
- * - quaternions
- * - \ref MatrixBase::cross() "cross product"
- * - \ref MatrixBase::unitOrthogonal() "orthognal vector generation"
- * - some linear components: parametrized-lines and hyperplanes
+ * - \link Quaternion quaternions \endlink
+ * - cross products (\ref MatrixBase::cross, \ref MatrixBase::cross3)
+ * - orthognal vector generation (\ref MatrixBase::unitOrthogonal)
+ * - some linear components: \link ParametrizedLine parametrized-lines \endlink and \link Hyperplane hyperplanes \endlink
+ * - \link AlignedBox axis aligned bounding boxes \endlink
+ * - \link umeyama least-square transformation fitting \endlink
*
* \code
* #include
diff --git a/Eigen/PaStiXSupport b/Eigen/PaStiXSupport
index 3411dface..de3a63b4d 100644
--- a/Eigen/PaStiXSupport
+++ b/Eigen/PaStiXSupport
@@ -12,7 +12,6 @@
#include "src/Core/util/DisableStupidWarnings.h"
-#include
extern "C" {
#include
#include
diff --git a/Eigen/QR b/Eigen/QR
index f74f365f1..25c781cc1 100644
--- a/Eigen/QR
+++ b/Eigen/QR
@@ -34,6 +34,7 @@
#include "src/QR/HouseholderQR.h"
#include "src/QR/FullPivHouseholderQR.h"
#include "src/QR/ColPivHouseholderQR.h"
+#include "src/QR/CompleteOrthogonalDecomposition.h"
#ifdef EIGEN_USE_LAPACKE
#include "src/QR/HouseholderQR_MKL.h"
#include "src/QR/ColPivHouseholderQR_MKL.h"
diff --git a/Eigen/SPQRSupport b/Eigen/SPQRSupport
index f9489dcd8..f70390c17 100644
--- a/Eigen/SPQRSupport
+++ b/Eigen/SPQRSupport
@@ -17,7 +17,7 @@
/** \ingroup Support_modules
* \defgroup SPQRSupport_Module SuiteSparseQR module
*
- * This module provides an interface to the SPQR library, which is part of the suitesparse package.
+ * This module provides an interface to the SPQR library, which is part of the suitesparse package.
*
* \code
* #include
diff --git a/Eigen/UmfPackSupport b/Eigen/UmfPackSupport
index 4a9f46a1e..00eec8087 100644
--- a/Eigen/UmfPackSupport
+++ b/Eigen/UmfPackSupport
@@ -19,7 +19,7 @@ extern "C" {
/** \ingroup Support_modules
* \defgroup UmfPackSupport_Module UmfPackSupport module
*
- * This module provides an interface to the UmfPack library which is part of the suitesparse package.
+ * This module provides an interface to the UmfPack library which is part of the suitesparse package.
* It provides the following factorization class:
* - class UmfPackLU: a multifrontal sequential LU factorization.
*
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
index 6fcae01f7..538aff956 100644
--- a/Eigen/src/Cholesky/LDLT.h
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -13,7 +13,7 @@
#ifndef EIGEN_LDLT_H
#define EIGEN_LDLT_H
-namespace Eigen {
+namespace Eigen {
namespace internal {
template struct LDLT_Traits;
@@ -28,8 +28,8 @@ namespace internal {
*
* \brief Robust Cholesky decomposition of a matrix with pivoting
*
- * \param MatrixType the type of the matrix of which to compute the LDL^T Cholesky decomposition
- * \param UpLo the triangular part that will be used for the decompositon: Lower (default) or Upper.
+ * \tparam _MatrixType the type of the matrix of which to compute the LDL^T Cholesky decomposition
+ * \tparam _UpLo the triangular part that will be used for the decompositon: Lower (default) or Upper.
* The other triangular part won't be read.
*
* Perform a robust Cholesky decomposition of a positive semidefinite or negative semidefinite
@@ -73,11 +73,11 @@ template class LDLT
* The default constructor is useful in cases in which the user intends to
* perform decompositions via LDLT::compute(const MatrixType&).
*/
- LDLT()
- : m_matrix(),
- m_transpositions(),
+ LDLT()
+ : m_matrix(),
+ m_transpositions(),
m_sign(internal::ZeroSign),
- m_isInitialized(false)
+ m_isInitialized(false)
{}
/** \brief Default Constructor with memory preallocation
@@ -168,7 +168,7 @@ template class LDLT
* \note_about_checking_solutions
*
* More precisely, this method solves \f$ A x = b \f$ using the decomposition \f$ A = P^T L D L^* P \f$
- * by solving the systems \f$ P^T y_1 = b \f$, \f$ L y_2 = y_1 \f$, \f$ D y_3 = y_2 \f$,
+ * by solving the systems \f$ P^T y_1 = b \f$, \f$ L y_2 = y_1 \f$, \f$ D y_3 = y_2 \f$,
* \f$ L^* y_4 = y_3 \f$ and \f$ P x = y_4 \f$ in succession. If the matrix \f$ A \f$ is singular, then
* \f$ D \f$ will also be singular (all the other matrices are invertible). In that case, the
* least-square solution of \f$ D y_3 = y_2 \f$ is computed. This does not mean that this function
@@ -192,6 +192,15 @@ template class LDLT
template
LDLT& compute(const EigenBase& matrix);
+ /** \returns an estimate of the reciprocal condition number of the matrix of
+ * which \c *this is the LDLT decomposition.
+ */
+ RealScalar rcond() const
+ {
+ eigen_assert(m_isInitialized && "LDLT is not initialized.");
+ return internal::rcond_estimate_helper(m_l1_norm, *this);
+ }
+
template
LDLT& rankUpdate(const MatrixBase& w, const RealScalar& alpha=1);
@@ -207,6 +216,13 @@ template class LDLT
MatrixType reconstructedMatrix() const;
+ /** \returns the adjoint of \c *this, that is, a const reference to the decomposition itself as the underlying matrix is self-adjoint.
+ *
+ * This method is provided for compatibility with other matrix decompositions, thus enabling generic code such as:
+ * \code x = decomposition.adjoint().solve(b) \endcode
+ */
+ const LDLT& adjoint() const { return *this; };
+
inline Index rows() const { return m_matrix.rows(); }
inline Index cols() const { return m_matrix.cols(); }
@@ -220,7 +236,7 @@ template class LDLT
eigen_assert(m_isInitialized && "LDLT is not initialized.");
return Success;
}
-
+
#ifndef EIGEN_PARSED_BY_DOXYGEN
template
EIGEN_DEVICE_FUNC
@@ -228,7 +244,7 @@ template class LDLT
#endif
protected:
-
+
static void check_template_parameters()
{
EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
@@ -241,6 +257,7 @@ template class LDLT
* is not stored), and the diagonal entries correspond to D.
*/
MatrixType m_matrix;
+ RealScalar m_l1_norm;
TranspositionType m_transpositions;
TmpMatrixType m_temporary;
internal::SignMatrix m_sign;
@@ -266,8 +283,8 @@ template<> struct ldlt_inplace
if (size <= 1)
{
transpositions.setIdentity();
- if (numext::real(mat.coeff(0,0)) > 0) sign = PositiveSemiDef;
- else if (numext::real(mat.coeff(0,0)) < 0) sign = NegativeSemiDef;
+ if (numext::real(mat.coeff(0,0)) > static_cast(0) ) sign = PositiveSemiDef;
+ else if (numext::real(mat.coeff(0,0)) < static_cast(0)) sign = NegativeSemiDef;
else sign = ZeroSign;
return true;
}
@@ -314,7 +331,7 @@ template<> struct ldlt_inplace
if(rs>0)
A21.noalias() -= A20 * temp.head(k);
}
-
+
// In some previous versions of Eigen (e.g., 3.2.1), the scaling was omitted if the pivot
// was smaller than the cutoff value. However, since LDLT is not rank-revealing
// we should only make sure that we do not introduce INF or NaN values.
@@ -324,12 +341,12 @@ template<> struct ldlt_inplace
A21 /= realAkk;
if (sign == PositiveSemiDef) {
- if (realAkk < 0) sign = Indefinite;
+ if (realAkk < static_cast(0)) sign = Indefinite;
} else if (sign == NegativeSemiDef) {
- if (realAkk > 0) sign = Indefinite;
+ if (realAkk > static_cast(0)) sign = Indefinite;
} else if (sign == ZeroSign) {
- if (realAkk > 0) sign = PositiveSemiDef;
- else if (realAkk < 0) sign = NegativeSemiDef;
+ if (realAkk > static_cast(0)) sign = PositiveSemiDef;
+ else if (realAkk < static_cast(0)) sign = NegativeSemiDef;
}
}
@@ -433,12 +450,25 @@ template
LDLT& LDLT::compute(const EigenBase& a)
{
check_template_parameters();
-
+
eigen_assert(a.rows()==a.cols());
const Index size = a.rows();
m_matrix = a.derived();
+ // Compute matrix L1 norm = max abs column sum.
+ m_l1_norm = RealScalar(0);
+ // TODO move this code to SelfAdjointView
+ for (Index col = 0; col < size; ++col) {
+ RealScalar abs_col_sum;
+ if (_UpLo == Lower)
+ abs_col_sum = m_matrix.col(col).tail(size - col).template lpNorm<1>() + m_matrix.row(col).head(col).template lpNorm<1>();
+ else
+ abs_col_sum = m_matrix.col(col).head(col).template lpNorm<1>() + m_matrix.row(col).tail(size - col).template lpNorm<1>();
+ if (abs_col_sum > m_l1_norm)
+ m_l1_norm = abs_col_sum;
+ }
+
m_transpositions.resize(size);
m_isInitialized = false;
m_temporary.resize(size);
@@ -466,7 +496,7 @@ LDLT& LDLT::rankUpdate(const MatrixBase::_solve_impl(const RhsType &rhs, DstType &dst) cons
// diagonal element is not well justified and leads to numerical issues in some cases.
// Moreover, Lapack's xSYTRS routines use 0 for the tolerance.
RealScalar tolerance = RealScalar(1) / NumTraits::highest();
-
+
for (Index i = 0; i < vecD.size(); ++i)
{
if(abs(vecD(i)) > tolerance)
diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
index 1f0091f3c..19578b216 100644
--- a/Eigen/src/Cholesky/LLT.h
+++ b/Eigen/src/Cholesky/LLT.h
@@ -10,7 +10,7 @@
#ifndef EIGEN_LLT_H
#define EIGEN_LLT_H
-namespace Eigen {
+namespace Eigen {
namespace internal{
template struct LLT_Traits;
@@ -22,8 +22,8 @@ template struct LLT_Traits;
*
* \brief Standard Cholesky decomposition (LL^T) of a matrix and associated features
*
- * \param MatrixType the type of the matrix of which we are computing the LL^T Cholesky decomposition
- * \param UpLo the triangular part that will be used for the decompositon: Lower (default) or Upper.
+ * \tparam _MatrixType the type of the matrix of which we are computing the LL^T Cholesky decomposition
+ * \tparam _UpLo the triangular part that will be used for the decompositon: Lower (default) or Upper.
* The other triangular part won't be read.
*
* This class performs a LL^T Cholesky decomposition of a symmetric, positive definite
@@ -40,7 +40,7 @@ template struct LLT_Traits;
*
* Example: \include LLT_example.cpp
* Output: \verbinclude LLT_example.out
- *
+ *
* \sa MatrixBase::llt(), SelfAdjointView::llt(), class LDLT
*/
/* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH)
@@ -135,6 +135,16 @@ template class LLT
template
LLT& compute(const EigenBase& matrix);
+ /** \returns an estimate of the reciprocal condition number of the matrix of
+ * which \c *this is the Cholesky decomposition.
+ */
+ RealScalar rcond() const
+ {
+ eigen_assert(m_isInitialized && "LLT is not initialized.");
+ eigen_assert(m_info == Success && "LLT failed because matrix appears to be negative");
+ return internal::rcond_estimate_helper(m_l1_norm, *this);
+ }
+
/** \returns the LLT decomposition matrix
*
* TODO: document the storage layout
@@ -159,12 +169,19 @@ template class LLT
return m_info;
}
+ /** \returns the adjoint of \c *this, that is, a const reference to the decomposition itself as the underlying matrix is self-adjoint.
+ *
+ * This method is provided for compatibility with other matrix decompositions, thus enabling generic code such as:
+ * \code x = decomposition.adjoint().solve(b) \endcode
+ */
+ const LLT& adjoint() const { return *this; };
+
inline Index rows() const { return m_matrix.rows(); }
inline Index cols() const { return m_matrix.cols(); }
template
LLT rankUpdate(const VectorType& vec, const RealScalar& sigma = 1);
-
+
#ifndef EIGEN_PARSED_BY_DOXYGEN
template
EIGEN_DEVICE_FUNC
@@ -172,17 +189,18 @@ template class LLT
#endif
protected:
-
+
static void check_template_parameters()
{
EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
}
-
+
/** \internal
* Used to compute and store L
* The strict upper part is not used and even not initialized.
*/
MatrixType m_matrix;
+ RealScalar m_l1_norm;
bool m_isInitialized;
ComputationInfo m_info;
};
@@ -268,7 +286,7 @@ template struct llt_inplace
static Index unblocked(MatrixType& mat)
{
using std::sqrt;
-
+
eigen_assert(mat.rows()==mat.cols());
const Index size = mat.rows();
for(Index k = 0; k < size; ++k)
@@ -328,7 +346,7 @@ template struct llt_inplace
return Eigen::internal::llt_rank_update_lower(mat, vec, sigma);
}
};
-
+
template struct llt_inplace
{
typedef typename NumTraits::Real RealScalar;
@@ -387,12 +405,25 @@ template
LLT& LLT::compute(const EigenBase& a)
{
check_template_parameters();
-
+
eigen_assert(a.rows()==a.cols());
const Index size = a.rows();
m_matrix.resize(size, size);
m_matrix = a.derived();
+ // Compute matrix L1 norm = max abs column sum.
+ m_l1_norm = RealScalar(0);
+ // TODO move this code to SelfAdjointView
+ for (Index col = 0; col < size; ++col) {
+ RealScalar abs_col_sum;
+ if (_UpLo == Lower)
+ abs_col_sum = m_matrix.col(col).tail(size - col).template lpNorm<1>() + m_matrix.row(col).head(col).template lpNorm<1>();
+ else
+ abs_col_sum = m_matrix.col(col).head(col).template lpNorm<1>() + m_matrix.row(col).tail(size - col).template lpNorm<1>();
+ if (abs_col_sum > m_l1_norm)
+ m_l1_norm = abs_col_sum;
+ }
+
m_isInitialized = true;
bool ok = Traits::inplace_decomposition(m_matrix);
m_info = ok ? Success : NumericalIssue;
@@ -419,7 +450,7 @@ LLT<_MatrixType,_UpLo> LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, c
return *this;
}
-
+
#ifndef EIGEN_PARSED_BY_DOXYGEN
template
template
@@ -431,15 +462,12 @@ void LLT<_MatrixType,_UpLo>::_solve_impl(const RhsType &rhs, DstType &dst) const
#endif
/** \internal use x = llt_object.solve(x);
- *
+ *
* This is the \em in-place version of solve().
*
* \param bAndX represents both the right-hand side matrix b and result x.
*
- * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD.
- *
- * This version avoids a copy when the right hand side matrix b is not
- * needed anymore.
+ * This version avoids a copy when the right hand side matrix b is not needed anymore.
*
* \sa LLT::solve(), MatrixBase::llt()
*/
@@ -486,7 +514,7 @@ SelfAdjointView::llt() const
return LLT(m_matrix);
}
#endif // __CUDACC__
-
+
} // end namespace Eigen
#endif // EIGEN_LLT_H
diff --git a/Eigen/src/Cholesky/LLT_MKL.h b/Eigen/src/Cholesky/LLT_MKL.h
index 0d42cb5bc..f5be3b240 100644
--- a/Eigen/src/Cholesky/LLT_MKL.h
+++ b/Eigen/src/Cholesky/LLT_MKL.h
@@ -53,11 +53,11 @@ template<> struct mkl_llt \
EIGTYPE* a; \
eigen_assert(m.rows()==m.cols()); \
/* Set up parameters for ?potrf */ \
- size = m.rows(); \
+ size = convert_index(m.rows()); \
StorageOrder = MatrixType::Flags&RowMajorBit?RowMajor:ColMajor; \
matrix_order = StorageOrder==RowMajor ? LAPACK_ROW_MAJOR : LAPACK_COL_MAJOR; \
a = &(m.coeffRef(0,0)); \
- lda = m.outerStride(); \
+ lda = convert_index(m.outerStride()); \
\
info = LAPACKE_##MKLPREFIX##potrf( matrix_order, uplo, size, (MKLTYPE*)a, lda ); \
info = (info==0) ? -1 : info>0 ? info-1 : size; \
diff --git a/Eigen/src/CholmodSupport/CholmodSupport.h b/Eigen/src/CholmodSupport/CholmodSupport.h
index 06421d5ed..b8020a92c 100644
--- a/Eigen/src/CholmodSupport/CholmodSupport.h
+++ b/Eigen/src/CholmodSupport/CholmodSupport.h
@@ -78,7 +78,7 @@ cholmod_sparse viewAsCholmod(SparseMatrix<_Scalar,_Options,_StorageIndex>& mat)
{
res.itype = CHOLMOD_INT;
}
- else if (internal::is_same<_StorageIndex,SuiteSparse_long>::value)
+ else if (internal::is_same<_StorageIndex,long>::value)
{
res.itype = CHOLMOD_LONG;
}
@@ -178,14 +178,14 @@ class CholmodBase : public SparseSolverBase
public:
CholmodBase()
- : m_cholmodFactor(0), m_info(Success)
+ : m_cholmodFactor(0), m_info(Success), m_factorizationIsOk(false), m_analysisIsOk(false)
{
m_shiftOffset[0] = m_shiftOffset[1] = RealScalar(0.0);
cholmod_start(&m_cholmod);
}
explicit CholmodBase(const MatrixType& matrix)
- : m_cholmodFactor(0), m_info(Success)
+ : m_cholmodFactor(0), m_info(Success), m_factorizationIsOk(false), m_analysisIsOk(false)
{
m_shiftOffset[0] = m_shiftOffset[1] = RealScalar(0.0);
cholmod_start(&m_cholmod);
@@ -273,9 +273,10 @@ class CholmodBase : public SparseSolverBase
const Index size = m_cholmodFactor->n;
EIGEN_UNUSED_VARIABLE(size);
eigen_assert(size==b.rows());
+
+ // Cholmod needs column-major stoarge without inner-stride, which corresponds to the default behavior of Ref.
+ Ref > b_ref(b.derived());
- // note: cd stands for Cholmod Dense
- Rhs& b_ref(b.const_cast_derived());
cholmod_dense b_cd = viewAsCholmod(b_ref);
cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod);
if(!x_cd)
@@ -327,6 +328,57 @@ class CholmodBase : public SparseSolverBase
return derived();
}
+ /** \returns the determinant of the underlying matrix from the current factorization */
+ Scalar determinant() const
+ {
+ using std::exp;
+ return exp(logDeterminant());
+ }
+
+ /** \returns the log determinant of the underlying matrix from the current factorization */
+ Scalar logDeterminant() const
+ {
+ using std::log;
+ using numext::real;
+ eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
+
+ RealScalar logDet = 0;
+ Scalar *x = static_cast(m_cholmodFactor->x);
+ if (m_cholmodFactor->is_super)
+ {
+ // Supernodal factorization stored as a packed list of dense column-major blocs,
+ // as described by the following structure:
+
+ // super[k] == index of the first column of the j-th super node
+ StorageIndex *super = static_cast(m_cholmodFactor->super);
+ // pi[k] == offset to the description of row indices
+ StorageIndex *pi = static_cast(m_cholmodFactor->pi);
+ // px[k] == offset to the respective dense block
+ StorageIndex *px = static_cast(m_cholmodFactor->px);
+
+ Index nb_super_nodes = m_cholmodFactor->nsuper;
+ for (Index k=0; k < nb_super_nodes; ++k)
+ {
+ StorageIndex ncols = super[k + 1] - super[k];
+ StorageIndex nrows = pi[k + 1] - pi[k];
+
+ Map, 0, InnerStride<> > sk(x + px[k], ncols, InnerStride<>(nrows+1));
+ logDet += sk.real().log().sum();
+ }
+ }
+ else
+ {
+ // Simplicial factorization stored as standard CSC matrix.
+ StorageIndex *p = static_cast(m_cholmodFactor->p);
+ Index size = m_cholmodFactor->n;
+ for (Index k=0; kis_ll)
+ logDet *= 2.0;
+ return logDet;
+ };
+
template
void dumpMemory(Stream& /*s*/)
{}
@@ -358,7 +410,7 @@ class CholmodBase : public SparseSolverBase
*
* This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed.
*
- * \sa \ref TutorialSparseDirectSolvers, class CholmodSupernodalLLT, class SimplicialLLT
+ * \sa \ref TutorialSparseSolverConcept, class CholmodSupernodalLLT, class SimplicialLLT
*/
template
class CholmodSimplicialLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT<_MatrixType, _UpLo> >
@@ -407,7 +459,7 @@ class CholmodSimplicialLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimpl
*
* This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed.
*
- * \sa \ref TutorialSparseDirectSolvers, class CholmodSupernodalLLT, class SimplicialLDLT
+ * \sa \ref TutorialSparseSolverConcept, class CholmodSupernodalLLT, class SimplicialLDLT
*/
template
class CholmodSimplicialLDLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT<_MatrixType, _UpLo> >
@@ -454,7 +506,7 @@ class CholmodSimplicialLDLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimp
*
* This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed.
*
- * \sa \ref TutorialSparseDirectSolvers
+ * \sa \ref TutorialSparseSolverConcept
*/
template
class CholmodSupernodalLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT<_MatrixType, _UpLo> >
@@ -503,7 +555,7 @@ class CholmodSupernodalLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSuper
*
* This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed.
*
- * \sa \ref TutorialSparseDirectSolvers
+ * \sa \ref TutorialSparseSolverConcept
*/
template
class CholmodDecomposition : public CholmodBase<_MatrixType, _UpLo, CholmodDecomposition<_MatrixType, _UpLo> >
diff --git a/Eigen/src/Core/Array.h b/Eigen/src/Core/Array.h
index e38eda72c..c0af4aa9d 100644
--- a/Eigen/src/Core/Array.h
+++ b/Eigen/src/Core/Array.h
@@ -12,7 +12,16 @@
namespace Eigen {
-/** \class Array
+namespace internal {
+template
+struct traits > : traits >
+{
+ typedef ArrayXpr XprKind;
+ typedef ArrayBase > XprBase;
+};
+}
+
+/** \class Array
* \ingroup Core_Module
*
* \brief General-purpose arrays with easy API for coefficient-wise operations
@@ -26,21 +35,12 @@ namespace Eigen {
*
* See documentation of class Matrix for detailed information on the template parameters
* storage layout.
- *
+ *
* This class can be extended with the help of the plugin mechanism described on the page
* \ref TopicCustomizingEigen by defining the preprocessor symbol \c EIGEN_ARRAY_PLUGIN.
*
- * \sa \ref TutorialArrayClass, \ref TopicClassHierarchy
+ * \sa \blank \ref TutorialArrayClass, \ref TopicClassHierarchy
*/
-namespace internal {
-template
-struct traits > : traits >
-{
- typedef ArrayXpr XprKind;
- typedef ArrayBase > XprBase;
-};
-}
-
template
class Array
: public PlainObjectBase >
@@ -147,7 +147,7 @@ class Array
}
#endif
-#ifdef EIGEN_HAVE_RVALUE_REFERENCES
+#if EIGEN_HAS_RVALUE_REFERENCES
EIGEN_DEVICE_FUNC
Array(Array&& other)
: Base(std::move(other))
diff --git a/Eigen/src/Core/ArrayBase.h b/Eigen/src/Core/ArrayBase.h
index b4c24a27a..57333af1a 100644
--- a/Eigen/src/Core/ArrayBase.h
+++ b/Eigen/src/Core/ArrayBase.h
@@ -103,7 +103,7 @@ template class ArrayBase
/** Special case of the template operator=, in order to prevent the compiler
* from generating a default operator= (issue hit with g++ 4.1)
*/
- EIGEN_DEVICE_FUNC
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
Derived& operator=(const ArrayBase& other)
{
internal::call_assignment(derived(), other.derived());
@@ -112,28 +112,28 @@ template class ArrayBase
/** Set all the entries to \a value.
* \sa DenseBase::setConstant(), DenseBase::fill() */
- EIGEN_DEVICE_FUNC
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
Derived& operator=(const Scalar &value)
{ Base::setConstant(value); return derived(); }
- EIGEN_DEVICE_FUNC
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
Derived& operator+=(const Scalar& scalar);
- EIGEN_DEVICE_FUNC
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
Derived& operator-=(const Scalar& scalar);
template
- EIGEN_DEVICE_FUNC
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
Derived& operator+=(const ArrayBase& other);
template
- EIGEN_DEVICE_FUNC
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
Derived& operator-=(const ArrayBase& other);
template
- EIGEN_DEVICE_FUNC
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
Derived& operator*=(const ArrayBase& other);
template
- EIGEN_DEVICE_FUNC
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
Derived& operator/=(const ArrayBase& other);
public:
@@ -217,7 +217,7 @@ template
EIGEN_STRONG_INLINE Derived &
ArrayBase::operator/=(const ArrayBase& other)
{
- call_assignment(derived(), other.derived(), internal::div_assign_op());
+ call_assignment(derived(), other.derived(), internal::div_assign_op());
return derived();
}
diff --git a/Eigen/src/Core/ArrayWrapper.h b/Eigen/src/Core/ArrayWrapper.h
index 4e484f290..6013d4d85 100644
--- a/Eigen/src/Core/ArrayWrapper.h
+++ b/Eigen/src/Core/ArrayWrapper.h
@@ -52,7 +52,7 @@ class ArrayWrapper : public ArrayBase >
const Scalar
>::type ScalarWithConstIfNotLvalue;
- typedef typename internal::ref_selector::type NestedExpressionType;
+ typedef typename internal::ref_selector::non_const_type NestedExpressionType;
EIGEN_DEVICE_FUNC
explicit EIGEN_STRONG_INLINE ArrayWrapper(ExpressionType& matrix) : m_expression(matrix) {}
@@ -67,7 +67,7 @@ class ArrayWrapper : public ArrayBase >
inline Index innerStride() const { return m_expression.innerStride(); }
EIGEN_DEVICE_FUNC
- inline ScalarWithConstIfNotLvalue* data() { return m_expression.const_cast_derived().data(); }
+ inline ScalarWithConstIfNotLvalue* data() { return m_expression.data(); }
EIGEN_DEVICE_FUNC
inline const Scalar* data() const { return m_expression.data(); }
@@ -80,13 +80,13 @@ class ArrayWrapper : public ArrayBase >
EIGEN_DEVICE_FUNC
inline Scalar& coeffRef(Index rowId, Index colId)
{
- return m_expression.const_cast_derived().coeffRef(rowId, colId);
+ return m_expression.coeffRef(rowId, colId);
}
EIGEN_DEVICE_FUNC
inline const Scalar& coeffRef(Index rowId, Index colId) const
{
- return m_expression.const_cast_derived().coeffRef(rowId, colId);
+ return m_expression.coeffRef(rowId, colId);
}
EIGEN_DEVICE_FUNC
@@ -98,13 +98,13 @@ class ArrayWrapper : public ArrayBase >
EIGEN_DEVICE_FUNC
inline Scalar& coeffRef(Index index)
{
- return m_expression.const_cast_derived().coeffRef(index);
+ return m_expression.coeffRef(index);
}
EIGEN_DEVICE_FUNC
inline const Scalar& coeffRef(Index index) const
{
- return m_expression.const_cast_derived().coeffRef(index);
+ return m_expression.coeffRef(index);
}
template
@@ -116,7 +116,7 @@ class ArrayWrapper : public ArrayBase >
template
inline void writePacket(Index rowId, Index colId, const PacketScalar& val)
{
- m_expression.const_cast_derived().template writePacket(rowId, colId, val);
+ m_expression.template writePacket(rowId, colId, val);
}
template
@@ -128,7 +128,7 @@ class ArrayWrapper : public ArrayBase >
template
inline void writePacket(Index index, const PacketScalar& val)
{
- m_expression.const_cast_derived().template writePacket(index, val);
+ m_expression.template writePacket(index, val);
}
template
@@ -145,11 +145,11 @@ class ArrayWrapper : public ArrayBase >
/** Forwards the resizing request to the nested expression
* \sa DenseBase::resize(Index) */
EIGEN_DEVICE_FUNC
- void resize(Index newSize) { m_expression.const_cast_derived().resize(newSize); }
+ void resize(Index newSize) { m_expression.resize(newSize); }
/** Forwards the resizing request to the nested expression
* \sa DenseBase::resize(Index,Index)*/
EIGEN_DEVICE_FUNC
- void resize(Index rows, Index cols) { m_expression.const_cast_derived().resize(rows,cols); }
+ void resize(Index rows, Index cols) { m_expression.resize(rows,cols); }
protected:
NestedExpressionType m_expression;
@@ -195,7 +195,7 @@ class MatrixWrapper : public MatrixBase >
const Scalar
>::type ScalarWithConstIfNotLvalue;
- typedef typename internal::ref_selector::type NestedExpressionType;
+ typedef typename internal::ref_selector::non_const_type NestedExpressionType;
EIGEN_DEVICE_FUNC
explicit inline MatrixWrapper(ExpressionType& matrix) : m_expression(matrix) {}
@@ -210,7 +210,7 @@ class MatrixWrapper : public MatrixBase >
inline Index innerStride() const { return m_expression.innerStride(); }
EIGEN_DEVICE_FUNC
- inline ScalarWithConstIfNotLvalue* data() { return m_expression.const_cast_derived().data(); }
+ inline ScalarWithConstIfNotLvalue* data() { return m_expression.data(); }
EIGEN_DEVICE_FUNC
inline const Scalar* data() const { return m_expression.data(); }
@@ -223,7 +223,7 @@ class MatrixWrapper : public MatrixBase >
EIGEN_DEVICE_FUNC
inline Scalar& coeffRef(Index rowId, Index colId)
{
- return m_expression.const_cast_derived().coeffRef(rowId, colId);
+ return m_expression.coeffRef(rowId, colId);
}
EIGEN_DEVICE_FUNC
@@ -241,13 +241,13 @@ class MatrixWrapper : public MatrixBase >
EIGEN_DEVICE_FUNC
inline Scalar& coeffRef(Index index)
{
- return m_expression.const_cast_derived().coeffRef(index);
+ return m_expression.coeffRef(index);
}
EIGEN_DEVICE_FUNC
inline const Scalar& coeffRef(Index index) const
{
- return m_expression.const_cast_derived().coeffRef(index);
+ return m_expression.coeffRef(index);
}
template
@@ -259,7 +259,7 @@ class MatrixWrapper : public MatrixBase >
template
inline void writePacket(Index rowId, Index colId, const PacketScalar& val)
{
- m_expression.const_cast_derived().template writePacket(rowId, colId, val);
+ m_expression.template writePacket(rowId, colId, val);
}
template
@@ -271,7 +271,7 @@ class MatrixWrapper : public MatrixBase >
template
inline void writePacket(Index index, const PacketScalar& val)
{
- m_expression.const_cast_derived().template writePacket(index, val);
+ m_expression.template writePacket(index, val);
}
EIGEN_DEVICE_FUNC
@@ -284,11 +284,11 @@ class MatrixWrapper : public MatrixBase >
/** Forwards the resizing request to the nested expression
* \sa DenseBase::resize(Index) */
EIGEN_DEVICE_FUNC
- void resize(Index newSize) { m_expression.const_cast_derived().resize(newSize); }
+ void resize(Index newSize) { m_expression.resize(newSize); }
/** Forwards the resizing request to the nested expression
* \sa DenseBase::resize(Index,Index)*/
EIGEN_DEVICE_FUNC
- void resize(Index rows, Index cols) { m_expression.const_cast_derived().resize(rows,cols); }
+ void resize(Index rows, Index cols) { m_expression.resize(rows,cols); }
protected:
NestedExpressionType m_expression;
diff --git a/Eigen/src/Core/AssignEvaluator.h b/Eigen/src/Core/AssignEvaluator.h
old mode 100755
new mode 100644
index 9dfffbcc4..4b914ac0c
--- a/Eigen/src/Core/AssignEvaluator.h
+++ b/Eigen/src/Core/AssignEvaluator.h
@@ -29,13 +29,10 @@ struct copy_using_evaluator_traits
{
typedef typename DstEvaluator::XprType Dst;
typedef typename Dst::Scalar DstScalar;
- // TODO distinguish between linear traversal and inner-traversals
- typedef typename find_best_packet::type PacketType;
enum {
DstFlags = DstEvaluator::Flags,
- SrcFlags = SrcEvaluator::Flags,
- RequiredAlignment = unpacket_traits::alignment
+ SrcFlags = SrcEvaluator::Flags
};
public:
@@ -55,28 +52,43 @@ private:
: int(DstFlags)&RowMajorBit ? int(Dst::MaxColsAtCompileTime)
: int(Dst::MaxRowsAtCompileTime),
OuterStride = int(outer_stride_at_compile_time::ret),
- MaxSizeAtCompileTime = Dst::SizeAtCompileTime,
- PacketSize = unpacket_traits::size
+ MaxSizeAtCompileTime = Dst::SizeAtCompileTime
};
+ // TODO distinguish between linear traversal and inner-traversals
+ typedef typename find_best_packet::type LinearPacketType;
+ typedef typename find_best_packet::type InnerPacketType;
+
+ enum {
+ LinearPacketSize = unpacket_traits::size,
+ InnerPacketSize = unpacket_traits::size
+ };
+
+public:
+ enum {
+ LinearRequiredAlignment = unpacket_traits::alignment,
+ InnerRequiredAlignment = unpacket_traits::alignment
+ };
+
+private:
enum {
DstIsRowMajor = DstFlags&RowMajorBit,
SrcIsRowMajor = SrcFlags&RowMajorBit,
StorageOrdersAgree = (int(DstIsRowMajor) == int(SrcIsRowMajor)),
- MightVectorize = StorageOrdersAgree
+ MightVectorize = bool(StorageOrdersAgree)
&& (int(DstFlags) & int(SrcFlags) & ActualPacketAccessBit)
- && (functor_traits::PacketAccess),
+ && bool(functor_traits::PacketAccess),
MayInnerVectorize = MightVectorize
- && int(InnerSize)!=Dynamic && int(InnerSize)%int(PacketSize)==0
- && int(OuterStride)!=Dynamic && int(OuterStride)%int(PacketSize)==0
- && int(JointAlignment)>=int(RequiredAlignment),
- MayLinearize = StorageOrdersAgree && (int(DstFlags) & int(SrcFlags) & LinearAccessBit),
- MayLinearVectorize = MightVectorize && MayLinearize && DstHasDirectAccess
- && ((int(DstAlignment)>=int(RequiredAlignment)) || MaxSizeAtCompileTime == Dynamic),
+ && int(InnerSize)!=Dynamic && int(InnerSize)%int(InnerPacketSize)==0
+ && int(OuterStride)!=Dynamic && int(OuterStride)%int(InnerPacketSize)==0
+ && (EIGEN_UNALIGNED_VECTORIZE || int(JointAlignment)>=int(InnerRequiredAlignment)),
+ MayLinearize = bool(StorageOrdersAgree) && (int(DstFlags) & int(SrcFlags) & LinearAccessBit),
+ MayLinearVectorize = bool(MightVectorize) && MayLinearize && DstHasDirectAccess
+ && (EIGEN_UNALIGNED_VECTORIZE || (int(DstAlignment)>=int(LinearRequiredAlignment)) || MaxSizeAtCompileTime == Dynamic),
/* If the destination isn't aligned, we have to do runtime checks and we don't unroll,
so it's only good for large enough sizes. */
- MaySliceVectorize = MightVectorize && DstHasDirectAccess
- && (int(InnerMaxSize)==Dynamic || int(InnerMaxSize)>=3*PacketSize)
+ MaySliceVectorize = bool(MightVectorize) && bool(DstHasDirectAccess)
+ && (int(InnerMaxSize)==Dynamic || int(InnerMaxSize)>=3*InnerPacketSize)
/* slice vectorization can be slow, so we only want it if the slices are big, which is
indicated by InnerMaxSize rather than InnerSize, think of the case of a dynamic block
in a fixed-size matrix */
@@ -84,7 +96,8 @@ private:
public:
enum {
- Traversal = int(MayInnerVectorize) ? int(InnerVectorizedTraversal)
+ Traversal = int(MayLinearVectorize) && (LinearPacketSize>InnerPacketSize) ? int(LinearVectorizedTraversal)
+ : int(MayInnerVectorize) ? int(InnerVectorizedTraversal)
: int(MayLinearVectorize) ? int(LinearVectorizedTraversal)
: int(MaySliceVectorize) ? int(SliceVectorizedTraversal)
: int(MayLinearize) ? int(LinearTraversal)
@@ -94,9 +107,14 @@ public:
|| int(Traversal) == SliceVectorizedTraversal
};
+ typedef typename conditional::type PacketType;
+
private:
enum {
- UnrollingLimit = EIGEN_UNROLLING_LIMIT * (Vectorized ? int(PacketSize) : 1),
+ ActualPacketSize = int(Traversal)==LinearVectorizedTraversal ? LinearPacketSize
+ : Vectorized ? InnerPacketSize
+ : 1,
+ UnrollingLimit = EIGEN_UNROLLING_LIMIT * ActualPacketSize,
MayUnrollCompletely = int(Dst::SizeAtCompileTime) != Dynamic
&& int(Dst::SizeAtCompileTime) * int(SrcEvaluator::CoeffReadCost) <= int(UnrollingLimit),
MayUnrollInner = int(InnerSize) != Dynamic
@@ -112,8 +130,9 @@ public:
: int(NoUnrolling)
)
: int(Traversal) == int(LinearVectorizedTraversal)
- ? ( bool(MayUnrollCompletely) && (int(DstAlignment)>=int(RequiredAlignment)) ? int(CompleteUnrolling)
- : int(NoUnrolling) )
+ ? ( bool(MayUnrollCompletely) && ( EIGEN_UNALIGNED_VECTORIZE || (int(DstAlignment)>=int(LinearRequiredAlignment)))
+ ? int(CompleteUnrolling)
+ : int(NoUnrolling) )
: int(Traversal) == int(LinearTraversal)
? ( bool(MayUnrollCompletely) ? int(CompleteUnrolling)
: int(NoUnrolling) )
@@ -131,11 +150,14 @@ public:
std::cerr.unsetf(std::ios::hex);
EIGEN_DEBUG_VAR(DstAlignment)
EIGEN_DEBUG_VAR(SrcAlignment)
- EIGEN_DEBUG_VAR(RequiredAlignment)
+ EIGEN_DEBUG_VAR(LinearRequiredAlignment)
+ EIGEN_DEBUG_VAR(InnerRequiredAlignment)
EIGEN_DEBUG_VAR(JointAlignment)
EIGEN_DEBUG_VAR(InnerSize)
EIGEN_DEBUG_VAR(InnerMaxSize)
- EIGEN_DEBUG_VAR(PacketSize)
+ EIGEN_DEBUG_VAR(LinearPacketSize)
+ EIGEN_DEBUG_VAR(InnerPacketSize)
+ EIGEN_DEBUG_VAR(ActualPacketSize)
EIGEN_DEBUG_VAR(StorageOrdersAgree)
EIGEN_DEBUG_VAR(MightVectorize)
EIGEN_DEBUG_VAR(MayLinearize)
@@ -236,12 +258,13 @@ struct copy_using_evaluator_innervec_CompleteUnrolling
enum {
outer = Index / DstXprType::InnerSizeAtCompileTime,
inner = Index % DstXprType::InnerSizeAtCompileTime,
- JointAlignment = Kernel::AssignmentTraits::JointAlignment
+ SrcAlignment = Kernel::AssignmentTraits::SrcAlignment,
+ DstAlignment = Kernel::AssignmentTraits::DstAlignment
};
EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE void run(Kernel &kernel)
{
- kernel.template assignPacketByOuterInner(outer, inner);
+ kernel.template assignPacketByOuterInner(outer, inner);
enum { NextIndex = Index + unpacket_traits::size };
copy_using_evaluator_innervec_CompleteUnrolling::run(kernel);
}
@@ -257,9 +280,13 @@ template
struct copy_using_evaluator_innervec_InnerUnrolling
{
typedef typename Kernel::PacketType PacketType;
+ enum {
+ SrcAlignment = Kernel::AssignmentTraits::SrcAlignment,
+ DstAlignment = Kernel::AssignmentTraits::DstAlignment
+ };
EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE void run(Kernel &kernel, Index outer)
{
- kernel.template assignPacketByOuterInner(outer, Index_);
+ kernel.template assignPacketByOuterInner(outer, Index_);
enum { NextIndex = Index_ + unpacket_traits::size };
copy_using_evaluator_innervec_InnerUnrolling::run(kernel, outer);
}
@@ -370,7 +397,7 @@ struct dense_assignment_loop
typedef typename Kernel::Scalar Scalar;
typedef typename Kernel::PacketType PacketType;
enum {
- requestedAlignment = Kernel::AssignmentTraits::RequiredAlignment,
+ requestedAlignment = Kernel::AssignmentTraits::LinearRequiredAlignment,
packetSize = unpacket_traits::size,
dstIsAligned = int(Kernel::AssignmentTraits::DstAlignment)>=int(requestedAlignment),
dstAlignment = packet_traits::AlignedOnScalar ? int(requestedAlignment)
@@ -413,6 +440,10 @@ template
struct dense_assignment_loop
{
typedef typename Kernel::PacketType PacketType;
+ enum {
+ SrcAlignment = Kernel::AssignmentTraits::SrcAlignment,
+ DstAlignment = Kernel::AssignmentTraits::DstAlignment
+ };
EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE void run(Kernel &kernel)
{
const Index innerSize = kernel.innerSize();
@@ -420,7 +451,7 @@ struct dense_assignment_loop
const Index packetSize = unpacket_traits::size;
for(Index outer = 0; outer < outerSize; ++outer)
for(Index inner = 0; inner < innerSize; inner+=packetSize)
- kernel.template assignPacketByOuterInner(outer, inner);
+ kernel.template assignPacketByOuterInner(outer, inner);
}
};
@@ -484,14 +515,14 @@ struct dense_assignment_loop
typedef typename Kernel::PacketType PacketType;
enum {
packetSize = unpacket_traits::size,
- requestedAlignment = int(Kernel::AssignmentTraits::RequiredAlignment),
+ requestedAlignment = int(Kernel::AssignmentTraits::InnerRequiredAlignment),
alignable = packet_traits::AlignedOnScalar || int(Kernel::AssignmentTraits::DstAlignment)>=sizeof(Scalar),
dstIsAligned = int(Kernel::AssignmentTraits::DstAlignment)>=int(requestedAlignment),
dstAlignment = alignable ? int(requestedAlignment)
: int(Kernel::AssignmentTraits::DstAlignment)
};
const Scalar *dst_ptr = &kernel.dstEvaluator().coeffRef(0,0);
- if((!bool(dstIsAligned)) && (size_t(dst_ptr) % sizeof(Scalar))>0)
+ if((!bool(dstIsAligned)) && (UIntPtr(dst_ptr) % sizeof(Scalar))>0)
{
// the pointer is not aligend-on scalar, so alignment is not possible
return dense_assignment_loop::run(kernel);
@@ -637,7 +668,7 @@ protected:
***************************************************************************/
template
-EIGEN_DEVICE_FUNC void call_dense_assignment_loop(const DstXprType& dst, const SrcXprType& src, const Functor &func)
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void call_dense_assignment_loop(const DstXprType& dst, const SrcXprType& src, const Functor &func)
{
eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols());
@@ -654,7 +685,7 @@ EIGEN_DEVICE_FUNC void call_dense_assignment_loop(const DstXprType& dst, const S
}
template
-EIGEN_DEVICE_FUNC void call_dense_assignment_loop(const DstXprType& dst, const SrcXprType& src)
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void call_dense_assignment_loop(const DstXprType& dst, const SrcXprType& src)
{
call_dense_assignment_loop(dst, src, internal::assign_op());
}
@@ -682,47 +713,53 @@ template< typename DstXprType, typename SrcXprType, typename Functor,
struct Assignment;
-// The only purpose of this call_assignment() function is to deal with noalias() / AssumeAliasing and automatic transposition.
-// Indeed, I (Gael) think that this concept of AssumeAliasing was a mistake, and it makes thing quite complicated.
-// So this intermediate function removes everything related to AssumeAliasing such that Assignment
+// The only purpose of this call_assignment() function is to deal with noalias() / "assume-aliasing" and automatic transposition.
+// Indeed, I (Gael) think that this concept of "assume-aliasing" was a mistake, and it makes thing quite complicated.
+// So this intermediate function removes everything related to "assume-aliasing" such that Assignment
// does not has to bother about these annoying details.
template
-EIGEN_DEVICE_FUNC void call_assignment(Dst& dst, const Src& src)
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
+void call_assignment(Dst& dst, const Src& src)
{
call_assignment(dst, src, internal::assign_op());
}
template
-EIGEN_DEVICE_FUNC void call_assignment(const Dst& dst, const Src& src)
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
+void call_assignment(const Dst& dst, const Src& src)
{
call_assignment(dst, src, internal::assign_op());
}
-// Deal with AssumeAliasing
+// Deal with "assume-aliasing"
template
-EIGEN_DEVICE_FUNC void call_assignment(Dst& dst, const Src& src, const Func& func, typename enable_if::AssumeAliasing==1, void*>::type = 0)
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
+void call_assignment(Dst& dst, const Src& src, const Func& func, typename enable_if< evaluator_assume_aliasing::value, void*>::type = 0)
{
typename plain_matrix_type::type tmp(src);
call_assignment_no_alias(dst, tmp, func);
}
template
-EIGEN_DEVICE_FUNC void call_assignment(Dst& dst, const Src& src, const Func& func, typename enable_if::AssumeAliasing==0, void*>::type = 0)
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
+void call_assignment(Dst& dst, const Src& src, const Func& func, typename enable_if::value, void*>::type = 0)
{
call_assignment_no_alias(dst, src, func);
}
-// by-pass AssumeAliasing
+// by-pass "assume-aliasing"
// When there is no aliasing, we require that 'dst' has been properly resized
template class StorageBase, typename Src, typename Func>
-EIGEN_DEVICE_FUNC void call_assignment(NoAlias& dst, const Src& src, const Func& func)
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
+void call_assignment(NoAlias& dst, const Src& src, const Func& func)
{
call_assignment_no_alias(dst.expression(), src, func);
}
template
-EIGEN_DEVICE_FUNC void call_assignment_no_alias(Dst& dst, const Src& src, const Func& func)
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
+void call_assignment_no_alias(Dst& dst, const Src& src, const Func& func)
{
enum {
NeedToTranspose = ( (int(Dst::RowsAtCompileTime) == 1 && int(Src::ColsAtCompileTime) == 1)
@@ -747,13 +784,15 @@ EIGEN_DEVICE_FUNC void call_assignment_no_alias(Dst& dst, const Src& src, const
Assignment::run(actualDst, src, func);
}
template
-EIGEN_DEVICE_FUNC void call_assignment_no_alias(Dst& dst, const Src& src)
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
+void call_assignment_no_alias(Dst& dst, const Src& src)
{
call_assignment_no_alias(dst, src, internal::assign_op());
}
template
-EIGEN_DEVICE_FUNC void call_assignment_no_alias_no_transpose(Dst& dst, const Src& src, const Func& func)
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
+void call_assignment_no_alias_no_transpose(Dst& dst, const Src& src, const Func& func)
{
Index dstRows = src.rows();
Index dstCols = src.cols();
@@ -767,7 +806,8 @@ EIGEN_DEVICE_FUNC void call_assignment_no_alias_no_transpose(Dst& dst, const Src
Assignment::run(dst, src, func);
}
template
-EIGEN_DEVICE_FUNC void call_assignment_no_alias_no_transpose(Dst& dst, const Src& src)
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
+void call_assignment_no_alias_no_transpose(Dst& dst, const Src& src)
{
call_assignment_no_alias_no_transpose(dst, src, internal::assign_op());
}
@@ -779,7 +819,8 @@ template void check_for_aliasing(const Dst &dst, con
template< typename DstXprType, typename SrcXprType, typename Functor, typename Scalar>
struct Assignment
{
- EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
{
eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols());
@@ -796,7 +837,8 @@ struct Assignment
template< typename DstXprType, typename SrcXprType, typename Functor, typename Scalar>
struct Assignment
{
- EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &/*func*/)
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &/*func*/)
{
eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols());
src.evalTo(dst);
diff --git a/Eigen/src/Core/Assign_MKL.h b/Eigen/src/Core/Assign_MKL.h
old mode 100755
new mode 100644
diff --git a/Eigen/src/Core/BandMatrix.h b/Eigen/src/Core/BandMatrix.h
index 87c124fdf..4978c9140 100644
--- a/Eigen/src/Core/BandMatrix.h
+++ b/Eigen/src/Core/BandMatrix.h
@@ -161,15 +161,15 @@ class BandMatrixBase : public EigenBase
*
* \brief Represents a rectangular matrix with a banded storage
*
- * \param _Scalar Numeric type, i.e. float, double, int
- * \param Rows Number of rows, or \b Dynamic
- * \param Cols Number of columns, or \b Dynamic
- * \param Supers Number of super diagonal
- * \param Subs Number of sub diagonal
- * \param _Options A combination of either \b #RowMajor or \b #ColMajor, and of \b #SelfAdjoint
- * The former controls \ref TopicStorageOrders "storage order", and defaults to
- * column-major. The latter controls whether the matrix represents a selfadjoint
- * matrix in which case either Supers of Subs have to be null.
+ * \tparam _Scalar Numeric type, i.e. float, double, int
+ * \tparam _Rows Number of rows, or \b Dynamic
+ * \tparam _Cols Number of columns, or \b Dynamic
+ * \tparam _Supers Number of super diagonal
+ * \tparam _Subs Number of sub diagonal
+ * \tparam _Options A combination of either \b #RowMajor or \b #ColMajor, and of \b #SelfAdjoint
+ * The former controls \ref TopicStorageOrders "storage order", and defaults to
+ * column-major. The latter controls whether the matrix represents a selfadjoint
+ * matrix in which case either Supers of Subs have to be null.
*
* \sa class TridiagonalMatrix
*/
@@ -302,9 +302,9 @@ class BandMatrixWrapper : public BandMatrixBase(Index,Index) and
- * most of the time this is the only way it is used.
- *
- * However, if you want to directly maniputate block expressions,
- * for instance if you want to write a function returning such an expression, you
- * will need to use this class.
- *
- * Here is an example illustrating the dynamic case:
- * \include class_Block.cpp
- * Output: \verbinclude class_Block.out
- *
- * \note Even though this expression has dynamic size, in the case where \a XprType
- * has fixed size, this expression inherits a fixed maximal size which means that evaluating
- * it does not cause a dynamic memory allocation.
- *
- * Here is an example illustrating the fixed-size case:
- * \include class_FixedBlock.cpp
- * Output: \verbinclude class_FixedBlock.out
- *
- * \sa DenseBase::block(Index,Index,Index,Index), DenseBase::block(Index,Index), class VectorBlock
- */
-
namespace internal {
template
struct traits > : traits
@@ -101,6 +66,40 @@ template class BlockImpl;
+/** \class Block
+ * \ingroup Core_Module
+ *
+ * \brief Expression of a fixed-size or dynamic-size block
+ *
+ * \tparam XprType the type of the expression in which we are taking a block
+ * \tparam BlockRows the number of rows of the block we are taking at compile time (optional)
+ * \tparam BlockCols the number of columns of the block we are taking at compile time (optional)
+ * \tparam InnerPanel is true, if the block maps to a set of rows of a row major matrix or
+ * to set of columns of a column major matrix (optional). The parameter allows to determine
+ * at compile time whether aligned access is possible on the block expression.
+ *
+ * This class represents an expression of either a fixed-size or dynamic-size block. It is the return
+ * type of DenseBase::block(Index,Index,Index,Index) and DenseBase::block(Index,Index) and
+ * most of the time this is the only way it is used.
+ *
+ * However, if you want to directly maniputate block expressions,
+ * for instance if you want to write a function returning such an expression, you
+ * will need to use this class.
+ *
+ * Here is an example illustrating the dynamic case:
+ * \include class_Block.cpp
+ * Output: \verbinclude class_Block.out
+ *
+ * \note Even though this expression has dynamic size, in the case where \a XprType
+ * has fixed size, this expression inherits a fixed maximal size which means that evaluating
+ * it does not cause a dynamic memory allocation.
+ *
+ * Here is an example illustrating the fixed-size case:
+ * \include class_FixedBlock.cpp
+ * Output: \verbinclude class_FixedBlock.out
+ *
+ * \sa DenseBase::block(Index,Index,Index,Index), DenseBase::block(Index,Index), class VectorBlock
+ */
template class Block
: public BlockImpl::StorageKind>
{
@@ -130,8 +129,8 @@ template class
: Impl(xpr, startRow, startCol)
{
EIGEN_STATIC_ASSERT(RowsAtCompileTime!=Dynamic && ColsAtCompileTime!=Dynamic,THIS_METHOD_IS_ONLY_FOR_FIXED_SIZE)
- eigen_assert(startRow >= 0 && BlockRows >= 1 && startRow + BlockRows <= xpr.rows()
- && startCol >= 0 && BlockCols >= 1 && startCol + BlockCols <= xpr.cols());
+ eigen_assert(startRow >= 0 && BlockRows >= 0 && startRow + BlockRows <= xpr.rows()
+ && startCol >= 0 && BlockCols >= 0 && startCol + BlockCols <= xpr.cols());
}
/** Dynamic-size constructor
@@ -174,6 +173,7 @@ template >::type
{
typedef Block BlockType;
+ typedef typename internal::ref_selector::non_const_type XprTypeNested;
public:
typedef typename internal::dense_xpr_base::type Base;
@@ -222,15 +222,13 @@ template
inline PacketScalar packet(Index rowId, Index colId) const
{
- return m_xpr.template packet
- (rowId + m_startRow.value(), colId + m_startCol.value());
+ return m_xpr.template packet(rowId + m_startRow.value(), colId + m_startCol.value());
}
template
inline void writePacket(Index rowId, Index colId, const PacketScalar& val)
{
- m_xpr.const_cast_derived().template writePacket
- (rowId + m_startRow.value(), colId + m_startCol.value(), val);
+ m_xpr.template writePacket(rowId + m_startRow.value(), colId + m_startCol.value(), val);
}
template
@@ -289,7 +282,7 @@ template
inline void writePacket(Index index, const PacketScalar& val)
{
- m_xpr.const_cast_derived().template writePacket
+ m_xpr.template writePacket
(m_startRow.value() + (RowsAtCompileTime == 1 ? 0 : index),
m_startCol.value() + (RowsAtCompileTime == 1 ? index : 0), val);
}
@@ -302,10 +295,13 @@ template::type& nestedExpression() const
+ const typename internal::remove_all::type& nestedExpression() const
{
return m_xpr;
}
+
+ EIGEN_DEVICE_FUNC
+ XprType& nestedExpression() { return m_xpr; }
EIGEN_DEVICE_FUNC
StorageIndex startRow() const
@@ -321,9 +317,9 @@ template m_startRow;
- const internal::variable_if_dynamic m_startCol;
+ XprTypeNested m_xpr;
+ const internal::variable_if_dynamic m_startRow;
+ const internal::variable_if_dynamic m_startCol;
const internal::variable_if_dynamic m_blockRows;
const internal::variable_if_dynamic m_blockCols;
};
@@ -334,6 +330,7 @@ class BlockImpl_dense
: public MapBase >
{
typedef Block BlockType;
+ typedef typename internal::ref_selector::non_const_type XprTypeNested;
enum {
XprTypeIsRowMajor = (int(traits::Flags)&RowMajorBit) != 0
};
@@ -351,7 +348,9 @@ class BlockImpl_dense
|| ((BlockRows==XprType::RowsAtCompileTime) && (BlockCols==1) && ( XprTypeIsRowMajor)) ? xpr.innerStride() : xpr.outerStride()),
BlockRows==1 ? 1 : xpr.rows(),
BlockCols==1 ? 1 : xpr.cols()),
- m_xpr(xpr)
+ m_xpr(xpr),
+ m_startRow( (BlockRows==1) && (BlockCols==XprType::ColsAtCompileTime) ? i : 0),
+ m_startCol( (BlockRows==XprType::RowsAtCompileTime) && (BlockCols==1) ? i : 0)
{
init();
}
@@ -361,7 +360,7 @@ class BlockImpl_dense
EIGEN_DEVICE_FUNC
inline BlockImpl_dense(XprType& xpr, Index startRow, Index startCol)
: Base(xpr.data()+xpr.innerStride()*(XprTypeIsRowMajor?startCol:startRow) + xpr.outerStride()*(XprTypeIsRowMajor?startRow:startCol)),
- m_xpr(xpr)
+ m_xpr(xpr), m_startRow(startRow), m_startCol(startCol)
{
init();
}
@@ -373,16 +372,19 @@ class BlockImpl_dense