Skip to content
Hide whitespace changes
Inline Side-by-side
external/eigen3/Eigen/src/Core/Inverse.h 0 → 100644
View file @ a394b22a
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_INVERSE_H
#define EIGEN_INVERSE_H
namespace Eigen {
template<typename XprType,typename StorageKind> class InverseImpl;
namespace internal {
template<typename XprType>
struct traits<Inverse<XprType> >
: traits<typename XprType::PlainObject>
{
typedef typename XprType::PlainObject PlainObject;
typedef traits<PlainObject> BaseTraits;
enum {
Flags = BaseTraits::Flags & RowMajorBit
};
};
} // end namespace internal
/** \class Inverse
*
* \brief Expression of the inverse of another expression
*
* \tparam XprType the type of the expression we are taking the inverse
*
* This class represents an abstract expression of A.inverse()
* and most of the time this is the only way it is used.
*
*/
template<typename XprType>
class Inverse : public InverseImpl<XprType,typename internal::traits<XprType>::StorageKind>
{
public:
typedef typename XprType::StorageIndex StorageIndex;
typedef typename XprType::PlainObject PlainObject;
typedef typename XprType::Scalar Scalar;
typedef typename internal::ref_selector<XprType>::type XprTypeNested;
typedef typename internal::remove_all<XprTypeNested>::type XprTypeNestedCleaned;
typedef typename internal::ref_selector<Inverse>::type Nested;
typedef typename internal::remove_all<XprType>::type NestedExpression;
explicit EIGEN_DEVICE_FUNC Inverse(const XprType &xpr)
: m_xpr(xpr)
{}
EIGEN_DEVICE_FUNC Index rows() const { return m_xpr.rows(); }
EIGEN_DEVICE_FUNC Index cols() const { return m_xpr.cols(); }
EIGEN_DEVICE_FUNC const XprTypeNestedCleaned& nestedExpression() const { return m_xpr; }
protected:
XprTypeNested m_xpr;
};
// Generic API dispatcher
template<typename XprType, typename StorageKind>
class InverseImpl
: public internal::generic_xpr_base<Inverse<XprType> >::type
{
public:
typedef typename internal::generic_xpr_base<Inverse<XprType> >::type Base;
typedef typename XprType::Scalar Scalar;
private:
Scalar coeff(Index row, Index col) const;
Scalar coeff(Index i) const;
};
namespace internal {
/** \internal
* \brief Default evaluator for Inverse expression.
*
* This default evaluator for Inverse expression simply evaluate the inverse into a temporary
* by a call to internal::call_assignment_no_alias.
* Therefore, inverse implementers only have to specialize Assignment<Dst,Inverse<...>, ...> for
* there own nested expression.
*
* \sa class Inverse
*/
template<typename ArgType>
struct unary_evaluator<Inverse<ArgType> >
: public evaluator<typename Inverse<ArgType>::PlainObject>
{
typedef Inverse<ArgType> InverseType;
typedef typename InverseType::PlainObject PlainObject;
typedef evaluator<PlainObject> Base;
enum { Flags = Base::Flags | EvalBeforeNestingBit };
unary_evaluator(const InverseType& inv_xpr)
: m_result(inv_xpr.rows(), inv_xpr.cols())
{
::new (static_cast<Base*>(this)) Base(m_result);
internal::call_assignment_no_alias(m_result, inv_xpr);
}
protected:
PlainObject m_result;
};
} // end namespace internal
} // end namespace Eigen
#endif // EIGEN_INVERSE_H
external/eigen3/Eigen/src/Core/Map.h 0 → 100644
View file @ a394b22a
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2007-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_MAP_H
#define EIGEN_MAP_H
namespace Eigen {
namespace internal {
template<typename PlainObjectType, int MapOptions, typename StrideType>
struct traits<Map<PlainObjectType, MapOptions, StrideType> >
: public traits<PlainObjectType>
{
typedef traits<PlainObjectType> TraitsBase;
enum {
InnerStrideAtCompileTime = StrideType::InnerStrideAtCompileTime == 0
? int(PlainObjectType::InnerStrideAtCompileTime)
: int(StrideType::InnerStrideAtCompileTime),
OuterStrideAtCompileTime = StrideType::OuterStrideAtCompileTime == 0
? int(PlainObjectType::OuterStrideAtCompileTime)
: int(StrideType::OuterStrideAtCompileTime),
Alignment = int(MapOptions)&int(AlignedMask),
Flags0 = TraitsBase::Flags & (~NestByRefBit),
Flags = is_lvalue<PlainObjectType>::value ? int(Flags0) : (int(Flags0) & ~LvalueBit)
};
private:
enum { Options }; // Expressions don't have Options
};
}
/** \class Map
* \ingroup Core_Module
*
* \brief A matrix or vector expression mapping an existing array of data.
*
* \tparam PlainObjectType the equivalent matrix type of the mapped data
* \tparam MapOptions specifies the pointer alignment in bytes. It can be: \c #Aligned128, , \c #Aligned64, \c #Aligned32, \c #Aligned16, \c #Aligned8 or \c #Unaligned.
* The default is \c #Unaligned.
* \tparam StrideType optionally specifies strides. By default, Map assumes the memory layout
* of an ordinary, contiguous array. This can be overridden by specifying strides.
* The type passed here must be a specialization of the Stride template, see examples below.
*
* This class represents a matrix or vector expression mapping an existing array of data.
* It can be used to let Eigen interface without any overhead with non-Eigen data structures,
* such as plain C arrays or structures from other libraries. By default, it assumes that the
* data is laid out contiguously in memory. You can however override this by explicitly specifying
* inner and outer strides.
*
* Here's an example of simply mapping a contiguous array as a \ref TopicStorageOrders "column-major" matrix:
* \include Map_simple.cpp
* Output: \verbinclude Map_simple.out
*
* If you need to map non-contiguous arrays, you can do so by specifying strides:
*
* Here's an example of mapping an array as a vector, specifying an inner stride, that is, the pointer
* increment between two consecutive coefficients. Here, we're specifying the inner stride as a compile-time
* fixed value.
* \include Map_inner_stride.cpp
* Output: \verbinclude Map_inner_stride.out
*
* Here's an example of mapping an array while specifying an outer stride. Here, since we're mapping
* as a column-major matrix, 'outer stride' means the pointer increment between two consecutive columns.
* Here, we're specifying the outer stride as a runtime parameter. Note that here \c OuterStride<> is
* a short version of \c OuterStride<Dynamic> because the default template parameter of OuterStride
* is \c Dynamic
* \include Map_outer_stride.cpp
* Output: \verbinclude Map_outer_stride.out
*
* For more details and for an example of specifying both an inner and an outer stride, see class Stride.
*
* \b Tip: to change the array of data mapped by a Map object, you can use the C++
* placement new syntax:
*
* Example: \include Map_placement_new.cpp
* Output: \verbinclude Map_placement_new.out
*
* This class is the return type of PlainObjectBase::Map() but can also be used directly.
*
* \sa PlainObjectBase::Map(), \ref TopicStorageOrders
*/
template<typename PlainObjectType, int MapOptions, typename StrideType> class Map
: public MapBase<Map<PlainObjectType, MapOptions, StrideType> >
{
public:
typedef MapBase<Map> Base;
EIGEN_DENSE_PUBLIC_INTERFACE(Map)
typedef typename Base::PointerType PointerType;
typedef PointerType PointerArgType;
EIGEN_DEVICE_FUNC
inline PointerType cast_to_pointer_type(PointerArgType ptr) { return ptr; }
EIGEN_DEVICE_FUNC
inline Index innerStride() const
{
return StrideType::InnerStrideAtCompileTime != 0 ? m_stride.inner() : 1;
}
EIGEN_DEVICE_FUNC
inline Index outerStride() const
{
return StrideType::OuterStrideAtCompileTime != 0 ? m_stride.outer()
: IsVectorAtCompileTime ? this->size()
: int(Flags)&RowMajorBit ? this->cols()
: this->rows();
}
/** Constructor in the fixed-size case.
*
* \param dataPtr pointer to the array to map
* \param stride optional Stride object, passing the strides.
*/
EIGEN_DEVICE_FUNC
explicit inline Map(PointerArgType dataPtr, const StrideType& stride = StrideType())
: Base(cast_to_pointer_type(dataPtr)), m_stride(stride)
{
PlainObjectType::Base::_check_template_params();
}
/** Constructor in the dynamic-size vector case.
*
* \param dataPtr pointer to the array to map
* \param size the size of the vector expression
* \param stride optional Stride object, passing the strides.
*/
EIGEN_DEVICE_FUNC
inline Map(PointerArgType dataPtr, Index size, const StrideType& stride = StrideType())
: Base(cast_to_pointer_type(dataPtr), size), m_stride(stride)
{
PlainObjectType::Base::_check_template_params();
}
/** Constructor in the dynamic-size matrix case.
*
* \param dataPtr pointer to the array to map
* \param rows the number of rows of the matrix expression
* \param cols the number of columns of the matrix expression
* \param stride optional Stride object, passing the strides.
*/
EIGEN_DEVICE_FUNC
inline Map(PointerArgType dataPtr, Index rows, Index cols, const StrideType& stride = StrideType())
: Base(cast_to_pointer_type(dataPtr), rows, cols), m_stride(stride)
{
PlainObjectType::Base::_check_template_params();
}
EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Map)
protected:
StrideType m_stride;
};
} // end namespace Eigen
#endif // EIGEN_MAP_H
external/eigen3/Eigen/src/Core/MapBase.h 0 → 100644
View file @ a394b22a
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2007-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_MAPBASE_H
#define EIGEN_MAPBASE_H
#define EIGEN_STATIC_ASSERT_INDEX_BASED_ACCESS(Derived) \
EIGEN_STATIC_ASSERT((int(internal::evaluator<Derived>::Flags) & LinearAccessBit) || Derived::IsVectorAtCompileTime, \
YOU_ARE_TRYING_TO_USE_AN_INDEX_BASED_ACCESSOR_ON_AN_EXPRESSION_THAT_DOES_NOT_SUPPORT_THAT)
namespace Eigen {
/** \ingroup Core_Module
*
* \brief Base class for dense Map and Block expression with direct access
*
* This base class provides the const low-level accessors (e.g. coeff, coeffRef) of dense
* Map and Block objects with direct access.
* Typical users do not have to directly deal with this class.
*
* This class can be extended by through the macro plugin \c EIGEN_MAPBASE_PLUGIN.
* See \link TopicCustomizing_Plugins customizing Eigen \endlink for details.
*
* The \c Derived class has to provide the following two methods describing the memory layout:
* \code Index innerStride() const; \endcode
* \code Index outerStride() const; \endcode
*
* \sa class Map, class Block
*/
template<typename Derived> class MapBase<Derived, ReadOnlyAccessors>
: public internal::dense_xpr_base<Derived>::type
{
public:
typedef typename internal::dense_xpr_base<Derived>::type Base;
enum {
RowsAtCompileTime = internal::traits<Derived>::RowsAtCompileTime,
ColsAtCompileTime = internal::traits<Derived>::ColsAtCompileTime,
SizeAtCompileTime = Base::SizeAtCompileTime
};
typedef typename internal::traits<Derived>::StorageKind StorageKind;
typedef typename internal::traits<Derived>::Scalar Scalar;
typedef typename internal::packet_traits<Scalar>::type PacketScalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef typename internal::conditional<
bool(internal::is_lvalue<Derived>::value),
Scalar *,
const Scalar *>::type
PointerType;
using Base::derived;
// using Base::RowsAtCompileTime;
// using Base::ColsAtCompileTime;
// using Base::SizeAtCompileTime;
using Base::MaxRowsAtCompileTime;
using Base::MaxColsAtCompileTime;
using Base::MaxSizeAtCompileTime;
using Base::IsVectorAtCompileTime;
using Base::Flags;
using Base::IsRowMajor;
using Base::rows;
using Base::cols;
using Base::size;
using Base::coeff;
using Base::coeffRef;
using Base::lazyAssign;
using Base::eval;
using Base::innerStride;
using Base::outerStride;
using Base::rowStride;
using Base::colStride;
// bug 217 - compile error on ICC 11.1
using Base::operator=;
typedef typename Base::CoeffReturnType CoeffReturnType;
/** \copydoc DenseBase::rows() */
EIGEN_DEVICE_FUNC inline Index rows() const { return m_rows.value(); }
/** \copydoc DenseBase::cols() */
EIGEN_DEVICE_FUNC inline Index cols() const { return m_cols.value(); }
/** Returns a pointer to the first coefficient of the matrix or vector.
*
* \note When addressing this data, make sure to honor the strides returned by innerStride() and outerStride().
*
* \sa innerStride(), outerStride()
*/
EIGEN_DEVICE_FUNC inline const Scalar* data() const { return m_data; }
/** \copydoc PlainObjectBase::coeff(Index,Index) const */
EIGEN_DEVICE_FUNC
inline const Scalar& coeff(Index rowId, Index colId) const
{
return m_data[colId * colStride() + rowId * rowStride()];
}
/** \copydoc PlainObjectBase::coeff(Index) const */
EIGEN_DEVICE_FUNC
inline const Scalar& coeff(Index index) const
{
EIGEN_STATIC_ASSERT_INDEX_BASED_ACCESS(Derived)
return m_data[index * innerStride()];
}
/** \copydoc PlainObjectBase::coeffRef(Index,Index) const */
EIGEN_DEVICE_FUNC
inline const Scalar& coeffRef(Index rowId, Index colId) const
{
return this->m_data[colId * colStride() + rowId * rowStride()];
}
/** \copydoc PlainObjectBase::coeffRef(Index) const */
EIGEN_DEVICE_FUNC
inline const Scalar& coeffRef(Index index) const
{
EIGEN_STATIC_ASSERT_INDEX_BASED_ACCESS(Derived)
return this->m_data[index * innerStride()];
}
/** \internal */
template<int LoadMode>
inline PacketScalar packet(Index rowId, Index colId) const
{
return internal::ploadt<PacketScalar, LoadMode>
(m_data + (colId * colStride() + rowId * rowStride()));
}
/** \internal */
template<int LoadMode>
inline PacketScalar packet(Index index) const
{
EIGEN_STATIC_ASSERT_INDEX_BASED_ACCESS(Derived)
return internal::ploadt<PacketScalar, LoadMode>(m_data + index * innerStride());
}
/** \internal Constructor for fixed size matrices or vectors */
EIGEN_DEVICE_FUNC
explicit inline MapBase(PointerType dataPtr) : m_data(dataPtr), m_rows(RowsAtCompileTime), m_cols(ColsAtCompileTime)
{
EIGEN_STATIC_ASSERT_FIXED_SIZE(Derived)
checkSanity<Derived>();
}
/** \internal Constructor for dynamically sized vectors */
EIGEN_DEVICE_FUNC
inline MapBase(PointerType dataPtr, Index vecSize)
: m_data(dataPtr),
m_rows(RowsAtCompileTime == Dynamic ? vecSize : Index(RowsAtCompileTime)),
m_cols(ColsAtCompileTime == Dynamic ? vecSize : Index(ColsAtCompileTime))
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
eigen_assert(vecSize >= 0);
eigen_assert(dataPtr == 0 || SizeAtCompileTime == Dynamic || SizeAtCompileTime == vecSize);
checkSanity<Derived>();
}
/** \internal Constructor for dynamically sized matrices */
EIGEN_DEVICE_FUNC
inline MapBase(PointerType dataPtr, Index rows, Index cols)
: m_data(dataPtr), m_rows(rows), m_cols(cols)
{
eigen_assert( (dataPtr == 0)
|| ( rows >= 0 && (RowsAtCompileTime == Dynamic || RowsAtCompileTime == rows)
&& cols >= 0 && (ColsAtCompileTime == Dynamic || ColsAtCompileTime == cols)));
checkSanity<Derived>();
}
#ifdef EIGEN_MAPBASE_PLUGIN
#include EIGEN_MAPBASE_PLUGIN
#endif
protected:
template<typename T>
EIGEN_DEVICE_FUNC
void checkSanity(typename internal::enable_if<(internal::traits<T>::Alignment>0),void*>::type = 0) const
{
#if EIGEN_MAX_ALIGN_BYTES>0
eigen_assert(( ((internal::UIntPtr(m_data) % internal::traits<Derived>::Alignment) == 0)
|| (cols() * rows() * innerStride() * sizeof(Scalar)) < internal::traits<Derived>::Alignment ) && "data is not aligned");
#endif
}
template<typename T>
EIGEN_DEVICE_FUNC
void checkSanity(typename internal::enable_if<internal::traits<T>::Alignment==0,void*>::type = 0) const
{}
PointerType m_data;
const internal::variable_if_dynamic<Index, RowsAtCompileTime> m_rows;
const internal::variable_if_dynamic<Index, ColsAtCompileTime> m_cols;
};
/** \ingroup Core_Module
*
* \brief Base class for non-const dense Map and Block expression with direct access
*
* This base class provides the non-const low-level accessors (e.g. coeff and coeffRef) of
* dense Map and Block objects with direct access.
* It inherits MapBase<Derived, ReadOnlyAccessors> which defines the const variant for reading specific entries.
*
* \sa class Map, class Block
*/
template<typename Derived> class MapBase<Derived, WriteAccessors>
: public MapBase<Derived, ReadOnlyAccessors>
{
typedef MapBase<Derived, ReadOnlyAccessors> ReadOnlyMapBase;
public:
typedef MapBase<Derived, ReadOnlyAccessors> Base;
typedef typename Base::Scalar Scalar;
typedef typename Base::PacketScalar PacketScalar;
typedef typename Base::StorageIndex StorageIndex;
typedef typename Base::PointerType PointerType;
using Base::derived;
using Base::rows;
using Base::cols;
using Base::size;
using Base::coeff;
using Base::coeffRef;
using Base::innerStride;
using Base::outerStride;
using Base::rowStride;
using Base::colStride;
typedef typename internal::conditional<
internal::is_lvalue<Derived>::value,
Scalar,
const Scalar
>::type ScalarWithConstIfNotLvalue;
EIGEN_DEVICE_FUNC
inline const Scalar* data() const { return this->m_data; }
EIGEN_DEVICE_FUNC
inline ScalarWithConstIfNotLvalue* data() { return this->m_data; } // no const-cast here so non-const-correct code will give a compile error
EIGEN_DEVICE_FUNC
inline ScalarWithConstIfNotLvalue& coeffRef(Index row, Index col)
{
return this->m_data[col * colStride() + row * rowStride()];
}
EIGEN_DEVICE_FUNC
inline ScalarWithConstIfNotLvalue& coeffRef(Index index)
{
EIGEN_STATIC_ASSERT_INDEX_BASED_ACCESS(Derived)
return this->m_data[index * innerStride()];
}
template<int StoreMode>
inline void writePacket(Index row, Index col, const PacketScalar& val)
{
internal::pstoret<Scalar, PacketScalar, StoreMode>
(this->m_data + (col * colStride() + row * rowStride()), val);
}
template<int StoreMode>
inline void writePacket(Index index, const PacketScalar& val)
{
EIGEN_STATIC_ASSERT_INDEX_BASED_ACCESS(Derived)
internal::pstoret<Scalar, PacketScalar, StoreMode>
(this->m_data + index * innerStride(), val);
}
EIGEN_DEVICE_FUNC explicit inline MapBase(PointerType dataPtr) : Base(dataPtr) {}
EIGEN_DEVICE_FUNC inline MapBase(PointerType dataPtr, Index vecSize) : Base(dataPtr, vecSize) {}
EIGEN_DEVICE_FUNC inline MapBase(PointerType dataPtr, Index rows, Index cols) : Base(dataPtr, rows, cols) {}
EIGEN_DEVICE_FUNC
Derived& operator=(const MapBase& other)
{
ReadOnlyMapBase::Base::operator=(other);
return derived();
}
// In theory we could simply refer to Base:Base::operator=, but MSVC does not like Base::Base,
// see bugs 821 and 920.
using ReadOnlyMapBase::Base::operator=;
};
#undef EIGEN_STATIC_ASSERT_INDEX_BASED_ACCESS
} // end namespace Eigen
#endif // EIGEN_MAPBASE_H
external/eigen3/Eigen/src/Core/MathFunctions.h 0 → 100644
View file @ a394b22a
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2006-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_MATHFUNCTIONS_H
#define EIGEN_MATHFUNCTIONS_H
// source: http://www.geom.uiuc.edu/~huberty/math5337/groupe/digits.html
// TODO this should better be moved to NumTraits
#define EIGEN_PI 3.141592653589793238462643383279502884197169399375105820974944592307816406L
namespace Eigen {
// On WINCE, std::abs is defined for int only, so let's defined our own overloads:
// This issue has been confirmed with MSVC 2008 only, but the issue might exist for more recent versions too.
#if EIGEN_OS_WINCE && EIGEN_COMP_MSVC && EIGEN_COMP_MSVC<=1500
long abs(long x) { return (labs(x)); }
double abs(double x) { return (fabs(x)); }
float abs(float x) { return (fabsf(x)); }
long double abs(long double x) { return (fabsl(x)); }
#endif
namespace internal {
/** \internal \class global_math_functions_filtering_base
*
* What it does:
* Defines a typedef 'type' as follows:
* - if type T has a member typedef Eigen_BaseClassForSpecializationOfGlobalMathFuncImpl, then
* global_math_functions_filtering_base<T>::type is a typedef for it.
* - otherwise, global_math_functions_filtering_base<T>::type is a typedef for T.
*
* How it's used:
* To allow to defined the global math functions (like sin...) in certain cases, like the Array expressions.
* When you do sin(array1+array2), the object array1+array2 has a complicated expression type, all what you want to know
* is that it inherits ArrayBase. So we implement a partial specialization of sin_impl for ArrayBase<Derived>.
* So we must make sure to use sin_impl<ArrayBase<Derived> > and not sin_impl<Derived>, otherwise our partial specialization
* won't be used. How does sin know that? That's exactly what global_math_functions_filtering_base tells it.
*
* How it's implemented:
* SFINAE in the style of enable_if. Highly susceptible of breaking compilers. With GCC, it sure does work, but if you replace
* the typename dummy by an integer template parameter, it doesn't work anymore!
*/
template<typename T, typename dummy = void>
struct global_math_functions_filtering_base
{
typedef T type;
};
template<typename T> struct always_void { typedef void type; };
template<typename T>
struct global_math_functions_filtering_base
<T,
typename always_void<typename T::Eigen_BaseClassForSpecializationOfGlobalMathFuncImpl>::type
>
{
typedef typename T::Eigen_BaseClassForSpecializationOfGlobalMathFuncImpl type;
};
#define EIGEN_MATHFUNC_IMPL(func, scalar) Eigen::internal::func##_impl<typename Eigen::internal::global_math_functions_filtering_base<scalar>::type>
#define EIGEN_MATHFUNC_RETVAL(func, scalar) typename Eigen::internal::func##_retval<typename Eigen::internal::global_math_functions_filtering_base<scalar>::type>::type
/****************************************************************************
* Implementation of real *
****************************************************************************/
template<typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
struct real_default_impl
{
typedef typename NumTraits<Scalar>::Real RealScalar;
EIGEN_DEVICE_FUNC
static inline RealScalar run(const Scalar& x)
{
return x;
}
};
template<typename Scalar>
struct real_default_impl<Scalar,true>
{
typedef typename NumTraits<Scalar>::Real RealScalar;
EIGEN_DEVICE_FUNC
static inline RealScalar run(const Scalar& x)
{
using std::real;
return real(x);
}
};
template<typename Scalar> struct real_impl : real_default_impl<Scalar> {};
#ifdef __CUDA_ARCH__
template<typename T>
struct real_impl<std::complex<T> >
{
typedef T RealScalar;
EIGEN_DEVICE_FUNC
static inline T run(const std::complex<T>& x)
{
return x.real();
}
};
#endif
template<typename Scalar>
struct real_retval
{
typedef typename NumTraits<Scalar>::Real type;
};
/****************************************************************************
* Implementation of imag *
****************************************************************************/
template<typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
struct imag_default_impl
{
typedef typename NumTraits<Scalar>::Real RealScalar;
EIGEN_DEVICE_FUNC
static inline RealScalar run(const Scalar&)
{
return RealScalar(0);
}
};
template<typename Scalar>
struct imag_default_impl<Scalar,true>
{
typedef typename NumTraits<Scalar>::Real RealScalar;
EIGEN_DEVICE_FUNC
static inline RealScalar run(const Scalar& x)
{
using std::imag;
return imag(x);
}
};
template<typename Scalar> struct imag_impl : imag_default_impl<Scalar> {};
#ifdef __CUDA_ARCH__
template<typename T>
struct imag_impl<std::complex<T> >
{
typedef T RealScalar;
EIGEN_DEVICE_FUNC
static inline T run(const std::complex<T>& x)
{
return x.imag();
}
};
#endif
template<typename Scalar>
struct imag_retval
{
typedef typename NumTraits<Scalar>::Real type;
};
/****************************************************************************
* Implementation of real_ref *
****************************************************************************/
template<typename Scalar>
struct real_ref_impl
{
typedef typename NumTraits<Scalar>::Real RealScalar;
EIGEN_DEVICE_FUNC
static inline RealScalar& run(Scalar& x)
{
return reinterpret_cast<RealScalar*>(&x)[0];
}
EIGEN_DEVICE_FUNC
static inline const RealScalar& run(const Scalar& x)
{
return reinterpret_cast<const RealScalar*>(&x)[0];
}
};
template<typename Scalar>
struct real_ref_retval
{
typedef typename NumTraits<Scalar>::Real & type;
};
/****************************************************************************
* Implementation of imag_ref *
****************************************************************************/
template<typename Scalar, bool IsComplex>
struct imag_ref_default_impl
{
typedef typename NumTraits<Scalar>::Real RealScalar;
EIGEN_DEVICE_FUNC
static inline RealScalar& run(Scalar& x)
{
return reinterpret_cast<RealScalar*>(&x)[1];
}
EIGEN_DEVICE_FUNC
static inline const RealScalar& run(const Scalar& x)
{
return reinterpret_cast<RealScalar*>(&x)[1];
}
};
template<typename Scalar>
struct imag_ref_default_impl<Scalar, false>
{
EIGEN_DEVICE_FUNC
static inline Scalar run(Scalar&)
{
return Scalar(0);
}
EIGEN_DEVICE_FUNC
static inline const Scalar run(const Scalar&)
{
return Scalar(0);
}
};
template<typename Scalar>
struct imag_ref_impl : imag_ref_default_impl<Scalar, NumTraits<Scalar>::IsComplex> {};
template<typename Scalar>
struct imag_ref_retval
{
typedef typename NumTraits<Scalar>::Real & type;
};
/****************************************************************************
* Implementation of conj *
****************************************************************************/
template<typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
struct conj_impl
{
EIGEN_DEVICE_FUNC
static inline Scalar run(const Scalar& x)
{
return x;
}
};
template<typename Scalar>
struct conj_impl<Scalar,true>
{
EIGEN_DEVICE_FUNC
static inline Scalar run(const Scalar& x)
{
using std::conj;
return conj(x);
}
};
template<typename Scalar>
struct conj_retval
{
typedef Scalar type;
};
/****************************************************************************
* Implementation of abs2 *
****************************************************************************/
template<typename Scalar,bool IsComplex>
struct abs2_impl_default
{
typedef typename NumTraits<Scalar>::Real RealScalar;
EIGEN_DEVICE_FUNC
static inline RealScalar run(const Scalar& x)
{
return x*x;
}
};
template<typename Scalar>
struct abs2_impl_default<Scalar, true> // IsComplex
{
typedef typename NumTraits<Scalar>::Real RealScalar;
EIGEN_DEVICE_FUNC
static inline RealScalar run(const Scalar& x)
{
return real(x)*real(x) + imag(x)*imag(x);
}
};
template<typename Scalar>
struct abs2_impl
{
typedef typename NumTraits<Scalar>::Real RealScalar;
EIGEN_DEVICE_FUNC
static inline RealScalar run(const Scalar& x)
{
return abs2_impl_default<Scalar,NumTraits<Scalar>::IsComplex>::run(x);
}
};
template<typename Scalar>
struct abs2_retval
{
typedef typename NumTraits<Scalar>::Real type;
};
/****************************************************************************
* Implementation of norm1 *
****************************************************************************/
template<typename Scalar, bool IsComplex>
struct norm1_default_impl
{
typedef typename NumTraits<Scalar>::Real RealScalar;
EIGEN_DEVICE_FUNC
static inline RealScalar run(const Scalar& x)
{
EIGEN_USING_STD_MATH(abs);
return abs(real(x)) + abs(imag(x));
}
};
template<typename Scalar>
struct norm1_default_impl<Scalar, false>
{
EIGEN_DEVICE_FUNC
static inline Scalar run(const Scalar& x)
{
EIGEN_USING_STD_MATH(abs);
return abs(x);
}
};
template<typename Scalar>
struct norm1_impl : norm1_default_impl<Scalar, NumTraits<Scalar>::IsComplex> {};
template<typename Scalar>
struct norm1_retval
{
typedef typename NumTraits<Scalar>::Real type;
};
/****************************************************************************
* Implementation of hypot *
****************************************************************************/
template<typename Scalar>
struct hypot_impl
{
typedef typename NumTraits<Scalar>::Real RealScalar;
static inline RealScalar run(const Scalar& x, const Scalar& y)
{
EIGEN_USING_STD_MATH(abs);
EIGEN_USING_STD_MATH(sqrt);
RealScalar _x = abs(x);
RealScalar _y = abs(y);
Scalar p, qp;
if(_x>_y)
{
p = _x;
qp = _y / p;
}
else
{
p = _y;
qp = _x / p;
}
if(p==RealScalar(0)) return RealScalar(0);
return p * sqrt(RealScalar(1) + qp*qp);
}
};
template<typename Scalar>
struct hypot_retval
{
typedef typename NumTraits<Scalar>::Real type;
};
/****************************************************************************
* Implementation of cast *
****************************************************************************/
template<typename OldType, typename NewType>
struct cast_impl
{
EIGEN_DEVICE_FUNC
static inline NewType run(const OldType& x)
{
return static_cast<NewType>(x);
}
};
// here, for once, we're plainly returning NewType: we don't want cast to do weird things.
template<typename OldType, typename NewType>
EIGEN_DEVICE_FUNC
inline NewType cast(const OldType& x)
{
return cast_impl<OldType, NewType>::run(x);
}
/****************************************************************************
* Implementation of round *
****************************************************************************/
#if EIGEN_HAS_CXX11_MATH
template<typename Scalar>
struct round_impl {
static inline Scalar run(const Scalar& x)
{
EIGEN_STATIC_ASSERT((!NumTraits<Scalar>::IsComplex), NUMERIC_TYPE_MUST_BE_REAL)
using std::round;
return round(x);
}
};
#else
template<typename Scalar>
struct round_impl
{
static inline Scalar run(const Scalar& x)
{
EIGEN_STATIC_ASSERT((!NumTraits<Scalar>::IsComplex), NUMERIC_TYPE_MUST_BE_REAL)
EIGEN_USING_STD_MATH(floor);
EIGEN_USING_STD_MATH(ceil);
return (x > Scalar(0)) ? floor(x + Scalar(0.5)) : ceil(x - Scalar(0.5));
}
};
#endif
template<typename Scalar>
struct round_retval
{
typedef Scalar type;
};
/****************************************************************************
* Implementation of arg *
****************************************************************************/
#if EIGEN_HAS_CXX11_MATH
template<typename Scalar>
struct arg_impl {
static inline Scalar run(const Scalar& x)
{
EIGEN_USING_STD_MATH(arg);
return arg(x);
}
};
#else
template<typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
struct arg_default_impl
{
typedef typename NumTraits<Scalar>::Real RealScalar;
EIGEN_DEVICE_FUNC
static inline RealScalar run(const Scalar& x)
{
return (x < Scalar(0)) ? Scalar(EIGEN_PI) : Scalar(0); }
};
template<typename Scalar>
struct arg_default_impl<Scalar,true>
{
typedef typename NumTraits<Scalar>::Real RealScalar;
EIGEN_DEVICE_FUNC
static inline RealScalar run(const Scalar& x)
{
EIGEN_USING_STD_MATH(arg);
return arg(x);
}
};
template<typename Scalar> struct arg_impl : arg_default_impl<Scalar> {};
#endif
template<typename Scalar>
struct arg_retval
{
typedef typename NumTraits<Scalar>::Real type;
};
/****************************************************************************
* Implementation of log1p *
****************************************************************************/
namespace std_fallback {
// fallback log1p implementation in case there is no log1p(Scalar) function in namespace of Scalar,
// or that there is no suitable std::log1p function available
template<typename Scalar>
EIGEN_DEVICE_FUNC inline Scalar log1p(const Scalar& x) {
EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
typedef typename NumTraits<Scalar>::Real RealScalar;
EIGEN_USING_STD_MATH(log);
Scalar x1p = RealScalar(1) + x;
return ( x1p == Scalar(1) ) ? x : x * ( log(x1p) / (x1p - RealScalar(1)) );
}
}
template<typename Scalar>
struct log1p_impl {
static inline Scalar run(const Scalar& x)
{
EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
#if EIGEN_HAS_CXX11_MATH
using std::log1p;
#endif
using std_fallback::log1p;
return log1p(x);
}
};
template<typename Scalar>
struct log1p_retval
{
typedef Scalar type;
};
/****************************************************************************
* Implementation of pow *
****************************************************************************/
template<typename ScalarX,typename ScalarY, bool IsInteger = NumTraits<ScalarX>::IsInteger&&NumTraits<ScalarY>::IsInteger>
struct pow_impl
{
//typedef Scalar retval;
typedef typename ScalarBinaryOpTraits<ScalarX,ScalarY,internal::scalar_pow_op<ScalarX,ScalarY> >::ReturnType result_type;
static EIGEN_DEVICE_FUNC inline result_type run(const ScalarX& x, const ScalarY& y)
{
EIGEN_USING_STD_MATH(pow);
return pow(x, y);
}
};
template<typename ScalarX,typename ScalarY>
struct pow_impl<ScalarX,ScalarY, true>
{
typedef ScalarX result_type;
static EIGEN_DEVICE_FUNC inline ScalarX run(ScalarX x, ScalarY y)
{
ScalarX res(1);
eigen_assert(!NumTraits<ScalarY>::IsSigned || y >= 0);
if(y & 1) res *= x;
y >>= 1;
while(y)
{
x *= x;
if(y&1) res *= x;
y >>= 1;
}
return res;
}
};
/****************************************************************************
* Implementation of random *
****************************************************************************/
template<typename Scalar,
bool IsComplex,
bool IsInteger>
struct random_default_impl {};
template<typename Scalar>
struct random_impl : random_default_impl<Scalar, NumTraits<Scalar>::IsComplex, NumTraits<Scalar>::IsInteger> {};
template<typename Scalar>
struct random_retval
{
typedef Scalar type;
};
template<typename Scalar> inline EIGEN_MATHFUNC_RETVAL(random, Scalar) random(const Scalar& x, const Scalar& y);
template<typename Scalar> inline EIGEN_MATHFUNC_RETVAL(random, Scalar) random();
template<typename Scalar>
struct random_default_impl<Scalar, false, false>
{
static inline Scalar run(const Scalar& x, const Scalar& y)
{
return x + (y-x) * Scalar(std::rand()) / Scalar(RAND_MAX);
}
static inline Scalar run()
{
return run(Scalar(NumTraits<Scalar>::IsSigned ? -1 : 0), Scalar(1));
}
};
enum {
meta_floor_log2_terminate,
meta_floor_log2_move_up,
meta_floor_log2_move_down,
meta_floor_log2_bogus
};
template<unsigned int n, int lower, int upper> struct meta_floor_log2_selector
{
enum { middle = (lower + upper) / 2,
value = (upper <= lower + 1) ? int(meta_floor_log2_terminate)
: (n < (1 << middle)) ? int(meta_floor_log2_move_down)
: (n==0) ? int(meta_floor_log2_bogus)
: int(meta_floor_log2_move_up)
};
};
template<unsigned int n,
int lower = 0,
int upper = sizeof(unsigned int) * CHAR_BIT - 1,
int selector = meta_floor_log2_selector<n, lower, upper>::value>
struct meta_floor_log2 {};
template<unsigned int n, int lower, int upper>
struct meta_floor_log2<n, lower, upper, meta_floor_log2_move_down>
{
enum { value = meta_floor_log2<n, lower, meta_floor_log2_selector<n, lower, upper>::middle>::value };
};
template<unsigned int n, int lower, int upper>
struct meta_floor_log2<n, lower, upper, meta_floor_log2_move_up>
{
enum { value = meta_floor_log2<n, meta_floor_log2_selector<n, lower, upper>::middle, upper>::value };
};
template<unsigned int n, int lower, int upper>
struct meta_floor_log2<n, lower, upper, meta_floor_log2_terminate>
{
enum { value = (n >= ((unsigned int)(1) << (lower+1))) ? lower+1 : lower };
};
template<unsigned int n, int lower, int upper>
struct meta_floor_log2<n, lower, upper, meta_floor_log2_bogus>
{
// no value, error at compile time
};
template<typename Scalar>
struct random_default_impl<Scalar, false, true>
{
static inline Scalar run(const Scalar& x, const Scalar& y)
{
typedef typename conditional<NumTraits<Scalar>::IsSigned,std::ptrdiff_t,std::size_t>::type ScalarX;
if(y<x)
return x;
// the following difference might overflow on a 32 bits system,
// but since y>=x the result converted to an unsigned long is still correct.
std::size_t range = ScalarX(y)-ScalarX(x);
std::size_t offset = 0;
// rejection sampling
std::size_t divisor = 1;
std::size_t multiplier = 1;
if(range<RAND_MAX) divisor = (std::size_t(RAND_MAX)+1)/(range+1);
else multiplier = 1 + range/(std::size_t(RAND_MAX)+1);
do {
offset = (std::size_t(std::rand()) * multiplier) / divisor;
} while (offset > range);
return Scalar(ScalarX(x) + offset);
}
static inline Scalar run()
{
#ifdef EIGEN_MAKING_DOCS
return run(Scalar(NumTraits<Scalar>::IsSigned ? -10 : 0), Scalar(10));
#else
enum { rand_bits = meta_floor_log2<(unsigned int)(RAND_MAX)+1>::value,
scalar_bits = sizeof(Scalar) * CHAR_BIT,
shift = EIGEN_PLAIN_ENUM_MAX(0, int(rand_bits) - int(scalar_bits)),
offset = NumTraits<Scalar>::IsSigned ? (1 << (EIGEN_PLAIN_ENUM_MIN(rand_bits,scalar_bits)-1)) : 0
};
return Scalar((std::rand() >> shift) - offset);
#endif
}
};
template<typename Scalar>
struct random_default_impl<Scalar, true, false>
{
static inline Scalar run(const Scalar& x, const Scalar& y)
{
return Scalar(random(real(x), real(y)),
random(imag(x), imag(y)));
}
static inline Scalar run()
{
typedef typename NumTraits<Scalar>::Real RealScalar;
return Scalar(random<RealScalar>(), random<RealScalar>());
}
};
template<typename Scalar>
inline EIGEN_MATHFUNC_RETVAL(random, Scalar) random(const Scalar& x, const Scalar& y)
{
return EIGEN_MATHFUNC_IMPL(random, Scalar)::run(x, y);
}
template<typename Scalar>
inline EIGEN_MATHFUNC_RETVAL(random, Scalar) random()
{
return EIGEN_MATHFUNC_IMPL(random, Scalar)::run();
}
// Implementatin of is* functions
// std::is* do not work with fast-math and gcc, std::is* are available on MSVC 2013 and newer, as well as in clang.
#if (EIGEN_HAS_CXX11_MATH && !(EIGEN_COMP_GNUC_STRICT && __FINITE_MATH_ONLY__)) || (EIGEN_COMP_MSVC>=1800) || (EIGEN_COMP_CLANG)
#define EIGEN_USE_STD_FPCLASSIFY 1
#else
#define EIGEN_USE_STD_FPCLASSIFY 0
#endif
template<typename T>
EIGEN_DEVICE_FUNC
typename internal::enable_if<internal::is_integral<T>::value,bool>::type
isnan_impl(const T&) { return false; }
template<typename T>
EIGEN_DEVICE_FUNC
typename internal::enable_if<internal::is_integral<T>::value,bool>::type
isinf_impl(const T&) { return false; }
template<typename T>
EIGEN_DEVICE_FUNC
typename internal::enable_if<internal::is_integral<T>::value,bool>::type
isfinite_impl(const T&) { return true; }
template<typename T>
EIGEN_DEVICE_FUNC
typename internal::enable_if<(!internal::is_integral<T>::value)&&(!NumTraits<T>::IsComplex),bool>::type
isfinite_impl(const T& x)
{
#ifdef __CUDA_ARCH__
return (::isfinite)(x);
#elif EIGEN_USE_STD_FPCLASSIFY
using std::isfinite;
return isfinite EIGEN_NOT_A_MACRO (x);
#else
return x<=NumTraits<T>::highest() && x>=NumTraits<T>::lowest();
#endif
}
template<typename T>
EIGEN_DEVICE_FUNC
typename internal::enable_if<(!internal::is_integral<T>::value)&&(!NumTraits<T>::IsComplex),bool>::type
isinf_impl(const T& x)
{
#ifdef __CUDA_ARCH__
return (::isinf)(x);
#elif EIGEN_USE_STD_FPCLASSIFY
using std::isinf;
return isinf EIGEN_NOT_A_MACRO (x);
#else
return x>NumTraits<T>::highest() || x<NumTraits<T>::lowest();
#endif
}
template<typename T>
EIGEN_DEVICE_FUNC
typename internal::enable_if<(!internal::is_integral<T>::value)&&(!NumTraits<T>::IsComplex),bool>::type
isnan_impl(const T& x)
{
#ifdef __CUDA_ARCH__
return (::isnan)(x);
#elif EIGEN_USE_STD_FPCLASSIFY
using std::isnan;
return isnan EIGEN_NOT_A_MACRO (x);
#else
return x != x;
#endif
}
#if (!EIGEN_USE_STD_FPCLASSIFY)
#if EIGEN_COMP_MSVC
template<typename T> EIGEN_DEVICE_FUNC bool isinf_msvc_helper(T x)
{
return _fpclass(x)==_FPCLASS_NINF || _fpclass(x)==_FPCLASS_PINF;
}
//MSVC defines a _isnan builtin function, but for double only
EIGEN_DEVICE_FUNC inline bool isnan_impl(const long double& x) { return _isnan(x)!=0; }
EIGEN_DEVICE_FUNC inline bool isnan_impl(const double& x) { return _isnan(x)!=0; }
EIGEN_DEVICE_FUNC inline bool isnan_impl(const float& x) { return _isnan(x)!=0; }
EIGEN_DEVICE_FUNC inline bool isinf_impl(const long double& x) { return isinf_msvc_helper(x); }
EIGEN_DEVICE_FUNC inline bool isinf_impl(const double& x) { return isinf_msvc_helper(x); }
EIGEN_DEVICE_FUNC inline bool isinf_impl(const float& x) { return isinf_msvc_helper(x); }
#elif (defined __FINITE_MATH_ONLY__ && __FINITE_MATH_ONLY__ && EIGEN_COMP_GNUC)
#if EIGEN_GNUC_AT_LEAST(5,0)
#define EIGEN_TMP_NOOPT_ATTRIB EIGEN_DEVICE_FUNC inline __attribute__((optimize("no-finite-math-only")))
#else
// NOTE the inline qualifier and noinline attribute are both needed: the former is to avoid linking issue (duplicate symbol),
// while the second prevent too aggressive optimizations in fast-math mode:
#define EIGEN_TMP_NOOPT_ATTRIB EIGEN_DEVICE_FUNC inline __attribute__((noinline,optimize("no-finite-math-only")))
#endif
template<> EIGEN_TMP_NOOPT_ATTRIB bool isnan_impl(const long double& x) { return __builtin_isnan(x); }
template<> EIGEN_TMP_NOOPT_ATTRIB bool isnan_impl(const double& x) { return __builtin_isnan(x); }
template<> EIGEN_TMP_NOOPT_ATTRIB bool isnan_impl(const float& x) { return __builtin_isnan(x); }
template<> EIGEN_TMP_NOOPT_ATTRIB bool isinf_impl(const double& x) { return __builtin_isinf(x); }
template<> EIGEN_TMP_NOOPT_ATTRIB bool isinf_impl(const float& x) { return __builtin_isinf(x); }
template<> EIGEN_TMP_NOOPT_ATTRIB bool isinf_impl(const long double& x) { return __builtin_isinf(x); }
#undef EIGEN_TMP_NOOPT_ATTRIB
#endif
#endif
// The following overload are defined at the end of this file
template<typename T> EIGEN_DEVICE_FUNC bool isfinite_impl(const std::complex<T>& x);
template<typename T> EIGEN_DEVICE_FUNC bool isnan_impl(const std::complex<T>& x);
template<typename T> EIGEN_DEVICE_FUNC bool isinf_impl(const std::complex<T>& x);
template<typename T> T generic_fast_tanh_float(const T& a_x);
} // end namespace internal
/****************************************************************************
* Generic math functions *
****************************************************************************/
namespace numext {
#ifndef __CUDA_ARCH__
template<typename T>
EIGEN_DEVICE_FUNC
EIGEN_ALWAYS_INLINE T mini(const T& x, const T& y)
{
EIGEN_USING_STD_MATH(min);
return min EIGEN_NOT_A_MACRO (x,y);
}
template<typename T>
EIGEN_DEVICE_FUNC
EIGEN_ALWAYS_INLINE T maxi(const T& x, const T& y)
{
EIGEN_USING_STD_MATH(max);
return max EIGEN_NOT_A_MACRO (x,y);
}
#else
template<typename T>
EIGEN_DEVICE_FUNC
EIGEN_ALWAYS_INLINE T mini(const T& x, const T& y)
{
return y < x ? y : x;
}
template<>
EIGEN_DEVICE_FUNC
EIGEN_ALWAYS_INLINE float mini(const float& x, const float& y)
{
return fminf(x, y);
}
template<typename T>
EIGEN_DEVICE_FUNC
EIGEN_ALWAYS_INLINE T maxi(const T& x, const T& y)
{
return x < y ? y : x;
}
template<>
EIGEN_DEVICE_FUNC
EIGEN_ALWAYS_INLINE float maxi(const float& x, const float& y)
{
return fmaxf(x, y);
}
#endif
template<typename Scalar>
EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(real, Scalar) real(const Scalar& x)
{
return EIGEN_MATHFUNC_IMPL(real, Scalar)::run(x);
}
template<typename Scalar>
EIGEN_DEVICE_FUNC
inline typename internal::add_const_on_value_type< EIGEN_MATHFUNC_RETVAL(real_ref, Scalar) >::type real_ref(const Scalar& x)
{
return internal::real_ref_impl<Scalar>::run(x);
}
template<typename Scalar>
EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(real_ref, Scalar) real_ref(Scalar& x)
{
return EIGEN_MATHFUNC_IMPL(real_ref, Scalar)::run(x);
}
template<typename Scalar>
EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(imag, Scalar) imag(const Scalar& x)
{
return EIGEN_MATHFUNC_IMPL(imag, Scalar)::run(x);
}
template<typename Scalar>
EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(arg, Scalar) arg(const Scalar& x)
{
return EIGEN_MATHFUNC_IMPL(arg, Scalar)::run(x);
}
template<typename Scalar>
EIGEN_DEVICE_FUNC
inline typename internal::add_const_on_value_type< EIGEN_MATHFUNC_RETVAL(imag_ref, Scalar) >::type imag_ref(const Scalar& x)
{
return internal::imag_ref_impl<Scalar>::run(x);
}
template<typename Scalar>
EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(imag_ref, Scalar) imag_ref(Scalar& x)
{
return EIGEN_MATHFUNC_IMPL(imag_ref, Scalar)::run(x);
}
template<typename Scalar>
EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(conj, Scalar) conj(const Scalar& x)
{
return EIGEN_MATHFUNC_IMPL(conj, Scalar)::run(x);
}
template<typename Scalar>
EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(abs2, Scalar) abs2(const Scalar& x)
{
return EIGEN_MATHFUNC_IMPL(abs2, Scalar)::run(x);
}
template<typename Scalar>
EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(norm1, Scalar) norm1(const Scalar& x)
{
return EIGEN_MATHFUNC_IMPL(norm1, Scalar)::run(x);
}
template<typename Scalar>
EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(hypot, Scalar) hypot(const Scalar& x, const Scalar& y)
{
return EIGEN_MATHFUNC_IMPL(hypot, Scalar)::run(x, y);
}
template<typename Scalar>
EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(log1p, Scalar) log1p(const Scalar& x)
{
return EIGEN_MATHFUNC_IMPL(log1p, Scalar)::run(x);
}
#ifdef __CUDACC__
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
float log1p(const float &x) { return ::log1pf(x); }
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
double log1p(const double &x) { return ::log1p(x); }
#endif
template<typename ScalarX,typename ScalarY>
EIGEN_DEVICE_FUNC
inline typename internal::pow_impl<ScalarX,ScalarY>::result_type pow(const ScalarX& x, const ScalarY& y)
{
return internal::pow_impl<ScalarX,ScalarY>::run(x, y);
}
template<typename T> EIGEN_DEVICE_FUNC bool (isnan) (const T &x) { return internal::isnan_impl(x); }
template<typename T> EIGEN_DEVICE_FUNC bool (isinf) (const T &x) { return internal::isinf_impl(x); }
template<typename T> EIGEN_DEVICE_FUNC bool (isfinite)(const T &x) { return internal::isfinite_impl(x); }
template<typename Scalar>
EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(round, Scalar) round(const Scalar& x)
{
return EIGEN_MATHFUNC_IMPL(round, Scalar)::run(x);
}
template<typename T>
EIGEN_DEVICE_FUNC
T (floor)(const T& x)
{
EIGEN_USING_STD_MATH(floor);
return floor(x);
}
#ifdef __CUDACC__
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
float floor(const float &x) { return ::floorf(x); }
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
double floor(const double &x) { return ::floor(x); }
#endif
template<typename T>
EIGEN_DEVICE_FUNC
T (ceil)(const T& x)
{
EIGEN_USING_STD_MATH(ceil);
return ceil(x);
}
#ifdef __CUDACC__
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
float ceil(const float &x) { return ::ceilf(x); }
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
double ceil(const double &x) { return ::ceil(x); }
#endif
/** Log base 2 for 32 bits positive integers.
* Conveniently returns 0 for x==0. */
inline int log2(int x)
{
eigen_assert(x>=0);
unsigned int v(x);
static const int table[32] = { 0, 9, 1, 10, 13, 21, 2, 29, 11, 14, 16, 18, 22, 25, 3, 30, 8, 12, 20, 28, 15, 17, 24, 7, 19, 27, 23, 6, 26, 5, 4, 31 };
v |= v >> 1;
v |= v >> 2;
v |= v >> 4;
v |= v >> 8;
v |= v >> 16;
return table[(v * 0x07C4ACDDU) >> 27];
}
/** \returns the square root of \a x.
*
* It is essentially equivalent to \code using std::sqrt; return sqrt(x); \endcode,
* but slightly faster for float/double and some compilers (e.g., gcc), thanks to
* specializations when SSE is enabled.
*
* It's usage is justified in performance critical functions, like norm/normalize.
*/
template<typename T>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
T sqrt(const T &x)
{
EIGEN_USING_STD_MATH(sqrt);
return sqrt(x);
}
template<typename T>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
T log(const T &x) {
EIGEN_USING_STD_MATH(log);
return log(x);
}
#ifdef __CUDACC__
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
float log(const float &x) { return ::logf(x); }
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
double log(const double &x) { return ::log(x); }
#endif
template<typename T>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
typename internal::enable_if<NumTraits<T>::IsSigned || NumTraits<T>::IsComplex,typename NumTraits<T>::Real>::type
abs(const T &x) {
EIGEN_USING_STD_MATH(abs);
return abs(x);
}
template<typename T>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
typename internal::enable_if<!(NumTraits<T>::IsSigned || NumTraits<T>::IsComplex),typename NumTraits<T>::Real>::type
abs(const T &x) {
return x;
}
#if defined(__SYCL_DEVICE_ONLY__)
EIGEN_ALWAYS_INLINE float abs(float x) { return cl::sycl::fabs(x); }
EIGEN_ALWAYS_INLINE double abs(double x) { return cl::sycl::fabs(x); }
#endif // defined(__SYCL_DEVICE_ONLY__)
#ifdef __CUDACC__
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
float abs(const float &x) { return ::fabsf(x); }
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
double abs(const double &x) { return ::fabs(x); }
template <> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
float abs(const std::complex<float>& x) {
return ::hypotf(x.real(), x.imag());
}
template <> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
double abs(const std::complex<double>& x) {
return ::hypot(x.real(), x.imag());
}
#endif
template<typename T>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
T exp(const T &x) {
EIGEN_USING_STD_MATH(exp);
return exp(x);
}
#ifdef __CUDACC__
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
float exp(const float &x) { return ::expf(x); }
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
double exp(const double &x) { return ::exp(x); }
#endif
template<typename T>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
T cos(const T &x) {
EIGEN_USING_STD_MATH(cos);
return cos(x);
}
#ifdef __CUDACC__
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
float cos(const float &x) { return ::cosf(x); }
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
double cos(const double &x) { return ::cos(x); }
#endif
template<typename T>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
T sin(const T &x) {
EIGEN_USING_STD_MATH(sin);
return sin(x);
}
#ifdef __CUDACC__
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
float sin(const float &x) { return ::sinf(x); }
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
double sin(const double &x) { return ::sin(x); }
#endif
template<typename T>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
T tan(const T &x) {
EIGEN_USING_STD_MATH(tan);
return tan(x);
}
#ifdef __CUDACC__
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
float tan(const float &x) { return ::tanf(x); }
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
double tan(const double &x) { return ::tan(x); }
#endif
template<typename T>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
T acos(const T &x) {
EIGEN_USING_STD_MATH(acos);
return acos(x);
}
#ifdef __CUDACC__
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
float acos(const float &x) { return ::acosf(x); }
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
double acos(const double &x) { return ::acos(x); }
#endif
template<typename T>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
T asin(const T &x) {
EIGEN_USING_STD_MATH(asin);
return asin(x);
}
#ifdef __CUDACC__
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
float asin(const float &x) { return ::asinf(x); }
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
double asin(const double &x) { return ::asin(x); }
#endif
template<typename T>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
T atan(const T &x) {
EIGEN_USING_STD_MATH(atan);
return atan(x);
}
#ifdef __CUDACC__
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
float atan(const float &x) { return ::atanf(x); }
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
double atan(const double &x) { return ::atan(x); }
#endif
template<typename T>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
T cosh(const T &x) {
EIGEN_USING_STD_MATH(cosh);
return cosh(x);
}
#ifdef __CUDACC__
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
float cosh(const float &x) { return ::coshf(x); }
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
double cosh(const double &x) { return ::cosh(x); }
#endif
template<typename T>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
T sinh(const T &x) {
EIGEN_USING_STD_MATH(sinh);
return sinh(x);
}
#ifdef __CUDACC__
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
float sinh(const float &x) { return ::sinhf(x); }
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
double sinh(const double &x) { return ::sinh(x); }
#endif
template<typename T>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
T tanh(const T &x) {
EIGEN_USING_STD_MATH(tanh);
return tanh(x);
}
#if (!defined(__CUDACC__)) && EIGEN_FAST_MATH
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
float tanh(float x) { return internal::generic_fast_tanh_float(x); }
#endif
#ifdef __CUDACC__
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
float tanh(const float &x) { return ::tanhf(x); }
template<> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
double tanh(const double &x) { return ::tanh(x); }
#endif
template <typename T>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
T fmod(const T& a, const T& b) {
EIGEN_USING_STD_MATH(fmod);
return fmod(a, b);
}
#ifdef __CUDACC__
template <>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
float fmod(const float& a, const float& b) {
return ::fmodf(a, b);
}
template <>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
double fmod(const double& a, const double& b) {
return ::fmod(a, b);
}
#endif
} // end namespace numext
namespace internal {
template<typename T>
EIGEN_DEVICE_FUNC bool isfinite_impl(const std::complex<T>& x)
{
return (numext::isfinite)(numext::real(x)) && (numext::isfinite)(numext::imag(x));
}
template<typename T>
EIGEN_DEVICE_FUNC bool isnan_impl(const std::complex<T>& x)
{
return (numext::isnan)(numext::real(x)) || (numext::isnan)(numext::imag(x));
}
template<typename T>
EIGEN_DEVICE_FUNC bool isinf_impl(const std::complex<T>& x)
{
return ((numext::isinf)(numext::real(x)) || (numext::isinf)(numext::imag(x))) && (!(numext::isnan)(x));
}
/****************************************************************************
* Implementation of fuzzy comparisons *
****************************************************************************/
template<typename Scalar,
bool IsComplex,
bool IsInteger>
struct scalar_fuzzy_default_impl {};
template<typename Scalar>
struct scalar_fuzzy_default_impl<Scalar, false, false>
{
typedef typename NumTraits<Scalar>::Real RealScalar;
template<typename OtherScalar> EIGEN_DEVICE_FUNC
static inline bool isMuchSmallerThan(const Scalar& x, const OtherScalar& y, const RealScalar& prec)
{
return numext::abs(x) <= numext::abs(y) * prec;
}
EIGEN_DEVICE_FUNC
static inline bool isApprox(const Scalar& x, const Scalar& y, const RealScalar& prec)
{
return numext::abs(x - y) <= numext::mini(numext::abs(x), numext::abs(y)) * prec;
}
EIGEN_DEVICE_FUNC
static inline bool isApproxOrLessThan(const Scalar& x, const Scalar& y, const RealScalar& prec)
{
return x <= y || isApprox(x, y, prec);
}
};
template<typename Scalar>
struct scalar_fuzzy_default_impl<Scalar, false, true>
{
typedef typename NumTraits<Scalar>::Real RealScalar;
template<typename OtherScalar> EIGEN_DEVICE_FUNC
static inline bool isMuchSmallerThan(const Scalar& x, const Scalar&, const RealScalar&)
{
return x == Scalar(0);
}
EIGEN_DEVICE_FUNC
static inline bool isApprox(const Scalar& x, const Scalar& y, const RealScalar&)
{
return x == y;
}
EIGEN_DEVICE_FUNC
static inline bool isApproxOrLessThan(const Scalar& x, const Scalar& y, const RealScalar&)
{
return x <= y;
}
};
template<typename Scalar>
struct scalar_fuzzy_default_impl<Scalar, true, false>
{
typedef typename NumTraits<Scalar>::Real RealScalar;
template<typename OtherScalar> EIGEN_DEVICE_FUNC
static inline bool isMuchSmallerThan(const Scalar& x, const OtherScalar& y, const RealScalar& prec)
{
return numext::abs2(x) <= numext::abs2(y) * prec * prec;
}
EIGEN_DEVICE_FUNC
static inline bool isApprox(const Scalar& x, const Scalar& y, const RealScalar& prec)
{
return numext::abs2(x - y) <= numext::mini(numext::abs2(x), numext::abs2(y)) * prec * prec;
}
};
template<typename Scalar>
struct scalar_fuzzy_impl : scalar_fuzzy_default_impl<Scalar, NumTraits<Scalar>::IsComplex, NumTraits<Scalar>::IsInteger> {};
template<typename Scalar, typename OtherScalar> EIGEN_DEVICE_FUNC
inline bool isMuchSmallerThan(const Scalar& x, const OtherScalar& y,
const typename NumTraits<Scalar>::Real &precision = NumTraits<Scalar>::dummy_precision())
{
return scalar_fuzzy_impl<Scalar>::template isMuchSmallerThan<OtherScalar>(x, y, precision);
}
template<typename Scalar> EIGEN_DEVICE_FUNC
inline bool isApprox(const Scalar& x, const Scalar& y,
const typename NumTraits<Scalar>::Real &precision = NumTraits<Scalar>::dummy_precision())
{
return scalar_fuzzy_impl<Scalar>::isApprox(x, y, precision);
}
template<typename Scalar> EIGEN_DEVICE_FUNC
inline bool isApproxOrLessThan(const Scalar& x, const Scalar& y,
const typename NumTraits<Scalar>::Real &precision = NumTraits<Scalar>::dummy_precision())
{
return scalar_fuzzy_impl<Scalar>::isApproxOrLessThan(x, y, precision);
}
/******************************************
*** The special case of the bool type ***
******************************************/
template<> struct random_impl<bool>
{
static inline bool run()
{
return random<int>(0,1)==0 ? false : true;
}
};
template<> struct scalar_fuzzy_impl<bool>
{
typedef bool RealScalar;
template<typename OtherScalar> EIGEN_DEVICE_FUNC
static inline bool isMuchSmallerThan(const bool& x, const bool&, const bool&)
{
return !x;
}
EIGEN_DEVICE_FUNC
static inline bool isApprox(bool x, bool y, bool)
{
return x == y;
}
EIGEN_DEVICE_FUNC
static inline bool isApproxOrLessThan(const bool& x, const bool& y, const bool&)
{
return (!x) || y;
}
};
} // end namespace internal
} // end namespace Eigen
#endif // EIGEN_MATHFUNCTIONS_H
external/eigen3/Eigen/src/Core/MathFunctionsImpl.h 0 → 100644
View file @ a394b22a
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2014 Pedro Gonnet (pedro.gonnet@gmail.com)
// Copyright (C) 2016 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_MATHFUNCTIONSIMPL_H
#define EIGEN_MATHFUNCTIONSIMPL_H
namespace Eigen {
namespace internal {
/** \internal \returns the hyperbolic tan of \a a (coeff-wise)
Doesn't do anything fancy, just a 13/6-degree rational interpolant which
is accurate up to a couple of ulp in the range [-9, 9], outside of which
the tanh(x) = +/-1.
This implementation works on both scalars and packets.
*/
template<typename T>
T generic_fast_tanh_float(const T& a_x)
{
// Clamp the inputs to the range [-9, 9] since anything outside
// this range is +/-1.0f in single-precision.
const T plus_9 = pset1<T>(9.f);
const T minus_9 = pset1<T>(-9.f);
// NOTE GCC prior to 6.3 might improperly optimize this max/min
// step such that if a_x is nan, x will be either 9 or -9,
// and tanh will return 1 or -1 instead of nan.
// This is supposed to be fixed in gcc6.3,
// see: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=72867
const T x = pmax(minus_9,pmin(plus_9,a_x));
// The monomial coefficients of the numerator polynomial (odd).
const T alpha_1 = pset1<T>(4.89352455891786e-03f);
const T alpha_3 = pset1<T>(6.37261928875436e-04f);
const T alpha_5 = pset1<T>(1.48572235717979e-05f);
const T alpha_7 = pset1<T>(5.12229709037114e-08f);
const T alpha_9 = pset1<T>(-8.60467152213735e-11f);
const T alpha_11 = pset1<T>(2.00018790482477e-13f);
const T alpha_13 = pset1<T>(-2.76076847742355e-16f);
// The monomial coefficients of the denominator polynomial (even).
const T beta_0 = pset1<T>(4.89352518554385e-03f);
const T beta_2 = pset1<T>(2.26843463243900e-03f);
const T beta_4 = pset1<T>(1.18534705686654e-04f);
const T beta_6 = pset1<T>(1.19825839466702e-06f);
// Since the polynomials are odd/even, we need x^2.
const T x2 = pmul(x, x);
// Evaluate the numerator polynomial p.
T p = pmadd(x2, alpha_13, alpha_11);
p = pmadd(x2, p, alpha_9);
p = pmadd(x2, p, alpha_7);
p = pmadd(x2, p, alpha_5);
p = pmadd(x2, p, alpha_3);
p = pmadd(x2, p, alpha_1);
p = pmul(x, p);
// Evaluate the denominator polynomial p.
T q = pmadd(x2, beta_6, beta_4);
q = pmadd(x2, q, beta_2);
q = pmadd(x2, q, beta_0);
// Divide the numerator by the denominator.
return pdiv(p, q);
}
} // end namespace internal
} // end namespace Eigen
#endif // EIGEN_MATHFUNCTIONSIMPL_H
external/eigen3/Eigen/src/Core/Matrix.h 0 → 100644
View file @ a394b22a
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2006-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_MATRIX_H
#define EIGEN_MATRIX_H
namespace Eigen {
namespace internal {
template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
struct traits<Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> >
{
private:
enum { size = internal::size_at_compile_time<_Rows,_Cols>::ret };
typedef typename find_best_packet<_Scalar,size>::type PacketScalar;
enum {
row_major_bit = _Options&RowMajor ? RowMajorBit : 0,
is_dynamic_size_storage = _MaxRows==Dynamic || _MaxCols==Dynamic,
max_size = is_dynamic_size_storage ? Dynamic : _MaxRows*_MaxCols,
default_alignment = compute_default_alignment<_Scalar,max_size>::value,
actual_alignment = ((_Options&DontAlign)==0) ? default_alignment : 0,
required_alignment = unpacket_traits<PacketScalar>::alignment,
packet_access_bit = (packet_traits<_Scalar>::Vectorizable && (EIGEN_UNALIGNED_VECTORIZE || (actual_alignment>=required_alignment))) ? PacketAccessBit : 0
};
public:
typedef _Scalar Scalar;
typedef Dense StorageKind;
typedef Eigen::Index StorageIndex;
typedef MatrixXpr XprKind;
enum {
RowsAtCompileTime = _Rows,
ColsAtCompileTime = _Cols,
MaxRowsAtCompileTime = _MaxRows,
MaxColsAtCompileTime = _MaxCols,
Flags = compute_matrix_flags<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::ret,
Options = _Options,
InnerStrideAtCompileTime = 1,
OuterStrideAtCompileTime = (Options&RowMajor) ? ColsAtCompileTime : RowsAtCompileTime,
// FIXME, the following flag in only used to define NeedsToAlign in PlainObjectBase
EvaluatorFlags = LinearAccessBit | DirectAccessBit | packet_access_bit | row_major_bit,
Alignment = actual_alignment
};
};
}
/** \class Matrix
* \ingroup Core_Module
*
* \brief The matrix class, also used for vectors and row-vectors
*
* The %Matrix class is the work-horse for all \em dense (\ref dense "note") matrices and vectors within Eigen.
* Vectors are matrices with one column, and row-vectors are matrices with one row.
*
* The %Matrix class encompasses \em both fixed-size and dynamic-size objects (\ref fixedsize "note").
*
* The first three template parameters are required:
* \tparam _Scalar Numeric type, e.g. float, double, int or std::complex<float>.
* User defined scalar types are supported as well (see \ref user_defined_scalars "here").
* \tparam _Rows Number of rows, or \b Dynamic
* \tparam _Cols Number of columns, or \b Dynamic
*
* The remaining template parameters are optional -- in most cases you don't have to worry about them.
* \tparam _Options A combination of either \b #RowMajor or \b #ColMajor, and of either
* \b #AutoAlign or \b #DontAlign.
* The former controls \ref TopicStorageOrders "storage order", and defaults to column-major. The latter controls alignment, which is required
* for vectorization. It defaults to aligning matrices except for fixed sizes that aren't a multiple of the packet size.
* \tparam _MaxRows Maximum number of rows. Defaults to \a _Rows (\ref maxrows "note").
* \tparam _MaxCols Maximum number of columns. Defaults to \a _Cols (\ref maxrows "note").
*
* Eigen provides a number of typedefs covering the usual cases. Here are some examples:
*
* \li \c Matrix2d is a 2x2 square matrix of doubles (\c Matrix<double, 2, 2>)
* \li \c Vector4f is a vector of 4 floats (\c Matrix<float, 4, 1>)
* \li \c RowVector3i is a row-vector of 3 ints (\c Matrix<int, 1, 3>)
*
* \li \c MatrixXf is a dynamic-size matrix of floats (\c Matrix<float, Dynamic, Dynamic>)
* \li \c VectorXf is a dynamic-size vector of floats (\c Matrix<float, Dynamic, 1>)
*
* \li \c Matrix2Xf is a partially fixed-size (dynamic-size) matrix of floats (\c Matrix<float, 2, Dynamic>)
* \li \c MatrixX3d is a partially dynamic-size (fixed-size) matrix of double (\c Matrix<double, Dynamic, 3>)
*
* See \link matrixtypedefs this page \endlink for a complete list of predefined \em %Matrix and \em Vector typedefs.
*
* You can access elements of vectors and matrices using normal subscripting:
*
* \code
* Eigen::VectorXd v(10);
* v[0] = 0.1;
* v[1] = 0.2;
* v(0) = 0.3;
* v(1) = 0.4;
*
* Eigen::MatrixXi m(10, 10);
* m(0, 1) = 1;
* m(0, 2) = 2;
* m(0, 3) = 3;
* \endcode
*
* This class can be extended with the help of the plugin mechanism described on the page
* \ref TopicCustomizing_Plugins by defining the preprocessor symbol \c EIGEN_MATRIX_PLUGIN.
*
* <i><b>Some notes:</b></i>
*
* <dl>
* <dt><b>\anchor dense Dense versus sparse:</b></dt>
* <dd>This %Matrix class handles dense, not sparse matrices and vectors. For sparse matrices and vectors, see the Sparse module.
*
* Dense matrices and vectors are plain usual arrays of coefficients. All the coefficients are stored, in an ordinary contiguous array.
* This is unlike Sparse matrices and vectors where the coefficients are stored as a list of nonzero coefficients.</dd>
*
* <dt><b>\anchor fixedsize Fixed-size versus dynamic-size:</b></dt>
* <dd>Fixed-size means that the numbers of rows and columns are known are compile-time. In this case, Eigen allocates the array
* of coefficients as a fixed-size array, as a class member. This makes sense for very small matrices, typically up to 4x4, sometimes up
* to 16x16. Larger matrices should be declared as dynamic-size even if one happens to know their size at compile-time.
*
* Dynamic-size means that the numbers of rows or columns are not necessarily known at compile-time. In this case they are runtime
* variables, and the array of coefficients is allocated dynamically on the heap.
*
* Note that \em dense matrices, be they Fixed-size or Dynamic-size, <em>do not</em> expand dynamically in the sense of a std::map.
* If you want this behavior, see the Sparse module.</dd>
*
* <dt><b>\anchor maxrows _MaxRows and _MaxCols:</b></dt>
* <dd>In most cases, one just leaves these parameters to the default values.
* These parameters mean the maximum size of rows and columns that the matrix may have. They are useful in cases
* when the exact numbers of rows and columns are not known are compile-time, but it is known at compile-time that they cannot
* exceed a certain value. This happens when taking dynamic-size blocks inside fixed-size matrices: in this case _MaxRows and _MaxCols
* are the dimensions of the original matrix, while _Rows and _Cols are Dynamic.</dd>
* </dl>
*
* <i><b>ABI and storage layout</b></i>
*
* The table below summarizes the ABI of some possible Matrix instances which is fixed thorough the lifetime of Eigen 3.
* <table class="manual">
* <tr><th>Matrix type</th><th>Equivalent C structure</th></tr>
* <tr><td>\code Matrix<T,Dynamic,Dynamic> \endcode</td><td>\code
* struct {
* T *data; // with (size_t(data)%EIGEN_MAX_ALIGN_BYTES)==0
* Eigen::Index rows, cols;
* };
* \endcode</td></tr>
* <tr class="alt"><td>\code
* Matrix<T,Dynamic,1>
* Matrix<T,1,Dynamic> \endcode</td><td>\code
* struct {
* T *data; // with (size_t(data)%EIGEN_MAX_ALIGN_BYTES)==0
* Eigen::Index size;
* };
* \endcode</td></tr>
* <tr><td>\code Matrix<T,Rows,Cols> \endcode</td><td>\code
* struct {
* T data[Rows*Cols]; // with (size_t(data)%A(Rows*Cols*sizeof(T)))==0
* };
* \endcode</td></tr>
* <tr class="alt"><td>\code Matrix<T,Dynamic,Dynamic,0,MaxRows,MaxCols> \endcode</td><td>\code
* struct {
* T data[MaxRows*MaxCols]; // with (size_t(data)%A(MaxRows*MaxCols*sizeof(T)))==0
* Eigen::Index rows, cols;
* };
* \endcode</td></tr>
* </table>
* Note that in this table Rows, Cols, MaxRows and MaxCols are all positive integers. A(S) is defined to the largest possible power-of-two
* smaller to EIGEN_MAX_STATIC_ALIGN_BYTES.
*
* \see MatrixBase for the majority of the API methods for matrices, \ref TopicClassHierarchy,
* \ref TopicStorageOrders
*/
template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
class Matrix
: public PlainObjectBase<Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> >
{
public:
/** \brief Base class typedef.
* \sa PlainObjectBase
*/
typedef PlainObjectBase<Matrix> Base;
enum { Options = _Options };
EIGEN_DENSE_PUBLIC_INTERFACE(Matrix)
typedef typename Base::PlainObject PlainObject;
using Base::base;
using Base::coeffRef;
/**
* \brief Assigns matrices to each other.
*
* \note This is a special case of the templated operator=. Its purpose is
* to prevent a default operator= from hiding the templated operator=.
*
* \callgraph
*/
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE Matrix& operator=(const Matrix& other)
{
return Base::_set(other);
}
/** \internal
* \brief Copies the value of the expression \a other into \c *this with automatic resizing.
*
* *this might be resized to match the dimensions of \a other. If *this was a null matrix (not already initialized),
* it will be initialized.
*
* Note that copying a row-vector into a vector (and conversely) is allowed.
* The resizing, if any, is then done in the appropriate way so that row-vectors
* remain row-vectors and vectors remain vectors.
*/
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE Matrix& operator=(const DenseBase<OtherDerived>& other)
{
return Base::_set(other);
}
/* Here, doxygen failed to copy the brief information when using \copydoc */
/**
* \brief Copies the generic expression \a other into *this.
* \copydetails DenseBase::operator=(const EigenBase<OtherDerived> &other)
*/
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE Matrix& operator=(const EigenBase<OtherDerived> &other)
{
return Base::operator=(other);
}
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE Matrix& operator=(const ReturnByValue<OtherDerived>& func)
{
return Base::operator=(func);
}
/** \brief Default constructor.
*
* For fixed-size matrices, does nothing.
*
* For dynamic-size matrices, creates an empty matrix of size 0. Does not allocate any array. Such a matrix
* is called a null matrix. This constructor is the unique way to create null matrices: resizing
* a matrix to 0 is not supported.
*
* \sa resize(Index,Index)
*/
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE Matrix() : Base()
{
Base::_check_template_params();
EIGEN_INITIALIZE_COEFFS_IF_THAT_OPTION_IS_ENABLED
}
// FIXME is it still needed
EIGEN_DEVICE_FUNC
explicit Matrix(internal::constructor_without_unaligned_array_assert)
: Base(internal::constructor_without_unaligned_array_assert())
{ Base::_check_template_params(); EIGEN_INITIALIZE_COEFFS_IF_THAT_OPTION_IS_ENABLED }
#if EIGEN_HAS_RVALUE_REFERENCES
EIGEN_DEVICE_FUNC
Matrix(Matrix&& other) EIGEN_NOEXCEPT_IF(std::is_nothrow_move_constructible<Scalar>::value)
: Base(std::move(other))
{
Base::_check_template_params();
if (RowsAtCompileTime!=Dynamic && ColsAtCompileTime!=Dynamic)
Base::_set_noalias(other);
}
EIGEN_DEVICE_FUNC
Matrix& operator=(Matrix&& other) EIGEN_NOEXCEPT_IF(std::is_nothrow_move_assignable<Scalar>::value)
{
other.swap(*this);
return *this;
}
#endif
#ifndef EIGEN_PARSED_BY_DOXYGEN
// This constructor is for both 1x1 matrices and dynamic vectors
template<typename T>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE explicit Matrix(const T& x)
{
Base::_check_template_params();
Base::template _init1<T>(x);
}
template<typename T0, typename T1>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE Matrix(const T0& x, const T1& y)
{
Base::_check_template_params();
Base::template _init2<T0,T1>(x, y);
}
#else
/** \brief Constructs a fixed-sized matrix initialized with coefficients starting at \a data */
EIGEN_DEVICE_FUNC
explicit Matrix(const Scalar *data);
/** \brief Constructs a vector or row-vector with given dimension. \only_for_vectors
*
* This is useful for dynamic-size vectors. For fixed-size vectors,
* it is redundant to pass these parameters, so one should use the default constructor
* Matrix() instead.
*
* \warning This constructor is disabled for fixed-size \c 1x1 matrices. For instance,
* calling Matrix<double,1,1>(1) will call the initialization constructor: Matrix(const Scalar&).
* For fixed-size \c 1x1 matrices it is therefore recommended to use the default
* constructor Matrix() instead, especially when using one of the non standard
* \c EIGEN_INITIALIZE_MATRICES_BY_{ZERO,\c NAN} macros (see \ref TopicPreprocessorDirectives).
*/
EIGEN_STRONG_INLINE explicit Matrix(Index dim);
/** \brief Constructs an initialized 1x1 matrix with the given coefficient */
Matrix(const Scalar& x);
/** \brief Constructs an uninitialized matrix with \a rows rows and \a cols columns.
*
* This is useful for dynamic-size matrices. For fixed-size matrices,
* it is redundant to pass these parameters, so one should use the default constructor
* Matrix() instead.
*
* \warning This constructor is disabled for fixed-size \c 1x2 and \c 2x1 vectors. For instance,
* calling Matrix2f(2,1) will call the initialization constructor: Matrix(const Scalar& x, const Scalar& y).
* For fixed-size \c 1x2 or \c 2x1 vectors it is therefore recommended to use the default
* constructor Matrix() instead, especially when using one of the non standard
* \c EIGEN_INITIALIZE_MATRICES_BY_{ZERO,\c NAN} macros (see \ref TopicPreprocessorDirectives).
*/
EIGEN_DEVICE_FUNC
Matrix(Index rows, Index cols);
/** \brief Constructs an initialized 2D vector with given coefficients */
Matrix(const Scalar& x, const Scalar& y);
#endif
/** \brief Constructs an initialized 3D vector with given coefficients */
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE Matrix(const Scalar& x, const Scalar& y, const Scalar& z)
{
Base::_check_template_params();
EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Matrix, 3)
m_storage.data()[0] = x;
m_storage.data()[1] = y;
m_storage.data()[2] = z;
}
/** \brief Constructs an initialized 4D vector with given coefficients */
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE Matrix(const Scalar& x, const Scalar& y, const Scalar& z, const Scalar& w)
{
Base::_check_template_params();
EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Matrix, 4)
m_storage.data()[0] = x;
m_storage.data()[1] = y;
m_storage.data()[2] = z;
m_storage.data()[3] = w;
}
/** \brief Copy constructor */
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE Matrix(const Matrix& other) : Base(other)
{ }
/** \brief Copy constructor for generic expressions.
* \sa MatrixBase::operator=(const EigenBase<OtherDerived>&)
*/
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE Matrix(const EigenBase<OtherDerived> &other)
: Base(other.derived())
{ }
EIGEN_DEVICE_FUNC inline Index innerStride() const { return 1; }
EIGEN_DEVICE_FUNC inline Index outerStride() const { return this->innerSize(); }
/////////// Geometry module ///////////
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
explicit Matrix(const RotationBase<OtherDerived,ColsAtCompileTime>& r);
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
Matrix& operator=(const RotationBase<OtherDerived,ColsAtCompileTime>& r);
// allow to extend Matrix outside Eigen
#ifdef EIGEN_MATRIX_PLUGIN
#include EIGEN_MATRIX_PLUGIN
#endif
protected:
template <typename Derived, typename OtherDerived, bool IsVector>
friend struct internal::conservative_resize_like_impl;
using Base::m_storage;
};
/** \defgroup matrixtypedefs Global matrix typedefs
*
* \ingroup Core_Module
*
* Eigen defines several typedef shortcuts for most common matrix and vector types.
*
* The general patterns are the following:
*
* \c MatrixSizeType where \c Size can be \c 2,\c 3,\c 4 for fixed size square matrices or \c X for dynamic size,
* and where \c Type can be \c i for integer, \c f for float, \c d for double, \c cf for complex float, \c cd
* for complex double.
*
* For example, \c Matrix3d is a fixed-size 3x3 matrix type of doubles, and \c MatrixXf is a dynamic-size matrix of floats.
*
* There are also \c VectorSizeType and \c RowVectorSizeType which are self-explanatory. For example, \c Vector4cf is
* a fixed-size vector of 4 complex floats.
*
* \sa class Matrix
*/
#define EIGEN_MAKE_TYPEDEFS(Type, TypeSuffix, Size, SizeSuffix) \
/** \ingroup matrixtypedefs */ \
typedef Matrix<Type, Size, Size> Matrix##SizeSuffix##TypeSuffix; \
/** \ingroup matrixtypedefs */ \
typedef Matrix<Type, Size, 1> Vector##SizeSuffix##TypeSuffix; \
/** \ingroup matrixtypedefs */ \
typedef Matrix<Type, 1, Size> RowVector##SizeSuffix##TypeSuffix;
#define EIGEN_MAKE_FIXED_TYPEDEFS(Type, TypeSuffix, Size) \
/** \ingroup matrixtypedefs */ \
typedef Matrix<Type, Size, Dynamic> Matrix##Size##X##TypeSuffix; \
/** \ingroup matrixtypedefs */ \
typedef Matrix<Type, Dynamic, Size> Matrix##X##Size##TypeSuffix;
#define EIGEN_MAKE_TYPEDEFS_ALL_SIZES(Type, TypeSuffix) \
EIGEN_MAKE_TYPEDEFS(Type, TypeSuffix, 2, 2) \
EIGEN_MAKE_TYPEDEFS(Type, TypeSuffix, 3, 3) \
EIGEN_MAKE_TYPEDEFS(Type, TypeSuffix, 4, 4) \
EIGEN_MAKE_TYPEDEFS(Type, TypeSuffix, Dynamic, X) \
EIGEN_MAKE_FIXED_TYPEDEFS(Type, TypeSuffix, 2) \
EIGEN_MAKE_FIXED_TYPEDEFS(Type, TypeSuffix, 3) \
EIGEN_MAKE_FIXED_TYPEDEFS(Type, TypeSuffix, 4)
EIGEN_MAKE_TYPEDEFS_ALL_SIZES(int, i)
EIGEN_MAKE_TYPEDEFS_ALL_SIZES(float, f)
EIGEN_MAKE_TYPEDEFS_ALL_SIZES(double, d)
EIGEN_MAKE_TYPEDEFS_ALL_SIZES(std::complex<float>, cf)
EIGEN_MAKE_TYPEDEFS_ALL_SIZES(std::complex<double>, cd)
#undef EIGEN_MAKE_TYPEDEFS_ALL_SIZES
#undef EIGEN_MAKE_TYPEDEFS
#undef EIGEN_MAKE_FIXED_TYPEDEFS
} // end namespace Eigen
#endif // EIGEN_MATRIX_H
external/eigen3/Eigen/src/Core/MatrixBase.h 0 → 100644
View file @ a394b22a
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2006-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_MATRIXBASE_H
#define EIGEN_MATRIXBASE_H
namespace Eigen {
/** \class MatrixBase
* \ingroup Core_Module
*
* \brief Base class for all dense matrices, vectors, and expressions
*
* This class is the base that is inherited by all matrix, vector, and related expression
* types. Most of the Eigen API is contained in this class, and its base classes. Other important
* classes for the Eigen API are Matrix, and VectorwiseOp.
*
* Note that some methods are defined in other modules such as the \ref LU_Module LU module
* for all functions related to matrix inversions.
*
* \tparam Derived is the derived type, e.g. a matrix type, or an expression, etc.
*
* When writing a function taking Eigen objects as argument, if you want your function
* to take as argument any matrix, vector, or expression, just let it take a
* MatrixBase argument. As an example, here is a function printFirstRow which, given
* a matrix, vector, or expression \a x, prints the first row of \a x.
*
* \code
template<typename Derived>
void printFirstRow(const Eigen::MatrixBase<Derived>& x)
{
cout << x.row(0) << endl;
}
* \endcode
*
* This class can be extended with the help of the plugin mechanism described on the page
* \ref TopicCustomizing_Plugins by defining the preprocessor symbol \c EIGEN_MATRIXBASE_PLUGIN.
*
* \sa \blank \ref TopicClassHierarchy
*/
template<typename Derived> class MatrixBase
: public DenseBase<Derived>
{
public:
#ifndef EIGEN_PARSED_BY_DOXYGEN
typedef MatrixBase StorageBaseType;
typedef typename internal::traits<Derived>::StorageKind StorageKind;
typedef typename internal::traits<Derived>::StorageIndex StorageIndex;
typedef typename internal::traits<Derived>::Scalar Scalar;
typedef typename internal::packet_traits<Scalar>::type PacketScalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef DenseBase<Derived> Base;
using Base::RowsAtCompileTime;
using Base::ColsAtCompileTime;
using Base::SizeAtCompileTime;
using Base::MaxRowsAtCompileTime;
using Base::MaxColsAtCompileTime;
using Base::MaxSizeAtCompileTime;
using Base::IsVectorAtCompileTime;
using Base::Flags;
using Base::derived;
using Base::const_cast_derived;
using Base::rows;
using Base::cols;
using Base::size;
using Base::coeff;
using Base::coeffRef;
using Base::lazyAssign;
using Base::eval;
using Base::operator+=;
using Base::operator-=;
using Base::operator*=;
using Base::operator/=;
typedef typename Base::CoeffReturnType CoeffReturnType;
typedef typename Base::ConstTransposeReturnType ConstTransposeReturnType;
typedef typename Base::RowXpr RowXpr;
typedef typename Base::ColXpr ColXpr;
#endif // not EIGEN_PARSED_BY_DOXYGEN
#ifndef EIGEN_PARSED_BY_DOXYGEN
/** type of the equivalent square matrix */
typedef Matrix<Scalar,EIGEN_SIZE_MAX(RowsAtCompileTime,ColsAtCompileTime),
EIGEN_SIZE_MAX(RowsAtCompileTime,ColsAtCompileTime)> SquareMatrixType;
#endif // not EIGEN_PARSED_BY_DOXYGEN
/** \returns the size of the main diagonal, which is min(rows(),cols()).
* \sa rows(), cols(), SizeAtCompileTime. */
EIGEN_DEVICE_FUNC
inline Index diagonalSize() const { return (numext::mini)(rows(),cols()); }
typedef typename Base::PlainObject PlainObject;
#ifndef EIGEN_PARSED_BY_DOXYGEN
/** \internal Represents a matrix with all coefficients equal to one another*/
typedef CwiseNullaryOp<internal::scalar_constant_op<Scalar>,PlainObject> ConstantReturnType;
/** \internal the return type of MatrixBase::adjoint() */
typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
CwiseUnaryOp<internal::scalar_conjugate_op<Scalar>, ConstTransposeReturnType>,
ConstTransposeReturnType
>::type AdjointReturnType;
/** \internal Return type of eigenvalues() */
typedef Matrix<std::complex<RealScalar>, internal::traits<Derived>::ColsAtCompileTime, 1, ColMajor> EigenvaluesReturnType;
/** \internal the return type of identity */
typedef CwiseNullaryOp<internal::scalar_identity_op<Scalar>,PlainObject> IdentityReturnType;
/** \internal the return type of unit vectors */
typedef Block<const CwiseNullaryOp<internal::scalar_identity_op<Scalar>, SquareMatrixType>,
internal::traits<Derived>::RowsAtCompileTime,
internal::traits<Derived>::ColsAtCompileTime> BasisReturnType;
#endif // not EIGEN_PARSED_BY_DOXYGEN
#define EIGEN_CURRENT_STORAGE_BASE_CLASS Eigen::MatrixBase
#define EIGEN_DOC_UNARY_ADDONS(X,Y)
# include "../plugins/CommonCwiseUnaryOps.h"
# include "../plugins/CommonCwiseBinaryOps.h"
# include "../plugins/MatrixCwiseUnaryOps.h"
# include "../plugins/MatrixCwiseBinaryOps.h"
# ifdef EIGEN_MATRIXBASE_PLUGIN
# include EIGEN_MATRIXBASE_PLUGIN
# endif
#undef EIGEN_CURRENT_STORAGE_BASE_CLASS
#undef EIGEN_DOC_UNARY_ADDONS
/** Special case of the template operator=, in order to prevent the compiler
* from generating a default operator= (issue hit with g++ 4.1)
*/
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
Derived& operator=(const MatrixBase& other);
// We cannot inherit here via Base::operator= since it is causing
// trouble with MSVC.
template <typename OtherDerived>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
Derived& operator=(const DenseBase<OtherDerived>& other);
template <typename OtherDerived>
EIGEN_DEVICE_FUNC
Derived& operator=(const EigenBase<OtherDerived>& other);
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
Derived& operator=(const ReturnByValue<OtherDerived>& other);
template<typename OtherDerived>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
Derived& operator+=(const MatrixBase<OtherDerived>& other);
template<typename OtherDerived>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
Derived& operator-=(const MatrixBase<OtherDerived>& other);
#ifdef __CUDACC__
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
const Product<Derived,OtherDerived,LazyProduct>
operator*(const MatrixBase<OtherDerived> &other) const
{ return this->lazyProduct(other); }
#else
template<typename OtherDerived>
const Product<Derived,OtherDerived>
operator*(const MatrixBase<OtherDerived> &other) const;
#endif
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
const Product<Derived,OtherDerived,LazyProduct>
lazyProduct(const MatrixBase<OtherDerived> &other) const;
template<typename OtherDerived>
Derived& operator*=(const EigenBase<OtherDerived>& other);
template<typename OtherDerived>
void applyOnTheLeft(const EigenBase<OtherDerived>& other);
template<typename OtherDerived>
void applyOnTheRight(const EigenBase<OtherDerived>& other);
template<typename DiagonalDerived>
EIGEN_DEVICE_FUNC
const Product<Derived, DiagonalDerived, LazyProduct>
operator*(const DiagonalBase<DiagonalDerived> &diagonal) const;
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
typename ScalarBinaryOpTraits<typename internal::traits<Derived>::Scalar,typename internal::traits<OtherDerived>::Scalar>::ReturnType
dot(const MatrixBase<OtherDerived>& other) const;
EIGEN_DEVICE_FUNC RealScalar squaredNorm() const;
EIGEN_DEVICE_FUNC RealScalar norm() const;
RealScalar stableNorm() const;
RealScalar blueNorm() const;
RealScalar hypotNorm() const;
EIGEN_DEVICE_FUNC const PlainObject normalized() const;
EIGEN_DEVICE_FUNC const PlainObject stableNormalized() const;
EIGEN_DEVICE_FUNC void normalize();
EIGEN_DEVICE_FUNC void stableNormalize();
EIGEN_DEVICE_FUNC const AdjointReturnType adjoint() const;
EIGEN_DEVICE_FUNC void adjointInPlace();
typedef Diagonal<Derived> DiagonalReturnType;
EIGEN_DEVICE_FUNC
DiagonalReturnType diagonal();
typedef typename internal::add_const<Diagonal<const Derived> >::type ConstDiagonalReturnType;
EIGEN_DEVICE_FUNC
ConstDiagonalReturnType diagonal() const;
template<int Index> struct DiagonalIndexReturnType { typedef Diagonal<Derived,Index> Type; };
template<int Index> struct ConstDiagonalIndexReturnType { typedef const Diagonal<const Derived,Index> Type; };
template<int Index>
EIGEN_DEVICE_FUNC
typename DiagonalIndexReturnType<Index>::Type diagonal();
template<int Index>
EIGEN_DEVICE_FUNC
typename ConstDiagonalIndexReturnType<Index>::Type diagonal() const;
typedef Diagonal<Derived,DynamicIndex> DiagonalDynamicIndexReturnType;
typedef typename internal::add_const<Diagonal<const Derived,DynamicIndex> >::type ConstDiagonalDynamicIndexReturnType;
EIGEN_DEVICE_FUNC
DiagonalDynamicIndexReturnType diagonal(Index index);
EIGEN_DEVICE_FUNC
ConstDiagonalDynamicIndexReturnType diagonal(Index index) const;
template<unsigned int Mode> struct TriangularViewReturnType { typedef TriangularView<Derived, Mode> Type; };
template<unsigned int Mode> struct ConstTriangularViewReturnType { typedef const TriangularView<const Derived, Mode> Type; };
template<unsigned int Mode>
EIGEN_DEVICE_FUNC
typename TriangularViewReturnType<Mode>::Type triangularView();
template<unsigned int Mode>
EIGEN_DEVICE_FUNC
typename ConstTriangularViewReturnType<Mode>::Type triangularView() const;
template<unsigned int UpLo> struct SelfAdjointViewReturnType { typedef SelfAdjointView<Derived, UpLo> Type; };
template<unsigned int UpLo> struct ConstSelfAdjointViewReturnType { typedef const SelfAdjointView<const Derived, UpLo> Type; };
template<unsigned int UpLo>
EIGEN_DEVICE_FUNC
typename SelfAdjointViewReturnType<UpLo>::Type selfadjointView();
template<unsigned int UpLo>
EIGEN_DEVICE_FUNC
typename ConstSelfAdjointViewReturnType<UpLo>::Type selfadjointView() const;
const SparseView<Derived> sparseView(const Scalar& m_reference = Scalar(0),
const typename NumTraits<Scalar>::Real& m_epsilon = NumTraits<Scalar>::dummy_precision()) const;
EIGEN_DEVICE_FUNC static const IdentityReturnType Identity();
EIGEN_DEVICE_FUNC static const IdentityReturnType Identity(Index rows, Index cols);
EIGEN_DEVICE_FUNC static const BasisReturnType Unit(Index size, Index i);
EIGEN_DEVICE_FUNC static const BasisReturnType Unit(Index i);
EIGEN_DEVICE_FUNC static const BasisReturnType UnitX();
EIGEN_DEVICE_FUNC static const BasisReturnType UnitY();
EIGEN_DEVICE_FUNC static const BasisReturnType UnitZ();
EIGEN_DEVICE_FUNC static const BasisReturnType UnitW();
EIGEN_DEVICE_FUNC
const DiagonalWrapper<const Derived> asDiagonal() const;
const PermutationWrapper<const Derived> asPermutation() const;
EIGEN_DEVICE_FUNC
Derived& setIdentity();
EIGEN_DEVICE_FUNC
Derived& setIdentity(Index rows, Index cols);
bool isIdentity(const RealScalar& prec = NumTraits<Scalar>::dummy_precision()) const;
bool isDiagonal(const RealScalar& prec = NumTraits<Scalar>::dummy_precision()) const;
bool isUpperTriangular(const RealScalar& prec = NumTraits<Scalar>::dummy_precision()) const;
bool isLowerTriangular(const RealScalar& prec = NumTraits<Scalar>::dummy_precision()) const;
template<typename OtherDerived>
bool isOrthogonal(const MatrixBase<OtherDerived>& other,
const RealScalar& prec = NumTraits<Scalar>::dummy_precision()) const;
bool isUnitary(const RealScalar& prec = NumTraits<Scalar>::dummy_precision()) const;
/** \returns true if each coefficients of \c *this and \a other are all exactly equal.
* \warning When using floating point scalar values you probably should rather use a
* fuzzy comparison such as isApprox()
* \sa isApprox(), operator!= */
template<typename OtherDerived>
EIGEN_DEVICE_FUNC inline bool operator==(const MatrixBase<OtherDerived>& other) const
{ return cwiseEqual(other).all(); }
/** \returns true if at least one pair of coefficients of \c *this and \a other are not exactly equal to each other.
* \warning When using floating point scalar values you probably should rather use a
* fuzzy comparison such as isApprox()
* \sa isApprox(), operator== */
template<typename OtherDerived>
EIGEN_DEVICE_FUNC inline bool operator!=(const MatrixBase<OtherDerived>& other) const
{ return cwiseNotEqual(other).any(); }
NoAlias<Derived,Eigen::MatrixBase > noalias();
// TODO forceAlignedAccess is temporarily disabled
// Need to find a nicer workaround.
inline const Derived& forceAlignedAccess() const { return derived(); }
inline Derived& forceAlignedAccess() { return derived(); }
template<bool Enable> inline const Derived& forceAlignedAccessIf() const { return derived(); }
template<bool Enable> inline Derived& forceAlignedAccessIf() { return derived(); }
EIGEN_DEVICE_FUNC Scalar trace() const;
template<int p> EIGEN_DEVICE_FUNC RealScalar lpNorm() const;
EIGEN_DEVICE_FUNC MatrixBase<Derived>& matrix() { return *this; }
EIGEN_DEVICE_FUNC const MatrixBase<Derived>& matrix() const { return *this; }
/** \returns an \link Eigen::ArrayBase Array \endlink expression of this matrix
* \sa ArrayBase::matrix() */
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE ArrayWrapper<Derived> array() { return ArrayWrapper<Derived>(derived()); }
/** \returns a const \link Eigen::ArrayBase Array \endlink expression of this matrix
* \sa ArrayBase::matrix() */
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const ArrayWrapper<const Derived> array() const { return ArrayWrapper<const Derived>(derived()); }
/////////// LU module ///////////
inline const FullPivLU<PlainObject> fullPivLu() const;
inline const PartialPivLU<PlainObject> partialPivLu() const;
inline const PartialPivLU<PlainObject> lu() const;
inline const Inverse<Derived> inverse() const;
template<typename ResultType>
inline void computeInverseAndDetWithCheck(
ResultType& inverse,
typename ResultType::Scalar& determinant,
bool& invertible,
const RealScalar& absDeterminantThreshold = NumTraits<Scalar>::dummy_precision()
) const;
template<typename ResultType>
inline void computeInverseWithCheck(
ResultType& inverse,
bool& invertible,
const RealScalar& absDeterminantThreshold = NumTraits<Scalar>::dummy_precision()
) const;
Scalar determinant() const;
/////////// Cholesky module ///////////
inline const LLT<PlainObject> llt() const;
inline const LDLT<PlainObject> ldlt() const;
/////////// QR module ///////////
inline const HouseholderQR<PlainObject> householderQr() const;
inline const ColPivHouseholderQR<PlainObject> colPivHouseholderQr() const;
inline const FullPivHouseholderQR<PlainObject> fullPivHouseholderQr() const;
inline const CompleteOrthogonalDecomposition<PlainObject> completeOrthogonalDecomposition() const;
/////////// Eigenvalues module ///////////
inline EigenvaluesReturnType eigenvalues() const;
inline RealScalar operatorNorm() const;
/////////// SVD module ///////////
inline JacobiSVD<PlainObject> jacobiSvd(unsigned int computationOptions = 0) const;
inline BDCSVD<PlainObject> bdcSvd(unsigned int computationOptions = 0) const;
/////////// Geometry module ///////////
#ifndef EIGEN_PARSED_BY_DOXYGEN
/// \internal helper struct to form the return type of the cross product
template<typename OtherDerived> struct cross_product_return_type {
typedef typename ScalarBinaryOpTraits<typename internal::traits<Derived>::Scalar,typename internal::traits<OtherDerived>::Scalar>::ReturnType Scalar;
typedef Matrix<Scalar,MatrixBase::RowsAtCompileTime,MatrixBase::ColsAtCompileTime> type;
};
#endif // EIGEN_PARSED_BY_DOXYGEN
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
#ifndef EIGEN_PARSED_BY_DOXYGEN
inline typename cross_product_return_type<OtherDerived>::type
#else
inline PlainObject
#endif
cross(const MatrixBase<OtherDerived>& other) const;
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
inline PlainObject cross3(const MatrixBase<OtherDerived>& other) const;
EIGEN_DEVICE_FUNC
inline PlainObject unitOrthogonal(void) const;
EIGEN_DEVICE_FUNC
inline Matrix<Scalar,3,1> eulerAngles(Index a0, Index a1, Index a2) const;
// put this as separate enum value to work around possible GCC 4.3 bug (?)
enum { HomogeneousReturnTypeDirection = ColsAtCompileTime==1&&RowsAtCompileTime==1 ? ((internal::traits<Derived>::Flags&RowMajorBit)==RowMajorBit ? Horizontal : Vertical)
: ColsAtCompileTime==1 ? Vertical : Horizontal };
typedef Homogeneous<Derived, HomogeneousReturnTypeDirection> HomogeneousReturnType;
EIGEN_DEVICE_FUNC
inline HomogeneousReturnType homogeneous() const;
enum {
SizeMinusOne = SizeAtCompileTime==Dynamic ? Dynamic : SizeAtCompileTime-1
};
typedef Block<const Derived,
internal::traits<Derived>::ColsAtCompileTime==1 ? SizeMinusOne : 1,
internal::traits<Derived>::ColsAtCompileTime==1 ? 1 : SizeMinusOne> ConstStartMinusOne;
typedef EIGEN_EXPR_BINARYOP_SCALAR_RETURN_TYPE(ConstStartMinusOne,Scalar,quotient) HNormalizedReturnType;
EIGEN_DEVICE_FUNC
inline const HNormalizedReturnType hnormalized() const;
////////// Householder module ///////////
void makeHouseholderInPlace(Scalar& tau, RealScalar& beta);
template<typename EssentialPart>
void makeHouseholder(EssentialPart& essential,
Scalar& tau, RealScalar& beta) const;
template<typename EssentialPart>
void applyHouseholderOnTheLeft(const EssentialPart& essential,
const Scalar& tau,
Scalar* workspace);
template<typename EssentialPart>
void applyHouseholderOnTheRight(const EssentialPart& essential,
const Scalar& tau,
Scalar* workspace);
///////// Jacobi module /////////
template<typename OtherScalar>
void applyOnTheLeft(Index p, Index q, const JacobiRotation<OtherScalar>& j);
template<typename OtherScalar>
void applyOnTheRight(Index p, Index q, const JacobiRotation<OtherScalar>& j);
///////// SparseCore module /////////
template<typename OtherDerived>
EIGEN_STRONG_INLINE const typename SparseMatrixBase<OtherDerived>::template CwiseProductDenseReturnType<Derived>::Type
cwiseProduct(const SparseMatrixBase<OtherDerived> &other) const
{
return other.cwiseProduct(derived());
}
///////// MatrixFunctions module /////////
typedef typename internal::stem_function<Scalar>::type StemFunction;
const MatrixExponentialReturnValue<Derived> exp() const;
const MatrixFunctionReturnValue<Derived> matrixFunction(StemFunction f) const;
const MatrixFunctionReturnValue<Derived> cosh() const;
const MatrixFunctionReturnValue<Derived> sinh() const;
const MatrixFunctionReturnValue<Derived> cos() const;
const MatrixFunctionReturnValue<Derived> sin() const;
const MatrixSquareRootReturnValue<Derived> sqrt() const;
const MatrixLogarithmReturnValue<Derived> log() const;
const MatrixPowerReturnValue<Derived> pow(const RealScalar& p) const;
const MatrixComplexPowerReturnValue<Derived> pow(const std::complex<RealScalar>& p) const;
protected:
EIGEN_DEVICE_FUNC MatrixBase() : Base() {}
private:
EIGEN_DEVICE_FUNC explicit MatrixBase(int);
EIGEN_DEVICE_FUNC MatrixBase(int,int);
template<typename OtherDerived> EIGEN_DEVICE_FUNC explicit MatrixBase(const MatrixBase<OtherDerived>&);
protected:
// mixing arrays and matrices is not legal
template<typename OtherDerived> Derived& operator+=(const ArrayBase<OtherDerived>& )
{EIGEN_STATIC_ASSERT(std::ptrdiff_t(sizeof(typename OtherDerived::Scalar))==-1,YOU_CANNOT_MIX_ARRAYS_AND_MATRICES); return *this;}
// mixing arrays and matrices is not legal
template<typename OtherDerived> Derived& operator-=(const ArrayBase<OtherDerived>& )
{EIGEN_STATIC_ASSERT(std::ptrdiff_t(sizeof(typename OtherDerived::Scalar))==-1,YOU_CANNOT_MIX_ARRAYS_AND_MATRICES); return *this;}
};
/***************************************************************************
* Implementation of matrix base methods
***************************************************************************/
/** replaces \c *this by \c *this * \a other.
*
* \returns a reference to \c *this
*
* Example: \include MatrixBase_applyOnTheRight.cpp
* Output: \verbinclude MatrixBase_applyOnTheRight.out
*/
template<typename Derived>
template<typename OtherDerived>
inline Derived&
MatrixBase<Derived>::operator*=(const EigenBase<OtherDerived> &other)
{
other.derived().applyThisOnTheRight(derived());
return derived();
}
/** replaces \c *this by \c *this * \a other. It is equivalent to MatrixBase::operator*=().
*
* Example: \include MatrixBase_applyOnTheRight.cpp
* Output: \verbinclude MatrixBase_applyOnTheRight.out
*/
template<typename Derived>
template<typename OtherDerived>
inline void MatrixBase<Derived>::applyOnTheRight(const EigenBase<OtherDerived> &other)
{
other.derived().applyThisOnTheRight(derived());
}
/** replaces \c *this by \a other * \c *this.
*
* Example: \include MatrixBase_applyOnTheLeft.cpp
* Output: \verbinclude MatrixBase_applyOnTheLeft.out
*/
template<typename Derived>
template<typename OtherDerived>
inline void MatrixBase<Derived>::applyOnTheLeft(const EigenBase<OtherDerived> &other)
{
other.derived().applyThisOnTheLeft(derived());
}
} // end namespace Eigen
#endif // EIGEN_MATRIXBASE_H
external/eigen3/Eigen/src/Core/NestByValue.h 0 → 100644
View file @ a394b22a
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_NESTBYVALUE_H
#define EIGEN_NESTBYVALUE_H
namespace Eigen {
namespace internal {
template<typename ExpressionType>
struct traits<NestByValue<ExpressionType> > : public traits<ExpressionType>
{};
}
/** \class NestByValue
* \ingroup Core_Module
*
* \brief Expression which must be nested by value
*
* \tparam ExpressionType the type of the object of which we are requiring nesting-by-value
*
* This class is the return type of MatrixBase::nestByValue()
* and most of the time this is the only way it is used.
*
* \sa MatrixBase::nestByValue()
*/
template<typename ExpressionType> class NestByValue
: public internal::dense_xpr_base< NestByValue<ExpressionType> >::type
{
public:
typedef typename internal::dense_xpr_base<NestByValue>::type Base;
EIGEN_DENSE_PUBLIC_INTERFACE(NestByValue)
EIGEN_DEVICE_FUNC explicit inline NestByValue(const ExpressionType& matrix) : m_expression(matrix) {}
EIGEN_DEVICE_FUNC inline Index rows() const { return m_expression.rows(); }
EIGEN_DEVICE_FUNC inline Index cols() const { return m_expression.cols(); }
EIGEN_DEVICE_FUNC inline Index outerStride() const { return m_expression.outerStride(); }
EIGEN_DEVICE_FUNC inline Index innerStride() const { return m_expression.innerStride(); }
EIGEN_DEVICE_FUNC inline const CoeffReturnType coeff(Index row, Index col) const
{
return m_expression.coeff(row, col);
}
EIGEN_DEVICE_FUNC inline Scalar& coeffRef(Index row, Index col)
{
return m_expression.const_cast_derived().coeffRef(row, col);
}
EIGEN_DEVICE_FUNC inline const CoeffReturnType coeff(Index index) const
{
return m_expression.coeff(index);
}
EIGEN_DEVICE_FUNC inline Scalar& coeffRef(Index index)
{
return m_expression.const_cast_derived().coeffRef(index);
}
template<int LoadMode>
inline const PacketScalar packet(Index row, Index col) const
{
return m_expression.template packet<LoadMode>(row, col);
}
template<int LoadMode>
inline void writePacket(Index row, Index col, const PacketScalar& x)
{
m_expression.const_cast_derived().template writePacket<LoadMode>(row, col, x);
}
template<int LoadMode>
inline const PacketScalar packet(Index index) const
{
return m_expression.template packet<LoadMode>(index);
}
template<int LoadMode>
inline void writePacket(Index index, const PacketScalar& x)
{
m_expression.const_cast_derived().template writePacket<LoadMode>(index, x);
}
EIGEN_DEVICE_FUNC operator const ExpressionType&() const { return m_expression; }
protected:
const ExpressionType m_expression;
};
/** \returns an expression of the temporary version of *this.
*/
template<typename Derived>
inline const NestByValue<Derived>
DenseBase<Derived>::nestByValue() const
{
return NestByValue<Derived>(derived());
}
} // end namespace Eigen
#endif // EIGEN_NESTBYVALUE_H
external/eigen3/Eigen/src/Core/NoAlias.h 0 → 100644
View file @ a394b22a
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_NOALIAS_H
#define EIGEN_NOALIAS_H
namespace Eigen {
/** \class NoAlias
* \ingroup Core_Module
*
* \brief Pseudo expression providing an operator = assuming no aliasing
*
* \tparam ExpressionType the type of the object on which to do the lazy assignment
*
* This class represents an expression with special assignment operators
* assuming no aliasing between the target expression and the source expression.
* More precisely it alloas to bypass the EvalBeforeAssignBit flag of the source expression.
* It is the return type of MatrixBase::noalias()
* and most of the time this is the only way it is used.
*
* \sa MatrixBase::noalias()
*/
template<typename ExpressionType, template <typename> class StorageBase>
class NoAlias
{
public:
typedef typename ExpressionType::Scalar Scalar;
explicit NoAlias(ExpressionType& expression) : m_expression(expression) {}
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE ExpressionType& operator=(const StorageBase<OtherDerived>& other)
{
call_assignment_no_alias(m_expression, other.derived(), internal::assign_op<Scalar,typename OtherDerived::Scalar>());
return m_expression;
}
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE ExpressionType& operator+=(const StorageBase<OtherDerived>& other)
{
call_assignment_no_alias(m_expression, other.derived(), internal::add_assign_op<Scalar,typename OtherDerived::Scalar>());
return m_expression;
}
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE ExpressionType& operator-=(const StorageBase<OtherDerived>& other)
{
call_assignment_no_alias(m_expression, other.derived(), internal::sub_assign_op<Scalar,typename OtherDerived::Scalar>());
return m_expression;
}
EIGEN_DEVICE_FUNC
ExpressionType& expression() const
{
return m_expression;
}
protected:
ExpressionType& m_expression;
};
/** \returns a pseudo expression of \c *this with an operator= assuming
* no aliasing between \c *this and the source expression.
*
* More precisely, noalias() allows to bypass the EvalBeforeAssignBit flag.
* Currently, even though several expressions may alias, only product
* expressions have this flag. Therefore, noalias() is only usefull when
* the source expression contains a matrix product.
*
* Here are some examples where noalias is usefull:
* \code
* D.noalias() = A * B;
* D.noalias() += A.transpose() * B;
* D.noalias() -= 2 * A * B.adjoint();
* \endcode
*
* On the other hand the following example will lead to a \b wrong result:
* \code
* A.noalias() = A * B;
* \endcode
* because the result matrix A is also an operand of the matrix product. Therefore,
* there is no alternative than evaluating A * B in a temporary, that is the default
* behavior when you write:
* \code
* A = A * B;
* \endcode
*
* \sa class NoAlias
*/
template<typename Derived>
NoAlias<Derived,MatrixBase> MatrixBase<Derived>::noalias()
{
return NoAlias<Derived, Eigen::MatrixBase >(derived());
}
} // end namespace Eigen
#endif // EIGEN_NOALIAS_H
external/eigen3/Eigen/src/Core/NumTraits.h 0 → 100644
View file @ a394b22a
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2006-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_NUMTRAITS_H
#define EIGEN_NUMTRAITS_H
namespace Eigen {
namespace internal {
// default implementation of digits10(), based on numeric_limits if specialized,
// 0 for integer types, and log10(epsilon()) otherwise.
template< typename T,
bool use_numeric_limits = std::numeric_limits<T>::is_specialized,
bool is_integer = NumTraits<T>::IsInteger>
struct default_digits10_impl
{
static int run() { return std::numeric_limits<T>::digits10; }
};
template<typename T>
struct default_digits10_impl<T,false,false> // Floating point
{
static int run() {
using std::log10;
using std::ceil;
typedef typename NumTraits<T>::Real Real;
return int(ceil(-log10(NumTraits<Real>::epsilon())));
}
};
template<typename T>
struct default_digits10_impl<T,false,true> // Integer
{
static int run() { return 0; }
};
} // end namespace internal
/** \class NumTraits
* \ingroup Core_Module
*
* \brief Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
*
* \tparam T the numeric type at hand
*
* This class stores enums, typedefs and static methods giving information about a numeric type.
*
* The provided data consists of:
* \li A typedef \c Real, giving the "real part" type of \a T. If \a T is already real,
* then \c Real is just a typedef to \a T. If \a T is \c std::complex<U> then \c Real
* is a typedef to \a U.
* \li A typedef \c NonInteger, giving the type that should be used for operations producing non-integral values,
* such as quotients, square roots, etc. If \a T is a floating-point type, then this typedef just gives
* \a T again. Note however that many Eigen functions such as internal::sqrt simply refuse to
* take integers. Outside of a few cases, Eigen doesn't do automatic type promotion. Thus, this typedef is
* only intended as a helper for code that needs to explicitly promote types.
* \li A typedef \c Literal giving the type to use for numeric literals such as "2" or "0.5". For instance, for \c std::complex<U>, Literal is defined as \c U.
* Of course, this type must be fully compatible with \a T. In doubt, just use \a T here.
* \li A typedef \a Nested giving the type to use to nest a value inside of the expression tree. If you don't know what
* this means, just use \a T here.
* \li An enum value \a IsComplex. It is equal to 1 if \a T is a \c std::complex
* type, and to 0 otherwise.
* \li An enum value \a IsInteger. It is equal to \c 1 if \a T is an integer type such as \c int,
* and to \c 0 otherwise.
* \li Enum values ReadCost, AddCost and MulCost representing a rough estimate of the number of CPU cycles needed
* to by move / add / mul instructions respectively, assuming the data is already stored in CPU registers.
* Stay vague here. No need to do architecture-specific stuff.
* \li An enum value \a IsSigned. It is equal to \c 1 if \a T is a signed type and to 0 if \a T is unsigned.
* \li An enum value \a RequireInitialization. It is equal to \c 1 if the constructor of the numeric type \a T must
* be called, and to 0 if it is safe not to call it. Default is 0 if \a T is an arithmetic type, and 1 otherwise.
* \li An epsilon() function which, unlike <a href="http://en.cppreference.com/w/cpp/types/numeric_limits/epsilon">std::numeric_limits::epsilon()</a>,
* it returns a \a Real instead of a \a T.
* \li A dummy_precision() function returning a weak epsilon value. It is mainly used as a default
* value by the fuzzy comparison operators.
* \li highest() and lowest() functions returning the highest and lowest possible values respectively.
* \li digits10() function returning the number of decimal digits that can be represented without change. This is
* the analogue of <a href="http://en.cppreference.com/w/cpp/types/numeric_limits/digits10">std::numeric_limits<T>::digits10</a>
* which is used as the default implementation if specialized.
*/
template<typename T> struct GenericNumTraits
{
enum {
IsInteger = std::numeric_limits<T>::is_integer,
IsSigned = std::numeric_limits<T>::is_signed,
IsComplex = 0,
RequireInitialization = internal::is_arithmetic<T>::value ? 0 : 1,
ReadCost = 1,
AddCost = 1,
MulCost = 1
};
typedef T Real;
typedef typename internal::conditional<
IsInteger,
typename internal::conditional<sizeof(T)<=2, float, double>::type,
T
>::type NonInteger;
typedef T Nested;
typedef T Literal;
EIGEN_DEVICE_FUNC
static inline Real epsilon()
{
return numext::numeric_limits<T>::epsilon();
}
EIGEN_DEVICE_FUNC
static inline int digits10()
{
return internal::default_digits10_impl<T>::run();
}
EIGEN_DEVICE_FUNC
static inline Real dummy_precision()
{
// make sure to override this for floating-point types
return Real(0);
}
EIGEN_DEVICE_FUNC
static inline T highest() {
return (numext::numeric_limits<T>::max)();
}
EIGEN_DEVICE_FUNC
static inline T lowest() {
return IsInteger ? (numext::numeric_limits<T>::min)() : (-(numext::numeric_limits<T>::max)());
}
EIGEN_DEVICE_FUNC
static inline T infinity() {
return numext::numeric_limits<T>::infinity();
}
EIGEN_DEVICE_FUNC
static inline T quiet_NaN() {
return numext::numeric_limits<T>::quiet_NaN();
}
};
template<typename T> struct NumTraits : GenericNumTraits<T>
{};
template<> struct NumTraits<float>
: GenericNumTraits<float>
{
EIGEN_DEVICE_FUNC
static inline float dummy_precision() { return 1e-5f; }
};
template<> struct NumTraits<double> : GenericNumTraits<double>
{
EIGEN_DEVICE_FUNC
static inline double dummy_precision() { return 1e-12; }
};
template<> struct NumTraits<long double>
: GenericNumTraits<long double>
{
static inline long double dummy_precision() { return 1e-15l; }
};
template<typename _Real> struct NumTraits<std::complex<_Real> >
: GenericNumTraits<std::complex<_Real> >
{
typedef _Real Real;
typedef typename NumTraits<_Real>::Literal Literal;
enum {
IsComplex = 1,
RequireInitialization = NumTraits<_Real>::RequireInitialization,
ReadCost = 2 * NumTraits<_Real>::ReadCost,
AddCost = 2 * NumTraits<Real>::AddCost,
MulCost = 4 * NumTraits<Real>::MulCost + 2 * NumTraits<Real>::AddCost
};
EIGEN_DEVICE_FUNC
static inline Real epsilon() { return NumTraits<Real>::epsilon(); }
EIGEN_DEVICE_FUNC
static inline Real dummy_precision() { return NumTraits<Real>::dummy_precision(); }
EIGEN_DEVICE_FUNC
static inline int digits10() { return NumTraits<Real>::digits10(); }
};
template<typename Scalar, int Rows, int Cols, int Options, int MaxRows, int MaxCols>
struct NumTraits<Array<Scalar, Rows, Cols, Options, MaxRows, MaxCols> >
{
typedef Array<Scalar, Rows, Cols, Options, MaxRows, MaxCols> ArrayType;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef Array<RealScalar, Rows, Cols, Options, MaxRows, MaxCols> Real;
typedef typename NumTraits<Scalar>::NonInteger NonIntegerScalar;
typedef Array<NonIntegerScalar, Rows, Cols, Options, MaxRows, MaxCols> NonInteger;
typedef ArrayType & Nested;
typedef typename NumTraits<Scalar>::Literal Literal;
enum {
IsComplex = NumTraits<Scalar>::IsComplex,
IsInteger = NumTraits<Scalar>::IsInteger,
IsSigned = NumTraits<Scalar>::IsSigned,
RequireInitialization = 1,
ReadCost = ArrayType::SizeAtCompileTime==Dynamic ? HugeCost : ArrayType::SizeAtCompileTime * NumTraits<Scalar>::ReadCost,
AddCost = ArrayType::SizeAtCompileTime==Dynamic ? HugeCost : ArrayType::SizeAtCompileTime * NumTraits<Scalar>::AddCost,
MulCost = ArrayType::SizeAtCompileTime==Dynamic ? HugeCost : ArrayType::SizeAtCompileTime * NumTraits<Scalar>::MulCost
};
EIGEN_DEVICE_FUNC
static inline RealScalar epsilon() { return NumTraits<RealScalar>::epsilon(); }
EIGEN_DEVICE_FUNC
static inline RealScalar dummy_precision() { return NumTraits<RealScalar>::dummy_precision(); }
static inline int digits10() { return NumTraits<Scalar>::digits10(); }
};
template<> struct NumTraits<std::string>
: GenericNumTraits<std::string>
{
enum {
RequireInitialization = 1,
ReadCost = HugeCost,
AddCost = HugeCost,
MulCost = HugeCost
};
static inline int digits10() { return 0; }
private:
static inline std::string epsilon();
static inline std::string dummy_precision();
static inline std::string lowest();
static inline std::string highest();
static inline std::string infinity();
static inline std::string quiet_NaN();
};
// Empty specialization for void to allow template specialization based on NumTraits<T>::Real with T==void and SFINAE.
template<> struct NumTraits<void> {};
} // end namespace Eigen
#endif // EIGEN_NUMTRAITS_H
external/eigen3/Eigen/src/Core/PermutationMatrix.h 0 → 100644
View file @ a394b22a
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2009-2015 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_PERMUTATIONMATRIX_H
#define EIGEN_PERMUTATIONMATRIX_H
namespace Eigen {
namespace internal {
enum PermPermProduct_t {PermPermProduct};
} // end namespace internal
/** \class PermutationBase
* \ingroup Core_Module
*
* \brief Base class for permutations
*
* \tparam Derived the derived class
*
* This class is the base class for all expressions representing a permutation matrix,
* internally stored as a vector of integers.
* The convention followed here is that if \f$ \sigma \f$ is a permutation, the corresponding permutation matrix
* \f$ P_\sigma \f$ is such that if \f$ (e_1,\ldots,e_p) \f$ is the canonical basis, we have:
* \f[ P_\sigma(e_i) = e_{\sigma(i)}. \f]
* This convention ensures that for any two permutations \f$ \sigma, \tau \f$, we have:
* \f[ P_{\sigma\circ\tau} = P_\sigma P_\tau. \f]
*
* Permutation matrices are square and invertible.
*
* Notice that in addition to the member functions and operators listed here, there also are non-member
* operator* to multiply any kind of permutation object with any kind of matrix expression (MatrixBase)
* on either side.
*
* \sa class PermutationMatrix, class PermutationWrapper
*/
template<typename Derived>
class PermutationBase : public EigenBase<Derived>
{
typedef internal::traits<Derived> Traits;
typedef EigenBase<Derived> Base;
public:
#ifndef EIGEN_PARSED_BY_DOXYGEN
typedef typename Traits::IndicesType IndicesType;
enum {
Flags = Traits::Flags,
RowsAtCompileTime = Traits::RowsAtCompileTime,
ColsAtCompileTime = Traits::ColsAtCompileTime,
MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
};
typedef typename Traits::StorageIndex StorageIndex;
typedef Matrix<StorageIndex,RowsAtCompileTime,ColsAtCompileTime,0,MaxRowsAtCompileTime,MaxColsAtCompileTime>
DenseMatrixType;
typedef PermutationMatrix<IndicesType::SizeAtCompileTime,IndicesType::MaxSizeAtCompileTime,StorageIndex>
PlainPermutationType;
typedef PlainPermutationType PlainObject;
using Base::derived;
typedef Inverse<Derived> InverseReturnType;
typedef void Scalar;
#endif
/** Copies the other permutation into *this */
template<typename OtherDerived>
Derived& operator=(const PermutationBase<OtherDerived>& other)
{
indices() = other.indices();
return derived();
}
/** Assignment from the Transpositions \a tr */
template<typename OtherDerived>
Derived& operator=(const TranspositionsBase<OtherDerived>& tr)
{
setIdentity(tr.size());
for(Index k=size()-1; k>=0; --k)
applyTranspositionOnTheRight(k,tr.coeff(k));
return derived();
}
#ifndef EIGEN_PARSED_BY_DOXYGEN
/** This is a special case of the templated operator=. Its purpose is to
* prevent a default operator= from hiding the templated operator=.
*/
Derived& operator=(const PermutationBase& other)
{
indices() = other.indices();
return derived();
}
#endif
/** \returns the number of rows */
inline Index rows() const { return Index(indices().size()); }
/** \returns the number of columns */
inline Index cols() const { return Index(indices().size()); }
/** \returns the size of a side of the respective square matrix, i.e., the number of indices */
inline Index size() const { return Index(indices().size()); }
#ifndef EIGEN_PARSED_BY_DOXYGEN
template<typename DenseDerived>
void evalTo(MatrixBase<DenseDerived>& other) const
{
other.setZero();
for (Index i=0; i<rows(); ++i)
other.coeffRef(indices().coeff(i),i) = typename DenseDerived::Scalar(1);
}
#endif
/** \returns a Matrix object initialized from this permutation matrix. Notice that it
* is inefficient to return this Matrix object by value. For efficiency, favor using
* the Matrix constructor taking EigenBase objects.
*/
DenseMatrixType toDenseMatrix() const
{
return derived();
}
/** const version of indices(). */
const IndicesType& indices() const { return derived().indices(); }
/** \returns a reference to the stored array representing the permutation. */
IndicesType& indices() { return derived().indices(); }
/** Resizes to given size.
*/
inline void resize(Index newSize)
{
indices().resize(newSize);
}
/** Sets *this to be the identity permutation matrix */
void setIdentity()
{
StorageIndex n = StorageIndex(size());
for(StorageIndex i = 0; i < n; ++i)
indices().coeffRef(i) = i;
}
/** Sets *this to be the identity permutation matrix of given size.
*/
void setIdentity(Index newSize)
{
resize(newSize);
setIdentity();
}
/** Multiplies *this by the transposition \f$(ij)\f$ on the left.
*
* \returns a reference to *this.
*
* \warning This is much slower than applyTranspositionOnTheRight(Index,Index):
* this has linear complexity and requires a lot of branching.
*
* \sa applyTranspositionOnTheRight(Index,Index)
*/
Derived& applyTranspositionOnTheLeft(Index i, Index j)
{
eigen_assert(i>=0 && j>=0 && i<size() && j<size());
for(Index k = 0; k < size(); ++k)
{
if(indices().coeff(k) == i) indices().coeffRef(k) = StorageIndex(j);
else if(indices().coeff(k) == j) indices().coeffRef(k) = StorageIndex(i);
}
return derived();
}
/** Multiplies *this by the transposition \f$(ij)\f$ on the right.
*
* \returns a reference to *this.
*
* This is a fast operation, it only consists in swapping two indices.
*
* \sa applyTranspositionOnTheLeft(Index,Index)
*/
Derived& applyTranspositionOnTheRight(Index i, Index j)
{
eigen_assert(i>=0 && j>=0 && i<size() && j<size());
std::swap(indices().coeffRef(i), indices().coeffRef(j));
return derived();
}
/** \returns the inverse permutation matrix.
*
* \note \blank \note_try_to_help_rvo
*/
inline InverseReturnType inverse() const
{ return InverseReturnType(derived()); }
/** \returns the tranpose permutation matrix.
*
* \note \blank \note_try_to_help_rvo
*/
inline InverseReturnType transpose() const
{ return InverseReturnType(derived()); }
/**** multiplication helpers to hopefully get RVO ****/
#ifndef EIGEN_PARSED_BY_DOXYGEN
protected:
template<typename OtherDerived>
void assignTranspose(const PermutationBase<OtherDerived>& other)
{
for (Index i=0; i<rows();++i) indices().coeffRef(other.indices().coeff(i)) = i;
}
template<typename Lhs,typename Rhs>
void assignProduct(const Lhs& lhs, const Rhs& rhs)
{
eigen_assert(lhs.cols() == rhs.rows());
for (Index i=0; i<rows();++i) indices().coeffRef(i) = lhs.indices().coeff(rhs.indices().coeff(i));
}
#endif
public:
/** \returns the product permutation matrix.
*
* \note \blank \note_try_to_help_rvo
*/
template<typename Other>
inline PlainPermutationType operator*(const PermutationBase<Other>& other) const
{ return PlainPermutationType(internal::PermPermProduct, derived(), other.derived()); }
/** \returns the product of a permutation with another inverse permutation.
*
* \note \blank \note_try_to_help_rvo
*/
template<typename Other>
inline PlainPermutationType operator*(const InverseImpl<Other,PermutationStorage>& other) const
{ return PlainPermutationType(internal::PermPermProduct, *this, other.eval()); }
/** \returns the product of an inverse permutation with another permutation.
*
* \note \blank \note_try_to_help_rvo
*/
template<typename Other> friend
inline PlainPermutationType operator*(const InverseImpl<Other, PermutationStorage>& other, const PermutationBase& perm)
{ return PlainPermutationType(internal::PermPermProduct, other.eval(), perm); }
/** \returns the determinant of the permutation matrix, which is either 1 or -1 depending on the parity of the permutation.
*
* This function is O(\c n) procedure allocating a buffer of \c n booleans.
*/
Index determinant() const
{
Index res = 1;
Index n = size();
Matrix<bool,RowsAtCompileTime,1,0,MaxRowsAtCompileTime> mask(n);
mask.fill(false);
Index r = 0;
while(r < n)
{
// search for the next seed
while(r<n && mask[r]) r++;
if(r>=n)
break;
// we got one, let's follow it until we are back to the seed
Index k0 = r++;
mask.coeffRef(k0) = true;
for(Index k=indices().coeff(k0); k!=k0; k=indices().coeff(k))
{
mask.coeffRef(k) = true;
res = -res;
}
}
return res;
}
protected:
};
namespace internal {
template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex>
struct traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex> >
: traits<Matrix<_StorageIndex,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
{
typedef PermutationStorage StorageKind;
typedef Matrix<_StorageIndex, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
typedef _StorageIndex StorageIndex;
typedef void Scalar;
};
}
/** \class PermutationMatrix
* \ingroup Core_Module
*
* \brief Permutation matrix
*
* \tparam SizeAtCompileTime the number of rows/cols, or Dynamic
* \tparam MaxSizeAtCompileTime the maximum number of rows/cols, or Dynamic. This optional parameter defaults to SizeAtCompileTime. Most of the time, you should not have to specify it.
* \tparam _StorageIndex the integer type of the indices
*
* This class represents a permutation matrix, internally stored as a vector of integers.
*
* \sa class PermutationBase, class PermutationWrapper, class DiagonalMatrix
*/
template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex>
class PermutationMatrix : public PermutationBase<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex> >
{
typedef PermutationBase<PermutationMatrix> Base;
typedef internal::traits<PermutationMatrix> Traits;
public:
typedef const PermutationMatrix& Nested;
#ifndef EIGEN_PARSED_BY_DOXYGEN
typedef typename Traits::IndicesType IndicesType;
typedef typename Traits::StorageIndex StorageIndex;
#endif
inline PermutationMatrix()
{}
/** Constructs an uninitialized permutation matrix of given size.
*/
explicit inline PermutationMatrix(Index size) : m_indices(size)
{
eigen_internal_assert(size <= NumTraits<StorageIndex>::highest());
}
/** Copy constructor. */
template<typename OtherDerived>
inline PermutationMatrix(const PermutationBase<OtherDerived>& other)
: m_indices(other.indices()) {}
#ifndef EIGEN_PARSED_BY_DOXYGEN
/** Standard copy constructor. Defined only to prevent a default copy constructor
* from hiding the other templated constructor */
inline PermutationMatrix(const PermutationMatrix& other) : m_indices(other.indices()) {}
#endif
/** Generic constructor from expression of the indices. The indices
* array has the meaning that the permutations sends each integer i to indices[i].
*
* \warning It is your responsibility to check that the indices array that you passes actually
* describes a permutation, i.e., each value between 0 and n-1 occurs exactly once, where n is the
* array's size.
*/
template<typename Other>
explicit inline PermutationMatrix(const MatrixBase<Other>& indices) : m_indices(indices)
{}
/** Convert the Transpositions \a tr to a permutation matrix */
template<typename Other>
explicit PermutationMatrix(const TranspositionsBase<Other>& tr)
: m_indices(tr.size())
{
*this = tr;
}
/** Copies the other permutation into *this */
template<typename Other>
PermutationMatrix& operator=(const PermutationBase<Other>& other)
{
m_indices = other.indices();
return *this;
}
/** Assignment from the Transpositions \a tr */
template<typename Other>
PermutationMatrix& operator=(const TranspositionsBase<Other>& tr)
{
return Base::operator=(tr.derived());
}
#ifndef EIGEN_PARSED_BY_DOXYGEN
/** This is a special case of the templated operator=. Its purpose is to
* prevent a default operator= from hiding the templated operator=.
*/
PermutationMatrix& operator=(const PermutationMatrix& other)
{
m_indices = other.m_indices;
return *this;
}
#endif
/** const version of indices(). */
const IndicesType& indices() const { return m_indices; }
/** \returns a reference to the stored array representing the permutation. */
IndicesType& indices() { return m_indices; }
/**** multiplication helpers to hopefully get RVO ****/
#ifndef EIGEN_PARSED_BY_DOXYGEN
template<typename Other>
PermutationMatrix(const InverseImpl<Other,PermutationStorage>& other)
: m_indices(other.derived().nestedExpression().size())
{
eigen_internal_assert(m_indices.size() <= NumTraits<StorageIndex>::highest());
StorageIndex end = StorageIndex(m_indices.size());
for (StorageIndex i=0; i<end;++i)
m_indices.coeffRef(other.derived().nestedExpression().indices().coeff(i)) = i;
}
template<typename Lhs,typename Rhs>
PermutationMatrix(internal::PermPermProduct_t, const Lhs& lhs, const Rhs& rhs)
: m_indices(lhs.indices().size())
{
Base::assignProduct(lhs,rhs);
}
#endif
protected:
IndicesType m_indices;
};
namespace internal {
template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex, int _PacketAccess>
struct traits<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex>,_PacketAccess> >
: traits<Matrix<_StorageIndex,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
{
typedef PermutationStorage StorageKind;
typedef Map<const Matrix<_StorageIndex, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1>, _PacketAccess> IndicesType;
typedef _StorageIndex StorageIndex;
typedef void Scalar;
};
}
template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex, int _PacketAccess>
class Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex>,_PacketAccess>
: public PermutationBase<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex>,_PacketAccess> >
{
typedef PermutationBase<Map> Base;
typedef internal::traits<Map> Traits;
public:
#ifndef EIGEN_PARSED_BY_DOXYGEN
typedef typename Traits::IndicesType IndicesType;
typedef typename IndicesType::Scalar StorageIndex;
#endif
inline Map(const StorageIndex* indicesPtr)
: m_indices(indicesPtr)
{}
inline Map(const StorageIndex* indicesPtr, Index size)
: m_indices(indicesPtr,size)
{}
/** Copies the other permutation into *this */
template<typename Other>
Map& operator=(const PermutationBase<Other>& other)
{ return Base::operator=(other.derived()); }
/** Assignment from the Transpositions \a tr */
template<typename Other>
Map& operator=(const TranspositionsBase<Other>& tr)
{ return Base::operator=(tr.derived()); }
#ifndef EIGEN_PARSED_BY_DOXYGEN
/** This is a special case of the templated operator=. Its purpose is to
* prevent a default operator= from hiding the templated operator=.
*/
Map& operator=(const Map& other)
{
m_indices = other.m_indices;
return *this;
}
#endif
/** const version of indices(). */
const IndicesType& indices() const { return m_indices; }
/** \returns a reference to the stored array representing the permutation. */
IndicesType& indices() { return m_indices; }
protected:
IndicesType m_indices;
};
template<typename _IndicesType> class TranspositionsWrapper;
namespace internal {
template<typename _IndicesType>
struct traits<PermutationWrapper<_IndicesType> >
{
typedef PermutationStorage StorageKind;
typedef void Scalar;
typedef typename _IndicesType::Scalar StorageIndex;
typedef _IndicesType IndicesType;
enum {
RowsAtCompileTime = _IndicesType::SizeAtCompileTime,
ColsAtCompileTime = _IndicesType::SizeAtCompileTime,
MaxRowsAtCompileTime = IndicesType::MaxSizeAtCompileTime,
MaxColsAtCompileTime = IndicesType::MaxSizeAtCompileTime,
Flags = 0
};
};
}
/** \class PermutationWrapper
* \ingroup Core_Module
*
* \brief Class to view a vector of integers as a permutation matrix
*
* \tparam _IndicesType the type of the vector of integer (can be any compatible expression)
*
* This class allows to view any vector expression of integers as a permutation matrix.
*
* \sa class PermutationBase, class PermutationMatrix
*/
template<typename _IndicesType>
class PermutationWrapper : public PermutationBase<PermutationWrapper<_IndicesType> >
{
typedef PermutationBase<PermutationWrapper> Base;
typedef internal::traits<PermutationWrapper> Traits;
public:
#ifndef EIGEN_PARSED_BY_DOXYGEN
typedef typename Traits::IndicesType IndicesType;
#endif
inline PermutationWrapper(const IndicesType& indices)
: m_indices(indices)
{}
/** const version of indices(). */
const typename internal::remove_all<typename IndicesType::Nested>::type&
indices() const { return m_indices; }
protected:
typename IndicesType::Nested m_indices;
};
/** \returns the matrix with the permutation applied to the columns.
*/
template<typename MatrixDerived, typename PermutationDerived>
EIGEN_DEVICE_FUNC
const Product<MatrixDerived, PermutationDerived, AliasFreeProduct>
operator*(const MatrixBase<MatrixDerived> &matrix,
const PermutationBase<PermutationDerived>& permutation)
{
return Product<MatrixDerived, PermutationDerived, AliasFreeProduct>
(matrix.derived(), permutation.derived());
}
/** \returns the matrix with the permutation applied to the rows.
*/
template<typename PermutationDerived, typename MatrixDerived>
EIGEN_DEVICE_FUNC
const Product<PermutationDerived, MatrixDerived, AliasFreeProduct>
operator*(const PermutationBase<PermutationDerived> &permutation,
const MatrixBase<MatrixDerived>& matrix)
{
return Product<PermutationDerived, MatrixDerived, AliasFreeProduct>
(permutation.derived(), matrix.derived());
}
template<typename PermutationType>
class InverseImpl<PermutationType, PermutationStorage>
: public EigenBase<Inverse<PermutationType> >
{
typedef typename PermutationType::PlainPermutationType PlainPermutationType;
typedef internal::traits<PermutationType> PermTraits;
protected:
InverseImpl() {}
public:
typedef Inverse<PermutationType> InverseType;
using EigenBase<Inverse<PermutationType> >::derived;
#ifndef EIGEN_PARSED_BY_DOXYGEN
typedef typename PermutationType::DenseMatrixType DenseMatrixType;
enum {
RowsAtCompileTime = PermTraits::RowsAtCompileTime,
ColsAtCompileTime = PermTraits::ColsAtCompileTime,
MaxRowsAtCompileTime = PermTraits::MaxRowsAtCompileTime,
MaxColsAtCompileTime = PermTraits::MaxColsAtCompileTime
};
#endif
#ifndef EIGEN_PARSED_BY_DOXYGEN
template<typename DenseDerived>
void evalTo(MatrixBase<DenseDerived>& other) const
{
other.setZero();
for (Index i=0; i<derived().rows();++i)
other.coeffRef(i, derived().nestedExpression().indices().coeff(i)) = typename DenseDerived::Scalar(1);
}
#endif
/** \return the equivalent permutation matrix */
PlainPermutationType eval() const { return derived(); }
DenseMatrixType toDenseMatrix() const { return derived(); }
/** \returns the matrix with the inverse permutation applied to the columns.
*/
template<typename OtherDerived> friend
const Product<OtherDerived, InverseType, AliasFreeProduct>
operator*(const MatrixBase<OtherDerived>& matrix, const InverseType& trPerm)
{
return Product<OtherDerived, InverseType, AliasFreeProduct>(matrix.derived(), trPerm.derived());
}
/** \returns the matrix with the inverse permutation applied to the rows.
*/
template<typename OtherDerived>
const Product<InverseType, OtherDerived, AliasFreeProduct>
operator*(const MatrixBase<OtherDerived>& matrix) const
{
return Product<InverseType, OtherDerived, AliasFreeProduct>(derived(), matrix.derived());
}
};
template<typename Derived>
const PermutationWrapper<const Derived> MatrixBase<Derived>::asPermutation() const
{
return derived();
}
namespace internal {
template<> struct AssignmentKind<DenseShape,PermutationShape> { typedef EigenBase2EigenBase Kind; };
} // end namespace internal
} // end namespace Eigen
#endif // EIGEN_PERMUTATIONMATRIX_H
external/eigen3/Eigen/src/Core/PlainObjectBase.h 0 → 100644
View file @ a394b22a
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_DENSESTORAGEBASE_H
#define EIGEN_DENSESTORAGEBASE_H
#if defined(EIGEN_INITIALIZE_MATRICES_BY_ZERO)
# define EIGEN_INITIALIZE_COEFFS
# define EIGEN_INITIALIZE_COEFFS_IF_THAT_OPTION_IS_ENABLED for(int i=0;i<base().size();++i) coeffRef(i)=Scalar(0);
#elif defined(EIGEN_INITIALIZE_MATRICES_BY_NAN)
# define EIGEN_INITIALIZE_COEFFS
# define EIGEN_INITIALIZE_COEFFS_IF_THAT_OPTION_IS_ENABLED for(int i=0;i<base().size();++i) coeffRef(i)=std::numeric_limits<Scalar>::quiet_NaN();
#else
# undef EIGEN_INITIALIZE_COEFFS
# define EIGEN_INITIALIZE_COEFFS_IF_THAT_OPTION_IS_ENABLED
#endif
namespace Eigen {
namespace internal {
template<int MaxSizeAtCompileTime> struct check_rows_cols_for_overflow {
template<typename Index>
EIGEN_DEVICE_FUNC
static EIGEN_ALWAYS_INLINE void run(Index, Index)
{
}
};
template<> struct check_rows_cols_for_overflow<Dynamic> {
template<typename Index>
EIGEN_DEVICE_FUNC
static EIGEN_ALWAYS_INLINE void run(Index rows, Index cols)
{
// http://hg.mozilla.org/mozilla-central/file/6c8a909977d3/xpcom/ds/CheckedInt.h#l242
// we assume Index is signed
Index max_index = (std::size_t(1) << (8 * sizeof(Index) - 1)) - 1; // assume Index is signed
bool error = (rows == 0 || cols == 0) ? false
: (rows > max_index / cols);
if (error)
throw_std_bad_alloc();
}
};
template <typename Derived,
typename OtherDerived = Derived,
bool IsVector = bool(Derived::IsVectorAtCompileTime) && bool(OtherDerived::IsVectorAtCompileTime)>
struct conservative_resize_like_impl;
template<typename MatrixTypeA, typename MatrixTypeB, bool SwapPointers> struct matrix_swap_impl;
} // end namespace internal
#ifdef EIGEN_PARSED_BY_DOXYGEN
namespace doxygen {
// This is a workaround to doxygen not being able to understand the inheritance logic
// when it is hidden by the dense_xpr_base helper struct.
// Moreover, doxygen fails to include members that are not documented in the declaration body of
// MatrixBase if we inherits MatrixBase<Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> >,
// this is why we simply inherits MatrixBase, though this does not make sense.
/** This class is just a workaround for Doxygen and it does not not actually exist. */
template<typename Derived> struct dense_xpr_base_dispatcher;
/** This class is just a workaround for Doxygen and it does not not actually exist. */
template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
struct dense_xpr_base_dispatcher<Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> >
: public MatrixBase {};
/** This class is just a workaround for Doxygen and it does not not actually exist. */
template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
struct dense_xpr_base_dispatcher<Array<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> >
: public ArrayBase {};
} // namespace doxygen
/** \class PlainObjectBase
* \ingroup Core_Module
* \brief %Dense storage base class for matrices and arrays.
*
* This class can be extended with the help of the plugin mechanism described on the page
* \ref TopicCustomizing_Plugins by defining the preprocessor symbol \c EIGEN_PLAINOBJECTBASE_PLUGIN.
*
* \tparam Derived is the derived type, e.g., a Matrix or Array
*
* \sa \ref TopicClassHierarchy
*/
template<typename Derived>
class PlainObjectBase : public doxygen::dense_xpr_base_dispatcher<Derived>
#else
template<typename Derived>
class PlainObjectBase : public internal::dense_xpr_base<Derived>::type
#endif
{
public:
enum { Options = internal::traits<Derived>::Options };
typedef typename internal::dense_xpr_base<Derived>::type Base;
typedef typename internal::traits<Derived>::StorageKind StorageKind;
typedef typename internal::traits<Derived>::Scalar Scalar;
typedef typename internal::packet_traits<Scalar>::type PacketScalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef Derived DenseType;
using Base::RowsAtCompileTime;
using Base::ColsAtCompileTime;
using Base::SizeAtCompileTime;
using Base::MaxRowsAtCompileTime;
using Base::MaxColsAtCompileTime;
using Base::MaxSizeAtCompileTime;
using Base::IsVectorAtCompileTime;
using Base::Flags;
template<typename PlainObjectType, int MapOptions, typename StrideType> friend class Eigen::Map;
friend class Eigen::Map<Derived, Unaligned>;
typedef Eigen::Map<Derived, Unaligned> MapType;
friend class Eigen::Map<const Derived, Unaligned>;
typedef const Eigen::Map<const Derived, Unaligned> ConstMapType;
#if EIGEN_MAX_ALIGN_BYTES>0
// for EIGEN_MAX_ALIGN_BYTES==0, AlignedMax==Unaligned, and many compilers generate warnings for friend-ing a class twice.
friend class Eigen::Map<Derived, AlignedMax>;
friend class Eigen::Map<const Derived, AlignedMax>;
#endif
typedef Eigen::Map<Derived, AlignedMax> AlignedMapType;
typedef const Eigen::Map<const Derived, AlignedMax> ConstAlignedMapType;
template<typename StrideType> struct StridedMapType { typedef Eigen::Map<Derived, Unaligned, StrideType> type; };
template<typename StrideType> struct StridedConstMapType { typedef Eigen::Map<const Derived, Unaligned, StrideType> type; };
template<typename StrideType> struct StridedAlignedMapType { typedef Eigen::Map<Derived, AlignedMax, StrideType> type; };
template<typename StrideType> struct StridedConstAlignedMapType { typedef Eigen::Map<const Derived, AlignedMax, StrideType> type; };
protected:
DenseStorage<Scalar, Base::MaxSizeAtCompileTime, Base::RowsAtCompileTime, Base::ColsAtCompileTime, Options> m_storage;
public:
enum { NeedsToAlign = (SizeAtCompileTime != Dynamic) && (internal::traits<Derived>::Alignment>0) };
EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(NeedsToAlign)
EIGEN_DEVICE_FUNC
Base& base() { return *static_cast<Base*>(this); }
EIGEN_DEVICE_FUNC
const Base& base() const { return *static_cast<const Base*>(this); }
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE Index rows() const { return m_storage.rows(); }
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE Index cols() const { return m_storage.cols(); }
/** This is an overloaded version of DenseCoeffsBase<Derived,ReadOnlyAccessors>::coeff(Index,Index) const
* provided to by-pass the creation of an evaluator of the expression, thus saving compilation efforts.
*
* See DenseCoeffsBase<Derived,ReadOnlyAccessors>::coeff(Index) const for details. */
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE const Scalar& coeff(Index rowId, Index colId) const
{
if(Flags & RowMajorBit)
return m_storage.data()[colId + rowId * m_storage.cols()];
else // column-major
return m_storage.data()[rowId + colId * m_storage.rows()];
}
/** This is an overloaded version of DenseCoeffsBase<Derived,ReadOnlyAccessors>::coeff(Index) const
* provided to by-pass the creation of an evaluator of the expression, thus saving compilation efforts.
*
* See DenseCoeffsBase<Derived,ReadOnlyAccessors>::coeff(Index) const for details. */
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE const Scalar& coeff(Index index) const
{
return m_storage.data()[index];
}
/** This is an overloaded version of DenseCoeffsBase<Derived,WriteAccessors>::coeffRef(Index,Index) const
* provided to by-pass the creation of an evaluator of the expression, thus saving compilation efforts.
*
* See DenseCoeffsBase<Derived,WriteAccessors>::coeffRef(Index,Index) const for details. */
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE Scalar& coeffRef(Index rowId, Index colId)
{
if(Flags & RowMajorBit)
return m_storage.data()[colId + rowId * m_storage.cols()];
else // column-major
return m_storage.data()[rowId + colId * m_storage.rows()];
}
/** This is an overloaded version of DenseCoeffsBase<Derived,WriteAccessors>::coeffRef(Index) const
* provided to by-pass the creation of an evaluator of the expression, thus saving compilation efforts.
*
* See DenseCoeffsBase<Derived,WriteAccessors>::coeffRef(Index) const for details. */
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE Scalar& coeffRef(Index index)
{
return m_storage.data()[index];
}
/** This is the const version of coeffRef(Index,Index) which is thus synonym of coeff(Index,Index).
* It is provided for convenience. */
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE const Scalar& coeffRef(Index rowId, Index colId) const
{
if(Flags & RowMajorBit)
return m_storage.data()[colId + rowId * m_storage.cols()];
else // column-major
return m_storage.data()[rowId + colId * m_storage.rows()];
}
/** This is the const version of coeffRef(Index) which is thus synonym of coeff(Index).
* It is provided for convenience. */
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE const Scalar& coeffRef(Index index) const
{
return m_storage.data()[index];
}
/** \internal */
template<int LoadMode>
EIGEN_STRONG_INLINE PacketScalar packet(Index rowId, Index colId) const
{
return internal::ploadt<PacketScalar, LoadMode>
(m_storage.data() + (Flags & RowMajorBit
? colId + rowId * m_storage.cols()
: rowId + colId * m_storage.rows()));
}
/** \internal */
template<int LoadMode>
EIGEN_STRONG_INLINE PacketScalar packet(Index index) const
{
return internal::ploadt<PacketScalar, LoadMode>(m_storage.data() + index);
}
/** \internal */
template<int StoreMode>
EIGEN_STRONG_INLINE void writePacket(Index rowId, Index colId, const PacketScalar& val)
{
internal::pstoret<Scalar, PacketScalar, StoreMode>
(m_storage.data() + (Flags & RowMajorBit
? colId + rowId * m_storage.cols()
: rowId + colId * m_storage.rows()), val);
}
/** \internal */
template<int StoreMode>
EIGEN_STRONG_INLINE void writePacket(Index index, const PacketScalar& val)
{
internal::pstoret<Scalar, PacketScalar, StoreMode>(m_storage.data() + index, val);
}
/** \returns a const pointer to the data array of this matrix */
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar *data() const
{ return m_storage.data(); }
/** \returns a pointer to the data array of this matrix */
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar *data()
{ return m_storage.data(); }
/** Resizes \c *this to a \a rows x \a cols matrix.
*
* This method is intended for dynamic-size matrices, although it is legal to call it on any
* matrix as long as fixed dimensions are left unchanged. If you only want to change the number
* of rows and/or of columns, you can use resize(NoChange_t, Index), resize(Index, NoChange_t).
*
* If the current number of coefficients of \c *this exactly matches the
* product \a rows * \a cols, then no memory allocation is performed and
* the current values are left unchanged. In all other cases, including
* shrinking, the data is reallocated and all previous values are lost.
*
* Example: \include Matrix_resize_int_int.cpp
* Output: \verbinclude Matrix_resize_int_int.out
*
* \sa resize(Index) for vectors, resize(NoChange_t, Index), resize(Index, NoChange_t)
*/
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
{
eigen_assert( EIGEN_IMPLIES(RowsAtCompileTime!=Dynamic,rows==RowsAtCompileTime)
&& EIGEN_IMPLIES(ColsAtCompileTime!=Dynamic,cols==ColsAtCompileTime)
&& EIGEN_IMPLIES(RowsAtCompileTime==Dynamic && MaxRowsAtCompileTime!=Dynamic,rows<=MaxRowsAtCompileTime)
&& EIGEN_IMPLIES(ColsAtCompileTime==Dynamic && MaxColsAtCompileTime!=Dynamic,cols<=MaxColsAtCompileTime)
&& rows>=0 && cols>=0 && "Invalid sizes when resizing a matrix or array.");
internal::check_rows_cols_for_overflow<MaxSizeAtCompileTime>::run(rows, cols);
#ifdef EIGEN_INITIALIZE_COEFFS
Index size = rows*cols;
bool size_changed = size != this->size();
m_storage.resize(size, rows, cols);
if(size_changed) EIGEN_INITIALIZE_COEFFS_IF_THAT_OPTION_IS_ENABLED
#else
m_storage.resize(rows*cols, rows, cols);
#endif
}
/** Resizes \c *this to a vector of length \a size
*
* \only_for_vectors. This method does not work for
* partially dynamic matrices when the static dimension is anything other
* than 1. For example it will not work with Matrix<double, 2, Dynamic>.
*
* Example: \include Matrix_resize_int.cpp
* Output: \verbinclude Matrix_resize_int.out
*
* \sa resize(Index,Index), resize(NoChange_t, Index), resize(Index, NoChange_t)
*/
EIGEN_DEVICE_FUNC
inline void resize(Index size)
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(PlainObjectBase)
eigen_assert(((SizeAtCompileTime == Dynamic && (MaxSizeAtCompileTime==Dynamic || size<=MaxSizeAtCompileTime)) || SizeAtCompileTime == size) && size>=0);
#ifdef EIGEN_INITIALIZE_COEFFS
bool size_changed = size != this->size();
#endif
if(RowsAtCompileTime == 1)
m_storage.resize(size, 1, size);
else
m_storage.resize(size, size, 1);
#ifdef EIGEN_INITIALIZE_COEFFS
if(size_changed) EIGEN_INITIALIZE_COEFFS_IF_THAT_OPTION_IS_ENABLED
#endif
}
/** Resizes the matrix, changing only the number of columns. For the parameter of type NoChange_t, just pass the special value \c NoChange
* as in the example below.
*
* Example: \include Matrix_resize_NoChange_int.cpp
* Output: \verbinclude Matrix_resize_NoChange_int.out
*
* \sa resize(Index,Index)
*/
EIGEN_DEVICE_FUNC
inline void resize(NoChange_t, Index cols)
{
resize(rows(), cols);
}
/** Resizes the matrix, changing only the number of rows. For the parameter of type NoChange_t, just pass the special value \c NoChange
* as in the example below.
*
* Example: \include Matrix_resize_int_NoChange.cpp
* Output: \verbinclude Matrix_resize_int_NoChange.out
*
* \sa resize(Index,Index)
*/
EIGEN_DEVICE_FUNC
inline void resize(Index rows, NoChange_t)
{
resize(rows, cols());
}
/** Resizes \c *this to have the same dimensions as \a other.
* Takes care of doing all the checking that's needed.
*
* Note that copying a row-vector into a vector (and conversely) is allowed.
* The resizing, if any, is then done in the appropriate way so that row-vectors
* remain row-vectors and vectors remain vectors.
*/
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE void resizeLike(const EigenBase<OtherDerived>& _other)
{
const OtherDerived& other = _other.derived();
internal::check_rows_cols_for_overflow<MaxSizeAtCompileTime>::run(other.rows(), other.cols());
const Index othersize = other.rows()*other.cols();
if(RowsAtCompileTime == 1)
{
eigen_assert(other.rows() == 1 || other.cols() == 1);
resize(1, othersize);
}
else if(ColsAtCompileTime == 1)
{
eigen_assert(other.rows() == 1 || other.cols() == 1);
resize(othersize, 1);
}
else resize(other.rows(), other.cols());
}
/** Resizes the matrix to \a rows x \a cols while leaving old values untouched.
*
* The method is intended for matrices of dynamic size. If you only want to change the number
* of rows and/or of columns, you can use conservativeResize(NoChange_t, Index) or
* conservativeResize(Index, NoChange_t).
*
* Matrices are resized relative to the top-left element. In case values need to be
* appended to the matrix they will be uninitialized.
*/
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE void conservativeResize(Index rows, Index cols)
{
internal::conservative_resize_like_impl<Derived>::run(*this, rows, cols);
}
/** Resizes the matrix to \a rows x \a cols while leaving old values untouched.
*
* As opposed to conservativeResize(Index rows, Index cols), this version leaves
* the number of columns unchanged.
*
* In case the matrix is growing, new rows will be uninitialized.
*/
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE void conservativeResize(Index rows, NoChange_t)
{
// Note: see the comment in conservativeResize(Index,Index)
conservativeResize(rows, cols());
}
/** Resizes the matrix to \a rows x \a cols while leaving old values untouched.
*
* As opposed to conservativeResize(Index rows, Index cols), this version leaves
* the number of rows unchanged.
*
* In case the matrix is growing, new columns will be uninitialized.
*/
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE void conservativeResize(NoChange_t, Index cols)
{
// Note: see the comment in conservativeResize(Index,Index)
conservativeResize(rows(), cols);
}
/** Resizes the vector to \a size while retaining old values.
*
* \only_for_vectors. This method does not work for
* partially dynamic matrices when the static dimension is anything other
* than 1. For example it will not work with Matrix<double, 2, Dynamic>.
*
* When values are appended, they will be uninitialized.
*/
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE void conservativeResize(Index size)
{
internal::conservative_resize_like_impl<Derived>::run(*this, size);
}
/** Resizes the matrix to \a rows x \a cols of \c other, while leaving old values untouched.
*
* The method is intended for matrices of dynamic size. If you only want to change the number
* of rows and/or of columns, you can use conservativeResize(NoChange_t, Index) or
* conservativeResize(Index, NoChange_t).
*
* Matrices are resized relative to the top-left element. In case values need to be
* appended to the matrix they will copied from \c other.
*/
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE void conservativeResizeLike(const DenseBase<OtherDerived>& other)
{
internal::conservative_resize_like_impl<Derived,OtherDerived>::run(*this, other);
}
/** This is a special case of the templated operator=. Its purpose is to
* prevent a default operator= from hiding the templated operator=.
*/
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE Derived& operator=(const PlainObjectBase& other)
{
return _set(other);
}
/** \sa MatrixBase::lazyAssign() */
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE Derived& lazyAssign(const DenseBase<OtherDerived>& other)
{
_resize_to_match(other);
return Base::lazyAssign(other.derived());
}
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE Derived& operator=(const ReturnByValue<OtherDerived>& func)
{
resize(func.rows(), func.cols());
return Base::operator=(func);
}
// Prevent user from trying to instantiate PlainObjectBase objects
// by making all its constructor protected. See bug 1074.
protected:
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE PlainObjectBase() : m_storage()
{
// _check_template_params();
// EIGEN_INITIALIZE_COEFFS_IF_THAT_OPTION_IS_ENABLED
}
#ifndef EIGEN_PARSED_BY_DOXYGEN
// FIXME is it still needed ?
/** \internal */
EIGEN_DEVICE_FUNC
explicit PlainObjectBase(internal::constructor_without_unaligned_array_assert)
: m_storage(internal::constructor_without_unaligned_array_assert())
{
// _check_template_params(); EIGEN_INITIALIZE_COEFFS_IF_THAT_OPTION_IS_ENABLED
}
#endif
#if EIGEN_HAS_RVALUE_REFERENCES
EIGEN_DEVICE_FUNC
PlainObjectBase(PlainObjectBase&& other) EIGEN_NOEXCEPT
: m_storage( std::move(other.m_storage) )
{
}
EIGEN_DEVICE_FUNC
PlainObjectBase& operator=(PlainObjectBase&& other) EIGEN_NOEXCEPT
{
using std::swap;
swap(m_storage, other.m_storage);
return *this;
}
#endif
/** Copy constructor */
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE PlainObjectBase(const PlainObjectBase& other)
: Base(), m_storage(other.m_storage) { }
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE PlainObjectBase(Index size, Index rows, Index cols)
: m_storage(size, rows, cols)
{
// _check_template_params();
// EIGEN_INITIALIZE_COEFFS_IF_THAT_OPTION_IS_ENABLED
}
/** \sa PlainObjectBase::operator=(const EigenBase<OtherDerived>&) */
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE PlainObjectBase(const DenseBase<OtherDerived> &other)
: m_storage()
{
_check_template_params();
resizeLike(other);
_set_noalias(other);
}
/** \sa PlainObjectBase::operator=(const EigenBase<OtherDerived>&) */
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE PlainObjectBase(const EigenBase<OtherDerived> &other)
: m_storage()
{
_check_template_params();
resizeLike(other);
*this = other.derived();
}
/** \brief Copy constructor with in-place evaluation */
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE PlainObjectBase(const ReturnByValue<OtherDerived>& other)
{
_check_template_params();
// FIXME this does not automatically transpose vectors if necessary
resize(other.rows(), other.cols());
other.evalTo(this->derived());
}
public:
/** \brief Copies the generic expression \a other into *this.
* \copydetails DenseBase::operator=(const EigenBase<OtherDerived> &other)
*/
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE Derived& operator=(const EigenBase<OtherDerived> &other)
{
_resize_to_match(other);
Base::operator=(other.derived());
return this->derived();
}
/** \name Map
* These are convenience functions returning Map objects. The Map() static functions return unaligned Map objects,
* while the AlignedMap() functions return aligned Map objects and thus should be called only with 16-byte-aligned
* \a data pointers.
*
* \see class Map
*/
//@{
static inline ConstMapType Map(const Scalar* data)
{ return ConstMapType(data); }
static inline MapType Map(Scalar* data)
{ return MapType(data); }
static inline ConstMapType Map(const Scalar* data, Index size)
{ return ConstMapType(data, size); }
static inline MapType Map(Scalar* data, Index size)
{ return MapType(data, size); }
static inline ConstMapType Map(const Scalar* data, Index rows, Index cols)
{ return ConstMapType(data, rows, cols); }
static inline MapType Map(Scalar* data, Index rows, Index cols)
{ return MapType(data, rows, cols); }
static inline ConstAlignedMapType MapAligned(const Scalar* data)
{ return ConstAlignedMapType(data); }
static inline AlignedMapType MapAligned(Scalar* data)
{ return AlignedMapType(data); }
static inline ConstAlignedMapType MapAligned(const Scalar* data, Index size)
{ return ConstAlignedMapType(data, size); }
static inline AlignedMapType MapAligned(Scalar* data, Index size)
{ return AlignedMapType(data, size); }
static inline ConstAlignedMapType MapAligned(const Scalar* data, Index rows, Index cols)
{ return ConstAlignedMapType(data, rows, cols); }
static inline AlignedMapType MapAligned(Scalar* data, Index rows, Index cols)
{ return AlignedMapType(data, rows, cols); }
template<int Outer, int Inner>
static inline typename StridedConstMapType<Stride<Outer, Inner> >::type Map(const Scalar* data, const Stride<Outer, Inner>& stride)
{ return typename StridedConstMapType<Stride<Outer, Inner> >::type(data, stride); }
template<int Outer, int Inner>
static inline typename StridedMapType<Stride<Outer, Inner> >::type Map(Scalar* data, const Stride<Outer, Inner>& stride)
{ return typename StridedMapType<Stride<Outer, Inner> >::type(data, stride); }
template<int Outer, int Inner>
static inline typename StridedConstMapType<Stride<Outer, Inner> >::type Map(const Scalar* data, Index size, const Stride<Outer, Inner>& stride)
{ return typename StridedConstMapType<Stride<Outer, Inner> >::type(data, size, stride); }
template<int Outer, int Inner>
static inline typename StridedMapType<Stride<Outer, Inner> >::type Map(Scalar* data, Index size, const Stride<Outer, Inner>& stride)
{ return typename StridedMapType<Stride<Outer, Inner> >::type(data, size, stride); }
template<int Outer, int Inner>
static inline typename StridedConstMapType<Stride<Outer, Inner> >::type Map(const Scalar* data, Index rows, Index cols, const Stride<Outer, Inner>& stride)
{ return typename StridedConstMapType<Stride<Outer, Inner> >::type(data, rows, cols, stride); }
template<int Outer, int Inner>
static inline typename StridedMapType<Stride<Outer, Inner> >::type Map(Scalar* data, Index rows, Index cols, const Stride<Outer, Inner>& stride)
{ return typename StridedMapType<Stride<Outer, Inner> >::type(data, rows, cols, stride); }
template<int Outer, int Inner>
static inline typename StridedConstAlignedMapType<Stride<Outer, Inner> >::type MapAligned(const Scalar* data, const Stride<Outer, Inner>& stride)
{ return typename StridedConstAlignedMapType<Stride<Outer, Inner> >::type(data, stride); }
template<int Outer, int Inner>
static inline typename StridedAlignedMapType<Stride<Outer, Inner> >::type MapAligned(Scalar* data, const Stride<Outer, Inner>& stride)
{ return typename StridedAlignedMapType<Stride<Outer, Inner> >::type(data, stride); }
template<int Outer, int Inner>
static inline typename StridedConstAlignedMapType<Stride<Outer, Inner> >::type MapAligned(const Scalar* data, Index size, const Stride<Outer, Inner>& stride)
{ return typename StridedConstAlignedMapType<Stride<Outer, Inner> >::type(data, size, stride); }
template<int Outer, int Inner>
static inline typename StridedAlignedMapType<Stride<Outer, Inner> >::type MapAligned(Scalar* data, Index size, const Stride<Outer, Inner>& stride)
{ return typename StridedAlignedMapType<Stride<Outer, Inner> >::type(data, size, stride); }
template<int Outer, int Inner>
static inline typename StridedConstAlignedMapType<Stride<Outer, Inner> >::type MapAligned(const Scalar* data, Index rows, Index cols, const Stride<Outer, Inner>& stride)
{ return typename StridedConstAlignedMapType<Stride<Outer, Inner> >::type(data, rows, cols, stride); }
template<int Outer, int Inner>
static inline typename StridedAlignedMapType<Stride<Outer, Inner> >::type MapAligned(Scalar* data, Index rows, Index cols, const Stride<Outer, Inner>& stride)
{ return typename StridedAlignedMapType<Stride<Outer, Inner> >::type(data, rows, cols, stride); }
//@}
using Base::setConstant;
EIGEN_DEVICE_FUNC Derived& setConstant(Index size, const Scalar& val);
EIGEN_DEVICE_FUNC Derived& setConstant(Index rows, Index cols, const Scalar& val);
using Base::setZero;
EIGEN_DEVICE_FUNC Derived& setZero(Index size);
EIGEN_DEVICE_FUNC Derived& setZero(Index rows, Index cols);
using Base::setOnes;
EIGEN_DEVICE_FUNC Derived& setOnes(Index size);
EIGEN_DEVICE_FUNC Derived& setOnes(Index rows, Index cols);
using Base::setRandom;
Derived& setRandom(Index size);
Derived& setRandom(Index rows, Index cols);
#ifdef EIGEN_PLAINOBJECTBASE_PLUGIN
#include EIGEN_PLAINOBJECTBASE_PLUGIN
#endif
protected:
/** \internal Resizes *this in preparation for assigning \a other to it.
* Takes care of doing all the checking that's needed.
*
* Note that copying a row-vector into a vector (and conversely) is allowed.
* The resizing, if any, is then done in the appropriate way so that row-vectors
* remain row-vectors and vectors remain vectors.
*/
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE void _resize_to_match(const EigenBase<OtherDerived>& other)
{
#ifdef EIGEN_NO_AUTOMATIC_RESIZING
eigen_assert((this->size()==0 || (IsVectorAtCompileTime ? (this->size() == other.size())
: (rows() == other.rows() && cols() == other.cols())))
&& "Size mismatch. Automatic resizing is disabled because EIGEN_NO_AUTOMATIC_RESIZING is defined");
EIGEN_ONLY_USED_FOR_DEBUG(other);
#else
resizeLike(other);
#endif
}
/**
* \brief Copies the value of the expression \a other into \c *this with automatic resizing.
*
* *this might be resized to match the dimensions of \a other. If *this was a null matrix (not already initialized),
* it will be initialized.
*
* Note that copying a row-vector into a vector (and conversely) is allowed.
* The resizing, if any, is then done in the appropriate way so that row-vectors
* remain row-vectors and vectors remain vectors.
*
* \sa operator=(const MatrixBase<OtherDerived>&), _set_noalias()
*
* \internal
*/
// aliasing is dealt once in internall::call_assignment
// so at this stage we have to assume aliasing... and resising has to be done later.
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE Derived& _set(const DenseBase<OtherDerived>& other)
{
internal::call_assignment(this->derived(), other.derived());
return this->derived();
}
/** \internal Like _set() but additionally makes the assumption that no aliasing effect can happen (which
* is the case when creating a new matrix) so one can enforce lazy evaluation.
*
* \sa operator=(const MatrixBase<OtherDerived>&), _set()
*/
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE Derived& _set_noalias(const DenseBase<OtherDerived>& other)
{
// I don't think we need this resize call since the lazyAssign will anyways resize
// and lazyAssign will be called by the assign selector.
//_resize_to_match(other);
// the 'false' below means to enforce lazy evaluation. We don't use lazyAssign() because
// it wouldn't allow to copy a row-vector into a column-vector.
internal::call_assignment_no_alias(this->derived(), other.derived(), internal::assign_op<Scalar,typename OtherDerived::Scalar>());
return this->derived();
}
template<typename T0, typename T1>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE void _init2(Index rows, Index cols, typename internal::enable_if<Base::SizeAtCompileTime!=2,T0>::type* = 0)
{
EIGEN_STATIC_ASSERT(bool(NumTraits<T0>::IsInteger) &&
bool(NumTraits<T1>::IsInteger),
FLOATING_POINT_ARGUMENT_PASSED__INTEGER_WAS_EXPECTED)
resize(rows,cols);
}
template<typename T0, typename T1>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE void _init2(const T0& val0, const T1& val1, typename internal::enable_if<Base::SizeAtCompileTime==2,T0>::type* = 0)
{
EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(PlainObjectBase, 2)
m_storage.data()[0] = Scalar(val0);
m_storage.data()[1] = Scalar(val1);
}
template<typename T0, typename T1>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE void _init2(const Index& val0, const Index& val1,
typename internal::enable_if< (!internal::is_same<Index,Scalar>::value)
&& (internal::is_same<T0,Index>::value)
&& (internal::is_same<T1,Index>::value)
&& Base::SizeAtCompileTime==2,T1>::type* = 0)
{
EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(PlainObjectBase, 2)
m_storage.data()[0] = Scalar(val0);
m_storage.data()[1] = Scalar(val1);
}
// The argument is convertible to the Index type and we either have a non 1x1 Matrix, or a dynamic-sized Array,
// then the argument is meant to be the size of the object.
template<typename T>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE void _init1(Index size, typename internal::enable_if< (Base::SizeAtCompileTime!=1 || !internal::is_convertible<T, Scalar>::value)
&& ((!internal::is_same<typename internal::traits<Derived>::XprKind,ArrayXpr>::value || Base::SizeAtCompileTime==Dynamic)),T>::type* = 0)
{
// NOTE MSVC 2008 complains if we directly put bool(NumTraits<T>::IsInteger) as the EIGEN_STATIC_ASSERT argument.
const bool is_integer = NumTraits<T>::IsInteger;
EIGEN_UNUSED_VARIABLE(is_integer);
EIGEN_STATIC_ASSERT(is_integer,
FLOATING_POINT_ARGUMENT_PASSED__INTEGER_WAS_EXPECTED)
resize(size);
}
// We have a 1x1 matrix/array => the argument is interpreted as the value of the unique coefficient (case where scalar type can be implicitely converted)
template<typename T>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE void _init1(const Scalar& val0, typename internal::enable_if<Base::SizeAtCompileTime==1 && internal::is_convertible<T, Scalar>::value,T>::type* = 0)
{
EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(PlainObjectBase, 1)
m_storage.data()[0] = val0;
}
// We have a 1x1 matrix/array => the argument is interpreted as the value of the unique coefficient (case where scalar type match the index type)
template<typename T>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE void _init1(const Index& val0,
typename internal::enable_if< (!internal::is_same<Index,Scalar>::value)
&& (internal::is_same<Index,T>::value)
&& Base::SizeAtCompileTime==1
&& internal::is_convertible<T, Scalar>::value,T*>::type* = 0)
{
EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(PlainObjectBase, 1)
m_storage.data()[0] = Scalar(val0);
}
// Initialize a fixed size matrix from a pointer to raw data
template<typename T>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE void _init1(const Scalar* data){
this->_set_noalias(ConstMapType(data));
}
// Initialize an arbitrary matrix from a dense expression
template<typename T, typename OtherDerived>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE void _init1(const DenseBase<OtherDerived>& other){
this->_set_noalias(other);
}
// Initialize an arbitrary matrix from an object convertible to the Derived type.
template<typename T>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE void _init1(const Derived& other){
this->_set_noalias(other);
}
// Initialize an arbitrary matrix from a generic Eigen expression
template<typename T, typename OtherDerived>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE void _init1(const EigenBase<OtherDerived>& other){
this->derived() = other;
}
template<typename T, typename OtherDerived>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE void _init1(const ReturnByValue<OtherDerived>& other)
{
resize(other.rows(), other.cols());
other.evalTo(this->derived());
}
template<typename T, typename OtherDerived, int ColsAtCompileTime>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE void _init1(const RotationBase<OtherDerived,ColsAtCompileTime>& r)
{
this->derived() = r;
}
// For fixed-size Array<Scalar,...>
template<typename T>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE void _init1(const Scalar& val0,
typename internal::enable_if< Base::SizeAtCompileTime!=Dynamic
&& Base::SizeAtCompileTime!=1
&& internal::is_convertible<T, Scalar>::value
&& internal::is_same<typename internal::traits<Derived>::XprKind,ArrayXpr>::value,T>::type* = 0)
{
Base::setConstant(val0);
}
// For fixed-size Array<Index,...>
template<typename T>
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE void _init1(const Index& val0,
typename internal::enable_if< (!internal::is_same<Index,Scalar>::value)
&& (internal::is_same<Index,T>::value)
&& Base::SizeAtCompileTime!=Dynamic
&& Base::SizeAtCompileTime!=1
&& internal::is_convertible<T, Scalar>::value
&& internal::is_same<typename internal::traits<Derived>::XprKind,ArrayXpr>::value,T*>::type* = 0)
{
Base::setConstant(val0);
}
template<typename MatrixTypeA, typename MatrixTypeB, bool SwapPointers>
friend struct internal::matrix_swap_impl;
public:
#ifndef EIGEN_PARSED_BY_DOXYGEN
/** \internal
* \brief Override DenseBase::swap() since for dynamic-sized matrices
* of same type it is enough to swap the data pointers.
*/
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
void swap(DenseBase<OtherDerived> & other)
{
enum { SwapPointers = internal::is_same<Derived, OtherDerived>::value && Base::SizeAtCompileTime==Dynamic };
internal::matrix_swap_impl<Derived, OtherDerived, bool(SwapPointers)>::run(this->derived(), other.derived());
}
/** \internal
* \brief const version forwarded to DenseBase::swap
*/
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
void swap(DenseBase<OtherDerived> const & other)
{ Base::swap(other.derived()); }
EIGEN_DEVICE_FUNC
static EIGEN_STRONG_INLINE void _check_template_params()
{
EIGEN_STATIC_ASSERT((EIGEN_IMPLIES(MaxRowsAtCompileTime==1 && MaxColsAtCompileTime!=1, (Options&RowMajor)==RowMajor)
&& EIGEN_IMPLIES(MaxColsAtCompileTime==1 && MaxRowsAtCompileTime!=1, (Options&RowMajor)==0)
&& ((RowsAtCompileTime == Dynamic) || (RowsAtCompileTime >= 0))
&& ((ColsAtCompileTime == Dynamic) || (ColsAtCompileTime >= 0))
&& ((MaxRowsAtCompileTime == Dynamic) || (MaxRowsAtCompileTime >= 0))
&& ((MaxColsAtCompileTime == Dynamic) || (MaxColsAtCompileTime >= 0))
&& (MaxRowsAtCompileTime == RowsAtCompileTime || RowsAtCompileTime==Dynamic)
&& (MaxColsAtCompileTime == ColsAtCompileTime || ColsAtCompileTime==Dynamic)
&& (Options & (DontAlign|RowMajor)) == Options),
INVALID_MATRIX_TEMPLATE_PARAMETERS)
}
enum { IsPlainObjectBase = 1 };
#endif
};
namespace internal {
template <typename Derived, typename OtherDerived, bool IsVector>
struct conservative_resize_like_impl
{
static void run(DenseBase<Derived>& _this, Index rows, Index cols)
{
if (_this.rows() == rows && _this.cols() == cols) return;
EIGEN_STATIC_ASSERT_DYNAMIC_SIZE(Derived)
if ( ( Derived::IsRowMajor && _this.cols() == cols) || // row-major and we change only the number of rows
(!Derived::IsRowMajor && _this.rows() == rows) ) // column-major and we change only the number of columns
{
internal::check_rows_cols_for_overflow<Derived::MaxSizeAtCompileTime>::run(rows, cols);
_this.derived().m_storage.conservativeResize(rows*cols,rows,cols);
}
else
{
// The storage order does not allow us to use reallocation.
typename Derived::PlainObject tmp(rows,cols);
const Index common_rows = numext::mini(rows, _this.rows());
const Index common_cols = numext::mini(cols, _this.cols());
tmp.block(0,0,common_rows,common_cols) = _this.block(0,0,common_rows,common_cols);
_this.derived().swap(tmp);
}
}
static void run(DenseBase<Derived>& _this, const DenseBase<OtherDerived>& other)
{
if (_this.rows() == other.rows() && _this.cols() == other.cols()) return;
// Note: Here is space for improvement. Basically, for conservativeResize(Index,Index),
// neither RowsAtCompileTime or ColsAtCompileTime must be Dynamic. If only one of the
// dimensions is dynamic, one could use either conservativeResize(Index rows, NoChange_t) or
// conservativeResize(NoChange_t, Index cols). For these methods new static asserts like
// EIGEN_STATIC_ASSERT_DYNAMIC_ROWS and EIGEN_STATIC_ASSERT_DYNAMIC_COLS would be good.
EIGEN_STATIC_ASSERT_DYNAMIC_SIZE(Derived)
EIGEN_STATIC_ASSERT_DYNAMIC_SIZE(OtherDerived)
if ( ( Derived::IsRowMajor && _this.cols() == other.cols()) || // row-major and we change only the number of rows
(!Derived::IsRowMajor && _this.rows() == other.rows()) ) // column-major and we change only the number of columns
{
const Index new_rows = other.rows() - _this.rows();
const Index new_cols = other.cols() - _this.cols();
_this.derived().m_storage.conservativeResize(other.size(),other.rows(),other.cols());
if (new_rows>0)
_this.bottomRightCorner(new_rows, other.cols()) = other.bottomRows(new_rows);
else if (new_cols>0)
_this.bottomRightCorner(other.rows(), new_cols) = other.rightCols(new_cols);
}
else
{
// The storage order does not allow us to use reallocation.
typename Derived::PlainObject tmp(other);
const Index common_rows = numext::mini(tmp.rows(), _this.rows());
const Index common_cols = numext::mini(tmp.cols(), _this.cols());
tmp.block(0,0,common_rows,common_cols) = _this.block(0,0,common_rows,common_cols);
_this.derived().swap(tmp);
}
}
};
// Here, the specialization for vectors inherits from the general matrix case
// to allow calling .conservativeResize(rows,cols) on vectors.
template <typename Derived, typename OtherDerived>
struct conservative_resize_like_impl<Derived,OtherDerived,true>
: conservative_resize_like_impl<Derived,OtherDerived,false>
{
using conservative_resize_like_impl<Derived,OtherDerived,false>::run;
static void run(DenseBase<Derived>& _this, Index size)
{
const Index new_rows = Derived::RowsAtCompileTime==1 ? 1 : size;
const Index new_cols = Derived::RowsAtCompileTime==1 ? size : 1;
_this.derived().m_storage.conservativeResize(size,new_rows,new_cols);
}
static void run(DenseBase<Derived>& _this, const DenseBase<OtherDerived>& other)
{
if (_this.rows() == other.rows() && _this.cols() == other.cols()) return;
const Index num_new_elements = other.size() - _this.size();
const Index new_rows = Derived::RowsAtCompileTime==1 ? 1 : other.rows();
const Index new_cols = Derived::RowsAtCompileTime==1 ? other.cols() : 1;
_this.derived().m_storage.conservativeResize(other.size(),new_rows,new_cols);
if (num_new_elements > 0)
_this.tail(num_new_elements) = other.tail(num_new_elements);
}
};
template<typename MatrixTypeA, typename MatrixTypeB, bool SwapPointers>
struct matrix_swap_impl
{
EIGEN_DEVICE_FUNC
static inline void run(MatrixTypeA& a, MatrixTypeB& b)
{
a.base().swap(b);
}
};
template<typename MatrixTypeA, typename MatrixTypeB>
struct matrix_swap_impl<MatrixTypeA, MatrixTypeB, true>
{
EIGEN_DEVICE_FUNC
static inline void run(MatrixTypeA& a, MatrixTypeB& b)
{
static_cast<typename MatrixTypeA::Base&>(a).m_storage.swap(static_cast<typename MatrixTypeB::Base&>(b).m_storage);
}
};
} // end namespace internal
} // end namespace Eigen
#endif // EIGEN_DENSESTORAGEBASE_H
external/eigen3/Eigen/src/Core/Product.h 0 → 100644
View file @ a394b22a
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_PRODUCT_H
#define EIGEN_PRODUCT_H
namespace Eigen {
template<typename Lhs, typename Rhs, int Option, typename StorageKind> class ProductImpl;
namespace internal {
template<typename Lhs, typename Rhs, int Option>
struct traits<Product<Lhs, Rhs, Option> >
{
typedef typename remove_all<Lhs>::type LhsCleaned;
typedef typename remove_all<Rhs>::type RhsCleaned;
typedef traits<LhsCleaned> LhsTraits;
typedef traits<RhsCleaned> RhsTraits;
typedef MatrixXpr XprKind;
typedef typename ScalarBinaryOpTraits<typename traits<LhsCleaned>::Scalar, typename traits<RhsCleaned>::Scalar>::ReturnType Scalar;
typedef typename product_promote_storage_type<typename LhsTraits::StorageKind,
typename RhsTraits::StorageKind,
internal::product_type<Lhs,Rhs>::ret>::ret StorageKind;
typedef typename promote_index_type<typename LhsTraits::StorageIndex,
typename RhsTraits::StorageIndex>::type StorageIndex;
enum {
RowsAtCompileTime = LhsTraits::RowsAtCompileTime,
ColsAtCompileTime = RhsTraits::ColsAtCompileTime,
MaxRowsAtCompileTime = LhsTraits::MaxRowsAtCompileTime,
MaxColsAtCompileTime = RhsTraits::MaxColsAtCompileTime,
// FIXME: only needed by GeneralMatrixMatrixTriangular
InnerSize = EIGEN_SIZE_MIN_PREFER_FIXED(LhsTraits::ColsAtCompileTime, RhsTraits::RowsAtCompileTime),
// The storage order is somewhat arbitrary here. The correct one will be determined through the evaluator.
Flags = (MaxRowsAtCompileTime==1 && MaxColsAtCompileTime!=1) ? RowMajorBit
: (MaxColsAtCompileTime==1 && MaxRowsAtCompileTime!=1) ? 0
: ( ((LhsTraits::Flags&NoPreferredStorageOrderBit) && (RhsTraits::Flags&RowMajorBit))
|| ((RhsTraits::Flags&NoPreferredStorageOrderBit) && (LhsTraits::Flags&RowMajorBit)) ) ? RowMajorBit
: NoPreferredStorageOrderBit
};
};
} // end namespace internal
/** \class Product
* \ingroup Core_Module
*
* \brief Expression of the product of two arbitrary matrices or vectors
*
* \tparam _Lhs the type of the left-hand side expression
* \tparam _Rhs the type of the right-hand side expression
*
* This class represents an expression of the product of two arbitrary matrices.
*
* The other template parameters are:
* \tparam Option can be DefaultProduct, AliasFreeProduct, or LazyProduct
*
*/
template<typename _Lhs, typename _Rhs, int Option>
class Product : public ProductImpl<_Lhs,_Rhs,Option,
typename internal::product_promote_storage_type<typename internal::traits<_Lhs>::StorageKind,
typename internal::traits<_Rhs>::StorageKind,
internal::product_type<_Lhs,_Rhs>::ret>::ret>
{
public:
typedef _Lhs Lhs;
typedef _Rhs Rhs;
typedef typename ProductImpl<
Lhs, Rhs, Option,
typename internal::product_promote_storage_type<typename internal::traits<Lhs>::StorageKind,
typename internal::traits<Rhs>::StorageKind,
internal::product_type<Lhs,Rhs>::ret>::ret>::Base Base;
EIGEN_GENERIC_PUBLIC_INTERFACE(Product)
typedef typename internal::ref_selector<Lhs>::type LhsNested;
typedef typename internal::ref_selector<Rhs>::type RhsNested;
typedef typename internal::remove_all<LhsNested>::type LhsNestedCleaned;
typedef typename internal::remove_all<RhsNested>::type RhsNestedCleaned;
EIGEN_DEVICE_FUNC Product(const Lhs& lhs, const Rhs& rhs) : m_lhs(lhs), m_rhs(rhs)
{
eigen_assert(lhs.cols() == rhs.rows()
&& "invalid matrix product"
&& "if you wanted a coeff-wise or a dot product use the respective explicit functions");
}
EIGEN_DEVICE_FUNC inline Index rows() const { return m_lhs.rows(); }
EIGEN_DEVICE_FUNC inline Index cols() const { return m_rhs.cols(); }
EIGEN_DEVICE_FUNC const LhsNestedCleaned& lhs() const { return m_lhs; }
EIGEN_DEVICE_FUNC const RhsNestedCleaned& rhs() const { return m_rhs; }
protected:
LhsNested m_lhs;
RhsNested m_rhs;
};
namespace internal {
template<typename Lhs, typename Rhs, int Option, int ProductTag = internal::product_type<Lhs,Rhs>::ret>
class dense_product_base
: public internal::dense_xpr_base<Product<Lhs,Rhs,Option> >::type
{};
/** Convertion to scalar for inner-products */
template<typename Lhs, typename Rhs, int Option>
class dense_product_base<Lhs, Rhs, Option, InnerProduct>
: public internal::dense_xpr_base<Product<Lhs,Rhs,Option> >::type
{
typedef Product<Lhs,Rhs,Option> ProductXpr;
typedef typename internal::dense_xpr_base<ProductXpr>::type Base;
public:
using Base::derived;
typedef typename Base::Scalar Scalar;
operator const Scalar() const
{
return internal::evaluator<ProductXpr>(derived()).coeff(0,0);
}
};
} // namespace internal
// Generic API dispatcher
template<typename Lhs, typename Rhs, int Option, typename StorageKind>
class ProductImpl : public internal::generic_xpr_base<Product<Lhs,Rhs,Option>, MatrixXpr, StorageKind>::type
{
public:
typedef typename internal::generic_xpr_base<Product<Lhs,Rhs,Option>, MatrixXpr, StorageKind>::type Base;
};
template<typename Lhs, typename Rhs, int Option>
class ProductImpl<Lhs,Rhs,Option,Dense>
: public internal::dense_product_base<Lhs,Rhs,Option>
{
typedef Product<Lhs, Rhs, Option> Derived;
public:
typedef typename internal::dense_product_base<Lhs, Rhs, Option> Base;
EIGEN_DENSE_PUBLIC_INTERFACE(Derived)
protected:
enum {
IsOneByOne = (RowsAtCompileTime == 1 || RowsAtCompileTime == Dynamic) &&
(ColsAtCompileTime == 1 || ColsAtCompileTime == Dynamic),
EnableCoeff = IsOneByOne || Option==LazyProduct
};
public:
EIGEN_DEVICE_FUNC Scalar coeff(Index row, Index col) const
{
EIGEN_STATIC_ASSERT(EnableCoeff, THIS_METHOD_IS_ONLY_FOR_INNER_OR_LAZY_PRODUCTS);
eigen_assert( (Option==LazyProduct) || (this->rows() == 1 && this->cols() == 1) );
return internal::evaluator<Derived>(derived()).coeff(row,col);
}
EIGEN_DEVICE_FUNC Scalar coeff(Index i) const
{
EIGEN_STATIC_ASSERT(EnableCoeff, THIS_METHOD_IS_ONLY_FOR_INNER_OR_LAZY_PRODUCTS);
eigen_assert( (Option==LazyProduct) || (this->rows() == 1 && this->cols() == 1) );
return internal::evaluator<Derived>(derived()).coeff(i);
}
};
} // end namespace Eigen
#endif // EIGEN_PRODUCT_H
external/eigen3/Eigen/src/Core/ProductEvaluators.h 0 → 100644
View file @ a394b22a
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2011 Jitse Niesen <jitse@maths.leeds.ac.uk>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_PRODUCTEVALUATORS_H
#define EIGEN_PRODUCTEVALUATORS_H
namespace Eigen {
namespace internal {
/** \internal
* Evaluator of a product expression.
* Since products require special treatments to handle all possible cases,
* we simply deffer the evaluation logic to a product_evaluator class
* which offers more partial specialization possibilities.
*
* \sa class product_evaluator
*/
template<typename Lhs, typename Rhs, int Options>
struct evaluator<Product<Lhs, Rhs, Options> >
: public product_evaluator<Product<Lhs, Rhs, Options> >
{
typedef Product<Lhs, Rhs, Options> XprType;
typedef product_evaluator<XprType> Base;
EIGEN_DEVICE_FUNC explicit evaluator(const XprType& xpr) : Base(xpr) {}
};
// Catch "scalar * ( A * B )" and transform it to "(A*scalar) * B"
// TODO we should apply that rule only if that's really helpful
template<typename Lhs, typename Rhs, typename Scalar1, typename Scalar2, typename Plain1>
struct evaluator_assume_aliasing<CwiseBinaryOp<internal::scalar_product_op<Scalar1,Scalar2>,
const CwiseNullaryOp<internal::scalar_constant_op<Scalar1>, Plain1>,
const Product<Lhs, Rhs, DefaultProduct> > >
{
static const bool value = true;
};
template<typename Lhs, typename Rhs, typename Scalar1, typename Scalar2, typename Plain1>
struct evaluator<CwiseBinaryOp<internal::scalar_product_op<Scalar1,Scalar2>,
const CwiseNullaryOp<internal::scalar_constant_op<Scalar1>, Plain1>,
const Product<Lhs, Rhs, DefaultProduct> > >
: public evaluator<Product<EIGEN_SCALAR_BINARYOP_EXPR_RETURN_TYPE(Scalar1,Lhs,product), Rhs, DefaultProduct> >
{
typedef CwiseBinaryOp<internal::scalar_product_op<Scalar1,Scalar2>,
const CwiseNullaryOp<internal::scalar_constant_op<Scalar1>, Plain1>,
const Product<Lhs, Rhs, DefaultProduct> > XprType;
typedef evaluator<Product<EIGEN_SCALAR_BINARYOP_EXPR_RETURN_TYPE(Scalar1,Lhs,product), Rhs, DefaultProduct> > Base;
EIGEN_DEVICE_FUNC explicit evaluator(const XprType& xpr)
: Base(xpr.lhs().functor().m_other * xpr.rhs().lhs() * xpr.rhs().rhs())
{}
};
template<typename Lhs, typename Rhs, int DiagIndex>
struct evaluator<Diagonal<const Product<Lhs, Rhs, DefaultProduct>, DiagIndex> >
: public evaluator<Diagonal<const Product<Lhs, Rhs, LazyProduct>, DiagIndex> >
{
typedef Diagonal<const Product<Lhs, Rhs, DefaultProduct>, DiagIndex> XprType;
typedef evaluator<Diagonal<const Product<Lhs, Rhs, LazyProduct>, DiagIndex> > Base;
EIGEN_DEVICE_FUNC explicit evaluator(const XprType& xpr)
: Base(Diagonal<const Product<Lhs, Rhs, LazyProduct>, DiagIndex>(
Product<Lhs, Rhs, LazyProduct>(xpr.nestedExpression().lhs(), xpr.nestedExpression().rhs()),
xpr.index() ))
{}
};
// Helper class to perform a matrix product with the destination at hand.
// Depending on the sizes of the factors, there are different evaluation strategies
// as controlled by internal::product_type.
template< typename Lhs, typename Rhs,
typename LhsShape = typename evaluator_traits<Lhs>::Shape,
typename RhsShape = typename evaluator_traits<Rhs>::Shape,
int ProductType = internal::product_type<Lhs,Rhs>::value>
struct generic_product_impl;
template<typename Lhs, typename Rhs>
struct evaluator_assume_aliasing<Product<Lhs, Rhs, DefaultProduct> > {
static const bool value = true;
};
// This is the default evaluator implementation for products:
// It creates a temporary and call generic_product_impl
template<typename Lhs, typename Rhs, int Options, int ProductTag, typename LhsShape, typename RhsShape>
struct product_evaluator<Product<Lhs, Rhs, Options>, ProductTag, LhsShape, RhsShape>
: public evaluator<typename Product<Lhs, Rhs, Options>::PlainObject>
{
typedef Product<Lhs, Rhs, Options> XprType;
typedef typename XprType::PlainObject PlainObject;
typedef evaluator<PlainObject> Base;
enum {
Flags = Base::Flags | EvalBeforeNestingBit
};
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
explicit product_evaluator(const XprType& xpr)
: m_result(xpr.rows(), xpr.cols())
{
::new (static_cast<Base*>(this)) Base(m_result);
// FIXME shall we handle nested_eval here?,
// if so, then we must take care at removing the call to nested_eval in the specializations (e.g., in permutation_matrix_product, transposition_matrix_product, etc.)
// typedef typename internal::nested_eval<Lhs,Rhs::ColsAtCompileTime>::type LhsNested;
// typedef typename internal::nested_eval<Rhs,Lhs::RowsAtCompileTime>::type RhsNested;
// typedef typename internal::remove_all<LhsNested>::type LhsNestedCleaned;
// typedef typename internal::remove_all<RhsNested>::type RhsNestedCleaned;
//
// const LhsNested lhs(xpr.lhs());
// const RhsNested rhs(xpr.rhs());
//
// generic_product_impl<LhsNestedCleaned, RhsNestedCleaned>::evalTo(m_result, lhs, rhs);
generic_product_impl<Lhs, Rhs, LhsShape, RhsShape, ProductTag>::evalTo(m_result, xpr.lhs(), xpr.rhs());
}
protected:
PlainObject m_result;
};
// The following three shortcuts are enabled only if the scalar types match excatly.
// TODO: we could enable them for different scalar types when the product is not vectorized.
// Dense = Product
template< typename DstXprType, typename Lhs, typename Rhs, int Options, typename Scalar>
struct Assignment<DstXprType, Product<Lhs,Rhs,Options>, internal::assign_op<Scalar,Scalar>, Dense2Dense,
typename enable_if<(Options==DefaultProduct || Options==AliasFreeProduct)>::type>
{
typedef Product<Lhs,Rhs,Options> SrcXprType;
static EIGEN_STRONG_INLINE
void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &)
{
Index dstRows = src.rows();
Index dstCols = src.cols();
if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
dst.resize(dstRows, dstCols);
// FIXME shall we handle nested_eval here?
generic_product_impl<Lhs, Rhs>::evalTo(dst, src.lhs(), src.rhs());
}
};
// Dense += Product
template< typename DstXprType, typename Lhs, typename Rhs, int Options, typename Scalar>
struct Assignment<DstXprType, Product<Lhs,Rhs,Options>, internal::add_assign_op<Scalar,Scalar>, Dense2Dense,
typename enable_if<(Options==DefaultProduct || Options==AliasFreeProduct)>::type>
{
typedef Product<Lhs,Rhs,Options> SrcXprType;
static EIGEN_STRONG_INLINE
void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op<Scalar,Scalar> &)
{
eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols());
// FIXME shall we handle nested_eval here?
generic_product_impl<Lhs, Rhs>::addTo(dst, src.lhs(), src.rhs());
}
};
// Dense -= Product
template< typename DstXprType, typename Lhs, typename Rhs, int Options, typename Scalar>
struct Assignment<DstXprType, Product<Lhs,Rhs,Options>, internal::sub_assign_op<Scalar,Scalar>, Dense2Dense,
typename enable_if<(Options==DefaultProduct || Options==AliasFreeProduct)>::type>
{
typedef Product<Lhs,Rhs,Options> SrcXprType;
static EIGEN_STRONG_INLINE
void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op<Scalar,Scalar> &)
{
eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols());
// FIXME shall we handle nested_eval here?
generic_product_impl<Lhs, Rhs>::subTo(dst, src.lhs(), src.rhs());
}
};
// Dense ?= scalar * Product
// TODO we should apply that rule if that's really helpful
// for instance, this is not good for inner products
template< typename DstXprType, typename Lhs, typename Rhs, typename AssignFunc, typename Scalar, typename ScalarBis, typename Plain>
struct Assignment<DstXprType, CwiseBinaryOp<internal::scalar_product_op<ScalarBis,Scalar>, const CwiseNullaryOp<internal::scalar_constant_op<ScalarBis>,Plain>,
const Product<Lhs,Rhs,DefaultProduct> >, AssignFunc, Dense2Dense>
{
typedef CwiseBinaryOp<internal::scalar_product_op<ScalarBis,Scalar>,
const CwiseNullaryOp<internal::scalar_constant_op<ScalarBis>,Plain>,
const Product<Lhs,Rhs,DefaultProduct> > SrcXprType;
static EIGEN_STRONG_INLINE
void run(DstXprType &dst, const SrcXprType &src, const AssignFunc& func)
{
call_assignment_no_alias(dst, (src.lhs().functor().m_other * src.rhs().lhs())*src.rhs().rhs(), func);
}
};
//----------------------------------------
// Catch "Dense ?= xpr + Product<>" expression to save one temporary
// FIXME we could probably enable these rules for any product, i.e., not only Dense and DefaultProduct
template<typename OtherXpr, typename Lhs, typename Rhs>
struct evaluator_assume_aliasing<CwiseBinaryOp<internal::scalar_sum_op<typename OtherXpr::Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, const OtherXpr,
const Product<Lhs,Rhs,DefaultProduct> >, DenseShape > {
static const bool value = true;
};
template<typename OtherXpr, typename Lhs, typename Rhs>
struct evaluator_assume_aliasing<CwiseBinaryOp<internal::scalar_difference_op<typename OtherXpr::Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, const OtherXpr,
const Product<Lhs,Rhs,DefaultProduct> >, DenseShape > {
static const bool value = true;
};
template<typename DstXprType, typename OtherXpr, typename ProductType, typename Func1, typename Func2>
struct assignment_from_xpr_op_product
{
template<typename SrcXprType, typename InitialFunc>
static EIGEN_STRONG_INLINE
void run(DstXprType &dst, const SrcXprType &src, const InitialFunc& /*func*/)
{
call_assignment_no_alias(dst, src.lhs(), Func1());
call_assignment_no_alias(dst, src.rhs(), Func2());
}
};
#define EIGEN_CATCH_ASSIGN_XPR_OP_PRODUCT(ASSIGN_OP,BINOP,ASSIGN_OP2) \
template< typename DstXprType, typename OtherXpr, typename Lhs, typename Rhs, typename DstScalar, typename SrcScalar, typename OtherScalar,typename ProdScalar> \
struct Assignment<DstXprType, CwiseBinaryOp<internal::BINOP<OtherScalar,ProdScalar>, const OtherXpr, \
const Product<Lhs,Rhs,DefaultProduct> >, internal::ASSIGN_OP<DstScalar,SrcScalar>, Dense2Dense> \
: assignment_from_xpr_op_product<DstXprType, OtherXpr, Product<Lhs,Rhs,DefaultProduct>, internal::ASSIGN_OP<DstScalar,OtherScalar>, internal::ASSIGN_OP2<DstScalar,ProdScalar> > \
{}
EIGEN_CATCH_ASSIGN_XPR_OP_PRODUCT(assign_op, scalar_sum_op,add_assign_op);
EIGEN_CATCH_ASSIGN_XPR_OP_PRODUCT(add_assign_op,scalar_sum_op,add_assign_op);
EIGEN_CATCH_ASSIGN_XPR_OP_PRODUCT(sub_assign_op,scalar_sum_op,sub_assign_op);
EIGEN_CATCH_ASSIGN_XPR_OP_PRODUCT(assign_op, scalar_difference_op,sub_assign_op);
EIGEN_CATCH_ASSIGN_XPR_OP_PRODUCT(add_assign_op,scalar_difference_op,sub_assign_op);
EIGEN_CATCH_ASSIGN_XPR_OP_PRODUCT(sub_assign_op,scalar_difference_op,add_assign_op);
//----------------------------------------
template<typename Lhs, typename Rhs>
struct generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,InnerProduct>
{
template<typename Dst>
static inline void evalTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
{
dst.coeffRef(0,0) = (lhs.transpose().cwiseProduct(rhs)).sum();
}
template<typename Dst>
static inline void addTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
{
dst.coeffRef(0,0) += (lhs.transpose().cwiseProduct(rhs)).sum();
}
template<typename Dst>
static void subTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
{ dst.coeffRef(0,0) -= (lhs.transpose().cwiseProduct(rhs)).sum(); }
};
/***********************************************************************
* Implementation of outer dense * dense vector product
***********************************************************************/
// Column major result
template<typename Dst, typename Lhs, typename Rhs, typename Func>
void outer_product_selector_run(Dst& dst, const Lhs &lhs, const Rhs &rhs, const Func& func, const false_type&)
{
evaluator<Rhs> rhsEval(rhs);
typename nested_eval<Lhs,Rhs::SizeAtCompileTime>::type actual_lhs(lhs);
// FIXME if cols is large enough, then it might be useful to make sure that lhs is sequentially stored
// FIXME not very good if rhs is real and lhs complex while alpha is real too
const Index cols = dst.cols();
for (Index j=0; j<cols; ++j)
func(dst.col(j), rhsEval.coeff(Index(0),j) * actual_lhs);
}
// Row major result
template<typename Dst, typename Lhs, typename Rhs, typename Func>
void outer_product_selector_run(Dst& dst, const Lhs &lhs, const Rhs &rhs, const Func& func, const true_type&)
{
evaluator<Lhs> lhsEval(lhs);
typename nested_eval<Rhs,Lhs::SizeAtCompileTime>::type actual_rhs(rhs);
// FIXME if rows is large enough, then it might be useful to make sure that rhs is sequentially stored
// FIXME not very good if lhs is real and rhs complex while alpha is real too
const Index rows = dst.rows();
for (Index i=0; i<rows; ++i)
func(dst.row(i), lhsEval.coeff(i,Index(0)) * actual_rhs);
}
template<typename Lhs, typename Rhs>
struct generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,OuterProduct>
{
template<typename T> struct is_row_major : internal::conditional<(int(T::Flags)&RowMajorBit), internal::true_type, internal::false_type>::type {};
typedef typename Product<Lhs,Rhs>::Scalar Scalar;
// TODO it would be nice to be able to exploit our *_assign_op functors for that purpose
struct set { template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() = src; } };
struct add { template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() += src; } };
struct sub { template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() -= src; } };
struct adds {
Scalar m_scale;
explicit adds(const Scalar& s) : m_scale(s) {}
template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const {
dst.const_cast_derived() += m_scale * src;
}
};
template<typename Dst>
static inline void evalTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
{
internal::outer_product_selector_run(dst, lhs, rhs, set(), is_row_major<Dst>());
}
template<typename Dst>
static inline void addTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
{
internal::outer_product_selector_run(dst, lhs, rhs, add(), is_row_major<Dst>());
}
template<typename Dst>
static inline void subTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
{
internal::outer_product_selector_run(dst, lhs, rhs, sub(), is_row_major<Dst>());
}
template<typename Dst>
static inline void scaleAndAddTo(Dst& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
{
internal::outer_product_selector_run(dst, lhs, rhs, adds(alpha), is_row_major<Dst>());
}
};
// This base class provides default implementations for evalTo, addTo, subTo, in terms of scaleAndAddTo
template<typename Lhs, typename Rhs, typename Derived>
struct generic_product_impl_base
{
typedef typename Product<Lhs,Rhs>::Scalar Scalar;
template<typename Dst>
static EIGEN_STRONG_INLINE void evalTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
{ dst.setZero(); scaleAndAddTo(dst, lhs, rhs, Scalar(1)); }
template<typename Dst>
static EIGEN_STRONG_INLINE void addTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
{ scaleAndAddTo(dst,lhs, rhs, Scalar(1)); }
template<typename Dst>
static EIGEN_STRONG_INLINE void subTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
{ scaleAndAddTo(dst, lhs, rhs, Scalar(-1)); }
template<typename Dst>
static EIGEN_STRONG_INLINE void scaleAndAddTo(Dst& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
{ Derived::scaleAndAddTo(dst,lhs,rhs,alpha); }
};
template<typename Lhs, typename Rhs>
struct generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,GemvProduct>
: generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,GemvProduct> >
{
typedef typename nested_eval<Lhs,1>::type LhsNested;
typedef typename nested_eval<Rhs,1>::type RhsNested;
typedef typename Product<Lhs,Rhs>::Scalar Scalar;
enum { Side = Lhs::IsVectorAtCompileTime ? OnTheLeft : OnTheRight };
typedef typename internal::remove_all<typename internal::conditional<int(Side)==OnTheRight,LhsNested,RhsNested>::type>::type MatrixType;
template<typename Dest>
static EIGEN_STRONG_INLINE void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
{
LhsNested actual_lhs(lhs);
RhsNested actual_rhs(rhs);
internal::gemv_dense_selector<Side,
(int(MatrixType::Flags)&RowMajorBit) ? RowMajor : ColMajor,
bool(internal::blas_traits<MatrixType>::HasUsableDirectAccess)
>::run(actual_lhs, actual_rhs, dst, alpha);
}
};
template<typename Lhs, typename Rhs>
struct generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,CoeffBasedProductMode>
{
typedef typename Product<Lhs,Rhs>::Scalar Scalar;
template<typename Dst>
static EIGEN_STRONG_INLINE void evalTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
{
// Same as: dst.noalias() = lhs.lazyProduct(rhs);
// but easier on the compiler side
call_assignment_no_alias(dst, lhs.lazyProduct(rhs), internal::assign_op<typename Dst::Scalar,Scalar>());
}
template<typename Dst>
static EIGEN_STRONG_INLINE void addTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
{
// dst.noalias() += lhs.lazyProduct(rhs);
call_assignment_no_alias(dst, lhs.lazyProduct(rhs), internal::add_assign_op<typename Dst::Scalar,Scalar>());
}
template<typename Dst>
static EIGEN_STRONG_INLINE void subTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
{
// dst.noalias() -= lhs.lazyProduct(rhs);
call_assignment_no_alias(dst, lhs.lazyProduct(rhs), internal::sub_assign_op<typename Dst::Scalar,Scalar>());
}
// template<typename Dst>
// static inline void scaleAndAddTo(Dst& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
// { dst.noalias() += alpha * lhs.lazyProduct(rhs); }
};
// This specialization enforces the use of a coefficient-based evaluation strategy
template<typename Lhs, typename Rhs>
struct generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,LazyCoeffBasedProductMode>
: generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,CoeffBasedProductMode> {};
// Case 2: Evaluate coeff by coeff
//
// This is mostly taken from CoeffBasedProduct.h
// The main difference is that we add an extra argument to the etor_product_*_impl::run() function
// for the inner dimension of the product, because evaluator object do not know their size.
template<int Traversal, int UnrollingIndex, typename Lhs, typename Rhs, typename RetScalar>
struct etor_product_coeff_impl;
template<int StorageOrder, int UnrollingIndex, typename Lhs, typename Rhs, typename Packet, int LoadMode>
struct etor_product_packet_impl;
template<typename Lhs, typename Rhs, int ProductTag>
struct product_evaluator<Product<Lhs, Rhs, LazyProduct>, ProductTag, DenseShape, DenseShape>
: evaluator_base<Product<Lhs, Rhs, LazyProduct> >
{
typedef Product<Lhs, Rhs, LazyProduct> XprType;
typedef typename XprType::Scalar Scalar;
typedef typename XprType::CoeffReturnType CoeffReturnType;
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
explicit product_evaluator(const XprType& xpr)
: m_lhs(xpr.lhs()),
m_rhs(xpr.rhs()),
m_lhsImpl(m_lhs), // FIXME the creation of the evaluator objects should result in a no-op, but check that!
m_rhsImpl(m_rhs), // Moreover, they are only useful for the packet path, so we could completely disable them when not needed,
// or perhaps declare them on the fly on the packet method... We have experiment to check what's best.
m_innerDim(xpr.lhs().cols())
{
EIGEN_INTERNAL_CHECK_COST_VALUE(NumTraits<Scalar>::MulCost);
EIGEN_INTERNAL_CHECK_COST_VALUE(NumTraits<Scalar>::AddCost);
EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost);
#if 0
std::cerr << "LhsOuterStrideBytes= " << LhsOuterStrideBytes << "\n";
std::cerr << "RhsOuterStrideBytes= " << RhsOuterStrideBytes << "\n";
std::cerr << "LhsAlignment= " << LhsAlignment << "\n";
std::cerr << "RhsAlignment= " << RhsAlignment << "\n";
std::cerr << "CanVectorizeLhs= " << CanVectorizeLhs << "\n";
std::cerr << "CanVectorizeRhs= " << CanVectorizeRhs << "\n";
std::cerr << "CanVectorizeInner= " << CanVectorizeInner << "\n";
std::cerr << "EvalToRowMajor= " << EvalToRowMajor << "\n";
std::cerr << "Alignment= " << Alignment << "\n";
std::cerr << "Flags= " << Flags << "\n";
#endif
}
// Everything below here is taken from CoeffBasedProduct.h
typedef typename internal::nested_eval<Lhs,Rhs::ColsAtCompileTime>::type LhsNested;
typedef typename internal::nested_eval<Rhs,Lhs::RowsAtCompileTime>::type RhsNested;
typedef typename internal::remove_all<LhsNested>::type LhsNestedCleaned;
typedef typename internal::remove_all<RhsNested>::type RhsNestedCleaned;
typedef evaluator<LhsNestedCleaned> LhsEtorType;
typedef evaluator<RhsNestedCleaned> RhsEtorType;
enum {
RowsAtCompileTime = LhsNestedCleaned::RowsAtCompileTime,
ColsAtCompileTime = RhsNestedCleaned::ColsAtCompileTime,
InnerSize = EIGEN_SIZE_MIN_PREFER_FIXED(LhsNestedCleaned::ColsAtCompileTime, RhsNestedCleaned::RowsAtCompileTime),
MaxRowsAtCompileTime = LhsNestedCleaned::MaxRowsAtCompileTime,
MaxColsAtCompileTime = RhsNestedCleaned::MaxColsAtCompileTime
};
typedef typename find_best_packet<Scalar,RowsAtCompileTime>::type LhsVecPacketType;
typedef typename find_best_packet<Scalar,ColsAtCompileTime>::type RhsVecPacketType;
enum {
LhsCoeffReadCost = LhsEtorType::CoeffReadCost,
RhsCoeffReadCost = RhsEtorType::CoeffReadCost,
CoeffReadCost = InnerSize==0 ? NumTraits<Scalar>::ReadCost
: InnerSize == Dynamic ? HugeCost
: InnerSize * (NumTraits<Scalar>::MulCost + LhsCoeffReadCost + RhsCoeffReadCost)
+ (InnerSize - 1) * NumTraits<Scalar>::AddCost,
Unroll = CoeffReadCost <= EIGEN_UNROLLING_LIMIT,
LhsFlags = LhsEtorType::Flags,
RhsFlags = RhsEtorType::Flags,
LhsRowMajor = LhsFlags & RowMajorBit,
RhsRowMajor = RhsFlags & RowMajorBit,
LhsVecPacketSize = unpacket_traits<LhsVecPacketType>::size,
RhsVecPacketSize = unpacket_traits<RhsVecPacketType>::size,
// Here, we don't care about alignment larger than the usable packet size.
LhsAlignment = EIGEN_PLAIN_ENUM_MIN(LhsEtorType::Alignment,LhsVecPacketSize*int(sizeof(typename LhsNestedCleaned::Scalar))),
RhsAlignment = EIGEN_PLAIN_ENUM_MIN(RhsEtorType::Alignment,RhsVecPacketSize*int(sizeof(typename RhsNestedCleaned::Scalar))),
SameType = is_same<typename LhsNestedCleaned::Scalar,typename RhsNestedCleaned::Scalar>::value,
CanVectorizeRhs = bool(RhsRowMajor) && (RhsFlags & PacketAccessBit) && (ColsAtCompileTime!=1),
CanVectorizeLhs = (!LhsRowMajor) && (LhsFlags & PacketAccessBit) && (RowsAtCompileTime!=1),
EvalToRowMajor = (MaxRowsAtCompileTime==1&&MaxColsAtCompileTime!=1) ? 1
: (MaxColsAtCompileTime==1&&MaxRowsAtCompileTime!=1) ? 0
: (bool(RhsRowMajor) && !CanVectorizeLhs),
Flags = ((unsigned int)(LhsFlags | RhsFlags) & HereditaryBits & ~RowMajorBit)
| (EvalToRowMajor ? RowMajorBit : 0)
// TODO enable vectorization for mixed types
| (SameType && (CanVectorizeLhs || CanVectorizeRhs) ? PacketAccessBit : 0)
| (XprType::IsVectorAtCompileTime ? LinearAccessBit : 0),
LhsOuterStrideBytes = int(LhsNestedCleaned::OuterStrideAtCompileTime) * int(sizeof(typename LhsNestedCleaned::Scalar)),
RhsOuterStrideBytes = int(RhsNestedCleaned::OuterStrideAtCompileTime) * int(sizeof(typename RhsNestedCleaned::Scalar)),
Alignment = bool(CanVectorizeLhs) ? (LhsOuterStrideBytes<=0 || (int(LhsOuterStrideBytes) % EIGEN_PLAIN_ENUM_MAX(1,LhsAlignment))!=0 ? 0 : LhsAlignment)
: bool(CanVectorizeRhs) ? (RhsOuterStrideBytes<=0 || (int(RhsOuterStrideBytes) % EIGEN_PLAIN_ENUM_MAX(1,RhsAlignment))!=0 ? 0 : RhsAlignment)
: 0,
/* CanVectorizeInner deserves special explanation. It does not affect the product flags. It is not used outside
* of Product. If the Product itself is not a packet-access expression, there is still a chance that the inner
* loop of the product might be vectorized. This is the meaning of CanVectorizeInner. Since it doesn't affect
* the Flags, it is safe to make this value depend on ActualPacketAccessBit, that doesn't affect the ABI.
*/
CanVectorizeInner = SameType
&& LhsRowMajor
&& (!RhsRowMajor)
&& (LhsFlags & RhsFlags & ActualPacketAccessBit)
&& (InnerSize % packet_traits<Scalar>::size == 0)
};
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const CoeffReturnType coeff(Index row, Index col) const
{
return (m_lhs.row(row).transpose().cwiseProduct( m_rhs.col(col) )).sum();
}
/* Allow index-based non-packet access. It is impossible though to allow index-based packed access,
* which is why we don't set the LinearAccessBit.
* TODO: this seems possible when the result is a vector
*/
EIGEN_DEVICE_FUNC const CoeffReturnType coeff(Index index) const
{
const Index row = (RowsAtCompileTime == 1 || MaxRowsAtCompileTime==1) ? 0 : index;
const Index col = (RowsAtCompileTime == 1 || MaxRowsAtCompileTime==1) ? index : 0;
return (m_lhs.row(row).transpose().cwiseProduct( m_rhs.col(col) )).sum();
}
template<int LoadMode, typename PacketType>
const PacketType packet(Index row, Index col) const
{
PacketType res;
typedef etor_product_packet_impl<bool(int(Flags)&RowMajorBit) ? RowMajor : ColMajor,
Unroll ? int(InnerSize) : Dynamic,
LhsEtorType, RhsEtorType, PacketType, LoadMode> PacketImpl;
PacketImpl::run(row, col, m_lhsImpl, m_rhsImpl, m_innerDim, res);
return res;
}
template<int LoadMode, typename PacketType>
const PacketType packet(Index index) const
{
const Index row = (RowsAtCompileTime == 1 || MaxRowsAtCompileTime==1) ? 0 : index;
const Index col = (RowsAtCompileTime == 1 || MaxRowsAtCompileTime==1) ? index : 0;
return packet<LoadMode,PacketType>(row,col);
}
protected:
typename internal::add_const_on_value_type<LhsNested>::type m_lhs;
typename internal::add_const_on_value_type<RhsNested>::type m_rhs;
LhsEtorType m_lhsImpl;
RhsEtorType m_rhsImpl;
// TODO: Get rid of m_innerDim if known at compile time
Index m_innerDim;
};
template<typename Lhs, typename Rhs>
struct product_evaluator<Product<Lhs, Rhs, DefaultProduct>, LazyCoeffBasedProductMode, DenseShape, DenseShape>
: product_evaluator<Product<Lhs, Rhs, LazyProduct>, CoeffBasedProductMode, DenseShape, DenseShape>
{
typedef Product<Lhs, Rhs, DefaultProduct> XprType;
typedef Product<Lhs, Rhs, LazyProduct> BaseProduct;
typedef product_evaluator<BaseProduct, CoeffBasedProductMode, DenseShape, DenseShape> Base;
enum {
Flags = Base::Flags | EvalBeforeNestingBit
};
EIGEN_DEVICE_FUNC explicit product_evaluator(const XprType& xpr)
: Base(BaseProduct(xpr.lhs(),xpr.rhs()))
{}
};
/****************************************
*** Coeff based product, Packet path ***
****************************************/
template<int UnrollingIndex, typename Lhs, typename Rhs, typename Packet, int LoadMode>
struct etor_product_packet_impl<RowMajor, UnrollingIndex, Lhs, Rhs, Packet, LoadMode>
{
static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, Packet &res)
{
etor_product_packet_impl<RowMajor, UnrollingIndex-1, Lhs, Rhs, Packet, LoadMode>::run(row, col, lhs, rhs, innerDim, res);
res = pmadd(pset1<Packet>(lhs.coeff(row, Index(UnrollingIndex-1))), rhs.template packet<LoadMode,Packet>(Index(UnrollingIndex-1), col), res);
}
};
template<int UnrollingIndex, typename Lhs, typename Rhs, typename Packet, int LoadMode>
struct etor_product_packet_impl<ColMajor, UnrollingIndex, Lhs, Rhs, Packet, LoadMode>
{
static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, Packet &res)
{
etor_product_packet_impl<ColMajor, UnrollingIndex-1, Lhs, Rhs, Packet, LoadMode>::run(row, col, lhs, rhs, innerDim, res);
res = pmadd(lhs.template packet<LoadMode,Packet>(row, Index(UnrollingIndex-1)), pset1<Packet>(rhs.coeff(Index(UnrollingIndex-1), col)), res);
}
};
template<typename Lhs, typename Rhs, typename Packet, int LoadMode>
struct etor_product_packet_impl<RowMajor, 1, Lhs, Rhs, Packet, LoadMode>
{
static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index /*innerDim*/, Packet &res)
{
res = pmul(pset1<Packet>(lhs.coeff(row, Index(0))),rhs.template packet<LoadMode,Packet>(Index(0), col));
}
};
template<typename Lhs, typename Rhs, typename Packet, int LoadMode>
struct etor_product_packet_impl<ColMajor, 1, Lhs, Rhs, Packet, LoadMode>
{
static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index /*innerDim*/, Packet &res)
{
res = pmul(lhs.template packet<LoadMode,Packet>(row, Index(0)), pset1<Packet>(rhs.coeff(Index(0), col)));
}
};
template<typename Lhs, typename Rhs, typename Packet, int LoadMode>
struct etor_product_packet_impl<RowMajor, 0, Lhs, Rhs, Packet, LoadMode>
{
static EIGEN_STRONG_INLINE void run(Index /*row*/, Index /*col*/, const Lhs& /*lhs*/, const Rhs& /*rhs*/, Index /*innerDim*/, Packet &res)
{
res = pset1<Packet>(typename unpacket_traits<Packet>::type(0));
}
};
template<typename Lhs, typename Rhs, typename Packet, int LoadMode>
struct etor_product_packet_impl<ColMajor, 0, Lhs, Rhs, Packet, LoadMode>
{
static EIGEN_STRONG_INLINE void run(Index /*row*/, Index /*col*/, const Lhs& /*lhs*/, const Rhs& /*rhs*/, Index /*innerDim*/, Packet &res)
{
res = pset1<Packet>(typename unpacket_traits<Packet>::type(0));
}
};
template<typename Lhs, typename Rhs, typename Packet, int LoadMode>
struct etor_product_packet_impl<RowMajor, Dynamic, Lhs, Rhs, Packet, LoadMode>
{
static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, Packet& res)
{
res = pset1<Packet>(typename unpacket_traits<Packet>::type(0));
for(Index i = 0; i < innerDim; ++i)
res = pmadd(pset1<Packet>(lhs.coeff(row, i)), rhs.template packet<LoadMode,Packet>(i, col), res);
}
};
template<typename Lhs, typename Rhs, typename Packet, int LoadMode>
struct etor_product_packet_impl<ColMajor, Dynamic, Lhs, Rhs, Packet, LoadMode>
{
static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, Packet& res)
{
res = pset1<Packet>(typename unpacket_traits<Packet>::type(0));
for(Index i = 0; i < innerDim; ++i)
res = pmadd(lhs.template packet<LoadMode,Packet>(row, i), pset1<Packet>(rhs.coeff(i, col)), res);
}
};
/***************************************************************************
* Triangular products
***************************************************************************/
template<int Mode, bool LhsIsTriangular,
typename Lhs, bool LhsIsVector,
typename Rhs, bool RhsIsVector>
struct triangular_product_impl;
template<typename Lhs, typename Rhs, int ProductTag>
struct generic_product_impl<Lhs,Rhs,TriangularShape,DenseShape,ProductTag>
: generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,TriangularShape,DenseShape,ProductTag> >
{
typedef typename Product<Lhs,Rhs>::Scalar Scalar;
template<typename Dest>
static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
{
triangular_product_impl<Lhs::Mode,true,typename Lhs::MatrixType,false,Rhs, Rhs::ColsAtCompileTime==1>
::run(dst, lhs.nestedExpression(), rhs, alpha);
}
};
template<typename Lhs, typename Rhs, int ProductTag>
struct generic_product_impl<Lhs,Rhs,DenseShape,TriangularShape,ProductTag>
: generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,DenseShape,TriangularShape,ProductTag> >
{
typedef typename Product<Lhs,Rhs>::Scalar Scalar;
template<typename Dest>
static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
{
triangular_product_impl<Rhs::Mode,false,Lhs,Lhs::RowsAtCompileTime==1, typename Rhs::MatrixType, false>::run(dst, lhs, rhs.nestedExpression(), alpha);
}
};
/***************************************************************************
* SelfAdjoint products
***************************************************************************/
template <typename Lhs, int LhsMode, bool LhsIsVector,
typename Rhs, int RhsMode, bool RhsIsVector>
struct selfadjoint_product_impl;
template<typename Lhs, typename Rhs, int ProductTag>
struct generic_product_impl<Lhs,Rhs,SelfAdjointShape,DenseShape,ProductTag>
: generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,SelfAdjointShape,DenseShape,ProductTag> >
{
typedef typename Product<Lhs,Rhs>::Scalar Scalar;
template<typename Dest>
static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
{
selfadjoint_product_impl<typename Lhs::MatrixType,Lhs::Mode,false,Rhs,0,Rhs::IsVectorAtCompileTime>::run(dst, lhs.nestedExpression(), rhs, alpha);
}
};
template<typename Lhs, typename Rhs, int ProductTag>
struct generic_product_impl<Lhs,Rhs,DenseShape,SelfAdjointShape,ProductTag>
: generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,DenseShape,SelfAdjointShape,ProductTag> >
{
typedef typename Product<Lhs,Rhs>::Scalar Scalar;
template<typename Dest>
static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
{
selfadjoint_product_impl<Lhs,0,Lhs::IsVectorAtCompileTime,typename Rhs::MatrixType,Rhs::Mode,false>::run(dst, lhs, rhs.nestedExpression(), alpha);
}
};
/***************************************************************************
* Diagonal products
***************************************************************************/
template<typename MatrixType, typename DiagonalType, typename Derived, int ProductOrder>
struct diagonal_product_evaluator_base
: evaluator_base<Derived>
{
typedef typename ScalarBinaryOpTraits<typename MatrixType::Scalar, typename DiagonalType::Scalar>::ReturnType Scalar;
public:
enum {
CoeffReadCost = NumTraits<Scalar>::MulCost + evaluator<MatrixType>::CoeffReadCost + evaluator<DiagonalType>::CoeffReadCost,
MatrixFlags = evaluator<MatrixType>::Flags,
DiagFlags = evaluator<DiagonalType>::Flags,
_StorageOrder = MatrixFlags & RowMajorBit ? RowMajor : ColMajor,
_ScalarAccessOnDiag = !((int(_StorageOrder) == ColMajor && int(ProductOrder) == OnTheLeft)
||(int(_StorageOrder) == RowMajor && int(ProductOrder) == OnTheRight)),
_SameTypes = is_same<typename MatrixType::Scalar, typename DiagonalType::Scalar>::value,
// FIXME currently we need same types, but in the future the next rule should be the one
//_Vectorizable = bool(int(MatrixFlags)&PacketAccessBit) && ((!_PacketOnDiag) || (_SameTypes && bool(int(DiagFlags)&PacketAccessBit))),
_Vectorizable = bool(int(MatrixFlags)&PacketAccessBit) && _SameTypes && (_ScalarAccessOnDiag || (bool(int(DiagFlags)&PacketAccessBit))),
_LinearAccessMask = (MatrixType::RowsAtCompileTime==1 || MatrixType::ColsAtCompileTime==1) ? LinearAccessBit : 0,
Flags = ((HereditaryBits|_LinearAccessMask) & (unsigned int)(MatrixFlags)) | (_Vectorizable ? PacketAccessBit : 0),
Alignment = evaluator<MatrixType>::Alignment
};
diagonal_product_evaluator_base(const MatrixType &mat, const DiagonalType &diag)
: m_diagImpl(diag), m_matImpl(mat)
{
EIGEN_INTERNAL_CHECK_COST_VALUE(NumTraits<Scalar>::MulCost);
EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost);
}
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar coeff(Index idx) const
{
return m_diagImpl.coeff(idx) * m_matImpl.coeff(idx);
}
protected:
template<int LoadMode,typename PacketType>
EIGEN_STRONG_INLINE PacketType packet_impl(Index row, Index col, Index id, internal::true_type) const
{
return internal::pmul(m_matImpl.template packet<LoadMode,PacketType>(row, col),
internal::pset1<PacketType>(m_diagImpl.coeff(id)));
}
template<int LoadMode,typename PacketType>
EIGEN_STRONG_INLINE PacketType packet_impl(Index row, Index col, Index id, internal::false_type) const
{
enum {
InnerSize = (MatrixType::Flags & RowMajorBit) ? MatrixType::ColsAtCompileTime : MatrixType::RowsAtCompileTime,
DiagonalPacketLoadMode = EIGEN_PLAIN_ENUM_MIN(LoadMode,((InnerSize%16) == 0) ? int(Aligned16) : int(evaluator<DiagonalType>::Alignment)) // FIXME hardcoded 16!!
};
return internal::pmul(m_matImpl.template packet<LoadMode,PacketType>(row, col),
m_diagImpl.template packet<DiagonalPacketLoadMode,PacketType>(id));
}
evaluator<DiagonalType> m_diagImpl;
evaluator<MatrixType> m_matImpl;
};
// diagonal * dense
template<typename Lhs, typename Rhs, int ProductKind, int ProductTag>
struct product_evaluator<Product<Lhs, Rhs, ProductKind>, ProductTag, DiagonalShape, DenseShape>
: diagonal_product_evaluator_base<Rhs, typename Lhs::DiagonalVectorType, Product<Lhs, Rhs, LazyProduct>, OnTheLeft>
{
typedef diagonal_product_evaluator_base<Rhs, typename Lhs::DiagonalVectorType, Product<Lhs, Rhs, LazyProduct>, OnTheLeft> Base;
using Base::m_diagImpl;
using Base::m_matImpl;
using Base::coeff;
typedef typename Base::Scalar Scalar;
typedef Product<Lhs, Rhs, ProductKind> XprType;
typedef typename XprType::PlainObject PlainObject;
enum {
StorageOrder = int(Rhs::Flags) & RowMajorBit ? RowMajor : ColMajor
};
EIGEN_DEVICE_FUNC explicit product_evaluator(const XprType& xpr)
: Base(xpr.rhs(), xpr.lhs().diagonal())
{
}
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar coeff(Index row, Index col) const
{
return m_diagImpl.coeff(row) * m_matImpl.coeff(row, col);
}
#ifndef __CUDACC__
template<int LoadMode,typename PacketType>
EIGEN_STRONG_INLINE PacketType packet(Index row, Index col) const
{
// FIXME: NVCC used to complain about the template keyword, but we have to check whether this is still the case.
// See also similar calls below.
return this->template packet_impl<LoadMode,PacketType>(row,col, row,
typename internal::conditional<int(StorageOrder)==RowMajor, internal::true_type, internal::false_type>::type());
}
template<int LoadMode,typename PacketType>
EIGEN_STRONG_INLINE PacketType packet(Index idx) const
{
return packet<LoadMode,PacketType>(int(StorageOrder)==ColMajor?idx:0,int(StorageOrder)==ColMajor?0:idx);
}
#endif
};
// dense * diagonal
template<typename Lhs, typename Rhs, int ProductKind, int ProductTag>
struct product_evaluator<Product<Lhs, Rhs, ProductKind>, ProductTag, DenseShape, DiagonalShape>
: diagonal_product_evaluator_base<Lhs, typename Rhs::DiagonalVectorType, Product<Lhs, Rhs, LazyProduct>, OnTheRight>
{
typedef diagonal_product_evaluator_base<Lhs, typename Rhs::DiagonalVectorType, Product<Lhs, Rhs, LazyProduct>, OnTheRight> Base;
using Base::m_diagImpl;
using Base::m_matImpl;
using Base::coeff;
typedef typename Base::Scalar Scalar;
typedef Product<Lhs, Rhs, ProductKind> XprType;
typedef typename XprType::PlainObject PlainObject;
enum { StorageOrder = int(Lhs::Flags) & RowMajorBit ? RowMajor : ColMajor };
EIGEN_DEVICE_FUNC explicit product_evaluator(const XprType& xpr)
: Base(xpr.lhs(), xpr.rhs().diagonal())
{
}
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar coeff(Index row, Index col) const
{
return m_matImpl.coeff(row, col) * m_diagImpl.coeff(col);
}
#ifndef __CUDACC__
template<int LoadMode,typename PacketType>
EIGEN_STRONG_INLINE PacketType packet(Index row, Index col) const
{
return this->template packet_impl<LoadMode,PacketType>(row,col, col,
typename internal::conditional<int(StorageOrder)==ColMajor, internal::true_type, internal::false_type>::type());
}
template<int LoadMode,typename PacketType>
EIGEN_STRONG_INLINE PacketType packet(Index idx) const
{
return packet<LoadMode,PacketType>(int(StorageOrder)==ColMajor?idx:0,int(StorageOrder)==ColMajor?0:idx);
}
#endif
};
/***************************************************************************
* Products with permutation matrices
***************************************************************************/
/** \internal
* \class permutation_matrix_product
* Internal helper class implementing the product between a permutation matrix and a matrix.
* This class is specialized for DenseShape below and for SparseShape in SparseCore/SparsePermutation.h
*/
template<typename ExpressionType, int Side, bool Transposed, typename ExpressionShape>
struct permutation_matrix_product;
template<typename ExpressionType, int Side, bool Transposed>
struct permutation_matrix_product<ExpressionType, Side, Transposed, DenseShape>
{
typedef typename nested_eval<ExpressionType, 1>::type MatrixType;
typedef typename remove_all<MatrixType>::type MatrixTypeCleaned;
template<typename Dest, typename PermutationType>
static inline void run(Dest& dst, const PermutationType& perm, const ExpressionType& xpr)
{
MatrixType mat(xpr);
const Index n = Side==OnTheLeft ? mat.rows() : mat.cols();
// FIXME we need an is_same for expression that is not sensitive to constness. For instance
// is_same_xpr<Block<const Matrix>, Block<Matrix> >::value should be true.
//if(is_same<MatrixTypeCleaned,Dest>::value && extract_data(dst) == extract_data(mat))
if(is_same_dense(dst, mat))
{
// apply the permutation inplace
Matrix<bool,PermutationType::RowsAtCompileTime,1,0,PermutationType::MaxRowsAtCompileTime> mask(perm.size());
mask.fill(false);
Index r = 0;
while(r < perm.size())
{
// search for the next seed
while(r<perm.size() && mask[r]) r++;
if(r>=perm.size())
break;
// we got one, let's follow it until we are back to the seed
Index k0 = r++;
Index kPrev = k0;
mask.coeffRef(k0) = true;
for(Index k=perm.indices().coeff(k0); k!=k0; k=perm.indices().coeff(k))
{
Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>(dst, k)
.swap(Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
(dst,((Side==OnTheLeft) ^ Transposed) ? k0 : kPrev));
mask.coeffRef(k) = true;
kPrev = k;
}
}
}
else
{
for(Index i = 0; i < n; ++i)
{
Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
(dst, ((Side==OnTheLeft) ^ Transposed) ? perm.indices().coeff(i) : i)
=
Block<const MatrixTypeCleaned,Side==OnTheLeft ? 1 : MatrixTypeCleaned::RowsAtCompileTime,Side==OnTheRight ? 1 : MatrixTypeCleaned::ColsAtCompileTime>
(mat, ((Side==OnTheRight) ^ Transposed) ? perm.indices().coeff(i) : i);
}
}
}
};
template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
struct generic_product_impl<Lhs, Rhs, PermutationShape, MatrixShape, ProductTag>
{
template<typename Dest>
static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
{
permutation_matrix_product<Rhs, OnTheLeft, false, MatrixShape>::run(dst, lhs, rhs);
}
};
template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
struct generic_product_impl<Lhs, Rhs, MatrixShape, PermutationShape, ProductTag>
{
template<typename Dest>
static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
{
permutation_matrix_product<Lhs, OnTheRight, false, MatrixShape>::run(dst, rhs, lhs);
}
};
template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
struct generic_product_impl<Inverse<Lhs>, Rhs, PermutationShape, MatrixShape, ProductTag>
{
template<typename Dest>
static void evalTo(Dest& dst, const Inverse<Lhs>& lhs, const Rhs& rhs)
{
permutation_matrix_product<Rhs, OnTheLeft, true, MatrixShape>::run(dst, lhs.nestedExpression(), rhs);
}
};
template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
struct generic_product_impl<Lhs, Inverse<Rhs>, MatrixShape, PermutationShape, ProductTag>
{
template<typename Dest>
static void evalTo(Dest& dst, const Lhs& lhs, const Inverse<Rhs>& rhs)
{
permutation_matrix_product<Lhs, OnTheRight, true, MatrixShape>::run(dst, rhs.nestedExpression(), lhs);
}
};
/***************************************************************************
* Products with transpositions matrices
***************************************************************************/
// FIXME could we unify Transpositions and Permutation into a single "shape"??
/** \internal
* \class transposition_matrix_product
* Internal helper class implementing the product between a permutation matrix and a matrix.
*/
template<typename ExpressionType, int Side, bool Transposed, typename ExpressionShape>
struct transposition_matrix_product
{
typedef typename nested_eval<ExpressionType, 1>::type MatrixType;
typedef typename remove_all<MatrixType>::type MatrixTypeCleaned;
template<typename Dest, typename TranspositionType>
static inline void run(Dest& dst, const TranspositionType& tr, const ExpressionType& xpr)
{
MatrixType mat(xpr);
typedef typename TranspositionType::StorageIndex StorageIndex;
const Index size = tr.size();
StorageIndex j = 0;
if(!is_same_dense(dst,mat))
dst = mat;
for(Index k=(Transposed?size-1:0) ; Transposed?k>=0:k<size ; Transposed?--k:++k)
if(Index(j=tr.coeff(k))!=k)
{
if(Side==OnTheLeft) dst.row(k).swap(dst.row(j));
else if(Side==OnTheRight) dst.col(k).swap(dst.col(j));
}
}
};
template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
struct generic_product_impl<Lhs, Rhs, TranspositionsShape, MatrixShape, ProductTag>
{
template<typename Dest>
static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
{
transposition_matrix_product<Rhs, OnTheLeft, false, MatrixShape>::run(dst, lhs, rhs);
}
};
template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
struct generic_product_impl<Lhs, Rhs, MatrixShape, TranspositionsShape, ProductTag>
{
template<typename Dest>
static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
{
transposition_matrix_product<Lhs, OnTheRight, false, MatrixShape>::run(dst, rhs, lhs);
}
};
template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
struct generic_product_impl<Transpose<Lhs>, Rhs, TranspositionsShape, MatrixShape, ProductTag>
{
template<typename Dest>
static void evalTo(Dest& dst, const Transpose<Lhs>& lhs, const Rhs& rhs)
{
transposition_matrix_product<Rhs, OnTheLeft, true, MatrixShape>::run(dst, lhs.nestedExpression(), rhs);
}
};
template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
struct generic_product_impl<Lhs, Transpose<Rhs>, MatrixShape, TranspositionsShape, ProductTag>
{
template<typename Dest>
static void evalTo(Dest& dst, const Lhs& lhs, const Transpose<Rhs>& rhs)
{
transposition_matrix_product<Lhs, OnTheRight, true, MatrixShape>::run(dst, rhs.nestedExpression(), lhs);
}
};
} // end namespace internal
} // end namespace Eigen
#endif // EIGEN_PRODUCT_EVALUATORS_H
external/eigen3/Eigen/src/Core/Random.h 0 → 100644
View file @ a394b22a
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_RANDOM_H
#define EIGEN_RANDOM_H
namespace Eigen {
namespace internal {
template<typename Scalar> struct scalar_random_op {
EIGEN_EMPTY_STRUCT_CTOR(scalar_random_op)
inline const Scalar operator() () const { return random<Scalar>(); }
};
template<typename Scalar>
struct functor_traits<scalar_random_op<Scalar> >
{ enum { Cost = 5 * NumTraits<Scalar>::MulCost, PacketAccess = false, IsRepeatable = false }; };
} // end namespace internal
/** \returns a random matrix expression
*
* Numbers are uniformly spread through their whole definition range for integer types,
* and in the [-1:1] range for floating point scalar types.
*
* The parameters \a rows and \a cols are the number of rows and of columns of
* the returned matrix. Must be compatible with this MatrixBase type.
*
* \not_reentrant
*
* This variant is meant to be used for dynamic-size matrix types. For fixed-size types,
* it is redundant to pass \a rows and \a cols as arguments, so Random() should be used
* instead.
*
*
* Example: \include MatrixBase_random_int_int.cpp
* Output: \verbinclude MatrixBase_random_int_int.out
*
* This expression has the "evaluate before nesting" flag so that it will be evaluated into
* a temporary matrix whenever it is nested in a larger expression. This prevents unexpected
* behavior with expressions involving random matrices.
*
* See DenseBase::NullaryExpr(Index, const CustomNullaryOp&) for an example using C++11 random generators.
*
* \sa DenseBase::setRandom(), DenseBase::Random(Index), DenseBase::Random()
*/
template<typename Derived>
inline const typename DenseBase<Derived>::RandomReturnType
DenseBase<Derived>::Random(Index rows, Index cols)
{
return NullaryExpr(rows, cols, internal::scalar_random_op<Scalar>());
}
/** \returns a random vector expression
*
* Numbers are uniformly spread through their whole definition range for integer types,
* and in the [-1:1] range for floating point scalar types.
*
* The parameter \a size is the size of the returned vector.
* Must be compatible with this MatrixBase type.
*
* \only_for_vectors
* \not_reentrant
*
* This variant is meant to be used for dynamic-size vector types. For fixed-size types,
* it is redundant to pass \a size as argument, so Random() should be used
* instead.
*
* Example: \include MatrixBase_random_int.cpp
* Output: \verbinclude MatrixBase_random_int.out
*
* This expression has the "evaluate before nesting" flag so that it will be evaluated into
* a temporary vector whenever it is nested in a larger expression. This prevents unexpected
* behavior with expressions involving random matrices.
*
* \sa DenseBase::setRandom(), DenseBase::Random(Index,Index), DenseBase::Random()
*/
template<typename Derived>
inline const typename DenseBase<Derived>::RandomReturnType
DenseBase<Derived>::Random(Index size)
{
return NullaryExpr(size, internal::scalar_random_op<Scalar>());
}
/** \returns a fixed-size random matrix or vector expression
*
* Numbers are uniformly spread through their whole definition range for integer types,
* and in the [-1:1] range for floating point scalar types.
*
* This variant is only for fixed-size MatrixBase types. For dynamic-size types, you
* need to use the variants taking size arguments.
*
* Example: \include MatrixBase_random.cpp
* Output: \verbinclude MatrixBase_random.out
*
* This expression has the "evaluate before nesting" flag so that it will be evaluated into
* a temporary matrix whenever it is nested in a larger expression. This prevents unexpected
* behavior with expressions involving random matrices.
*
* \not_reentrant
*
* \sa DenseBase::setRandom(), DenseBase::Random(Index,Index), DenseBase::Random(Index)
*/
template<typename Derived>
inline const typename DenseBase<Derived>::RandomReturnType
DenseBase<Derived>::Random()
{
return NullaryExpr(RowsAtCompileTime, ColsAtCompileTime, internal::scalar_random_op<Scalar>());
}
/** Sets all coefficients in this expression to random values.
*
* Numbers are uniformly spread through their whole definition range for integer types,
* and in the [-1:1] range for floating point scalar types.
*
* \not_reentrant
*
* Example: \include MatrixBase_setRandom.cpp
* Output: \verbinclude MatrixBase_setRandom.out
*
* \sa class CwiseNullaryOp, setRandom(Index), setRandom(Index,Index)
*/
template<typename Derived>
inline Derived& DenseBase<Derived>::setRandom()
{
return *this = Random(rows(), cols());
}
/** Resizes to the given \a newSize, and sets all coefficients in this expression to random values.
*
* Numbers are uniformly spread through their whole definition range for integer types,
* and in the [-1:1] range for floating point scalar types.
*
* \only_for_vectors
* \not_reentrant
*
* Example: \include Matrix_setRandom_int.cpp
* Output: \verbinclude Matrix_setRandom_int.out
*
* \sa DenseBase::setRandom(), setRandom(Index,Index), class CwiseNullaryOp, DenseBase::Random()
*/
template<typename Derived>
EIGEN_STRONG_INLINE Derived&
PlainObjectBase<Derived>::setRandom(Index newSize)
{
resize(newSize);
return setRandom();
}
/** Resizes to the given size, and sets all coefficients in this expression to random values.
*
* Numbers are uniformly spread through their whole definition range for integer types,
* and in the [-1:1] range for floating point scalar types.
*
* \not_reentrant
*
* \param rows the new number of rows
* \param cols the new number of columns
*
* Example: \include Matrix_setRandom_int_int.cpp
* Output: \verbinclude Matrix_setRandom_int_int.out
*
* \sa DenseBase::setRandom(), setRandom(Index), class CwiseNullaryOp, DenseBase::Random()
*/
template<typename Derived>
EIGEN_STRONG_INLINE Derived&
PlainObjectBase<Derived>::setRandom(Index rows, Index cols)
{
resize(rows, cols);
return setRandom();
}
} // end namespace Eigen
#endif // EIGEN_RANDOM_H
external/eigen3/Eigen/src/Core/Redux.h 0 → 100644
View file @ a394b22a
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_REDUX_H
#define EIGEN_REDUX_H
namespace Eigen {
namespace internal {
// TODO
// * implement other kind of vectorization
// * factorize code
/***************************************************************************
* Part 1 : the logic deciding a strategy for vectorization and unrolling
***************************************************************************/
template<typename Func, typename Derived>
struct redux_traits
{
public:
typedef typename find_best_packet<typename Derived::Scalar,Derived::SizeAtCompileTime>::type PacketType;
enum {
PacketSize = unpacket_traits<PacketType>::size,
InnerMaxSize = int(Derived::IsRowMajor)
? Derived::MaxColsAtCompileTime
: Derived::MaxRowsAtCompileTime
};
enum {
MightVectorize = (int(Derived::Flags)&ActualPacketAccessBit)
&& (functor_traits<Func>::PacketAccess),
MayLinearVectorize = bool(MightVectorize) && (int(Derived::Flags)&LinearAccessBit),
MaySliceVectorize = bool(MightVectorize) && int(InnerMaxSize)>=3*PacketSize
};
public:
enum {
Traversal = int(MayLinearVectorize) ? int(LinearVectorizedTraversal)
: int(MaySliceVectorize) ? int(SliceVectorizedTraversal)
: int(DefaultTraversal)
};
public:
enum {
Cost = Derived::SizeAtCompileTime == Dynamic ? HugeCost
: Derived::SizeAtCompileTime * Derived::CoeffReadCost + (Derived::SizeAtCompileTime-1) * functor_traits<Func>::Cost,
UnrollingLimit = EIGEN_UNROLLING_LIMIT * (int(Traversal) == int(DefaultTraversal) ? 1 : int(PacketSize))
};
public:
enum {
Unrolling = Cost <= UnrollingLimit ? CompleteUnrolling : NoUnrolling
};
#ifdef EIGEN_DEBUG_ASSIGN
static void debug()
{
std::cerr << "Xpr: " << typeid(typename Derived::XprType).name() << std::endl;
std::cerr.setf(std::ios::hex, std::ios::basefield);
EIGEN_DEBUG_VAR(Derived::Flags)
std::cerr.unsetf(std::ios::hex);
EIGEN_DEBUG_VAR(InnerMaxSize)
EIGEN_DEBUG_VAR(PacketSize)
EIGEN_DEBUG_VAR(MightVectorize)
EIGEN_DEBUG_VAR(MayLinearVectorize)
EIGEN_DEBUG_VAR(MaySliceVectorize)
EIGEN_DEBUG_VAR(Traversal)
EIGEN_DEBUG_VAR(UnrollingLimit)
EIGEN_DEBUG_VAR(Unrolling)
std::cerr << std::endl;
}
#endif
};
/***************************************************************************
* Part 2 : unrollers
***************************************************************************/
/*** no vectorization ***/
template<typename Func, typename Derived, int Start, int Length>
struct redux_novec_unroller
{
enum {
HalfLength = Length/2
};
typedef typename Derived::Scalar Scalar;
EIGEN_DEVICE_FUNC
static EIGEN_STRONG_INLINE Scalar run(const Derived &mat, const Func& func)
{
return func(redux_novec_unroller<Func, Derived, Start, HalfLength>::run(mat,func),
redux_novec_unroller<Func, Derived, Start+HalfLength, Length-HalfLength>::run(mat,func));
}
};
template<typename Func, typename Derived, int Start>
struct redux_novec_unroller<Func, Derived, Start, 1>
{
enum {
outer = Start / Derived::InnerSizeAtCompileTime,
inner = Start % Derived::InnerSizeAtCompileTime
};
typedef typename Derived::Scalar Scalar;
EIGEN_DEVICE_FUNC
static EIGEN_STRONG_INLINE Scalar run(const Derived &mat, const Func&)
{
return mat.coeffByOuterInner(outer, inner);
}
};
// 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 redux_novec_unroller<Func, Derived, Start, 0>
{
typedef typename Derived::Scalar Scalar;
EIGEN_DEVICE_FUNC
static EIGEN_STRONG_INLINE Scalar run(const Derived&, const Func&) { return Scalar(); }
};
/*** vectorization ***/
template<typename Func, typename Derived, int Start, int Length>
struct redux_vec_unroller
{
enum {
PacketSize = redux_traits<Func, Derived>::PacketSize,
HalfLength = Length/2
};
typedef typename Derived::Scalar Scalar;
typedef typename redux_traits<Func, Derived>::PacketType PacketScalar;
static EIGEN_STRONG_INLINE PacketScalar run(const Derived &mat, const Func& func)
{
return func.packetOp(
redux_vec_unroller<Func, Derived, Start, HalfLength>::run(mat,func),
redux_vec_unroller<Func, Derived, Start+HalfLength, Length-HalfLength>::run(mat,func) );
}
};
template<typename Func, typename Derived, int Start>
struct redux_vec_unroller<Func, Derived, Start, 1>
{
enum {
index = Start * redux_traits<Func, Derived>::PacketSize,
outer = index / int(Derived::InnerSizeAtCompileTime),
inner = index % int(Derived::InnerSizeAtCompileTime),
alignment = Derived::Alignment
};
typedef typename Derived::Scalar Scalar;
typedef typename redux_traits<Func, Derived>::PacketType PacketScalar;
static EIGEN_STRONG_INLINE PacketScalar run(const Derived &mat, const Func&)
{
return mat.template packetByOuterInner<alignment,PacketScalar>(outer, inner);
}
};
/***************************************************************************
* Part 3 : implementation of all cases
***************************************************************************/
template<typename Func, typename Derived,
int Traversal = redux_traits<Func, Derived>::Traversal,
int Unrolling = redux_traits<Func, Derived>::Unrolling
>
struct redux_impl;
template<typename Func, typename Derived>
struct redux_impl<Func, Derived, DefaultTraversal, NoUnrolling>
{
typedef typename Derived::Scalar Scalar;
EIGEN_DEVICE_FUNC
static EIGEN_STRONG_INLINE Scalar run(const Derived &mat, const Func& func)
{
eigen_assert(mat.rows()>0 && mat.cols()>0 && "you are using an empty matrix");
Scalar res;
res = mat.coeffByOuterInner(0, 0);
for(Index i = 1; i < mat.innerSize(); ++i)
res = func(res, mat.coeffByOuterInner(0, i));
for(Index i = 1; i < mat.outerSize(); ++i)
for(Index j = 0; j < mat.innerSize(); ++j)
res = func(res, mat.coeffByOuterInner(i, j));
return res;
}
};
template<typename Func, typename Derived>
struct redux_impl<Func,Derived, DefaultTraversal, CompleteUnrolling>
: public redux_novec_unroller<Func,Derived, 0, Derived::SizeAtCompileTime>
{};
template<typename Func, typename Derived>
struct redux_impl<Func, Derived, LinearVectorizedTraversal, NoUnrolling>
{
typedef typename Derived::Scalar Scalar;
typedef typename redux_traits<Func, Derived>::PacketType PacketScalar;
static Scalar run(const Derived &mat, const Func& func)
{
const Index size = mat.size();
const Index packetSize = redux_traits<Func, Derived>::PacketSize;
const int packetAlignment = unpacket_traits<PacketScalar>::alignment;
enum {
alignment0 = (bool(Derived::Flags & DirectAccessBit) && bool(packet_traits<Scalar>::AlignedOnScalar)) ? int(packetAlignment) : int(Unaligned),
alignment = EIGEN_PLAIN_ENUM_MAX(alignment0, Derived::Alignment)
};
const Index alignedStart = internal::first_default_aligned(mat.nestedExpression());
const Index alignedSize2 = ((size-alignedStart)/(2*packetSize))*(2*packetSize);
const Index alignedSize = ((size-alignedStart)/(packetSize))*(packetSize);
const Index alignedEnd2 = alignedStart + alignedSize2;
const Index alignedEnd = alignedStart + alignedSize;
Scalar res;
if(alignedSize)
{
PacketScalar packet_res0 = mat.template packet<alignment,PacketScalar>(alignedStart);
if(alignedSize>packetSize) // we have at least two packets to partly unroll the loop
{
PacketScalar packet_res1 = mat.template packet<alignment,PacketScalar>(alignedStart+packetSize);
for(Index index = alignedStart + 2*packetSize; index < alignedEnd2; index += 2*packetSize)
{
packet_res0 = func.packetOp(packet_res0, mat.template packet<alignment,PacketScalar>(index));
packet_res1 = func.packetOp(packet_res1, mat.template packet<alignment,PacketScalar>(index+packetSize));
}
packet_res0 = func.packetOp(packet_res0,packet_res1);
if(alignedEnd>alignedEnd2)
packet_res0 = func.packetOp(packet_res0, mat.template packet<alignment,PacketScalar>(alignedEnd2));
}
res = func.predux(packet_res0);
for(Index index = 0; index < alignedStart; ++index)
res = func(res,mat.coeff(index));
for(Index index = alignedEnd; index < size; ++index)
res = func(res,mat.coeff(index));
}
else // too small to vectorize anything.
// since this is dynamic-size hence inefficient anyway for such small sizes, don't try to optimize.
{
res = mat.coeff(0);
for(Index index = 1; index < size; ++index)
res = func(res,mat.coeff(index));
}
return res;
}
};
// NOTE: for SliceVectorizedTraversal we simply bypass unrolling
template<typename Func, typename Derived, int Unrolling>
struct redux_impl<Func, Derived, SliceVectorizedTraversal, Unrolling>
{
typedef typename Derived::Scalar Scalar;
typedef typename redux_traits<Func, Derived>::PacketType PacketType;
EIGEN_DEVICE_FUNC static Scalar run(const Derived &mat, const Func& func)
{
eigen_assert(mat.rows()>0 && mat.cols()>0 && "you are using an empty matrix");
const Index innerSize = mat.innerSize();
const Index outerSize = mat.outerSize();
enum {
packetSize = redux_traits<Func, Derived>::PacketSize
};
const Index packetedInnerSize = ((innerSize)/packetSize)*packetSize;
Scalar res;
if(packetedInnerSize)
{
PacketType packet_res = mat.template packet<Unaligned,PacketType>(0,0);
for(Index j=0; j<outerSize; ++j)
for(Index i=(j==0?packetSize:0); i<packetedInnerSize; i+=Index(packetSize))
packet_res = func.packetOp(packet_res, mat.template packetByOuterInner<Unaligned,PacketType>(j,i));
res = func.predux(packet_res);
for(Index j=0; j<outerSize; ++j)
for(Index i=packetedInnerSize; i<innerSize; ++i)
res = func(res, mat.coeffByOuterInner(j,i));
}
else // too small to vectorize anything.
// since this is dynamic-size hence inefficient anyway for such small sizes, don't try to optimize.
{
res = redux_impl<Func, Derived, DefaultTraversal, NoUnrolling>::run(mat, func);
}
return res;
}
};
template<typename Func, typename Derived>
struct redux_impl<Func, Derived, LinearVectorizedTraversal, CompleteUnrolling>
{
typedef typename Derived::Scalar Scalar;
typedef typename redux_traits<Func, Derived>::PacketType PacketScalar;
enum {
PacketSize = redux_traits<Func, Derived>::PacketSize,
Size = Derived::SizeAtCompileTime,
VectorizedSize = (Size / PacketSize) * PacketSize
};
EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Scalar run(const Derived &mat, const Func& func)
{
eigen_assert(mat.rows()>0 && mat.cols()>0 && "you are using an empty matrix");
if (VectorizedSize > 0) {
Scalar res = func.predux(redux_vec_unroller<Func, Derived, 0, Size / PacketSize>::run(mat,func));
if (VectorizedSize != Size)
res = func(res,redux_novec_unroller<Func, Derived, VectorizedSize, Size-VectorizedSize>::run(mat,func));
return res;
}
else {
return redux_novec_unroller<Func, Derived, 0, Size>::run(mat,func);
}
}
};
// evaluator adaptor
template<typename _XprType>
class redux_evaluator
{
public:
typedef _XprType XprType;
EIGEN_DEVICE_FUNC explicit redux_evaluator(const XprType &xpr) : m_evaluator(xpr), m_xpr(xpr) {}
typedef typename XprType::Scalar Scalar;
typedef typename XprType::CoeffReturnType CoeffReturnType;
typedef typename XprType::PacketScalar PacketScalar;
typedef typename XprType::PacketReturnType PacketReturnType;
enum {
MaxRowsAtCompileTime = XprType::MaxRowsAtCompileTime,
MaxColsAtCompileTime = XprType::MaxColsAtCompileTime,
// TODO we should not remove DirectAccessBit and rather find an elegant way to query the alignment offset at runtime from the evaluator
Flags = evaluator<XprType>::Flags & ~DirectAccessBit,
IsRowMajor = XprType::IsRowMajor,
SizeAtCompileTime = XprType::SizeAtCompileTime,
InnerSizeAtCompileTime = XprType::InnerSizeAtCompileTime,
CoeffReadCost = evaluator<XprType>::CoeffReadCost,
Alignment = evaluator<XprType>::Alignment
};
EIGEN_DEVICE_FUNC Index rows() const { return m_xpr.rows(); }
EIGEN_DEVICE_FUNC Index cols() const { return m_xpr.cols(); }
EIGEN_DEVICE_FUNC Index size() const { return m_xpr.size(); }
EIGEN_DEVICE_FUNC Index innerSize() const { return m_xpr.innerSize(); }
EIGEN_DEVICE_FUNC Index outerSize() const { return m_xpr.outerSize(); }
EIGEN_DEVICE_FUNC
CoeffReturnType coeff(Index row, Index col) const
{ return m_evaluator.coeff(row, col); }
EIGEN_DEVICE_FUNC
CoeffReturnType coeff(Index index) const
{ return m_evaluator.coeff(index); }
template<int LoadMode, typename PacketType>
PacketType packet(Index row, Index col) const
{ return m_evaluator.template packet<LoadMode,PacketType>(row, col); }
template<int LoadMode, typename PacketType>
PacketType packet(Index index) const
{ return m_evaluator.template packet<LoadMode,PacketType>(index); }
EIGEN_DEVICE_FUNC
CoeffReturnType coeffByOuterInner(Index outer, Index inner) const
{ return m_evaluator.coeff(IsRowMajor ? outer : inner, IsRowMajor ? inner : outer); }
template<int LoadMode, typename PacketType>
PacketType packetByOuterInner(Index outer, Index inner) const
{ return m_evaluator.template packet<LoadMode,PacketType>(IsRowMajor ? outer : inner, IsRowMajor ? inner : outer); }
const XprType & nestedExpression() const { return m_xpr; }
protected:
internal::evaluator<XprType> m_evaluator;
const XprType &m_xpr;
};
} // end namespace internal
/***************************************************************************
* Part 4 : public API
***************************************************************************/
/** \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 associative operator. Both current C++98 and C++11 functor styles are handled.
*
* \sa DenseBase::sum(), DenseBase::minCoeff(), DenseBase::maxCoeff(), MatrixBase::colwise(), MatrixBase::rowwise()
*/
template<typename Derived>
template<typename Func>
typename internal::traits<Derived>::Scalar
DenseBase<Derived>::redux(const Func& func) const
{
eigen_assert(this->rows()>0 && this->cols()>0 && "you are using an empty matrix");
typedef typename internal::redux_evaluator<Derived> ThisEvaluator;
ThisEvaluator thisEval(derived());
return internal::redux_impl<Func, ThisEvaluator>::run(thisEval, func);
}
/** \returns the minimum of all coefficients of \c *this.
* \warning the result is undefined if \c *this contains NaN.
*/
template<typename Derived>
EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
DenseBase<Derived>::minCoeff() const
{
return derived().redux(Eigen::internal::scalar_min_op<Scalar,Scalar>());
}
/** \returns the maximum of all coefficients of \c *this.
* \warning the result is undefined if \c *this contains NaN.
*/
template<typename Derived>
EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
DenseBase<Derived>::maxCoeff() const
{
return derived().redux(Eigen::internal::scalar_max_op<Scalar,Scalar>());
}
/** \returns the sum of all coefficients of \c *this
*
* If \c *this is empty, then the value 0 is returned.
*
* \sa trace(), prod(), mean()
*/
template<typename Derived>
EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
DenseBase<Derived>::sum() const
{
if(SizeAtCompileTime==0 || (SizeAtCompileTime==Dynamic && size()==0))
return Scalar(0);
return derived().redux(Eigen::internal::scalar_sum_op<Scalar,Scalar>());
}
/** \returns the mean of all coefficients of *this
*
* \sa trace(), prod(), sum()
*/
template<typename Derived>
EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
DenseBase<Derived>::mean() const
{
#ifdef __INTEL_COMPILER
#pragma warning push
#pragma warning ( disable : 2259 )
#endif
return Scalar(derived().redux(Eigen::internal::scalar_sum_op<Scalar,Scalar>())) / Scalar(this->size());
#ifdef __INTEL_COMPILER
#pragma warning pop
#endif
}
/** \returns the product of all coefficients of *this
*
* Example: \include MatrixBase_prod.cpp
* Output: \verbinclude MatrixBase_prod.out
*
* \sa sum(), mean(), trace()
*/
template<typename Derived>
EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
DenseBase<Derived>::prod() const
{
if(SizeAtCompileTime==0 || (SizeAtCompileTime==Dynamic && size()==0))
return Scalar(1);
return derived().redux(Eigen::internal::scalar_product_op<Scalar>());
}
/** \returns the trace of \c *this, i.e. the sum of the coefficients on the main diagonal.
*
* \c *this can be any matrix, not necessarily square.
*
* \sa diagonal(), sum()
*/
template<typename Derived>
EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
MatrixBase<Derived>::trace() const
{
return derived().diagonal().sum();
}
} // end namespace Eigen
#endif // EIGEN_REDUX_H
external/eigen3/Eigen/src/Core/Ref.h 0 → 100644
View file @ a394b22a
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_REF_H
#define EIGEN_REF_H
namespace Eigen {
namespace internal {
template<typename _PlainObjectType, int _Options, typename _StrideType>
struct traits<Ref<_PlainObjectType, _Options, _StrideType> >
: public traits<Map<_PlainObjectType, _Options, _StrideType> >
{
typedef _PlainObjectType PlainObjectType;
typedef _StrideType StrideType;
enum {
Options = _Options,
Flags = traits<Map<_PlainObjectType, _Options, _StrideType> >::Flags | NestByRefBit,
Alignment = traits<Map<_PlainObjectType, _Options, _StrideType> >::Alignment
};
template<typename Derived> struct match {
enum {
HasDirectAccess = internal::has_direct_access<Derived>::ret,
StorageOrderMatch = PlainObjectType::IsVectorAtCompileTime || Derived::IsVectorAtCompileTime || ((PlainObjectType::Flags&RowMajorBit)==(Derived::Flags&RowMajorBit)),
InnerStrideMatch = int(StrideType::InnerStrideAtCompileTime)==int(Dynamic)
|| int(StrideType::InnerStrideAtCompileTime)==int(Derived::InnerStrideAtCompileTime)
|| (int(StrideType::InnerStrideAtCompileTime)==0 && int(Derived::InnerStrideAtCompileTime)==1),
OuterStrideMatch = Derived::IsVectorAtCompileTime
|| int(StrideType::OuterStrideAtCompileTime)==int(Dynamic) || int(StrideType::OuterStrideAtCompileTime)==int(Derived::OuterStrideAtCompileTime),
// NOTE, this indirection of evaluator<Derived>::Alignment is needed
// to workaround a very strange bug in MSVC related to the instantiation
// of has_*ary_operator in evaluator<CwiseNullaryOp>.
// This line is surprisingly very sensitive. For instance, simply adding parenthesis
// as "DerivedAlignment = (int(evaluator<Derived>::Alignment))," will make MSVC fail...
DerivedAlignment = int(evaluator<Derived>::Alignment),
AlignmentMatch = (int(traits<PlainObjectType>::Alignment)==int(Unaligned)) || (DerivedAlignment >= int(Alignment)), // FIXME the first condition is not very clear, it should be replaced by the required alignment
ScalarTypeMatch = internal::is_same<typename PlainObjectType::Scalar, typename Derived::Scalar>::value,
MatchAtCompileTime = HasDirectAccess && StorageOrderMatch && InnerStrideMatch && OuterStrideMatch && AlignmentMatch && ScalarTypeMatch
};
typedef typename internal::conditional<MatchAtCompileTime,internal::true_type,internal::false_type>::type type;
};
};
template<typename Derived>
struct traits<RefBase<Derived> > : public traits<Derived> {};
}
template<typename Derived> class RefBase
: public MapBase<Derived>
{
typedef typename internal::traits<Derived>::PlainObjectType PlainObjectType;
typedef typename internal::traits<Derived>::StrideType StrideType;
public:
typedef MapBase<Derived> Base;
EIGEN_DENSE_PUBLIC_INTERFACE(RefBase)
EIGEN_DEVICE_FUNC inline Index innerStride() const
{
return StrideType::InnerStrideAtCompileTime != 0 ? m_stride.inner() : 1;
}
EIGEN_DEVICE_FUNC inline Index outerStride() const
{
return StrideType::OuterStrideAtCompileTime != 0 ? m_stride.outer()
: IsVectorAtCompileTime ? this->size()
: int(Flags)&RowMajorBit ? this->cols()
: this->rows();
}
EIGEN_DEVICE_FUNC RefBase()
: Base(0,RowsAtCompileTime==Dynamic?0:RowsAtCompileTime,ColsAtCompileTime==Dynamic?0:ColsAtCompileTime),
// Stride<> does not allow default ctor for Dynamic strides, so let' initialize it with dummy values:
m_stride(StrideType::OuterStrideAtCompileTime==Dynamic?0:StrideType::OuterStrideAtCompileTime,
StrideType::InnerStrideAtCompileTime==Dynamic?0:StrideType::InnerStrideAtCompileTime)
{}
EIGEN_INHERIT_ASSIGNMENT_OPERATORS(RefBase)
protected:
typedef Stride<StrideType::OuterStrideAtCompileTime,StrideType::InnerStrideAtCompileTime> StrideBase;
template<typename Expression>
EIGEN_DEVICE_FUNC void construct(Expression& expr)
{
if(PlainObjectType::RowsAtCompileTime==1)
{
eigen_assert(expr.rows()==1 || expr.cols()==1);
::new (static_cast<Base*>(this)) Base(expr.data(), 1, expr.size());
}
else if(PlainObjectType::ColsAtCompileTime==1)
{
eigen_assert(expr.rows()==1 || expr.cols()==1);
::new (static_cast<Base*>(this)) Base(expr.data(), expr.size(), 1);
}
else
::new (static_cast<Base*>(this)) Base(expr.data(), expr.rows(), expr.cols());
if(Expression::IsVectorAtCompileTime && (!PlainObjectType::IsVectorAtCompileTime) && ((Expression::Flags&RowMajorBit)!=(PlainObjectType::Flags&RowMajorBit)))
::new (&m_stride) StrideBase(expr.innerStride(), StrideType::InnerStrideAtCompileTime==0?0:1);
else
::new (&m_stride) StrideBase(StrideType::OuterStrideAtCompileTime==0?0:expr.outerStride(),
StrideType::InnerStrideAtCompileTime==0?0:expr.innerStride());
}
StrideBase m_stride;
};
/** \class Ref
* \ingroup Core_Module
*
* \brief A matrix or vector expression mapping an existing expression
*
* \tparam PlainObjectType the equivalent matrix type of the mapped data
* \tparam Options specifies the pointer alignment in bytes. It can be: \c #Aligned128, , \c #Aligned64, \c #Aligned32, \c #Aligned16, \c #Aligned8 or \c #Unaligned.
* The default is \c #Unaligned.
* \tparam StrideType optionally specifies strides. By default, Ref implies a contiguous storage along the inner dimension (inner stride==1),
* but accepts a variable outer stride (leading dimension).
* This can be overridden by specifying strides.
* The type passed here must be a specialization of the Stride template, see examples below.
*
* This class provides a way to write non-template functions taking Eigen objects as parameters while limiting the number of copies.
* A Ref<> object can represent either a const expression or a l-value:
* \code
* // in-out argument:
* void foo1(Ref<VectorXf> x);
*
* // read-only const argument:
* void foo2(const Ref<const VectorXf>& x);
* \endcode
*
* In the in-out case, the input argument must satisfy the constraints of the actual Ref<> type, otherwise a compilation issue will be triggered.
* By default, a Ref<VectorXf> can reference any dense vector expression of float having a contiguous memory layout.
* Likewise, a Ref<MatrixXf> can reference any column-major dense matrix expression of float whose column's elements are contiguously stored with
* the possibility to have a constant space in-between each column, i.e. the inner stride must be equal to 1, but the outer stride (or leading dimension)
* can be greater than the number of rows.
*
* In the const case, if the input expression does not match the above requirement, then it is evaluated into a temporary before being passed to the function.
* Here are some examples:
* \code
* MatrixXf A;
* VectorXf a;
* foo1(a.head()); // OK
* foo1(A.col()); // OK
* foo1(A.row()); // Compilation error because here innerstride!=1
* foo2(A.row()); // Compilation error because A.row() is a 1xN object while foo2 is expecting a Nx1 object
* foo2(A.row().transpose()); // The row is copied into a contiguous temporary
* foo2(2*a); // The expression is evaluated into a temporary
* foo2(A.col().segment(2,4)); // No temporary
* \endcode
*
* The range of inputs that can be referenced without temporary can be enlarged using the last two template parameters.
* Here is an example accepting an innerstride!=1:
* \code
* // in-out argument:
* void foo3(Ref<VectorXf,0,InnerStride<> > x);
* foo3(A.row()); // OK
* \endcode
* The downside here is that the function foo3 might be significantly slower than foo1 because it won't be able to exploit vectorization, and will involve more
* expensive address computations even if the input is contiguously stored in memory. To overcome this issue, one might propose to overload internally calling a
* template function, e.g.:
* \code
* // in the .h:
* void foo(const Ref<MatrixXf>& A);
* void foo(const Ref<MatrixXf,0,Stride<> >& A);
*
* // in the .cpp:
* template<typename TypeOfA> void foo_impl(const TypeOfA& A) {
* ... // crazy code goes here
* }
* void foo(const Ref<MatrixXf>& A) { foo_impl(A); }
* void foo(const Ref<MatrixXf,0,Stride<> >& A) { foo_impl(A); }
* \endcode
*
*
* \sa PlainObjectBase::Map(), \ref TopicStorageOrders
*/
template<typename PlainObjectType, int Options, typename StrideType> class Ref
: public RefBase<Ref<PlainObjectType, Options, StrideType> >
{
private:
typedef internal::traits<Ref> Traits;
template<typename Derived>
EIGEN_DEVICE_FUNC inline Ref(const PlainObjectBase<Derived>& expr,
typename internal::enable_if<bool(Traits::template match<Derived>::MatchAtCompileTime),Derived>::type* = 0);
public:
typedef RefBase<Ref> Base;
EIGEN_DENSE_PUBLIC_INTERFACE(Ref)
#ifndef EIGEN_PARSED_BY_DOXYGEN
template<typename Derived>
EIGEN_DEVICE_FUNC inline Ref(PlainObjectBase<Derived>& expr,
typename internal::enable_if<bool(Traits::template match<Derived>::MatchAtCompileTime),Derived>::type* = 0)
{
EIGEN_STATIC_ASSERT(bool(Traits::template match<Derived>::MatchAtCompileTime), STORAGE_LAYOUT_DOES_NOT_MATCH);
Base::construct(expr.derived());
}
template<typename Derived>
EIGEN_DEVICE_FUNC inline Ref(const DenseBase<Derived>& expr,
typename internal::enable_if<bool(Traits::template match<Derived>::MatchAtCompileTime),Derived>::type* = 0)
#else
/** Implicit constructor from any dense expression */
template<typename Derived>
inline Ref(DenseBase<Derived>& expr)
#endif
{
EIGEN_STATIC_ASSERT(bool(internal::is_lvalue<Derived>::value), THIS_EXPRESSION_IS_NOT_A_LVALUE__IT_IS_READ_ONLY);
EIGEN_STATIC_ASSERT(bool(Traits::template match<Derived>::MatchAtCompileTime), STORAGE_LAYOUT_DOES_NOT_MATCH);
EIGEN_STATIC_ASSERT(!Derived::IsPlainObjectBase,THIS_EXPRESSION_IS_NOT_A_LVALUE__IT_IS_READ_ONLY);
Base::construct(expr.const_cast_derived());
}
EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Ref)
};
// this is the const ref version
template<typename TPlainObjectType, int Options, typename StrideType> class Ref<const TPlainObjectType, Options, StrideType>
: public RefBase<Ref<const TPlainObjectType, Options, StrideType> >
{
typedef internal::traits<Ref> Traits;
public:
typedef RefBase<Ref> Base;
EIGEN_DENSE_PUBLIC_INTERFACE(Ref)
template<typename Derived>
EIGEN_DEVICE_FUNC inline Ref(const DenseBase<Derived>& expr,
typename internal::enable_if<bool(Traits::template match<Derived>::ScalarTypeMatch),Derived>::type* = 0)
{
// std::cout << match_helper<Derived>::HasDirectAccess << "," << match_helper<Derived>::OuterStrideMatch << "," << match_helper<Derived>::InnerStrideMatch << "\n";
// std::cout << int(StrideType::OuterStrideAtCompileTime) << " - " << int(Derived::OuterStrideAtCompileTime) << "\n";
// std::cout << int(StrideType::InnerStrideAtCompileTime) << " - " << int(Derived::InnerStrideAtCompileTime) << "\n";
construct(expr.derived(), typename Traits::template match<Derived>::type());
}
EIGEN_DEVICE_FUNC inline Ref(const Ref& other) : Base(other) {
// copy constructor shall not copy the m_object, to avoid unnecessary malloc and copy
}
template<typename OtherRef>
EIGEN_DEVICE_FUNC inline Ref(const RefBase<OtherRef>& other) {
construct(other.derived(), typename Traits::template match<OtherRef>::type());
}
protected:
template<typename Expression>
EIGEN_DEVICE_FUNC void construct(const Expression& expr,internal::true_type)
{
Base::construct(expr);
}
template<typename Expression>
EIGEN_DEVICE_FUNC void construct(const Expression& expr, internal::false_type)
{
internal::call_assignment_no_alias(m_object,expr,internal::assign_op<Scalar,Scalar>());
Base::construct(m_object);
}
protected:
TPlainObjectType m_object;
};
} // end namespace Eigen
#endif // EIGEN_REF_H
external/eigen3/Eigen/src/Core/Replicate.h 0 → 100644
View file @ a394b22a
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_REPLICATE_H
#define EIGEN_REPLICATE_H
namespace Eigen {
namespace internal {
template<typename MatrixType,int RowFactor,int ColFactor>
struct traits<Replicate<MatrixType,RowFactor,ColFactor> >
: traits<MatrixType>
{
typedef typename MatrixType::Scalar Scalar;
typedef typename traits<MatrixType>::StorageKind StorageKind;
typedef typename traits<MatrixType>::XprKind XprKind;
typedef typename ref_selector<MatrixType>::type MatrixTypeNested;
typedef typename remove_reference<MatrixTypeNested>::type _MatrixTypeNested;
enum {
RowsAtCompileTime = RowFactor==Dynamic || int(MatrixType::RowsAtCompileTime)==Dynamic
? Dynamic
: RowFactor * MatrixType::RowsAtCompileTime,
ColsAtCompileTime = ColFactor==Dynamic || int(MatrixType::ColsAtCompileTime)==Dynamic
? Dynamic
: ColFactor * MatrixType::ColsAtCompileTime,
//FIXME we don't propagate the max sizes !!!
MaxRowsAtCompileTime = RowsAtCompileTime,
MaxColsAtCompileTime = ColsAtCompileTime,
IsRowMajor = MaxRowsAtCompileTime==1 && MaxColsAtCompileTime!=1 ? 1
: MaxColsAtCompileTime==1 && MaxRowsAtCompileTime!=1 ? 0
: (MatrixType::Flags & RowMajorBit) ? 1 : 0,
// FIXME enable DirectAccess with negative strides?
Flags = IsRowMajor ? RowMajorBit : 0
};
};
}
/**
* \class Replicate
* \ingroup Core_Module
*
* \brief Expression of the multiple replication of a matrix or vector
*
* \tparam MatrixType the type of the object we are replicating
* \tparam RowFactor number of repetitions at compile time along the vertical direction, can be Dynamic.
* \tparam ColFactor number of repetitions at compile time along the horizontal direction, can be Dynamic.
*
* This class represents an expression of the multiple replication of a matrix or vector.
* It is the return type of DenseBase::replicate() and most of the time
* this is the only way it is used.
*
* \sa DenseBase::replicate()
*/
template<typename MatrixType,int RowFactor,int ColFactor> class Replicate
: public internal::dense_xpr_base< Replicate<MatrixType,RowFactor,ColFactor> >::type
{
typedef typename internal::traits<Replicate>::MatrixTypeNested MatrixTypeNested;
typedef typename internal::traits<Replicate>::_MatrixTypeNested _MatrixTypeNested;
public:
typedef typename internal::dense_xpr_base<Replicate>::type Base;
EIGEN_DENSE_PUBLIC_INTERFACE(Replicate)
typedef typename internal::remove_all<MatrixType>::type NestedExpression;
template<typename OriginalMatrixType>
EIGEN_DEVICE_FUNC
inline explicit Replicate(const OriginalMatrixType& matrix)
: m_matrix(matrix), m_rowFactor(RowFactor), m_colFactor(ColFactor)
{
EIGEN_STATIC_ASSERT((internal::is_same<typename internal::remove_const<MatrixType>::type,OriginalMatrixType>::value),
THE_MATRIX_OR_EXPRESSION_THAT_YOU_PASSED_DOES_NOT_HAVE_THE_EXPECTED_TYPE)
eigen_assert(RowFactor!=Dynamic && ColFactor!=Dynamic);
}
template<typename OriginalMatrixType>
EIGEN_DEVICE_FUNC
inline Replicate(const OriginalMatrixType& matrix, Index rowFactor, Index colFactor)
: m_matrix(matrix), m_rowFactor(rowFactor), m_colFactor(colFactor)
{
EIGEN_STATIC_ASSERT((internal::is_same<typename internal::remove_const<MatrixType>::type,OriginalMatrixType>::value),
THE_MATRIX_OR_EXPRESSION_THAT_YOU_PASSED_DOES_NOT_HAVE_THE_EXPECTED_TYPE)
}
EIGEN_DEVICE_FUNC
inline Index rows() const { return m_matrix.rows() * m_rowFactor.value(); }
EIGEN_DEVICE_FUNC
inline Index cols() const { return m_matrix.cols() * m_colFactor.value(); }
EIGEN_DEVICE_FUNC
const _MatrixTypeNested& nestedExpression() const
{
return m_matrix;
}
protected:
MatrixTypeNested m_matrix;
const internal::variable_if_dynamic<Index, RowFactor> m_rowFactor;
const internal::variable_if_dynamic<Index, ColFactor> m_colFactor;
};
/**
* \return an expression of the replication of \c *this
*
* Example: \include MatrixBase_replicate.cpp
* Output: \verbinclude MatrixBase_replicate.out
*
* \sa VectorwiseOp::replicate(), DenseBase::replicate(Index,Index), class Replicate
*/
template<typename Derived>
template<int RowFactor, int ColFactor>
const Replicate<Derived,RowFactor,ColFactor>
DenseBase<Derived>::replicate() const
{
return Replicate<Derived,RowFactor,ColFactor>(derived());
}
/**
* \return an expression of the replication of each column (or row) of \c *this
*
* Example: \include DirectionWise_replicate_int.cpp
* Output: \verbinclude DirectionWise_replicate_int.out
*
* \sa VectorwiseOp::replicate(), DenseBase::replicate(), class Replicate
*/
template<typename ExpressionType, int Direction>
const typename VectorwiseOp<ExpressionType,Direction>::ReplicateReturnType
VectorwiseOp<ExpressionType,Direction>::replicate(Index factor) const
{
return typename VectorwiseOp<ExpressionType,Direction>::ReplicateReturnType
(_expression(),Direction==Vertical?factor:1,Direction==Horizontal?factor:1);
}
} // end namespace Eigen
#endif // EIGEN_REPLICATE_H
external/eigen3/Eigen/src/Core/ReturnByValue.h 0 → 100644
View file @ a394b22a
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2009-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_RETURNBYVALUE_H
#define EIGEN_RETURNBYVALUE_H
namespace Eigen {
namespace internal {
template<typename Derived>
struct traits<ReturnByValue<Derived> >
: public traits<typename traits<Derived>::ReturnType>
{
enum {
// We're disabling the DirectAccess because e.g. the constructor of
// the Block-with-DirectAccess expression requires to have a coeffRef method.
// Also, we don't want to have to implement the stride stuff.
Flags = (traits<typename traits<Derived>::ReturnType>::Flags
| EvalBeforeNestingBit) & ~DirectAccessBit
};
};
/* The ReturnByValue object doesn't even have a coeff() method.
* So the only way that nesting it in an expression can work, is by evaluating it into a plain matrix.
* So internal::nested always gives the plain return matrix type.
*
* FIXME: I don't understand why we need this specialization: isn't this taken care of by the EvalBeforeNestingBit ??
* Answer: EvalBeforeNestingBit should be deprecated since we have the evaluators
*/
template<typename Derived,int n,typename PlainObject>
struct nested_eval<ReturnByValue<Derived>, n, PlainObject>
{
typedef typename traits<Derived>::ReturnType type;
};
} // end namespace internal
/** \class ReturnByValue
* \ingroup Core_Module
*
*/
template<typename Derived> class ReturnByValue
: public internal::dense_xpr_base< ReturnByValue<Derived> >::type, internal::no_assignment_operator
{
public:
typedef typename internal::traits<Derived>::ReturnType ReturnType;
typedef typename internal::dense_xpr_base<ReturnByValue>::type Base;
EIGEN_DENSE_PUBLIC_INTERFACE(ReturnByValue)
template<typename Dest>
EIGEN_DEVICE_FUNC
inline void evalTo(Dest& dst) const
{ static_cast<const Derived*>(this)->evalTo(dst); }
EIGEN_DEVICE_FUNC inline Index rows() const { return static_cast<const Derived*>(this)->rows(); }
EIGEN_DEVICE_FUNC inline Index cols() const { return static_cast<const Derived*>(this)->cols(); }
#ifndef EIGEN_PARSED_BY_DOXYGEN
#define Unusable YOU_ARE_TRYING_TO_ACCESS_A_SINGLE_COEFFICIENT_IN_A_SPECIAL_EXPRESSION_WHERE_THAT_IS_NOT_ALLOWED_BECAUSE_THAT_WOULD_BE_INEFFICIENT
class Unusable{
Unusable(const Unusable&) {}
Unusable& operator=(const Unusable&) {return *this;}
};
const Unusable& coeff(Index) const { return *reinterpret_cast<const Unusable*>(this); }
const Unusable& coeff(Index,Index) const { return *reinterpret_cast<const Unusable*>(this); }
Unusable& coeffRef(Index) { return *reinterpret_cast<Unusable*>(this); }
Unusable& coeffRef(Index,Index) { return *reinterpret_cast<Unusable*>(this); }
#undef Unusable
#endif
};
template<typename Derived>
template<typename OtherDerived>
Derived& DenseBase<Derived>::operator=(const ReturnByValue<OtherDerived>& other)
{
other.evalTo(derived());
return derived();
}
namespace internal {
// Expression is evaluated in a temporary; default implementation of Assignment is bypassed so that
// when a ReturnByValue expression is assigned, the evaluator is not constructed.
// TODO: Finalize port to new regime; ReturnByValue should not exist in the expression world
template<typename Derived>
struct evaluator<ReturnByValue<Derived> >
: public evaluator<typename internal::traits<Derived>::ReturnType>
{
typedef ReturnByValue<Derived> XprType;
typedef typename internal::traits<Derived>::ReturnType PlainObject;
typedef evaluator<PlainObject> Base;
EIGEN_DEVICE_FUNC explicit evaluator(const XprType& xpr)
: m_result(xpr.rows(), xpr.cols())
{
::new (static_cast<Base*>(this)) Base(m_result);
xpr.evalTo(m_result);
}
protected:
PlainObject m_result;
};
} // end namespace internal
} // end namespace Eigen
#endif // EIGEN_RETURNBYVALUE_H
external/eigen3/Eigen/src/Core/Reverse.h 0 → 100644
View file @ a394b22a
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2009 Ricard Marxer <email@ricardmarxer.com>
// Copyright (C) 2009-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_REVERSE_H
#define EIGEN_REVERSE_H
namespace Eigen {
namespace internal {
template<typename MatrixType, int Direction>
struct traits<Reverse<MatrixType, Direction> >
: traits<MatrixType>
{
typedef typename MatrixType::Scalar Scalar;
typedef typename traits<MatrixType>::StorageKind StorageKind;
typedef typename traits<MatrixType>::XprKind XprKind;
typedef typename ref_selector<MatrixType>::type MatrixTypeNested;
typedef typename remove_reference<MatrixTypeNested>::type _MatrixTypeNested;
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
Flags = _MatrixTypeNested::Flags & (RowMajorBit | LvalueBit)
};
};
template<typename PacketType, bool ReversePacket> struct reverse_packet_cond
{
static inline PacketType run(const PacketType& x) { return preverse(x); }
};
template<typename PacketType> struct reverse_packet_cond<PacketType,false>
{
static inline PacketType run(const PacketType& x) { return x; }
};
} // end namespace internal
/** \class Reverse
* \ingroup Core_Module
*
* \brief Expression of the reverse of a vector or matrix
*
* \tparam MatrixType the type of the object of which we are taking the reverse
* \tparam Direction defines the direction of the reverse operation, can be Vertical, Horizontal, or BothDirections
*
* This class represents an expression of the reverse of a vector.
* It is the return type of MatrixBase::reverse() and VectorwiseOp::reverse()
* and most of the time this is the only way it is used.
*
* \sa MatrixBase::reverse(), VectorwiseOp::reverse()
*/
template<typename MatrixType, int Direction> class Reverse
: public internal::dense_xpr_base< Reverse<MatrixType, Direction> >::type
{
public:
typedef typename internal::dense_xpr_base<Reverse>::type Base;
EIGEN_DENSE_PUBLIC_INTERFACE(Reverse)
typedef typename internal::remove_all<MatrixType>::type NestedExpression;
using Base::IsRowMajor;
protected:
enum {
PacketSize = internal::packet_traits<Scalar>::size,
IsColMajor = !IsRowMajor,
ReverseRow = (Direction == Vertical) || (Direction == BothDirections),
ReverseCol = (Direction == Horizontal) || (Direction == BothDirections),
OffsetRow = ReverseRow && IsColMajor ? PacketSize : 1,
OffsetCol = ReverseCol && IsRowMajor ? PacketSize : 1,
ReversePacket = (Direction == BothDirections)
|| ((Direction == Vertical) && IsColMajor)
|| ((Direction == Horizontal) && IsRowMajor)
};
typedef internal::reverse_packet_cond<PacketScalar,ReversePacket> reverse_packet;
public:
EIGEN_DEVICE_FUNC explicit inline Reverse(const MatrixType& matrix) : m_matrix(matrix) { }
EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Reverse)
EIGEN_DEVICE_FUNC inline Index rows() const { return m_matrix.rows(); }
EIGEN_DEVICE_FUNC inline Index cols() const { return m_matrix.cols(); }
EIGEN_DEVICE_FUNC inline Index innerStride() const
{
return -m_matrix.innerStride();
}
EIGEN_DEVICE_FUNC const typename internal::remove_all<typename MatrixType::Nested>::type&
nestedExpression() const
{
return m_matrix;
}
protected:
typename MatrixType::Nested m_matrix;
};
/** \returns an expression of the reverse of *this.
*
* Example: \include MatrixBase_reverse.cpp
* Output: \verbinclude MatrixBase_reverse.out
*
*/
template<typename Derived>
inline typename DenseBase<Derived>::ReverseReturnType
DenseBase<Derived>::reverse()
{
return ReverseReturnType(derived());
}
//reverse const overload moved DenseBase.h due to a CUDA compiler bug
/** This is the "in place" version of reverse: it reverses \c *this.
*
* In most cases it is probably better to simply use the reversed expression
* of a matrix. However, when reversing the matrix data itself is really needed,
* then this "in-place" version is probably the right choice because it provides
* the following additional benefits:
* - less error prone: doing the same operation with .reverse() requires special care:
* \code m = m.reverse().eval(); \endcode
* - this API enables reverse operations without the need for a temporary
* - it allows future optimizations (cache friendliness, etc.)
*
* \sa VectorwiseOp::reverseInPlace(), reverse() */
template<typename Derived>
inline void DenseBase<Derived>::reverseInPlace()
{
if(cols()>rows())
{
Index half = cols()/2;
leftCols(half).swap(rightCols(half).reverse());
if((cols()%2)==1)
{
Index half2 = rows()/2;
col(half).head(half2).swap(col(half).tail(half2).reverse());
}
}
else
{
Index half = rows()/2;
topRows(half).swap(bottomRows(half).reverse());
if((rows()%2)==1)
{
Index half2 = cols()/2;
row(half).head(half2).swap(row(half).tail(half2).reverse());
}
}
}
namespace internal {
template<int Direction>
struct vectorwise_reverse_inplace_impl;
template<>
struct vectorwise_reverse_inplace_impl<Vertical>
{
template<typename ExpressionType>
static void run(ExpressionType &xpr)
{
Index half = xpr.rows()/2;
xpr.topRows(half).swap(xpr.bottomRows(half).colwise().reverse());
}
};
template<>
struct vectorwise_reverse_inplace_impl<Horizontal>
{
template<typename ExpressionType>
static void run(ExpressionType &xpr)
{
Index half = xpr.cols()/2;
xpr.leftCols(half).swap(xpr.rightCols(half).rowwise().reverse());
}
};
} // end namespace internal
/** This is the "in place" version of VectorwiseOp::reverse: it reverses each column or row of \c *this.
*
* In most cases it is probably better to simply use the reversed expression
* of a matrix. However, when reversing the matrix data itself is really needed,
* then this "in-place" version is probably the right choice because it provides
* the following additional benefits:
* - less error prone: doing the same operation with .reverse() requires special care:
* \code m = m.reverse().eval(); \endcode
* - this API enables reverse operations without the need for a temporary
*
* \sa DenseBase::reverseInPlace(), reverse() */
template<typename ExpressionType, int Direction>
void VectorwiseOp<ExpressionType,Direction>::reverseInPlace()
{
internal::vectorwise_reverse_inplace_impl<Direction>::run(_expression().const_cast_derived());
}
} // end namespace Eigen
#endif // EIGEN_REVERSE_H