GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-12 03:39:01 +08:00
Code Issues Packages Projects Releases Wiki Activity

merge

Browse Source
This commit is contained in:
Benoit Jacob 2009-10-15 16:09:43 -04:00
commit 3c4a025a54
83 changed files with 2390 additions and 768 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
.hgignore
View File

@ -20,4 +20,5 @@ CMakeCache.txt
tags
.*.swp
activity.png
gmon.out
*.out
*.php*

59
CMakeLists.txt
View File

@ -33,20 +33,11 @@ include(CheckCXXCompilerFlag)
set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)
option(EIGEN_BUILD_TESTS "Build tests" OFF)
option(EIGEN_BUILD_DEMOS "Build demos" OFF)
if(NOT WIN32)
option(EIGEN_BUILD_LIB "Build the binary shared library" OFF)
endif(NOT WIN32)
option(EIGEN_BUILD_BTL "Build benchmark suite" OFF)
if(NOT WIN32)
option(EIGEN_BUILD_PKGCONFIG "Build pkg-config .pc file for Eigen" ON)
endif(NOT WIN32)
if(EIGEN_BUILD_LIB)
option(EIGEN_TEST_LIB "Build the unit tests using the library (disable -pedantic)" OFF)
endif(EIGEN_BUILD_LIB)
set(CMAKE_INCLUDE_CURRENT_DIR ON)
if(CMAKE_COMPILER_IS_GNUCXX)
@ -95,7 +86,7 @@ if(MSVC)
option(EIGEN_TEST_SSE2 "Enable/Disable SSE2 in tests/examples" OFF)
if(EIGEN_TEST_SSE2)
if(NOT CMAKE_CL_64)
# arch is not supported on 64 bit systems, SSE is enabled automatically.
# arch is not supported on 64 bit systems, SSE is enabled automatically.
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /arch:SSE2")
endif(NOT CMAKE_CL_64)
message("Enabling SSE2 in tests/examples")
@ -108,6 +99,10 @@ if(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION)
message("Disabling vectorization in tests/examples")
endif(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION)
option(EIGEN_TEST_C++0x "Enables all C++0x features." OFF)
option(EIGEN_TEST_RVALUE_REF_SUPPORT "Enable rvalue references for unit tests." OFF)
include_directories(${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR})
set(INCLUDE_INSTALL_DIR
@ -125,23 +120,43 @@ endif(EIGEN_BUILD_PKGCONFIG)
add_subdirectory(Eigen)
add_subdirectory(doc)
add_subdirectory(doc EXCLUDE_FROM_ALL)
if(EIGEN_BUILD_TESTS)
include(CTest)
add_subdirectory(test)
endif(EIGEN_BUILD_TESTS)
include(CTest)
enable_testing() # must be called from the root CMakeLists, see man page
if(EIGEN_CMAKE_RUN_FROM_CTEST)
add_subdirectory(test) # can't do EXCLUDE_FROM_ALL here, breaks CTest
else(EIGEN_CMAKE_RUN_FROM_CTEST)
add_subdirectory(test EXCLUDE_FROM_ALL)
endif(EIGEN_CMAKE_RUN_FROM_CTEST)
add_subdirectory(unsupported)
if(EIGEN_BUILD_DEMOS)
add_subdirectory(demos)
endif(EIGEN_BUILD_DEMOS)
add_subdirectory(demos EXCLUDE_FROM_ALL)
add_subdirectory(blas EXCLUDE_FROM_ALL)
# TODO: consider also replacing EIGEN_BUILD_BTL by a custom target "make btl"?
if(EIGEN_BUILD_BTL)
add_subdirectory(bench/btl)
add_subdirectory(bench/btl EXCLUDE_FROM_ALL)
endif(EIGEN_BUILD_BTL)
if(EIGEN_BUILD_TESTS)
ei_testing_print_summary()
endif(EIGEN_BUILD_TESTS)
ei_testing_print_summary()
if(NOT MSVC_IDE)
message("")
message("Configured Eigen ${EIGEN_VERSION_NUMBER}")
message("You can now do the following:")
message("--------------+----------------------------------------------------------------")
message("Command | Description")
message("--------------+----------------------------------------------------------------")
message("make install | Install to ${CMAKE_INSTALL_PREFIX}")
message(" | * To change that: cmake . -DCMAKE_INSTALL_PREFIX=yourpath")
message("make btest | Build the unit tests")
message(" | * That takes lots of memory! Easy on the -j option")
message("make test | Build and run the unit tests (using CTest)")
message("make test_qr | Build a specific test, here test_qr. To run it: test/test_qr")
message("make debug_qr | Build a test with full debug info. To run it: test/debug_qr")
message("make blas | Build BLAS library (not the same thing as Eigen)")
message("make doc | Generate the API documentation, requires Doxygen & LaTeX")
message("--------------+----------------------------------------------------------------")
endif(NOT MSVC_IDE)

2
Doxyfile
View File

@ -5,7 +5,7 @@
#---------------------------------------------------------------------------
DOXYFILE_ENCODING = UTF-8
PROJECT_NAME = Eigen
PROJECT_NUMBER = 2.0
PROJECT_NUMBER = you-got-it-wrong
OUTPUT_DIRECTORY = ./
CREATE_SUBDIRS = NO
OUTPUT_LANGUAGE = English

26
Eigen/src/Array/Replicate.h
View File

@ -45,14 +45,10 @@ struct ei_traits<Replicate<MatrixType,RowFactor,ColFactor> >
typedef typename ei_nested<MatrixType>::type MatrixTypeNested;
typedef typename ei_unref<MatrixTypeNested>::type _MatrixTypeNested;
enum {
RowsPlusOne = (MatrixType::RowsAtCompileTime != Dynamic) ?
int(MatrixType::RowsAtCompileTime) + 1 : Dynamic,
ColsPlusOne = (MatrixType::ColsAtCompileTime != Dynamic) ?
int(MatrixType::ColsAtCompileTime) + 1 : Dynamic,
RowsAtCompileTime = RowFactor==Dynamic || MatrixType::RowsAtCompileTime==Dynamic
RowsAtCompileTime = RowFactor==Dynamic || int(MatrixType::RowsAtCompileTime)==Dynamic
? Dynamic
: RowFactor * MatrixType::RowsAtCompileTime,
ColsAtCompileTime = ColFactor==Dynamic || MatrixType::ColsAtCompileTime==Dynamic
ColsAtCompileTime = ColFactor==Dynamic || int(MatrixType::ColsAtCompileTime)==Dynamic
? Dynamic
: ColFactor * MatrixType::ColsAtCompileTime,
MaxRowsAtCompileTime = RowsAtCompileTime,
@ -69,15 +65,22 @@ template<typename MatrixType,int RowFactor,int ColFactor> class Replicate
EIGEN_GENERIC_PUBLIC_INTERFACE(Replicate)
inline Replicate(const MatrixType& matrix)
template<typename OriginalMatrixType>
inline explicit Replicate(const OriginalMatrixType& matrix)
: m_matrix(matrix), m_rowFactor(RowFactor), m_colFactor(ColFactor)
{
EIGEN_STATIC_ASSERT((ei_is_same_type<MatrixType,OriginalMatrixType>::ret),
THE_MATRIX_OR_EXPRESSION_THAT_YOU_PASSED_DOES_NOT_HAVE_THE_EXPECTED_TYPE)
ei_assert(RowFactor!=Dynamic && ColFactor!=Dynamic);
}
inline Replicate(const MatrixType& matrix, int rowFactor, int colFactor)
template<typename OriginalMatrixType>
inline Replicate(const OriginalMatrixType& matrix, int rowFactor, int colFactor)
: m_matrix(matrix), m_rowFactor(rowFactor), m_colFactor(colFactor)
{}
{
EIGEN_STATIC_ASSERT((ei_is_same_type<MatrixType,OriginalMatrixType>::ret),
THE_MATRIX_OR_EXPRESSION_THAT_YOU_PASSED_DOES_NOT_HAVE_THE_EXPECTED_TYPE)
}
inline int rows() const { return m_matrix.rows() * m_rowFactor.value(); }
inline int cols() const { return m_matrix.cols() * m_colFactor.value(); }
@ -91,6 +94,9 @@ template<typename MatrixType,int RowFactor,int ColFactor> class Replicate
const typename MatrixType::Nested m_matrix;
const ei_int_if_dynamic<RowFactor> m_rowFactor;
const ei_int_if_dynamic<ColFactor> m_colFactor;
private:
Replicate& operator=(const Replicate&);
};
/** \nonstableyet
@ -106,7 +112,7 @@ template<int RowFactor, int ColFactor>
inline const Replicate<Derived,RowFactor,ColFactor>
MatrixBase<Derived>::replicate() const
{
return derived();
return Replicate<Derived,RowFactor,ColFactor>(derived());
}
/** \nonstableyet

11
Eigen/src/Array/VectorwiseOp.h
View File

@ -95,6 +95,14 @@ class PartialReduxExpr : ei_no_assignment_operator,
return m_functor(m_matrix.row(i));
}
const Scalar coeff(int index) const
{
if (Direction==Vertical)
return m_functor(m_matrix.col(index));
else
return m_functor(m_matrix.row(index));
}
protected:
const MatrixTypeNested m_matrix;
const MemberOp m_functor;
@ -442,6 +450,9 @@ template<typename ExpressionType, int Direction> class VectorwiseOp
protected:
ExpressionTypeNested m_matrix;
private:
VectorwiseOp& operator=(const VectorwiseOp&);
};
/** \array_module

9
Eigen/src/Core/CMakeLists.txt
View File

@ -5,13 +5,6 @@ INSTALL(FILES
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Core COMPONENT Devel
)
FILE(GLOB Eigen_Core_Product_SRCS "products/*.h")
INSTALL(FILES
${Eigen_Core_Product_SRCS}
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Core/products COMPONENT Devel
)
ADD_SUBDIRECTORY(products)
ADD_SUBDIRECTORY(util)
ADD_SUBDIRECTORY(arch)

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

@ -116,6 +116,9 @@ struct CommaInitializer
int m_row; // current row id
int m_col; // current col id
int m_currentBlockRows; // current block height
private:
CommaInitializer& operator=(const CommaInitializer&);
};
/** \anchor MatrixBaseCommaInitRef

3
Eigen/src/Core/Cwise.h
View File

@ -178,6 +178,9 @@ template<typename ExpressionType> class Cwise
protected:
ExpressionTypeNested m_matrix;
private:
Cwise& operator=(const Cwise&);
};
/** \returns a Cwise wrapper of *this providing additional coefficient-wise operations

1
Eigen/src/Core/CwiseNullaryOp.h
View File

@ -147,7 +147,6 @@ EIGEN_STRONG_INLINE const CwiseNullaryOp<CustomNullaryOp, Derived>
MatrixBase<Derived>::NullaryExpr(int size, const CustomNullaryOp& func)
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
ei_assert(IsVectorAtCompileTime);
if(RowsAtCompileTime == 1) return CwiseNullaryOp<CustomNullaryOp, Derived>(1, size, func);
else return CwiseNullaryOp<CustomNullaryOp, Derived>(size, 1, func);
}

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

@ -351,6 +351,8 @@ struct ei_scalar_multiple_op {
EIGEN_STRONG_INLINE const PacketScalar packetOp(const PacketScalar& a) const
{ return ei_pmul(a, ei_pset1(m_other)); }
const Scalar m_other;
private:
ei_scalar_multiple_op& operator=(const ei_scalar_multiple_op&);
};
template<typename Scalar>
struct ei_functor_traits<ei_scalar_multiple_op<Scalar> >
@ -378,6 +380,8 @@ struct ei_scalar_quotient1_impl {
EIGEN_STRONG_INLINE const PacketScalar packetOp(const PacketScalar& a) const
{ return ei_pmul(a, ei_pset1(m_other)); }
const Scalar m_other;
private:
ei_scalar_quotient1_impl& operator=(const ei_scalar_quotient1_impl&);
};
template<typename Scalar>
struct ei_functor_traits<ei_scalar_quotient1_impl<Scalar,true> >
@ -423,6 +427,8 @@ struct ei_scalar_constant_op {
EIGEN_STRONG_INLINE const Scalar operator() (int, int = 0) const { return m_other; }
EIGEN_STRONG_INLINE const PacketScalar packetOp() const { return ei_pset1(m_other); }
const Scalar m_other;
private:
ei_scalar_constant_op& operator=(const ei_scalar_constant_op&);
};
template<typename Scalar>
struct ei_functor_traits<ei_scalar_constant_op<Scalar> >

30
Eigen/src/Core/MathFunctions.h
View File

@ -317,4 +317,34 @@ inline bool ei_isApproxOrLessThan(long double a, long double b, long double prec
return a <= b || ei_isApprox(a, b, prec);
}
/**************
*** bool ***
**************/
template<> inline bool precision<bool>() { return 0; }
inline bool ei_real(bool x) { return x; }
inline bool& ei_real_ref(bool& x) { return x; }
inline bool ei_imag(bool) { return 0; }
inline bool ei_conj(bool x) { return x; }
inline bool ei_abs(bool x) { return x; }
inline bool ei_abs2(bool x) { return x; }
inline bool ei_sqrt(bool x) { return x; }
template<> inline bool ei_random()
{
return (ei_random<int>(0,1) == 1);
}
inline bool ei_isMuchSmallerThan(bool a, bool, bool = precision<bool>())
{
return !a;
}
inline bool ei_isApprox(bool a, bool b, bool = precision<bool>())
{
return a == b;
}
inline bool ei_isApproxOrLessThan(bool a, bool b, bool = precision<bool>())
{
return int(a) <= int(b);
}
#endif // EIGEN_MATHFUNCTIONS_H

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

@ -538,7 +538,7 @@ class Matrix
* data pointers.
*/
template<typename OtherDerived>
void swap(const MatrixBase<OtherDerived>& other);
void swap(MatrixBase<OtherDerived> EIGEN_REF_TO_TEMPORARY other);
/** \name Map
* These are convenience functions returning Map objects. The Map() static functions return unaligned Map objects,
@ -707,6 +707,8 @@ struct ei_conservative_resize_like_impl
{
static void run(MatrixBase<Derived>& _this, const MatrixBase<OtherDerived>& other)
{
if (_this.rows() == other.rows() && _this.cols() == other.cols()) return;
// Note: Here is space for improvement. Basically, for conservativeResize(int,int),
// neither RowsAtCompileTime or ColsAtCompileTime must be Dynamic. If only one of the
// dimensions is dynamic, one could use either conservativeResize(int rows, NoChange_t) or
@ -728,6 +730,8 @@ struct ei_conservative_resize_like_impl<Derived,OtherDerived,true>
{
static void run(MatrixBase<Derived>& _this, const MatrixBase<OtherDerived>& other)
{
if (_this.rows() == other.rows() && _this.cols() == other.cols()) return;
// segment(...) will check whether Derived/OtherDerived are vectors!
typename MatrixBase<Derived>::PlainMatrixType tmp(other);
const int common_size = std::min<int>(_this.size(),tmp.size());
@ -756,7 +760,7 @@ struct ei_matrix_swap_impl<MatrixType, OtherDerived, true>
template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
template<typename OtherDerived>
inline void Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::swap(const MatrixBase<OtherDerived>& other)
inline void Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::swap(MatrixBase<OtherDerived> EIGEN_REF_TO_TEMPORARY other)
{
enum { SwapPointers = ei_is_same_type<Matrix, OtherDerived>::ret && Base::SizeAtCompileTime==Dynamic };
ei_matrix_swap_impl<Matrix, OtherDerived, bool(SwapPointers)>::run(*this, *const_cast<MatrixBase<OtherDerived>*>(&other));

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

@ -592,7 +592,7 @@ template<typename Derived> class MatrixBase
{ return typename ei_eval<Derived>::type(derived()); }
template<typename OtherDerived>
void swap(const MatrixBase<OtherDerived>& other);
void swap(MatrixBase<OtherDerived> EIGEN_REF_TO_TEMPORARY other);
template<unsigned int Added>
const Flagged<Derived, Added, 0> marked() const;
@ -702,8 +702,10 @@ template<typename Derived> class MatrixBase
const LU<PlainMatrixType> lu() const;
const PartialLU<PlainMatrixType> partialLu() const;
const PlainMatrixType inverse() const;
void computeInverse(PlainMatrixType *result) const;
bool computeInverseWithCheck( PlainMatrixType *result ) const;
template<typename ResultType>
void computeInverse(ResultType *result) const;
template<typename ResultType>
bool computeInverseWithCheck(ResultType *result ) const;
Scalar determinant() const;
/////////// Cholesky module ///////////

42
Eigen/src/Core/MatrixStorage.h
View File

@ -29,32 +29,38 @@
struct ei_constructor_without_unaligned_array_assert {};
/** \internal
* Static array automatically aligned if the total byte size is a multiple of 16 and the matrix options require auto alignment
* Static array. If the MatrixOptions require auto-alignment, the array will be automatically aligned:
* to 16 bytes boundary if the total size is a multiple of 16 bytes.
*/
template <typename T, int Size, int MatrixOptions,
bool Align = (!(MatrixOptions&DontAlign)) && (((Size*sizeof(T))&0xf)==0)
> struct ei_matrix_array
{
EIGEN_ALIGN_128 T array[Size];
ei_matrix_array()
{
#ifndef EIGEN_DISABLE_UNALIGNED_ARRAY_ASSERT
ei_assert((reinterpret_cast<size_t>(array) & 0xf) == 0
&& "this assertion is explained here: http://eigen.tuxfamily.org/dox/UnalignedArrayAssert.html **** READ THIS WEB PAGE !!! ****");
#endif
}
ei_matrix_array(ei_constructor_without_unaligned_array_assert) {}
};
template <typename T, int Size, int MatrixOptions> struct ei_matrix_array<T,Size,MatrixOptions,false>
int Alignment = (MatrixOptions&DontAlign) ? 0
: (((Size*sizeof(T))%16)==0) ? 16
: 0 >
struct ei_matrix_array
{
T array[Size];
ei_matrix_array() {}
ei_matrix_array(ei_constructor_without_unaligned_array_assert) {}
};
#ifdef EIGEN_DISABLE_UNALIGNED_ARRAY_ASSERT
#define EIGEN_MAKE_UNALIGNED_ARRAY_ASSERT(sizemask)
#else
#define EIGEN_MAKE_UNALIGNED_ARRAY_ASSERT(sizemask) \
ei_assert((reinterpret_cast<size_t>(array) & sizemask) == 0 \
&& "this assertion is explained here: " \
"http://eigen.tuxfamily.org/dox/UnalignedArrayAssert.html" \
" **** READ THIS WEB PAGE !!! ****");
#endif
template <typename T, int Size, int MatrixOptions>
struct ei_matrix_array<T, Size, MatrixOptions, 16>
{
EIGEN_ALIGN16 T array[Size];
ei_matrix_array() { EIGEN_MAKE_UNALIGNED_ARRAY_ASSERT(0xf) }
ei_matrix_array(ei_constructor_without_unaligned_array_assert) {}
};
/** \internal
*
* \class ei_matrix_storage

3
Eigen/src/Core/NestByValue.h
View File

@ -102,6 +102,9 @@ template<typename ExpressionType> class NestByValue
protected:
const ExpressionType m_expression;
private:
NestByValue& operator=(const NestByValue&);
};
/** \returns an expression of the temporary version of *this.

3
Eigen/src/Core/NoAlias.h
View File

@ -73,6 +73,9 @@ class NoAlias
protected:
ExpressionType& m_expression;
private:
NoAlias& operator=(const NoAlias&);
};
/** \returns a pseudo expression of \c *this with an operator= assuming

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

@ -135,7 +135,7 @@ struct ProductReturnType<Lhs,Rhs,UnrolledProduct>
{
typedef typename ei_nested<Lhs,Rhs::ColsAtCompileTime>::type LhsNested;
typedef typename ei_nested<Rhs,Lhs::RowsAtCompileTime>::type RhsNested;
typedef GeneralProduct<Lhs, Rhs, UnrolledProduct> Type;
typedef GeneralProduct<LhsNested, RhsNested, UnrolledProduct> Type;
};
@ -211,6 +211,9 @@ class GeneralProduct<Lhs, Rhs, OuterProduct>
{
ei_outer_product_selector<(int(Dest::Flags)&RowMajorBit) ? RowMajor : ColMajor>::run(*this, dest, alpha);
}
private:
GeneralProduct& operator=(const GeneralProduct&);
};
template<> struct ei_outer_product_selector<ColMajor> {
@ -276,6 +279,9 @@ class GeneralProduct<Lhs, Rhs, GemvProduct>
ei_gemv_selector<Side,(int(MatrixType::Flags)&RowMajorBit) ? RowMajor : ColMajor,
bool(ei_blas_traits<MatrixType>::ActualAccess)>::run(*this, dst, alpha);
}
private:
GeneralProduct& operator=(const GeneralProduct&);
};
// The vector is on the left => transposition

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

@ -137,6 +137,8 @@ class ProductBase : public MatrixBase<Derived>
void coeffRef(int,int);
void coeff(int) const;
void coeffRef(int);
ProductBase& operator=(const ProductBase&);
};
template<typename NestedProduct>

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

@ -112,6 +112,16 @@ struct ei_redux_novec_unroller<Func, Derived, Start, 1>
}
};
// This is actually dead code and will never be called. It is required
// to prevent false warnings regarding failed inlining though
// for 0 length run() will never be called at all.
template<typename Func, typename Derived, int Start>
struct ei_redux_novec_unroller<Func, Derived, Start, 0>
{
typedef typename Derived::Scalar Scalar;
EIGEN_STRONG_INLINE static Scalar run(const Derived&, const Func&) { return Scalar(); }
};
/*** vectorization ***/
template<typename Func, typename Derived, int Start, int Length>
@ -297,7 +307,7 @@ struct ei_redux_impl<Func, Derived, LinearVectorization, CompleteUnrolling>
/** \returns the result of a full redux operation on the whole matrix or vector using \a func
*
* The template parameter \a BinaryOp is the type of the functor \a func which must be
* an assiociative operator. Both current STL and TR1 functor styles are handled.
* an associative operator. Both current STL and TR1 functor styles are handled.
*
* \sa MatrixBase::sum(), MatrixBase::minCoeff(), MatrixBase::maxCoeff(), MatrixBase::colwise(), MatrixBase::rowwise()
*/

14
Eigen/src/Core/StableNorm.h
View File

@ -108,21 +108,15 @@ MatrixBase<Derived>::blueNorm() const
iemax = std::numeric_limits<RealScalar>::max_exponent; // maximum exponent
rbig = std::numeric_limits<RealScalar>::max(); // largest floating-point number
// Check the basic machine-dependent constants.
if(iemin > 1 - 2*it || 1+it>iemax || (it==2 && ibeta<5)
|| (it<=4 && ibeta <= 3 ) || it<2)
{
ei_assert(false && "the algorithm cannot be guaranteed on this computer");
}
iexp = -((1-iemin)/2);
b1 = RealScalar(std::pow(double(ibeta),iexp)); // lower boundary of midrange
b1 = RealScalar(std::pow(RealScalar(ibeta),RealScalar(iexp))); // lower boundary of midrange
iexp = (iemax + 1 - it)/2;
b2 = RealScalar(std::pow(double(ibeta),iexp)); // upper boundary of midrange
b2 = RealScalar(std::pow(RealScalar(ibeta),RealScalar(iexp))); // upper boundary of midrange
iexp = (2-iemin)/2;
s1m = RealScalar(std::pow(double(ibeta),iexp)); // scaling factor for lower range
s1m = RealScalar(std::pow(RealScalar(ibeta),RealScalar(iexp))); // scaling factor for lower range
iexp = - ((iemax+it)/2);
s2m = RealScalar(std::pow(double(ibeta),iexp)); // scaling factor for upper range
s2m = RealScalar(std::pow(RealScalar(ibeta),RealScalar(iexp))); // scaling factor for upper range
overfl = rbig*s2m; // overfow boundary for abig
eps = RealScalar(std::pow(double(ibeta), 1-it));

11
Eigen/src/Core/Swap.h
View File

@ -117,6 +117,9 @@ template<typename ExpressionType> class SwapWrapper
protected:
ExpressionType& m_expression;
private:
SwapWrapper& operator=(const SwapWrapper&);
};
/** swaps *this with the expression \a other.
@ -128,15 +131,9 @@ template<typename ExpressionType> class SwapWrapper
*/
template<typename Derived>
template<typename OtherDerived>
void MatrixBase<Derived>::swap(const MatrixBase<OtherDerived>& other)
void MatrixBase<Derived>::swap(MatrixBase<OtherDerived> EIGEN_REF_TO_TEMPORARY other)
{
(SwapWrapper<Derived>(derived())).lazyAssign(other);
}
#endif // EIGEN_SWAP_H

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

@ -300,13 +300,13 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView
}
template<typename OtherDerived>
void swap(const TriangularBase<OtherDerived>& other)
void swap(TriangularBase<OtherDerived> EIGEN_REF_TO_TEMPORARY other)
{
TriangularView<SwapWrapper<MatrixType>,Mode>(const_cast<MatrixType&>(m_matrix)).lazyAssign(other.derived());
}
template<typename OtherDerived>
void swap(const MatrixBase<OtherDerived>& other)
void swap(MatrixBase<OtherDerived> EIGEN_REF_TO_TEMPORARY other)
{
TriangularView<SwapWrapper<MatrixType>,Mode>(const_cast<MatrixType&>(m_matrix)).lazyAssign(other.derived());
}

6
Eigen/src/Core/arch/AltiVec/PacketMath.h
View File

@ -265,14 +265,14 @@ template<> inline void ei_pstoreu(int* to , const v4i& from )
template<> inline float ei_pfirst(const v4f& a)
{
float EIGEN_ALIGN_128 af[4];
float EIGEN_ALIGN16 af[4];
vec_st(a, 0, af);
return af[0];
}
template<> inline int ei_pfirst(const v4i& a)
{
int EIGEN_ALIGN_128 ai[4];
int EIGEN_ALIGN16 ai[4];
vec_st(a, 0, ai);
return ai[0];
}
@ -373,7 +373,7 @@ inline float ei_predux_mul(const v4f& a)
inline int ei_predux_mul(const v4i& a)
{
EIGEN_ALIGN_128 int aux[4];
EIGEN_ALIGN16 int aux[4];
ei_pstore(aux, a);
return aux[0] * aux[1] * aux[2] * aux[3];
}

6
Eigen/src/Core/arch/SSE/PacketMath.h
View File

@ -359,7 +359,7 @@ template<> EIGEN_STRONG_INLINE int ei_predux_mul<Packet4i>(const Packet4i& a)
// after some experiments, it is seems this is the fastest way to implement it
// for GCC (eg., reusing ei_pmul is very slow !)
// TODO try to call _mm_mul_epu32 directly
EIGEN_ALIGN_128 int aux[4];
EIGEN_ALIGN16 int aux[4];
ei_pstore(aux, a);
return (aux[0] * aux[1]) * (aux[2] * aux[3]);;
}
@ -378,7 +378,7 @@ template<> EIGEN_STRONG_INLINE int ei_predux_min<Packet4i>(const Packet4i& a)
{
// after some experiments, it is seems this is the fastest way to implement it
// for GCC (eg., it does not like using std::min after the ei_pstore !!)
EIGEN_ALIGN_128 int aux[4];
EIGEN_ALIGN16 int aux[4];
ei_pstore(aux, a);
register int aux0 = aux[0]<aux[1] ? aux[0] : aux[1];
register int aux2 = aux[2]<aux[3] ? aux[2] : aux[3];
@ -399,7 +399,7 @@ template<> EIGEN_STRONG_INLINE int ei_predux_max<Packet4i>(const Packet4i& a)
{
// after some experiments, it is seems this is the fastest way to implement it
// for GCC (eg., it does not like using std::min after the ei_pstore !!)
EIGEN_ALIGN_128 int aux[4];
EIGEN_ALIGN16 int aux[4];
ei_pstore(aux, a);
register int aux0 = aux[0]>aux[1] ? aux[0] : aux[1];
register int aux2 = aux[2]>aux[3] ? aux[2] : aux[3];

6
Eigen/src/Core/products/CMakeLists.txt Normal file
View File

@ -0,0 +1,6 @@
FILE(GLOB Eigen_Core_Product_SRCS "*.h")
INSTALL(FILES
${Eigen_Core_Product_SRCS}
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Core/products COMPONENT Devel
)

35
Eigen/src/Core/products/GeneralMatrixVector.h
View File

@ -57,8 +57,7 @@ void ei_cache_friendly_product_colmajor_times_vector(
if(ConjugateRhs)
alpha = ei_conj(alpha);
// std::cerr << "prod " << size << " " << rhs.size() << "\n";
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef typename ei_packet_traits<Scalar>::type Packet;
const int PacketSize = sizeof(Packet)/sizeof(Scalar);
@ -69,9 +68,9 @@ void ei_cache_friendly_product_colmajor_times_vector(
const int PeelAlignedMask = PacketSize*peels-1;
// How many coeffs of the result do we have to skip to be aligned.
// Here we assume data are at least aligned on the base scalar type that is mandatory anyway.
const int alignedStart = ei_alignmentOffset(res,size);
const int alignedSize = PacketSize>1 ? alignedStart + ((size-alignedStart) & ~PacketAlignedMask) : 0;
// Here we assume data are at least aligned on the base scalar type.
int alignedStart = ei_alignmentOffset(res,size);
int alignedSize = PacketSize>1 ? alignedStart + ((size-alignedStart) & ~PacketAlignedMask) : 0;
const int peeledSize = peels>1 ? alignedStart + ((alignedSize-alignedStart) & ~PeelAlignedMask) : alignedStart;
const int alignmentStep = PacketSize>1 ? (PacketSize - lhsStride % PacketSize) & PacketAlignedMask : 0;
@ -84,12 +83,18 @@ void ei_cache_friendly_product_colmajor_times_vector(
// find how many columns do we have to skip to be aligned with the result (if possible)
int skipColumns = 0;
if (PacketSize>1)
// if the data cannot be aligned (TODO add some compile time tests when possible, e.g. for floats)
if( (size_t(lhs)%sizeof(RealScalar)) || (size_t(res)%sizeof(RealScalar)) )
{
alignedSize = 0;
alignedStart = 0;
}
else if (PacketSize>1)
{
ei_internal_assert(size_t(lhs+lhsAlignmentOffset)%sizeof(Packet)==0 || size<PacketSize);
while (skipColumns<PacketSize &&
alignedStart != ((lhsAlignmentOffset + alignmentStep*skipColumns)%PacketSize))
alignedStart != ((lhsAlignmentOffset + alignmentStep*skipColumns)%PacketSize))
++skipColumns;
if (skipColumns==PacketSize)
{
@ -263,6 +268,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
ei_conj_helper<ConjugateLhs,ConjugateRhs> cj;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef typename ei_packet_traits<Scalar>::type Packet;
const int PacketSize = sizeof(Packet)/sizeof(Scalar);
@ -274,9 +280,10 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
const int size = rhsSize;
// How many coeffs of the result do we have to skip to be aligned.
// Here we assume data are at least aligned on the base scalar type that is mandatory anyway.
const int alignedStart = ei_alignmentOffset(rhs, size);
const int alignedSize = PacketSize>1 ? alignedStart + ((size-alignedStart) & ~PacketAlignedMask) : 0;
// Here we assume data are at least aligned on the base scalar type
// if that's not the case then vectorization is discarded, see below.
int alignedStart = ei_alignmentOffset(rhs, size);
int alignedSize = PacketSize>1 ? alignedStart + ((size-alignedStart) & ~PacketAlignedMask) : 0;
const int peeledSize = peels>1 ? alignedStart + ((alignedSize-alignedStart) & ~PeelAlignedMask) : alignedStart;
const int alignmentStep = PacketSize>1 ? (PacketSize - lhsStride % PacketSize) & PacketAlignedMask : 0;
@ -289,7 +296,13 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
// find how many rows do we have to skip to be aligned with rhs (if possible)
int skipRows = 0;
if (PacketSize>1)
// if the data cannot be aligned (TODO add some compile time tests when possible, e.g. for floats)
if( (size_t(lhs)%sizeof(RealScalar)) || (size_t(rhs)%sizeof(RealScalar)) )
{
alignedSize = 0;
alignedStart = 0;
}
else if (PacketSize>1)
{
ei_internal_assert(size_t(lhs+lhsAlignmentOffset)%sizeof(Packet)==0 || size<PacketSize);

2
Eigen/src/Core/products/SelfadjointRank2Update.h
View File

@ -38,10 +38,8 @@ struct ei_selfadjoint_rank2_update_selector<Scalar,UType,VType,LowerTriangular>
static void run(Scalar* mat, int stride, const UType& u, const VType& v, Scalar alpha)
{
const int size = u.size();
// std::cerr << "lower \n" << u.transpose() << "\n" << v.transpose() << "\n\n";
for (int i=0; i<size; ++i)
{
// std::cerr <<
Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i+i, size-i) +=
(alpha * ei_conj(u.coeff(i))) * v.end(size-i)
+ (alpha * ei_conj(v.coeff(i))) * u.end(size-i);

26
Eigen/src/Core/util/Macros.h
View File

@ -29,8 +29,8 @@
#undef minor
#define EIGEN_WORLD_VERSION 2
#define EIGEN_MAJOR_VERSION 0
#define EIGEN_MINOR_VERSION 52
#define EIGEN_MAJOR_VERSION 90
#define EIGEN_MINOR_VERSION 0
#define EIGEN_VERSION_AT_LEAST(x,y,z) (EIGEN_WORLD_VERSION>x || (EIGEN_WORLD_VERSION>=x && \
(EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \
@ -202,25 +202,28 @@ using Eigen::ei_cos;
#define EIGEN_ASM_COMMENT(X)
#endif
/* EIGEN_ALIGN_128 forces data to be 16-byte aligned, EVEN if vectorization (EIGEN_VECTORIZE) is disabled,
/* EIGEN_ALIGN_TO_BOUNDARY(n) forces data to be n-byte aligned. This is used to satisfy SIMD requirements.
* However, we do that EVEN if vectorization (EIGEN_VECTORIZE) is disabled,
* so that vectorization doesn't affect binary compatibility.
*
* If we made alignment depend on whether or not EIGEN_VECTORIZE is defined, it would be impossible to link
* vectorized and non-vectorized code.
*/
#if !EIGEN_ALIGN
#define EIGEN_ALIGN_128
#define EIGEN_ALIGN_TO_BOUNDARY(n)
#elif (defined __GNUC__)
#define EIGEN_ALIGN_128 __attribute__((aligned(16)))
#define EIGEN_ALIGN_TO_BOUNDARY(n) __attribute__((aligned(n)))
#elif (defined _MSC_VER)
#define EIGEN_ALIGN_128 __declspec(align(16))
#define EIGEN_ALIGN_TO_BOUNDARY(n) __declspec(align(n))
#elif (defined __SUNPRO_CC)
// FIXME not sure about this one:
#define EIGEN_ALIGN_128 __attribute__((aligned(16)))
#define EIGEN_ALIGN_TO_BOUNDARY(n) __attribute__((aligned(n)))
#else
#error Please tell me what is the equivalent of __attribute__((aligned(16))) for your compiler
#error Please tell me what is the equivalent of __attribute__((aligned(n))) for your compiler
#endif
#define EIGEN_ALIGN16 EIGEN_ALIGN_TO_BOUNDARY(16)
#ifdef EIGEN_DONT_USE_RESTRICT_KEYWORD
#define EIGEN_RESTRICT
#endif
@ -251,6 +254,13 @@ using Eigen::ei_cos;
// needed to define it here as escaping characters in CMake add_definition's argument seems very problematic.
#define EIGEN_DOCS_IO_FORMAT IOFormat(3, 0, " ", "\n", "", "")
// C++0x features
#if defined(__GXX_EXPERIMENTAL_CXX0X__) || (defined(_MSC_VER) && (_MSC_VER >= 1600))
#define EIGEN_REF_TO_TEMPORARY &&
#else
#define EIGEN_REF_TO_TEMPORARY const &
#endif
#ifdef _MSC_VER
#define EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Derived) \
using Base::operator =; \

5
Eigen/src/Core/util/StaticAssert.h
View File

@ -41,7 +41,7 @@
#ifndef EIGEN_NO_STATIC_ASSERT
#ifdef __GXX_EXPERIMENTAL_CXX0X__
#if defined(__GXX_EXPERIMENTAL_CXX0X__) || (defined(_MSC_VER) && (_MSC_VER >= 1600))
// if native static_assert is enabled, let's use it
#define EIGEN_STATIC_ASSERT(X,MSG) static_assert(X,#MSG);
@ -77,7 +77,8 @@
THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES,
INVALID_MATRIX_TEMPLATE_PARAMETERS,
BOTH_MATRICES_MUST_HAVE_THE_SAME_STORAGE_ORDER,
THIS_METHOD_IS_ONLY_FOR_DIAGONAL_MATRIX
THIS_METHOD_IS_ONLY_FOR_DIAGONAL_MATRIX,
THE_MATRIX_OR_EXPRESSION_THAT_YOU_PASSED_DOES_NOT_HAVE_THE_EXPECTED_TYPE
};
};

297
Eigen/src/Eigenvalues/ComplexEigenSolver.h
View File

@ -1,148 +1,149 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Claire Maurice
// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H
#define EIGEN_COMPLEX_EIGEN_SOLVER_H
/** \eigenvalues_module \ingroup Eigenvalues_Module
* \nonstableyet
*
* \class ComplexEigenSolver
*
* \brief Eigen values/vectors solver for general complex matrices
*
* \param MatrixType the type of the matrix of which we are computing the eigen decomposition
*
* \sa class EigenSolver, class SelfAdjointEigenSolver
*/
template<typename _MatrixType> class ComplexEigenSolver
{
public:
typedef _MatrixType MatrixType;
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef std::complex<RealScalar> Complex;
typedef Matrix<Complex, MatrixType::ColsAtCompileTime,1> EigenvalueType;
typedef Matrix<Complex, MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime> EigenvectorType;
/**
* \brief Default Constructor.
*
* The default constructor is useful in cases in which the user intends to
* perform decompositions via ComplexEigenSolver::compute(const MatrixType&).
*/
ComplexEigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false)
{}
ComplexEigenSolver(const MatrixType& matrix)
: m_eivec(matrix.rows(),matrix.cols()),
m_eivalues(matrix.cols()),
m_isInitialized(false)
{
compute(matrix);
}
EigenvectorType eigenvectors(void) const
{
ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
return m_eivec;
}
EigenvalueType eigenvalues() const
{
ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
return m_eivalues;
}
void compute(const MatrixType& matrix);
protected:
MatrixType m_eivec;
EigenvalueType m_eivalues;
bool m_isInitialized;
};
template<typename MatrixType>
void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
{
// this code is inspired from Jampack
assert(matrix.cols() == matrix.rows());
int n = matrix.cols();
m_eivalues.resize(n,1);
RealScalar eps = epsilon<RealScalar>();
// Reduce to complex Schur form
ComplexSchur<MatrixType> schur(matrix);
m_eivalues = schur.matrixT().diagonal();
m_eivec.setZero();
Scalar d2, z;
RealScalar norm = matrix.norm();
// compute the (normalized) eigenvectors
for(int k=n-1 ; k>=0 ; k--)
{
d2 = schur.matrixT().coeff(k,k);
m_eivec.coeffRef(k,k) = Scalar(1.0,0.0);
for(int i=k-1 ; i>=0 ; i--)
{
m_eivec.coeffRef(i,k) = -schur.matrixT().coeff(i,k);
if(k-i-1>0)
m_eivec.coeffRef(i,k) -= (schur.matrixT().row(i).segment(i+1,k-i-1) * m_eivec.col(k).segment(i+1,k-i-1)).value();
z = schur.matrixT().coeff(i,i) - d2;
if(z==Scalar(0))
ei_real_ref(z) = eps * norm;
m_eivec.coeffRef(i,k) = m_eivec.coeff(i,k) / z;
}
m_eivec.col(k).normalize();
}
m_eivec = schur.matrixU() * m_eivec;
m_isInitialized = true;
// sort the eigenvalues
{
for (int i=0; i<n; i++)
{
int k;
m_eivalues.cwise().abs().end(n-i).minCoeff(&k);
if (k != 0)
{
k += i;
std::swap(m_eivalues[k],m_eivalues[i]);
m_eivec.col(i).swap(m_eivec.col(k));
}
}
}
}
#endif // EIGEN_COMPLEX_EIGEN_SOLVER_H
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Claire Maurice
// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H
#define EIGEN_COMPLEX_EIGEN_SOLVER_H
/** \eigenvalues_module \ingroup Eigenvalues_Module
* \nonstableyet
*
* \class ComplexEigenSolver
*
* \brief Eigen values/vectors solver for general complex matrices
*
* \param MatrixType the type of the matrix of which we are computing the eigen decomposition
*
* \sa class EigenSolver, class SelfAdjointEigenSolver
*/
template<typename _MatrixType> class ComplexEigenSolver
{
public:
typedef _MatrixType MatrixType;
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef std::complex<RealScalar> Complex;
typedef Matrix<Complex, MatrixType::ColsAtCompileTime,1> EigenvalueType;
typedef Matrix<Complex, MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime> EigenvectorType;
/**
* \brief Default Constructor.
*
* The default constructor is useful in cases in which the user intends to
* perform decompositions via ComplexEigenSolver::compute(const MatrixType&).
*/
ComplexEigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false)
{}
ComplexEigenSolver(const MatrixType& matrix)
: m_eivec(matrix.rows(),matrix.cols()),
m_eivalues(matrix.cols()),
m_isInitialized(false)
{
compute(matrix);
}
EigenvectorType eigenvectors(void) const
{
ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
return m_eivec;
}
EigenvalueType eigenvalues() const
{
ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
return m_eivalues;
}
void compute(const MatrixType& matrix);
protected:
MatrixType m_eivec;
EigenvalueType m_eivalues;
bool m_isInitialized;
};
template<typename MatrixType>
void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
{
// this code is inspired from Jampack
assert(matrix.cols() == matrix.rows());
int n = matrix.cols();
m_eivalues.resize(n,1);
m_eivec.resize(n,n);
RealScalar eps = epsilon<RealScalar>();
// Reduce to complex Schur form
ComplexSchur<MatrixType> schur(matrix);
m_eivalues = schur.matrixT().diagonal();
m_eivec.setZero();
Scalar d2, z;
RealScalar norm = matrix.norm();
// compute the (normalized) eigenvectors
for(int k=n-1 ; k>=0 ; k--)
{
d2 = schur.matrixT().coeff(k,k);
m_eivec.coeffRef(k,k) = Scalar(1.0,0.0);
for(int i=k-1 ; i>=0 ; i--)
{
m_eivec.coeffRef(i,k) = -schur.matrixT().coeff(i,k);
if(k-i-1>0)
m_eivec.coeffRef(i,k) -= (schur.matrixT().row(i).segment(i+1,k-i-1) * m_eivec.col(k).segment(i+1,k-i-1)).value();
z = schur.matrixT().coeff(i,i) - d2;
if(z==Scalar(0))
ei_real_ref(z) = eps * norm;
m_eivec.coeffRef(i,k) = m_eivec.coeff(i,k) / z;
}
m_eivec.col(k).normalize();
}
m_eivec = schur.matrixU() * m_eivec;
m_isInitialized = true;
// sort the eigenvalues
{
for (int i=0; i<n; i++)
{
int k;
m_eivalues.cwise().abs().end(n-i).minCoeff(&k);
if (k != 0)
{
k += i;
std::swap(m_eivalues[k],m_eivalues[i]);
m_eivec.col(i).swap(m_eivec.col(k));
}
}
}
}
#endif // EIGEN_COMPLEX_EIGEN_SOLVER_H

1
Eigen/src/Eigenvalues/EigenSolver.h
View File

@ -194,6 +194,7 @@ EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matr
assert(matrix.cols() == matrix.rows());
int n = matrix.cols();
m_eivalues.resize(n,1);
m_eivec.resize(n,n);
MatrixType matH = matrix;
RealVectorType ort(n);

1
Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
View File

@ -168,6 +168,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::compute(
assert(matrix.cols() == matrix.rows());
int n = matrix.cols();
m_eivalues.resize(n,1);
m_eivec.resize(n,n);
if(n==1)
{

4
Eigen/src/Householder/HouseholderSequence.h
View File

@ -153,10 +153,12 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence
}
protected:
typename VectorsType::Nested m_vectors;
typename CoeffsType::Nested m_coeffs;
bool m_trans;
private:
HouseholderSequence& operator=(const HouseholderSequence&);
};
template<typename VectorsType, typename CoeffsType>

74
Eigen/src/LU/Inverse.h
View File

@ -95,8 +95,8 @@ bool ei_compute_inverse_size3(const XprType& matrix, MatrixType* result)
return true;
}
template<typename MatrixType>
bool ei_compute_inverse_size4_helper(const MatrixType& matrix, MatrixType* result)
template<typename MatrixType, typename ResultType>
bool ei_compute_inverse_size4_helper(const MatrixType& matrix, ResultType* result)
{
/* Let's split M into four 2x2 blocks:
* (P Q)
@ -195,47 +195,47 @@ bool ei_compute_inverse_size4_with_check(const XprType& matrix, MatrixType* resu
*** Part 2 : selector and MatrixBase methods ***
***********************************************/
template<typename MatrixType, int Size = MatrixType::RowsAtCompileTime>
template<typename MatrixType, typename ResultType, int Size = MatrixType::RowsAtCompileTime>
struct ei_compute_inverse
{
static inline void run(const MatrixType& matrix, MatrixType* result)
static inline void run(const MatrixType& matrix, ResultType* result)
{
result = matrix.partialLu().inverse();
}
};
template<typename MatrixType>
struct ei_compute_inverse<MatrixType, 1>
template<typename MatrixType, typename ResultType>
struct ei_compute_inverse<MatrixType, ResultType, 1>
{
static inline void run(const MatrixType& matrix, MatrixType* result)
static inline void run(const MatrixType& matrix, ResultType* result)
{
typedef typename MatrixType::Scalar Scalar;
result->coeffRef(0,0) = Scalar(1) / matrix.coeff(0,0);
}
};
template<typename MatrixType>
struct ei_compute_inverse<MatrixType, 2>
template<typename MatrixType, typename ResultType>
struct ei_compute_inverse<MatrixType, ResultType, 2>
{
static inline void run(const MatrixType& matrix, MatrixType* result)
static inline void run(const MatrixType& matrix, ResultType* result)
{
ei_compute_inverse_size2(matrix, result);
}
};
template<typename MatrixType>
struct ei_compute_inverse<MatrixType, 3>
template<typename MatrixType, typename ResultType>
struct ei_compute_inverse<MatrixType, ResultType, 3>
{
static inline void run(const MatrixType& matrix, MatrixType* result)
static inline void run(const MatrixType& matrix, ResultType* result)
{
ei_compute_inverse_size3<false, MatrixType, MatrixType>(matrix, result);
ei_compute_inverse_size3<false, MatrixType, ResultType>(matrix, result);
}
};
template<typename MatrixType>
struct ei_compute_inverse<MatrixType, 4>
template<typename MatrixType, typename ResultType>
struct ei_compute_inverse<MatrixType, ResultType, 4>
{
static inline void run(const MatrixType& matrix, MatrixType* result)
static inline void run(const MatrixType& matrix, ResultType* result)
{
ei_compute_inverse_size4_with_check(matrix, result);
}
@ -256,11 +256,12 @@ struct ei_compute_inverse<MatrixType, 4>
* \sa inverse(), computeInverseWithCheck()
*/
template<typename Derived>
inline void MatrixBase<Derived>::computeInverse(PlainMatrixType *result) const
template<typename ResultType>
inline void MatrixBase<Derived>::computeInverse(ResultType *result) const
{
ei_assert(rows() == cols());
EIGEN_STATIC_ASSERT(NumTraits<Scalar>::HasFloatingPoint,NUMERIC_TYPE_MUST_BE_FLOATING_POINT)
ei_compute_inverse<PlainMatrixType>::run(eval(), result);
ei_compute_inverse<PlainMatrixType, ResultType>::run(eval(), result);
}
/** \lu_module
@ -289,10 +290,10 @@ inline const typename MatrixBase<Derived>::PlainMatrixType MatrixBase<Derived>::
* Compute inverse with invertibility check *
*******************************************/
template<typename MatrixType, int Size = MatrixType::RowsAtCompileTime>
template<typename MatrixType, typename ResultType, int Size = MatrixType::RowsAtCompileTime>
struct ei_compute_inverse_with_check
{
static inline bool run(const MatrixType& matrix, MatrixType* result)
static inline bool run(const MatrixType& matrix, ResultType* result)
{
typedef typename MatrixType::Scalar Scalar;
LU<MatrixType> lu( matrix );
@ -302,10 +303,10 @@ struct ei_compute_inverse_with_check
}
};
template<typename MatrixType>
struct ei_compute_inverse_with_check<MatrixType, 1>
template<typename MatrixType, typename ResultType>
struct ei_compute_inverse_with_check<MatrixType, ResultType, 1>
{
static inline bool run(const MatrixType& matrix, MatrixType* result)
static inline bool run(const MatrixType& matrix, ResultType* result)
{
typedef typename MatrixType::Scalar Scalar;
if( matrix.coeff(0,0) == Scalar(0) ) return false;
@ -314,28 +315,28 @@ struct ei_compute_inverse_with_check<MatrixType, 1>
}
};
template<typename MatrixType>
struct ei_compute_inverse_with_check<MatrixType, 2>
template<typename MatrixType, typename ResultType>
struct ei_compute_inverse_with_check<MatrixType, ResultType, 2>
{
static inline bool run(const MatrixType& matrix, MatrixType* result)
static inline bool run(const MatrixType& matrix, ResultType* result)
{
return ei_compute_inverse_size2_with_check(matrix, result);
}
};
template<typename MatrixType>
struct ei_compute_inverse_with_check<MatrixType, 3>
template<typename MatrixType, typename ResultType>
struct ei_compute_inverse_with_check<MatrixType, ResultType, 3>
{
static inline bool run(const MatrixType& matrix, MatrixType* result)
static inline bool run(const MatrixType& matrix, ResultType* result)
{
return ei_compute_inverse_size3<true, MatrixType, MatrixType>(matrix, result);
return ei_compute_inverse_size3<true, MatrixType, ResultType>(matrix, result);
}
};
template<typename MatrixType>
struct ei_compute_inverse_with_check<MatrixType, 4>
template<typename MatrixType, typename ResultType>
struct ei_compute_inverse_with_check<MatrixType, ResultType, 4>
{
static inline bool run(const MatrixType& matrix, MatrixType* result)
static inline bool run(const MatrixType& matrix, ResultType* result)
{
return ei_compute_inverse_size4_with_check(matrix, result);
}
@ -352,11 +353,12 @@ struct ei_compute_inverse_with_check<MatrixType, 4>
* \sa inverse(), computeInverse()
*/
template<typename Derived>
inline bool MatrixBase<Derived>::computeInverseWithCheck(PlainMatrixType *result) const
template<typename ResultType>
inline bool MatrixBase<Derived>::computeInverseWithCheck(ResultType *result) const
{
ei_assert(rows() == cols());
EIGEN_STATIC_ASSERT(NumTraits<Scalar>::HasFloatingPoint,NUMERIC_TYPE_MUST_BE_FLOATING_POINT)
return ei_compute_inverse_with_check<PlainMatrixType>::run(eval(), result);
return ei_compute_inverse_with_check<PlainMatrixType, ResultType>::run(eval(), result);
}

3
Eigen/src/QR/ColPivotingHouseholderQR.h
View File

@ -62,7 +62,7 @@ template<typename MatrixType> class ColPivotingHouseholderQR
typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType;
typedef Matrix<RealScalar, 1, ColsAtCompileTime> RealRowVectorType;
typedef typename HouseholderSequence<MatrixQType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
/**
* \brief Default Constructor.
@ -351,7 +351,6 @@ bool ColPivotingHouseholderQR<MatrixType>::solve(
}
const int rows = m_qr.rows();
const int cols = b.cols();
ei_assert(b.rows() == rows);
typename OtherDerived::PlainMatrixType c(b);

20
Eigen/src/QR/HouseholderQR.h
View File

@ -59,10 +59,10 @@ template<typename MatrixType> class HouseholderQR
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType;
typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, AutoAlign | (ei_traits<MatrixType>::Flags&RowMajorBit ? RowMajor : ColMajor)> MatrixQType;
typedef Matrix<Scalar, DiagSizeAtCompileTime, 1> HCoeffsType;
typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
typedef typename HouseholderSequence<MatrixQType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
/**
* \brief Default Constructor.
@ -206,18 +206,22 @@ void HouseholderQR<MatrixType>::solve(
) const
{
ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
result->derived().resize(m_qr.cols(), b.cols());
const int rows = m_qr.rows();
const int cols = b.cols();
const int rank = std::min(m_qr.rows(), m_qr.cols());
ei_assert(b.rows() == rows);
result->resize(rows, cols);
*result = b;
result->applyOnTheLeft(matrixQAsHouseholderSequence().inverse());
typename OtherDerived::PlainMatrixType c(b);
// Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
c.applyOnTheLeft(makeHouseholderSequence(m_qr.corner(TopLeft,rows,rank), m_hCoeffs.start(rank)).transpose());
const int rank = std::min(result->rows(), result->cols());
m_qr.corner(TopLeft, rank, rank)
.template triangularView<UpperTriangular>()
.solveInPlace(result->corner(TopLeft, rank, result->cols()));
.solveInPlace(c.corner(TopLeft, rank, c.cols()));
result->corner(TopLeft, rank, c.cols()) = c.corner(TopLeft,rank, c.cols());
result->corner(BottomLeft, result->rows()-rank, c.cols()).setZero();
}
/** \returns the matrix Q */

3
Eigen/src/SVD/JacobiSVD.h
View File

@ -26,8 +26,9 @@
#define EIGEN_JACOBISVD_H
// forward declarations (needed by ICC)
// the empty bodies are required by VC
template<typename MatrixType, unsigned int Options, bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex>
struct ei_svd_precondition_2x2_block_to_be_real;
struct ei_svd_precondition_2x2_block_to_be_real {};
template<typename MatrixType, unsigned int Options,
bool PossiblyMoreRowsThanCols = (Options & AtLeastAsManyColsAsRows) == 0

2
Eigen/src/Sparse/SparseMatrixBase.h
View File

@ -515,7 +515,7 @@ template<typename Derived> class SparseMatrixBase : public AnyMatrixBase<Derived
{ return typename ei_eval<Derived>::type(derived()); }
// template<typename OtherDerived>
// void swap(const MatrixBase<OtherDerived>& other);
// void swap(MatrixBase<OtherDerived> EIGEN_REF_TO_TEMPORARY other);
template<unsigned int Added>
const SparseFlagged<Derived, Added, 0> marked() const;

35
INSTALL Normal file
View File

@ -0,0 +1,35 @@
Installation instructions for Eigen
***********************************
Explanation before starting
***************************
Eigen consists only of header files, hence there is nothing to compile
before you can use it. Moreover, these header files do not depend on your
platform, they are the same for everybody.
Method 1. Installing without using CMake
****************************************
You can use right away the headers in the Eigen/ subdirectory. In order
to install, just copy this Eigen/ subdirectory to your favorite location.
If you also want the unsupported features, copy the unsupported/
subdirectory too.
Method 2. Installing using CMake
********************************
Let's call this directory 'source_dir' (where this INSTALL file is).
Before starting, create another directory which we will call 'build_dir'.
Do:
cd build_dir
cmake source_dir
make install
The "make install" step may require administrator privileges.
You can adjust the installation destination (the "prefix")
by passing the -DCMAKE_INSTALL_PREFIX=myprefix option to cmake, as is
explained in the message that cmake prints at the end.

118
Mainpage.dox
View File

@ -1,127 +1,19 @@
// Please don't remove the following lines:
// this is the only way to specify doxygen options
// to api.kde.org's scripts
// DOXYGEN_SET_PROJECT_NAME = Eigen2
// DOXYGEN_SET_PROJECT_NUMBER = "2.0 - trunk"
// DOXYGEN_SET_CREATE_SUBDIRS = NO
// DOXYGEN_SET_BRIEF_MEMBER_DESC = YES
// DOXYGEN_SET_REPEAT_BRIEF = YES
// DOXYGEN_SET_ALWAYS_DETAILED_SEC = NO
// DOXYGEN_SET_INLINE_INHERITED_MEMB = NO
// DOXYGEN_SET_FULL_PATH_NAMES = NO
// DOXYGEN_SET_SHORT_NAMES = NO
// DOXYGEN_SET_JAVADOC_AUTOBRIEF = NO
// DOXYGEN_SET_QT_AUTOBRIEF = NO
// DOXYGEN_SET_MULTILINE_CPP_IS_BRIEF = NO
// DOXYGEN_SET_DETAILS_AT_TOP = YES
// DOXYGEN_SET_INHERIT_DOCS = YES
// DOXYGEN_SET_ALIASES = "only_for_vectors=This is only for vectors (either row-vectors or column-vectors), i.e. matrices which are known at compile-time to have either one row or one column." "array_module=This is defined in the %Array module. \code #include <Eigen/Array> \endcode" "lu_module=This is defined in the %LU module. \code #include <Eigen/LU> \endcode" "cholesky_module=This is defined in the %Cholesky module. \code #include <Eigen/Cholesky> \endcode" "qr_module=This is defined in the %QR module. \code #include <Eigen/QR> \endcode" "svd_module=This is defined in the %SVD module. \code #include <Eigen/SVD> \endcode" "geometry_module=This is defined in the %Geometry module. \code #include <Eigen/Geometry> \endcode" "leastsquares_module=This is defined in the %LeastSquares module. \code #include <Eigen/LeastSquares> \endcode" "addexample=\anchor" "label=\bug" "redstar=<a href='#warningarraymodule' style='color:red;text-decoration: none;'>*</a>" "nonstableyet=\warning This is not considered to be part of the stable public API yet. Changes may happen in future releases. See \ref Experimental \"Experimental parts of Eigen\""
// DOXYGEN_SET_DISTRIBUTE_GROUP_DOC = NO
// DOXYGEN_SET_SUBGROUPING = YES
// DOXYGEN_SET_TYPEDEF_HIDES_STRUCT = NO
// DOXYGEN_SET_EXTRACT_ALL = NO
// DOXYGEN_SET_EXTRACT_PRIVATE = NO
// DOXYGEN_SET_EXTRACT_STATIC = NO
// DOXYGEN_SET_EXTRACT_LOCAL_CLASSES = NO
// DOXYGEN_SET_EXTRACT_LOCAL_METHODS = NO
// DOXYGEN_SET_EXTRACT_ANON_NSPACES = NO
// DOXYGEN_SET_HIDE_UNDOC_MEMBERS = NO
// DOXYGEN_SET_HIDE_UNDOC_CLASSES = YES
// DOXYGEN_SET_HIDE_FRIEND_COMPOUNDS = YES
// DOXYGEN_SET_HIDE_IN_BODY_DOCS = NO
// DOXYGEN_SET_INTERNAL_DOCS = NO
// DOXYGEN_SET_CASE_SENSE_NAMES = YES
// DOXYGEN_SET_HIDE_SCOPE_NAMES = YES
// DOXYGEN_SET_SHOW_INCLUDE_FILES = YES
// DOXYGEN_SET_INLINE_INFO = YES
// DOXYGEN_SET_SORT_MEMBER_DOCS = YES
// DOXYGEN_SET_SORT_BRIEF_DOCS = YES
// DOXYGEN_SET_SORT_GROUP_NAMES = NO
// DOXYGEN_SET_SORT_BY_SCOPE_NAME = NO
// DOXYGEN_SET_GENERATE_TODOLIST = NO
// DOXYGEN_SET_GENERATE_TESTLIST = NO
// DOXYGEN_SET_GENERATE_BUGLIST = NO
// DOXYGEN_SET_GENERATE_DEPRECATEDLIST= NO
// DOXYGEN_SET_SHOW_USED_FILES = YES
// DOXYGEN_SET_SHOW_DIRECTORIES = NO
// DOXYGEN_SET_SHOW_FILES = YES
// DOXYGEN_SET_SHOW_NAMESPACES = NO
// DOXYGEN_SET_WARN_IF_UNDOCUMENTED = NO
// DOXYGEN_SET_WARN_NO_PARAMDOC = NO
// DOXYGEN_SET_INPUT = @topdir@/eigen2/Eigen @topdir@/eigen2/doc @topdir@/eigen2/build/doc @topdir@/eigen2/unsupported/Eigen
// DOXYGEN_SET_EXCLUDE = *.sh *.in
// DOXYGEN_SET_EXAMPLE_PATH = @topdir@/eigen2/doc/snippets/ @topdir@/eigen2/doc/examples/ @topdir@/eigen2/build/doc/examples/ @topdir@/eigen2/build/doc/snippets/
// DOXYGEN_SET_FILE_PATTERNS = *
// DOXYGEN_SET_RECURSIVE = NO
// DOXYGEN_SET_FILTER_SOURCE_FILES = YES
// DOXYGEN_EXCLUDE_SYMBOLS = MatrixBase<* MapBase<* RotationBase<* Matrix<*
// DOXYGEN_SET_SOURCE_BROWSER = NO
// DOXYGEN_SET_INLINE_SOURCES = NO
// DOXYGEN_SET_STRIP_CODE_COMMENTS = YES
// DOXYGEN_SET_REFERENCED_BY_RELATION = YES
// DOXYGEN_SET_REFERENCES_RELATION = YES
// DOXYGEN_SET_REFERENCES_LINK_SOURCE = YES
// DOXYGEN_SET_VERBATIM_HEADERS = YES
// DOXYGEN_SET_ALPHABETICAL_INDEX = NO
// DOXYGEN_SET_HTML_ALIGN_MEMBERS = YES
// DOXYGEN_SET_GENERATE_TREEVIEW = NO
// DOXYGEN_SET_FORMULA_FONTSIZE = 12
// DOXYGEN_SET_GENERATE_LATEX = NO
// DOXYGEN_SET_EXTRA_PACKAGES = amssymb
// DOXYGEN_SET_ENABLE_PREPROCESSING = YES
// DOXYGEN_SET_MACRO_EXPANSION = YES
// DOXYGEN_SET_EXPAND_ONLY_PREDEF = YES
// DOXYGEN_SET_SEARCH_INCLUDES = YES
// DOXYGEN_SET_PREDEFINED = EIGEN_EMPTY_STRUCT EIGEN_PARSED_BY_DOXYGEN EIGEN_VECTORIZE EIGEN_QT_SUPPORT EIGEN_STRONG_INLINE=inline
// DOXYGEN_SET_EXPAND_AS_DEFINED = EIGEN_MAKE_SCALAR_OPS EIGEN_MAKE_TYPEDEFS EIGEN_MAKE_TYPEDEFS_ALL_SIZES EIGEN_CWISE_UNOP_RETURN_TYPE EIGEN_CWISE_BINOP_RETURN_TYPE
// DOXYGEN_SET_SKIP_FUNCTION_MACROS = YES
// DOXYGEN_SET_CLASS_DIAGRAMS = NO
// DOXYGEN_SET_HIDE_UNDOC_RELATIONS = NO
// DOXYGEN_SET_HAVE_DOT = NO
// DOXYGEN_SET_CLASS_GRAPH = NO
// DOXYGEN_SET_COLLABORATION_GRAPH = NO
// DOXYGEN_SET_GROUP_GRAPHS = NO
// DOXYGEN_SET_UML_LOOK = NO
// DOXYGEN_SET_TEMPLATE_RELATIONS = NO
// DOXYGEN_SET_INCLUDE_GRAPH = NO
// DOXYGEN_SET_INCLUDED_BY_GRAPH = NO
// DOXYGEN_SET_CALL_GRAPH = NO
// DOXYGEN_SET_CALLER_GRAPH = NO
// DOXYGEN_SET_GRAPHICAL_HIERARCHY = NO
// DOXYGEN_SET_DIRECTORY_GRAPH = NO
o /** \mainpage Eigen
<h3>If you see this page, then you have not properly generated the documentation. Namely, you have run doxygen from the source directory, which is not appropriate for generating the documentation of Eigen.</h3>
In order to generate the documentation of Eigen, please follow these steps:
<ul>
<li>make sure you have the required software installed: cmake, doxygen, and a C++ compiler.
<li>make sure you have the required software installed: CMake, Doxygen, LaTeX, and a C++ compiler.
<li>create a new directory, which we will call the "build directory", outside of the Eigen source directory.</li>
<li>enter the build directory</li>
<li>configure the project: <pre>cmake /path/to/source/directory</pre></li>
<li>now generate the documentaion: <pre>make doc</pre> or, if you have two CPUs, <pre>make doc -j2</pre> Note that this will compile the examples, run them, and integrate their output into the documentation, which can take some time.</li>
<li>now generate the documentaion: <pre>make doc</pre> or, if you have two CPUs, <pre>make doc -j3</pre> Note that this will compile the examples, run them, and integrate their output into the documentation, which can take some time.</li>
</ul>
After doing that, you will find the HTML documentation in the doc/html/ subdirectory of the build directory.
<h2>Note however that the documentation is available online here:
<a href="http://eigen.tuxfamily.org/dox/">http://eigen.tuxfamily.org/dox</a></h2>
<h2>Note however that the documentation is available online here:</h2>
<a href="http://eigen.tuxfamily.org/dox/">http://eigen.tuxfamily.org/dox</a> for the stable version and
<a href="http://eigen.tuxfamily.org/dox-devel/">http://eigen.tuxfamily.org/dox-devel</a> for the development branch version.
*/

12
bench/BenchSparseUtil.h
View File

@ -27,23 +27,23 @@ typedef SparseMatrix<Scalar> EigenSparseMatrix;
void fillMatrix(float density, int rows, int cols, EigenSparseMatrix& dst)
{
dst.startFill(rows*cols*density);
dst.reserve(rows*cols*density);
for(int j = 0; j < cols; j++)
{
for(int i = 0; i < rows; i++)
{
Scalar v = (ei_random<float>(0,1) < density) ? ei_random<Scalar>() : 0;
if (v!=0)
dst.fillrand(i,j) = v;
dst.insert(i,j) = v;
}
}
dst.endFill();
dst.finalize();
}
void fillMatrix2(int nnzPerCol, int rows, int cols, EigenSparseMatrix& dst)
{
std::cout << "alloc " << nnzPerCol*cols << "\n";
dst.startFill(nnzPerCol*cols);
dst.reserve(nnzPerCol*cols);
for(int j = 0; j < cols; j++)
{
std::set<int> aux;
@ -54,10 +54,10 @@ void fillMatrix2(int nnzPerCol, int rows, int cols, EigenSparseMatrix& dst)
k = ei_random<int>(0,rows-1);
aux.insert(k);
dst.fillrand(k,j) = ei_random<Scalar>();
dst.insert(k,j) = ei_random<Scalar>();
}
}
dst.endFill();
dst.finalize();
}
void eiToDense(const EigenSparseMatrix& src, DenseMatrix& dst)

320
bench/sparse_setter.cpp
View File

@ -12,7 +12,15 @@
#endif
#ifndef REPEAT
#define REPEAT 1
#define REPEAT 2
#endif
#ifndef NBTRIES
#define NBTRIES 2
#endif
#ifndef KK
#define KK 10
#endif
#ifndef NOGOOGLE
@ -22,7 +30,7 @@
#include "BenchSparseUtil.h"
#define CHECK_MEM
#define CHECK_MEM
// #define CHECK_MEM std/**/::cout << "check mem\n"; getchar();
#define BENCH(X) \
@ -37,9 +45,13 @@ typedef std::vector<Vector2i> Coordinates;
typedef std::vector<float> Values;
EIGEN_DONT_INLINE Scalar* setinnerrand_eigen(const Coordinates& coords, const Values& vals);
EIGEN_DONT_INLINE Scalar* setrand_eigen_dynamic(const Coordinates& coords, const Values& vals);
EIGEN_DONT_INLINE Scalar* setrand_eigen_compact(const Coordinates& coords, const Values& vals);
EIGEN_DONT_INLINE Scalar* setrand_eigen_sumeq(const Coordinates& coords, const Values& vals);
EIGEN_DONT_INLINE Scalar* setrand_eigen_gnu_hash(const Coordinates& coords, const Values& vals);
EIGEN_DONT_INLINE Scalar* setrand_eigen_google_dense(const Coordinates& coords, const Values& vals);
EIGEN_DONT_INLINE Scalar* setrand_eigen_google_sparse(const Coordinates& coords, const Values& vals);
EIGEN_DONT_INLINE Scalar* setrand_scipy(const Coordinates& coords, const Values& vals);
EIGEN_DONT_INLINE Scalar* setrand_ublas_mapped(const Coordinates& coords, const Values& vals);
EIGEN_DONT_INLINE Scalar* setrand_ublas_coord(const Coordinates& coords, const Values& vals);
EIGEN_DONT_INLINE Scalar* setrand_ublas_compressed(const Coordinates& coords, const Values& vals);
@ -50,17 +62,36 @@ int main(int argc, char *argv[])
{
int rows = SIZE;
int cols = SIZE;
bool fullyrand = false;
//float density = float(NBPERROW)/float(SIZE);
bool fullyrand = true;
BenchTimer timer;
Coordinates coords;
Values values;
if(fullyrand)
{
for (int i=0; i<cols*NBPERROW; ++i)
Coordinates pool;
pool.reserve(cols*NBPERROW);
std::cerr << "fill pool" << "\n";
for (int i=0; i<cols*NBPERROW; )
{
coords.push_back(Vector2i(ei_random<int>(0,rows-1),ei_random<int>(0,cols-1)));
// DynamicSparseMatrix<int> stencil(SIZE,SIZE);
Vector2i ij(ei_random<int>(0,rows-1),ei_random<int>(0,cols-1));
// if(stencil.coeffRef(ij.x(), ij.y())==0)
{
// stencil.coeffRef(ij.x(), ij.y()) = 1;
pool.push_back(ij);
}
++i;
}
std::cerr << "pool ok" << "\n";
int n = cols*NBPERROW*KK;
coords.reserve(n);
values.reserve(n);
for (int i=0; i<n; ++i)
{
int i = ei_random<int>(0,pool.size());
coords.push_back(pool[i]);
values.push_back(ei_random<Scalar>());
}
}
@ -79,67 +110,55 @@ int main(int argc, char *argv[])
// dense matrices
#ifdef DENSEMATRIX
{
timer.reset();
timer.start();
for (int k=0; k<REPEAT; ++k)
setrand_eigen_dense(coords,values);
timer.stop();
BENCH(setrand_eigen_dense(coords,values);)
std::cout << "Eigen Dense\t" << timer.value() << "\n";
}
#endif
// eigen sparse matrices
if (!fullyrand)
// if (!fullyrand)
// {
// BENCH(setinnerrand_eigen(coords,values);)
// std::cout << "Eigen fillrand\t" << timer.value() << "\n";
// }
{
timer.reset();
timer.start();
for (int k=0; k<REPEAT; ++k)
setinnerrand_eigen(coords,values);
timer.stop();
std::cout << "Eigen fillrand\t" << timer.value() << "\n";
BENCH(setrand_eigen_dynamic(coords,values);)
std::cout << "Eigen dynamic\t" << timer.value() << "\n";
}
// {
// BENCH(setrand_eigen_compact(coords,values);)
// std::cout << "Eigen compact\t" << timer.value() << "\n";
// }
{
BENCH(setrand_eigen_sumeq(coords,values);)
std::cout << "Eigen sumeq\t" << timer.value() << "\n";
}
{
timer.reset();
timer.start();
for (int k=0; k<REPEAT; ++k)
setrand_eigen_gnu_hash(coords,values);
timer.stop();
std::cout << "Eigen std::map\t" << timer.value() << "\n";
// BENCH(setrand_eigen_gnu_hash(coords,values);)
// std::cout << "Eigen std::map\t" << timer.value() << "\n";
}
{
BENCH(setrand_scipy(coords,values);)
std::cout << "scipy\t" << timer.value() << "\n";
}
#ifndef NOGOOGLE
{
timer.reset();
timer.start();
for (int k=0; k<REPEAT; ++k)
setrand_eigen_google_dense(coords,values);
timer.stop();
BENCH(setrand_eigen_google_dense(coords,values);)
std::cout << "Eigen google dense\t" << timer.value() << "\n";
}
{
timer.reset();
timer.start();
for (int k=0; k<REPEAT; ++k)
setrand_eigen_google_sparse(coords,values);
timer.stop();
BENCH(setrand_eigen_google_sparse(coords,values);)
std::cout << "Eigen google sparse\t" << timer.value() << "\n";
}
#endif
#ifndef NOUBLAS
{
timer.reset();
timer.start();
for (int k=0; k<REPEAT; ++k)
setrand_ublas_mapped(coords,values);
timer.stop();
std::cout << "ublas mapped\t" << timer.value() << "\n";
// BENCH(setrand_ublas_mapped(coords,values);)
// std::cout << "ublas mapped\t" << timer.value() << "\n";
}
{
timer.reset();
timer.start();
for (int k=0; k<REPEAT; ++k)
setrand_ublas_genvec(coords,values);
timer.stop();
BENCH(setrand_ublas_genvec(coords,values);)
std::cout << "ublas vecofvec\t" << timer.value() << "\n";
}
/*{
@ -159,16 +178,12 @@ int main(int argc, char *argv[])
std::cout << "ublas coord\t" << timer.value() << "\n";
}*/
#endif
// MTL4
#ifndef NOMTL
{
timer.reset();
timer.start();
for (int k=0; k<REPEAT; ++k)
setrand_mtl(coords,values);
timer.stop();
BENCH(setrand_mtl(coords,values));
std::cout << "MTL\t" << timer.value() << "\n";
}
#endif
@ -180,16 +195,63 @@ EIGEN_DONT_INLINE Scalar* setinnerrand_eigen(const Coordinates& coords, const Va
{
using namespace Eigen;
SparseMatrix<Scalar> mat(SIZE,SIZE);
mat.startFill(2000000/*coords.size()*/);
//mat.startFill(2000000/*coords.size()*/);
for (int i=0; i<coords.size(); ++i)
{
mat.fillrand(coords[i].x(), coords[i].y()) = vals[i];
mat.insert(coords[i].x(), coords[i].y()) = vals[i];
}
mat.endFill();
mat.finalize();
CHECK_MEM;
return 0;
}
EIGEN_DONT_INLINE Scalar* setrand_eigen_dynamic(const Coordinates& coords, const Values& vals)
{
using namespace Eigen;
DynamicSparseMatrix<Scalar> mat(SIZE,SIZE);
mat.reserve(coords.size()/10);
for (int i=0; i<coords.size(); ++i)
{
mat.coeffRef(coords[i].x(), coords[i].y()) += vals[i];
}
mat.finalize();
CHECK_MEM;
return &mat.coeffRef(coords[0].x(), coords[0].y());
}
EIGEN_DONT_INLINE Scalar* setrand_eigen_sumeq(const Coordinates& coords, const Values& vals)
{
using namespace Eigen;
int n = coords.size()/KK;
DynamicSparseMatrix<Scalar> mat(SIZE,SIZE);
for (int j=0; j<KK; ++j)
{
DynamicSparseMatrix<Scalar> aux(SIZE,SIZE);
mat.reserve(n);
for (int i=j*n; i<(j+1)*n; ++i)
{
aux.insert(coords[i].x(), coords[i].y()) += vals[i];
}
aux.finalize();
mat += aux;
}
return &mat.coeffRef(coords[0].x(), coords[0].y());
}
EIGEN_DONT_INLINE Scalar* setrand_eigen_compact(const Coordinates& coords, const Values& vals)
{
using namespace Eigen;
DynamicSparseMatrix<Scalar> setter(SIZE,SIZE);
setter.reserve(coords.size()/10);
for (int i=0; i<coords.size(); ++i)
{
setter.coeffRef(coords[i].x(), coords[i].y()) += vals[i];
}
SparseMatrix<Scalar> mat = setter;
CHECK_MEM;
return &mat.coeffRef(coords[0].x(), coords[0].y());
}
EIGEN_DONT_INLINE Scalar* setrand_eigen_gnu_hash(const Coordinates& coords, const Values& vals)
{
using namespace Eigen;
@ -198,11 +260,11 @@ EIGEN_DONT_INLINE Scalar* setrand_eigen_gnu_hash(const Coordinates& coords, cons
RandomSetter<SparseMatrix<Scalar>, StdMapTraits > setter(mat);
for (int i=0; i<coords.size(); ++i)
{
setter(coords[i].x(), coords[i].y()) = vals[i];
setter(coords[i].x(), coords[i].y()) += vals[i];
}
CHECK_MEM;
}
return 0;//&mat.coeffRef(coords[0].x(), coords[0].y());
return &mat.coeffRef(coords[0].x(), coords[0].y());
}
#ifndef NOGOOGLE
@ -213,10 +275,10 @@ EIGEN_DONT_INLINE Scalar* setrand_eigen_google_dense(const Coordinates& coords,
{
RandomSetter<SparseMatrix<Scalar>, GoogleDenseHashMapTraits> setter(mat);
for (int i=0; i<coords.size(); ++i)
setter(coords[i].x(), coords[i].y()) = vals[i];
setter(coords[i].x(), coords[i].y()) += vals[i];
CHECK_MEM;
}
return 0;//&mat.coeffRef(coords[0].x(), coords[0].y());
return &mat.coeffRef(coords[0].x(), coords[0].y());
}
EIGEN_DONT_INLINE Scalar* setrand_eigen_google_sparse(const Coordinates& coords, const Values& vals)
@ -226,13 +288,139 @@ EIGEN_DONT_INLINE Scalar* setrand_eigen_google_sparse(const Coordinates& coords,
{
RandomSetter<SparseMatrix<Scalar>, GoogleSparseHashMapTraits> setter(mat);
for (int i=0; i<coords.size(); ++i)
setter(coords[i].x(), coords[i].y()) = vals[i];
setter(coords[i].x(), coords[i].y()) += vals[i];
CHECK_MEM;
}
return 0;//&mat.coeffRef(coords[0].x(), coords[0].y());
return &mat.coeffRef(coords[0].x(), coords[0].y());
}
#endif
template <class T>
void coo_tocsr(const int n_row,
const int n_col,
const int nnz,
const Coordinates Aij,
const Values Ax,
int Bp[],
int Bj[],
T Bx[])
{
//compute number of non-zero entries per row of A coo_tocsr
std::fill(Bp, Bp + n_row, 0);
for (int n = 0; n < nnz; n++){
Bp[Aij[n].x()]++;
}
//cumsum the nnz per row to get Bp[]
for(int i = 0, cumsum = 0; i < n_row; i++){
int temp = Bp[i];
Bp[i] = cumsum;
cumsum += temp;
}
Bp[n_row] = nnz;
//write Aj,Ax into Bj,Bx
for(int n = 0; n < nnz; n++){
int row = Aij[n].x();
int dest = Bp[row];
Bj[dest] = Aij[n].y();
Bx[dest] = Ax[n];
Bp[row]++;
}
for(int i = 0, last = 0; i <= n_row; i++){
int temp = Bp[i];
Bp[i] = last;
last = temp;
}
//now Bp,Bj,Bx form a CSR representation (with possible duplicates)
}
template< class T1, class T2 >
bool kv_pair_less(const std::pair<T1,T2>& x, const std::pair<T1,T2>& y){
return x.first < y.first;
}
template<class I, class T>
void csr_sort_indices(const I n_row,
const I Ap[],
I Aj[],
T Ax[])
{
std::vector< std::pair<I,T> > temp;
for(I i = 0; i < n_row; i++){
I row_start = Ap[i];
I row_end = Ap[i+1];
temp.clear();
for(I jj = row_start; jj < row_end; jj++){
temp.push_back(std::make_pair(Aj[jj],Ax[jj]));
}
std::sort(temp.begin(),temp.end(),kv_pair_less<I,T>);
for(I jj = row_start, n = 0; jj < row_end; jj++, n++){
Aj[jj] = temp[n].first;
Ax[jj] = temp[n].second;
}
}
}
template <class I, class T>
void csr_sum_duplicates(const I n_row,
const I n_col,
I Ap[],
I Aj[],
T Ax[])
{
I nnz = 0;
I row_end = 0;
for(I i = 0; i < n_row; i++){
I jj = row_end;
row_end = Ap[i+1];
while( jj < row_end ){
I j = Aj[jj];
T x = Ax[jj];
jj++;
while( jj < row_end && Aj[jj] == j ){
x += Ax[jj];
jj++;
}
Aj[nnz] = j;
Ax[nnz] = x;
nnz++;
}
Ap[i+1] = nnz;
}
}
EIGEN_DONT_INLINE Scalar* setrand_scipy(const Coordinates& coords, const Values& vals)
{
using namespace Eigen;
SparseMatrix<Scalar> mat(SIZE,SIZE);
mat.resizeNonZeros(coords.size());
// std::cerr << "setrand_scipy...\n";
coo_tocsr<Scalar>(SIZE,SIZE, coords.size(), coords, vals, mat._outerIndexPtr(), mat._innerIndexPtr(), mat._valuePtr());
// std::cerr << "coo_tocsr ok\n";
csr_sort_indices(SIZE, mat._outerIndexPtr(), mat._innerIndexPtr(), mat._valuePtr());
csr_sum_duplicates(SIZE, SIZE, mat._outerIndexPtr(), mat._innerIndexPtr(), mat._valuePtr());
mat.resizeNonZeros(mat._outerIndexPtr()[SIZE]);
return &mat.coeffRef(coords[0].x(), coords[0].y());
}
#ifndef NOUBLAS
EIGEN_DONT_INLINE Scalar* setrand_ublas_mapped(const Coordinates& coords, const Values& vals)
{
@ -242,7 +430,7 @@ EIGEN_DONT_INLINE Scalar* setrand_ublas_mapped(const Coordinates& coords, const
mapped_matrix<Scalar> aux(SIZE,SIZE);
for (int i=0; i<coords.size(); ++i)
{
aux(coords[i].x(), coords[i].y()) = vals[i];
aux(coords[i].x(), coords[i].y()) += vals[i];
}
CHECK_MEM;
compressed_matrix<Scalar> mat(aux);
@ -278,12 +466,12 @@ EIGEN_DONT_INLINE Scalar* setrand_ublas_genvec(const Coordinates& coords, const
using namespace boost;
using namespace boost::numeric;
using namespace boost::numeric::ublas;
// ublas::vector<coordinate_vector<Scalar> > foo;
generalized_vector_of_vector<Scalar, row_major, ublas::vector<coordinate_vector<Scalar> > > aux(SIZE,SIZE);
for (int i=0; i<coords.size(); ++i)
{
aux(coords[i].x(), coords[i].y()) = vals[i];
aux(coords[i].x(), coords[i].y()) += vals[i];
}
CHECK_MEM;
compressed_matrix<Scalar,row_major> mat(aux);

14
blas/CMakeLists.txt Normal file
View File

@ -0,0 +1,14 @@
project(EigenBlas)
add_custom_target(blas)
set(EigenBlas_SRCS single.cpp double.cpp complex_single.cpp complex_double.cpp)
add_library(eigen_blas SHARED ${EigenBlas_SRCS})
add_dependencies(blas eigen_blas)
install(TARGETS eigen_blas
RUNTIME DESTINATION bin
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib)

7
blas/README.txt Normal file
View File

@ -0,0 +1,7 @@
This directory contains a BLAS library built on top of Eigen.
This is currently a work in progress which is far to be ready for use,
but feel free to contribute to it if you wish.
If you want to compile it, set the cmake variable EIGEN_BUILD_BLAS to "on".

115
blas/common.h Normal file
View File

@ -0,0 +1,115 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifndef EIGEN_BLAS_COMMON_H
#define EIGEN_BLAS_COMMON_H
#ifndef SCALAR
#error the token SCALAR must be defined to compile this file
#endif
#ifdef __cplusplus
extern "C"
{
#endif
#include <blas.h>
#ifdef __cplusplus
}
#endif
#define NOTR 0
#define TR 1
#define ADJ 2
#define LEFT 0
#define RIGHT 1
#define UP 0
#define LO 1
#define NUNIT 0
#define UNIT 1
#define OP(X) ( ((X)=='N' || (X)=='n') ? NOTR \
: ((X)=='T' || (X)=='t') ? TR \
: ((X)=='C' || (X)=='c') ? ADJ \
: 0xff)
#define SIDE(X) ( ((X)=='L' || (X)=='l') ? LEFT \
: ((X)=='R' || (X)=='r') ? RIGHT \
: 0xff)
#define UPLO(X) ( ((X)=='U' || (X)=='u') ? UP \
: ((X)=='L' || (X)=='l') ? LO \
: 0xff)
#define DIAG(X) ( ((X)=='N' || (X)=='N') ? NUNIT \
: ((X)=='U' || (X)=='u') ? UNIT \
: 0xff)
#include <Eigen/Core>
#include <Eigen/Jacobi>
using namespace Eigen;
template<typename T>
Block<NestByValue<Map<Matrix<T,Dynamic,Dynamic> > >, Dynamic, Dynamic>
matrix(T* data, int rows, int cols, int stride)
{
return Map<Matrix<T,Dynamic,Dynamic> >(data, stride, cols).nestByValue().block(0,0,rows,cols);
}
template<typename T>
Block<NestByValue<Map<Matrix<T,Dynamic,Dynamic,RowMajor> > >, Dynamic, 1>
vector(T* data, int size, int incr)
{
return Map<Matrix<T,Dynamic,Dynamic,RowMajor> >(data, size, incr).nestByValue().col(0);
}
template<typename T>
Map<Matrix<T,Dynamic,1> >
vector(T* data, int size)
{
return Map<Matrix<T,Dynamic,1> >(data, size);
}
typedef SCALAR Scalar;
typedef NumTraits<Scalar>::Real RealScalar;
typedef std::complex<RealScalar> Complex;
enum
{
IsComplex = Eigen::NumTraits<SCALAR>::IsComplex,
Conj = IsComplex
};
typedef Block<NestByValue<Map<Matrix<Scalar,Dynamic,Dynamic> > >, Dynamic, Dynamic> MatrixType;
typedef Block<NestByValue<Map<Matrix<Scalar,Dynamic,Dynamic, RowMajor> > >, Dynamic, 1> StridedVectorType;
typedef Map<Matrix<Scalar,Dynamic,1> > CompactVectorType;
#define EIGEN_BLAS_FUNC(X) EIGEN_CAT(SCALAR_SUFFIX,X##_)
#endif // EIGEN_BLAS_COMMON_H

21
Eigen/src/QR/QrInstantiations.cpp → blas/complex_double.cpp
View File

@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@ -22,17 +22,10 @@
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifndef EIGEN_EXTERN_INSTANTIATIONS
#define EIGEN_EXTERN_INSTANTIATIONS
#endif
#include "../../Core"
#undef EIGEN_EXTERN_INSTANTIATIONS
#define SCALAR std::complex<double>
#define SCALAR_SUFFIX c
#define ISCOMPLEX 1
#include "../../QR"
namespace Eigen
{
EIGEN_QR_MODULE_INSTANTIATE();
}
#include "level1_impl.h"
#include "level2_impl.h"
#include "level3_impl.h"

18
Eigen/src/Cholesky/CholeskyInstantiations.cpp → blas/complex_single.cpp
View File

@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@ -22,14 +22,10 @@
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifndef EIGEN_EXTERN_INSTANTIATIONS
#define EIGEN_EXTERN_INSTANTIATIONS
#endif
#include "../../Core"
#undef EIGEN_EXTERN_INSTANTIATIONS
#define SCALAR std::complex<float>
#define SCALAR_SUFFIX z
#define ISCOMPLEX 1
#include "../../Cholesky"
namespace Eigen {
EIGEN_CHOLESKY_MODULE_INSTANTIATE();
}
#include "level1_impl.h"
#include "level2_impl.h"
#include "level3_impl.h"

30
Eigen/src/Core/CoreInstantiations.cpp → blas/double.cpp
View File

@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@ -22,26 +22,10 @@
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifdef EIGEN_EXTERN_INSTANTIATIONS
#undef EIGEN_EXTERN_INSTANTIATIONS
#endif
#define SCALAR double
#define SCALAR_SUFFIX d
#define ISCOMPLEX 0
#include "../../Core"
namespace Eigen
{
#define EIGEN_INSTANTIATE_PRODUCT(TYPE) \
template static void ei_cache_friendly_product<TYPE>( \
int _rows, int _cols, int depth, \
bool _lhsRowMajor, const TYPE* _lhs, int _lhsStride, \
bool _rhsRowMajor, const TYPE* _rhs, int _rhsStride, \
bool resRowMajor, TYPE* res, int resStride)
EIGEN_INSTANTIATE_PRODUCT(float);
EIGEN_INSTANTIATE_PRODUCT(double);
EIGEN_INSTANTIATE_PRODUCT(int);
EIGEN_INSTANTIATE_PRODUCT(std::complex<float>);
EIGEN_INSTANTIATE_PRODUCT(std::complex<double>);
}
#include "level1_impl.h"
#include "level2_impl.h"
#include "level3_impl.h"

225
blas/level1_impl.h Normal file
View File

@ -0,0 +1,225 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#include "common.h"
int EIGEN_BLAS_FUNC(axpy)(int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy)
{
Scalar* x = reinterpret_cast<Scalar*>(px);
Scalar* y = reinterpret_cast<Scalar*>(py);
Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
if(*incx==1 && *incy==1)
vector(y,*n) += alpha * vector(x,*n);
else
vector(y,*n,*incy) += alpha * vector(x,*n,*incx);
return 1;
}
// computes the sum of magnitudes of all vector elements or, for a complex vector x, the sum
// res = |Rex1| + |Imx1| + |Rex2| + |Imx2| + ... + |Rexn| + |Imxn|, where x is a vector of order n
RealScalar EIGEN_BLAS_FUNC(asum)(int *n, RealScalar *px, int *incx)
{
int size = IsComplex ? 2* *n : *n;
if(*incx==1)
return vector(px,size).cwise().abs().sum();
else
return vector(px,size,*incx).cwise().abs().sum();
return 1;
}
int EIGEN_BLAS_FUNC(copy)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
{
int size = IsComplex ? 2* *n : *n;
if(*incx==1 && *incy==1)
vector(py,size) = vector(px,size);
else
vector(py,size,*incy) = vector(px,size,*incx);
return 1;
}
// computes a vector-vector dot product.
Scalar EIGEN_BLAS_FUNC(dot)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
{
Scalar* x = reinterpret_cast<Scalar*>(px);
Scalar* y = reinterpret_cast<Scalar*>(py);
if(*incx==1 && *incy==1)
return (vector(x,*n).cwise()*vector(y,*n)).sum();
return (vector(x,*n,*incx).cwise()*vector(y,*n,*incy)).sum();
}
/*
// computes a vector-vector dot product with extended precision.
Scalar EIGEN_BLAS_FUNC(sdot)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
{
// TODO
Scalar* x = reinterpret_cast<Scalar*>(px);
Scalar* y = reinterpret_cast<Scalar*>(py);
if(*incx==1 && *incy==1)
return vector(x,*n).dot(vector(y,*n));
return vector(x,*n,*incx).dot(vector(y,*n,*incy));
}
*/
#if ISCOMPLEX
// computes a dot product of a conjugated vector with another vector.
Scalar EIGEN_BLAS_FUNC(dotc)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
{
Scalar* x = reinterpret_cast<Scalar*>(px);
Scalar* y = reinterpret_cast<Scalar*>(py);
if(*incx==1 && *incy==1)
return vector(x,*n).dot(vector(y,*n));
return vector(x,*n,*incx).dot(vector(y,*n,*incy));
}
// computes a vector-vector dot product without complex conjugation.
Scalar EIGEN_BLAS_FUNC(dotu)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
{
Scalar* x = reinterpret_cast<Scalar*>(px);
Scalar* y = reinterpret_cast<Scalar*>(py);
if(*incx==1 && *incy==1)
return (vector(x,*n).cwise()*vector(y,*n)).sum();
return (vector(x,*n,*incx).cwise()*vector(y,*n,*incy)).sum();
}
#endif // ISCOMPLEX
// computes the Euclidean norm of a vector.
Scalar EIGEN_BLAS_FUNC(nrm2)(int *n, RealScalar *px, int *incx)
{
Scalar* x = reinterpret_cast<Scalar*>(px);
if(*incx==1)
return vector(x,*n).norm();
return vector(x,*n,*incx).norm();
}
int EIGEN_BLAS_FUNC(rot)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pc, RealScalar *ps)
{
Scalar* x = reinterpret_cast<Scalar*>(px);
Scalar* y = reinterpret_cast<Scalar*>(py);
Scalar c = *reinterpret_cast<Scalar*>(pc);
Scalar s = *reinterpret_cast<Scalar*>(ps);
StridedVectorType vx(vector(x,*n,*incx));
StridedVectorType vy(vector(y,*n,*incy));
ei_apply_rotation_in_the_plane(vx, vy, PlanarRotation<Scalar>(c,s));
return 1;
}
int EIGEN_BLAS_FUNC(rotg)(RealScalar *pa, RealScalar *pb, RealScalar *pc, RealScalar *ps)
{
Scalar a = *reinterpret_cast<Scalar*>(pa);
Scalar b = *reinterpret_cast<Scalar*>(pb);
Scalar* c = reinterpret_cast<Scalar*>(pc);
Scalar* s = reinterpret_cast<Scalar*>(ps);
PlanarRotation<Scalar> r;
r.makeGivens(a,b);
*c = r.c();
*s = r.s();
return 1;
}
#if !ISCOMPLEX
/*
// performs rotation of points in the modified plane.
int EIGEN_BLAS_FUNC(rotm)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *param)
{
Scalar* x = reinterpret_cast<Scalar*>(px);
Scalar* y = reinterpret_cast<Scalar*>(py);
// TODO
return 0;
}
// computes the modified parameters for a Givens rotation.
int EIGEN_BLAS_FUNC(rotmg)(RealScalar *d1, RealScalar *d2, RealScalar *x1, RealScalar *x2, RealScalar *param)
{
// TODO
return 0;
}
*/
#endif // !ISCOMPLEX
int EIGEN_BLAS_FUNC(scal)(int *n, RealScalar *px, int *incx, RealScalar *palpha)
{
Scalar* x = reinterpret_cast<Scalar*>(px);
Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
if(*incx==1)
vector(x,*n) *= alpha;
vector(x,*n,*incx) *= alpha;
return 1;
}
int EIGEN_BLAS_FUNC(swap)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
{
int size = IsComplex ? 2* *n : *n;
if(*incx==1 && *incy==1)
vector(py,size).swap(vector(px,size));
else
vector(py,size,*incy).swap(vector(px,size,*incx));
return 1;
}
#if !ISCOMPLEX
RealScalar EIGEN_BLAS_FUNC(casum)(int *n, RealScalar *px, int *incx)
{
Complex* x = reinterpret_cast<Complex*>(px);
if(*incx==1)
return vector(x,*n).cwise().abs().sum();
else
return vector(x,*n,*incx).cwise().abs().sum();
return 1;
}
#endif // ISCOMPLEX

214
blas/level2_impl.h Normal file
View File

@ -0,0 +1,214 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#include "common.h"
int EIGEN_BLAS_FUNC(gemv)(char *opa, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *incb, RealScalar *pbeta, RealScalar *pc, int *incc)
{
Scalar* a = reinterpret_cast<Scalar*>(pa);
Scalar* b = reinterpret_cast<Scalar*>(pb);
Scalar* c = reinterpret_cast<Scalar*>(pc);
Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
Scalar beta = *reinterpret_cast<Scalar*>(pbeta);
if(beta!=Scalar(1))
vector(c, *m, *incc) *= beta;
if(OP(*opa)==NOTR)
if(*incc==1)
vector(c,*m) += alpha * matrix(a,*m,*n,*lda) * vector(b,*n,*incb);
else
vector(c,*m,*incc) += alpha * matrix(a,*m,*n,*lda) * vector(b,*n,*incb);
else if(OP(*opa)==TR)
if(*incb==1)
vector(c,*m,*incc) += alpha * matrix(a,*n,*m,*lda).transpose() * vector(b,*n);
else
vector(c,*m,*incc) += alpha * matrix(a,*n,*m,*lda).transpose() * vector(b,*n,*incb);
else if(OP(*opa)==TR)
if(*incb==1)
vector(c,*m,*incc) += alpha * matrix(a,*n,*m,*lda).adjoint() * vector(b,*n);
else
vector(c,*m,*incc) += alpha * matrix(a,*n,*m,*lda).adjoint() * vector(b,*n,*incb);
else
return 0;
return 1;
}
/*
int EIGEN_BLAS_FUNC(trsv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pa, int *lda, RealScalar *pb, int *incb)
{
typedef void (*functype)(int, const Scalar *, int, Scalar *, int);
functype func[16];
static bool init = false;
if(!init)
{
for(int k=0; k<16; ++k)
func[k] = 0;
// func[NOTR | (UP << 2) | (NUNIT << 3)] = (ei_triangular_solve_vector<Scalar, UpperTriangular|0, false,ColMajor,ColMajor>::run);
// func[TR | (UP << 2) | (NUNIT << 3)] = (ei_triangular_solve_vector<Scalar, UpperTriangular|0, false,RowMajor,ColMajor>::run);
// func[ADJ | (UP << 2) | (NUNIT << 3)] = (ei_triangular_solve_vector<Scalar, UpperTriangular|0, Conj, RowMajor,ColMajor>::run);
//
// func[NOTR | (LO << 2) | (NUNIT << 3)] = (ei_triangular_solve_vector<Scalar, LowerTriangular|0, false,ColMajor,ColMajor>::run);
// func[TR | (LO << 2) | (NUNIT << 3)] = (ei_triangular_solve_vector<Scalar, LowerTriangular|0, false,RowMajor,ColMajor>::run);
// func[ADJ | (LO << 2) | (NUNIT << 3)] = (ei_triangular_solve_vector<Scalar, LowerTriangular|0, Conj, RowMajor,ColMajor>::run);
//
// func[NOTR | (UP << 3) | (UNIT << 3)] = (ei_triangular_solve_vector<Scalar, UpperTriangular|UnitDiagBit,false,ColMajor,ColMajor>::run);
// func[TR | (UP << 2) | (UNIT << 3)] = (ei_triangular_solve_vector<Scalar, UpperTriangular|UnitDiagBit,false,RowMajor,ColMajor>::run);
// func[ADJ | (UP << 2) | (UNIT << 3)] = (ei_triangular_solve_vector<Scalar, UpperTriangular|UnitDiagBit,Conj, RowMajor,ColMajor>::run);
//
// func[NOTR | (LO << 2) | (UNIT << 3)] = (ei_triangular_solve_vector<Scalar, LowerTriangular|UnitDiagBit,false,ColMajor,ColMajor>::run);
// func[TR | (LO << 2) | (UNIT << 3)] = (ei_triangular_solve_vector<Scalar, LowerTriangular|UnitDiagBit,false,RowMajor,ColMajor>::run);
// func[ADJ | (LO << 2) | (UNIT << 3)] = (ei_triangular_solve_vector<Scalar, LowerTriangular|UnitDiagBit,Conj, RowMajor,ColMajor>::run);
init = true;
}
Scalar* a = reinterpret_cast<Scalar*>(pa);
Scalar* b = reinterpret_cast<Scalar*>(pb);
int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
if(code>=16 || func[code]==0)
return 0;
func[code](*n, a, *lda, b, *incb);
return 1;
}
*/
/*
int EIGEN_BLAS_FUNC(trmv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pa, int *lda, RealScalar *pb, int *incb)
{
// TODO
typedef void (*functype)(int, const Scalar *, int, const Scalar *, int, Scalar *, int);
functype func[16];
static bool init = false;
if(!init)
{
for(int k=0; k<16; ++k)
func[k] = 0;
// func[NOTR | (UP << 2) | (NUNIT << 3)] = (ei_product_triangular_matrix_vector<Scalar,UpperTriangular|0, true, ColMajor,false,ColMajor,false,ColMajor>::run);
// func[TR | (UP << 2) | (NUNIT << 3)] = (ei_product_triangular_matrix_vector<Scalar,UpperTriangular|0, true, RowMajor,false,ColMajor,false,ColMajor>::run);
// func[ADJ | (UP << 2) | (NUNIT << 3)] = (ei_product_triangular_matrix_vector<Scalar,UpperTriangular|0, true, RowMajor,Conj, ColMajor,false,ColMajor>::run);
//
// func[NOTR | (LO << 2) | (NUNIT << 3)] = (ei_product_triangular_matrix_vector<Scalar,LowerTriangular|0, true, ColMajor,false,ColMajor,false,ColMajor>::run);
// func[TR | (LO << 2) | (NUNIT << 3)] = (ei_product_triangular_matrix_vector<Scalar,LowerTriangular|0, true, RowMajor,false,ColMajor,false,ColMajor>::run);
// func[ADJ | (LO << 2) | (NUNIT << 3)] = (ei_product_triangular_matrix_vector<Scalar,LowerTriangular|0, true, RowMajor,Conj, ColMajor,false,ColMajor>::run);
//
// func[NOTR | (UP << 2) | (UNIT << 3)] = (ei_product_triangular_matrix_vector<Scalar,UpperTriangular|UnitDiagBit,true, ColMajor,false,ColMajor,false,ColMajor>::run);
// func[TR | (UP << 2) | (UNIT << 3)] = (ei_product_triangular_matrix_vector<Scalar,UpperTriangular|UnitDiagBit,true, RowMajor,false,ColMajor,false,ColMajor>::run);
// func[ADJ | (UP << 2) | (UNIT << 3)] = (ei_product_triangular_matrix_vector<Scalar,UpperTriangular|UnitDiagBit,true, RowMajor,Conj, ColMajor,false,ColMajor>::run);
//
// func[NOTR | (LO << 2) | (UNIT << 3)] = (ei_product_triangular_matrix_vector<Scalar,LowerTriangular|UnitDiagBit,true, ColMajor,false,ColMajor,false,ColMajor>::run);
// func[TR | (LO << 2) | (UNIT << 3)] = (ei_product_triangular_matrix_vector<Scalar,LowerTriangular|UnitDiagBit,true, RowMajor,false,ColMajor,false,ColMajor>::run);
// func[ADJ | (LO << 2) | (UNIT << 3)] = (ei_product_triangular_matrix_vector<Scalar,LowerTriangular|UnitDiagBit,true, RowMajor,Conj, ColMajor,false,ColMajor>::run);
init = true;
}
Scalar* a = reinterpret_cast<Scalar*>(pa);
Scalar* b = reinterpret_cast<Scalar*>(pb);
int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
if(code>=16 || func[code]==0)
return 0;
func[code](*n, a, *lda, b, *incb, b, *incb);
return 1;
}
*/
/*
int EIGEN_BLAS_FUNC(syr)(char *uplo, int *n, RealScalar *palpha, RealScalar *pa, int *inca, RealScalar *pc, int *ldc)
{
// TODO
typedef void (*functype)(int, const Scalar *, int, Scalar *, int, Scalar);
functype func[2];
static bool init = false;
if(!init)
{
for(int k=0; k<2; ++k)
func[k] = 0;
// func[UP] = (ei_selfadjoint_product<Scalar,ColMajor,ColMajor,false,UpperTriangular>::run);
// func[LO] = (ei_selfadjoint_product<Scalar,ColMajor,ColMajor,false,LowerTriangular>::run);
init = true;
}
Scalar* a = reinterpret_cast<Scalar*>(pa);
Scalar* c = reinterpret_cast<Scalar*>(pc);
Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
int code = UPLO(*uplo);
if(code>=2 || func[code]==0)
return 0;
func[code](*n, a, *inca, c, *ldc, alpha);
return 1;
}
*/
/*
int EIGEN_BLAS_FUNC(syr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *pa, int *inca, RealScalar *pb, int *incb, RealScalar *pc, int *ldc)
{
// TODO
typedef void (*functype)(int, const Scalar *, int, const Scalar *, int, Scalar *, int, Scalar);
functype func[2];
static bool init = false;
if(!init)
{
for(int k=0; k<2; ++k)
func[k] = 0;
// func[UP] = (ei_selfadjoint_product<Scalar,ColMajor,ColMajor,false,UpperTriangular>::run);
// func[LO] = (ei_selfadjoint_product<Scalar,ColMajor,ColMajor,false,LowerTriangular>::run);
init = true;
}
Scalar* a = reinterpret_cast<Scalar*>(pa);
Scalar* b = reinterpret_cast<Scalar*>(pb);
Scalar* c = reinterpret_cast<Scalar*>(pc);
Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
int code = UPLO(*uplo);
if(code>=2 || func[code]==0)
return 0;
func[code](*n, a, *inca, b, *incb, c, *ldc, alpha);
return 1;
}
*/
#if ISCOMPLEX
#endif // ISCOMPLEX

365
blas/level3_impl.h Normal file
View File

@ -0,0 +1,365 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#include "common.h"
int EIGEN_BLAS_FUNC(gemm)(char *opa, char *opb, int *m, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc)
{
typedef void (*functype)(int, int, int, const Scalar *, int, const Scalar *, int, Scalar *, int, Scalar);
functype func[12];
static bool init = false;
if(!init)
{
for(int k=0; k<12; ++k)
func[k] = 0;
func[NOTR | (NOTR << 2)] = (ei_general_matrix_matrix_product<Scalar,ColMajor,false,ColMajor,false,ColMajor>::run);
func[TR | (NOTR << 2)] = (ei_general_matrix_matrix_product<Scalar,RowMajor,false,ColMajor,false,ColMajor>::run);
func[ADJ | (NOTR << 2)] = (ei_general_matrix_matrix_product<Scalar,RowMajor,Conj, ColMajor,false,ColMajor>::run);
func[NOTR | (TR << 2)] = (ei_general_matrix_matrix_product<Scalar,ColMajor,false,RowMajor,false,ColMajor>::run);
func[TR | (TR << 2)] = (ei_general_matrix_matrix_product<Scalar,RowMajor,false,RowMajor,false,ColMajor>::run);
func[ADJ | (TR << 2)] = (ei_general_matrix_matrix_product<Scalar,RowMajor,Conj, RowMajor,false,ColMajor>::run);
func[NOTR | (ADJ << 2)] = (ei_general_matrix_matrix_product<Scalar,ColMajor,false,RowMajor,Conj, ColMajor>::run);
func[TR | (ADJ << 2)] = (ei_general_matrix_matrix_product<Scalar,RowMajor,false,RowMajor,Conj, ColMajor>::run);
func[ADJ | (ADJ << 2)] = (ei_general_matrix_matrix_product<Scalar,RowMajor,Conj, RowMajor,Conj, ColMajor>::run);
init = true;
}
Scalar* a = reinterpret_cast<Scalar*>(pa);
Scalar* b = reinterpret_cast<Scalar*>(pb);
Scalar* c = reinterpret_cast<Scalar*>(pc);
Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
Scalar beta = *reinterpret_cast<Scalar*>(pbeta);
if(beta!=Scalar(1))
matrix(c, *m, *n, *ldc) *= beta;
int code = OP(*opa) | (OP(*opb) << 2);
if(code>=12 || func[code]==0)
return 0;
func[code](*m, *n, *k, a, *lda, b, *ldb, c, *ldc, alpha);
return 1;
}
int EIGEN_BLAS_FUNC(trsm)(char *side, char *uplo, char *opa, char *diag, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb)
{
typedef void (*functype)(int, int, const Scalar *, int, Scalar *, int);
functype func[32];
static bool init = false;
if(!init)
{
for(int k=0; k<32; ++k)
func[k] = 0;
func[NOTR | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, UpperTriangular|0, false,ColMajor,ColMajor>::run);
func[TR | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, UpperTriangular|0, false,RowMajor,ColMajor>::run);
func[ADJ | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, UpperTriangular|0, Conj, RowMajor,ColMajor>::run);
func[NOTR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,UpperTriangular|0, false,ColMajor,ColMajor>::run);
func[TR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,UpperTriangular|0, false,RowMajor,ColMajor>::run);
func[ADJ | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,UpperTriangular|0, Conj, RowMajor,ColMajor>::run);
func[NOTR | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, LowerTriangular|0, false,ColMajor,ColMajor>::run);
func[TR | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, LowerTriangular|0, false,RowMajor,ColMajor>::run);
func[ADJ | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, LowerTriangular|0, Conj, RowMajor,ColMajor>::run);
func[NOTR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,LowerTriangular|0, false,ColMajor,ColMajor>::run);
func[TR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,LowerTriangular|0, false,RowMajor,ColMajor>::run);
func[ADJ | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,LowerTriangular|0, Conj, RowMajor,ColMajor>::run);
func[NOTR | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, UpperTriangular|UnitDiagBit,false,ColMajor,ColMajor>::run);
func[TR | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, UpperTriangular|UnitDiagBit,false,RowMajor,ColMajor>::run);
func[ADJ | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, UpperTriangular|UnitDiagBit,Conj, RowMajor,ColMajor>::run);
func[NOTR | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,UpperTriangular|UnitDiagBit,false,ColMajor,ColMajor>::run);
func[TR | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,UpperTriangular|UnitDiagBit,false,RowMajor,ColMajor>::run);
func[ADJ | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,UpperTriangular|UnitDiagBit,Conj, RowMajor,ColMajor>::run);
func[NOTR | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, LowerTriangular|UnitDiagBit,false,ColMajor,ColMajor>::run);
func[TR | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, LowerTriangular|UnitDiagBit,false,RowMajor,ColMajor>::run);
func[ADJ | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, LowerTriangular|UnitDiagBit,Conj, RowMajor,ColMajor>::run);
func[NOTR | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,LowerTriangular|UnitDiagBit,false,ColMajor,ColMajor>::run);
func[TR | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,LowerTriangular|UnitDiagBit,false,RowMajor,ColMajor>::run);
func[ADJ | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,LowerTriangular|UnitDiagBit,Conj, RowMajor,ColMajor>::run);
init = true;
}
Scalar* a = reinterpret_cast<Scalar*>(pa);
Scalar* b = reinterpret_cast<Scalar*>(pb);
Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
// TODO handle alpha
int code = OP(*opa) | (SIDE(*side) << 2) | (UPLO(*uplo) << 3) | (DIAG(*diag) << 4);
if(code>=32 || func[code]==0)
return 0;
func[code](*m, *n, a, *lda, b, *ldb);
return 1;
}
// b = alpha*op(a)*b for side = 'L'or'l'
// b = alpha*b*op(a) for side = 'R'or'r'
int EIGEN_BLAS_FUNC(trmm)(char *side, char *uplo, char *opa, char *diag, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb)
{
typedef void (*functype)(int, int, const Scalar *, int, const Scalar *, int, Scalar *, int, Scalar);
functype func[32];
static bool init = false;
if(!init)
{
for(int k=0; k<32; ++k)
func[k] = 0;
func[NOTR | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,UpperTriangular|0, true, ColMajor,false,ColMajor,false,ColMajor>::run);
func[TR | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,UpperTriangular|0, true, RowMajor,false,ColMajor,false,ColMajor>::run);
func[ADJ | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,UpperTriangular|0, true, RowMajor,Conj, ColMajor,false,ColMajor>::run);
func[NOTR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,UpperTriangular|0, false,ColMajor,false,ColMajor,false,ColMajor>::run);
func[TR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,UpperTriangular|0, false,ColMajor,false,RowMajor,false,ColMajor>::run);
func[ADJ | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,UpperTriangular|0, false,ColMajor,false,RowMajor,Conj, ColMajor>::run);
func[NOTR | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,LowerTriangular|0, true, ColMajor,false,ColMajor,false,ColMajor>::run);
func[TR | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,LowerTriangular|0, true, RowMajor,false,ColMajor,false,ColMajor>::run);
func[ADJ | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,LowerTriangular|0, true, RowMajor,Conj, ColMajor,false,ColMajor>::run);
func[NOTR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,LowerTriangular|0, false,ColMajor,false,ColMajor,false,ColMajor>::run);
func[TR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,LowerTriangular|0, false,ColMajor,false,RowMajor,false,ColMajor>::run);
func[ADJ | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,LowerTriangular|0, false,ColMajor,false,RowMajor,Conj, ColMajor>::run);
func[NOTR | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,UpperTriangular|UnitDiagBit,true, ColMajor,false,ColMajor,false,ColMajor>::run);
func[TR | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,UpperTriangular|UnitDiagBit,true, RowMajor,false,ColMajor,false,ColMajor>::run);
func[ADJ | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,UpperTriangular|UnitDiagBit,true, RowMajor,Conj, ColMajor,false,ColMajor>::run);
func[NOTR | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,UpperTriangular|UnitDiagBit,false,ColMajor,false,ColMajor,false,ColMajor>::run);
func[TR | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,UpperTriangular|UnitDiagBit,false,ColMajor,false,RowMajor,false,ColMajor>::run);
func[ADJ | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,UpperTriangular|UnitDiagBit,false,ColMajor,false,RowMajor,Conj, ColMajor>::run);
func[NOTR | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,LowerTriangular|UnitDiagBit,true, ColMajor,false,ColMajor,false,ColMajor>::run);
func[TR | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,LowerTriangular|UnitDiagBit,true, RowMajor,false,ColMajor,false,ColMajor>::run);
func[ADJ | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,LowerTriangular|UnitDiagBit,true, RowMajor,Conj, ColMajor,false,ColMajor>::run);
func[NOTR | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,LowerTriangular|UnitDiagBit,false,ColMajor,false,ColMajor,false,ColMajor>::run);
func[TR | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,LowerTriangular|UnitDiagBit,false,ColMajor,false,RowMajor,false,ColMajor>::run);
func[ADJ | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,LowerTriangular|UnitDiagBit,false,ColMajor,false,RowMajor,Conj, ColMajor>::run);
init = true;
}
Scalar* a = reinterpret_cast<Scalar*>(pa);
Scalar* b = reinterpret_cast<Scalar*>(pb);
Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
int code = OP(*opa) | (SIDE(*side) << 2) | (UPLO(*uplo) << 3) | (DIAG(*diag) << 4);
if(code>=32 || func[code]==0)
return 0;
func[code](*m, *n, a, *lda, b, *ldb, b, *ldb, alpha);
return 1;
}
// c = alpha*a*b + beta*c for side = 'L'or'l'
// c = alpha*b*a + beta*c for side = 'R'or'r
int EIGEN_BLAS_FUNC(symm)(char *side, char *uplo, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc)
{
Scalar* a = reinterpret_cast<Scalar*>(pa);
Scalar* b = reinterpret_cast<Scalar*>(pb);
Scalar* c = reinterpret_cast<Scalar*>(pc);
Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
Scalar beta = *reinterpret_cast<Scalar*>(pbeta);
if(beta!=Scalar(1))
matrix(c, *m, *n, *ldc) *= beta;
if(SIDE(*side)==LEFT)
if(UPLO(*uplo)==UP)
ei_product_selfadjoint_matrix<Scalar, RowMajor,true,false, ColMajor,false,false, ColMajor>::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha);
else if(UPLO(*uplo)==LO)
ei_product_selfadjoint_matrix<Scalar, ColMajor,true,false, ColMajor,false,false, ColMajor>::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha);
else
return 0;
else if(SIDE(*side)==RIGHT)
if(UPLO(*uplo)==UP)
ei_product_selfadjoint_matrix<Scalar, ColMajor,false,false, RowMajor,true,false, ColMajor>::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha);
else if(UPLO(*uplo)==LO)
ei_product_selfadjoint_matrix<Scalar, ColMajor,false,false, ColMajor,true,false, ColMajor>::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha);
else
return 0;
else
return 0;
return 1;
}
// c = alpha*a*a' + beta*c for op = 'N'or'n'
// c = alpha*a'*a + beta*c for op = 'T'or't','C'or'c'
int EIGEN_BLAS_FUNC(syrk)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pbeta, RealScalar *pc, int *ldc)
{
typedef void (*functype)(int, int, const Scalar *, int, Scalar *, int, Scalar);
functype func[8];
static bool init = false;
if(!init)
{
for(int k=0; k<8; ++k)
func[k] = 0;
func[NOTR | (UP << 2)] = (ei_selfadjoint_product<Scalar,ColMajor,ColMajor,true, UpperTriangular>::run);
func[TR | (UP << 2)] = (ei_selfadjoint_product<Scalar,RowMajor,ColMajor,false,UpperTriangular>::run);
func[ADJ | (UP << 2)] = (ei_selfadjoint_product<Scalar,RowMajor,ColMajor,false,UpperTriangular>::run);
func[NOTR | (LO << 2)] = (ei_selfadjoint_product<Scalar,ColMajor,ColMajor,true, LowerTriangular>::run);
func[TR | (LO << 2)] = (ei_selfadjoint_product<Scalar,RowMajor,ColMajor,false,LowerTriangular>::run);
func[ADJ | (LO << 2)] = (ei_selfadjoint_product<Scalar,RowMajor,ColMajor,false,LowerTriangular>::run);
init = true;
}
Scalar* a = reinterpret_cast<Scalar*>(pa);
Scalar* c = reinterpret_cast<Scalar*>(pc);
Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
Scalar beta = *reinterpret_cast<Scalar*>(pbeta);
int code = OP(*op) | (UPLO(*uplo) << 2);
if(code>=8 || func[code]==0)
return 0;
if(beta!=Scalar(1))
matrix(c, *n, *n, *ldc) *= beta;
func[code](*n, *k, a, *lda, c, *ldc, alpha);
return 1;
}
// c = alpha*a*b' + alpha*b*a' + beta*c for op = 'N'or'n'
// c = alpha*a'*b + alpha*b'*a + beta*c for op = 'T'or't'
int EIGEN_BLAS_FUNC(syr2k)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc)
{
Scalar* a = reinterpret_cast<Scalar*>(pa);
Scalar* b = reinterpret_cast<Scalar*>(pb);
Scalar* c = reinterpret_cast<Scalar*>(pc);
Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
Scalar beta = *reinterpret_cast<Scalar*>(pbeta);
// TODO
return 0;
}
#if ISCOMPLEX
// c = alpha*a*b + beta*c for side = 'L'or'l'
// c = alpha*b*a + beta*c for side = 'R'or'r
int EIGEN_BLAS_FUNC(hemm)(char *side, char *uplo, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc)
{
Scalar* a = reinterpret_cast<Scalar*>(pa);
Scalar* b = reinterpret_cast<Scalar*>(pb);
Scalar* c = reinterpret_cast<Scalar*>(pc);
Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
Scalar beta = *reinterpret_cast<Scalar*>(pbeta);
if(beta!=Scalar(1))
matrix(c, *m, *n, *ldc) *= beta;
if(SIDE(*side)==LEFT)
if(UPLO(*uplo)==UP)
ei_product_selfadjoint_matrix<Scalar, RowMajor,true,Conj, ColMajor,false,false, ColMajor>::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha);
else if(UPLO(*uplo)==LO)
ei_product_selfadjoint_matrix<Scalar, ColMajor,true,false, ColMajor,false,false, ColMajor>::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha);
else
return 0;
else if(SIDE(*side)==RIGHT)
if(UPLO(*uplo)==UP)
ei_product_selfadjoint_matrix<Scalar, ColMajor,false,false, RowMajor,true,Conj, ColMajor>::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha);
else if(UPLO(*uplo)==LO)
ei_product_selfadjoint_matrix<Scalar, ColMajor,false,false, ColMajor,true,false, ColMajor>::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha);
else
return 0;
else
return 0;
return 1;
}
// c = alpha*a*conj(a') + beta*c for op = 'N'or'n'
// c = alpha*conj(a')*a + beta*c for op = 'C'or'c'
int EIGEN_BLAS_FUNC(herk)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pbeta, RealScalar *pc, int *ldc)
{
typedef void (*functype)(int, int, const Scalar *, int, Scalar *, int, Scalar);
functype func[8];
static bool init = false;
if(!init)
{
for(int k=0; k<8; ++k)
func[k] = 0;
func[NOTR | (UP << 2)] = (ei_selfadjoint_product<Scalar,ColMajor,ColMajor,true, UpperTriangular>::run);
func[ADJ | (UP << 2)] = (ei_selfadjoint_product<Scalar,RowMajor,ColMajor,false,UpperTriangular>::run);
func[NOTR | (LO << 2)] = (ei_selfadjoint_product<Scalar,ColMajor,ColMajor,true, LowerTriangular>::run);
func[ADJ | (LO << 2)] = (ei_selfadjoint_product<Scalar,RowMajor,ColMajor,false,LowerTriangular>::run);
init = true;
}
Scalar* a = reinterpret_cast<Scalar*>(pa);
Scalar* c = reinterpret_cast<Scalar*>(pc);
Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
Scalar beta = *reinterpret_cast<Scalar*>(pbeta);
int code = OP(*op) | (UPLO(*uplo) << 2);
if(code>=8 || func[code]==0)
return 0;
if(beta!=Scalar(1))
matrix(c, *n, *n, *ldc) *= beta;
func[code](*n, *k, a, *lda, c, *ldc, alpha);
return 1;
}
// c = alpha*a*conj(b') + conj(alpha)*b*conj(a') + beta*c, for op = 'N'or'n'
// c = alpha*conj(b')*a + conj(alpha)*conj(a')*b + beta*c, for op = 'C'or'c'
int EIGEN_BLAS_FUNC(her2k)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc)
{
Scalar* a = reinterpret_cast<Scalar*>(pa);
Scalar* b = reinterpret_cast<Scalar*>(pb);
Scalar* c = reinterpret_cast<Scalar*>(pc);
Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
Scalar beta = *reinterpret_cast<Scalar*>(pbeta);
// TODO
return 0;
}
#endif // ISCOMPLEX

31
blas/single.cpp Normal file
View File

@ -0,0 +1,31 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#define SCALAR float
#define SCALAR_SUFFIX s
#define ISCOMPLEX 0
#include "level1_impl.h"
#include "level2_impl.h"
#include "level3_impl.h"

65
cmake/EigenTesting.cmake
View File

@ -27,7 +27,11 @@ endmacro(ei_add_property)
# void test_<testname>() { ... }
#
# this macro add an executable test_<testname> as well as a ctest test
# named <testname>
# named <testname>.
#
# it also adds another executable debug_<testname> that compiles in full debug mode
# and is not added to the test target. The idea is that when a test fails you want
# a quick way of rebuilding this specific test in full debug mode.
#
# On platforms with bash simply run:
# "ctest -V" or "ctest -V -R <testname>"
@ -36,38 +40,65 @@ endmacro(ei_add_property)
macro(ei_add_test testname)
set(targetname test_${testname})
if(NOT MSVC_IDE)
set(debug_targetname debug_${testname})
endif(NOT MSVC_IDE)
set(filename ${testname}.cpp)
add_executable(${targetname} ${filename})
add_dependencies(btest ${targetname})
if(NOT MSVC_IDE)
add_executable(${debug_targetname} EXCLUDE_FROM_ALL ${filename})
endif(NOT MSVC_IDE)
if(NOT EIGEN_NO_ASSERTION_CHECKING)
if(MSVC)
set_target_properties(${targetname} PROPERTIES COMPILE_FLAGS "/EHsc")
if(NOT MSVC_IDE)
set_target_properties(${debug_targetname} PROPERTIES COMPILE_FLAGS "/EHsc")
endif(NOT MSVC_IDE)
else(MSVC)
set_target_properties(${targetname} PROPERTIES COMPILE_FLAGS "-fexceptions")
set_target_properties(${debug_targetname} PROPERTIES COMPILE_FLAGS "-fexceptions")
endif(MSVC)
option(EIGEN_DEBUG_ASSERTS "Enable debuging of assertions" OFF)
if(EIGEN_DEBUG_ASSERTS)
ei_add_target_property(${targetname} COMPILE_FLAGS "-DEIGEN_DEBUG_ASSERTS=1")
if(NOT MSVC_IDE)
ei_add_target_property(${debug_targetname} COMPILE_FLAGS "-DEIGEN_DEBUG_ASSERTS=1")
endif(NOT MSVC_IDE)
endif(EIGEN_DEBUG_ASSERTS)
else(NOT EIGEN_NO_ASSERTION_CHECKING)
ei_add_target_property(${targetname} COMPILE_FLAGS "-DEIGEN_NO_ASSERTION_CHECKING=1")
if(NOT MSVC_IDE)
ei_add_target_property(${debug_targetname} COMPILE_FLAGS "-DEIGEN_NO_ASSERTION_CHECKING=1")
endif(NOT MSVC_IDE)
endif(NOT EIGEN_NO_ASSERTION_CHECKING)
# let the user pass e.g. optimization flags, but don't apply them to the debug target
if(${ARGC} GREATER 1)
ei_add_target_property(${targetname} COMPILE_FLAGS "${ARGV1}")
endif(${ARGC} GREATER 1)
ei_add_target_property(${targetname} COMPILE_FLAGS "-DEIGEN_TEST_FUNC=${testname}")
# for the debug target, add full debug options
if(CMAKE_COMPILER_IS_GNUCXX)
# O0 is in principle redundant here, but doesn't hurt
ei_add_target_property(${debug_targetname} COMPILE_FLAGS "-O0 -g3")
elseif(MSVC)
if(NOT MSVC_IDE)
ei_add_target_property(${debug_targetname} COMPILE_FLAGS "/Od /Zi")
endif(NOT MSVC_IDE)
endif(CMAKE_COMPILER_IS_GNUCXX)
if(TEST_LIB)
target_link_libraries(${targetname} Eigen2)
endif(TEST_LIB)
ei_add_target_property(${targetname} COMPILE_FLAGS "-DEIGEN_TEST_FUNC=${testname}")
if(NOT MSVC_IDE)
ei_add_target_property(${debug_targetname} COMPILE_FLAGS "-DEIGEN_TEST_FUNC=${testname}")
endif(NOT MSVC_IDE)
target_link_libraries(${targetname} ${EXTERNAL_LIBS})
if(${ARGC} GREATER 2)
@ -75,11 +106,18 @@ macro(ei_add_test testname)
string(LENGTH "${ARGV2_stripped}" ARGV2_stripped_length)
if(${ARGV2_stripped_length} GREATER 0)
target_link_libraries(${targetname} ${ARGV2})
if(NOT MSVC_IDE)
target_link_libraries(${debug_targetname} ${ARGV2})
endif(NOT MSVC_IDE)
endif(${ARGV2_stripped_length} GREATER 0)
endif(${ARGC} GREATER 2)
if(WIN32)
add_test(${testname} "${targetname}")
if(CYGWIN)
add_test(${testname} "${Eigen_SOURCE_DIR}/test/runtest.sh" "${testname}")
else(CYGWIN)
add_test(${testname} "${targetname}")
endif(CYGWIN)
else(WIN32)
add_test(${testname} "${Eigen_SOURCE_DIR}/test/runtest.sh" "${testname}")
endif(WIN32)
@ -102,31 +140,31 @@ macro(ei_testing_print_summary)
if(EIGEN_TEST_SSE2)
message("SSE2: ON")
else(EIGEN_TEST_SSE2)
message("SSE2: AUTO")
message("SSE2: Using architecture defaults")
endif(EIGEN_TEST_SSE2)
if(EIGEN_TEST_SSE3)
message("SSE3: ON")
else(EIGEN_TEST_SSE3)
message("SSE3: AUTO")
message("SSE3: Using architecture defaults")
endif(EIGEN_TEST_SSE3)
if(EIGEN_TEST_SSSE3)
message("SSSE3: ON")
else(EIGEN_TEST_SSSE3)
message("SSSE3: AUTO")
message("SSSE3: Using architecture defaults")
endif(EIGEN_TEST_SSSE3)
if(EIGEN_TEST_ALTIVEC)
message("Altivec: ON")
message("Altivec: Using architecture defaults")
else(EIGEN_TEST_ALTIVEC)
message("Altivec: AUTO")
message("Altivec: Using architecture defaults")
endif(EIGEN_TEST_ALTIVEC)
if(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION)
message("Explicit vec: OFF")
else(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION)
message("Explicit vec: AUTO")
message("Explicit vec: Using architecture defaults")
endif(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION)
message("\n${EIGEN_TESTING_SUMMARY}")
@ -159,6 +197,9 @@ if(CMAKE_COMPILER_IS_GNUCXX)
else(EIGEN_COVERAGE_TESTING)
set(COVERAGE_FLAGS "")
endif(EIGEN_COVERAGE_TESTING)
if(EIGEN_TEST_RVALUE_REF_SUPPORT OR EIGEN_TEST_C++0x)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++0x")
endif(EIGEN_TEST_RVALUE_REF_SUPPORT OR EIGEN_TEST_C++0x)
if(CMAKE_SYSTEM_NAME MATCHES Linux)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${COVERAGE_FLAGS} -g2")
set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "${CMAKE_CXX_FLAGS_RELWITHDEBINFO} ${COVERAGE_FLAGS} -O2 -g2")

12
demos/CMakeLists.txt
View File

@ -1,3 +1,11 @@
project(EigenDemos)
add_subdirectory(mandelbrot)
add_subdirectory(opengl)
add_custom_target(demos)
find_package(Qt4)
if(QT4_FOUND)
add_subdirectory(mandelbrot)
add_subdirectory(opengl)
else(QT4_FOUND)
message(STATUS "Qt4 not found, so disabling the mandelbrot and opengl demos")
endif(QT4_FOUND)

1
demos/mandelbrot/CMakeLists.txt
View File

@ -16,5 +16,6 @@ set(mandelbrot_SRCS
qt4_automoc(${mandelbrot_SRCS})
add_executable(mandelbrot ${mandelbrot_SRCS})
add_dependencies(demos mandelbrot)
target_link_libraries(mandelbrot ${QT_QTCORE_LIBRARY} ${QT_QTGUI_LIBRARY})

4
demos/mix_eigen_and_c/README
View File

@ -1,7 +1,9 @@
This is an example of how one can wrap some of Eigen into a C library.
To try this with GCC, do:
g++ -c binary_library.cpp -O2 -msse2 -I ../..
gcc example.c binary_library.o -o example -lstdc++
./example
This is an example of how one can wrap some of Eigen into a C library.
TODO: add CMakeLists, add more explanations here

2
demos/opengl/CMakeLists.txt
View File

@ -1,4 +1,3 @@
find_package(Qt4 REQUIRED)
find_package(OpenGL REQUIRED)
@ -14,6 +13,7 @@ set(quaternion_demo_SRCS gpuhelper.cpp icosphere.cpp camera.cpp trackball.cpp q
qt4_automoc(${quaternion_demo_SRCS})
add_executable(quaternion_demo ${quaternion_demo_SRCS})
add_dependencies(demos quaternion_demo)
target_link_libraries(quaternion_demo
${QT_QTCORE_LIBRARY} ${QT_QTGUI_LIBRARY}

2
demos/opengl/camera.cpp
View File

@ -260,7 +260,7 @@ void Camera::activateGL(void)
Vector3f Camera::unProject(const Vector2f& uv, float depth) const
{
Matrix4f inv = mViewMatrix.inverse();
Matrix4f inv = mViewMatrix.inverse().matrix();
return unProject(uv, depth, inv);
}

3
doc/AsciiQuickReference.txt
View File

@ -118,6 +118,9 @@ s = R.maxCoeff(&r, &c) // [aa, bb] = max(R); [cc, dd] = max(aa);
R.sum() // sum(R(:))
R.colwise.sum() // sum(R)
R.rowwise.sum() // sum(R, 2) or sum(R')'
R.prod() // prod(R(:))
R.colwise.prod() // prod(R)
R.rowwise.prod() // prod(R, 2) or prod(R')'
R.trace() // trace(R)
R.all() // all(R(:))
R.colwise().all() // all(R)

67
doc/C07_TutorialSparse.dox
View File

@ -129,7 +129,7 @@ SparseVector<double> vec(size);
for (SparseVector<double>::InnerIterator it(vec); it; ++it)
{
it.value(); // == vec[ it.index() ]
it.index();
it.index();
}
\endcode
</td></tr>
@ -138,58 +138,51 @@ for (SparseVector<double>::InnerIterator it(vec); it; ++it)
\section TutorialSparseFilling Filling a sparse matrix
A DynamicSparseMatrix object can be set and updated just like any dense matrix using the coeffRef(row,col) method. If the coefficient is not stored yet, then it will be inserted in the matrix. Here is an example:
Owing to the special storage scheme of a SparseMatrix, it is obvious that for performance reasons a sparse matrix cannot be filled as easily as a dense matrix. For instance the cost of a purely random insertion into a SparseMatrix is in O(nnz) where nnz is the current number of non zeros. In order to cover all uses cases with best efficiency, Eigen provides various mechanisms, from the easiest but slowest, to the fastest but restrictive one.
If you don't have any prior knowledge about the order your matrix will be filled, then the best choice is to use a DynamicSparseMatrix. With a DynamicSparseMatrix, you can add or modify any coefficients at any time using the coeffRef(row,col) method. Here is an example:
\code
DynamicSparseMatrix<float> aux(1000,1000);
aux.reserve(estimated_number_of_non_zero); // optional
for (...)
for each i
for each j interacting with i
aux.coeffRef(i,j) += foo(o1,o2);
SparseMatrix<float> mat(aux); // convert the DynamicSparseMatrix to a SparseMatrix
for each j // the j can be random
for each i interacting with j // the i can be random
aux.coeffRef(i,j) += foo(i,j);
\endcode
Then the DynamicSparseMatrix object can be converted to a compact SparseMatrix to be used, e.g., by one of our supported solver:
\code
SparseMatrix<float> mat(aux);
\endcode
Sometimes, however, we simply want to set all the coefficients of a matrix before using it through standard matrix operations (addition, product, etc.). In that case it faster to use the low-level startFill()/fill()/fillrand()/endFill() interface. Even though this interface is availabe for both sparse matrix types, their respective restrictions slightly differ from one representation to the other. In all case, a call to startFill() set the matrix to zero, and the fill*() functions will fail if the coefficient already exist.
In order to optimize this process, instead of the generic coeffRef(i,j) method one can also use:
- \code m.insert(i,j) = value; \endcode which assumes the coefficient of coordinate (row,col) does not already exist (otherwise this is a programming error and your program will stop).
- \code m.insertBack(i,j) = value; \endcode which, in addition to the requirements of insert(), also assumes that the coefficient of coordinate (row,col) will be inserted at the end of the target inner-vector. More precisely, if the matrix m is column major, then the row index of the last non zero coefficient of the j-th column must be smaller than i.
As a first difference, for SparseMatrix, the fill*() functions can only be called inside a startFill()/endFill() pair, and no other member functions are allowed during the filling process, i.e., until endFill() has been called. On the other hand, a DynamicSparseMatrix is always in a stable state, and the startFill()/endFill() functions are only for compatibility purpose.
Another difference is that the fill*() functions must be called with increasing outer indices for a SparseMatrix, while they can be random for a DynamicSparseMatrix.
Actually, the SparseMatrix class also supports random insertion via the insert() method. However, its uses should be reserved in cases where the inserted non zero is nearly the last one of the compact storage array. In practice, this means it should be used only to perform random (or sorted) insertion into the current inner-vector while filling the inner-vectors in an increasing order. Moreover, with a SparseMatrix an insertion session must be closed by a call to finalize() before any use of the matrix. Here is an example for a column major matrix:
Finally, the fill() function assumes the coefficient are inserted in a sorted order per inner vector, while the fillrand() variante allows random insertions (the outer indices must still be sorted for SparseMatrix).
Some examples:
1 - If you can set the coefficients in exactly the same order that the storage order, then the matrix can be filled directly and very efficiently. Here is an example initializing a random, row-major sparse matrix:
\code
SparseMatrix<double,RowMajor> m(rows,cols);
m.startFill(rows*cols*percent_of_non_zero); // estimate of the number of nonzeros (optional)
for (int i=0; i\<rows; ++i)
for (int j=0; j\<cols; ++j)
if (rand()\<percent_of_non_zero)
m.fill(i,j) = rand();
m.endFill();
SparseMatrix<float> mat(1000,1000);
mat.reserve(estimated_number_of_non_zero); // optional
for each j // should be in increasing order for performance reasons
for each i interacting with j // the i can be random
mat.insert(i,j) = foo(i,j); // optional for a DynamicSparseMatrix
mat.finalize();
\endcode
2 - If you can set each outer vector in a consistent order, but do not have sorted data for each inner vector, then you can use fillrand() instead of fill():
\code
SparseMatrix<double,RowMajor> m(rows,cols);
m.startFill(rows*cols*percent_of_non_zero); // estimate of the number of nonzeros (optional)
for (int i=0; i\<rows; ++i)
for (int k=0; k\<cols*percent_of_non_zero; ++k)
m.fillrand(i,rand(0,cols)) = rand();
m.endFill();
\endcode
The fillrand() function performs a sorted insertion into an array sequentially stored in memory and requires a copy of all coefficients stored after its target position. This method is therefore reserved for matrices having only a few elements per row/column (up to 50) and works better if the insertion are almost sorted.
Finally, the fastest way to fill a SparseMatrix object is to insert the elements in a purely coherence order (increasing inner index per increasing outer index). To this end, Eigen provides a very low but optimal API and illustrated below:
3 - Eventually, if none of the above solution is practicable for you, then you have to use a RandomSetter which temporarily wraps the matrix into a more flexible hash map allowing complete random accesses:
\code
SparseMatrix<double,RowMajor> m(rows,cols);
SparseMatrix<float> mat(1000,1000);
mat.reserve(estimated_number_of_non_zero); // optional
for(int j=0; j<1000; ++j)
{
RandomSetter<SparseMatrix<double,RowMajor> > setter(m);
for (int k=0; k\<cols*rows*percent_of_non_zero; ++k)
setter(rand(0,rows), rand(0,cols)) = rand();
mat.startVec(j); // optional for a DynamicSparseMatrix
for each i interacting with j // with increasing i
mat.insertBack(i,j) = foo(i,j);
}
mat.finalize(); // optional for a DynamicSparseMatrix
\endcode
The matrix \c m is set at the destruction of the setter, hence the use of a nested block. This imposed syntax has the advantage to emphasize the critical section where m is not valid and cannot be used.
\section TutorialSparseFeatureSet Supported operators and functions

1
doc/CMakeLists.txt
View File

@ -1,3 +1,4 @@
project(EigenDoc)
set_directory_properties(PROPERTIES EXCLUDE_FROM_ALL TRUE)

2
doc/I00_CustomizingEigen.dox
View File

@ -125,7 +125,7 @@ In order to add support for a custom type \c T you need:
3 - define a couple of math functions for your type such as: ei_sqrt, ei_abs, etc...
(see the file Eigen/src/Core/MathFunctions.h)
Here is a concrete example adding support for the Adolc's \c adouble type. <a href="http://www.math.tu-dresden.de/~adol-c/">Adolc</a> is an automatic differentiation library. The type \c adouble is basically a real value tracking the values of any number of partial derivatives.
Here is a concrete example adding support for the Adolc's \c adouble type. <a href="https://projects.coin-or.org/ADOL-C">Adolc</a> is an automatic differentiation library. The type \c adouble is basically a real value tracking the values of any number of partial derivatives.
\code
#ifndef ADLOCSUPPORT_H

23
scripts/eigen_gen_credits Executable file
View File

@ -0,0 +1,23 @@
#!/bin/sh
# this script must be run from the eigen2/ directory.
# when running hg churn from the scripts/ subdir, i hit a divide-by-zero error.
#
# like this:
# cd eigen2
# USER=yourtuxfamilyuser scripts/eigen_gen_credits
rm -f eigen_gen_credits.log
hg pull >> eigen_gen_credits.log
wget http://eigen.tuxfamily.org/index.php?title=ContributorsInfo -O online-info.out -a eigen_gen_credits.log
hg churn -r 37: --changesets -t {author} > churn-changesets.out
hg churn -r 37: -t {author} > churn-changedlines.out
g++ scripts/eigen_gen_credits.cpp -o e
./e > credits.out
rsync credits.out $USER@ssh.tuxfamily.org:eigen/eigen.tuxfamily.org-web/htdocs/credits.out || (echo "upload failed"; exit 1)
ssh $USER@ssh.tuxfamily.org "cd eigen/eigen.tuxfamily.org-web/htdocs; chmod 664 credits.out; echo Main_Page | /usr/bin/php maintenance/purgeList.php"

213
scripts/eigen_gen_credits.cpp Normal file
View File

@ -0,0 +1,213 @@
#include <string>
#include <sstream>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <map>
#include <list>
using namespace std;
// this function takes a line that may contain a name and/or email address,
// and returns just the name, while fixing the "bad cases".
std::string contributor_name(const std::string& line)
{
string result;
size_t position_of_email_address = line.find_first_of('<');
if(position_of_email_address != string::npos)
{
// there is an e-mail address.
// Hauke once committed as "John Smith", fix that.
if(line.find("hauke.heibel") != string::npos)
result = "Hauke Heibel";
else
{
// just remove the e-mail address
result = line.substr(0, position_of_email_address);
}
}
else
{
// there is no e-mail address.
if(line.find("convert-repo") != string::npos)
result = "";
else
result = line;
}
// remove trailing spaces
size_t length = result.length();
while(length >= 1 && result[length-1] == ' ') result.erase(--length);
return result;
}
// parses hg churn output to generate a contributors map.
map<string,int> contributors_map_from_churn_output(const char *filename)
{
map<string,int> contributors_map;
string line;
ifstream churn_out;
churn_out.open(filename, ios::in);
while(!getline(churn_out,line).eof())
{
// remove the histograms "******" that hg churn may draw at the end of some lines
size_t first_star = line.find_first_of('*');
if(first_star != string::npos) line.erase(first_star);
// remove trailing spaces
size_t length = line.length();
while(length >= 1 && line[length-1] == ' ') line.erase(--length);
// now the last space indicates where the number starts
size_t last_space = line.find_last_of(' ');
// get the number (of changesets or of modified lines for each contributor)
int number;
istringstream(line.substr(last_space+1)) >> number;
// get the name of the contributor
line.erase(last_space);
string name = contributor_name(line);
map<string,int>::iterator it = contributors_map.find(name);
// if new contributor, insert
if(it == contributors_map.end())
contributors_map.insert(pair<string,int>(name, number));
// if duplicate, just add the number
else
it->second += number;
}
churn_out.close();
return contributors_map;
}
// find the last name, i.e. the last word.
// for "van den Schbling" types of last names, that's not a problem, that's actually what we want.
string lastname(const string& name)
{
size_t last_space = name.find_last_of(' ');
if(last_space >= name.length()-1) return name;
else return name.substr(last_space+1);
}
struct contributor
{
string name;
int changedlines;
int changesets;
string url;
string misc;
contributor() : changedlines(0), changesets(0) {}
bool operator < (const contributor& other)
{
return lastname(name).compare(lastname(other.name)) < 0;
}
};
void add_online_info_into_contributors_list(list<contributor>& contributors_list, const char *filename)
{
string line;
ifstream online_info;
online_info.open(filename, ios::in);
while(!getline(online_info,line).eof())
{
string hgname, realname, url, misc;
size_t last_bar = line.find_last_of('|');
if(last_bar == string::npos) continue;
if(last_bar < line.length())
misc = line.substr(last_bar+1);
line.erase(last_bar);
last_bar = line.find_last_of('|');
if(last_bar == string::npos) continue;
if(last_bar < line.length())
url = line.substr(last_bar+1);
line.erase(last_bar);
last_bar = line.find_last_of('|');
if(last_bar == string::npos) continue;
if(last_bar < line.length())
realname = line.substr(last_bar+1);
line.erase(last_bar);
hgname = line;
// remove the example line
if(hgname.find("MercurialName") != string::npos) continue;
list<contributor>::iterator it;
for(it=contributors_list.begin(); it != contributors_list.end() && it->name != hgname; ++it)
{}
if(it == contributors_list.end())
{
contributor c;
c.name = realname;
c.url = url;
c.misc = misc;
contributors_list.push_back(c);
}
else
{
it->name = realname;
it->url = url;
it->misc = misc;
}
}
}
int main()
{
// parse the hg churn output files
map<string,int> contributors_map_for_changedlines = contributors_map_from_churn_output("churn-changedlines.out");
//map<string,int> contributors_map_for_changesets = contributors_map_from_churn_output("churn-changesets.out");
// merge into the contributors list
list<contributor> contributors_list;
map<string,int>::iterator it;
for(it=contributors_map_for_changedlines.begin(); it != contributors_map_for_changedlines.end(); ++it)
{
contributor c;
c.name = it->first;
c.changedlines = it->second;
c.changesets = 0; //contributors_map_for_changesets.find(it->first)->second;
contributors_list.push_back(c);
}
add_online_info_into_contributors_list(contributors_list, "online-info.out");
contributors_list.sort();
cout << "{| cellpadding=\"5\"\n";
cout << "!\n";
cout << "! Lines changed\n";
cout << "!\n";
list<contributor>::iterator itc;
int i = 0;
for(itc=contributors_list.begin(); itc != contributors_list.end(); ++itc)
{
if(itc->name.length() == 0) continue;
if(i%2) cout << "|-\n";
else cout << "|- style=\"background:#FFFFD0\"\n";
if(itc->url.length())
cout << "| [" << itc->url << " " << itc->name << "]\n";
else
cout << "| " << itc->name << "\n";
if(itc->changedlines)
cout << "| " << itc->changedlines << "\n";
else
cout << "| (no information)\n";
cout << "| " << itc->misc << "\n";
i++;
}
cout << "|}" << endl;
}

11
scripts/eigen_gen_docs
View File

@ -1,18 +1,17 @@
#!/bin/sh
# configuration
USER='orzel'
# You should call this script with USER set as you want, else some default
# will be used
USER=${USER:-'orzel'}
# step 1 : update
hg pull -u || (echo "update failed"; exit 1)
# step 2 : build
# step 1 : build
# todo if 'build is not there, create one:
#mkdir build
(cd build && cmake .. && make -j3 doc) || (echo "make failed"; exit 1)
#todo: n+1 where n = number of cpus
#step 3 : upload
#step 2 : upload
BRANCH=`hg branch`
if [ $BRANCH == "default" ]
then

4
test/CMakeLists.txt
View File

@ -1,6 +1,6 @@
project(EigenTesting)
add_custom_target(btest)
include(EigenTesting)
enable_testing()
ei_init_testing()
find_package(GSL)

10
test/array_replicate.cpp
View File

@ -42,9 +42,9 @@ template<typename MatrixType> void replicate(const MatrixType& m)
MatrixType m1 = MatrixType::Random(rows, cols),
m2 = MatrixType::Random(rows, cols);
VectorType v1 = VectorType::Random(rows);
MatrixX x1, x2;
VectorX vx1;
@ -56,17 +56,17 @@ template<typename MatrixType> void replicate(const MatrixType& m)
for(int i=0; i<f1; i++)
x1.block(i*rows,j*cols,rows,cols) = m1;
VERIFY_IS_APPROX(x1, m1.replicate(f1,f2));
x2.resize(2*rows,3*cols);
x2 << m2, m2, m2,
m2, m2, m2;
VERIFY_IS_APPROX(x2, (m2.template replicate<2,3>()));
x2.resize(rows,f1);
for (int j=0; j<f1; ++j)
x2.col(j) = v1;
VERIFY_IS_APPROX(x2, v1.rowwise().replicate(f1));
vx1.resize(rows*f2);
for (int j=0; j<f2; ++j)
vx1.segment(j*rows,rows) = v1;

14
test/main.h
View File

@ -248,18 +248,22 @@ template<typename MatrixType>
void createRandomMatrixOfRank(int desired_rank, int rows, int cols, MatrixType& m)
{
typedef typename ei_traits<MatrixType>::Scalar Scalar;
typedef Matrix<Scalar, Dynamic, 1> VectorType;
enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime };
MatrixType a = MatrixType::Random(rows,rows);
typedef Matrix<Scalar, Dynamic, 1> VectorType;
typedef Matrix<Scalar, Rows, Rows> MatrixAType;
typedef Matrix<Scalar, Cols, Cols> MatrixBType;
MatrixAType a = MatrixAType::Random(rows,rows);
MatrixType d = MatrixType::Identity(rows,cols);
MatrixType b = MatrixType::Random(cols,cols);
MatrixBType b = MatrixBType::Random(cols,cols);
// set the diagonal such that only desired_rank non-zero entries reamain
const int diag_size = std::min(d.rows(),d.cols());
d.diagonal().segment(desired_rank, diag_size-desired_rank) = VectorType::Zero(diag_size-desired_rank);
HouseholderQR<MatrixType> qra(a);
HouseholderQR<MatrixType> qrb(b);
HouseholderQR<MatrixAType> qra(a);
HouseholderQR<MatrixBType> qrb(b);
m = qra.matrixQ() * d * qrb.matrixQ();
}

14
test/packetmath.cpp
View File

@ -99,10 +99,10 @@ template<typename Scalar> void packetmath()
const int PacketSize = ei_packet_traits<Scalar>::size;
const int size = PacketSize*4;
EIGEN_ALIGN_128 Scalar data1[ei_packet_traits<Scalar>::size*4];
EIGEN_ALIGN_128 Scalar data2[ei_packet_traits<Scalar>::size*4];
EIGEN_ALIGN_128 Packet packets[PacketSize*2];
EIGEN_ALIGN_128 Scalar ref[ei_packet_traits<Scalar>::size*4];
EIGEN_ALIGN16 Scalar data1[ei_packet_traits<Scalar>::size*4];
EIGEN_ALIGN16 Scalar data2[ei_packet_traits<Scalar>::size*4];
EIGEN_ALIGN16 Packet packets[PacketSize*2];
EIGEN_ALIGN16 Scalar ref[ei_packet_traits<Scalar>::size*4];
for (int i=0; i<size; ++i)
{
data1[i] = ei_random<Scalar>();
@ -202,9 +202,9 @@ template<typename Scalar> void packetmath_real()
const int PacketSize = ei_packet_traits<Scalar>::size;
const int size = PacketSize*4;
EIGEN_ALIGN_128 Scalar data1[ei_packet_traits<Scalar>::size*4];
EIGEN_ALIGN_128 Scalar data2[ei_packet_traits<Scalar>::size*4];
EIGEN_ALIGN_128 Scalar ref[ei_packet_traits<Scalar>::size*4];
EIGEN_ALIGN16 Scalar data1[ei_packet_traits<Scalar>::size*4];
EIGEN_ALIGN16 Scalar data2[ei_packet_traits<Scalar>::size*4];
EIGEN_ALIGN16 Scalar ref[ei_packet_traits<Scalar>::size*4];
for (int i=0; i<size; ++i)
{

1
test/product_extra.cpp
View File

@ -114,7 +114,6 @@ template<typename MatrixType> void product_extra(const MatrixType& m)
VERIFY_IS_APPROX(m1.col(j2).adjoint() * m1.block(0,j,m1.rows(),c), m1.col(j2).adjoint().eval() * m1.block(0,j,m1.rows(),c).eval());
VERIFY_IS_APPROX(m1.block(i,0,r,m1.cols()) * m1.row(i2).adjoint(), m1.block(i,0,r,m1.cols()).eval() * m1.row(i2).adjoint().eval());
}
void test_product_extra()

6
test/product_small.cpp
View File

@ -34,4 +34,10 @@ void test_product_small()
CALL_SUBTEST( product(Matrix4d()) );
CALL_SUBTEST( product(Matrix4f()) );
}
{
// test compilation of (outer_product) * vector
Vector3f v = Vector3f::Random();
VERIFY_IS_APPROX( (v * v.transpose()) * v, (v * v.transpose()).eval() * v);
}
}

44
test/qr.cpp
View File

@ -31,30 +31,40 @@ template<typename MatrixType> void qr(const MatrixType& m)
int cols = m.cols();
typedef typename MatrixType::Scalar Scalar;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> SquareMatrixType;
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
MatrixType a = MatrixType::Random(rows,cols);
HouseholderQR<MatrixType> qrOfA(a);
MatrixType r = qrOfA.matrixQR();
MatrixQType q = qrOfA.matrixQ();
VERIFY_IS_UNITARY(q);
// FIXME need better way to construct trapezoid
for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0);
VERIFY_IS_APPROX(a, qrOfA.matrixQ() * r);
}
SquareMatrixType b = a.adjoint() * a;
template<typename MatrixType, int Cols2> void qr_fixedsize()
{
enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime };
typedef typename MatrixType::Scalar Scalar;
Matrix<Scalar,Rows,Cols> m1 = Matrix<Scalar,Rows,Cols>::Random();
HouseholderQR<Matrix<Scalar,Rows,Cols> > qr(m1);
// check tridiagonalization
Tridiagonalization<SquareMatrixType> tridiag(b);
VERIFY_IS_APPROX(b, tridiag.matrixQ() * tridiag.matrixT() * tridiag.matrixQ().adjoint());
Matrix<Scalar,Rows,Cols> r = qr.matrixQR();
// FIXME need better way to construct trapezoid
for(int i = 0; i < Rows; i++) for(int j = 0; j < Cols; j++) if(i>j) r(i,j) = Scalar(0);
// check hessenberg decomposition
HessenbergDecomposition<SquareMatrixType> hess(b);
VERIFY_IS_APPROX(b, hess.matrixQ() * hess.matrixH() * hess.matrixQ().adjoint());
VERIFY_IS_APPROX(tridiag.matrixT(), hess.matrixH());
b = SquareMatrixType::Random(cols,cols);
hess.compute(b);
VERIFY_IS_APPROX(b, hess.matrixQ() * hess.matrixH() * hess.matrixQ().adjoint());
VERIFY_IS_APPROX(m1, qr.matrixQ() * r);
Matrix<Scalar,Cols,Cols2> m2 = Matrix<Scalar,Cols,Cols2>::Random(Cols,Cols2);
Matrix<Scalar,Rows,Cols2> m3 = m1*m2;
m2 = Matrix<Scalar,Cols,Cols2>::Random(Cols,Cols2);
qr.solve(m3, &m2);
VERIFY_IS_APPROX(m3, m1*m2);
}
template<typename MatrixType> void qr_invertible()
@ -105,11 +115,11 @@ template<typename MatrixType> void qr_verify_assert()
void test_qr()
{
for(int i = 0; i < 1; i++) {
// FIXME : very weird bug here
// CALL_SUBTEST( qr(Matrix2f()) );
CALL_SUBTEST( qr(Matrix4d()) );
CALL_SUBTEST( qr(MatrixXf(47,40)) );
CALL_SUBTEST( qr(MatrixXcd(17,7)) );
CALL_SUBTEST( qr(MatrixXf(47,40)) );
CALL_SUBTEST( qr(MatrixXcd(17,7)) );
CALL_SUBTEST(( qr_fixedsize<Matrix<float,3,4>, 2 >() ));
CALL_SUBTEST(( qr_fixedsize<Matrix<double,6,2>, 4 >() ));
CALL_SUBTEST(( qr_fixedsize<Matrix<double,2,5>, 7 >() ));
}
for(int i = 0; i < g_repeat; i++) {

42
test/qr_colpivoting.cpp
View File

@ -32,7 +32,7 @@ template<typename MatrixType> void qr()
int rank = ei_random<int>(1, std::min(rows, cols)-1);
typedef typename MatrixType::Scalar Scalar;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> SquareMatrixType;
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
MatrixType m1;
createRandomMatrixOfRank(rank,rows,cols,m1);
@ -44,6 +44,10 @@ template<typename MatrixType> void qr()
VERIFY(!qr.isSurjective());
MatrixType r = qr.matrixQR();
MatrixQType q = qr.matrixQ();
VERIFY_IS_UNITARY(q);
// FIXME need better way to construct trapezoid
for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0);
@ -63,6 +67,40 @@ template<typename MatrixType> void qr()
VERIFY(!qr.solve(m3, &m2));
}
template<typename MatrixType, int Cols2> void qr_fixedsize()
{
enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime };
typedef typename MatrixType::Scalar Scalar;
int rank = ei_random<int>(1, std::min(int(Rows), int(Cols))-1);
Matrix<Scalar,Rows,Cols> m1;
createRandomMatrixOfRank(rank,Rows,Cols,m1);
ColPivotingHouseholderQR<Matrix<Scalar,Rows,Cols> > qr(m1);
VERIFY_IS_APPROX(rank, qr.rank());
VERIFY(Cols - qr.rank() == qr.dimensionOfKernel());
VERIFY(!qr.isInjective());
VERIFY(!qr.isInvertible());
VERIFY(!qr.isSurjective());
Matrix<Scalar,Rows,Cols> r = qr.matrixQR();
// FIXME need better way to construct trapezoid
for(int i = 0; i < Rows; i++) for(int j = 0; j < Cols; j++) if(i>j) r(i,j) = Scalar(0);
Matrix<Scalar,Rows,Cols> b = qr.matrixQ() * r;
Matrix<Scalar,Rows,Cols> c = MatrixType::Zero(Rows,Cols);
for(int i = 0; i < Cols; ++i) c.col(qr.colsPermutation().coeff(i)) = b.col(i);
VERIFY_IS_APPROX(m1, c);
Matrix<Scalar,Cols,Cols2> m2 = Matrix<Scalar,Cols,Cols2>::Random(Cols,Cols2);
Matrix<Scalar,Rows,Cols2> m3 = m1*m2;
m2 = Matrix<Scalar,Cols,Cols2>::Random(Cols,Cols2);
VERIFY(qr.solve(m3, &m2));
VERIFY_IS_APPROX(m3, m1*m2);
m3 = Matrix<Scalar,Rows,Cols2>::Random(Rows,Cols2);
VERIFY(!qr.solve(m3, &m2));
}
template<typename MatrixType> void qr_invertible()
{
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
@ -120,6 +158,8 @@ void test_qr_colpivoting()
CALL_SUBTEST( qr<MatrixXf>() );
CALL_SUBTEST( qr<MatrixXd>() );
CALL_SUBTEST( qr<MatrixXcd>() );
CALL_SUBTEST(( qr_fixedsize<Matrix<float,3,5>, 4 >() ));
CALL_SUBTEST(( qr_fixedsize<Matrix<double,6,2>, 3 >() ));
}
for(int i = 0; i < g_repeat; i++) {

6
test/qr_fullpivoting.cpp
View File

@ -32,7 +32,7 @@ template<typename MatrixType> void qr()
int rank = ei_random<int>(1, std::min(rows, cols)-1);
typedef typename MatrixType::Scalar Scalar;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> SquareMatrixType;
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
MatrixType m1;
createRandomMatrixOfRank(rank,rows,cols,m1);
@ -44,6 +44,10 @@ template<typename MatrixType> void qr()
VERIFY(!qr.isSurjective());
MatrixType r = qr.matrixQR();
MatrixQType q = qr.matrixQ();
VERIFY_IS_UNITARY(q);
// FIXME need better way to construct trapezoid
for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0);

14
test/stable_norm.cpp
View File

@ -33,6 +33,20 @@ template<typename MatrixType> void stable_norm(const MatrixType& m)
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
// Check the basic machine-dependent constants.
{
int ibeta, it, iemin, iemax;
ibeta = std::numeric_limits<RealScalar>::radix; // base for floating-point numbers
it = std::numeric_limits<RealScalar>::digits; // number of base-beta digits in mantissa
iemin = std::numeric_limits<RealScalar>::min_exponent; // minimum exponent
iemax = std::numeric_limits<RealScalar>::max_exponent; // maximum exponent
VERIFY( (!(iemin > 1 - 2*it || 1+it>iemax || (it==2 && ibeta<5) || (it<=4 && ibeta <= 3 ) || it<2))
&& "the stable norm algorithm cannot be guaranteed on this computer");
}
int rows = m.rows();
int cols = m.cols();

6
test/testsuite.cmake
View File

@ -148,7 +148,7 @@ endif(NOT EIGEN_NO_UPDATE)
# which ctest command to use for running the dashboard
SET (CTEST_COMMAND "${EIGEN_CMAKE_DIR}ctest -D ${EIGEN_MODE}")
# what cmake command to use for configuring this dashboard
SET (CTEST_CMAKE_COMMAND "${EIGEN_CMAKE_DIR}cmake -DEIGEN_BUILD_TESTS=on ")
SET (CTEST_CMAKE_COMMAND "${EIGEN_CMAKE_DIR}cmake -DEIGEN_CMAKE_RUN_FROM_CTEST=ON")
####################################################################
# The values in this section are optional you can either
@ -175,9 +175,9 @@ if(WIN32 AND NOT UNIX)
else(EIGEN_GENERATOR_TYPE)
set(CTEST_CMAKE_COMMAND "${CTEST_CMAKE_COMMAND} -G \"NMake Makefiles\" -DCMAKE_MAKE_PROGRAM=nmake")
SET (CTEST_INITIAL_CACHE "
MAKECOMMAND:STRING=nmake -i
MAKECOMMAND:STRING=nmake /i
CMAKE_MAKE_PROGRAM:FILEPATH=nmake
CMAKE_GENERATOR:INTERNAL=NMake Makefiles
CMAKE_GENERATOR:INTERNAL=NMake Makefiles
CMAKE_BUILD_TYPE:STRING=Release
BUILDNAME:STRING=${EIGEN_BUILD_STRING}
SITE:STRING=${EIGEN_SITE}

94
test/unalignedassert.cpp
View File

@ -24,52 +24,38 @@
#include "main.h"
struct Good1
struct TestNew1
{
MatrixXd m; // good: m will allocate its own array, taking care of alignment.
Good1() : m(20,20) {}
TestNew1() : m(20,20) {}
};
struct Good2
struct TestNew2
{
Matrix3d m; // good: m's size isn't a multiple of 16 bytes, so m doesn't have to be aligned
Matrix3d m; // good: m's size isn't a multiple of 16 bytes, so m doesn't have to be 16-byte aligned,
// 8-byte alignment is good enough here, which we'll get automatically
};
struct Good3
struct TestNew3
{
Vector2f m; // good: same reason
Vector2f m; // good: m's size isn't a multiple of 16 bytes, so m doesn't have to be 16-byte aligned
};
struct Bad4
{
Vector2d m; // bad: sizeof(m)%16==0 so alignment is required
};
struct Bad5
{
Matrix<float, 2, 6> m; // bad: same reason
};
struct Bad6
{
Matrix<double, 3, 4> m; // bad: same reason
};
struct Good7
struct TestNew4
{
EIGEN_MAKE_ALIGNED_OPERATOR_NEW
Vector2d m;
float f; // make the struct have sizeof%16!=0 to make it a little more tricky when we allow an array of 2 such objects
};
struct Good8
struct TestNew5
{
EIGEN_MAKE_ALIGNED_OPERATOR_NEW
float f; // try the f at first -- the EIGEN_ALIGN_128 attribute of m should make that still work
float f; // try the f at first -- the EIGEN_ALIGN16 attribute of m should make that still work
Matrix4f m;
};
struct Good9
struct TestNew6
{
Matrix<float,2,2,DontAlign> m; // good: no alignment requested
float f;
@ -94,34 +80,58 @@ void check_unalignedassert_good()
#if EIGEN_ALIGN
template<typename T>
void check_unalignedassert_bad()
void construct_at_boundary(int boundary)
{
float buf[sizeof(T)+16];
float *unaligned = buf;
while((reinterpret_cast<size_t>(unaligned)&0xf)==0) ++unaligned; // make sure unaligned is really unaligned
T *x = ::new(static_cast<void*>(unaligned)) T;
char buf[sizeof(T)+256];
size_t _buf = reinterpret_cast<size_t>(buf);
_buf += (16 - (_buf % 16)); // make 16-byte aligned
_buf += boundary; // make exact boundary-aligned
T *x = ::new(reinterpret_cast<void*>(_buf)) T;
x->~T();
}
#endif
void unalignedassert()
{
check_unalignedassert_good<Good1>();
check_unalignedassert_good<Good2>();
check_unalignedassert_good<Good3>();
#if EIGEN_ALIGN
VERIFY_RAISES_ASSERT(check_unalignedassert_bad<Bad4>());
VERIFY_RAISES_ASSERT(check_unalignedassert_bad<Bad5>());
VERIFY_RAISES_ASSERT(check_unalignedassert_bad<Bad6>());
#endif
#if EIGEN_ALIGN
construct_at_boundary<Vector2f>(4);
construct_at_boundary<Vector3f>(4);
construct_at_boundary<Vector4f>(16);
construct_at_boundary<Matrix2f>(16);
construct_at_boundary<Matrix3f>(4);
construct_at_boundary<Matrix4f>(16);
check_unalignedassert_good<Good7>();
check_unalignedassert_good<Good8>();
check_unalignedassert_good<Good9>();
construct_at_boundary<Vector2d>(16);
construct_at_boundary<Vector3d>(4);
construct_at_boundary<Vector4d>(16);
construct_at_boundary<Matrix2d>(16);
construct_at_boundary<Matrix3d>(4);
construct_at_boundary<Matrix4d>(16);
construct_at_boundary<Vector2cf>(16);
construct_at_boundary<Vector3cf>(4);
construct_at_boundary<Vector2cd>(16);
construct_at_boundary<Vector3cd>(16);
#endif
check_unalignedassert_good<TestNew1>();
check_unalignedassert_good<TestNew2>();
check_unalignedassert_good<TestNew3>();
check_unalignedassert_good<TestNew4>();
check_unalignedassert_good<TestNew5>();
check_unalignedassert_good<TestNew6>();
check_unalignedassert_good<Depends<true> >();
#if EIGEN_ALIGN
VERIFY_RAISES_ASSERT(check_unalignedassert_bad<Depends<false> >());
VERIFY_RAISES_ASSERT(construct_at_boundary<Vector4f>(8));
VERIFY_RAISES_ASSERT(construct_at_boundary<Matrix4f>(8));
VERIFY_RAISES_ASSERT(construct_at_boundary<Vector2d>(8));
VERIFY_RAISES_ASSERT(construct_at_boundary<Vector4d>(8));
VERIFY_RAISES_ASSERT(construct_at_boundary<Matrix2d>(8));
VERIFY_RAISES_ASSERT(construct_at_boundary<Matrix4d>(8));
VERIFY_RAISES_ASSERT(construct_at_boundary<Vector2cf>(8));
VERIFY_RAISES_ASSERT(construct_at_boundary<Vector2cd>(8));
#endif
}

4
test/visitor.cpp
View File

@ -40,7 +40,7 @@ template<typename MatrixType> void matrixVisitor(const MatrixType& p)
m(i) = ei_random<Scalar>();
Scalar minc = Scalar(1000), maxc = Scalar(-1000);
int minrow,mincol,maxrow,maxcol;
int minrow=0,mincol=0,maxrow=0,maxcol=0;
for(int j = 0; j < cols; j++)
for(int i = 0; i < rows; i++)
{
@ -86,7 +86,7 @@ template<typename VectorType> void vectorVisitor(const VectorType& w)
v(i) = ei_random<Scalar>();
Scalar minc = Scalar(1000), maxc = Scalar(-1000);
int minidx,maxidx;
int minidx=0,maxidx=0;
for(int i = 0; i < size; i++)
{
if(v(i) < minc)

10
unsupported/CMakeLists.txt
View File

@ -1,9 +1,3 @@
add_subdirectory(Eigen)
add_subdirectory(doc)
if(EIGEN_BUILD_TESTS)
add_subdirectory(test)
endif(EIGEN_BUILD_TESTS)
add_subdirectory(doc EXCLUDE_FROM_ALL)
add_subdirectory(test EXCLUDE_FROM_ALL)

4
unsupported/Eigen/AdolcForward
View File

@ -29,7 +29,7 @@
//
// This file provides support for adolc's adouble type in forward mode.
// ADOL-C is a C++ automatic differentiation library,
// see http://www.math.tu-dresden.de/~adol-c/ for more information.
// see https://projects.coin-or.org/ADOL-C for more information.
//
// Note that the maximal number of directions is controlled by
// the preprocessor token NUMBER_DIRECTIONS. The default is 2.
@ -63,7 +63,7 @@ namespace Eigen {
* \defgroup AdolcForward_Module Adolc forward module
* This module provides support for adolc's adouble type in forward mode.
* ADOL-C is a C++ automatic differentiation library,
* see http://www.math.tu-dresden.de/~adol-c/ for more information.
* see https://projects.coin-or.org/ADOL-C for more information.
* It mainly consists in:
* - a struct Eigen::NumTraits<adtl::adouble> specialization
* - overloads of ei_* math function for adtl::adouble type.

25
unsupported/Eigen/src/AutoDiff/AutoDiffJacobian.h
View File

@ -46,17 +46,18 @@ public:
InputsAtCompileTime = Functor::InputsAtCompileTime,
ValuesAtCompileTime = Functor::ValuesAtCompileTime
};
typedef typename Functor::InputType InputType;
typedef typename Functor::ValueType ValueType;
typedef typename Functor::JacobianType JacobianType;
typedef AutoDiffScalar<Matrix<double,InputsAtCompileTime,1> > ActiveScalar;
typedef Matrix<double,InputsAtCompileTime,1> DerivativeType;
typedef AutoDiffScalar<DerivativeType> ActiveScalar;
typedef Matrix<ActiveScalar, InputsAtCompileTime, 1> ActiveInput;
typedef Matrix<ActiveScalar, ValuesAtCompileTime, 1> ActiveValue;
void operator() (const InputType& x, ValueType* v, JacobianType* _jac) const
void operator() (const InputType& x, ValueType* v, JacobianType* _jac=0) const
{
ei_assert(v!=0);
if (!_jac)
@ -69,26 +70,20 @@ public:
ActiveInput ax = x.template cast<ActiveScalar>();
ActiveValue av(jac.rows());
if(InputsAtCompileTime==Dynamic)
{
for (int j=0; j<jac.cols(); j++)
ax[j].derivatives().resize(this->inputs());
for (int j=0; j<jac.rows(); j++)
av[j].derivatives().resize(this->inputs());
}
for (int j=0; j<jac.cols(); j++)
for (int i=0; i<jac.cols(); i++)
ax[i].derivatives().coeffRef(j) = i==j ? 1 : 0;
for (int i=0; i<jac.cols(); i++)
ax[i].derivatives() = DerivativeType::Unit(this->inputs(),i);
Functor::operator()(ax, &av);
for (int i=0; i<jac.rows(); i++)
{
(*v)[i] = av[i].value();
for (int j=0; j<jac.cols(); j++)
jac.coeffRef(i,j) = av[i].derivatives().coeff(j);
jac.row(i) = av[i].derivatives();
}
}
protected:

132
unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h
View File

@ -27,6 +27,18 @@
namespace Eigen {
template<typename A, typename B>
struct ei_make_coherent_impl {
static void run(A& a, B& b) {}
};
// resize a to match b is a.size()==0, and conversely.
template<typename A, typename B>
void ei_make_coherent(const A& a, const B&b)
{
ei_make_coherent_impl<A,B>::run(a.const_cast_derived(), b.const_cast_derived());
}
/** \class AutoDiffScalar
* \brief A scalar type replacement with automatic differentation capability
*
@ -35,7 +47,7 @@ namespace Eigen {
* This class represents a scalar value while tracking its respective derivatives.
*
* It supports the following list of global math function:
* - std::abs, std::sqrt, std::pow, std::exp, std::log, std::sin, std::cos,
* - std::abs, std::sqrt, std::pow, std::exp, std::log, std::sin, std::cos,
* - ei_abs, ei_sqrt, ei_pow, ei_exp, ei_log, ei_sin, ei_cos,
* - ei_conj, ei_real, ei_imag, ei_abs2.
*
@ -49,29 +61,29 @@ class AutoDiffScalar
{
public:
typedef typename ei_traits<DerType>::Scalar Scalar;
inline AutoDiffScalar() {}
inline AutoDiffScalar(const Scalar& value)
: m_value(value)
{
if(m_derivatives.size()>0)
m_derivatives.setZero();
}
inline AutoDiffScalar(const Scalar& value, const DerType& der)
: m_value(value), m_derivatives(der)
{}
template<typename OtherDerType>
inline AutoDiffScalar(const AutoDiffScalar<OtherDerType>& other)
: m_value(other.value()), m_derivatives(other.derivatives())
{}
inline AutoDiffScalar(const AutoDiffScalar& other)
: m_value(other.value()), m_derivatives(other.derivatives())
{}
template<typename OtherDerType>
inline AutoDiffScalar& operator=(const AutoDiffScalar<OtherDerType>& other)
{
@ -79,32 +91,33 @@ class AutoDiffScalar
m_derivatives = other.derivatives();
return *this;
}
inline AutoDiffScalar& operator=(const AutoDiffScalar& other)
{
m_value = other.value();
m_derivatives = other.derivatives();
return *this;
}
// inline operator const Scalar& () const { return m_value; }
// inline operator Scalar& () { return m_value; }
inline const Scalar& value() const { return m_value; }
inline Scalar& value() { return m_value; }
inline const DerType& derivatives() const { return m_derivatives; }
inline DerType& derivatives() { return m_derivatives; }
template<typename OtherDerType>
inline const AutoDiffScalar<CwiseBinaryOp<ei_scalar_sum_op<Scalar>,DerType,OtherDerType> >
operator+(const AutoDiffScalar<OtherDerType>& other) const
{
ei_make_coherent(m_derivatives, other.derivatives());
return AutoDiffScalar<CwiseBinaryOp<ei_scalar_sum_op<Scalar>,DerType,OtherDerType> >(
m_value + other.value(),
m_derivatives + other.derivatives());
}
template<typename OtherDerType>
inline AutoDiffScalar&
operator+=(const AutoDiffScalar<OtherDerType>& other)
@ -112,16 +125,17 @@ class AutoDiffScalar
(*this) = (*this) + other;
return *this;
}
template<typename OtherDerType>
inline const AutoDiffScalar<CwiseBinaryOp<ei_scalar_difference_op<Scalar>, DerType,OtherDerType> >
operator-(const AutoDiffScalar<OtherDerType>& other) const
{
ei_make_coherent(m_derivatives, other.derivatives());
return AutoDiffScalar<CwiseBinaryOp<ei_scalar_difference_op<Scalar>, DerType,OtherDerType> >(
m_value - other.value(),
m_derivatives - other.derivatives());
}
template<typename OtherDerType>
inline AutoDiffScalar&
operator-=(const AutoDiffScalar<OtherDerType>& other)
@ -129,7 +143,7 @@ class AutoDiffScalar
*this = *this - other;
return *this;
}
template<typename OtherDerType>
inline const AutoDiffScalar<CwiseUnaryOp<ei_scalar_opposite_op<Scalar>, DerType> >
operator-() const
@ -138,7 +152,7 @@ class AutoDiffScalar
-m_value,
-m_derivatives);
}
inline const AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, DerType> >
operator*(const Scalar& other) const
{
@ -146,7 +160,7 @@ class AutoDiffScalar
m_value * other,
(m_derivatives * other));
}
friend inline const AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, DerType> >
operator*(const Scalar& other, const AutoDiffScalar& a)
{
@ -154,7 +168,7 @@ class AutoDiffScalar
a.value() * other,
a.derivatives() * other);
}
inline const AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, DerType> >
operator/(const Scalar& other) const
{
@ -162,7 +176,7 @@ class AutoDiffScalar
m_value / other,
(m_derivatives * (Scalar(1)/other)));
}
friend inline const AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, DerType> >
operator/(const Scalar& other, const AutoDiffScalar& a)
{
@ -170,7 +184,7 @@ class AutoDiffScalar
other / a.value(),
a.derivatives() * (-Scalar(1)/other));
}
template<typename OtherDerType>
inline const AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>,
NestByValue<CwiseBinaryOp<ei_scalar_difference_op<Scalar>,
@ -178,6 +192,7 @@ class AutoDiffScalar
NestByValue<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, OtherDerType> > > > > >
operator/(const AutoDiffScalar<OtherDerType>& other) const
{
ei_make_coherent(m_derivatives, other.derivatives());
return AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>,
NestByValue<CwiseBinaryOp<ei_scalar_difference_op<Scalar>,
NestByValue<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, DerType> >,
@ -186,45 +201,91 @@ class AutoDiffScalar
((m_derivatives * other.value()).nestByValue() - (m_value * other.derivatives()).nestByValue()).nestByValue()
* (Scalar(1)/(other.value()*other.value())));
}
template<typename OtherDerType>
inline const AutoDiffScalar<CwiseBinaryOp<ei_scalar_sum_op<Scalar>,
NestByValue<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, DerType> >,
NestByValue<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, OtherDerType> > > >
operator*(const AutoDiffScalar<OtherDerType>& other) const
{
ei_make_coherent(m_derivatives, other.derivatives());
return AutoDiffScalar<CwiseBinaryOp<ei_scalar_sum_op<Scalar>,
NestByValue<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, DerType> >,
NestByValue<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, OtherDerType> > > >(
m_value * other.value(),
(m_derivatives * other.value()).nestByValue() + (m_value * other.derivatives()).nestByValue());
}
inline AutoDiffScalar& operator*=(const Scalar& other)
{
*this = *this * other;
return *this;
}
template<typename OtherDerType>
inline AutoDiffScalar& operator*=(const AutoDiffScalar<OtherDerType>& other)
{
*this = *this * other;
return *this;
}
protected:
Scalar m_value;
DerType m_derivatives;
};
template<typename A_Scalar, int A_Rows, int A_Cols, int A_Options, int A_MaxRows, int A_MaxCols, typename B>
struct ei_make_coherent_impl<Matrix<A_Scalar, A_Rows, A_Cols, A_Options, A_MaxRows, A_MaxCols>, B> {
typedef Matrix<A_Scalar, A_Rows, A_Cols, A_Options, A_MaxRows, A_MaxCols> A;
static void run(A& a, B& b) {
if((A_Rows==Dynamic || A_Cols==Dynamic) && (a.size()==0))
{
a.resize(b.size());
a.setZero();
}
}
};
template<typename A, typename B_Scalar, int B_Rows, int B_Cols, int B_Options, int B_MaxRows, int B_MaxCols>
struct ei_make_coherent_impl<A, Matrix<B_Scalar, B_Rows, B_Cols, B_Options, B_MaxRows, B_MaxCols> > {
typedef Matrix<B_Scalar, B_Rows, B_Cols, B_Options, B_MaxRows, B_MaxCols> B;
static void run(A& a, B& b) {
if((B_Rows==Dynamic || B_Cols==Dynamic) && (b.size()==0))
{
b.resize(a.size());
b.setZero();
}
}
};
template<typename A_Scalar, int A_Rows, int A_Cols, int A_Options, int A_MaxRows, int A_MaxCols,
typename B_Scalar, int B_Rows, int B_Cols, int B_Options, int B_MaxRows, int B_MaxCols>
struct ei_make_coherent_impl<Matrix<A_Scalar, A_Rows, A_Cols, A_Options, A_MaxRows, A_MaxCols>,
Matrix<B_Scalar, B_Rows, B_Cols, B_Options, B_MaxRows, B_MaxCols> > {
typedef Matrix<A_Scalar, A_Rows, A_Cols, A_Options, A_MaxRows, A_MaxCols> A;
typedef Matrix<B_Scalar, B_Rows, B_Cols, B_Options, B_MaxRows, B_MaxCols> B;
static void run(A& a, B& b) {
if((A_Rows==Dynamic || A_Cols==Dynamic) && (a.size()==0))
{
a.resize(b.size());
a.setZero();
}
else if((B_Rows==Dynamic || B_Cols==Dynamic) && (b.size()==0))
{
b.resize(a.size());
b.setZero();
}
}
};
}
#define EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(FUNC,CODE) \
template<typename DerType> \
inline const AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<typename ei_traits<DerType>::Scalar>, DerType> > \
FUNC(const AutoDiffScalar<DerType>& x) { \
inline const Eigen::AutoDiffScalar<Eigen::CwiseUnaryOp<Eigen::ei_scalar_multiple_op<typename Eigen::ei_traits<DerType>::Scalar>, DerType> > \
FUNC(const Eigen::AutoDiffScalar<DerType>& x) { \
using namespace Eigen; \
typedef typename ei_traits<DerType>::Scalar Scalar; \
typedef AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, DerType> > ReturnType; \
CODE; \
@ -234,34 +295,35 @@ namespace std
{
EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(abs,
return ReturnType(std::abs(x.value()), x.derivatives() * (sign(x.value())));)
EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(sqrt,
Scalar sqrtx = std::sqrt(x.value());
return ReturnType(sqrtx,x.derivatives() * (Scalar(0.5) / sqrtx));)
EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(cos,
return ReturnType(std::cos(x.value()), x.derivatives() * (-std::sin(x.value())));)
EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(sin,
return ReturnType(std::sin(x.value()),x.derivatives() * std::cos(x.value()));)
EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(exp,
Scalar expx = std::exp(x.value());
return ReturnType(expx,x.derivatives() * expx);)
EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(ei_log,
return ReturnType(std::log(x.value),x.derivatives() * (Scalar(1).x.value()));)
template<typename DerType>
inline const AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<typename ei_traits<DerType>::Scalar>, DerType> >
pow(const AutoDiffScalar<DerType>& x, typename ei_traits<DerType>::Scalar y)
inline const Eigen::AutoDiffScalar<Eigen::CwiseUnaryOp<Eigen::ei_scalar_multiple_op<typename Eigen::ei_traits<DerType>::Scalar>, DerType> >
pow(const Eigen::AutoDiffScalar<DerType>& x, typename Eigen::ei_traits<DerType>::Scalar y)
{
using namespace Eigen;
typedef typename ei_traits<DerType>::Scalar Scalar;
return AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, DerType> >(
std::pow(x.value(),y),
x.derivatives() * (y * std::pow(x.value(),y-1)));
}
}
namespace Eigen {

38
unsupported/test/autodiff.cpp
View File

@ -46,12 +46,12 @@ struct TestFunc1
typedef Matrix<Scalar,InputsAtCompileTime,1> InputType;
typedef Matrix<Scalar,ValuesAtCompileTime,1> ValueType;
typedef Matrix<Scalar,ValuesAtCompileTime,InputsAtCompileTime> JacobianType;
int m_inputs, m_values;
TestFunc1() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {}
TestFunc1(int inputs, int values) : m_inputs(inputs), m_values(values) {}
int inputs() const { return m_inputs; }
int values() const { return m_values; }
@ -111,7 +111,7 @@ struct TestFunc1
}
};
template<typename Func> void adolc_forward_jacobian(const Func& f)
template<typename Func> void forward_jacobian(const Func& f)
{
typename Func::InputType x = Func::InputType::Random(f.inputs());
typename Func::ValueType y(f.values()), yref(f.values());
@ -134,21 +134,29 @@ template<typename Func> void adolc_forward_jacobian(const Func& f)
VERIFY_IS_APPROX(j, jref);
}
void test_autodiff()
void test_autodiff_scalar()
{
std::cerr << foo<float>(1,2) << "\n";
AutoDiffScalar<Vector2f> ax(1,Vector2f::UnitX());
AutoDiffScalar<Vector2f> ay(2,Vector2f::UnitY());
std::cerr << foo<AutoDiffScalar<Vector2f> >(ax,ay).value() << " <> "
<< foo<AutoDiffScalar<Vector2f> >(ax,ay).derivatives().transpose() << "\n\n";
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,2,2>()) ));
CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,2,3>()) ));
CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,3,2>()) ));
CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,3,3>()) ));
CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double>(3,3)) ));
}
// exit(1);
}
void test_autodiff_jacobian()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST(( forward_jacobian(TestFunc1<double,2,2>()) ));
CALL_SUBTEST(( forward_jacobian(TestFunc1<double,2,3>()) ));
CALL_SUBTEST(( forward_jacobian(TestFunc1<double,3,2>()) ));
CALL_SUBTEST(( forward_jacobian(TestFunc1<double,3,3>()) ));
CALL_SUBTEST(( forward_jacobian(TestFunc1<double>(3,3)) ));
}
}
void test_autodiff()
{
test_autodiff_scalar();
test_autodiff_jacobian();
}