GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-11 19:29:02 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merged latest changes from the parent

Browse Source
This commit is contained in:
Benoit Steiner 2014-03-18 12:58:08 -07:00
commit 8a0845ebd7
36 changed files with 373 additions and 175 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
CTestConfig.cmake
View File

@ -4,10 +4,14 @@
## # The following are required to uses Dart and the Cdash dashboard
## ENABLE_TESTING()
## INCLUDE(CTest)
set(CTEST_PROJECT_NAME "Eigen3.2")
set(CTEST_PROJECT_NAME "Eigen")
set(CTEST_NIGHTLY_START_TIME "00:00:00 UTC")
set(CTEST_DROP_METHOD "http")
set(CTEST_DROP_SITE "manao.inria.fr")
set(CTEST_DROP_LOCATION "/CDash/submit.php?project=Eigen3.2")
set(CTEST_DROP_LOCATION "/CDash/submit.php?project=Eigen")
set(CTEST_DROP_SITE_CDASH TRUE)
set(CTEST_PROJECT_SUBPROJECTS
Official
Unsupported
)

8
Eigen/src/Cholesky/LDLT.h
View File

@ -291,13 +291,6 @@ template<> struct ldlt_inplace<Lower>
cutoff = abs(NumTraits<Scalar>::epsilon() * biggest_in_corner);
}
// Finish early if the matrix is not full rank.
if(biggest_in_corner < cutoff)
{
for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i;
break;
}
transpositions.coeffRef(k) = index_of_biggest_in_corner;
if(k != index_of_biggest_in_corner)
{
@ -333,6 +326,7 @@ template<> struct ldlt_inplace<Lower>
if(rs>0)
A21.noalias() -= A20 * temp.head(k);
}
if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff))
A21 /= mat.coeffRef(k,k);

2
Eigen/src/Core/ArrayWrapper.h
View File

@ -49,7 +49,7 @@ class ArrayWrapper : public ArrayBase<ArrayWrapper<ExpressionType> >
typedef typename internal::nested<ExpressionType>::type NestedExpressionType;
EIGEN_DEVICE_FUNC
inline ArrayWrapper(ExpressionType& matrix) : m_expression(matrix) {}
EIGEN_STRONG_INLINE ArrayWrapper(ExpressionType& matrix) : m_expression(matrix) {}
EIGEN_DEVICE_FUNC
inline Index rows() const { return m_expression.rows(); }

4
Eigen/src/Core/AssignEvaluator.h
View File

@ -574,13 +574,13 @@ public:
template<int StoreMode, int LoadMode>
void assignPacket(Index row, Index col)
{
m_functor.assignPacket<StoreMode>(&m_dst.coeffRef(row,col), m_src.template packet<LoadMode>(row,col));
m_functor.template assignPacket<StoreMode>(&m_dst.coeffRef(row,col), m_src.template packet<LoadMode>(row,col));
}
template<int StoreMode, int LoadMode>
void assignPacket(Index index)
{
m_functor.assignPacket<StoreMode>(&m_dst.coeffRef(index), m_src.template packet<LoadMode>(index));
m_functor.template assignPacket<StoreMode>(&m_dst.coeffRef(index), m_src.template packet<LoadMode>(index));
}
template<int StoreMode, int LoadMode>

12
Eigen/src/Core/CommaInitializer.h
View File

@ -45,6 +45,18 @@ struct CommaInitializer
m_xpr.block(0, 0, other.rows(), other.cols()) = other;
}
/* Copy/Move constructor which transfers ownership. This is crucial in
* absence of return value optimization to avoid assertions during destruction. */
// FIXME in C++11 mode this could be replaced by a proper RValue constructor
EIGEN_DEVICE_FUNC
inline CommaInitializer(const CommaInitializer& o)
: m_xpr(o.m_xpr), m_row(o.m_row), m_col(o.m_col), m_currentBlockRows(o.m_currentBlockRows) {
// Mark original object as finished. In absence of R-value references we need to const_cast:
const_cast<CommaInitializer&>(o).m_row = m_xpr.rows();
const_cast<CommaInitializer&>(o).m_col = m_xpr.cols();
const_cast<CommaInitializer&>(o).m_currentBlockRows = 0;
}
/* inserts a scalar value in the target matrix */
EIGEN_DEVICE_FUNC
CommaInitializer& operator,(const Scalar& s)

8
Eigen/src/Core/MatrixBase.h
View File

@ -356,8 +356,6 @@ template<typename Derived> class MatrixBase
Scalar trace() const;
/////////// Array module ///////////
template<int p> EIGEN_DEVICE_FUNC RealScalar lpNorm() const;
EIGEN_DEVICE_FUNC MatrixBase<Derived>& matrix() { return *this; }
@ -365,8 +363,10 @@ template<typename Derived> class MatrixBase
/** \returns an \link Eigen::ArrayBase Array \endlink expression of this matrix
* \sa ArrayBase::matrix() */
EIGEN_DEVICE_FUNC ArrayWrapper<Derived> array() { return derived(); }
EIGEN_DEVICE_FUNC const ArrayWrapper<const Derived> array() const { return derived(); }
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE ArrayWrapper<Derived> array() { return derived(); }
/** \returns a const \link Eigen::ArrayBase Array \endlink expression of this matrix
* \sa ArrayBase::matrix() */
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const ArrayWrapper<const Derived> array() const { return derived(); }
/////////// LU module ///////////

6
Eigen/src/Core/Ref.h
View File

@ -101,7 +101,7 @@ struct traits<Ref<_PlainObjectType, _Options, _StrideType> >
template<typename Derived> struct match {
enum {
HasDirectAccess = internal::has_direct_access<Derived>::ret,
StorageOrderMatch = PlainObjectType::IsVectorAtCompileTime || ((PlainObjectType::Flags&RowMajorBit)==(Derived::Flags&RowMajorBit)),
StorageOrderMatch = PlainObjectType::IsVectorAtCompileTime || Derived::IsVectorAtCompileTime || ((PlainObjectType::Flags&RowMajorBit)==(Derived::Flags&RowMajorBit)),
InnerStrideMatch = int(StrideType::InnerStrideAtCompileTime)==int(Dynamic)
|| int(StrideType::InnerStrideAtCompileTime)==int(Derived::InnerStrideAtCompileTime)
|| (int(StrideType::InnerStrideAtCompileTime)==0 && int(Derived::InnerStrideAtCompileTime)==1),
@ -172,6 +172,10 @@ protected:
}
else
::new (static_cast<Base*>(this)) Base(expr.data(), expr.rows(), expr.cols());
if(Expression::IsVectorAtCompileTime && (!PlainObjectType::IsVectorAtCompileTime) && ((Expression::Flags&RowMajorBit)!=(PlainObjectType::Flags&RowMajorBit)))
::new (&m_stride) StrideBase(expr.innerStride(), StrideType::InnerStrideAtCompileTime==0?0:1);
else
::new (&m_stride) StrideBase(StrideType::OuterStrideAtCompileTime==0?0:expr.outerStride(),
StrideType::InnerStrideAtCompileTime==0?0:expr.innerStride());
}

7
Eigen/src/Core/products/GeneralBlockPanelKernel.h
View File

@ -1277,6 +1277,7 @@ EIGEN_DONT_INLINE void gemm_pack_rhs<Scalar, Index, nr, ColMajor, Conjugate, Pan
template<typename Scalar, typename Index, int nr, bool Conjugate, bool PanelMode>
struct gemm_pack_rhs<Scalar, Index, nr, RowMajor, Conjugate, PanelMode>
{
typedef typename packet_traits<Scalar>::type Packet;
enum { PacketSize = packet_traits<Scalar>::size };
EIGEN_DONT_INLINE void operator()(Scalar* blockB, const Scalar* rhs, Index rhsStride, Index depth, Index cols, Index stride=0, Index offset=0);
};
@ -1298,6 +1299,11 @@ EIGEN_DONT_INLINE void gemm_pack_rhs<Scalar, Index, nr, RowMajor, Conjugate, Pan
if(PanelMode) count += nr * offset;
for(Index k=0; k<depth; k++)
{
if (nr == PacketSize) {
Packet A = ploadu<Packet>(&rhs[k*rhsStride + j2]);
pstoreu(blockB+count, cj.pconj(A));
count += PacketSize;
} else {
const Scalar* b0 = &rhs[k*rhsStride + j2];
blockB[count+0] = cj(b0[0]);
blockB[count+1] = cj(b0[1]);
@ -1305,6 +1311,7 @@ EIGEN_DONT_INLINE void gemm_pack_rhs<Scalar, Index, nr, RowMajor, Conjugate, Pan
if(nr==4) blockB[count+3] = cj(b0[3]);
count += nr;
}
}
// skip what we have after
if(PanelMode) count += nr * (stride-offset-depth);
}

6
Eigen/src/Core/products/SelfadjointMatrixVector.h
View File

@ -113,9 +113,9 @@ EIGEN_DONT_INLINE void selfadjoint_matrix_vector_product<Scalar,Index,StorageOrd
for (size_t i=starti; i<alignedStart; ++i)
{
res[i] += t0 * A0[i] + t1 * A1[i];
t2 += numext::conj(A0[i]) * rhs[i];
t3 += numext::conj(A1[i]) * rhs[i];
res[i] += cj0.pmul(A0[i], t0) + cj0.pmul(A1[i],t1);
t2 += cj1.pmul(A0[i], rhs[i]);
t3 += cj1.pmul(A1[i], rhs[i]);
}
// Yes this an optimization for gcc 4.3 and 4.4 (=> huge speed up)
// gcc 4.2 does this optimization automatically.

8
Eigen/src/Core/util/Memory.h
View File

@ -274,12 +274,12 @@ inline void* aligned_realloc(void *ptr, size_t new_size, size_t old_size)
// The defined(_mm_free) is just here to verify that this MSVC version
// implements _mm_malloc/_mm_free based on the corresponding _aligned_
// functions. This may not always be the case and we just try to be safe.
#if defined(_MSC_VER) && defined(_mm_free)
#if defined(_MSC_VER) && (!defined(_WIN32_WCE)) && defined(_mm_free)
result = _aligned_realloc(ptr,new_size,EIGEN_ALIGN_BYTES);
#else
result = generic_aligned_realloc(ptr,new_size,old_size);
#endif
#elif defined(_MSC_VER)
#elif defined(_MSC_VER) && (!defined(_WIN32_WCE))
result = _aligned_realloc(ptr,new_size,EIGEN_ALIGN_BYTES);
#else
result = handmade_aligned_realloc(ptr,new_size,old_size);
@ -464,7 +464,7 @@ template<typename T, bool Align> inline void conditional_aligned_delete_auto(T *
* There is also the variant first_aligned(const MatrixBase&) defined in DenseCoeffsBase.h.
*/
template<typename Scalar, typename Index>
static inline Index first_aligned(const Scalar* array, Index size)
inline Index first_aligned(const Scalar* array, Index size)
{
enum { PacketSize = packet_traits<Scalar>::size,
PacketAlignedMask = PacketSize-1
@ -492,7 +492,7 @@ static inline Index first_aligned(const Scalar* array, Index size)
/** \internal Returns the smallest integer multiple of \a base and greater or equal to \a size
*/
template<typename Index>
inline static Index first_multiple(Index size, Index base)
inline Index first_multiple(Index size, Index base)
{
return ((size+base-1)/base)*base;
}

21
Eigen/src/Geometry/Quaternion.h
View File

@ -34,8 +34,9 @@ struct quaternionbase_assign_impl;
template<class Derived>
class QuaternionBase : public RotationBase<Derived, 3>
{
public:
typedef RotationBase<Derived, 3> Base;
public:
using Base::operator*;
using Base::derived;
@ -203,6 +204,8 @@ public:
* \li \c Quaternionf for \c float
* \li \c Quaterniond for \c double
*
* \warning Operations interpreting the quaternion as rotation have undefined behavior if the quaternion is not normalized.
*
* \sa class AngleAxis, class Transform
*/
@ -223,10 +226,10 @@ struct traits<Quaternion<_Scalar,_Options> >
template<typename _Scalar, int _Options>
class Quaternion : public QuaternionBase<Quaternion<_Scalar,_Options> >
{
public:
typedef QuaternionBase<Quaternion<_Scalar,_Options> > Base;
enum { IsAligned = internal::traits<Quaternion>::IsAligned };
public:
typedef _Scalar Scalar;
EIGEN_INHERIT_ASSIGNMENT_EQUAL_OPERATOR(Quaternion)
@ -334,9 +337,9 @@ template<typename _Scalar, int _Options>
class Map<const Quaternion<_Scalar>, _Options >
: public QuaternionBase<Map<const Quaternion<_Scalar>, _Options> >
{
public:
typedef QuaternionBase<Map<const Quaternion<_Scalar>, _Options> > Base;
public:
typedef _Scalar Scalar;
typedef typename internal::traits<Map>::Coefficients Coefficients;
EIGEN_INHERIT_ASSIGNMENT_EQUAL_OPERATOR(Map)
@ -344,7 +347,7 @@ class Map<const Quaternion<_Scalar>, _Options >
/** Constructs a Mapped Quaternion object from the pointer \a coeffs
*
* The pointer \a coeffs must reference the four coeffecients of Quaternion in the following order:
* The pointer \a coeffs must reference the four coefficients of Quaternion in the following order:
* \code *coeffs == {x, y, z, w} \endcode
*
* If the template parameter _Options is set to #Aligned, then the pointer coeffs must be aligned. */
@ -371,9 +374,9 @@ template<typename _Scalar, int _Options>
class Map<Quaternion<_Scalar>, _Options >
: public QuaternionBase<Map<Quaternion<_Scalar>, _Options> >
{
public:
typedef QuaternionBase<Map<Quaternion<_Scalar>, _Options> > Base;
public:
typedef _Scalar Scalar;
typedef typename internal::traits<Map>::Coefficients Coefficients;
EIGEN_INHERIT_ASSIGNMENT_EQUAL_OPERATOR(Map)
@ -464,7 +467,7 @@ QuaternionBase<Derived>::_transformVector(Vector3 v) const
// Note that this algorithm comes from the optimization by hand
// of the conversion to a Matrix followed by a Matrix/Vector product.
// It appears to be much faster than the common algorithm found
// in the litterature (30 versus 39 flops). It also requires two
// in the literature (30 versus 39 flops). It also requires two
// Vector3 as temporaries.
Vector3 uv = this->vec().cross(v);
uv += uv;
@ -667,10 +670,10 @@ QuaternionBase<Derived>::angularDistance(const QuaternionBase<OtherDerived>& oth
{
using std::acos;
using std::abs;
double d = abs(this->dot(other));
if (d>=1.0)
Scalar d = abs(this->dot(other));
if (d>=Scalar(1))
return Scalar(0);
return static_cast<Scalar>(2 * acos(d));
return Scalar(2) * acos(d);
}

2
Eigen/src/Geometry/Scaling.h
View File

@ -65,7 +65,7 @@ public:
Transform<Scalar,Dim,(int(Mode)==int(Isometry)?Affine:Mode)> res = t;
res.prescale(factor());
return res;
}
}
/** Concatenates a uniform scaling and a linear transformation matrix */
// TODO returns an expression

4
Eigen/src/Geometry/Transform.h
View File

@ -530,9 +530,9 @@ public:
inline Transform& operator=(const UniformScaling<Scalar>& t);
inline Transform& operator*=(const UniformScaling<Scalar>& s) { return scale(s.factor()); }
inline Transform<Scalar,Dim,(int(Mode)==int(Isometry)?Affine:Mode)> operator*(const UniformScaling<Scalar>& s) const
inline TransformTimeDiagonalReturnType operator*(const UniformScaling<Scalar>& s) const
{
Transform<Scalar,Dim,(int(Mode)==int(Isometry)?Affine:Mode),Options> res = *this;
TransformTimeDiagonalReturnType res = *this;
res.scale(s.factor());
return res;
}

12
Eigen/src/PaStiXSupport/PaStiXSupport.h
View File

@ -12,6 +12,14 @@
namespace Eigen {
#if defined(DCOMPLEX)
#define PASTIX_COMPLEX COMPLEX
#define PASTIX_DCOMPLEX DCOMPLEX
#else
#define PASTIX_COMPLEX std::complex<float>
#define PASTIX_DCOMPLEX std::complex<double>
#endif
/** \ingroup PaStiXSupport_Module
* \brief Interface to the PaStix solver
*
@ -74,14 +82,14 @@ namespace internal
{
if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
if (nbrhs == 0) {x = NULL; nbrhs=1;}
c_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<COMPLEX*>(vals), perm, invp, reinterpret_cast<COMPLEX*>(x), nbrhs, iparm, dparm);
c_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<PASTIX_COMPLEX*>(vals), perm, invp, reinterpret_cast<PASTIX_COMPLEX*>(x), nbrhs, iparm, dparm);
}
void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, std::complex<double> *vals, int *perm, int * invp, std::complex<double> *x, int nbrhs, int *iparm, double *dparm)
{
if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
if (nbrhs == 0) {x = NULL; nbrhs=1;}
z_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<DCOMPLEX*>(vals), perm, invp, reinterpret_cast<DCOMPLEX*>(x), nbrhs, iparm, dparm);
z_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<PASTIX_DCOMPLEX*>(vals), perm, invp, reinterpret_cast<PASTIX_DCOMPLEX*>(x), nbrhs, iparm, dparm);
}
// Convert the matrix to Fortran-style Numbering

2
Eigen/src/QR/ColPivHouseholderQR.h
View File

@ -350,7 +350,7 @@ template<typename _MatrixType> class ColPivHouseholderQR
return m_usePrescribedThreshold ? m_prescribedThreshold
// this formula comes from experimenting (see "LU precision tuning" thread on the list)
// and turns out to be identical to Higham's formula used already in LDLt.
: NumTraits<Scalar>::epsilon() * m_qr.diagonalSize();
: NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize());
}
/** \returns the number of nonzero pivots in the QR decomposition.

4
Eigen/src/QR/FullPivHouseholderQR.h
View File

@ -346,7 +346,7 @@ template<typename _MatrixType> class FullPivHouseholderQR
return m_usePrescribedThreshold ? m_prescribedThreshold
// this formula comes from experimenting (see "LU precision tuning" thread on the list)
// and turns out to be identical to Higham's formula used already in LDLt.
: NumTraits<Scalar>::epsilon() * m_qr.diagonalSize();
: NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize());
}
/** \returns the number of nonzero pivots in the QR decomposition.
@ -417,7 +417,7 @@ FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons
m_temp.resize(cols);
m_precision = NumTraits<Scalar>::epsilon() * size;
m_precision = NumTraits<Scalar>::epsilon() * RealScalar(size);
m_rows_transpositions.resize(size);
m_cols_transpositions.resize(size);

5
Eigen/src/SparseCore/SparseBlock.h
View File

@ -338,7 +338,10 @@ const Block<const Derived,Dynamic,Dynamic,true> SparseMatrixBase<Derived>::inner
namespace internal {
template< typename XprType, int BlockRows, int BlockCols, bool InnerPanel,
bool OuterVector = (BlockCols==1 && XprType::IsRowMajor) || (BlockRows==1 && !XprType::IsRowMajor)>
bool OuterVector = (BlockCols==1 && XprType::IsRowMajor)
| // FIXME | instead of || to please GCC 4.4.0 stupid warning "suggest parentheses around &&".
// revert to || as soon as not needed anymore.
(BlockRows==1 && !XprType::IsRowMajor)>
class GenericSparseBlockInnerIteratorImpl;
}

8
bench/btl/CMakeLists.txt
View File

@ -48,6 +48,12 @@ include_directories(
# set(DEFAULT_LIBRARIES ${MKL_LIBRARIES})
# endif (MKL_FOUND)
find_library(EIGEN_BTL_RT_LIBRARY rt)
# if we cannot find it easily, then we don't need it!
if(NOT EIGEN_BTL_RT_LIBRARY)
set(EIGEN_BTL_RT_LIBRARY "")
endif()
MACRO(BTL_ADD_BENCH targetname)
foreach(_current_var ${ARGN})
@ -70,7 +76,7 @@ MACRO(BTL_ADD_BENCH targetname)
IF(BUILD_${targetname})
ADD_EXECUTABLE(${targetname} ${_sources})
ADD_TEST(${targetname} "${targetname}")
target_link_libraries(${targetname} ${DEFAULT_LIBRARIES} rt)
target_link_libraries(${targetname} ${DEFAULT_LIBRARIES} ${EIGEN_BTL_RT_LIBRARY})
ENDIF(BUILD_${targetname})
ENDMACRO(BTL_ADD_BENCH)

4
bench/btl/generic_bench/bench.hh
View File

@ -102,8 +102,8 @@ BTL_DONT_INLINE void bench( int size_min, int size_max, int nb_point )
// merge the two data
std::vector<int> newSizes;
std::vector<double> newFlops;
int i=0;
int j=0;
unsigned int i=0;
unsigned int j=0;
while (i<tab_sizes.size() && j<oldSizes.size())
{
if (tab_sizes[i] == oldSizes[j])

2
bench/btl/generic_bench/btl.hh
View File

@ -46,7 +46,7 @@
#if (defined __GNUC__) && (!defined __INTEL_COMPILER) && !defined(__arm__) && !defined(__powerpc__)
#define BTL_DISABLE_SSE_EXCEPTIONS() { \
int aux; \
int aux = 0; \
asm( \
"stmxcsr %[aux] \n\t" \
"orl $32832, %[aux] \n\t" \

10
bench/btl/generic_bench/init/init_matrix.hh
View File

@ -29,7 +29,7 @@ BTL_DONT_INLINE void init_row(Vector & X, int size, int row){
X.resize(size);
for (int j=0;j<X.size();j++){
for (unsigned int j=0;j<X.size();j++){
X[j]=typename Vector::value_type(init_function(row,j));
}
}
@ -42,7 +42,7 @@ BTL_DONT_INLINE void init_row(Vector & X, int size, int row){
template<double init_function(int,int),class Vector>
BTL_DONT_INLINE void init_matrix(Vector & A, int size){
A.resize(size);
for (int row=0; row<A.size() ; row++){
for (unsigned int row=0; row<A.size() ; row++){
init_row<init_function>(A[row],size,row);
}
}
@ -50,11 +50,11 @@ BTL_DONT_INLINE void init_matrix(Vector & A, int size){
template<double init_function(int,int),class Matrix>
BTL_DONT_INLINE void init_matrix_symm(Matrix& A, int size){
A.resize(size);
for (int row=0; row<A.size() ; row++)
for (unsigned int row=0; row<A.size() ; row++)
A[row].resize(size);
for (int row=0; row<A.size() ; row++){
for (unsigned int row=0; row<A.size() ; row++){
A[row][row] = init_function(row,row);
for (int col=0; col<row ; col++){
for (unsigned int col=0; col<row ; col++){
double x = init_function(row,col);
A[row][col] = A[col][row] = x;
}

2
bench/btl/generic_bench/init/init_vector.hh
View File

@ -29,7 +29,7 @@ void init_vector(Vector & X, int size){
X.resize(size);
for (int i=0;i<X.size();i++){
for (unsigned int i=0;i<X.size();i++){
X[i]=typename Vector::value_type(init_function(i));
}
}

2
bench/btl/generic_bench/timers/portable_perf_analyzer.hh
View File

@ -78,7 +78,7 @@ public:
// time measurement
action.calculate();
_chronos.start();
for (int ii=0;ii<_nb_calc;ii++)
for (unsigned int ii=0;ii<_nb_calc;ii++)
{
action.calculate();
}

46
bench/btl/generic_bench/timers/portable_timer.hh
View File

@ -34,7 +34,7 @@
// timer -------------------------------------------------------------------//
// A timer object measures CPU time.
#ifdef _MSC_VER
#if defined(_MSC_VER)
#define NOMINMAX
#include <windows.h>
@ -87,6 +87,48 @@
}; // Portable_Timer
#elif defined(__APPLE__)
#include <CoreServices/CoreServices.h>
#include <mach/mach_time.h>
class Portable_Timer
{
public:
Portable_Timer()
{
}
void start()
{
m_start_time = double(mach_absolute_time())*1e-9;;
}
void stop()
{
m_stop_time = double(mach_absolute_time())*1e-9;;
}
double elapsed()
{
return user_time();
}
double user_time()
{
return m_stop_time - m_start_time;
}
private:
double m_stop_time, m_start_time;
}; // Portable_Timer (Apple)
#else
#include <sys/time.h>
@ -138,7 +180,7 @@ private:
int m_clkid;
double m_stop_time, m_start_time;
}; // Portable_Timer
}; // Portable_Timer (Linux)
#endif

8
bench/btl/libs/eigen3/eigen3_interface.hh
View File

@ -52,8 +52,8 @@ public :
static BTL_DONT_INLINE void matrix_from_stl(gene_matrix & A, stl_matrix & A_stl){
A.resize(A_stl[0].size(), A_stl.size());
for (int j=0; j<A_stl.size() ; j++){
for (int i=0; i<A_stl[j].size() ; i++){
for (unsigned int j=0; j<A_stl.size() ; j++){
for (unsigned int i=0; i<A_stl[j].size() ; i++){
A.coeffRef(i,j) = A_stl[j][i];
}
}
@ -62,13 +62,13 @@ public :
static BTL_DONT_INLINE void vector_from_stl(gene_vector & B, stl_vector & B_stl){
B.resize(B_stl.size(),1);
for (int i=0; i<B_stl.size() ; i++){
for (unsigned int i=0; i<B_stl.size() ; i++){
B.coeffRef(i) = B_stl[i];
}
}
static BTL_DONT_INLINE void vector_to_stl(gene_vector & B, stl_vector & B_stl){
for (int i=0; i<B_stl.size() ; i++){
for (unsigned int i=0; i<B_stl.size() ; i++){
B_stl[i] = B.coeff(i);
}
}

2
doc/TopicLinearAlgebraDecompositions.dox
View File

@ -249,7 +249,7 @@ For an introduction on linear solvers and decompositions, check this \link Tutor
<dt><b>Implicit Multi Threading (MT)</b></dt>
<dd>Means the algorithm can take advantage of multicore processors via OpenMP. "Implicit" means the algortihm itself is not parallelized, but that it relies on parallelized matrix-matrix product rountines.</dd>
<dt><b>Explicit Multi Threading (MT)</b></dt>
<dd>Means the algorithm is explicitely parallelized to take advantage of multicore processors via OpenMP.</dd>
<dd>Means the algorithm is explicitly parallelized to take advantage of multicore processors via OpenMP.</dd>
<dt><b>Meta-unroller</b></dt>
<dd>Means the algorithm is automatically and explicitly unrolled for very small fixed size matrices.</dd>
<dt><b></b></dt>

2
doc/TopicMultithreading.dox
View File

@ -39,7 +39,7 @@ int main(int argc, char** argv)
}
\endcode
\warning note that all functions generating random matrices are \b not re-entrant nor thread-safe. Those include DenseBase::Random(), and DenseBase::setRandom() despite a call to Eigen::initParallel(). This is because these functions are based on std::rand which is not re-entrant. For thread-safe random generator, we recommend the use of boost::random of c++11 random feature.
\warning note that all functions generating random matrices are \b not re-entrant nor thread-safe. Those include DenseBase::Random(), and DenseBase::setRandom() despite a call to Eigen::initParallel(). This is because these functions are based on std::rand which is not re-entrant. For thread-safe random generator, we recommend the use of boost::random or c++11 random feature.
In the case your application is parallelized with OpenMP, you might want to disable Eigen's own parallization as detailed in the previous section.

2
doc/snippets/EigenSolver_eigenvectors.cpp
View File

@ -1,4 +1,4 @@
MatrixXd ones = MatrixXd::Ones(3,3);
EigenSolver<MatrixXd> es(ones);
cout << "The first eigenvector of the 3x3 matrix of ones is:"
<< endl << es.eigenvectors().col(1) << endl;
<< endl << es.eigenvectors().col(0) << endl;

55
test/CMakeLists.txt
View File

@ -13,11 +13,26 @@ if(NOT EXISTS ${CMAKE_CURRENT_BINARY_DIR}/split_test_helper.h)
endforeach()
endif()
# check if we have a Fortran compiler
include("../cmake/language_support.cmake")
workaround_9220(Fortran EIGEN_Fortran_COMPILER_WORKS)
if(EIGEN_Fortran_COMPILER_WORKS)
enable_language(Fortran OPTIONAL)
if(NOT CMAKE_Fortran_COMPILER)
set(EIGEN_Fortran_COMPILER_WORKS OFF)
endif()
endif()
if(NOT EIGEN_Fortran_COMPILER_WORKS)
# search for a default Lapack library to complete Eigen's one
find_package(LAPACK)
endif()
# configure blas/lapack (use Eigen's ones)
set(BLAS_FOUND TRUE)
set(LAPACK_FOUND TRUE)
set(BLAS_LIBRARIES eigen_blas)
set(LAPACK_LIBRARIES eigen_lapack)
set(EIGEN_BLAS_LIBRARIES eigen_blas)
set(EIGEN_LAPACK_LIBRARIES eigen_lapack)
set(EIGEN_TEST_MATRIX_DIR "" CACHE STRING "Enable testing of realword sparse matrices contained in the specified path")
if(EIGEN_TEST_MATRIX_DIR)
@ -32,33 +47,33 @@ endif(EIGEN_TEST_MATRIX_DIR)
set(SPARSE_LIBS " ")
find_package(Cholmod)
if(CHOLMOD_FOUND AND BLAS_FOUND AND LAPACK_FOUND)
if(CHOLMOD_FOUND)
add_definitions("-DEIGEN_CHOLMOD_SUPPORT")
include_directories(${CHOLMOD_INCLUDES})
set(SPARSE_LIBS ${SPARSE_LIBS} ${CHOLMOD_LIBRARIES} ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES})
set(CHOLMOD_ALL_LIBS ${CHOLMOD_LIBRARIES} ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES})
set(SPARSE_LIBS ${SPARSE_LIBS} ${CHOLMOD_LIBRARIES} ${EIGEN_BLAS_LIBRARIES} ${EIGEN_LAPACK_LIBRARIES})
set(CHOLMOD_ALL_LIBS ${CHOLMOD_LIBRARIES} ${EIGEN_BLAS_LIBRARIES} ${EIGEN_LAPACK_LIBRARIES})
ei_add_property(EIGEN_TESTED_BACKENDS "Cholmod, ")
else()
ei_add_property(EIGEN_MISSING_BACKENDS "Cholmod, ")
endif()
find_package(Umfpack)
if(UMFPACK_FOUND AND BLAS_FOUND)
if(UMFPACK_FOUND)
add_definitions("-DEIGEN_UMFPACK_SUPPORT")
include_directories(${UMFPACK_INCLUDES})
set(SPARSE_LIBS ${SPARSE_LIBS} ${UMFPACK_LIBRARIES} ${BLAS_LIBRARIES})
set(UMFPACK_ALL_LIBS ${UMFPACK_LIBRARIES} ${BLAS_LIBRARIES})
set(SPARSE_LIBS ${SPARSE_LIBS} ${UMFPACK_LIBRARIES} ${EIGEN_BLAS_LIBRARIES})
set(UMFPACK_ALL_LIBS ${UMFPACK_LIBRARIES} ${EIGEN_BLAS_LIBRARIES})
ei_add_property(EIGEN_TESTED_BACKENDS "UmfPack, ")
else()
ei_add_property(EIGEN_MISSING_BACKENDS "UmfPack, ")
endif()
find_package(SuperLU)
if(SUPERLU_FOUND AND BLAS_FOUND)
if(SUPERLU_FOUND)
add_definitions("-DEIGEN_SUPERLU_SUPPORT")
include_directories(${SUPERLU_INCLUDES})
set(SPARSE_LIBS ${SPARSE_LIBS} ${SUPERLU_LIBRARIES} ${BLAS_LIBRARIES})
set(SUPERLU_ALL_LIBS ${SUPERLU_LIBRARIES} ${BLAS_LIBRARIES})
set(SPARSE_LIBS ${SPARSE_LIBS} ${SUPERLU_LIBRARIES} ${EIGEN_BLAS_LIBRARIES})
set(SUPERLU_ALL_LIBS ${SUPERLU_LIBRARIES} ${EIGEN_BLAS_LIBRARIES})
ei_add_property(EIGEN_TESTED_BACKENDS "SuperLU, ")
else()
ei_add_property(EIGEN_MISSING_BACKENDS "SuperLU, ")
@ -68,7 +83,7 @@ endif()
find_package(Pastix)
find_package(Scotch)
find_package(Metis)
if(PASTIX_FOUND AND BLAS_FOUND)
if(PASTIX_FOUND)
add_definitions("-DEIGEN_PASTIX_SUPPORT")
include_directories(${PASTIX_INCLUDES})
if(SCOTCH_FOUND)
@ -80,8 +95,8 @@ if(PASTIX_FOUND AND BLAS_FOUND)
else(SCOTCH_FOUND)
ei_add_property(EIGEN_MISSING_BACKENDS "PaStiX, ")
endif(SCOTCH_FOUND)
set(SPARSE_LIBS ${SPARSE_LIBS} ${PASTIX_LIBRARIES} ${ORDERING_LIBRARIES} ${BLAS_LIBRARIES})
set(PASTIX_ALL_LIBS ${PASTIX_LIBRARIES} ${BLAS_LIBRARIES})
set(SPARSE_LIBS ${SPARSE_LIBS} ${PASTIX_LIBRARIES} ${ORDERING_LIBRARIES} ${EIGEN_BLAS_LIBRARIES})
set(PASTIX_ALL_LIBS ${PASTIX_LIBRARIES} ${EIGEN_BLAS_LIBRARIES})
ei_add_property(EIGEN_TESTED_BACKENDS "PaStiX, ")
else()
ei_add_property(EIGEN_MISSING_BACKENDS "PaStiX, ")
@ -96,16 +111,14 @@ else()
endif()
find_package(SPQR)
if(SPQR_FOUND AND BLAS_FOUND AND LAPACK_FOUND)
if(CHOLMOD_FOUND)
if(SPQR_FOUND AND CHOLMOD_FOUND AND (EIGEN_Fortran_COMPILER_WORKS OR LAPACK_FOUND) )
add_definitions("-DEIGEN_SPQR_SUPPORT")
include_directories(${SPQR_INCLUDES})
set(SPQR_ALL_LIBS ${SPQR_LIBRARIES} ${CHOLMOD_LIBRARIES} ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES})
set(SPQR_ALL_LIBS ${SPQR_LIBRARIES} ${CHOLMOD_LIBRARIES} ${EIGEN_LAPACK_LIBRARIES} ${EIGEN_BLAS_LIBRARIES} ${LAPACK_LIBRARIES})
set(SPARSE_LIBS ${SPARSE_LIBS} ${SPQR_ALL_LIBS})
ei_add_property(EIGEN_TESTED_BACKENDS "SPQR, ")
else(CHOLMOD_FOUND)
else()
ei_add_property(EIGEN_MISSING_BACKENDS "SPQR, ")
endif(CHOLMOD_FOUND)
endif()
option(EIGEN_TEST_NOQT "Disable Qt support in unit tests" OFF)

23
test/block.cpp
View File

@ -10,6 +10,26 @@
#define EIGEN_NO_STATIC_ASSERT // otherwise we fail at compile time on unused paths
#include "main.h"
template<typename MatrixType, typename Index, typename Scalar>
typename Eigen::internal::enable_if<!NumTraits<typename MatrixType::Scalar>::IsComplex,typename MatrixType::Scalar>::type
block_real_only(const MatrixType &m1, Index r1, Index r2, Index c1, Index c2, const Scalar& s1) {
// check cwise-Functions:
VERIFY_IS_APPROX(m1.row(r1).cwiseMax(s1), m1.cwiseMax(s1).row(r1));
VERIFY_IS_APPROX(m1.col(c1).cwiseMin(s1), m1.cwiseMin(s1).col(c1));
VERIFY_IS_APPROX(m1.block(r1,c1,r2-r1+1,c2-c1+1).cwiseMin(s1), m1.cwiseMin(s1).block(r1,c1,r2-r1+1,c2-c1+1));
VERIFY_IS_APPROX(m1.block(r1,c1,r2-r1+1,c2-c1+1).cwiseMax(s1), m1.cwiseMax(s1).block(r1,c1,r2-r1+1,c2-c1+1));
return Scalar(0);
}
template<typename MatrixType, typename Index, typename Scalar>
typename Eigen::internal::enable_if<NumTraits<typename MatrixType::Scalar>::IsComplex,typename MatrixType::Scalar>::type
block_real_only(const MatrixType &, Index, Index, Index, Index, const Scalar&) {
return Scalar(0);
}
template<typename MatrixType> void block(const MatrixType& m)
{
typedef typename MatrixType::Index Index;
@ -37,6 +57,8 @@ template<typename MatrixType> void block(const MatrixType& m)
Index c1 = internal::random<Index>(0,cols-1);
Index c2 = internal::random<Index>(c1,cols-1);
block_real_only(m1, r1, r2, c1, c1, s1);
//check row() and col()
VERIFY_IS_EQUAL(m1.col(c1).transpose(), m1.transpose().row(c1));
//check operator(), both constant and non-constant, on row() and col()
@ -52,6 +74,7 @@ template<typename MatrixType> void block(const MatrixType& m)
m1.col(c1).col(0) += s1 * m1_copy.col(c2);
VERIFY_IS_APPROX(m1.col(c1), m1_copy.col(c1) + Scalar(2) * s1 * m1_copy.col(c2));
//check block()
Matrix<Scalar,Dynamic,Dynamic> b1(1,1); b1(0,0) = m1(r1,c1);

32
test/cholesky.cpp
View File

@ -179,6 +179,38 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
// restore
if(sign == -1)
symm = -symm;
// check matrices coming from linear constraints with Lagrange multipliers
if(rows>=3)
{
SquareMatrixType A = symm;
int c = internal::random<int>(0,rows-2);
A.bottomRightCorner(c,c).setZero();
// Make sure a solution exists:
vecX.setRandom();
vecB = A * vecX;
vecX.setZero();
ldltlo.compute(A);
VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix());
vecX = ldltlo.solve(vecB);
VERIFY_IS_APPROX(A * vecX, vecB);
}
// check non-full rank matrices
if(rows>=3)
{
int r = internal::random<int>(1,rows-1);
Matrix<Scalar,Dynamic,Dynamic> a = Matrix<Scalar,Dynamic,Dynamic>::Random(rows,r);
SquareMatrixType A = a * a.adjoint();
// Make sure a solution exists:
vecX.setRandom();
vecB = A * vecX;
vecX.setZero();
ldltlo.compute(A);
VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix());
vecX = ldltlo.solve(vecB);
VERIFY_IS_APPROX(A * vecX, vecB);
}
}
// update/downdate

104
test/ref.cpp
View File

@ -154,59 +154,79 @@ template<typename PlainObjectType> void check_const_correctness(const PlainObjec
VERIFY( !(Ref<ConstPlainObjectType, Aligned>::Flags & LvalueBit) );
}
EIGEN_DONT_INLINE void call_ref_1(Ref<VectorXf> ) { }
EIGEN_DONT_INLINE void call_ref_2(const Ref<const VectorXf>& ) { }
EIGEN_DONT_INLINE void call_ref_3(Ref<VectorXf,0,InnerStride<> > ) { }
EIGEN_DONT_INLINE void call_ref_4(const Ref<const VectorXf,0,InnerStride<> >& ) { }
EIGEN_DONT_INLINE void call_ref_5(Ref<MatrixXf,0,OuterStride<> > ) { }
EIGEN_DONT_INLINE void call_ref_6(const Ref<const MatrixXf,0,OuterStride<> >& ) { }
template<typename B>
EIGEN_DONT_INLINE void call_ref_1(Ref<VectorXf> a, const B &b) { VERIFY_IS_EQUAL(a,b); }
template<typename B>
EIGEN_DONT_INLINE void call_ref_2(const Ref<const VectorXf>& a, const B &b) { VERIFY_IS_EQUAL(a,b); }
template<typename B>
EIGEN_DONT_INLINE void call_ref_3(Ref<VectorXf,0,InnerStride<> > a, const B &b) { VERIFY_IS_EQUAL(a,b); }
template<typename B>
EIGEN_DONT_INLINE void call_ref_4(const Ref<const VectorXf,0,InnerStride<> >& a, const B &b) { VERIFY_IS_EQUAL(a,b); }
template<typename B>
EIGEN_DONT_INLINE void call_ref_5(Ref<MatrixXf,0,OuterStride<> > a, const B &b) { VERIFY_IS_EQUAL(a,b); }
template<typename B>
EIGEN_DONT_INLINE void call_ref_6(const Ref<const MatrixXf,0,OuterStride<> >& a, const B &b) { VERIFY_IS_EQUAL(a,b); }
template<typename B>
EIGEN_DONT_INLINE void call_ref_7(Ref<Matrix<float,Dynamic,3> > a, const B &b) { VERIFY_IS_EQUAL(a,b); }
void call_ref()
{
VectorXcf ca(10);
VectorXf a(10);
VectorXcf ca = VectorXcf::Random(10);
VectorXf a = VectorXf::Random(10);
RowVectorXf b = RowVectorXf::Random(10);
MatrixXf A = MatrixXf::Random(10,10);
RowVector3f c = RowVector3f::Random();
const VectorXf& ac(a);
VectorBlock<VectorXf> ab(a,0,3);
MatrixXf A(10,10);
const VectorBlock<VectorXf> abc(a,0,3);
VERIFY_EVALUATION_COUNT( call_ref_1(a), 0);
//call_ref_1(ac); // does not compile because ac is const
VERIFY_EVALUATION_COUNT( call_ref_1(ab), 0);
VERIFY_EVALUATION_COUNT( call_ref_1(a.head(4)), 0);
VERIFY_EVALUATION_COUNT( call_ref_1(abc), 0);
VERIFY_EVALUATION_COUNT( call_ref_1(A.col(3)), 0);
// call_ref_1(A.row(3)); // does not compile because innerstride!=1
VERIFY_EVALUATION_COUNT( call_ref_3(A.row(3)), 0);
VERIFY_EVALUATION_COUNT( call_ref_4(A.row(3)), 0);
//call_ref_1(a+a); // does not compile for obvious reason
VERIFY_EVALUATION_COUNT( call_ref_2(A*A.col(1)), 1); // evaluated into a temp
VERIFY_EVALUATION_COUNT( call_ref_2(ac.head(5)), 0);
VERIFY_EVALUATION_COUNT( call_ref_2(ac), 0);
VERIFY_EVALUATION_COUNT( call_ref_2(a), 0);
VERIFY_EVALUATION_COUNT( call_ref_2(ab), 0);
VERIFY_EVALUATION_COUNT( call_ref_2(a.head(4)), 0);
VERIFY_EVALUATION_COUNT( call_ref_2(a+a), 1); // evaluated into a temp
VERIFY_EVALUATION_COUNT( call_ref_2(ca.imag()), 1); // evaluated into a temp
VERIFY_EVALUATION_COUNT( call_ref_1(a,a), 0);
VERIFY_EVALUATION_COUNT( call_ref_1(b,b.transpose()), 0);
// call_ref_1(ac); // does not compile because ac is const
VERIFY_EVALUATION_COUNT( call_ref_1(ab,ab), 0);
VERIFY_EVALUATION_COUNT( call_ref_1(a.head(4),a.head(4)), 0);
VERIFY_EVALUATION_COUNT( call_ref_1(abc,abc), 0);
VERIFY_EVALUATION_COUNT( call_ref_1(A.col(3),A.col(3)), 0);
// call_ref_1(A.row(3)); // does not compile because innerstride!=1
VERIFY_EVALUATION_COUNT( call_ref_3(A.row(3),A.row(3).transpose()), 0);
VERIFY_EVALUATION_COUNT( call_ref_4(A.row(3),A.row(3).transpose()), 0);
// call_ref_1(a+a); // does not compile for obvious reason
VERIFY_EVALUATION_COUNT( call_ref_4(ac.head(5)), 0);
VERIFY_EVALUATION_COUNT( call_ref_4(a+a), 1); // evaluated into a temp
VERIFY_EVALUATION_COUNT( call_ref_4(ca.imag()), 0);
MatrixXf tmp = A*A.col(1);
VERIFY_EVALUATION_COUNT( call_ref_2(A*A.col(1), tmp), 1); // evaluated into a temp
VERIFY_EVALUATION_COUNT( call_ref_2(ac.head(5),ac.head(5)), 0);
VERIFY_EVALUATION_COUNT( call_ref_2(ac,ac), 0);
VERIFY_EVALUATION_COUNT( call_ref_2(a,a), 0);
VERIFY_EVALUATION_COUNT( call_ref_2(ab,ab), 0);
VERIFY_EVALUATION_COUNT( call_ref_2(a.head(4),a.head(4)), 0);
tmp = a+a;
VERIFY_EVALUATION_COUNT( call_ref_2(a+a,tmp), 1); // evaluated into a temp
VERIFY_EVALUATION_COUNT( call_ref_2(ca.imag(),ca.imag()), 1); // evaluated into a temp
VERIFY_EVALUATION_COUNT( call_ref_5(a), 0);
VERIFY_EVALUATION_COUNT( call_ref_5(a.head(3)), 0);
VERIFY_EVALUATION_COUNT( call_ref_5(A), 0);
// call_ref_5(A.transpose()); // does not compile
VERIFY_EVALUATION_COUNT( call_ref_5(A.block(1,1,2,2)), 0);
VERIFY_EVALUATION_COUNT( call_ref_4(ac.head(5),ac.head(5)), 0);
tmp = a+a;
VERIFY_EVALUATION_COUNT( call_ref_4(a+a,tmp), 1); // evaluated into a temp
VERIFY_EVALUATION_COUNT( call_ref_4(ca.imag(),ca.imag()), 0);
VERIFY_EVALUATION_COUNT( call_ref_6(a), 0);
VERIFY_EVALUATION_COUNT( call_ref_6(a.head(3)), 0);
VERIFY_EVALUATION_COUNT( call_ref_6(A.row(3)), 1); // evaluated into a temp thouth it could be avoided by viewing it as a 1xn matrix
VERIFY_EVALUATION_COUNT( call_ref_6(A+A), 1); // evaluated into a temp
VERIFY_EVALUATION_COUNT( call_ref_6(A), 0);
VERIFY_EVALUATION_COUNT( call_ref_6(A.transpose()), 1); // evaluated into a temp because the storage orders do not match
VERIFY_EVALUATION_COUNT( call_ref_6(A.block(1,1,2,2)), 0);
VERIFY_EVALUATION_COUNT( call_ref_5(a,a), 0);
VERIFY_EVALUATION_COUNT( call_ref_5(a.head(3),a.head(3)), 0);
VERIFY_EVALUATION_COUNT( call_ref_5(A,A), 0);
// call_ref_5(A.transpose()); // does not compile
VERIFY_EVALUATION_COUNT( call_ref_5(A.block(1,1,2,2),A.block(1,1,2,2)), 0);
VERIFY_EVALUATION_COUNT( call_ref_5(b,b), 0); // storage order do not match, but this is a degenerate case that should work
VERIFY_EVALUATION_COUNT( call_ref_5(a.row(3),a.row(3)), 0);
VERIFY_EVALUATION_COUNT( call_ref_6(a,a), 0);
VERIFY_EVALUATION_COUNT( call_ref_6(a.head(3),a.head(3)), 0);
VERIFY_EVALUATION_COUNT( call_ref_6(A.row(3),A.row(3)), 1); // evaluated into a temp thouth it could be avoided by viewing it as a 1xn matrix
tmp = A+A;
VERIFY_EVALUATION_COUNT( call_ref_6(A+A,tmp), 1); // evaluated into a temp
VERIFY_EVALUATION_COUNT( call_ref_6(A,A), 0);
VERIFY_EVALUATION_COUNT( call_ref_6(A.transpose(),A.transpose()), 1); // evaluated into a temp because the storage orders do not match
VERIFY_EVALUATION_COUNT( call_ref_6(A.block(1,1,2,2),A.block(1,1,2,2)), 0);
VERIFY_EVALUATION_COUNT( call_ref_7(c,c), 0);
}
void test_ref()

2
test/sizeof.cpp
View File

@ -21,6 +21,8 @@ template<typename MatrixType> void verifySizeOf(const MatrixType&)
void test_sizeof()
{
CALL_SUBTEST(verifySizeOf(Matrix<float, 1, 1>()) );
CALL_SUBTEST(verifySizeOf(Vector2d()) );
CALL_SUBTEST(verifySizeOf(Vector4f()) );
CALL_SUBTEST(verifySizeOf(Matrix4d()) );
CALL_SUBTEST(verifySizeOf(Matrix<double, 4, 2>()) );
CALL_SUBTEST(verifySizeOf(Matrix<bool, 7, 5>()) );

49
unsupported/Eigen/src/IterativeSolvers/MINRES.h
View File

@ -37,22 +37,31 @@ namespace Eigen {
typedef typename Dest::Scalar Scalar;
typedef Matrix<Scalar,Dynamic,1> VectorType;
// Check for zero rhs
const RealScalar rhsNorm2(rhs.squaredNorm());
if(rhsNorm2 == 0)
{
x.setZero();
iters = 0;
tol_error = 0;
return;
}
// initialize
const int maxIters(iters); // initialize maxIters to iters
const int N(mat.cols()); // the size of the matrix
const RealScalar rhsNorm2(rhs.squaredNorm());
const RealScalar threshold2(tol_error*tol_error*rhsNorm2); // convergence threshold (compared to residualNorm2)
// Initialize preconditioned Lanczos
// VectorType v_old(N); // will be initialized inside loop
VectorType v_old(N); // will be initialized inside loop
VectorType v( VectorType::Zero(N) ); //initialize v
VectorType v_new(rhs-mat*x); //initialize v_new
RealScalar residualNorm2(v_new.squaredNorm());
// VectorType w(N); // will be initialized inside loop
VectorType w(N); // will be initialized inside loop
VectorType w_new(precond.solve(v_new)); // initialize w_new
// RealScalar beta; // will be initialized inside loop
RealScalar beta_new2(v_new.dot(w_new));
eigen_assert(beta_new2 >= 0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE");
eigen_assert(beta_new2 >= 0.0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE");
RealScalar beta_new(sqrt(beta_new2));
const RealScalar beta_one(beta_new);
v_new /= beta_new;
@ -62,14 +71,14 @@ namespace Eigen {
RealScalar c_old(1.0);
RealScalar s(0.0); // the sine of the Givens rotation
RealScalar s_old(0.0); // the sine of the Givens rotation
// VectorType p_oold(N); // will be initialized in loop
VectorType p_oold(N); // will be initialized in loop
VectorType p_old(VectorType::Zero(N)); // initialize p_old=0
VectorType p(p_old); // initialize p=0
RealScalar eta(1.0);
iters = 0; // reset iters
while ( iters < maxIters ){
while ( iters < maxIters )
{
// Preconditioned Lanczos
/* Note that there are 4 variants on the Lanczos algorithm. These are
* described in Paige, C. C. (1972). Computational variants of
@ -81,17 +90,17 @@ namespace Eigen {
* A. Greenbaum, Iterative Methods for Solving Linear Systems, SIAM (1987).
*/
const RealScalar beta(beta_new);
// v_old = v; // update: at first time step, this makes v_old = 0 so value of beta doesn't matter
const VectorType v_old(v); // NOT SURE IF CREATING v_old EVERY ITERATION IS EFFICIENT
v_old = v; // update: at first time step, this makes v_old = 0 so value of beta doesn't matter
// const VectorType v_old(v); // NOT SURE IF CREATING v_old EVERY ITERATION IS EFFICIENT
v = v_new; // update
// w = w_new; // update
const VectorType w(w_new); // NOT SURE IF CREATING w EVERY ITERATION IS EFFICIENT
w = w_new; // update
// const VectorType w(w_new); // NOT SURE IF CREATING w EVERY ITERATION IS EFFICIENT
v_new.noalias() = mat*w - beta*v_old; // compute v_new
const RealScalar alpha = v_new.dot(w);
v_new -= alpha*v; // overwrite v_new
w_new = precond.solve(v_new); // overwrite w_new
beta_new2 = v_new.dot(w_new); // compute beta_new
eigen_assert(beta_new2 >= 0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE");
eigen_assert(beta_new2 >= 0.0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE");
beta_new = sqrt(beta_new2); // compute beta_new
v_new /= beta_new; // overwrite v_new for next iteration
w_new /= beta_new; // overwrite w_new for next iteration
@ -107,28 +116,34 @@ namespace Eigen {
s=beta_new/r1; // new sine
// Update solution
// p_oold = p_old;
const VectorType p_oold(p_old); // NOT SURE IF CREATING p_oold EVERY ITERATION IS EFFICIENT
p_oold = p_old;
// const VectorType p_oold(p_old); // NOT SURE IF CREATING p_oold EVERY ITERATION IS EFFICIENT
p_old = p;
p.noalias()=(w-r2*p_old-r3*p_oold) /r1; // IS NOALIAS REQUIRED?
x += beta_one*c*eta*p;
/* Update the squared residual. Note that this is the estimated residual.
The real residual |Ax-b|^2 may be slightly larger */
residualNorm2 *= s*s;
if ( residualNorm2 < threshold2){
if ( residualNorm2 < threshold2)
{
break;
}
eta=-s*eta; // update eta
iters++; // increment iteration number (for output purposes)
}
tol_error = std::sqrt(residualNorm2 / rhsNorm2); // return error. Note that this is the estimated error. The real error |Ax-b|/|b| may be slightly larger
/* Compute error. Note that this is the estimated error. The real
error |Ax-b|/|b| may be slightly larger */
tol_error = std::sqrt(residualNorm2 / rhsNorm2);
}
}
template< typename _MatrixType, int _UpLo=Lower,
typename _Preconditioner = IdentityPreconditioner>
// typename _Preconditioner = IdentityPreconditioner<typename _MatrixType::Scalar> > // preconditioner must be positive definite
class MINRES;
namespace internal {

2
unsupported/Eigen/src/Splines/Spline.h
View File

@ -57,7 +57,7 @@ namespace Eigen
**/
Spline()
: m_knots(1, (Degree==Dynamic ? 2 : 2*Degree+2))
, m_ctrls(ControlPointVectorType::Zero(2,(Degree==Dynamic ? 1 : Degree+1)))
, m_ctrls(ControlPointVectorType::Zero(Dimension,(Degree==Dynamic ? 1 : Degree+1)))
{
// in theory this code can go to the initializer list but it will get pretty
// much unreadable ...

30
unsupported/test/minres.cpp
View File

@ -1,8 +1,8 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2011 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2012 Giacomo Po <gpo@ucla.edu>
// Copyright (C) 2011 Gael Guennebaud <g.gael@free.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
@ -14,19 +14,29 @@
template<typename T> void test_minres_T()
{
MINRES<SparseMatrix<T>, Lower, DiagonalPreconditioner<T> > minres_colmajor_diag;
MINRES<SparseMatrix<T>, Lower, IdentityPreconditioner > minres_colmajor_I;
// MINRES<SparseMatrix<T>, Lower, IncompleteLUT<T> > minres_colmajor_ilut;
//minres<SparseMatrix<T>, SSORPreconditioner<T> > minres_colmajor_ssor;
// Identity preconditioner
MINRES<SparseMatrix<T>, Lower, IdentityPreconditioner > minres_colmajor_lower_I;
MINRES<SparseMatrix<T>, Upper, IdentityPreconditioner > minres_colmajor_upper_I;
// Diagonal preconditioner
MINRES<SparseMatrix<T>, Lower, DiagonalPreconditioner<T> > minres_colmajor_lower_diag;
MINRES<SparseMatrix<T>, Upper, DiagonalPreconditioner<T> > minres_colmajor_upper_diag;
// call tests for SPD matrix
CALL_SUBTEST( check_sparse_spd_solving(minres_colmajor_lower_I) );
CALL_SUBTEST( check_sparse_spd_solving(minres_colmajor_upper_I) );
CALL_SUBTEST( check_sparse_spd_solving(minres_colmajor_lower_diag) );
CALL_SUBTEST( check_sparse_spd_solving(minres_colmajor_upper_diag) );
// TO DO: symmetric semi-definite matrix
// TO DO: symmetric indefinite matrix
CALL_SUBTEST( check_sparse_square_solving(minres_colmajor_diag) );
CALL_SUBTEST( check_sparse_spd_solving(minres_colmajor_I) );
// CALL_SUBTEST( check_sparse_square_solving(minres_colmajor_ilut) );
//CALL_SUBTEST( check_sparse_square_solving(minres_colmajor_ssor) );
}
void test_minres()
{
CALL_SUBTEST_1(test_minres_T<double>());
// CALL_SUBTEST_2(test_minres_T<std::complex<double> >());
// CALL_SUBTEST_2(test_minres_T<std::compex<double> >());
}