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

merge with tip

Browse Source
This commit is contained in:
Thomas Capricelli 2009-09-10 02:37:08 +02:00
commit 72746838ad
91 changed files with 2602 additions and 905 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Eigen/CMakeLists.txt
View File

@ -1,4 +1,4 @@
set(Eigen_HEADERS Core LU Cholesky QR Geometry Sparse Array SVD LeastSquares QtAlignedMalloc StdVector Householder Jacobi)
set(Eigen_HEADERS Core LU Cholesky QR Geometry Sparse Array SVD LeastSquares QtAlignedMalloc StdVector Householder Jacobi Eigenvalues)
if(EIGEN_BUILD_LIB)
set(Eigen_SRCS

74
Eigen/Eigenvalues Normal file
View File

@ -0,0 +1,74 @@
#ifndef EIGEN_EIGENVALUES_MODULE_H
#define EIGEN_EIGENVALUES_MODULE_H
#include "Core"
#include "src/Core/util/DisableMSVCWarnings.h"
#include "Cholesky"
#include "Jacobi"
#include "Householder"
// Note that EIGEN_HIDE_HEAVY_CODE has to be defined per module
#if (defined EIGEN_EXTERN_INSTANTIATIONS) && (EIGEN_EXTERN_INSTANTIATIONS>=2)
#ifndef EIGEN_HIDE_HEAVY_CODE
#define EIGEN_HIDE_HEAVY_CODE
#endif
#elif defined EIGEN_HIDE_HEAVY_CODE
#undef EIGEN_HIDE_HEAVY_CODE
#endif
namespace Eigen {
/** \defgroup Eigenvalues_Module Eigenvalues module
*
* \nonstableyet
*
* This module mainly provides various eigenvalue solvers.
* This module also provides some MatrixBase methods, including:
* - MatrixBase::eigenvalues(),
* - MatrixBase::operatorNorm()
*
* \code
* #include <Eigen/Eigenvalues>
* \endcode
*/
#include "src/Eigenvalues/Tridiagonalization.h"
#include "src/Eigenvalues/EigenSolver.h"
#include "src/Eigenvalues/SelfAdjointEigenSolver.h"
#include "src/Eigenvalues/HessenbergDecomposition.h"
#include "src/Eigenvalues/ComplexSchur.h"
#include "src/Eigenvalues/ComplexEigenSolver.h"
// declare all classes for a given matrix type
#define EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(MATRIXTYPE,PREFIX) \
PREFIX template class Tridiagonalization<MATRIXTYPE>; \
PREFIX template class HessenbergDecomposition<MATRIXTYPE>; \
PREFIX template class SelfAdjointEigenSolver<MATRIXTYPE>
// removed because it does not support complex yet
// PREFIX template class EigenSolver<MATRIXTYPE>
// declare all class for all types
#define EIGEN_EIGENVALUES_MODULE_INSTANTIATE(PREFIX) \
EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(Matrix2f,PREFIX); \
EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(Matrix2d,PREFIX); \
EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(Matrix3f,PREFIX); \
EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(Matrix3d,PREFIX); \
EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(Matrix4f,PREFIX); \
EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(Matrix4d,PREFIX); \
EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(MatrixXf,PREFIX); \
EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(MatrixXd,PREFIX); \
EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(MatrixXcf,PREFIX); \
EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(MatrixXcd,PREFIX)
#ifdef EIGEN_EXTERN_INSTANTIATIONS
EIGEN_EIGENVALUES_MODULE_INSTANTIATE(extern);
#endif // EIGEN_EXTERN_INSTANTIATIONS
} // namespace Eigen
#include "src/Core/util/EnableMSVCWarnings.h"
#endif // EIGEN_EIGENVALUES_MODULE_H

6
Eigen/Jacobi
View File

@ -8,11 +8,15 @@
namespace Eigen {
/** \defgroup Jacobi_Module Jacobi module
* This module provides Jacobi rotations.
* This module provides Jacobi and Givens rotations.
*
* \code
* #include <Eigen/Jacobi>
* \endcode
*
* In addition to listed classes, it defines the two following MatrixBase methods to apply a Jacobi or Givens rotation:
* - MatrixBase::applyOnTheLeft()
* - MatrixBase::applyOnTheRight().
*/
#include "src/Jacobi/Jacobi.h"

2
Eigen/LeastSquares
View File

@ -5,7 +5,7 @@
#include "src/Core/util/DisableMSVCWarnings.h"
#include "QR"
#include "Eigenvalues"
#include "Geometry"
namespace Eigen {

17
Eigen/QR
View File

@ -24,11 +24,9 @@ namespace Eigen {
*
* \nonstableyet
*
* This module mainly provides QR decomposition and an eigen value solver.
* This module provides various QR decompositions
* This module also provides some MatrixBase methods, including:
* - MatrixBase::qr(),
* - MatrixBase::eigenvalues(),
* - MatrixBase::operatorNorm()
*
* \code
* #include <Eigen/QR>
@ -38,20 +36,10 @@ namespace Eigen {
#include "src/QR/HouseholderQR.h"
#include "src/QR/FullPivotingHouseholderQR.h"
#include "src/QR/ColPivotingHouseholderQR.h"
#include "src/QR/Tridiagonalization.h"
#include "src/QR/EigenSolver.h"
#include "src/QR/SelfAdjointEigenSolver.h"
#include "src/QR/HessenbergDecomposition.h"
// declare all classes for a given matrix type
#define EIGEN_QR_MODULE_INSTANTIATE_TYPE(MATRIXTYPE,PREFIX) \
PREFIX template class HouseholderQR<MATRIXTYPE>; \
PREFIX template class Tridiagonalization<MATRIXTYPE>; \
PREFIX template class HessenbergDecomposition<MATRIXTYPE>; \
PREFIX template class SelfAdjointEigenSolver<MATRIXTYPE>
// removed because it does not support complex yet
// PREFIX template class EigenSolver<MATRIXTYPE>
// declare all class for all types
#define EIGEN_QR_MODULE_INSTANTIATE(PREFIX) \
@ -74,4 +62,7 @@ namespace Eigen {
#include "src/Core/util/EnableMSVCWarnings.h"
// FIXME for compatibility we include Eigenvalues here:
#include "Eigenvalues"
#endif // EIGEN_QR_MODULE_H

6
Eigen/QtAlignedMalloc
View File

@ -6,8 +6,10 @@
#if (!EIGEN_MALLOC_ALREADY_ALIGNED)
#include "src/Core/util/DisableMSVCWarnings.h"
void *qMalloc(size_t size)
{std::cerr << "ok\n";
{
return Eigen::ei_aligned_malloc(size);
}
@ -24,6 +26,8 @@ void *qRealloc(void *ptr, size_t size)
return newPtr;
}
#include "src/Core/util/EnableMSVCWarnings.h"
#endif
#endif // EIGEN_QTMALLOC_MODULE_H

4
Eigen/SVD
View File

@ -1,7 +1,7 @@
#ifndef EIGEN_SVD_MODULE_H
#define EIGEN_SVD_MODULE_H
#include "Core"
#include "QR"
#include "Householder"
#include "Jacobi"
@ -23,7 +23,7 @@ namespace Eigen {
*/
#include "src/SVD/SVD.h"
#include "src/SVD/JacobiSquareSVD.h"
#include "src/SVD/JacobiSVD.h"
} // namespace Eigen

8
Eigen/src/Array/BooleanRedux.h
View File

@ -78,10 +78,8 @@ struct ei_any_unroller<Derived, Dynamic>
};
/** \array_module
*
* \returns true if all coefficients are true
*
* \addexample CwiseAll \label How to check whether a point is inside a box (using operator< and all())
* \returns true if all coefficients are true
*
* Example: \include MatrixBase_all.cpp
* Output: \verbinclude MatrixBase_all.out
@ -107,7 +105,7 @@ inline bool MatrixBase<Derived>::all() const
}
/** \array_module
*
*
* \returns true if at least one coefficient is true
*
* \sa MatrixBase::all()
@ -131,7 +129,7 @@ inline bool MatrixBase<Derived>::any() const
}
/** \array_module
*
*
* \returns the number of coefficients which evaluate to true
*
* \sa MatrixBase::all(), MatrixBase::any()

10
Eigen/src/Array/Random.h
View File

@ -34,7 +34,7 @@ struct ei_functor_traits<ei_scalar_random_op<Scalar> >
{ enum { Cost = 5 * NumTraits<Scalar>::MulCost, PacketAccess = false, IsRepeatable = false }; };
/** \array_module
*
*
* \returns a random matrix expression
*
* The parameters \a rows and \a cols are the number of rows and of columns of
@ -44,8 +44,6 @@ struct ei_functor_traits<ei_scalar_random_op<Scalar> >
* it is redundant to pass \a rows and \a cols as arguments, so Random() should be used
* instead.
*
* \addexample RandomExample \label How to create a matrix with random coefficients
*
* Example: \include MatrixBase_random_int_int.cpp
* Output: \verbinclude MatrixBase_random_int_int.out
*
@ -63,7 +61,7 @@ MatrixBase<Derived>::Random(int rows, int cols)
}
/** \array_module
*
*
* \returns a random vector expression
*
* The parameter \a size is the size of the returned vector.
@ -92,7 +90,7 @@ MatrixBase<Derived>::Random(int size)
}
/** \array_module
*
*
* \returns a fixed-size random matrix or vector expression
*
* This variant is only for fixed-size MatrixBase types. For dynamic-size types, you
@ -115,7 +113,7 @@ MatrixBase<Derived>::Random()
}
/** \array_module
*
*
* Sets all coefficients in this expression to random values.
*
* Example: \include MatrixBase_setRandom.cpp

1
Eigen/src/CMakeLists.txt
View File

@ -9,3 +9,4 @@ ADD_SUBDIRECTORY(LeastSquares)
ADD_SUBDIRECTORY(Sparse)
ADD_SUBDIRECTORY(Jacobi)
ADD_SUBDIRECTORY(Householder)
ADD_SUBDIRECTORY(Eigenvalues)

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

@ -218,8 +218,8 @@ LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a)
int endSize = size - j - 1;
if (endSize > 0) {
_temporary.end(endSize) = ( m_matrix.block(j+1,0, endSize, j)
* m_matrix.col(j).start(j).conjugate() ).lazy();
_temporary.end(endSize).noalias() = m_matrix.block(j+1,0, endSize, j)
* m_matrix.col(j).start(j).conjugate();
m_matrix.row(j).end(endSize) = m_matrix.row(j).end(endSize).conjugate()
- _temporary.end(endSize).transpose();

2
Eigen/src/Cholesky/LLT.h
View File

@ -141,7 +141,7 @@ template<> struct ei_llt_inplace<LowerTriangular>
if (x<=RealScalar(0))
return false;
mat.coeffRef(k,k) = x = ei_sqrt(x);
if (k>0 && rs>0) A21 -= (A20 * A10.adjoint()).lazy();
if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
if (rs>0) A21 *= RealScalar(1)/x;
}
return true;

26
Eigen/src/Core/Block.h
View File

@ -201,6 +201,13 @@ template<typename MatrixType, int BlockRows, int BlockCols, int PacketAccess, in
m_startCol.value() + (RowsAtCompileTime == 1 ? index : 0), x);
}
#ifdef EIGEN_PARSED_BY_DOXYGEN
/** \sa MapBase::data() */
inline const Scalar* data() const;
/** \sa MapBase::stride() */
inline int stride() const;
#endif
protected:
const typename MatrixType::Nested m_matrix;
@ -269,7 +276,14 @@ class Block<MatrixType,BlockRows,BlockCols,PacketAccess,HasDirectAccess>
&& startCol >= 0 && blockCols >= 0 && startCol + blockCols <= matrix.cols());
}
inline int stride(void) const { return m_matrix.stride(); }
/** \sa MapBase::stride() */
inline int stride() const
{
return ((!Base::IsVectorAtCompileTime)
|| (BlockRows==1 && ((Flags&RowMajorBit)==0))
|| (BlockCols==1 && ((Flags&RowMajorBit)==RowMajorBit)))
? m_matrix.stride() : 1;
}
#ifndef __SUNPRO_CC
// FIXME sunstudio is not friendly with the above friend...
@ -294,8 +308,6 @@ class Block<MatrixType,BlockRows,BlockCols,PacketAccess,HasDirectAccess>
* \param blockRows the number of rows in the block
* \param blockCols the number of columns in the block
*
* \addexample BlockIntIntIntInt \label How to reference a sub-matrix (dynamic-size)
*
* Example: \include MatrixBase_block_int_int_int_int.cpp
* Output: \verbinclude MatrixBase_block_int_int_int_int.out
*
@ -327,8 +339,6 @@ inline const typename BlockReturnType<Derived>::Type MatrixBase<Derived>
* \param cRows the number of rows in the corner
* \param cCols the number of columns in the corner
*
* \addexample BlockCornerDynamicSize \label How to reference a sub-corner of a matrix
*
* Example: \include MatrixBase_corner_enum_int_int.cpp
* Output: \verbinclude MatrixBase_corner_enum_int_int.out
*
@ -438,8 +448,6 @@ MatrixBase<Derived>::corner(CornerType type) const
* \param startRow the first row in the block
* \param startCol the first column in the block
*
* \addexample BlockSubMatrixFixedSize \label How to reference a sub-matrix (fixed-size)
*
* Example: \include MatrixBase_block_int_int.cpp
* Output: \verbinclude MatrixBase_block_int_int.out
*
@ -466,8 +474,6 @@ MatrixBase<Derived>::block(int startRow, int startCol) const
}
/** \returns an expression of the \a i-th column of *this. Note that the numbering starts at 0.
*
* \addexample BlockColumn \label How to reference a single column of a matrix
*
* Example: \include MatrixBase_col.cpp
* Output: \verbinclude MatrixBase_col.out
@ -489,8 +495,6 @@ MatrixBase<Derived>::col(int i) const
}
/** \returns an expression of the \a i-th row of *this. Note that the numbering starts at 0.
*
* \addexample BlockRow \label How to reference a single row of a matrix
*
* Example: \include MatrixBase_row.cpp
* Output: \verbinclude MatrixBase_row.out

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

@ -124,8 +124,6 @@ struct CommaInitializer
* The coefficients must be provided in a row major order and exactly match
* the size of the matrix. Otherwise an assertion is raised.
*
* \addexample CommaInit \label How to easily set all the coefficients of a matrix
*
* Example: \include MatrixBase_set.cpp
* Output: \verbinclude MatrixBase_set.out
*

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

@ -291,8 +291,6 @@ Cwise<ExpressionType>::max(const MatrixBase<OtherDerived> &other) const
* The template parameter \a CustomBinaryOp is the type of the functor
* of the custom operator (see class CwiseBinaryOp for an example)
*
* \addexample CustomCwiseBinaryFunctors \label How to use custom coeff wise binary functors
*
* Here is an example illustrating the use of custom functors:
* \include class_CwiseBinaryOp.cpp
* Output: \verbinclude class_CwiseBinaryOp.out

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

@ -317,8 +317,6 @@ Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::setConstant(int row
* it is redundant to pass \a rows and \a cols as arguments, so Zero() should be used
* instead.
*
* \addexample Zero \label How to take get a zero matrix
*
* Example: \include MatrixBase_zero_int_int.cpp
* Output: \verbinclude MatrixBase_zero_int_int.out
*
@ -448,8 +446,6 @@ Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::setZero(int rows, i
* it is redundant to pass \a rows and \a cols as arguments, so Ones() should be used
* instead.
*
* \addexample One \label How to get a matrix with all coefficients equal one
*
* Example: \include MatrixBase_ones_int_int.cpp
* Output: \verbinclude MatrixBase_ones_int_int.out
*
@ -576,8 +572,6 @@ Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::setOnes(int rows, i
* it is redundant to pass \a rows and \a cols as arguments, so Identity() should be used
* instead.
*
* \addexample Identity \label How to get an identity matrix
*
* Example: \include MatrixBase_identity_int_int.cpp
* Output: \verbinclude MatrixBase_identity_int_int.out
*

16
Eigen/src/Core/CwiseUnaryOp.h
View File

@ -109,8 +109,6 @@ class CwiseUnaryOp : ei_no_assignment_operator,
* The template parameter \a CustomUnaryOp is the type of the functor
* of the custom unary operator.
*
* \addexample CustomCwiseUnaryFunctors \label How to use custom coeff wise unary functors
*
* Example:
* \include class_CwiseUnaryOp.cpp
* Output: \verbinclude class_CwiseUnaryOp.out
@ -234,7 +232,7 @@ Cwise<ExpressionType>::log() const
}
/** \relates MatrixBase */
/** \returns an expression of \c *this scaled by the scalar factor \a scalar */
template<typename Derived>
EIGEN_STRONG_INLINE const typename MatrixBase<Derived>::ScalarMultipleReturnType
MatrixBase<Derived>::operator*(const Scalar& scalar) const
@ -243,7 +241,17 @@ MatrixBase<Derived>::operator*(const Scalar& scalar) const
(derived(), ei_scalar_multiple_op<Scalar>(scalar));
}
/** \relates MatrixBase */
/** Overloaded for efficient real matrix times complex scalar value */
template<typename Derived>
EIGEN_STRONG_INLINE const CwiseUnaryOp<ei_scalar_multiple2_op<typename ei_traits<Derived>::Scalar,
std::complex<typename ei_traits<Derived>::Scalar> >, Derived>
MatrixBase<Derived>::operator*(const std::complex<Scalar>& scalar) const
{
return CwiseUnaryOp<ei_scalar_multiple2_op<Scalar,std::complex<Scalar> >, Derived>
(*static_cast<const Derived*>(this), ei_scalar_multiple2_op<Scalar,std::complex<Scalar> >(scalar));
}
/** \returns an expression of \c *this divided by the scalar value \a scalar */
template<typename Derived>
EIGEN_STRONG_INLINE const CwiseUnaryOp<ei_scalar_quotient1_op<typename ei_traits<Derived>::Scalar>, Derived>
MatrixBase<Derived>::operator/(const Scalar& scalar) const

5
Eigen/src/Core/CwiseUnaryView.h
View File

@ -53,8 +53,7 @@ struct ei_traits<CwiseUnaryView<ViewOp, MatrixType> >
};
template<typename ViewOp, typename MatrixType>
class CwiseUnaryView : ei_no_assignment_operator,
public MatrixBase<CwiseUnaryView<ViewOp, MatrixType> >
class CwiseUnaryView : public MatrixBase<CwiseUnaryView<ViewOp, MatrixType> >
{
public:
@ -99,8 +98,6 @@ class CwiseUnaryView : ei_no_assignment_operator,
* The template parameter \a CustomUnaryOp is the type of the functor
* of the custom unary operator.
*
* \addexample CustomCwiseUnaryFunctors \label How to use custom coeff wise unary functors
*
* Example:
* \include class_CwiseUnaryOp.cpp
* Output: \verbinclude class_CwiseUnaryOp.out

26
Eigen/src/Core/DiagonalMatrix.h
View File

@ -32,7 +32,7 @@ class DiagonalBase : public AnyMatrixBase<Derived>
public:
typedef typename ei_traits<Derived>::DiagonalVectorType DiagonalVectorType;
typedef typename DiagonalVectorType::Scalar Scalar;
enum {
RowsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
ColsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
@ -41,14 +41,14 @@ class DiagonalBase : public AnyMatrixBase<Derived>
IsVectorAtCompileTime = 0,
Flags = 0
};
typedef Matrix<Scalar, RowsAtCompileTime, ColsAtCompileTime, 0, MaxRowsAtCompileTime, MaxColsAtCompileTime> DenseMatrixType;
#ifndef EIGEN_PARSED_BY_DOXYGEN
inline const Derived& derived() const { return *static_cast<const Derived*>(this); }
inline Derived& derived() { return *static_cast<Derived*>(this); }
#endif // not EIGEN_PARSED_BY_DOXYGEN
DenseMatrixType toDenseMatrix() const { return derived(); }
template<typename DenseDerived>
void evalToDense(MatrixBase<DenseDerived> &other) const;
@ -64,7 +64,7 @@ class DiagonalBase : public AnyMatrixBase<Derived>
inline int rows() const { return diagonal().size(); }
inline int cols() const { return diagonal().size(); }
template<typename MatrixDerived>
const DiagonalProduct<MatrixDerived, Derived, DiagonalOnTheLeft>
operator*(const MatrixBase<MatrixDerived> &matrix) const;
@ -100,20 +100,20 @@ class DiagonalMatrix
: public DiagonalBase<DiagonalMatrix<_Scalar,_Size> >
{
public:
typedef typename ei_traits<DiagonalMatrix>::DiagonalVectorType DiagonalVectorType;
typedef const DiagonalMatrix& Nested;
typedef _Scalar Scalar;
protected:
DiagonalVectorType m_diagonal;
public:
inline const DiagonalVectorType& diagonal() const { return m_diagonal; }
inline DiagonalVectorType& diagonal() { return m_diagonal; }
/** Default constructor without initialization */
inline DiagonalMatrix() {}
@ -152,7 +152,7 @@ class DiagonalMatrix
m_diagonal = other.m_diagonal();
return *this;
}
inline void resize(int size) { m_diagonal.resize(size); }
inline void setZero() { m_diagonal.setZero(); }
inline void setZero(int size) { m_diagonal.setZero(size); }
@ -194,10 +194,10 @@ class DiagonalWrapper
public:
typedef _DiagonalVectorType DiagonalVectorType;
typedef DiagonalWrapper Nested;
inline DiagonalWrapper(const DiagonalVectorType& diagonal) : m_diagonal(diagonal) {}
const DiagonalVectorType& diagonal() const { return m_diagonal; }
protected:
const typename DiagonalVectorType::Nested m_diagonal;
};
@ -207,8 +207,6 @@ class DiagonalWrapper
*
* \only_for_vectors
*
* \addexample AsDiagonalExample \label How to build a diagonal matrix from a vector
*
* Example: \include MatrixBase_asDiagonal.cpp
* Output: \verbinclude MatrixBase_asDiagonal.out
*

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

@ -29,7 +29,7 @@
template<typename MatrixType, typename DiagonalType, int ProductOrder>
struct ei_traits<DiagonalProduct<MatrixType, DiagonalType, ProductOrder> >
{
typedef typename MatrixType::Scalar Scalar;
typedef typename ei_scalar_product_traits<typename MatrixType::Scalar, typename DiagonalType::Scalar>::ReturnType Scalar;
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
@ -62,7 +62,7 @@ class DiagonalProduct : ei_no_assignment_operator,
{
return m_diagonal.diagonal().coeff(ProductOrder == DiagonalOnTheLeft ? row : col) * m_matrix.coeff(row, col);
}
template<int LoadMode>
EIGEN_STRONG_INLINE PacketScalar packet(int row, int col) const
{
@ -72,7 +72,7 @@ class DiagonalProduct : ei_no_assignment_operator,
DiagonalVectorPacketLoadMode = (LoadMode == Aligned && ((InnerSize%16) == 0)) ? Aligned : Unaligned
};
const int indexInDiagonalVector = ProductOrder == DiagonalOnTheLeft ? row : col;
if((int(StorageOrder) == RowMajor && int(ProductOrder) == DiagonalOnTheLeft)
||(int(StorageOrder) == ColMajor && int(ProductOrder) == DiagonalOnTheRight))
{

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

@ -62,7 +62,21 @@ template<typename Derived> class MapBase
inline int rows() const { return m_rows.value(); }
inline int cols() const { return m_cols.value(); }
/** Returns the leading dimension (for matrices) or the increment (for vectors) to be used with data().
*
* More precisely:
* - for a column major matrix it returns the number of elements between two successive columns
* - for a row major matrix it returns the number of elements between two successive rows
* - for a vector it returns the number of elements between two successive coefficients
* This function has to be used together with the MapBase::data() function.
*
* \sa MapBase::data() */
inline int stride() const { return derived().stride(); }
/** Returns a pointer to the first coefficient of the matrix or vector.
* This function has to be used together with the stride() function.
*
* \sa MapBase::stride() */
inline const Scalar* data() const { return m_data; }
template<bool IsForceAligned,typename Dummy> struct force_aligned_impl {

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

@ -116,6 +116,7 @@ inline float ei_imag(float) { return 0.f; }
inline float ei_conj(float x) { return x; }
inline float ei_abs(float x) { return std::abs(x); }
inline float ei_abs2(float x) { return x*x; }
inline float ei_norm1(float x) { return ei_abs(x); }
inline float ei_sqrt(float x) { return std::sqrt(x); }
inline float ei_exp(float x) { return std::exp(x); }
inline float ei_log(float x) { return std::log(x); }
@ -164,6 +165,7 @@ inline double ei_imag(double) { return 0.; }
inline double ei_conj(double x) { return x; }
inline double ei_abs(double x) { return std::abs(x); }
inline double ei_abs2(double x) { return x*x; }
inline double ei_norm1(double x) { return ei_abs(x); }
inline double ei_sqrt(double x) { return std::sqrt(x); }
inline double ei_exp(double x) { return std::exp(x); }
inline double ei_log(double x) { return std::log(x); }
@ -212,6 +214,7 @@ inline float& ei_imag_ref(std::complex<float>& x) { return reinterpret_cast<floa
inline std::complex<float> ei_conj(const std::complex<float>& x) { return std::conj(x); }
inline float ei_abs(const std::complex<float>& x) { return std::abs(x); }
inline float ei_abs2(const std::complex<float>& x) { return std::norm(x); }
inline float ei_norm1(const std::complex<float> &x) { return(ei_abs(x.real()) + ei_abs(x.imag())); }
inline std::complex<float> ei_exp(std::complex<float> x) { return std::exp(x); }
inline std::complex<float> ei_sin(std::complex<float> x) { return std::sin(x); }
inline std::complex<float> ei_cos(std::complex<float> x) { return std::cos(x); }
@ -248,6 +251,7 @@ inline double& ei_imag_ref(std::complex<double>& x) { return reinterpret_cast<do
inline std::complex<double> ei_conj(const std::complex<double>& x) { return std::conj(x); }
inline double ei_abs(const std::complex<double>& x) { return std::abs(x); }
inline double ei_abs2(const std::complex<double>& x) { return std::norm(x); }
inline double ei_norm1(const std::complex<double> &x) { return(ei_abs(x.real()) + ei_abs(x.imag())); }
inline std::complex<double> ei_exp(std::complex<double> x) { return std::exp(x); }
inline std::complex<double> ei_sin(std::complex<double> x) { return std::sin(x); }
inline std::complex<double> ei_cos(std::complex<double> x) { return std::cos(x); }

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

@ -25,6 +25,7 @@
#ifndef EIGEN_MATRIX_H
#define EIGEN_MATRIX_H
template <typename Derived, typename OtherDerived, bool IsVector = static_cast<bool>(Derived::IsVectorAtCompileTime)> struct ei_conservative_resize_like_impl;
/** \class Matrix
*
@ -126,6 +127,7 @@ class Matrix
EIGEN_GENERIC_PUBLIC_INTERFACE(Matrix)
enum { Options = _Options };
typedef typename Base::PlainMatrixType PlainMatrixType;
friend class Eigen::Map<Matrix, Unaligned>;
typedef class Eigen::Map<Matrix, Unaligned> UnalignedMapType;
friend class Eigen::Map<Matrix, Aligned>;
@ -139,15 +141,27 @@ class Matrix
&& SizeAtCompileTime!=Dynamic && ((sizeof(Scalar)*SizeAtCompileTime)%16)==0 };
EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(NeedsToAlign)
Base& base() { return *static_cast<Base*>(this); }
const Base& base() const { return *static_cast<const Base*>(this); }
EIGEN_STRONG_INLINE int rows() const { return m_storage.rows(); }
EIGEN_STRONG_INLINE int cols() const { return m_storage.cols(); }
EIGEN_STRONG_INLINE int stride(void) const
/** Returns the leading dimension (for matrices) or the increment (for vectors) to be used with data().
*
* More precisely:
* - for a column major matrix it returns the number of elements between two successive columns
* - for a row major matrix it returns the number of elements between two successive rows
* - for a vector it returns the number of elements between two successive coefficients
* This function has to be used together with the MapBase::data() function.
*
* \sa Matrix::data() */
EIGEN_STRONG_INLINE int stride() const
{
if(Flags & RowMajorBit)
return m_storage.cols();
if(IsVectorAtCompileTime)
return 1;
else
return m_storage.rows();
return (Flags & RowMajorBit) ? m_storage.cols() : m_storage.rows();
}
EIGEN_STRONG_INLINE const Scalar& coeff(int row, int col) const
@ -295,7 +309,7 @@ class Matrix
*/
template<typename OtherDerived>
EIGEN_STRONG_INLINE void resizeLike(const MatrixBase<OtherDerived>& other)
{
{
if(RowsAtCompileTime == 1)
{
ei_assert(other.isVector());
@ -309,6 +323,51 @@ class Matrix
else resize(other.rows(), other.cols());
}
/** Resizes \c *this to a \a rows x \a cols matrix while leaving old values of *this untouched.
*
* This method is intended for dynamic-size matrices. If you only want to change the number
* of rows and/or of columns, you can use conservativeResize(NoChange_t, int),
* conservativeResize(int, NoChange_t).
*
* The top-left part of the resized matrix will be the same as the overlapping top-left corner
* of *this. In case values need to be appended to the matrix they will be uninitialized.
*/
EIGEN_STRONG_INLINE void conservativeResize(int rows, int cols)
{
conservativeResizeLike(PlainMatrixType(rows, cols));
}
EIGEN_STRONG_INLINE void conservativeResize(int rows, NoChange_t)
{
// Note: see the comment in conservativeResize(int,int)
conservativeResize(rows, cols());
}
EIGEN_STRONG_INLINE void conservativeResize(NoChange_t, int cols)
{
// Note: see the comment in conservativeResize(int,int)
conservativeResize(rows(), cols);
}
/** Resizes \c *this to a vector of length \a size while retaining old values of *this.
*
* \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_STRONG_INLINE void conservativeResize(int size)
{
conservativeResizeLike(PlainMatrixType(size));
}
template<typename OtherDerived>
EIGEN_STRONG_INLINE void conservativeResizeLike(const MatrixBase<OtherDerived>& other)
{
ei_conservative_resize_like_impl<Matrix, OtherDerived>::run(*this, other);
}
/** 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),
@ -329,7 +388,15 @@ class Matrix
*/
EIGEN_STRONG_INLINE Matrix& operator=(const Matrix& other)
{
return _set(other);
return _set(other);
}
/** \sa MatrixBase::lazyAssign() */
template<typename OtherDerived>
EIGEN_STRONG_INLINE Matrix& lazyAssign(const MatrixBase<OtherDerived>& other)
{
_resize_to_match(other);
return Base::lazyAssign(other.derived());
}
template<typename OtherDerived,typename OtherEvalType>
@ -470,13 +537,8 @@ class Matrix
/** Override MatrixBase::swap() since for dynamic-sized matrices of same type it is enough to swap the
* data pointers.
*/
inline void swap(Matrix& other)
{
if (Base::SizeAtCompileTime==Dynamic)
m_storage.swap(other.m_storage);
else
this->Base::swap(other);
}
template<typename OtherDerived>
void swap(const MatrixBase<OtherDerived>& other);
/** \name Map
* These are convenience functions returning Map objects. The Map() static functions return unaligned Map objects,
@ -635,8 +697,71 @@ class Matrix
m_storage.data()[0] = x;
m_storage.data()[1] = y;
}
template<typename MatrixType, typename OtherDerived, bool SwapPointers>
friend struct ei_matrix_swap_impl;
};
template <typename Derived, typename OtherDerived, bool IsVector>
struct ei_conservative_resize_like_impl
{
static void run(MatrixBase<Derived>& _this, const MatrixBase<OtherDerived>& other)
{
// Note: Here is space for improvement. Basically, for conservativeResize(int,int),
// neither RowsAtCompileTime or ColsAtCompileTime must be Dynamic. If only one of the
// dimensions is dynamic, one could use either conservativeResize(int rows, NoChange_t) or
// conservativeResize(NoChange_t, int 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)
typename MatrixBase<Derived>::PlainMatrixType tmp(other);
const int common_rows = std::min(tmp.rows(), _this.rows());
const int common_cols = std::min(tmp.cols(), _this.cols());
tmp.block(0,0,common_rows,common_cols) = _this.block(0,0,common_rows,common_cols);
_this.derived().swap(tmp);
}
};
template <typename Derived, typename OtherDerived>
struct ei_conservative_resize_like_impl<Derived,OtherDerived,true>
{
static void run(MatrixBase<Derived>& _this, const MatrixBase<OtherDerived>& other)
{
// segment(...) will check whether Derived/OtherDerived are vectors!
typename MatrixBase<Derived>::PlainMatrixType tmp(other);
const int common_size = std::min<int>(_this.size(),tmp.size());
tmp.segment(0,common_size) = _this.segment(0,common_size);
_this.derived().swap(tmp);
}
};
template<typename MatrixType, typename OtherDerived, bool SwapPointers>
struct ei_matrix_swap_impl
{
static inline void run(MatrixType& matrix, MatrixBase<OtherDerived>& other)
{
matrix.base().swap(other);
}
};
template<typename MatrixType, typename OtherDerived>
struct ei_matrix_swap_impl<MatrixType, OtherDerived, true>
{
static inline void run(MatrixType& matrix, MatrixBase<OtherDerived>& other)
{
matrix.m_storage.swap(other.derived().m_storage);
}
};
template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
template<typename OtherDerived>
inline void Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::swap(const MatrixBase<OtherDerived>& other)
{
enum { SwapPointers = ei_is_same_type<Matrix, OtherDerived>::ret && Base::SizeAtCompileTime==Dynamic };
ei_matrix_swap_impl<Matrix, OtherDerived, bool(SwapPointers)>::run(*this, *const_cast<MatrixBase<OtherDerived>*>(&other));
}
/** \defgroup matrixtypedefs Global matrix typedefs
*
* \ingroup Core_Module

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

@ -91,9 +91,8 @@ template<typename Derived> struct AnyMatrixBase
*/
template<typename Derived> class MatrixBase
#ifndef EIGEN_PARSED_BY_DOXYGEN
: public AnyMatrixBase<Derived>,
public ei_special_scalar_op_base<Derived,typename ei_traits<Derived>::Scalar,
typename NumTraits<typename ei_traits<Derived>::Scalar>::Real>
: public ei_special_scalar_op_base<Derived,typename ei_traits<Derived>::Scalar,
typename NumTraits<typename ei_traits<Derived>::Scalar>::Real>
#endif // not EIGEN_PARSED_BY_DOXYGEN
{
public:
@ -419,10 +418,17 @@ template<typename Derived> class MatrixBase
const CwiseUnaryOp<ei_scalar_quotient1_op<typename ei_traits<Derived>::Scalar>, Derived>
operator/(const Scalar& scalar) const;
inline friend const CwiseUnaryOp<ei_scalar_multiple_op<typename ei_traits<Derived>::Scalar>, Derived>
const CwiseUnaryOp<ei_scalar_multiple2_op<Scalar,std::complex<Scalar> >, Derived>
operator*(const std::complex<Scalar>& scalar) const;
inline friend const ScalarMultipleReturnType
operator*(const Scalar& scalar, const MatrixBase& matrix)
{ return matrix*scalar; }
inline friend const CwiseUnaryOp<ei_scalar_multiple2_op<Scalar,std::complex<Scalar> >, Derived>
operator*(const std::complex<Scalar>& scalar, const MatrixBase& matrix)
{ return matrix*scalar; }
template<typename OtherDerived>
const typename ProductReturnType<Derived,OtherDerived>::Type
operator*(const MatrixBase<OtherDerived> &other) const;
@ -803,11 +809,10 @@ template<typename Derived> class MatrixBase
///////// Jacobi module /////////
void applyJacobiOnTheLeft(int p, int q, Scalar c, Scalar s);
void applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s);
bool makeJacobi(int p, int q, Scalar *c, Scalar *s) const;
bool makeJacobiForAtA(int p, int q, Scalar *c, Scalar *s) const;
bool makeJacobiForAAt(int p, int q, Scalar *c, Scalar *s) const;
template<typename OtherScalar>
void applyOnTheLeft(int p, int q, const PlanarRotation<OtherScalar>& j);
template<typename OtherScalar>
void applyOnTheRight(int p, int q, const PlanarRotation<OtherScalar>& j);
#ifdef EIGEN_MATRIXBASE_PLUGIN
#include EIGEN_MATRIXBASE_PLUGIN

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

@ -63,7 +63,7 @@ template<typename Lhs, typename Rhs> struct ei_product_type
Cols = Rhs::ColsAtCompileTime,
Depth = EIGEN_ENUM_MIN(Lhs::ColsAtCompileTime,Rhs::RowsAtCompileTime)
};
// the splitting into different lines of code here, introducing the _select enums and the typedef below,
// is to work around an internal compiler error with gcc 4.1 and 4.2.
private:
@ -73,7 +73,7 @@ private:
depth_select = Depth>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD ? Large : (Depth==1 ? 1 : Small)
};
typedef ei_product_type_selector<rows_select, cols_select, depth_select> product_type_selector;
public:
enum {
value = product_type_selector::ret
@ -84,18 +84,18 @@ public:
* based on the three dimensions of the product.
* This is a compile time mapping from {1,Small,Large}^3 -> {product types} */
// FIXME I'm not sure the current mapping is the ideal one.
template<int Rows, int Cols> struct ei_product_type_selector<Rows,Cols,1> { enum { ret = OuterProduct }; };
template<int Depth> struct ei_product_type_selector<1,1,Depth> { enum { ret = InnerProduct }; };
template<> struct ei_product_type_selector<1,1,1> { enum { ret = InnerProduct }; };
template<> struct ei_product_type_selector<Small,1,Small> { enum { ret = UnrolledProduct }; };
template<> struct ei_product_type_selector<1,Small,Small> { enum { ret = UnrolledProduct }; };
template<int Rows, int Cols> struct ei_product_type_selector<Rows, Cols, 1> { enum { ret = OuterProduct }; };
template<int Depth> struct ei_product_type_selector<1, 1, Depth> { enum { ret = InnerProduct }; };
template<> struct ei_product_type_selector<1, 1, 1> { enum { ret = InnerProduct }; };
template<> struct ei_product_type_selector<Small,1, Small> { enum { ret = UnrolledProduct }; };
template<> struct ei_product_type_selector<1, Small,Small> { enum { ret = UnrolledProduct }; };
template<> struct ei_product_type_selector<Small,Small,Small> { enum { ret = UnrolledProduct }; };
template<> struct ei_product_type_selector<1,Large,Small> { enum { ret = GemvProduct }; };
template<> struct ei_product_type_selector<1,Large,Large> { enum { ret = GemvProduct }; };
template<> struct ei_product_type_selector<1,Small,Large> { enum { ret = GemvProduct }; };
template<> struct ei_product_type_selector<Large,1,Small> { enum { ret = GemvProduct }; };
template<> struct ei_product_type_selector<Large,1,Large> { enum { ret = GemvProduct }; };
template<> struct ei_product_type_selector<Small,1,Large> { enum { ret = GemvProduct }; };
template<> struct ei_product_type_selector<1, Large,Small> { enum { ret = GemvProduct }; };
template<> struct ei_product_type_selector<1, Large,Large> { enum { ret = GemvProduct }; };
template<> struct ei_product_type_selector<1, Small,Large> { enum { ret = GemvProduct }; };
template<> struct ei_product_type_selector<Large,1, Small> { enum { ret = GemvProduct }; };
template<> struct ei_product_type_selector<Large,1, Large> { enum { ret = GemvProduct }; };
template<> struct ei_product_type_selector<Small,1, Large> { enum { ret = GemvProduct }; };
template<> struct ei_product_type_selector<Small,Small,Large> { enum { ret = GemmProduct }; };
template<> struct ei_product_type_selector<Large,Small,Large> { enum { ret = GemmProduct }; };
template<> struct ei_product_type_selector<Small,Large,Large> { enum { ret = GemmProduct }; };
@ -144,7 +144,7 @@ struct ProductReturnType<Lhs,Rhs,UnrolledProduct>
***********************************************************************/
// FIXME : maybe the "inner product" could return a Scalar
// instead of a 1x1 matrix ??
// instead of a 1x1 matrix ??
// Pro: more natural for the user
// Cons: this could be a problem if in a meta unrolled algorithm a matrix-matrix
// product ends up to a row-vector times col-vector product... To tackle this use
@ -162,7 +162,11 @@ class GeneralProduct<Lhs, Rhs, InnerProduct>
public:
EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct)
GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {}
GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
{
EIGEN_STATIC_ASSERT((ei_is_same_type<typename Lhs::RealScalar, typename Rhs::RealScalar>::ret),
YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
}
EIGEN_STRONG_INLINE Scalar value() const
{
@ -197,11 +201,15 @@ class GeneralProduct<Lhs, Rhs, OuterProduct>
public:
EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct)
GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {}
GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
{
EIGEN_STATIC_ASSERT((ei_is_same_type<typename Lhs::RealScalar, typename Rhs::RealScalar>::ret),
YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
}
template<typename Dest> void scaleAndAddTo(Dest& dest, Scalar alpha) const
{
ei_outer_product_selector<Dest::Flags&RowMajorBit>::run(*this, dest, alpha);
ei_outer_product_selector<(int(Dest::Flags)&RowMajorBit) ? RowMajor : ColMajor>::run(*this, dest, alpha);
}
};
@ -209,6 +217,7 @@ template<> struct ei_outer_product_selector<ColMajor> {
template<typename ProductType, typename Dest>
EIGEN_DONT_INLINE static void run(const ProductType& prod, Dest& dest, typename ProductType::Scalar alpha) {
// FIXME make sure lhs is sequentially stored
// FIXME not very good if rhs is real and lhs complex while alpha is real too
const int cols = dest.cols();
for (int j=0; j<cols; ++j)
dest.col(j) += (alpha * prod.rhs().coeff(j)) * prod.lhs();
@ -219,6 +228,7 @@ template<> struct ei_outer_product_selector<RowMajor> {
template<typename ProductType, typename Dest>
EIGEN_DONT_INLINE static void run(const ProductType& prod, Dest& dest, typename ProductType::Scalar alpha) {
// FIXME make sure rhs is sequentially stored
// FIXME not very good if lhs is real and rhs complex while alpha is real too
const int rows = dest.rows();
for (int i=0; i<rows; ++i)
dest.row(i) += (alpha * prod.lhs().coeff(i)) * prod.rhs();
@ -251,7 +261,11 @@ class GeneralProduct<Lhs, Rhs, GemvProduct>
public:
EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct)
GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {}
GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
{
EIGEN_STATIC_ASSERT((ei_is_same_type<typename Lhs::Scalar, typename Rhs::Scalar>::ret),
YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
}
enum { Side = Lhs::IsVectorAtCompileTime ? OnTheLeft : OnTheRight };
typedef typename ei_meta_if<int(Side)==OnTheRight,_LhsNested,_RhsNested>::ret MatrixType;
@ -259,8 +273,8 @@ class GeneralProduct<Lhs, Rhs, GemvProduct>
template<typename Dest> void scaleAndAddTo(Dest& dst, Scalar alpha) const
{
ei_assert(m_lhs.rows() == dst.rows() && m_rhs.cols() == dst.cols());
ei_gemv_selector<Side,int(MatrixType::Flags)&RowMajorBit,
ei_blas_traits<MatrixType>::ActualAccess>::run(*this, dst, alpha);
ei_gemv_selector<Side,(int(MatrixType::Flags)&RowMajorBit) ? RowMajor : ColMajor,
bool(ei_blas_traits<MatrixType>::ActualAccess)>::run(*this, dst, alpha);
}
};
@ -339,7 +353,7 @@ template<> struct ei_gemv_selector<OnTheRight,RowMajor,true>
Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(prod.lhs())
* RhsBlasTraits::extractScalarFactor(prod.rhs());
enum {
DirectlyUseRhs = ((ei_packet_traits<Scalar>::size==1) || (_ActualRhsType::Flags&ActualPacketAccessBit))
&& (!(_ActualRhsType::Flags & RowMajorBit))

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

@ -33,7 +33,7 @@ struct ei_traits<ProductBase<Derived,_Lhs,_Rhs> >
{
typedef typename ei_cleantype<_Lhs>::type Lhs;
typedef typename ei_cleantype<_Rhs>::type Rhs;
typedef typename ei_traits<Lhs>::Scalar Scalar;
typedef typename ei_scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType Scalar;
enum {
RowsAtCompileTime = ei_traits<Lhs>::RowsAtCompileTime,
ColsAtCompileTime = ei_traits<Rhs>::ColsAtCompileTime,
@ -146,7 +146,7 @@ class ScaledProduct;
// functions of ProductBase, because, otherwise we would have to
// define all overloads defined in MatrixBase. Furthermore, Using
// "using Base::operator*" would not work with MSVC.
//
//
// Also note that here we accept any compatible scalar types
template<typename Derived,typename Lhs,typename Rhs>
const ScaledProduct<Derived>

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

@ -30,7 +30,7 @@ template<typename Lhs, typename Rhs,
int Side, // can be OnTheLeft/OnTheRight
int Unrolling = Rhs::IsVectorAtCompileTime && Rhs::SizeAtCompileTime <= 8 // FIXME
? CompleteUnrolling : NoUnrolling,
int StorageOrder = int(Lhs::Flags) & RowMajorBit,
int StorageOrder = (int(Lhs::Flags) & RowMajorBit) ? RowMajor : ColMajor,
int RhsCols = Rhs::ColsAtCompileTime
>
struct ei_triangular_solver_selector;
@ -157,7 +157,7 @@ struct ei_triangular_solver_selector<Lhs,Rhs,Side,Mode,NoUnrolling,StorageOrder,
{
const ActualLhsType actualLhs = LhsProductTraits::extract(lhs);
ei_triangular_solve_matrix<Scalar,Side,Mode,LhsProductTraits::NeedToConjugate,StorageOrder,
Rhs::Flags&RowMajorBit>
(Rhs::Flags&RowMajorBit) ? RowMajor : ColMajor>
::run(lhs.rows(), Side==OnTheLeft? rhs.cols() : rhs.rows(), &actualLhs.coeff(0,0), actualLhs.stride(), &rhs.coeffRef(0,0), rhs.stride());
}
};

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

@ -88,7 +88,7 @@ template<typename VectorType, int Size, int PacketAccess> class VectorBlock
using Base::operator-=;
using Base::operator*=;
using Base::operator/=;
/** Dynamic-size constructor
*/
inline VectorBlock(const VectorType& vector, int start, int size)
@ -96,7 +96,7 @@ template<typename VectorType, int Size, int PacketAccess> class VectorBlock
IsColVector ? start : 0, IsColVector ? 0 : start,
IsColVector ? size : 1, IsColVector ? 1 : size)
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorBlock);
}
@ -114,8 +114,6 @@ template<typename VectorType, int Size, int PacketAccess> class VectorBlock
*
* \only_for_vectors
*
* \addexample VectorBlockIntInt \label How to reference a sub-vector (dynamic size)
*
* \param start the first coefficient in the segment
* \param size the number of coefficients in the segment
*
@ -151,8 +149,6 @@ MatrixBase<Derived>::segment(int start, int size) const
*
* \param size the number of coefficients in the block
*
* \addexample BlockInt \label How to reference a sub-vector (fixed-size)
*
* Example: \include MatrixBase_start_int.cpp
* Output: \verbinclude MatrixBase_start_int.out
*
@ -185,8 +181,6 @@ MatrixBase<Derived>::start(int size) const
*
* \param size the number of coefficients in the block
*
* \addexample BlockEnd \label How to reference the end of a vector (fixed-size)
*
* Example: \include MatrixBase_end_int.cpp
* Output: \verbinclude MatrixBase_end_int.out
*
@ -251,8 +245,6 @@ MatrixBase<Derived>::segment(int start) const
*
* The template parameter \a Size is the number of coefficients in the block
*
* \addexample BlockStart \label How to reference the start of a vector (fixed-size)
*
* Example: \include MatrixBase_template_int_start.cpp
* Output: \verbinclude MatrixBase_template_int_start.out
*

12
Eigen/src/Core/products/GeneralMatrixMatrix.h
View File

@ -135,7 +135,11 @@ class GeneralProduct<Lhs, Rhs, GemmProduct>
public:
EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct)
GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {}
GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
{
EIGEN_STATIC_ASSERT((ei_is_same_type<typename Lhs::Scalar, typename Rhs::Scalar>::ret),
YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
}
template<typename Dest> void scaleAndAddTo(Dest& dst, Scalar alpha) const
{
@ -149,9 +153,9 @@ class GeneralProduct<Lhs, Rhs, GemmProduct>
ei_general_matrix_matrix_product<
Scalar,
(_ActualLhsType::Flags&RowMajorBit)?RowMajor:ColMajor, bool(LhsBlasTraits::NeedToConjugate),
(_ActualRhsType::Flags&RowMajorBit)?RowMajor:ColMajor, bool(RhsBlasTraits::NeedToConjugate),
(Dest::Flags&RowMajorBit)?RowMajor:ColMajor>
(_ActualLhsType::Flags&RowMajorBit) ? RowMajor : ColMajor, bool(LhsBlasTraits::NeedToConjugate),
(_ActualRhsType::Flags&RowMajorBit) ? RowMajor : ColMajor, bool(RhsBlasTraits::NeedToConjugate),
(Dest::Flags&RowMajorBit) ? RowMajor : ColMajor>
::run(
this->rows(), this->cols(), lhs.cols(),
(const Scalar*)&(lhs.const_cast_derived().coeffRef(0,0)), lhs.stride(),

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

@ -41,7 +41,7 @@ static EIGEN_DONT_INLINE void ei_product_selfadjoint_vector(
const int PacketSize = sizeof(Packet)/sizeof(Scalar);
enum {
IsRowMajor = StorageOrder==RowMajorBit ? 1 : 0,
IsRowMajor = StorageOrder==RowMajor ? 1 : 0,
IsLower = UpLo == LowerTriangularBit ? 1 : 0,
FirstTriangular = IsRowMajor == IsLower
};
@ -185,14 +185,14 @@ struct SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,0,true>
Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(m_lhs)
* RhsBlasTraits::extractScalarFactor(m_rhs);
ei_assert((&dst.coeff(1))-(&dst.coeff(0))==1 && "not implemented yet");
ei_product_selfadjoint_vector<Scalar, ei_traits<_ActualLhsType>::Flags&RowMajorBit, int(LhsUpLo), bool(LhsBlasTraits::NeedToConjugate), bool(RhsBlasTraits::NeedToConjugate)>
ei_assert(dst.stride()==1 && "not implemented yet");
ei_product_selfadjoint_vector<Scalar, (ei_traits<_ActualLhsType>::Flags&RowMajorBit) ? RowMajor : ColMajor, int(LhsUpLo), bool(LhsBlasTraits::NeedToConjugate), bool(RhsBlasTraits::NeedToConjugate)>
(
lhs.rows(), // size
&lhs.coeff(0,0), lhs.stride(), // lhs info
&rhs.coeff(0), (&rhs.coeff(1))-(&rhs.coeff(0)), // rhs info
&dst.coeffRef(0), // result info
actualAlpha // scale factor
lhs.rows(), // size
&lhs.coeff(0,0), lhs.stride(), // lhs info
&rhs.coeff(0), rhs.stride(), // rhs info
&dst.coeffRef(0), // result info
actualAlpha // scale factor
);
}
};

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

@ -132,7 +132,7 @@ SelfAdjointView<MatrixType,UpLo>& SelfAdjointView<MatrixType,UpLo>
Scalar actualAlpha = alpha * UBlasTraits::extractScalarFactor(u.derived());
enum { IsRowMajor = (ei_traits<MatrixType>::Flags&RowMajorBit)?1:0 };
enum { IsRowMajor = (ei_traits<MatrixType>::Flags&RowMajorBit) ? 1 : 0 };
ei_selfadjoint_product<Scalar,
_ActualUType::Flags&RowMajorBit ? RowMajor : ColMajor,

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

@ -85,7 +85,7 @@ SelfAdjointView<MatrixType,UpLo>& SelfAdjointView<MatrixType,UpLo>
Scalar actualAlpha = alpha * UBlasTraits::extractScalarFactor(u.derived())
* VBlasTraits::extractScalarFactor(v.derived());
enum { IsRowMajor = (ei_traits<MatrixType>::Flags&RowMajorBit)?1:0 };
enum { IsRowMajor = (ei_traits<MatrixType>::Flags&RowMajorBit) ? 1 : 0 };
ei_selfadjoint_rank2_update_selector<Scalar,
typename ei_cleantype<typename ei_conj_expr_if<IsRowMajor ^ UBlasTraits::NeedToConjugate,_ActualUType>::ret>::type,
typename ei_cleantype<typename ei_conj_expr_if<IsRowMajor ^ VBlasTraits::NeedToConjugate,_ActualVType>::ret>::type,

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

@ -145,7 +145,7 @@ struct TriangularProduct<Mode,true,Lhs,false,Rhs,true>
Mode,
LhsBlasTraits::NeedToConjugate,
RhsBlasTraits::NeedToConjugate,
ei_traits<Lhs>::Flags&RowMajorBit>
(int(ei_traits<Lhs>::Flags)&RowMajorBit) ? RowMajor : ColMajor>
::run(lhs,rhs,dst,actualAlpha);
}
};

9
Eigen/src/Core/util/Constants.h
View File

@ -238,6 +238,15 @@ enum {
OnTheRight = 2
};
// options for SVD decomposition
enum {
SkipU = 0x1,
SkipV = 0x2,
AtLeastAsManyRowsAsCols = 0x4,
AtLeastAsManyColsAsRows = 0x8,
Square = AtLeastAsManyRowsAsCols | AtLeastAsManyColsAsRows
};
/* the following could as well be written:
* enum NoChange_t { NoChange };
* but it feels dangerous to disambiguate overloaded functions on enum/integer types.

3
Eigen/src/Core/util/DisableMSVCWarnings.h
View File

@ -1,5 +1,6 @@
#ifdef _MSC_VER
// 4273 - QtAlignedMalloc, inconsistent dll linkage
#pragma warning( push )
#pragma warning( disable : 4181 4244 4127 4211 4717 )
#pragma warning( disable : 4181 4244 4127 4211 4273 4522 4717 )
#endif

2
Eigen/src/Core/util/ForwardDeclarations.h
View File

@ -120,8 +120,10 @@ template<typename MatrixType> class HouseholderQR;
template<typename MatrixType> class ColPivotingHouseholderQR;
template<typename MatrixType> class FullPivotingHouseholderQR;
template<typename MatrixType> class SVD;
template<typename MatrixType, unsigned int Options = 0> class JacobiSVD;
template<typename MatrixType, int UpLo = LowerTriangular> class LLT;
template<typename MatrixType> class LDLT;
template<typename Scalar> class PlanarRotation;
// Geometry module:
template<typename Derived, int _Dim> class RotationBase;

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

@ -246,6 +246,14 @@ using Eigen::ei_cos;
// needed to define it here as escaping characters in CMake add_definition's argument seems very problematic.
#define EIGEN_DOCS_IO_FORMAT IOFormat(3, 0, " ", "\n", "", "")
#ifdef _MSC_VER
#define EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Derived) \
using Base::operator =; \
using Base::operator +=; \
using Base::operator -=; \
using Base::operator *=; \
using Base::operator /=;
#else
#define EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Derived) \
using Base::operator =; \
using Base::operator +=; \
@ -256,6 +264,7 @@ EIGEN_STRONG_INLINE Derived& operator=(const Derived& other) \
{ \
return Base::operator=(other); \
}
#endif
#define _EIGEN_GENERIC_PUBLIC_INTERFACE(Derived, BaseClass) \
typedef BaseClass Base; \
@ -278,6 +287,9 @@ enum { RowsAtCompileTime = Eigen::ei_traits<Derived>::RowsAtCompileTime, \
_EIGEN_GENERIC_PUBLIC_INTERFACE(Derived, Eigen::MatrixBase<Derived>)
#define EIGEN_ENUM_MIN(a,b) (((int)a <= (int)b) ? (int)a : (int)b)
#define EIGEN_SIZE_MIN(a,b) (((int)a == 1 || (int)b == 1) ? 1 \
: ((int)a == Dynamic || (int)b == Dynamic) ? Dynamic \
: ((int)a <= (int)b) ? (int)a : (int)b)
#define EIGEN_ENUM_MAX(a,b) (((int)a >= (int)b) ? (int)a : (int)b)
#define EIGEN_LOGICAL_XOR(a,b) (((a) || (b)) && !((a) && (b)))

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

@ -62,6 +62,7 @@
THIS_METHOD_IS_ONLY_FOR_MATRICES_OF_A_SPECIFIC_SIZE,
YOU_MADE_A_PROGRAMMING_MISTAKE,
YOU_CALLED_A_FIXED_SIZE_METHOD_ON_A_DYNAMIC_SIZE_MATRIX_OR_VECTOR,
YOU_CALLED_A_DYNAMIC_SIZE_METHOD_ON_A_FIXED_SIZE_MATRIX_OR_VECTOR,
UNALIGNED_LOAD_AND_STORE_OPERATIONS_UNIMPLEMENTED_ON_ALTIVEC,
NUMERIC_TYPE_MUST_BE_FLOATING_POINT,
NUMERIC_TYPE_MUST_BE_REAL,
@ -114,6 +115,11 @@
EIGEN_STATIC_ASSERT(TYPE::SizeAtCompileTime!=Eigen::Dynamic, \
YOU_CALLED_A_FIXED_SIZE_METHOD_ON_A_DYNAMIC_SIZE_MATRIX_OR_VECTOR)
// static assertion failing if the type \a TYPE is not dynamic-size
#define EIGEN_STATIC_ASSERT_DYNAMIC_SIZE(TYPE) \
EIGEN_STATIC_ASSERT(TYPE::SizeAtCompileTime==Eigen::Dynamic, \
YOU_CALLED_A_DYNAMIC_SIZE_METHOD_ON_A_FIXED_SIZE_MATRIX_OR_VECTOR)
// static assertion failing if the type \a TYPE is not a vector type of the given size
#define EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(TYPE, SIZE) \
EIGEN_STATIC_ASSERT(TYPE::IsVectorAtCompileTime && TYPE::SizeAtCompileTime==SIZE, \

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

@ -217,7 +217,7 @@ template<unsigned int Flags> struct ei_are_flags_consistent
* overloads for complex types */
template<typename Derived,typename Scalar,typename OtherScalar,
bool EnableIt = !ei_is_same_type<Scalar,OtherScalar>::ret >
struct ei_special_scalar_op_base
struct ei_special_scalar_op_base : public AnyMatrixBase<Derived>
{
// dummy operator* so that the
// "using ei_special_scalar_op_base::operator*" compiles
@ -225,7 +225,7 @@ struct ei_special_scalar_op_base
};
template<typename Derived,typename Scalar,typename OtherScalar>
struct ei_special_scalar_op_base<Derived,Scalar,OtherScalar,true>
struct ei_special_scalar_op_base<Derived,Scalar,OtherScalar,true> : public AnyMatrixBase<Derived>
{
const CwiseUnaryOp<ei_scalar_multiple2_op<Scalar,OtherScalar>, Derived>
operator*(const OtherScalar& scalar) const
@ -233,6 +233,10 @@ struct ei_special_scalar_op_base<Derived,Scalar,OtherScalar,true>
return CwiseUnaryOp<ei_scalar_multiple2_op<Scalar,OtherScalar>, Derived>
(*static_cast<const Derived*>(this), ei_scalar_multiple2_op<Scalar,OtherScalar>(scalar));
}
inline friend const CwiseUnaryOp<ei_scalar_multiple2_op<Scalar,OtherScalar>, Derived>
operator*(const OtherScalar& scalar, const Derived& matrix)
{ return matrix*scalar; }
};
/** \internal Gives the type of a sub-matrix or sub-vector of a matrix of type \a ExpressionType and size \a Size

6
Eigen/src/Eigenvalues/CMakeLists.txt Normal file
View File

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

148
Eigen/src/Eigenvalues/ComplexEigenSolver.h Normal file
View File

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

237
Eigen/src/Eigenvalues/ComplexSchur.h Normal file
View File

@ -0,0 +1,237 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Claire Maurice
// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifndef EIGEN_COMPLEX_SCHUR_H
#define EIGEN_COMPLEX_SCHUR_H
/** \eigenvalues_module \ingroup Eigenvalues_Module
* \nonstableyet
*
* \class ComplexShur
*
* \brief Performs a complex Shur decomposition of a real or complex square matrix
*
*/
template<typename _MatrixType> class ComplexSchur
{
public:
typedef _MatrixType MatrixType;
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef std::complex<RealScalar> Complex;
typedef Matrix<Complex, MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime> ComplexMatrixType;
/**
* \brief Default Constructor.
*
* The default constructor is useful in cases in which the user intends to
* perform decompositions via ComplexSchur::compute(const MatrixType&).
*/
ComplexSchur() : m_matT(), m_matU(), m_isInitialized(false)
{}
ComplexSchur(const MatrixType& matrix)
: m_matT(matrix.rows(),matrix.cols()),
m_matU(matrix.rows(),matrix.cols()),
m_isInitialized(false)
{
compute(matrix);
}
ComplexMatrixType matrixU() const
{
ei_assert(m_isInitialized && "ComplexSchur is not initialized.");
return m_matU;
}
ComplexMatrixType matrixT() const
{
ei_assert(m_isInitialized && "ComplexShur is not initialized.");
return m_matT;
}
void compute(const MatrixType& matrix);
protected:
ComplexMatrixType m_matT, m_matU;
bool m_isInitialized;
};
/** Computes the principal value of the square root of the complex \a z. */
template<typename RealScalar>
std::complex<RealScalar> ei_sqrt(const std::complex<RealScalar> &z)
{
RealScalar t, tre, tim;
t = ei_abs(z);
if (ei_abs(ei_real(z)) <= ei_abs(ei_imag(z)))
{
// No cancellation in these formulas
tre = ei_sqrt(0.5*(t + ei_real(z)));
tim = ei_sqrt(0.5*(t - ei_real(z)));
}
else
{
// Stable computation of the above formulas
if (z.real() > 0)
{
tre = t + z.real();
tim = ei_abs(ei_imag(z))*ei_sqrt(0.5/tre);
tre = ei_sqrt(0.5*tre);
}
else
{
tim = t - z.real();
tre = ei_abs(ei_imag(z))*ei_sqrt(0.5/tim);
tim = ei_sqrt(0.5*tim);
}
}
if(z.imag() < 0)
tim = -tim;
return (std::complex<RealScalar>(tre,tim));
}
template<typename MatrixType>
void ComplexSchur<MatrixType>::compute(const MatrixType& matrix)
{
// this code is inspired from Jampack
assert(matrix.cols() == matrix.rows());
int n = matrix.cols();
// Reduce to Hessenberg form
HessenbergDecomposition<MatrixType> hess(matrix);
m_matT = hess.matrixH();
m_matU = hess.matrixQ();
int iu = m_matT.cols() - 1;
int il;
RealScalar d,sd,sf;
Complex c,b,disc,r1,r2,kappa;
RealScalar eps = epsilon<RealScalar>();
int iter = 0;
while(true)
{
//locate the range in which to iterate
while(iu > 0)
{
d = ei_norm1(m_matT.coeffRef(iu,iu)) + ei_norm1(m_matT.coeffRef(iu-1,iu-1));
sd = ei_norm1(m_matT.coeffRef(iu,iu-1));
if(sd >= eps * d) break; // FIXME : precision criterion ??
m_matT.coeffRef(iu,iu-1) = Complex(0);
iter = 0;
--iu;
}
if(iu==0) break;
iter++;
if(iter >= 30)
{
// FIXME : what to do when iter==MAXITER ??
std::cerr << "MAXITER" << std::endl;
return;
}
il = iu-1;
while( il > 0 )
{
// check if the current 2x2 block on the diagonal is upper triangular
d = ei_norm1(m_matT.coeffRef(il,il)) + ei_norm1(m_matT.coeffRef(il-1,il-1));
sd = ei_norm1(m_matT.coeffRef(il,il-1));
if(sd < eps * d) break; // FIXME : precision criterion ??
--il;
}
if( il != 0 ) m_matT.coeffRef(il,il-1) = Complex(0);
// compute the shift (the normalization by sf is to avoid under/overflow)
Matrix<Scalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1);
sf = t.cwise().abs().sum();
t /= sf;
c = t.determinant();
b = t.diagonal().sum();
disc = ei_sqrt(b*b - RealScalar(4)*c);
r1 = (b+disc)/RealScalar(2);
r2 = (b-disc)/RealScalar(2);
if(ei_norm1(r1) > ei_norm1(r2))
r2 = c/r1;
else
r1 = c/r2;
if(ei_norm1(r1-t.coeff(1,1)) < ei_norm1(r2-t.coeff(1,1)))
kappa = sf * r1;
else
kappa = sf * r2;
// perform the QR step using Givens rotations
PlanarRotation<Complex> rot;
rot.makeGivens(m_matT.coeff(il,il) - kappa, m_matT.coeff(il+1,il));
for(int i=il ; i<iu ; i++)
{
m_matT.block(0,i,n,n-i).applyOnTheLeft(i, i+1, rot.adjoint());
m_matT.block(0,0,std::min(i+2,iu)+1,n).applyOnTheRight(i, i+1, rot);
m_matU.applyOnTheRight(i, i+1, rot);
if(i != iu-1)
{
int i1 = i+1;
int i2 = i+2;
rot.makeGivens(m_matT.coeffRef(i1,i), m_matT.coeffRef(i2,i), &m_matT.coeffRef(i1,i));
m_matT.coeffRef(i2,i) = Complex(0);
}
}
}
// FIXME : is it necessary ?
/*
for(int i=0 ; i<n ; i++)
for(int j=0 ; j<n ; j++)
{
if(ei_abs(ei_real(m_matT.coeff(i,j))) < eps)
ei_real_ref(m_matT.coeffRef(i,j)) = 0;
if(ei_imag(ei_abs(m_matT.coeff(i,j))) < eps)
ei_imag_ref(m_matT.coeffRef(i,j)) = 0;
}
*/
m_isInitialized = true;
}
#endif // EIGEN_COMPLEX_SCHUR_H

28
Eigen/src/QR/EigenSolver.h → Eigen/src/Eigenvalues/EigenSolver.h
View File

@ -25,7 +25,7 @@
#ifndef EIGEN_EIGENSOLVER_H
#define EIGEN_EIGENSOLVER_H
/** \ingroup QR_Module
/** \eigenvalues_module \ingroup Eigenvalues_Module
* \nonstableyet
*
* \class EigenSolver
@ -53,7 +53,7 @@ template<typename _MatrixType> class EigenSolver
typedef Matrix<RealScalar, MatrixType::ColsAtCompileTime, 1> RealVectorType;
typedef Matrix<RealScalar, Dynamic, 1> RealVectorTypeX;
/**
/**
* \brief Default Constructor.
*
* The default constructor is useful in cases in which the user intends to
@ -103,19 +103,19 @@ template<typename _MatrixType> class EigenSolver
*
* \sa pseudoEigenvalueMatrix()
*/
const MatrixType& pseudoEigenvectors() const
{
const MatrixType& pseudoEigenvectors() const
{
ei_assert(m_isInitialized && "EigenSolver is not initialized.");
return m_eivec;
return m_eivec;
}
MatrixType pseudoEigenvalueMatrix() const;
/** \returns the eigenvalues as a column vector */
EigenvalueType eigenvalues() const
{
EigenvalueType eigenvalues() const
{
ei_assert(m_isInitialized && "EigenSolver is not initialized.");
return m_eivalues;
return m_eivalues;
}
EigenSolver& compute(const MatrixType& matrix);
@ -244,11 +244,11 @@ void EigenSolver<MatrixType>::orthes(MatrixType& matH, RealVectorType& ort)
// Apply Householder similarity transformation
// H = (I-u*u'/h)*H*(I-u*u')/h)
int bSize = high-m+1;
matH.block(m, m, bSize, n-m) -= ((ort.segment(m, bSize)/h)
* (ort.segment(m, bSize).transpose() * matH.block(m, m, bSize, n-m))).lazy();
matH.block(m, m, bSize, n-m).noalias() -= ((ort.segment(m, bSize)/h)
* (ort.segment(m, bSize).transpose() * matH.block(m, m, bSize, n-m)));
matH.block(0, m, high+1, bSize) -= ((matH.block(0, m, high+1, bSize) * ort.segment(m, bSize))
* (ort.segment(m, bSize)/h).transpose()).lazy();
matH.block(0, m, high+1, bSize).noalias() -= ((matH.block(0, m, high+1, bSize) * ort.segment(m, bSize))
* (ort.segment(m, bSize)/h).transpose());
ort.coeffRef(m) = scale*ort.coeff(m);
matH.coeffRef(m,m-1) = scale*g;
@ -265,8 +265,8 @@ void EigenSolver<MatrixType>::orthes(MatrixType& matH, RealVectorType& ort)
ort.segment(m+1, high-m) = matH.col(m-1).segment(m+1, high-m);
int bSize = high-m+1;
m_eivec.block(m, m, bSize, bSize) += ( (ort.segment(m, bSize) / (matH.coeff(m,m-1) * ort.coeff(m)))
* (ort.segment(m, bSize).transpose() * m_eivec.block(m, m, bSize, bSize)) ).lazy();
m_eivec.block(m, m, bSize, bSize).noalias() += ( (ort.segment(m, bSize) / (matH.coeff(m,m-1) * ort.coeff(m)))
* (ort.segment(m, bSize).transpose() * m_eivec.block(m, m, bSize, bSize)) );
}
}
}

4
Eigen/src/QR/HessenbergDecomposition.h → Eigen/src/Eigenvalues/HessenbergDecomposition.h
View File

@ -25,7 +25,7 @@
#ifndef EIGEN_HESSENBERGDECOMPOSITION_H
#define EIGEN_HESSENBERGDECOMPOSITION_H
/** \ingroup QR_Module
/** \eigenvalues_module \ingroup Eigenvalues_Module
* \nonstableyet
*
* \class HessenbergDecomposition
@ -156,7 +156,7 @@ void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVector
// Apply similarity transformation to remaining columns,
// i.e., compute A = H A H'
// A = H A
matA.corner(BottomRight, remainingSize, remainingSize)
.applyHouseholderOnTheLeft(matA.col(i).end(remainingSize-1), h, &temp.coeffRef(0));

53
Eigen/src/QR/SelfAdjointEigenSolver.h → Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
View File

@ -25,7 +25,7 @@
#ifndef EIGEN_SELFADJOINTEIGENSOLVER_H
#define EIGEN_SELFADJOINTEIGENSOLVER_H
/** \qr_module \ingroup QR_Module
/** \eigenvalues_module \ingroup Eigenvalues_Module
* \nonstableyet
*
* \class SelfAdjointEigenSolver
@ -135,31 +135,9 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
#ifndef EIGEN_HIDE_HEAVY_CODE
// from Golub's "Matrix Computations", algorithm 5.1.3
template<typename Scalar>
static void ei_givens_rotation(Scalar a, Scalar b, Scalar& c, Scalar& s)
{
if (b==0)
{
c = 1; s = 0;
}
else if (ei_abs(b)>ei_abs(a))
{
Scalar t = -a/b;
s = Scalar(1)/ei_sqrt(1+t*t);
c = s * t;
}
else
{
Scalar t = -b/a;
c = Scalar(1)/ei_sqrt(1+t*t);
s = c * t;
}
}
/** \internal
*
* \qr_module
* \eigenvalues_module \ingroup Eigenvalues_Module
*
* Performs a QR step on a tridiagonal symmetric matrix represented as a
* pair of two vectors \a diag and \a subdiag.
@ -288,7 +266,7 @@ compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors
#endif // EIGEN_HIDE_HEAVY_CODE
/** \qr_module
/** \eigenvalues_module
*
* \returns a vector listing the eigenvalues of this matrix.
*/
@ -329,7 +307,7 @@ template<typename Derived> struct ei_operatorNorm_selector<Derived, false>
}
};
/** \qr_module
/** \eigenvalues_module
*
* \returns the matrix norm of this matrix.
*/
@ -353,34 +331,33 @@ static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int st
for (int k = start; k < end; ++k)
{
RealScalar c, s;
ei_givens_rotation(x, z, c, s);
PlanarRotation<RealScalar> rot;
rot.makeGivens(x, z);
// do T = G' T G
RealScalar sdk = s * diag[k] + c * subdiag[k];
RealScalar dkp1 = s * subdiag[k] + c * diag[k+1];
RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k];
RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k+1];
diag[k] = c * (c * diag[k] - s * subdiag[k]) - s * (c * subdiag[k] - s * diag[k+1]);
diag[k+1] = s * sdk + c * dkp1;
subdiag[k] = c * sdk - s * dkp1;
diag[k] = rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k+1]);
diag[k+1] = rot.s() * sdk + rot.c() * dkp1;
subdiag[k] = rot.c() * sdk - rot.s() * dkp1;
if (k > start)
subdiag[k - 1] = c * subdiag[k-1] - s * z;
subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z;
x = subdiag[k];
if (k < end - 1)
{
z = -s * subdiag[k+1];
subdiag[k + 1] = c * subdiag[k+1];
z = -rot.s() * subdiag[k+1];
subdiag[k + 1] = rot.c() * subdiag[k+1];
}
// apply the givens rotation to the unit matrix Q = Q * G
// G only modifies the two columns k and k+1
if (matrixQ)
{
Map<Matrix<Scalar,Dynamic,Dynamic> > q(matrixQ,n,n);
q.applyJacobiOnTheRight(k,k+1,c,s);
q.applyOnTheRight(k,k+1,rot);
}
}
}

2
Eigen/src/QR/Tridiagonalization.h → Eigen/src/Eigenvalues/Tridiagonalization.h
View File

@ -25,7 +25,7 @@
#ifndef EIGEN_TRIDIAGONALIZATION_H
#define EIGEN_TRIDIAGONALIZATION_H
/** \ingroup QR_Module
/** \eigenvalues_module \ingroup Eigenvalues_Module
* \nonstableyet
*
* \class Tridiagonalization

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

@ -39,8 +39,6 @@
* \li \c AngleAxisf for \c float
* \li \c AngleAxisd for \c double
*
* \addexample AngleAxisForEuler \label How to define a rotation from Euler-angles
*
* Combined with MatrixBase::Unit{X,Y,Z}, AngleAxis can be used to easily
* mimic Euler-angles. Here is an example:
* \include AngleAxis_mimic_euler.cpp
@ -85,7 +83,7 @@ public:
AngleAxis() {}
/** Constructs and initialize the angle-axis rotation from an \a angle in radian
* and an \a axis which \b must \b be \b normalized.
*
*
* \warning If the \a axis vector is not normalized, then the angle-axis object
* represents an invalid rotation. */
template<typename Derived>

10
Eigen/src/Geometry/Umeyama.h
View File

@ -144,7 +144,7 @@ umeyama(const MatrixBase<Derived>& src, const MatrixBase<OtherDerived>& dst, boo
// const Scalar dst_var = dst_demean.rowwise().squaredNorm().sum() * one_over_n;
// Eq. (38)
const MatrixType sigma = (one_over_n * dst_demean * src_demean.transpose()).lazy();
const MatrixType sigma = one_over_n * dst_demean * src_demean.transpose();
SVD<MatrixType> svd(sigma);
@ -160,14 +160,14 @@ umeyama(const MatrixBase<Derived>& src, const MatrixBase<OtherDerived>& dst, boo
int rank = 0; for (int i=0; i<m; ++i) if (!ei_isMuchSmallerThan(d.coeff(i),d.coeff(0))) ++rank;
if (rank == m-1) {
if ( svd.matrixU().determinant() * svd.matrixV().determinant() > 0 ) {
Rt.block(0,0,m,m) = (svd.matrixU()*svd.matrixV().transpose()).lazy();
Rt.block(0,0,m,m).noalias() = svd.matrixU()*svd.matrixV().transpose();
} else {
const Scalar s = S(m-1); S(m-1) = -1;
Rt.block(0,0,m,m) = (svd.matrixU() * S.asDiagonal() * svd.matrixV().transpose()).lazy();
Rt.block(0,0,m,m).noalias() = svd.matrixU() * S.asDiagonal() * svd.matrixV().transpose();
S(m-1) = s;
}
} else {
Rt.block(0,0,m,m) = (svd.matrixU() * S.asDiagonal() * svd.matrixV().transpose()).lazy();
Rt.block(0,0,m,m).noalias() = svd.matrixU() * S.asDiagonal() * svd.matrixV().transpose();
}
// Eq. (42)
@ -177,7 +177,7 @@ umeyama(const MatrixBase<Derived>& src, const MatrixBase<OtherDerived>& dst, boo
// Note that we first assign dst_mean to the destination so that there no need
// for a temporary.
Rt.col(m).start(m) = dst_mean;
Rt.col(m).start(m) -= (c*Rt.corner(TopLeft,m,m)*src_mean).lazy();
Rt.col(m).start(m).noalias() -= c*Rt.corner(TopLeft,m,m)*src_mean;
if (with_scaling) Rt.block(0,0,m,m) *= c;

364
Eigen/src/Jacobi/Jacobi.h
View File

@ -26,108 +26,282 @@
#ifndef EIGEN_JACOBI_H
#define EIGEN_JACOBI_H
/** Applies the counter clock wise 2D rotation of angle \c theta given by its
* cosine \a c and sine \a s to the set of 2D vectors of cordinates \a x and \a y:
* \f$ x = c x - s' y \f$
* \f$ y = s x + c y \f$
/** \ingroup Jacobi_Module
* \jacobi_module
* \class PlanarRotation
* \brief Represents a rotation in the plane from a cosine-sine pair.
*
* \sa MatrixBase::applyJacobiOnTheLeft(), MatrixBase::applyJacobiOnTheRight()
* This class represents a Jacobi or Givens rotation.
* This is a 2D rotation in the plane \c J of angle \f$ \theta \f$ defined by
* its cosine \c c and sine \c s as follow:
* \f$ J = \left ( \begin{array}{cc} c & \overline s \\ -s & \overline c \end{array} \right ) \f$
*
* You can apply the respective counter-clockwise rotation to a column vector \c v by
* applying its adjoint on the left: \f$ v = J^* v \f$ that translates to the following Eigen code:
* \code
* v.applyOnTheLeft(J.adjoint());
* \endcode
*
* \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight()
*/
template<typename VectorX, typename VectorY>
void ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, typename VectorX::Scalar c, typename VectorY::Scalar s);
/** Applies a rotation in the plane defined by \a c, \a s to the rows \a p and \a q of \c *this.
* More precisely, it computes B = J' * B, with J = [c s ; -s' c] and B = [ *this.row(p) ; *this.row(q) ]
* \sa MatrixBase::applyJacobiOnTheRight(), ei_apply_rotation_in_the_plane()
*/
template<typename Derived>
inline void MatrixBase<Derived>::applyJacobiOnTheLeft(int p, int q, Scalar c, Scalar s)
template<typename Scalar> class PlanarRotation
{
RowXpr x(row(p));
RowXpr y(row(q));
ei_apply_rotation_in_the_plane(x, y, ei_conj(c), ei_conj(s));
}
public:
typedef typename NumTraits<Scalar>::Real RealScalar;
/** Applies a rotation in the plane defined by \a c, \a s to the columns \a p and \a q of \c *this.
* More precisely, it computes B = B * J, with J = [c s ; -s' c] and B = [ *this.col(p) ; *this.col(q) ]
* \sa MatrixBase::applyJacobiOnTheLeft(), ei_apply_rotation_in_the_plane()
*/
template<typename Derived>
inline void MatrixBase<Derived>::applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s)
{
ColXpr x(col(p));
ColXpr y(col(q));
ei_apply_rotation_in_the_plane(x, y, c, s);
}
/** Default constructor without any initialization. */
PlanarRotation() {}
/** Computes the cosine-sine pair (\a c, \a s) such that its associated
* rotation \f$ J = ( \begin{array}{cc} c & s \\ -s' c \end{array} )\f$
* applied to both the right and left of the 2x2 matrix
* \f$ B = ( \begin{array}{cc} x & y \\ * & z \end{array} )\f$ yields
* a diagonal matrix A: \f$ A = J' B J \f$
/** Construct a planar rotation from a cosine-sine pair (\a c, \c s). */
PlanarRotation(const Scalar& c, const Scalar& s) : m_c(c), m_s(s) {}
Scalar& c() { return m_c; }
Scalar c() const { return m_c; }
Scalar& s() { return m_s; }
Scalar s() const { return m_s; }
/** Concatenates two planar rotation */
PlanarRotation operator*(const PlanarRotation& other)
{
return PlanarRotation(m_c * other.m_c - ei_conj(m_s) * other.m_s,
ei_conj(m_c * ei_conj(other.m_s) + ei_conj(m_s) * ei_conj(other.m_c)));
}
/** Returns the transposed transformation */
PlanarRotation transpose() const { return PlanarRotation(m_c, -ei_conj(m_s)); }
/** Returns the adjoint transformation */
PlanarRotation adjoint() const { return PlanarRotation(ei_conj(m_c), -m_s); }
template<typename Derived>
bool makeJacobi(const MatrixBase<Derived>&, int p, int q);
bool makeJacobi(RealScalar x, Scalar y, RealScalar z);
void makeGivens(const Scalar& p, const Scalar& q, Scalar* z=0);
protected:
void makeGivens(const Scalar& p, const Scalar& q, Scalar* z, ei_meta_true);
void makeGivens(const Scalar& p, const Scalar& q, Scalar* z, ei_meta_false);
Scalar m_c, m_s;
};
/** Makes \c *this as a Jacobi rotation \a J such that applying \a J on both the right and left sides of the selfadjoint 2x2 matrix
* \f$ B = \left ( \begin{array}{cc} x & y \\ \overline y & z \end{array} \right )\f$ yields a diagonal matrix \f$ A = J^* B J \f$
*
* \sa MatrixBase::makeJacobi(const MatrixBase<Derived>&, int, int), MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight()
*/
template<typename Scalar>
bool ei_makeJacobi(Scalar x, Scalar y, Scalar z, Scalar *c, Scalar *s)
bool PlanarRotation<Scalar>::makeJacobi(RealScalar x, Scalar y, RealScalar z)
{
if(y == 0)
typedef typename NumTraits<Scalar>::Real RealScalar;
if(y == Scalar(0))
{
*c = Scalar(1);
*s = Scalar(0);
m_c = Scalar(1);
m_s = Scalar(0);
return false;
}
else
{
Scalar tau = (z - x) / (2 * y);
Scalar w = ei_sqrt(1 + ei_abs2(tau));
Scalar t;
RealScalar tau = (x-z)/(RealScalar(2)*ei_abs(y));
RealScalar w = ei_sqrt(ei_abs2(tau) + 1);
RealScalar t;
if(tau>0)
t = Scalar(1) / (tau + w);
{
t = RealScalar(1) / (tau + w);
}
else
t = Scalar(1) / (tau - w);
*c = Scalar(1) / ei_sqrt(1 + ei_abs2(t));
*s = *c * t;
{
t = RealScalar(1) / (tau - w);
}
RealScalar sign_t = t > 0 ? 1 : -1;
RealScalar n = RealScalar(1) / ei_sqrt(ei_abs2(t)+1);
m_s = - sign_t * (ei_conj(y) / ei_abs(y)) * ei_abs(t) * n;
m_c = n;
return true;
}
}
template<typename Derived>
inline bool MatrixBase<Derived>::makeJacobi(int p, int q, Scalar *c, Scalar *s) const
{
return ei_makeJacobi(coeff(p,p), coeff(p,q), coeff(q,q), c, s);
}
template<typename Derived>
inline bool MatrixBase<Derived>::makeJacobiForAtA(int p, int q, Scalar *c, Scalar *s) const
{
return ei_makeJacobi(ei_abs2(coeff(p,p)) + ei_abs2(coeff(q,p)),
ei_conj(coeff(p,p))*coeff(p,q) + ei_conj(coeff(q,p))*coeff(q,q),
ei_abs2(coeff(p,q)) + ei_abs2(coeff(q,q)),
c,s);
}
template<typename Derived>
inline bool MatrixBase<Derived>::makeJacobiForAAt(int p, int q, Scalar *c, Scalar *s) const
{
return ei_makeJacobi(ei_abs2(coeff(p,p)) + ei_abs2(coeff(p,q)),
ei_conj(coeff(q,p))*coeff(p,p) + ei_conj(coeff(q,q))*coeff(p,q),
ei_abs2(coeff(q,p)) + ei_abs2(coeff(q,q)),
c,s);
}
/** Makes \c *this as a Jacobi rotation \c J such that applying \a J on both the right and left sides of the 2x2 selfadjoint matrix
* \f$ B = \left ( \begin{array}{cc} \text{this}_{pp} & \text{this}_{pq} \\ (\text{this}_{pq})^* & \text{this}_{qq} \end{array} \right )\f$ yields
* a diagonal matrix \f$ A = J^* B J \f$
*
* Example: \include Jacobi_makeJacobi.cpp
* Output: \verbinclude Jacobi_makeJacobi.out
*
* \sa PlanarRotation::makeJacobi(RealScalar, Scalar, RealScalar), MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight()
*/
template<typename Scalar>
inline void ei_normalizeJacobi(Scalar *c, Scalar *s, const Scalar& x, const Scalar& y)
template<typename Derived>
inline bool PlanarRotation<Scalar>::makeJacobi(const MatrixBase<Derived>& m, int p, int q)
{
Scalar a = x * *c - y * *s;
Scalar b = x * *s + y * *c;
if(ei_abs(b)>ei_abs(a)) {
Scalar x = *c;
*c = -*s;
*s = x;
return makeJacobi(ei_real(m.coeff(p,p)), m.coeff(p,q), ei_real(m.coeff(q,q)));
}
/** Makes \c *this as a Givens rotation \c G such that applying \f$ G^* \f$ to the left of the vector
* \f$ V = \left ( \begin{array}{c} p \\ q \end{array} \right )\f$ yields:
* \f$ G^* V = \left ( \begin{array}{c} r \\ 0 \end{array} \right )\f$.
*
* The value of \a z is returned if \a z is not null (the default is null).
* Also note that G is built such that the cosine is always real.
*
* Example: \include Jacobi_makeGivens.cpp
* Output: \verbinclude Jacobi_makeGivens.out
*
* This function implements the continuous Givens rotation generation algorithm
* found in Anderson (2000), Discontinuous Plane Rotations and the Symmetric Eigenvalue Problem.
* LAPACK Working Note 150, University of Tennessee, UT-CS-00-454, December 4, 2000.
*
* \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight()
*/
template<typename Scalar>
void PlanarRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar* z)
{
makeGivens(p, q, z, typename ei_meta_if<NumTraits<Scalar>::IsComplex, ei_meta_true, ei_meta_false>::ret());
}
// specialization for complexes
template<typename Scalar>
void PlanarRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar* r, ei_meta_true)
{
if(q==Scalar(0))
{
m_c = ei_real(p)<0 ? Scalar(-1) : Scalar(1);
m_s = 0;
if(r) *r = m_c * p;
}
else if(p==Scalar(0))
{
m_c = 0;
m_s = -q/ei_abs(q);
if(r) *r = ei_abs(q);
}
else
{
RealScalar p1 = ei_norm1(p);
RealScalar q1 = ei_norm1(q);
if(p1>=q1)
{
Scalar ps = p / p1;
RealScalar p2 = ei_abs2(ps);
Scalar qs = q / p1;
RealScalar q2 = ei_abs2(qs);
RealScalar u = ei_sqrt(RealScalar(1) + q2/p2);
if(ei_real(p)<RealScalar(0))
u = -u;
m_c = Scalar(1)/u;
m_s = -qs*ei_conj(ps)*(m_c/p2);
if(r) *r = p * u;
}
else
{
Scalar ps = p / q1;
RealScalar p2 = ei_abs2(ps);
Scalar qs = q / q1;
RealScalar q2 = ei_abs2(qs);
RealScalar u = q1 * ei_sqrt(p2 + q2);
if(ei_real(p)<RealScalar(0))
u = -u;
p1 = ei_abs(p);
ps = p/p1;
m_c = p1/u;
m_s = -ei_conj(ps) * (q/u);
if(r) *r = ps * u;
}
}
}
template<typename VectorX, typename VectorY>
void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, typename VectorX::Scalar c, typename VectorY::Scalar s)
// specialization for reals
template<typename Scalar>
void PlanarRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar* r, ei_meta_false)
{
if(q==0)
{
m_c = p<Scalar(0) ? Scalar(-1) : Scalar(1);
m_s = 0;
if(r) *r = ei_abs(p);
}
else if(p==0)
{
m_c = 0;
m_s = q<Scalar(0) ? Scalar(1) : Scalar(-1);
if(r) *r = ei_abs(q);
}
else if(ei_abs(p) > ei_abs(q))
{
Scalar t = q/p;
Scalar u = ei_sqrt(Scalar(1) + ei_abs2(t));
if(p<Scalar(0))
u = -u;
m_c = Scalar(1)/u;
m_s = -t * m_c;
if(r) *r = p * u;
}
else
{
Scalar t = p/q;
Scalar u = ei_sqrt(Scalar(1) + ei_abs2(t));
if(q<Scalar(0))
u = -u;
m_s = -Scalar(1)/u;
m_c = -t * m_s;
if(r) *r = q * u;
}
}
/****************************************************************************************
* Implementation of MatrixBase methods
****************************************************************************************/
/** \jacobi_module
* Applies the clock wise 2D rotation \a j to the set of 2D vectors of cordinates \a x and \a y:
* \f$ \left ( \begin{array}{cc} x \\ y \end{array} \right ) = J \left ( \begin{array}{cc} x \\ y \end{array} \right ) \f$
*
* \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight()
*/
template<typename VectorX, typename VectorY, typename OtherScalar>
void ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, const PlanarRotation<OtherScalar>& j);
/** \jacobi_module
* Applies the rotation in the plane \a j to the rows \a p and \a q of \c *this, i.e., it computes B = J * B,
* with \f$ B = \left ( \begin{array}{cc} \text{*this.row}(p) \\ \text{*this.row}(q) \end{array} \right ) \f$.
*
* \sa class PlanarRotation, MatrixBase::applyOnTheRight(), ei_apply_rotation_in_the_plane()
*/
template<typename Derived>
template<typename OtherScalar>
inline void MatrixBase<Derived>::applyOnTheLeft(int p, int q, const PlanarRotation<OtherScalar>& j)
{
RowXpr x(row(p));
RowXpr y(row(q));
ei_apply_rotation_in_the_plane(x, y, j);
}
/** \ingroup Jacobi_Module
* Applies the rotation in the plane \a j to the columns \a p and \a q of \c *this, i.e., it computes B = B * J
* with \f$ B = \left ( \begin{array}{cc} \text{*this.col}(p) & \text{*this.col}(q) \end{array} \right ) \f$.
*
* \sa class PlanarRotation, MatrixBase::applyOnTheLeft(), ei_apply_rotation_in_the_plane()
*/
template<typename Derived>
template<typename OtherScalar>
inline void MatrixBase<Derived>::applyOnTheRight(int p, int q, const PlanarRotation<OtherScalar>& j)
{
ColXpr x(col(p));
ColXpr y(col(q));
ei_apply_rotation_in_the_plane(x, y, j.transpose());
}
template<typename VectorX, typename VectorY, typename OtherScalar>
void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, const PlanarRotation<OtherScalar>& j)
{
typedef typename VectorX::Scalar Scalar;
ei_assert(_x.size() == _y.size());
@ -138,7 +312,7 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY&
Scalar* EIGEN_RESTRICT x = &_x.coeffRef(0);
Scalar* EIGEN_RESTRICT y = &_y.coeffRef(0);
if (incrx==1 && incry==1)
if((VectorX::Flags & VectorY::Flags & PacketAccessBit) && incrx==1 && incry==1)
{
// both vectors are sequentially stored in memory => vectorization
typedef typename ei_packet_traits<Scalar>::type Packet;
@ -147,16 +321,16 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY&
int alignedStart = ei_alignmentOffset(y, size);
int alignedEnd = alignedStart + ((size-alignedStart)/PacketSize)*PacketSize;
const Packet pc = ei_pset1(c);
const Packet ps = ei_pset1(s);
const Packet pc = ei_pset1(Scalar(j.c()));
const Packet ps = ei_pset1(Scalar(j.s()));
ei_conj_helper<NumTraits<Scalar>::IsComplex,false> cj;
for(int i=0; i<alignedStart; ++i)
{
Scalar xi = x[i];
Scalar yi = y[i];
x[i] = c * xi - ei_conj(s) * yi;
y[i] = s * xi + c * yi;
x[i] = j.c() * xi + ei_conj(j.s()) * yi;
y[i] = -j.s() * xi + ei_conj(j.c()) * yi;
}
Scalar* px = x + alignedStart;
@ -168,8 +342,8 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY&
{
Packet xi = ei_pload(px);
Packet yi = ei_pload(py);
ei_pstore(px, ei_psub(ei_pmul(pc,xi),cj.pmul(ps,yi)));
ei_pstore(py, ei_padd(ei_pmul(ps,xi),ei_pmul(pc,yi)));
ei_pstore(px, ei_padd(ei_pmul(pc,xi),cj.pmul(ps,yi)));
ei_pstore(py, ei_psub(ei_pmul(pc,yi),ei_pmul(ps,xi)));
px += PacketSize;
py += PacketSize;
}
@ -183,10 +357,10 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY&
Packet xi1 = ei_ploadu(px+PacketSize);
Packet yi = ei_pload (py);
Packet yi1 = ei_pload (py+PacketSize);
ei_pstoreu(px, ei_psub(ei_pmul(pc,xi),cj.pmul(ps,yi)));
ei_pstoreu(px+PacketSize, ei_psub(ei_pmul(pc,xi1),cj.pmul(ps,yi1)));
ei_pstore (py, ei_padd(ei_pmul(ps,xi),ei_pmul(pc,yi)));
ei_pstore (py+PacketSize, ei_padd(ei_pmul(ps,xi1),ei_pmul(pc,yi1)));
ei_pstoreu(px, ei_padd(ei_pmul(pc,xi),cj.pmul(ps,yi)));
ei_pstoreu(px+PacketSize, ei_padd(ei_pmul(pc,xi1),cj.pmul(ps,yi1)));
ei_pstore (py, ei_psub(ei_pmul(pc,yi),ei_pmul(ps,xi)));
ei_pstore (py+PacketSize, ei_psub(ei_pmul(pc,yi1),ei_pmul(ps,xi1)));
px += Peeling*PacketSize;
py += Peeling*PacketSize;
}
@ -194,8 +368,8 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY&
{
Packet xi = ei_ploadu(x+peelingEnd);
Packet yi = ei_pload (y+peelingEnd);
ei_pstoreu(x+peelingEnd, ei_psub(ei_pmul(pc,xi),cj.pmul(ps,yi)));
ei_pstore (y+peelingEnd, ei_padd(ei_pmul(ps,xi),ei_pmul(pc,yi)));
ei_pstoreu(x+peelingEnd, ei_padd(ei_pmul(pc,xi),cj.pmul(ps,yi)));
ei_pstore (y+peelingEnd, ei_psub(ei_pmul(pc,yi),ei_pmul(ps,xi)));
}
}
@ -203,8 +377,8 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY&
{
Scalar xi = x[i];
Scalar yi = y[i];
x[i] = c * xi - ei_conj(s) * yi;
y[i] = s * xi + c * yi;
x[i] = j.c() * xi + ei_conj(j.s()) * yi;
y[i] = -j.s() * xi + ei_conj(j.c()) * yi;
}
}
else
@ -213,8 +387,8 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY&
{
Scalar xi = *x;
Scalar yi = *y;
*x = c * xi - ei_conj(s) * yi;
*y = s * xi + c * yi;
*x = j.c() * xi + ei_conj(j.s()) * yi;
*y = -j.s() * xi + ei_conj(j.c()) * yi;
x += incrx;
y += incry;
}

8
Eigen/src/QR/FullPivotingHouseholderQR.h
View File

@ -67,12 +67,10 @@ template<typename MatrixType> class FullPivotingHouseholderQR
* The default constructor is useful in cases in which the user intends to
* perform decompositions via FullPivotingHouseholderQR::compute(const MatrixType&).
*/
FullPivotingHouseholderQR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {}
FullPivotingHouseholderQR() : m_isInitialized(false) {}
FullPivotingHouseholderQR(const MatrixType& matrix)
: m_qr(matrix.rows(), matrix.cols()),
m_hCoeffs(std::min(matrix.rows(),matrix.cols())),
m_isInitialized(false)
: m_isInitialized(false)
{
compute(matrix);
}
@ -286,7 +284,7 @@ FullPivotingHouseholderQR<MatrixType>& FullPivotingHouseholderQR<MatrixType>::co
m_cols_permutation.resize(matrix.cols());
int number_of_transpositions = 0;
RealScalar biggest;
RealScalar biggest(0);
for (int k = 0; k < size; ++k)
{

5
Eigen/src/QR/QrInstantiations.cpp
View File

@ -33,11 +33,6 @@
namespace Eigen
{
template static void ei_tridiagonal_qr_step(float* , float* , int, int, float* , int);
template static void ei_tridiagonal_qr_step(double* , double* , int, int, double* , int);
template static void ei_tridiagonal_qr_step(float* , float* , int, int, std::complex<float>* , int);
template static void ei_tridiagonal_qr_step(double* , double* , int, int, std::complex<double>* , int);
EIGEN_QR_MODULE_INSTANTIATE();
}

352
Eigen/src/SVD/JacobiSVD.h Normal file
View File

@ -0,0 +1,352 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifndef EIGEN_JACOBISVD_H
#define EIGEN_JACOBISVD_H
/** \ingroup SVD_Module
* \nonstableyet
*
* \class JacobiSVD
*
* \brief Jacobi SVD decomposition of a square matrix
*
* \param MatrixType the type of the matrix of which we are computing the SVD decomposition
* \param Options a bit field of flags offering the following options: \c SkipU and \c SkipV allow to skip the computation of
* the unitaries \a U and \a V respectively; \c AtLeastAsManyRowsAsCols and \c AtLeastAsManyRowsAsCols allows
* to hint the compiler to only generate the corresponding code paths; \c Square is equivantent to the combination of
* the latter two bits and is useful when you know that the matrix is square. Note that when this information can
* be automatically deduced from the matrix type (e.g. a Matrix3f is always square), Eigen does it for you.
*
* \sa MatrixBase::jacobiSvd()
*/
template<typename MatrixType, unsigned int Options> class JacobiSVD
{
private:
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
enum {
ComputeU = (Options & SkipU) == 0,
ComputeV = (Options & SkipV) == 0,
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
DiagSizeAtCompileTime = EIGEN_SIZE_MIN(RowsAtCompileTime,ColsAtCompileTime),
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN(MaxRowsAtCompileTime,MaxColsAtCompileTime),
MatrixOptions = MatrixType::Options
};
typedef Matrix<Scalar, Dynamic, Dynamic, MatrixOptions> DummyMatrixType;
typedef typename ei_meta_if<ComputeU,
Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime,
MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime>,
DummyMatrixType>::ret MatrixUType;
typedef typename ei_meta_if<ComputeV,
Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime,
MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime>,
DummyMatrixType>::ret MatrixVType;
typedef Matrix<RealScalar, DiagSizeAtCompileTime, 1,
MatrixOptions, MaxDiagSizeAtCompileTime, 1> SingularValuesType;
typedef Matrix<Scalar, 1, RowsAtCompileTime, MatrixOptions, 1, MaxRowsAtCompileTime> RowType;
typedef Matrix<Scalar, RowsAtCompileTime, 1, MatrixOptions, MaxRowsAtCompileTime, 1> ColType;
typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime,
MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime>
WorkMatrixType;
public:
JacobiSVD() : m_isInitialized(false) {}
JacobiSVD(const MatrixType& matrix) : m_isInitialized(false)
{
compute(matrix);
}
JacobiSVD& compute(const MatrixType& matrix);
const MatrixUType& matrixU() const
{
ei_assert(m_isInitialized && "JacobiSVD is not initialized.");
return m_matrixU;
}
const SingularValuesType& singularValues() const
{
ei_assert(m_isInitialized && "JacobiSVD is not initialized.");
return m_singularValues;
}
const MatrixVType& matrixV() const
{
ei_assert(m_isInitialized && "JacobiSVD is not initialized.");
return m_matrixV;
}
protected:
MatrixUType m_matrixU;
MatrixVType m_matrixV;
SingularValuesType m_singularValues;
bool m_isInitialized;
template<typename _MatrixType, unsigned int _Options, bool _IsComplex>
friend struct ei_svd_precondition_2x2_block_to_be_real;
template<typename _MatrixType, unsigned int _Options, bool _PossiblyMoreRowsThanCols>
friend struct ei_svd_precondition_if_more_rows_than_cols;
template<typename _MatrixType, unsigned int _Options, bool _PossiblyMoreRowsThanCols>
friend struct ei_svd_precondition_if_more_cols_than_rows;
};
template<typename MatrixType, unsigned int Options, bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex>
struct ei_svd_precondition_2x2_block_to_be_real
{
typedef JacobiSVD<MatrixType, Options> SVD;
static void run(typename SVD::WorkMatrixType&, JacobiSVD<MatrixType, Options>&, int, int) {}
};
template<typename MatrixType, unsigned int Options>
struct ei_svd_precondition_2x2_block_to_be_real<MatrixType, Options, true>
{
typedef JacobiSVD<MatrixType, Options> SVD;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
enum { ComputeU = SVD::ComputeU, ComputeV = SVD::ComputeV };
static void run(typename SVD::WorkMatrixType& work_matrix, JacobiSVD<MatrixType, Options>& svd, int p, int q)
{
Scalar z;
PlanarRotation<Scalar> rot;
RealScalar n = ei_sqrt(ei_abs2(work_matrix.coeff(p,p)) + ei_abs2(work_matrix.coeff(q,p)));
if(n==0)
{
z = ei_abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
work_matrix.row(p) *= z;
if(ComputeU) svd.m_matrixU.col(p) *= ei_conj(z);
z = ei_abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
work_matrix.row(q) *= z;
if(ComputeU) svd.m_matrixU.col(q) *= ei_conj(z);
}
else
{
rot.c() = ei_conj(work_matrix.coeff(p,p)) / n;
rot.s() = work_matrix.coeff(q,p) / n;
work_matrix.applyOnTheLeft(p,q,rot);
if(ComputeU) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint());
if(work_matrix.coeff(p,q) != Scalar(0))
{
Scalar z = ei_abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
work_matrix.col(q) *= z;
if(ComputeV) svd.m_matrixV.col(q) *= z;
}
if(work_matrix.coeff(q,q) != Scalar(0))
{
z = ei_abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
work_matrix.row(q) *= z;
if(ComputeU) svd.m_matrixU.col(q) *= ei_conj(z);
}
}
}
};
template<typename MatrixType, typename RealScalar>
void ei_real_2x2_jacobi_svd(const MatrixType& matrix, int p, int q,
PlanarRotation<RealScalar> *j_left,
PlanarRotation<RealScalar> *j_right)
{
Matrix<RealScalar,2,2> m;
m << ei_real(matrix.coeff(p,p)), ei_real(matrix.coeff(p,q)),
ei_real(matrix.coeff(q,p)), ei_real(matrix.coeff(q,q));
PlanarRotation<RealScalar> rot1;
RealScalar t = m.coeff(0,0) + m.coeff(1,1);
RealScalar d = m.coeff(1,0) - m.coeff(0,1);
if(t == RealScalar(0))
{
rot1.c() = 0;
rot1.s() = d > 0 ? 1 : -1;
}
else
{
RealScalar u = d / t;
rot1.c() = RealScalar(1) / ei_sqrt(1 + ei_abs2(u));
rot1.s() = rot1.c() * u;
}
m.applyOnTheLeft(0,1,rot1);
j_right->makeJacobi(m,0,1);
*j_left = rot1 * j_right->transpose();
}
template<typename MatrixType, unsigned int Options,
bool PossiblyMoreRowsThanCols = (Options & AtLeastAsManyColsAsRows) == 0
&& (MatrixType::RowsAtCompileTime==Dynamic
|| MatrixType::RowsAtCompileTime>MatrixType::ColsAtCompileTime)>
struct ei_svd_precondition_if_more_rows_than_cols
{
typedef JacobiSVD<MatrixType, Options> SVD;
static bool run(const MatrixType&, typename SVD::WorkMatrixType&, JacobiSVD<MatrixType, Options>&) { return false; }
};
template<typename MatrixType, unsigned int Options>
struct ei_svd_precondition_if_more_rows_than_cols<MatrixType, Options, true>
{
typedef JacobiSVD<MatrixType, Options> SVD;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
enum { ComputeU = SVD::ComputeU, ComputeV = SVD::ComputeV };
static bool run(const MatrixType& matrix, typename SVD::WorkMatrixType& work_matrix, SVD& svd)
{
int rows = matrix.rows();
int cols = matrix.cols();
int diagSize = cols;
if(rows > cols)
{
FullPivotingHouseholderQR<MatrixType> qr(matrix);
work_matrix = qr.matrixQR().block(0,0,diagSize,diagSize).template triangularView<UpperTriangular>();
if(ComputeU) svd.m_matrixU = qr.matrixQ();
if(ComputeV)
for(int i = 0; i < cols; i++)
svd.m_matrixV.coeffRef(qr.colsPermutation().coeff(i),i) = Scalar(1);
return true;
}
else return false;
}
};
template<typename MatrixType, unsigned int Options,
bool PossiblyMoreColsThanRows = (Options & AtLeastAsManyRowsAsCols) == 0
&& (MatrixType::ColsAtCompileTime==Dynamic
|| MatrixType::ColsAtCompileTime>MatrixType::RowsAtCompileTime)>
struct ei_svd_precondition_if_more_cols_than_rows
{
typedef JacobiSVD<MatrixType, Options> SVD;
static bool run(const MatrixType&, typename SVD::WorkMatrixType&, JacobiSVD<MatrixType, Options>&) { return false; }
};
template<typename MatrixType, unsigned int Options>
struct ei_svd_precondition_if_more_cols_than_rows<MatrixType, Options, true>
{
typedef JacobiSVD<MatrixType, Options> SVD;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
enum {
ComputeU = SVD::ComputeU,
ComputeV = SVD::ComputeV,
RowsAtCompileTime = SVD::RowsAtCompileTime,
ColsAtCompileTime = SVD::ColsAtCompileTime,
MaxRowsAtCompileTime = SVD::MaxRowsAtCompileTime,
MaxColsAtCompileTime = SVD::MaxColsAtCompileTime,
MatrixOptions = SVD::MatrixOptions
};
static bool run(const MatrixType& matrix, typename SVD::WorkMatrixType& work_matrix, SVD& svd)
{
int rows = matrix.rows();
int cols = matrix.cols();
int diagSize = rows;
if(cols > rows)
{
typedef Matrix<Scalar,ColsAtCompileTime,RowsAtCompileTime,
MatrixOptions,MaxColsAtCompileTime,MaxRowsAtCompileTime>
TransposeTypeWithSameStorageOrder;
FullPivotingHouseholderQR<TransposeTypeWithSameStorageOrder> qr(matrix.adjoint());
work_matrix = qr.matrixQR().block(0,0,diagSize,diagSize).template triangularView<UpperTriangular>().adjoint();
if(ComputeV) svd.m_matrixV = qr.matrixQ();
if(ComputeU)
for(int i = 0; i < rows; i++)
svd.m_matrixU.coeffRef(qr.colsPermutation().coeff(i),i) = Scalar(1);
return true;
}
else return false;
}
};
template<typename MatrixType, unsigned int Options>
JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute(const MatrixType& matrix)
{
WorkMatrixType work_matrix;
int rows = matrix.rows();
int cols = matrix.cols();
int diagSize = std::min(rows, cols);
if(ComputeU) m_matrixU = MatrixUType::Zero(rows,rows);
if(ComputeV) m_matrixV = MatrixVType::Zero(cols,cols);
m_singularValues.resize(diagSize);
const RealScalar precision = 2 * epsilon<Scalar>();
if(!ei_svd_precondition_if_more_rows_than_cols<MatrixType, Options>::run(matrix, work_matrix, *this)
&& !ei_svd_precondition_if_more_cols_than_rows<MatrixType, Options>::run(matrix, work_matrix, *this))
{
work_matrix = matrix.block(0,0,diagSize,diagSize);
if(ComputeU) m_matrixU.diagonal().setOnes();
if(ComputeV) m_matrixV.diagonal().setOnes();
}
bool finished = false;
while(!finished)
{
finished = true;
for(int p = 1; p < diagSize; ++p)
{
for(int q = 0; q < p; ++q)
{
if(std::max(ei_abs(work_matrix.coeff(p,q)),ei_abs(work_matrix.coeff(q,p)))
> std::max(ei_abs(work_matrix.coeff(p,p)),ei_abs(work_matrix.coeff(q,q)))*precision)
{
finished = false;
ei_svd_precondition_2x2_block_to_be_real<MatrixType, Options>::run(work_matrix, *this, p, q);
PlanarRotation<RealScalar> j_left, j_right;
ei_real_2x2_jacobi_svd(work_matrix, p, q, &j_left, &j_right);
work_matrix.applyOnTheLeft(p,q,j_left);
if(ComputeU) m_matrixU.applyOnTheRight(p,q,j_left.transpose());
work_matrix.applyOnTheRight(p,q,j_right);
if(ComputeV) m_matrixV.applyOnTheRight(p,q,j_right);
}
}
}
}
for(int i = 0; i < diagSize; ++i)
{
RealScalar a = ei_abs(work_matrix.coeff(i,i));
m_singularValues.coeffRef(i) = a;
if(ComputeU && (a!=RealScalar(0))) m_matrixU.col(i) *= work_matrix.coeff(i,i)/a;
}
for(int i = 0; i < diagSize; i++)
{
int pos;
m_singularValues.end(diagSize-i).maxCoeff(&pos);
if(pos)
{
pos += i;
std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos));
if(ComputeU) m_matrixU.col(pos).swap(m_matrixU.col(i));
if(ComputeV) m_matrixV.col(pos).swap(m_matrixV.col(i));
}
}
m_isInitialized = true;
return *this;
}
#endif // EIGEN_JACOBISVD_H

169
Eigen/src/SVD/JacobiSquareSVD.h
View File

@ -1,169 +0,0 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifndef EIGEN_JACOBISQUARESVD_H
#define EIGEN_JACOBISQUARESVD_H
/** \ingroup SVD_Module
* \nonstableyet
*
* \class JacobiSquareSVD
*
* \brief Jacobi SVD decomposition of a square matrix
*
* \param MatrixType the type of the matrix of which we are computing the SVD decomposition
* \param ComputeU whether the U matrix should be computed
* \param ComputeV whether the V matrix should be computed
*
* \sa MatrixBase::jacobiSvd()
*/
template<typename MatrixType, bool ComputeU, bool ComputeV> class JacobiSquareSVD
{
private:
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
Options = MatrixType::Options
};
typedef Matrix<Scalar, Dynamic, Dynamic, Options> DummyMatrixType;
typedef typename ei_meta_if<ComputeU,
Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime,
Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime>,
DummyMatrixType>::ret MatrixUType;
typedef typename Diagonal<MatrixType,0>::PlainMatrixType SingularValuesType;
typedef Matrix<Scalar, 1, RowsAtCompileTime, Options, 1, MaxRowsAtCompileTime> RowType;
typedef Matrix<Scalar, RowsAtCompileTime, 1, Options, MaxRowsAtCompileTime, 1> ColType;
public:
JacobiSquareSVD() : m_isInitialized(false) {}
JacobiSquareSVD(const MatrixType& matrix) : m_isInitialized(false)
{
compute(matrix);
}
JacobiSquareSVD& compute(const MatrixType& matrix);
const MatrixUType& matrixU() const
{
ei_assert(m_isInitialized && "SVD is not initialized.");
return m_matrixU;
}
const SingularValuesType& singularValues() const
{
ei_assert(m_isInitialized && "SVD is not initialized.");
return m_singularValues;
}
const MatrixUType& matrixV() const
{
ei_assert(m_isInitialized && "SVD is not initialized.");
return m_matrixV;
}
protected:
MatrixUType m_matrixU;
MatrixUType m_matrixV;
SingularValuesType m_singularValues;
bool m_isInitialized;
};
template<typename MatrixType, bool ComputeU, bool ComputeV>
JacobiSquareSVD<MatrixType, ComputeU, ComputeV>& JacobiSquareSVD<MatrixType, ComputeU, ComputeV>::compute(const MatrixType& matrix)
{
MatrixType work_matrix(matrix);
int size = matrix.rows();
if(ComputeU) m_matrixU = MatrixUType::Identity(size,size);
if(ComputeV) m_matrixV = MatrixUType::Identity(size,size);
m_singularValues.resize(size);
const RealScalar precision = 2 * epsilon<Scalar>();
sweep_again:
for(int p = 1; p < size; ++p)
{
for(int q = 0; q < p; ++q)
{
Scalar c, s;
while(std::max(ei_abs(work_matrix.coeff(p,q)),ei_abs(work_matrix.coeff(q,p)))
> std::max(ei_abs(work_matrix.coeff(p,p)),ei_abs(work_matrix.coeff(q,q)))*precision)
{
if(work_matrix.makeJacobiForAtA(p,q,&c,&s))
{
work_matrix.applyJacobiOnTheRight(p,q,c,s);
if(ComputeV) m_matrixV.applyJacobiOnTheRight(p,q,c,s);
}
if(work_matrix.makeJacobiForAAt(p,q,&c,&s))
{
ei_normalizeJacobi(&c, &s, work_matrix.coeff(p,p), work_matrix.coeff(q,p)),
work_matrix.applyJacobiOnTheLeft(p,q,c,s);
if(ComputeU) m_matrixU.applyJacobiOnTheRight(p,q,c,s);
}
}
}
}
RealScalar biggestOnDiag = work_matrix.diagonal().cwise().abs().maxCoeff();
RealScalar maxAllowedOffDiag = biggestOnDiag * precision;
for(int p = 0; p < size; ++p)
{
for(int q = 0; q < p; ++q)
if(ei_abs(work_matrix.coeff(p,q)) > maxAllowedOffDiag)
goto sweep_again;
for(int q = p+1; q < size; ++q)
if(ei_abs(work_matrix.coeff(p,q)) > maxAllowedOffDiag)
goto sweep_again;
}
m_singularValues = work_matrix.diagonal().cwise().abs();
RealScalar biggestSingularValue = m_singularValues.maxCoeff();
for(int i = 0; i < size; ++i)
{
RealScalar a = ei_abs(work_matrix.coeff(i,i));
m_singularValues.coeffRef(i) = a;
if(ComputeU && !ei_isMuchSmallerThan(a, biggestSingularValue)) m_matrixU.col(i) *= work_matrix.coeff(i,i)/a;
}
for(int i = 0; i < size; i++)
{
int pos;
m_singularValues.end(size-i).maxCoeff(&pos);
if(pos)
{
pos += i;
std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos));
if(ComputeU) m_matrixU.col(pos).swap(m_matrixU.col(i));
if(ComputeV) m_matrixV.col(pos).swap(m_matrixV.col(i));
}
}
m_isInitialized = true;
return *this;
}
#endif // EIGEN_JACOBISQUARESVD_H

10
Eigen/src/SVD/SVD.h
View File

@ -309,7 +309,7 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
h = Scalar(1.0)/h;
c = g*h;
s = -f*h;
V.applyJacobiOnTheRight(i,nm,c,s);
V.applyOnTheRight(i,nm,PlanarRotation<Scalar>(c,s));
}
}
z = W[k];
@ -342,7 +342,7 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
y = W[i];
h = s*g;
g = c*g;
z = pythag(f,h);
rv1[j] = z;
c = f/z;
@ -351,8 +351,8 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
g = g*c - x*s;
h = y*s;
y *= c;
V.applyJacobiOnTheRight(i,j,c,s);
V.applyOnTheRight(i,j,PlanarRotation<Scalar>(c,s));
z = pythag(f,h);
W[j] = z;
// Rotation can be arbitrary if z = 0.
@ -364,7 +364,7 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
}
f = c*g + s*y;
x = c*y - s*g;
A.applyJacobiOnTheRight(i,j,c,s);
A.applyOnTheRight(i,j,PlanarRotation<Scalar>(c,s));
}
rv1[l] = 0.0;
rv1[k] = f;

4
Eigen/src/Sparse/AmbiVector.h
View File

@ -41,7 +41,7 @@ template<typename _Scalar> class AmbiVector
resize(size);
}
void init(RealScalar estimatedDensity);
void init(double estimatedDensity);
void init(int mode);
int nonZeros() const;
@ -143,7 +143,7 @@ int AmbiVector<Scalar>::nonZeros() const
}
template<typename Scalar>
void AmbiVector<Scalar>::init(RealScalar estimatedDensity)
void AmbiVector<Scalar>::init(double estimatedDensity)
{
if (estimatedDensity>0.1)
init(IsDense);

15
Eigen/src/Sparse/CholmodSupport.h
View File

@ -99,7 +99,7 @@ cholmod_dense ei_cholmod_map_eigen_to_dense(MatrixBase<Derived>& mat)
res.nrow = mat.rows();
res.ncol = mat.cols();
res.nzmax = res.nrow * res.ncol;
res.d = mat.derived().stride();
res.d = Derived::IsVectorAtCompileTime ? mat.derived().size() : mat.derived().stride();
res.x = mat.derived().data();
res.z = 0;
@ -157,7 +157,7 @@ class SparseLLT<MatrixType,Cholmod> : public SparseLLT<MatrixType>
inline const typename Base::CholMatrixType& matrixL(void) const;
template<typename Derived>
void solveInPlace(MatrixBase<Derived> &b) const;
bool solveInPlace(MatrixBase<Derived> &b) const;
void compute(const MatrixType& matrix);
@ -216,7 +216,7 @@ SparseLLT<MatrixType,Cholmod>::matrixL() const
template<typename MatrixType>
template<typename Derived>
void SparseLLT<MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const
bool SparseLLT<MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const
{
const int size = m_cholmodFactor->n;
ei_assert(size==b.rows());
@ -228,9 +228,16 @@ void SparseLLT<MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const
// as long as our own triangular sparse solver is not fully optimal,
// let's use CHOLMOD's one:
cholmod_dense cdb = ei_cholmod_map_eigen_to_dense(b);
cholmod_dense* x = cholmod_solve(CHOLMOD_LDLt, m_cholmodFactor, &cdb, &m_cholmod);
//cholmod_dense* x = cholmod_solve(CHOLMOD_LDLt, m_cholmodFactor, &cdb, &m_cholmod);
cholmod_dense* x = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &cdb, &m_cholmod);
if(!x)
{
std::cerr << "Eigen: cholmod_solve failed\n";
return false;
}
b = Matrix<typename Base::Scalar,Dynamic,1>::Map(reinterpret_cast<typename Base::Scalar*>(x->x),b.rows());
cholmod_free_dense(&x, &m_cholmod);
return true;
}
#endif // EIGEN_CHOLMODSUPPORT_H

124
Eigen/src/Sparse/SuperLUSupport.h
View File

@ -48,6 +48,37 @@ DECL_GSSVX(SuperLU_C,cgssvx,float,std::complex<float>)
DECL_GSSVX(SuperLU_D,dgssvx,double,double)
DECL_GSSVX(SuperLU_Z,zgssvx,double,std::complex<double>)
#ifdef MILU_ALPHA
#define EIGEN_SUPERLU_HAS_ILU
#endif
#ifdef EIGEN_SUPERLU_HAS_ILU
// similarly for the incomplete factorization using gsisx
#define DECL_GSISX(NAMESPACE,FNAME,FLOATTYPE,KEYTYPE) \
inline float SuperLU_gsisx(superlu_options_t *options, SuperMatrix *A, \
int *perm_c, int *perm_r, int *etree, char *equed, \
FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \
SuperMatrix *U, void *work, int lwork, \
SuperMatrix *B, SuperMatrix *X, \
FLOATTYPE *recip_pivot_growth, \
FLOATTYPE *rcond, \
SuperLUStat_t *stats, int *info, KEYTYPE) { \
using namespace NAMESPACE; \
mem_usage_t mem_usage; \
NAMESPACE::FNAME(options, A, perm_c, perm_r, etree, equed, R, C, L, \
U, work, lwork, B, X, recip_pivot_growth, rcond, \
&mem_usage, stats, info); \
return mem_usage.for_lu; /* bytes used by the factor storage */ \
}
DECL_GSISX(SuperLU_S,sgsisx,float,float)
DECL_GSISX(SuperLU_C,cgsisx,float,std::complex<float>)
DECL_GSISX(SuperLU_D,dgsisx,double,double)
DECL_GSISX(SuperLU_Z,zgsisx,double,std::complex<double>)
#endif
template<typename MatrixType>
struct SluMatrixMapHelper;
@ -71,7 +102,7 @@ struct SluMatrix : SuperMatrix
Store = &storage;
storage = other.storage;
}
SluMatrix& operator=(const SluMatrix& other)
{
SuperMatrix::operator=(static_cast<const SuperMatrix&>(other));
@ -130,7 +161,7 @@ struct SluMatrix : SuperMatrix
res.nrow = mat.rows();
res.ncol = mat.cols();
res.storage.lda = mat.stride();
res.storage.lda = MatrixType::IsVectorAtCompileTime ? mat.size() : mat.stride();
res.storage.values = mat.data();
return res;
}
@ -373,7 +404,7 @@ void SparseLU<MatrixType,SuperLU>::compute(const MatrixType& a)
m_sluA = m_matrix.asSluMatrix();
memset(&m_sluL,0,sizeof m_sluL);
memset(&m_sluU,0,sizeof m_sluU);
m_sluEqued = 'B';
//m_sluEqued = 'B';
int info = 0;
m_p.resize(size);
@ -395,14 +426,44 @@ void SparseLU<MatrixType,SuperLU>::compute(const MatrixType& a)
m_sluX = m_sluB;
StatInit(&m_sluStat);
SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
&m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
&m_sluL, &m_sluU,
NULL, 0,
&m_sluB, &m_sluX,
&recip_pivot_gross, &rcond,
&ferr, &berr,
&m_sluStat, &info, Scalar());
if (m_flags&IncompleteFactorization)
{
#ifdef EIGEN_SUPERLU_HAS_ILU
ilu_set_default_options(&m_sluOptions);
// no attempt to preserve column sum
m_sluOptions.ILU_MILU = SILU;
// only basic ILU(k) support -- no direct control over memory consumption
// better to use ILU_DropRule = DROP_BASIC | DROP_AREA
// and set ILU_FillFactor to max memory growth
m_sluOptions.ILU_DropRule = DROP_BASIC;
m_sluOptions.ILU_DropTol = Base::m_precision;
SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
&m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
&m_sluL, &m_sluU,
NULL, 0,
&m_sluB, &m_sluX,
&recip_pivot_gross, &rcond,
&m_sluStat, &info, Scalar());
#else
std::cerr << "Incomplete factorization is only available in SuperLU v4\n";
Base::m_succeeded = false;
return;
#endif
}
else
{
SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
&m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
&m_sluL, &m_sluU,
NULL, 0,
&m_sluB, &m_sluX,
&recip_pivot_gross, &rcond,
&ferr, &berr,
&m_sluStat, &info, Scalar());
}
StatFree(&m_sluStat);
m_extractedDataAreDirty = true;
@ -440,17 +501,36 @@ bool SparseLU<MatrixType,SuperLU>::solve(const MatrixBase<BDerived> &b,
StatInit(&m_sluStat);
int info = 0;
RealScalar recip_pivot_gross, rcond;
SuperLU_gssvx(
&m_sluOptions, &m_sluA,
m_q.data(), m_p.data(),
&m_sluEtree[0], &m_sluEqued,
&m_sluRscale[0], &m_sluCscale[0],
&m_sluL, &m_sluU,
NULL, 0,
&m_sluB, &m_sluX,
&recip_pivot_gross, &rcond,
&m_sluFerr[0], &m_sluBerr[0],
&m_sluStat, &info, Scalar());
if (m_flags&IncompleteFactorization)
{
#ifdef EIGEN_SUPERLU_HAS_ILU
SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
&m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
&m_sluL, &m_sluU,
NULL, 0,
&m_sluB, &m_sluX,
&recip_pivot_gross, &rcond,
&m_sluStat, &info, Scalar());
#else
std::cerr << "Incomplete factorization is only available in SuperLU v4\n";
return false;
#endif
}
else
{
SuperLU_gssvx(
&m_sluOptions, &m_sluA,
m_q.data(), m_p.data(),
&m_sluEtree[0], &m_sluEqued,
&m_sluRscale[0], &m_sluCscale[0],
&m_sluL, &m_sluU,
NULL, 0,
&m_sluB, &m_sluX,
&recip_pivot_gross, &rcond,
&m_sluFerr[0], &m_sluBerr[0],
&m_sluStat, &info, Scalar());
}
StatFree(&m_sluStat);
// reset to previous state

14
cmake/EigenTesting.cmake
View File

@ -153,11 +153,17 @@ macro(ei_init_testing)
endmacro(ei_init_testing)
if(CMAKE_COMPILER_IS_GNUCXX)
option(EIGEN_COVERAGE_TESTING "Enable/disable gcov" OFF)
if(EIGEN_COVERAGE_TESTING)
set(COVERAGE_FLAGS "-fprofile-arcs -ftest-coverage")
else(EIGEN_COVERAGE_TESTING)
set(COVERAGE_FLAGS "")
endif(EIGEN_COVERAGE_TESTING)
if(CMAKE_SYSTEM_NAME MATCHES Linux)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g2")
set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "${CMAKE_CXX_FLAGS_RELWITHDEBINFO} -O2 -g2")
set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -fno-inline-functions")
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -O0 -g2")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${COVERAGE_FLAGS} -g2")
set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "${CMAKE_CXX_FLAGS_RELWITHDEBINFO} ${COVERAGE_FLAGS} -O2 -g2")
set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} ${COVERAGE_FLAGS} -fno-inline-functions")
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} ${COVERAGE_FLAGS} -O0 -g2")
endif(CMAKE_SYSTEM_NAME MATCHES Linux)
set(EI_OFLAG "-O2")
elseif(MSVC)

14
doc/C01_QuickStartGuide.dox
View File

@ -129,7 +129,7 @@ The default constructor leaves coefficients uninitialized. Any dynamic size is s
Matrix3f A; // construct 3x3 matrix with uninitialized coefficients
A(0,0) = 5; // OK
MatrixXf B; // construct 0x0 matrix without allocating anything
A(0,0) = 5; // Error, B is uninitialized, doesn't have any coefficients to address
B(0,0) = 5; // Error, B is uninitialized, doesn't have any coefficients to address
\endcode
In the above example, B is an uninitialized matrix. What to do with such a matrix? You can call resize() on it, or you can assign another matrix to it. Like this:
@ -261,7 +261,7 @@ v = 6 6 6
\subsection TutorialCasting Casting
In Eigen, any matrices of same size and same scalar type are all naturally compatible. The scalar type can be explicitely casted to another one using the template MatrixBase::cast() function:
In Eigen, any matrices of same size and same scalar type are all naturally compatible. The scalar type can be explicitly casted to another one using the template MatrixBase::cast() function:
\code
Matrix3d md(1,2,3);
Matrix3f mf = md.cast<float>();
@ -328,7 +328,7 @@ In short, all arithmetic operators can be used right away as in the following ex
mat4 -= mat1*1.5 + mat2 * (mat3/4);
\endcode
which includes two matrix scalar products ("mat1*1.5" and "mat3/4"), a matrix-matrix product ("mat2 * (mat3/4)"),
a matrix addition ("+") and substraction with assignment ("-=").
a matrix addition ("+") and subtraction with assignment ("-=").
<table class="tutorial_code">
<tr><td>
@ -464,7 +464,7 @@ mat = 2 7 8
Also note that maxCoeff and minCoeff can takes optional arguments returning the coordinates of the respective min/max coeff: \link MatrixBase::maxCoeff(int*,int*) const maxCoeff(int* i, int* j) \endlink, \link MatrixBase::minCoeff(int*,int*) const minCoeff(int* i, int* j) \endlink.
<span class="note">\b Side \b note: The all() and any() functions are especially useful in combinaison with coeff-wise comparison operators (\ref CwiseAll "example").</span>
<span class="note">\b Side \b note: The all() and any() functions are especially useful in combination with coeff-wise comparison operators (\ref CwiseAll "example").</span>
@ -578,7 +578,7 @@ vec1.normalize();\endcode
<a href="#" class="top">top</a>\section TutorialCoreTriangularMatrix Dealing with triangular matrices
Currently, Eigen does not provide any explcit triangular matrix, with storage class. Instead, we
Currently, Eigen does not provide any explicit triangular matrix, with storage class. Instead, we
can reference a triangular part of a square matrix or expression to perform special treatment on it.
This is achieved by the class TriangularView and the MatrixBase::triangularView template function.
Note that the opposite triangular part of the matrix is never referenced, and so it can, e.g., store
@ -595,12 +595,12 @@ m.triangularView<Eigen::LowerTriangular>()
m.triangularView<Eigen::UnitLowerTriangular>()\endcode
</td></tr>
<tr><td>
Writting to a specific triangular part:\n (only the referenced triangular part is evaluated)
Writing to a specific triangular part:\n (only the referenced triangular part is evaluated)
</td><td>\code
m1.triangularView<Eigen::LowerTriangular>() = m2 + m3 \endcode
</td></tr>
<tr><td>
Convertion to a dense matrix setting the opposite triangular part to zero:
Conversion to a dense matrix setting the opposite triangular part to zero:
</td><td>\code
m2 = m1.triangularView<Eigen::UnitUpperTriangular>()\endcode
</td></tr>

3
doc/Doxyfile.in
View File

@ -208,10 +208,11 @@ ALIASES = "only_for_vectors=This is only for vectors (either row-
"lu_module=This is defined in the %LU module. \code #include <Eigen/LU> \endcode" \
"cholesky_module=This is defined in the %Cholesky module. \code #include <Eigen/Cholesky> \endcode" \
"qr_module=This is defined in the %QR module. \code #include <Eigen/QR> \endcode" \
"jacobi_module=This is defined in the %Jacobi module. \code #include <Eigen/Jacobi> \endcode" \
"svd_module=This is defined in the %SVD module. \code #include <Eigen/SVD> \endcode" \
"geometry_module=This is defined in the %Geometry module. \code #include <Eigen/Geometry> \endcode" \
"leastsquares_module=This is defined in the %LeastSquares module. \code #include <Eigen/LeastSquares> \endcode" \
"addexample=\anchor" \
"eigenvalues_module=This is defined in the %Eigenvalues module. \code #include <Eigen/Eigenvalues> \endcode" \
"label=\bug" \
"redstar=<a href='#warningarraymodule' style='color:red;text-decoration: none;'>*</a>" \
"nonstableyet=\warning This is not considered to be part of the stable public API yet. Changes may happen in future releases. See \ref Experimental \"Experimental parts of Eigen\""

6
doc/snippets/Jacobi_makeGivens.cpp Normal file
View File

@ -0,0 +1,6 @@
Vector2f v = Vector2f::Random();
PlanarRotation<float> G;
G.makeGivens(v.x(), v.y());
cout << "Here is the vector v:" << endl << v << endl;
v.applyOnTheLeft(0, 1, G.adjoint());
cout << "Here is the vector J' * v:" << endl << v << endl;

8
doc/snippets/Jacobi_makeJacobi.cpp Normal file
View File

@ -0,0 +1,8 @@
Matrix2f m = Matrix2f::Random();
m = (m + m.adjoint()).eval();
PlanarRotation<float> J;
J.makeJacobi(m, 0, 1);
cout << "Here is the matrix m:" << endl << m << endl;
m.applyOnTheLeft(0, 1, J.adjoint());
m.applyOnTheRight(0, 1, J);
cout << "Here is the matrix J' * m * J:" << endl << m << endl;

1
doc/snippets/compile_snippet.cpp.in
View File

@ -4,6 +4,7 @@
#include <Eigen/QR>
#include <Eigen/Cholesky>
#include <Eigen/Geometry>
#include <Eigen/Jacobi>
using namespace Eigen;
using namespace std;

5
test/CMakeLists.txt
View File

@ -115,6 +115,7 @@ ei_add_test(product_trmv ${EI_OFLAG})
ei_add_test(product_trmm ${EI_OFLAG})
ei_add_test(product_trsm ${EI_OFLAG})
ei_add_test(product_notemporary ${EI_OFLAG})
ei_add_test(stable_norm)
ei_add_test(bandmatrix)
ei_add_test(cholesky " " "${GSL_LIBRARIES}")
ei_add_test(lu ${EI_OFLAG})
@ -125,7 +126,9 @@ ei_add_test(qr_colpivoting)
ei_add_test(qr_fullpivoting)
ei_add_test(eigensolver_selfadjoint " " "${GSL_LIBRARIES}")
ei_add_test(eigensolver_generic " " "${GSL_LIBRARIES}")
ei_add_test(eigensolver_complex)
ei_add_test(svd)
ei_add_test(jacobisvd ${EI_OFLAG})
ei_add_test(geo_orthomethods)
ei_add_test(geo_homogeneous)
ei_add_test(geo_quaternion)
@ -146,6 +149,8 @@ ei_add_test(sparse_product)
ei_add_test(sparse_solvers " " "${SPARSE_LIBS}")
ei_add_test(umeyama)
ei_add_test(householder)
ei_add_test(swap)
ei_add_test(conservative_resize)
ei_add_property(EIGEN_TESTING_SUMMARY "CXX: ${CMAKE_CXX_COMPILER}\n")
if(CMAKE_COMPILER_IS_GNUCXX)

9
test/adjoint.cpp
View File

@ -72,13 +72,6 @@ template<typename MatrixType> void adjoint(const MatrixType& m)
if(NumTraits<Scalar>::HasFloatingPoint)
VERIFY_IS_APPROX(v1.squaredNorm(), v1.norm() * v1.norm());
VERIFY_IS_MUCH_SMALLER_THAN(ei_abs(vzero.dot(v1)), static_cast<RealScalar>(1));
if(NumTraits<Scalar>::HasFloatingPoint)
{
VERIFY_IS_MUCH_SMALLER_THAN(vzero.norm(), static_cast<RealScalar>(1));
VERIFY_IS_APPROX(v1.norm(), v1.stableNorm());
VERIFY_IS_APPROX(v1.blueNorm(), v1.stableNorm());
VERIFY_IS_APPROX(v1.hypotNorm(), v1.stableNorm());
}
// check compatibility of dot and adjoint
VERIFY(ei_isApprox(v1.dot(square * v2), (square.adjoint() * v1).dot(v2), largerEps));
@ -124,7 +117,7 @@ void test_adjoint()
}
// test a large matrix only once
CALL_SUBTEST( adjoint(Matrix<float, 100, 100>()) );
{
MatrixXcf a(10,10), b(10,10);
VERIFY_RAISES_ASSERT(a = a.transpose());

12
test/bandmatrix.cpp
View File

@ -44,21 +44,21 @@ template<typename MatrixType> void bandmatrix(const MatrixType& _m)
dm1.diagonal().setConstant(123);
for (int i=1; i<=m.supers();++i)
{
m.diagonal(i).setConstant(i);
dm1.diagonal(i).setConstant(i);
m.diagonal(i).setConstant(static_cast<RealScalar>(i));
dm1.diagonal(i).setConstant(static_cast<RealScalar>(i));
}
for (int i=1; i<=m.subs();++i)
{
m.diagonal(-i).setConstant(-i);
dm1.diagonal(-i).setConstant(-i);
m.diagonal(-i).setConstant(-static_cast<RealScalar>(i));
dm1.diagonal(-i).setConstant(-static_cast<RealScalar>(i));
}
//std::cerr << m.m_data << "\n\n" << m.toDense() << "\n\n" << dm1 << "\n\n\n\n";
VERIFY_IS_APPROX(dm1,m.toDense());
for (int i=0; i<cols; ++i)
{
m.col(i).setConstant(i+1);
dm1.col(i).setConstant(i+1);
m.col(i).setConstant(static_cast<RealScalar>(i+1));
dm1.col(i).setConstant(static_cast<RealScalar>(i+1));
}
int d = std::min(rows,cols);
int a = std::max(0,cols-d-supers);

9
test/basicstuff.cpp
View File

@ -99,15 +99,6 @@ template<typename MatrixType> void basicStuff(const MatrixType& m)
MatrixType m4;
VERIFY_IS_APPROX(m4 = m1,m1);
// test swap
m3 = m1;
m1.swap(m2);
VERIFY_IS_APPROX(m3, m2);
if(rows*cols>=3)
{
VERIFY_IS_NOT_APPROX(m3, m1);
}
m3.real() = m1.real();
VERIFY_IS_APPROX(static_cast<const MatrixType&>(m3).real(), static_cast<const MatrixType&>(m1).real());
VERIFY_IS_APPROX(static_cast<const MatrixType&>(m3).real(), m1.real());

18
test/cholesky.cpp
View File

@ -82,7 +82,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
// // test gsl itself !
// VERIFY_IS_APPROX(vecB, _vecB);
// VERIFY_IS_APPROX(vecX, _vecX);
//
//
// Gsl::free(gMatA);
// Gsl::free(gSymm);
// Gsl::free(gVecB);
@ -149,16 +149,16 @@ void test_cholesky()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( cholesky(Matrix<double,1,1>()) );
// CALL_SUBTEST( cholesky(MatrixXd(1,1)) );
// CALL_SUBTEST( cholesky(Matrix2d()) );
// CALL_SUBTEST( cholesky(Matrix3f()) );
// CALL_SUBTEST( cholesky(Matrix4d()) );
CALL_SUBTEST( cholesky(MatrixXd(1,1)) );
CALL_SUBTEST( cholesky(Matrix2d()) );
CALL_SUBTEST( cholesky(Matrix3f()) );
CALL_SUBTEST( cholesky(Matrix4d()) );
CALL_SUBTEST( cholesky(MatrixXd(200,200)) );
CALL_SUBTEST( cholesky(MatrixXcd(100,100)) );
}
// CALL_SUBTEST( cholesky_verify_assert<Matrix3f>() );
// CALL_SUBTEST( cholesky_verify_assert<Matrix3d>() );
// CALL_SUBTEST( cholesky_verify_assert<MatrixXf>() );
// CALL_SUBTEST( cholesky_verify_assert<MatrixXd>() );
CALL_SUBTEST( cholesky_verify_assert<Matrix3f>() );
CALL_SUBTEST( cholesky_verify_assert<Matrix3d>() );
CALL_SUBTEST( cholesky_verify_assert<MatrixXf>() );
CALL_SUBTEST( cholesky_verify_assert<MatrixXd>() );
}

129
test/conservative_resize.cpp Normal file
View File

@ -0,0 +1,129 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra. Eigen itself is part of the KDE project.
//
// Copyright (C) 2009 Hauke Heibel <hauke.heibel@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or1 FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#include "main.h"
#include <Eigen/Core>
#include <Eigen/Array>
using namespace Eigen;
template <typename Scalar, int Storage>
void run_matrix_tests()
{
typedef Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Storage> MatrixType;
MatrixType m, n;
// boundary cases ...
m = n = MatrixType::Random(50,50);
m.conservativeResize(1,50);
VERIFY_IS_APPROX(m, n.block(0,0,1,50));
m = n = MatrixType::Random(50,50);
m.conservativeResize(50,1);
VERIFY_IS_APPROX(m, n.block(0,0,50,1));
m = n = MatrixType::Random(50,50);
m.conservativeResize(50,50);
VERIFY_IS_APPROX(m, n.block(0,0,50,50));
// random shrinking ...
for (int i=0; i<25; ++i)
{
const int rows = ei_random<int>(1,50);
const int cols = ei_random<int>(1,50);
m = n = MatrixType::Random(50,50);
m.conservativeResize(rows,cols);
VERIFY_IS_APPROX(m, n.block(0,0,rows,cols));
}
// random growing with zeroing ...
for (int i=0; i<25; ++i)
{
const int rows = ei_random<int>(50,75);
const int cols = ei_random<int>(50,75);
m = n = MatrixType::Random(50,50);
m.conservativeResizeLike(MatrixType::Zero(rows,cols));
VERIFY_IS_APPROX(m.block(0,0,n.rows(),n.cols()), n);
VERIFY( rows<=50 || m.block(50,0,rows-50,cols).sum() == Scalar(0) );
VERIFY( cols<=50 || m.block(0,50,rows,cols-50).sum() == Scalar(0) );
}
}
template <typename Scalar>
void run_vector_tests()
{
typedef Matrix<Scalar, 1, Eigen::Dynamic> MatrixType;
MatrixType m, n;
// boundary cases ...
m = n = MatrixType::Random(50);
m.conservativeResize(1);
VERIFY_IS_APPROX(m, n.segment(0,1));
m = n = MatrixType::Random(50);
m.conservativeResize(50);
VERIFY_IS_APPROX(m, n.segment(0,50));
// random shrinking ...
for (int i=0; i<50; ++i)
{
const int size = ei_random<int>(1,50);
m = n = MatrixType::Random(50);
m.conservativeResize(size);
VERIFY_IS_APPROX(m, n.segment(0,size));
}
// random growing with zeroing ...
for (int i=0; i<50; ++i)
{
const int size = ei_random<int>(50,100);
m = n = MatrixType::Random(50);
m.conservativeResizeLike(MatrixType::Zero(size));
VERIFY_IS_APPROX(m.segment(0,50), n);
VERIFY( size<=50 || m.segment(50,size-50).sum() == Scalar(0) );
}
}
void test_conservative_resize()
{
run_matrix_tests<int, Eigen::RowMajor>();
run_matrix_tests<int, Eigen::ColMajor>();
run_matrix_tests<float, Eigen::RowMajor>();
run_matrix_tests<float, Eigen::ColMajor>();
run_matrix_tests<double, Eigen::RowMajor>();
run_matrix_tests<double, Eigen::ColMajor>();
run_matrix_tests<std::complex<float>, Eigen::RowMajor>();
run_matrix_tests<std::complex<float>, Eigen::ColMajor>();
run_matrix_tests<std::complex<double>, Eigen::RowMajor>();
run_matrix_tests<std::complex<double>, Eigen::ColMajor>();
run_vector_tests<int>();
run_vector_tests<float>();
run_vector_tests<double>();
run_vector_tests<std::complex<float> >();
run_vector_tests<std::complex<double> >();
}

61
test/eigensolver_complex.cpp Normal file
View File

@ -0,0 +1,61 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra. Eigen itself is part of the KDE project.
//
// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#include "main.h"
#include <Eigen/Eigenvalues>
#include <Eigen/LU>
template<typename MatrixType> void eigensolver(const MatrixType& m)
{
/* this test covers the following files:
ComplexEigenSolver.h, and indirectly ComplexSchur.h
*/
int rows = m.rows();
int cols = m.cols();
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealVectorType;
typedef typename std::complex<typename NumTraits<typename MatrixType::Scalar>::Real> Complex;
MatrixType a = MatrixType::Random(rows,cols);
MatrixType symmA = a.adjoint() * a;
ComplexEigenSolver<MatrixType> ei0(symmA);
VERIFY_IS_APPROX(symmA * ei0.eigenvectors(), ei0.eigenvectors() * ei0.eigenvalues().asDiagonal());
ComplexEigenSolver<MatrixType> ei1(a);
VERIFY_IS_APPROX(a * ei1.eigenvectors(), ei1.eigenvectors() * ei1.eigenvalues().asDiagonal());
}
void test_eigensolver_complex()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( eigensolver(Matrix4cf()) );
CALL_SUBTEST( eigensolver(MatrixXcd(14,14)) );
}
}

2
test/eigensolver_generic.cpp
View File

@ -23,7 +23,7 @@
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#include "main.h"
#include <Eigen/QR>
#include <Eigen/Eigenvalues>
#ifdef HAS_GSL
#include "gsl_helper.h"

2
test/eigensolver_selfadjoint.cpp
View File

@ -23,7 +23,7 @@
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#include "main.h"
#include <Eigen/QR>
#include <Eigen/Eigenvalues>
#ifdef HAS_GSL
#include "gsl_helper.h"

6
test/geo_orthomethods.cpp
View File

@ -86,10 +86,10 @@ template<typename Scalar, int Size> void orthomethods(int size=Size)
VERIFY_IS_MUCH_SMALLER_THAN(v0.unitOrthogonal().dot(v0), Scalar(1));
VERIFY_IS_APPROX(v0.unitOrthogonal().norm(), RealScalar(1));
if (size>3)
if (size>=3)
{
v0.template start<3>().setZero();
v0.end(size-3).setRandom();
v0.template start<2>().setZero();
v0.end(size-2).setRandom();
VERIFY_IS_MUCH_SMALLER_THAN(v0.unitOrthogonal().dot(v0), Scalar(1));
VERIFY_IS_APPROX(v0.unitOrthogonal().norm(), RealScalar(1));

106
test/jacobisvd.cpp Normal file
View File

@ -0,0 +1,106 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#include "main.h"
#include <Eigen/SVD>
#include <Eigen/LU>
template<typename MatrixType, unsigned int Options> void svd(const MatrixType& m = MatrixType(), bool pickrandom = true)
{
int rows = m.rows();
int cols = m.cols();
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime
};
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixUType;
typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime> MatrixVType;
typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType;
typedef Matrix<Scalar, ColsAtCompileTime, 1> InputVectorType;
MatrixType a;
if(pickrandom) a = MatrixType::Random(rows,cols);
else a = m;
JacobiSVD<MatrixType,Options> svd(a);
MatrixType sigma = MatrixType::Zero(rows,cols);
sigma.diagonal() = svd.singularValues().template cast<Scalar>();
MatrixUType u = svd.matrixU();
MatrixVType v = svd.matrixV();
VERIFY_IS_APPROX(a, u * sigma * v.adjoint());
VERIFY_IS_UNITARY(u);
VERIFY_IS_UNITARY(v);
}
template<typename MatrixType> void svd_verify_assert()
{
MatrixType tmp;
SVD<MatrixType> svd;
//VERIFY_RAISES_ASSERT(svd.solve(tmp, &tmp))
VERIFY_RAISES_ASSERT(svd.matrixU())
VERIFY_RAISES_ASSERT(svd.singularValues())
VERIFY_RAISES_ASSERT(svd.matrixV())
/*VERIFY_RAISES_ASSERT(svd.computeUnitaryPositive(&tmp,&tmp))
VERIFY_RAISES_ASSERT(svd.computePositiveUnitary(&tmp,&tmp))
VERIFY_RAISES_ASSERT(svd.computeRotationScaling(&tmp,&tmp))
VERIFY_RAISES_ASSERT(svd.computeScalingRotation(&tmp,&tmp))*/
}
void test_jacobisvd()
{
for(int i = 0; i < g_repeat; i++) {
Matrix2cd m;
m << 0, 1,
0, 1;
CALL_SUBTEST(( svd<Matrix2cd,0>(m, false) ));
m << 1, 0,
1, 0;
CALL_SUBTEST(( svd<Matrix2cd,0>(m, false) ));
Matrix2d n;
n << 1, 1,
1, -1;
CALL_SUBTEST(( svd<Matrix2d,0>(n, false) ));
CALL_SUBTEST(( svd<Matrix3f,0>() ));
CALL_SUBTEST(( svd<Matrix4d,Square>() ));
CALL_SUBTEST(( svd<Matrix<float,3,5> , AtLeastAsManyColsAsRows>() ));
CALL_SUBTEST(( svd<Matrix<double,Dynamic,2> , AtLeastAsManyRowsAsCols>(Matrix<double,Dynamic,2>(10,2)) ));
CALL_SUBTEST(( svd<MatrixXf,Square>(MatrixXf(50,50)) ));
CALL_SUBTEST(( svd<MatrixXcd,AtLeastAsManyRowsAsCols>(MatrixXcd(14,7)) ));
}
CALL_SUBTEST(( svd<MatrixXf,0>(MatrixXf(300,200)) ));
CALL_SUBTEST(( svd<MatrixXcd,AtLeastAsManyColsAsRows>(MatrixXcd(100,150)) ));
CALL_SUBTEST(( svd_verify_assert<Matrix3f>() ));
CALL_SUBTEST(( svd_verify_assert<Matrix3d>() ));
CALL_SUBTEST(( svd_verify_assert<MatrixXf>() ));
CALL_SUBTEST(( svd_verify_assert<MatrixXd>() ));
}

12
test/main.h
View File

@ -30,6 +30,10 @@
#include <vector>
#include <typeinfo>
#ifdef NDEBUG
#undef NDEBUG
#endif
#ifndef EIGEN_TEST_FUNC
#error EIGEN_TEST_FUNC must be defined
#endif
@ -153,6 +157,8 @@ namespace Eigen
#define VERIFY_IS_APPROX_OR_LESS_THAN(a, b) VERIFY(test_ei_isApproxOrLessThan(a, b))
#define VERIFY_IS_NOT_APPROX_OR_LESS_THAN(a, b) VERIFY(!test_ei_isApproxOrLessThan(a, b))
#define VERIFY_IS_UNITARY(a) VERIFY(test_isUnitary(a))
#define CALL_SUBTEST(FUNC) do { \
g_test_stack.push_back(EI_PP_MAKE_STRING(FUNC)); \
FUNC; \
@ -227,6 +233,12 @@ inline bool test_ei_isMuchSmallerThan(const MatrixBase<Derived>& m,
return m.isMuchSmallerThan(s, test_precision<typename ei_traits<Derived>::Scalar>());
}
template<typename Derived>
inline bool test_isUnitary(const MatrixBase<Derived>& m)
{
return m.isUnitary(test_precision<typename ei_traits<Derived>::Scalar>());
}
template<typename MatrixType>
void createRandomMatrixOfRank(int desired_rank, int rows, int cols, MatrixType& m)
{

127
test/mixingtypes.cpp
View File

@ -33,6 +33,7 @@
#include "main.h"
using namespace std;
template<int SizeAtCompileType> void mixingtypes(int size = SizeAtCompileType)
{
@ -45,6 +46,104 @@ template<int SizeAtCompileType> void mixingtypes(int size = SizeAtCompileType)
typedef Matrix<std::complex<float>, SizeAtCompileType, 1> Vec_cf;
typedef Matrix<std::complex<double>, SizeAtCompileType, 1> Vec_cd;
Mat_f mf = Mat_f::Random(size,size);
Mat_d md = mf.template cast<double>();
Mat_cf mcf = Mat_cf::Random(size,size);
Mat_cd mcd = mcf.template cast<complex<double> >();
Vec_f vf = Vec_f::Random(size,1);
Vec_d vd = vf.template cast<double>();
Vec_cf vcf = Vec_cf::Random(size,1);
Vec_cd vcd = vcf.template cast<complex<double> >();
float sf = ei_random<float>();
double sd = ei_random<double>();
complex<float> scf = ei_random<complex<float> >();
complex<double> scd = ei_random<complex<double> >();
mf+mf;
VERIFY_RAISES_ASSERT(mf+md);
VERIFY_RAISES_ASSERT(mf+mcf);
VERIFY_RAISES_ASSERT(vf=vd);
VERIFY_RAISES_ASSERT(vf+=vd);
VERIFY_RAISES_ASSERT(mcd=md);
// check scalar products
VERIFY_IS_APPROX(vcf * sf , vcf * complex<float>(sf));
VERIFY_IS_APPROX(sd * vcd, complex<double>(sd) * vcd);
VERIFY_IS_APPROX(vf * scf , vf.template cast<complex<float> >() * scf);
VERIFY_IS_APPROX(scd * vd, scd * vd.template cast<complex<double> >());
// check dot product
vf.dot(vf);
VERIFY_RAISES_ASSERT(vd.dot(vf));
VERIFY_RAISES_ASSERT(vcf.dot(vf)); // yeah eventually we should allow this but i'm too lazy to make that change now in Dot.h
// especially as that might be rewritten as cwise product .sum() which would make that automatic.
// check diagonal product
VERIFY_IS_APPROX(vf.asDiagonal() * mcf, vf.template cast<complex<float> >().asDiagonal() * mcf);
VERIFY_IS_APPROX(vcd.asDiagonal() * md, vcd.asDiagonal() * md.template cast<complex<double> >());
VERIFY_IS_APPROX(mcf * vf.asDiagonal(), mcf * vf.template cast<complex<float> >().asDiagonal());
VERIFY_IS_APPROX(md * vcd.asDiagonal(), md.template cast<complex<double> >() * vcd.asDiagonal());
// vd.asDiagonal() * mf; // does not even compile
// vcd.asDiagonal() * mf; // does not even compile
// check inner product
VERIFY_IS_APPROX((vf.transpose() * vcf).value(), (vf.template cast<complex<float> >().transpose() * vcf).value());
// check outer product
VERIFY_IS_APPROX((vf * vcf.transpose()).eval(), (vf.template cast<complex<float> >() * vcf.transpose()).eval());
}
void mixingtypes_large(int size)
{
static const int SizeAtCompileType = Dynamic;
typedef Matrix<float, SizeAtCompileType, SizeAtCompileType> Mat_f;
typedef Matrix<double, SizeAtCompileType, SizeAtCompileType> Mat_d;
typedef Matrix<std::complex<float>, SizeAtCompileType, SizeAtCompileType> Mat_cf;
typedef Matrix<std::complex<double>, SizeAtCompileType, SizeAtCompileType> Mat_cd;
typedef Matrix<float, SizeAtCompileType, 1> Vec_f;
typedef Matrix<double, SizeAtCompileType, 1> Vec_d;
typedef Matrix<std::complex<float>, SizeAtCompileType, 1> Vec_cf;
typedef Matrix<std::complex<double>, SizeAtCompileType, 1> Vec_cd;
Mat_f mf(size,size);
Mat_d md(size,size);
Mat_cf mcf(size,size);
Mat_cd mcd(size,size);
Vec_f vf(size,1);
Vec_d vd(size,1);
Vec_cf vcf(size,1);
Vec_cd vcd(size,1);
mf*mf;
// FIXME large products does not allow mixing types
VERIFY_RAISES_ASSERT(md*mcd);
VERIFY_RAISES_ASSERT(mcd*md);
VERIFY_RAISES_ASSERT(mf*vcf);
VERIFY_RAISES_ASSERT(mcf*vf);
VERIFY_RAISES_ASSERT(mcf *= mf);
// VERIFY_RAISES_ASSERT(vcd = md*vcd); // does not even compile (cannot convert complex to double)
VERIFY_RAISES_ASSERT(vcf = mcf*vf);
// VERIFY_RAISES_ASSERT(mf*md); // does not even compile
// VERIFY_RAISES_ASSERT(mcf*mcd); // does not even compile
// VERIFY_RAISES_ASSERT(mcf*vcd); // does not even compile
VERIFY_RAISES_ASSERT(vcf = mf*vf);
}
template<int SizeAtCompileType> void mixingtypes_small()
{
int size = SizeAtCompileType;
typedef Matrix<float, SizeAtCompileType, SizeAtCompileType> Mat_f;
typedef Matrix<double, SizeAtCompileType, SizeAtCompileType> Mat_d;
typedef Matrix<std::complex<float>, SizeAtCompileType, SizeAtCompileType> Mat_cf;
typedef Matrix<std::complex<double>, SizeAtCompileType, SizeAtCompileType> Mat_cd;
typedef Matrix<float, SizeAtCompileType, 1> Vec_f;
typedef Matrix<double, SizeAtCompileType, 1> Vec_d;
typedef Matrix<std::complex<float>, SizeAtCompileType, 1> Vec_cf;
typedef Matrix<std::complex<double>, SizeAtCompileType, 1> Vec_cd;
Mat_f mf(size,size);
Mat_d md(size,size);
Mat_cf mcf(size,size);
@ -54,14 +153,11 @@ template<int SizeAtCompileType> void mixingtypes(int size = SizeAtCompileType)
Vec_cf vcf(size,1);
Vec_cd vcd(size,1);
mf+mf;
VERIFY_RAISES_ASSERT(mf+md);
VERIFY_RAISES_ASSERT(mf+mcf);
VERIFY_RAISES_ASSERT(vf=vd);
VERIFY_RAISES_ASSERT(vf+=vd);
VERIFY_RAISES_ASSERT(mcd=md);
mf*mf;
// FIXME shall we discard those products ?
// 1) currently they work only if SizeAtCompileType is small enough
// 2) in case we vectorize complexes this might be difficult to still allow that
md*mcd;
mcd*md;
mf*vcf;
@ -69,20 +165,19 @@ template<int SizeAtCompileType> void mixingtypes(int size = SizeAtCompileType)
mcf *= mf;
vcd = md*vcd;
vcf = mcf*vf;
VERIFY_RAISES_ASSERT(mf*md);
VERIFY_RAISES_ASSERT(mcf*mcd);
VERIFY_RAISES_ASSERT(mcf*vcd);
// VERIFY_RAISES_ASSERT(mf*md); // does not even compile
// VERIFY_RAISES_ASSERT(mcf*mcd); // does not even compile
// VERIFY_RAISES_ASSERT(mcf*vcd); // does not even compile
VERIFY_RAISES_ASSERT(vcf = mf*vf);
vf.dot(vf);
VERIFY_RAISES_ASSERT(vd.dot(vf));
VERIFY_RAISES_ASSERT(vcf.dot(vf)); // yeah eventually we should allow this but i'm too lazy to make that change now in Dot.h
} // especially as that might be rewritten as cwise product .sum() which would make that automatic.
}
void test_mixingtypes()
{
// check that our operator new is indeed called:
CALL_SUBTEST(mixingtypes<3>(3));
CALL_SUBTEST(mixingtypes<4>(4));
CALL_SUBTEST(mixingtypes<3>());
CALL_SUBTEST(mixingtypes<4>());
CALL_SUBTEST(mixingtypes<Dynamic>(20));
CALL_SUBTEST(mixingtypes_small<4>());
CALL_SUBTEST(mixingtypes_large(20));
}

12
test/product.h
View File

@ -81,7 +81,7 @@ template<typename MatrixType> void product(const MatrixType& m)
m3 = m1;
m3 *= m1.transpose() * m2;
VERIFY_IS_APPROX(m3, m1 * (m1.transpose()*m2));
VERIFY_IS_APPROX(m3, m1.lazy() * (m1.transpose()*m2));
VERIFY_IS_APPROX(m3, m1 * (m1.transpose()*m2));
// continue testing Product.h: distributivity
VERIFY_IS_APPROX(square*(m1 + m2), square*m1+square*m2);
@ -109,26 +109,26 @@ template<typename MatrixType> void product(const MatrixType& m)
// test optimized operator+= path
res = square;
res += (m1 * m2.transpose()).lazy();
res.noalias() += m1 * m2.transpose();
VERIFY_IS_APPROX(res, square + m1 * m2.transpose());
if (NumTraits<Scalar>::HasFloatingPoint && std::min(rows,cols)>1)
{
VERIFY(areNotApprox(res,square + m2 * m1.transpose()));
}
vcres = vc2;
vcres += (m1.transpose() * v1).lazy();
vcres.noalias() += m1.transpose() * v1;
VERIFY_IS_APPROX(vcres, vc2 + m1.transpose() * v1);
// test optimized operator-= path
res = square;
res -= (m1 * m2.transpose()).lazy();
res.noalias() -= m1 * m2.transpose();
VERIFY_IS_APPROX(res, square - (m1 * m2.transpose()));
if (NumTraits<Scalar>::HasFloatingPoint && std::min(rows,cols)>1)
{
VERIFY(areNotApprox(res,square - m2 * m1.transpose()));
}
vcres = vc2;
vcres -= (m1.transpose() * v1).lazy();
vcres.noalias() -= m1.transpose() * v1;
VERIFY_IS_APPROX(vcres, vc2 - m1.transpose() * v1);
tm1 = m1;
@ -145,7 +145,7 @@ template<typename MatrixType> void product(const MatrixType& m)
VERIFY_IS_APPROX(res, m1 * m2.transpose());
res2 = square2;
res2 += (m1.transpose() * m2).lazy();
res2.noalias() += m1.transpose() * m2;
VERIFY_IS_APPROX(res2, square2 + m1.transpose() * m2);
if (NumTraits<Scalar>::HasFloatingPoint && std::min(rows,cols)>1)
{

36
test/product_extra.cpp
View File

@ -59,16 +59,13 @@ template<typename MatrixType> void product_extra(const MatrixType& m)
// r0 = ei_random<int>(0,rows/2-1),
// r1 = ei_random<int>(rows/2,rows);
// all the expressions in this test should be compiled as a single matrix product
// TODO: add internal checks to verify that
VERIFY_IS_APPROX(m3 = (m1 * m2.adjoint()).lazy(), m1 * m2.adjoint().eval());
VERIFY_IS_APPROX(m3 = (m1.adjoint() * square.adjoint()).lazy(), m1.adjoint().eval() * square.adjoint().eval());
VERIFY_IS_APPROX(m3 = (m1.adjoint() * m2).lazy(), m1.adjoint().eval() * m2);
VERIFY_IS_APPROX(m3 = ((s1 * m1.adjoint()) * m2).lazy(), (s1 * m1.adjoint()).eval() * m2);
VERIFY_IS_APPROX(m3 = ((- m1.adjoint() * s1) * (s3 * m2)).lazy(), (- m1.adjoint() * s1).eval() * (s3 * m2).eval());
VERIFY_IS_APPROX(m3 = ((s2 * m1.adjoint() * s1) * m2).lazy(), (s2 * m1.adjoint() * s1).eval() * m2);
VERIFY_IS_APPROX(m3 = ((-m1*s2) * s1*m2.adjoint()).lazy(), (-m1*s2).eval() * (s1*m2.adjoint()).eval());
VERIFY_IS_APPROX(m3.noalias() = m1 * m2.adjoint(), m1 * m2.adjoint().eval());
VERIFY_IS_APPROX(m3.noalias() = m1.adjoint() * square.adjoint(), m1.adjoint().eval() * square.adjoint().eval());
VERIFY_IS_APPROX(m3.noalias() = m1.adjoint() * m2, m1.adjoint().eval() * m2);
VERIFY_IS_APPROX(m3.noalias() = (s1 * m1.adjoint()) * m2, (s1 * m1.adjoint()).eval() * m2);
VERIFY_IS_APPROX(m3.noalias() = (- m1.adjoint() * s1) * (s3 * m2), (- m1.adjoint() * s1).eval() * (s3 * m2).eval());
VERIFY_IS_APPROX(m3.noalias() = (s2 * m1.adjoint() * s1) * m2, (s2 * m1.adjoint() * s1).eval() * m2);
VERIFY_IS_APPROX(m3.noalias() = (-m1*s2) * s1*m2.adjoint(), (-m1*s2).eval() * (s1*m2.adjoint()).eval());
// a very tricky case where a scale factor has to be automatically conjugated:
VERIFY_IS_APPROX( m1.adjoint() * (s1*m2).conjugate(), (m1.adjoint()).eval() * ((s1*m2).conjugate()).eval());
@ -76,7 +73,6 @@ template<typename MatrixType> void product_extra(const MatrixType& m)
// test all possible conjugate combinations for the four matrix-vector product cases:
// std::cerr << "a\n";
VERIFY_IS_APPROX((-m1.conjugate() * s2) * (s1 * vc2),
(-m1.conjugate()*s2).eval() * (s1 * vc2).eval());
VERIFY_IS_APPROX((-m1 * s2) * (s1 * vc2.conjugate()),
@ -84,7 +80,6 @@ template<typename MatrixType> void product_extra(const MatrixType& m)
VERIFY_IS_APPROX((-m1.conjugate() * s2) * (s1 * vc2.conjugate()),
(-m1.conjugate()*s2).eval() * (s1 * vc2.conjugate()).eval());
// std::cerr << "b\n";
VERIFY_IS_APPROX((s1 * vc2.transpose()) * (-m1.adjoint() * s2),
(s1 * vc2.transpose()).eval() * (-m1.adjoint()*s2).eval());
VERIFY_IS_APPROX((s1 * vc2.adjoint()) * (-m1.transpose() * s2),
@ -92,7 +87,6 @@ template<typename MatrixType> void product_extra(const MatrixType& m)
VERIFY_IS_APPROX((s1 * vc2.adjoint()) * (-m1.adjoint() * s2),
(s1 * vc2.adjoint()).eval() * (-m1.adjoint()*s2).eval());
// std::cerr << "c\n";
VERIFY_IS_APPROX((-m1.adjoint() * s2) * (s1 * v1.transpose()),
(-m1.adjoint()*s2).eval() * (s1 * v1.transpose()).eval());
VERIFY_IS_APPROX((-m1.transpose() * s2) * (s1 * v1.adjoint()),
@ -100,7 +94,6 @@ template<typename MatrixType> void product_extra(const MatrixType& m)
VERIFY_IS_APPROX((-m1.adjoint() * s2) * (s1 * v1.adjoint()),
(-m1.adjoint()*s2).eval() * (s1 * v1.adjoint()).eval());
// std::cerr << "d\n";
VERIFY_IS_APPROX((s1 * v1) * (-m1.conjugate() * s2),
(s1 * v1).eval() * (-m1.conjugate()*s2).eval());
VERIFY_IS_APPROX((s1 * v1.conjugate()) * (-m1 * s2),
@ -111,13 +104,24 @@ template<typename MatrixType> void product_extra(const MatrixType& m)
VERIFY_IS_APPROX((-m1.adjoint() * s2) * (s1 * v1.adjoint()),
(-m1.adjoint()*s2).eval() * (s1 * v1.adjoint()).eval());
// test the vector-matrix product with non aligned starts
int i = ei_random<int>(0,m1.rows()-2);
int j = ei_random<int>(0,m1.cols()-2);
int r = ei_random<int>(1,m1.rows()-i);
int c = ei_random<int>(1,m1.cols()-j);
int i2 = ei_random<int>(0,m1.rows()-1);
int j2 = ei_random<int>(0,m1.cols()-1);
VERIFY_IS_APPROX(m1.col(j2).adjoint() * m1.block(0,j,m1.rows(),c), m1.col(j2).adjoint().eval() * m1.block(0,j,m1.rows(),c).eval());
VERIFY_IS_APPROX(m1.block(i,0,r,m1.cols()) * m1.row(i2).adjoint(), m1.block(i,0,r,m1.cols()).eval() * m1.row(i2).adjoint().eval());
}
void test_product_extra()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( product_extra(MatrixXf(ei_random<int>(1,320), ei_random<int>(1,320))) );
CALL_SUBTEST( product_extra(MatrixXf(ei_random<int>(2,320), ei_random<int>(2,320))) );
CALL_SUBTEST( product_extra(MatrixXcf(ei_random<int>(50,50), ei_random<int>(50,50))) );
CALL_SUBTEST( product_extra(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(ei_random<int>(1,50), ei_random<int>(1,50))) );
CALL_SUBTEST( product_extra(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(ei_random<int>(2,50), ei_random<int>(2,50))) );
}
}

44
test/product_notemporary.cpp
View File

@ -69,53 +69,47 @@ template<typename MatrixType> void product_notemporary(const MatrixType& m)
r1 = ei_random<int>(8,rows-r0);
VERIFY_EVALUATION_COUNT( m3 = (m1 * m2.adjoint()), 1);
VERIFY_EVALUATION_COUNT( m3 = (m1 * m2.adjoint()).lazy(), 0);
VERIFY_EVALUATION_COUNT( m3.noalias() = m1 * m2.adjoint(), 0);
// NOTE in this case the slow product is used:
// FIXME:
VERIFY_EVALUATION_COUNT( m3.noalias() = s1 * (m1 * m2.transpose()), 0);
VERIFY_EVALUATION_COUNT( m3 = (s1 * m1 * s2 * m2.adjoint()).lazy(), 0);
VERIFY_EVALUATION_COUNT( m3 = (s1 * m1 * s2 * (m1*s3+m2*s2).adjoint()).lazy(), 1);
VERIFY_EVALUATION_COUNT( m3 = ((s1 * m1).adjoint() * s2 * m2).lazy(), 0);
VERIFY_EVALUATION_COUNT( m3 += (s1 * (-m1*s3).adjoint() * (s2 * m2 * s3)).lazy(), 0);
VERIFY_EVALUATION_COUNT( m3.noalias() = s1 * m1 * s2 * m2.adjoint(), 0);
VERIFY_EVALUATION_COUNT( m3.noalias() = s1 * m1 * s2 * (m1*s3+m2*s2).adjoint(), 1);
VERIFY_EVALUATION_COUNT( m3.noalias() = (s1 * m1).adjoint() * s2 * m2, 0);
VERIFY_EVALUATION_COUNT( m3.noalias() += s1 * (-m1*s3).adjoint() * (s2 * m2 * s3), 0);
VERIFY_EVALUATION_COUNT( m3 -= (s1 * (m1.transpose() * m2)).lazy(), 0);
VERIFY_EVALUATION_COUNT( m3.noalias() -= s1 * (m1.transpose() * m2), 0);
VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1) += (-m1.block(r0,c0,r1,c1) * (s2*m2.block(r0,c0,r1,c1)).adjoint()).lazy() ), 0);
VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1) -= (s1 * m1.block(r0,c0,r1,c1) * m2.block(c0,r0,c1,r1)).lazy() ), 0);
VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1).noalias() += -m1.block(r0,c0,r1,c1) * (s2*m2.block(r0,c0,r1,c1)).adjoint() ), 0);
VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1).noalias() -= s1 * m1.block(r0,c0,r1,c1) * m2.block(c0,r0,c1,r1) ), 0);
// NOTE this is because the Block expression is not handled yet by our expression analyser
VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1) = (s1 * m1.block(r0,c0,r1,c1) * (s1*m2).block(c0,r0,c1,r1)).lazy() ), 1);
VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1).noalias() = s1 * m1.block(r0,c0,r1,c1) * (s1*m2).block(c0,r0,c1,r1) ), 1);
VERIFY_EVALUATION_COUNT( m3 -= ((s1 * m1).template triangularView<LowerTriangular>() * m2).lazy(), 0);
VERIFY_EVALUATION_COUNT( rm3 = ((s1 * m1.adjoint()).template triangularView<UpperTriangular>() * (m2+m2)).lazy(), 1);
VERIFY_EVALUATION_COUNT( rm3 = ((s1 * m1.adjoint()).template triangularView<UnitUpperTriangular>() * m2.adjoint()).lazy(), 0);
VERIFY_EVALUATION_COUNT( m3.noalias() -= (s1 * m1).template triangularView<LowerTriangular>() * m2, 0);
VERIFY_EVALUATION_COUNT( rm3.noalias() = (s1 * m1.adjoint()).template triangularView<UpperTriangular>() * (m2+m2), 1);
VERIFY_EVALUATION_COUNT( rm3.noalias() = (s1 * m1.adjoint()).template triangularView<UnitUpperTriangular>() * m2.adjoint(), 0);
VERIFY_EVALUATION_COUNT( rm3.col(c0) = ((s1 * m1.adjoint()).template triangularView<UnitUpperTriangular>() * (s2*m2.row(c0)).adjoint()).lazy(), 0);
VERIFY_EVALUATION_COUNT( rm3.col(c0).noalias() = (s1 * m1.adjoint()).template triangularView<UnitUpperTriangular>() * (s2*m2.row(c0)).adjoint(), 0);
VERIFY_EVALUATION_COUNT( m1.template triangularView<LowerTriangular>().solveInPlace(m3), 0);
VERIFY_EVALUATION_COUNT( m1.adjoint().template triangularView<LowerTriangular>().solveInPlace(m3.transpose()), 0);
VERIFY_EVALUATION_COUNT( m3 -= ((s1 * m1).adjoint().template selfadjointView<LowerTriangular>() * (-m2*s3).adjoint()).lazy(), 0);
VERIFY_EVALUATION_COUNT( m3 = (s2 * m2.adjoint() * (s1 * m1.adjoint()).template selfadjointView<UpperTriangular>()).lazy(), 0);
VERIFY_EVALUATION_COUNT( rm3 = ((s1 * m1.adjoint()).template selfadjointView<LowerTriangular>() * m2.adjoint()).lazy(), 0);
VERIFY_EVALUATION_COUNT( m3.noalias() -= (s1 * m1).adjoint().template selfadjointView<LowerTriangular>() * (-m2*s3).adjoint(), 0);
VERIFY_EVALUATION_COUNT( m3.noalias() = s2 * m2.adjoint() * (s1 * m1.adjoint()).template selfadjointView<UpperTriangular>(), 0);
VERIFY_EVALUATION_COUNT( rm3.noalias() = (s1 * m1.adjoint()).template selfadjointView<LowerTriangular>() * m2.adjoint(), 0);
VERIFY_EVALUATION_COUNT( m3.col(c0) = ((s1 * m1).adjoint().template selfadjointView<LowerTriangular>() * (-m2.row(c0)*s3).adjoint()).lazy(), 0);
VERIFY_EVALUATION_COUNT( m3.col(c0) -= ((s1 * m1).adjoint().template selfadjointView<UpperTriangular>() * (-m2.row(c0)*s3).adjoint()).lazy(), 0);
VERIFY_EVALUATION_COUNT( m3.col(c0).noalias() = (s1 * m1).adjoint().template selfadjointView<LowerTriangular>() * (-m2.row(c0)*s3).adjoint(), 0);
VERIFY_EVALUATION_COUNT( m3.col(c0).noalias() -= (s1 * m1).adjoint().template selfadjointView<UpperTriangular>() * (-m2.row(c0)*s3).adjoint(), 0);
VERIFY_EVALUATION_COUNT( m3.block(r0,c0,r1,c1) += ((m1.block(r0,r0,r1,r1).template selfadjointView<UpperTriangular>() * (s1*m2.block(r0,c0,r1,c1)) )).lazy(), 0);
VERIFY_EVALUATION_COUNT( m3.block(r0,c0,r1,c1) = ((m1.block(r0,r0,r1,r1).template selfadjointView<UpperTriangular>() * m2.block(r0,c0,r1,c1) )).lazy(), 0);
VERIFY_EVALUATION_COUNT( m3.block(r0,c0,r1,c1).noalias() += m1.block(r0,r0,r1,r1).template selfadjointView<UpperTriangular>() * (s1*m2.block(r0,c0,r1,c1)), 0);
VERIFY_EVALUATION_COUNT( m3.block(r0,c0,r1,c1).noalias() = m1.block(r0,r0,r1,r1).template selfadjointView<UpperTriangular>() * m2.block(r0,c0,r1,c1), 0);
VERIFY_EVALUATION_COUNT( m3.template selfadjointView<LowerTriangular>().rankUpdate(m2.adjoint()), 0);
m3.resize(1,1);
VERIFY_EVALUATION_COUNT( m3 = ((m1.block(r0,r0,r1,r1).template selfadjointView<LowerTriangular>() * m2.block(r0,c0,r1,c1) )).lazy(), 0);
VERIFY_EVALUATION_COUNT( m3.noalias() = m1.block(r0,r0,r1,r1).template selfadjointView<LowerTriangular>() * m2.block(r0,c0,r1,c1), 0);
m3.resize(1,1);
VERIFY_EVALUATION_COUNT( m3 = ((m1.block(r0,r0,r1,r1).template triangularView<UnitUpperTriangular>() * m2.block(r0,c0,r1,c1) )).lazy(), 0);
VERIFY_EVALUATION_COUNT( m3.noalias() = m1.block(r0,r0,r1,r1).template triangularView<UnitUpperTriangular>() * m2.block(r0,c0,r1,c1), 0);
}
void test_product_notemporary()

2
test/product_symm.cpp
View File

@ -96,7 +96,7 @@ template<typename Scalar, int Size, int OtherSize> void symm(int size = Size, in
m2 = m1.template triangularView<UpperTriangular>(); rhs13 = rhs12;
VERIFY_IS_APPROX(rhs12 += (s1 * ((m2.adjoint()).template selfadjointView<LowerTriangular>() * (s2*rhs3).conjugate())).lazy(),
VERIFY_IS_APPROX(rhs12.noalias() += s1 * ((m2.adjoint()).template selfadjointView<LowerTriangular>() * (s2*rhs3).conjugate()),
rhs13 += (s1*m1.adjoint()) * (s2*rhs3).conjugate());
// test matrix * selfadjoint

4
test/product_trsm.cpp
View File

@ -40,8 +40,8 @@ template<typename Scalar> void trsm(int size,int cols)
Matrix<Scalar,Dynamic,Dynamic,ColMajor> cmRhs(size,cols), ref(size,cols);
Matrix<Scalar,Dynamic,Dynamic,RowMajor> rmRhs(size,cols);
cmLhs.setRandom(); cmLhs *= 0.1; cmLhs.diagonal().cwise() += 1;
rmLhs.setRandom(); rmLhs *= 0.1; rmLhs.diagonal().cwise() += 1;
cmLhs.setRandom(); cmLhs *= static_cast<RealScalar>(0.1); cmLhs.diagonal().cwise() += static_cast<RealScalar>(1);
rmLhs.setRandom(); rmLhs *= static_cast<RealScalar>(0.1); rmLhs.diagonal().cwise() += static_cast<RealScalar>(1);
VERIFY_TRSM(cmLhs.conjugate().template triangularView<LowerTriangular>(), cmRhs);
VERIFY_TRSM(cmLhs .template triangularView<UpperTriangular>(), cmRhs);

79
test/stable_norm.cpp Normal file
View File

@ -0,0 +1,79 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#include "main.h"
template<typename MatrixType> void stable_norm(const MatrixType& m)
{
/* this test covers the following files:
StableNorm.h
*/
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
int rows = m.rows();
int cols = m.cols();
Scalar big = ei_random<Scalar>() * std::numeric_limits<RealScalar>::max() * RealScalar(1e-4);
Scalar small = static_cast<RealScalar>(1)/big;
MatrixType vzero = MatrixType::Zero(rows, cols),
vrand = MatrixType::Random(rows, cols),
vbig(rows, cols),
vsmall(rows,cols);
vbig.fill(big);
vsmall.fill(small);
VERIFY_IS_MUCH_SMALLER_THAN(vzero.norm(), static_cast<RealScalar>(1));
VERIFY_IS_APPROX(vrand.stableNorm(), vrand.norm());
VERIFY_IS_APPROX(vrand.blueNorm(), vrand.norm());
VERIFY_IS_APPROX(vrand.hypotNorm(), vrand.norm());
RealScalar size = static_cast<RealScalar>(m.size());
// test overflow
VERIFY_IS_NOT_APPROX(static_cast<Scalar>(vbig.norm()), ei_sqrt(size)*big); // here the default norm must fail
VERIFY_IS_APPROX(static_cast<Scalar>(vbig.stableNorm()), ei_sqrt(size)*big);
VERIFY_IS_APPROX(static_cast<Scalar>(vbig.blueNorm()), ei_sqrt(size)*big);
VERIFY_IS_APPROX(static_cast<Scalar>(vbig.hypotNorm()), ei_sqrt(size)*big);
// test underflow
VERIFY_IS_NOT_APPROX(static_cast<Scalar>(vsmall.norm()), ei_sqrt(size)*small); // here the default norm must fail
VERIFY_IS_APPROX(static_cast<Scalar>(vsmall.stableNorm()), ei_sqrt(size)*small);
VERIFY_IS_APPROX(static_cast<Scalar>(vsmall.blueNorm()), ei_sqrt(size)*small);
VERIFY_IS_APPROX(static_cast<Scalar>(vsmall.hypotNorm()), ei_sqrt(size)*small);
}
void test_stable_norm()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( stable_norm(Matrix<float, 1, 1>()) );
CALL_SUBTEST( stable_norm(Vector4d()) );
CALL_SUBTEST( stable_norm(VectorXd(ei_random<int>(10,2000))) );
CALL_SUBTEST( stable_norm(VectorXf(ei_random<int>(10,2000))) );
CALL_SUBTEST( stable_norm(VectorXcd(ei_random<int>(10,2000))) );
}
}

45
test/submatrices.cpp
View File

@ -170,6 +170,48 @@ template<typename MatrixType> void submatrices(const MatrixType& m)
VERIFY(ei_real(ones.row(r1).dot(ones.row(r2))) == RealScalar(cols));
}
template<typename MatrixType>
void compare_using_data_and_stride(const MatrixType& m)
{
int rows = m.rows();
int cols = m.cols();
int size = m.size();
int stride = m.stride();
const typename MatrixType::Scalar* data = m.data();
for(int j=0;j<cols;++j)
for(int i=0;i<rows;++i)
VERIFY_IS_APPROX(m.coeff(i,j), data[(MatrixType::Flags&RowMajorBit) ? i*stride+j : j*stride + i]);
if(MatrixType::IsVectorAtCompileTime)
{
VERIFY_IS_APPROX(stride, int((&m.coeff(1))-(&m.coeff(0))));
for (int i=0;i<size;++i)
VERIFY_IS_APPROX(m.coeff(i), data[i*stride]);
}
}
template<typename MatrixType>
void data_and_stride(const MatrixType& m)
{
int rows = m.rows();
int cols = m.cols();
int r1 = ei_random<int>(0,rows-1);
int r2 = ei_random<int>(r1,rows-1);
int c1 = ei_random<int>(0,cols-1);
int c2 = ei_random<int>(c1,cols-1);
MatrixType m1 = MatrixType::Random(rows, cols);
compare_using_data_and_stride(m1.block(r1, c1, r2-r1+1, c2-c1+1));
compare_using_data_and_stride(m1.transpose().block(c1, r1, c2-c1+1, r2-r1+1));
compare_using_data_and_stride(m1.row(r1));
compare_using_data_and_stride(m1.col(c1));
compare_using_data_and_stride(m1.row(r1).transpose());
compare_using_data_and_stride(m1.col(c1).transpose());
}
void test_submatrices()
{
for(int i = 0; i < g_repeat; i++) {
@ -179,5 +221,8 @@ void test_submatrices()
CALL_SUBTEST( submatrices(MatrixXi(8, 12)) );
CALL_SUBTEST( submatrices(MatrixXcd(20, 20)) );
CALL_SUBTEST( submatrices(MatrixXf(20, 20)) );
CALL_SUBTEST( data_and_stride(MatrixXf(ei_random(5,50), ei_random(5,50))) );
CALL_SUBTEST( data_and_stride(Matrix<int,Dynamic,Dynamic,RowMajor>(ei_random(5,50), ei_random(5,50))) );
}
}

6
test/svd.cpp
View File

@ -41,15 +41,11 @@ template<typename MatrixType> void svd(const MatrixType& m)
Matrix<Scalar, MatrixType::RowsAtCompileTime, 1>::Random(rows,1);
Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> x(cols,1), x2(cols,1);
RealScalar largerEps = test_precision<RealScalar>();
if (ei_is_same_type<RealScalar,float>::ret)
largerEps = 1e-3f;
{
SVD<MatrixType> svd(a);
MatrixType sigma = MatrixType::Zero(rows,cols);
MatrixType matU = MatrixType::Zero(rows,rows);
sigma.block(0,0,cols,cols) = svd.singularValues().asDiagonal();
sigma.diagonal() = svd.singularValues();
matU = svd.matrixU();
VERIFY_IS_APPROX(a, matU * sigma * svd.matrixV().transpose());
}

98
test/swap.cpp Normal file
View File

@ -0,0 +1,98 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#define EIGEN_NO_STATIC_ASSERT
#include "main.h"
template<typename T>
struct other_matrix_type
{
typedef int type;
};
template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
struct other_matrix_type<Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> >
{
typedef Matrix<_Scalar, _Rows, _Cols, _Options^RowMajor, _MaxRows, _MaxCols> type;
};
template<typename MatrixType> void swap(const MatrixType& m)
{
typedef typename other_matrix_type<MatrixType>::type OtherMatrixType;
typedef typename MatrixType::Scalar Scalar;
ei_assert((!ei_is_same_type<MatrixType,OtherMatrixType>::ret));
int rows = m.rows();
int cols = m.cols();
// construct 3 matrix guaranteed to be distinct
MatrixType m1 = MatrixType::Random(rows,cols);
MatrixType m2 = MatrixType::Random(rows,cols) + Scalar(100) * MatrixType::Identity(rows,cols);
OtherMatrixType m3 = OtherMatrixType::Random(rows,cols) + Scalar(200) * OtherMatrixType::Identity(rows,cols);
MatrixType m1_copy = m1;
MatrixType m2_copy = m2;
OtherMatrixType m3_copy = m3;
// test swapping 2 matrices of same type
m1.swap(m2);
VERIFY_IS_APPROX(m1,m2_copy);
VERIFY_IS_APPROX(m2,m1_copy);
m1 = m1_copy;
m2 = m2_copy;
// test swapping 2 matrices of different types
m1.swap(m3);
VERIFY_IS_APPROX(m1,m3_copy);
VERIFY_IS_APPROX(m3,m1_copy);
m1 = m1_copy;
m3 = m3_copy;
// test swapping matrix with expression
m1.swap(m2.block(0,0,rows,cols));
VERIFY_IS_APPROX(m1,m2_copy);
VERIFY_IS_APPROX(m2,m1_copy);
m1 = m1_copy;
m2 = m2_copy;
// test swapping two expressions of different types
m1.transpose().swap(m3.transpose());
VERIFY_IS_APPROX(m1,m3_copy);
VERIFY_IS_APPROX(m3,m1_copy);
m1 = m1_copy;
m3 = m3_copy;
// test assertion on mismatching size -- matrix case
VERIFY_RAISES_ASSERT(m1.swap(m1.row(0)));
// test assertion on mismatching size -- xpr case
VERIFY_RAISES_ASSERT(m1.row(0).swap(m1));
}
void test_swap()
{
CALL_SUBTEST( swap(Matrix3f()) ); // fixed size, no vectorization
CALL_SUBTEST( swap(Matrix4d()) ); // fixed size, possible vectorization
CALL_SUBTEST( swap(MatrixXd(3,3)) ); // dyn size, no vectorization
CALL_SUBTEST( swap(MatrixXf(30,30)) ); // dyn size, possible vectorization
}

38
test/testsuite.cmake
View File

@ -27,6 +27,13 @@
# - EIGEN_WORK_DIR: directory used to download the source files and make the builds
# default: folder which contains this script
# - EIGEN_CMAKE_ARGS: additional arguments passed to cmake
# - EIGEN_GENERATOR_TYPE: allows to overwrite the generator type
# default: nmake (windows
# See http://www.cmake.org/cmake/help/cmake2.6docs.html#section_Generators for a complete
# list of supported generators.
# - EIGEN_NO_UPDATE: allows to submit dash boards from local repositories
# This might be interesting in case you want to submit dashboards
# including local changes.
# - CTEST_SOURCE_DIRECTORY: path to eigen's src (use a new and empty folder, not the one you are working on)
# default: <EIGEN_WORK_DIR>/src
# - CTEST_BINARY_DIRECTORY: build directory
@ -132,11 +139,11 @@ endif(NOT EIGEN_MODE)
## mandatory variables (the default should be ok in most cases):
if(NOT IGNORE_CVS)
if(NOT EIGEN_NO_UPDATE)
SET (CTEST_CVS_COMMAND "hg")
SET (CTEST_CVS_CHECKOUT "${CTEST_CVS_COMMAND} clone http://bitbucket.org/eigen/eigen2 \"${CTEST_SOURCE_DIRECTORY}\"")
SET(CTEST_BACKUP_AND_RESTORE TRUE) # the backup is CVS related ...
endif(NOT IGNORE_CVS)
endif(NOT EIGEN_NO_UPDATE)
# which ctest command to use for running the dashboard
SET (CTEST_COMMAND "${EIGEN_CMAKE_DIR}ctest -D ${EIGEN_MODE}")
@ -158,15 +165,24 @@ SET(CTEST_START_WITH_EMPTY_BINARY_DIRECTORY TRUE)
# any quotes inside of this string if you use it
if(WIN32 AND NOT UNIX)
#message(SEND_ERROR "win32")
set(CTEST_CMAKE_COMMAND "${CTEST_CMAKE_COMMAND} -G \"NMake Makefiles\" -DCMAKE_MAKE_PROGRAM=nmake")
SET (CTEST_INITIAL_CACHE "
MAKECOMMAND:STRING=nmake -i
CMAKE_MAKE_PROGRAM:FILEPATH=nmake
CMAKE_GENERATOR:INTERNAL=NMake Makefiles
CMAKE_BUILD_TYPE:STRING=Release
BUILDNAME:STRING=${EIGEN_BUILD_STRING}
SITE:STRING=${EIGEN_SITE}
")
if(EIGEN_GENERATOR_TYPE)
set(CTEST_CMAKE_COMMAND "${CTEST_CMAKE_COMMAND} -G \"${EIGEN_GENERATOR_TYPE}\"")
SET (CTEST_INITIAL_CACHE "
CMAKE_BUILD_TYPE:STRING=Release
BUILDNAME:STRING=${EIGEN_BUILD_STRING}
SITE:STRING=${EIGEN_SITE}
")
else(EIGEN_GENERATOR_TYPE)
set(CTEST_CMAKE_COMMAND "${CTEST_CMAKE_COMMAND} -G \"NMake Makefiles\" -DCMAKE_MAKE_PROGRAM=nmake")
SET (CTEST_INITIAL_CACHE "
MAKECOMMAND:STRING=nmake -i
CMAKE_MAKE_PROGRAM:FILEPATH=nmake
CMAKE_GENERATOR:INTERNAL=NMake Makefiles
CMAKE_BUILD_TYPE:STRING=Release
BUILDNAME:STRING=${EIGEN_BUILD_STRING}
SITE:STRING=${EIGEN_SITE}
")
endif(EIGEN_GENERATOR_TYPE)
else(WIN32 AND NOT UNIX)
SET (CTEST_INITIAL_CACHE "
BUILDNAME:STRING=${EIGEN_BUILD_STRING}

4
test/umeyama.cpp
View File

@ -127,7 +127,7 @@ void run_test(int dim, int num_elements)
MatrixX src = MatrixX::Random(dim+1, num_elements);
src.row(dim) = Matrix<Scalar, 1, Dynamic>::Constant(num_elements, Scalar(1));
MatrixX dst = (cR_t*src).lazy();
MatrixX dst = cR_t*src;
MatrixX cR_t_umeyama = umeyama(src.block(0,0,dim,num_elements), dst.block(0,0,dim,num_elements));
@ -160,7 +160,7 @@ void run_fixed_size_test(int num_elements)
MatrixX src = MatrixX::Random(dim+1, num_elements);
src.row(dim) = Matrix<Scalar, 1, Dynamic>::Constant(num_elements, Scalar(1));
MatrixX dst = (cR_t*src).lazy();
MatrixX dst = cR_t*src;
Block<MatrixX, Dimension, Dynamic> src_block(src,0,0,dim,num_elements);
Block<MatrixX, Dimension, Dynamic> dst_block(dst,0,0,dim,num_elements);

2
unsupported/Eigen/BVH
View File

@ -94,8 +94,6 @@ namespace Eigen {
* responsibility of the intersectObject function to keep track of the results in whatever manner is appropriate.
* The cartesian product intersection and the BVMinimize queries are similar--see their individual documentation.
*
* \addexample BVH_Example \label How to use a BVH to find the closest pair between two point sets
*
* The following is a simple but complete example for how to use the BVH to accelerate the search for a closest red-blue point pair:
* \include BVH_Example.cpp
* Output: \verbinclude BVH_Example.out

464
unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h
View File

@ -25,6 +25,10 @@
#ifndef EIGEN_MATRIX_EXPONENTIAL
#define EIGEN_MATRIX_EXPONENTIAL
#ifdef _MSC_VER
template <typename Scalar> Scalar log2(Scalar v) { return std::log(v)/std::log(Scalar(2)); }
#endif
/** \brief Compute the matrix exponential.
*
* \param M matrix whose exponential is to be computed.
@ -61,257 +65,243 @@ template <typename Derived>
EIGEN_STRONG_INLINE void ei_matrix_exponential(const MatrixBase<Derived> &M,
typename MatrixBase<Derived>::PlainMatrixType* result);
/** \brief Class for computing the matrix exponential.*/
template <typename MatrixType>
class MatrixExponential {
/** \internal \brief Internal helper functions for computing the
* matrix exponential.
*/
namespace MatrixExponentialInternal {
public:
/** \brief Compute the matrix exponential.
*
* \param M matrix whose exponential is to be computed.
* \param result pointer to the matrix in which to store the result.
*/
MatrixExponential(const MatrixType &M, MatrixType *result);
#ifdef _MSC_VER
template <typename Scalar> Scalar log2(Scalar v) { return std::log(v)/std::log(Scalar(2)); }
#endif
private:
// Prevent copying
MatrixExponential(const MatrixExponential&);
MatrixExponential& operator=(const MatrixExponential&);
/** \brief Compute the (3,3)-Pad&eacute; approximant to the exponential.
*
* After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
* approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
*
* \param A Argument of matrix exponential
*/
void pade3(const MatrixType &A);
/** \brief Compute the (5,5)-Pad&eacute; approximant to the exponential.
*
* After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
* approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
*
* \param A Argument of matrix exponential
*/
void pade5(const MatrixType &A);
/** \brief Compute the (7,7)-Pad&eacute; approximant to the exponential.
*
* After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
* approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
*
* \param A Argument of matrix exponential
*/
void pade7(const MatrixType &A);
/** \brief Compute the (9,9)-Pad&eacute; approximant to the exponential.
*
* After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
* approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
*
* \param A Argument of matrix exponential
*/
void pade9(const MatrixType &A);
/** \brief Compute the (13,13)-Pad&eacute; approximant to the exponential.
*
* After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
* approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
*
* \param A Argument of matrix exponential
*/
void pade13(const MatrixType &A);
/** \brief Compute Pad&eacute; approximant to the exponential.
*
* Computes \c m_U, \c m_V and \c m_squarings such that
* \f$ (V+U)(V-U)^{-1} \f$ is a Pad&eacute; of
* \f$ \exp(2^{-\mbox{squarings}}M) \f$ around \f$ M = 0 \f$. The
* degree of the Pad&eacute; approximant and the value of
* squarings are chosen such that the approximation error is no
* more than the round-off error.
*
* The argument of this function should correspond with the (real
* part of) the entries of \c m_M. It is used to select the
* correct implementation using overloading.
*/
void computeUV(double);
/** \brief Compute Pad&eacute; approximant to the exponential.
*
* \sa computeUV(double);
*/
void computeUV(float);
/** \internal \brief Compute the (3,3)-Pad&eacute; approximant to
* the exponential.
*
* After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
* approximant of \f$ \exp(M) \f$ around \f$ M = 0 \f$.
*
* \param M Argument of matrix exponential
* \param Id Identity matrix of same size as M
* \param tmp Temporary storage, to be provided by the caller
* \param M2 Temporary storage, to be provided by the caller
* \param U Even-degree terms in numerator of Pad&eacute; approximant
* \param V Odd-degree terms in numerator of Pad&eacute; approximant
*/
template <typename MatrixType>
EIGEN_STRONG_INLINE void pade3(const MatrixType &M, const MatrixType& Id, MatrixType& tmp,
MatrixType& M2, MatrixType& U, MatrixType& V)
{
typedef typename ei_traits<MatrixType>::Scalar Scalar;
const Scalar b[] = {120., 60., 12., 1.};
M2 = (M * M).lazy();
tmp = b[3]*M2 + b[1]*Id;
U = (M * tmp).lazy();
V = b[2]*M2 + b[0]*Id;
}
/** \internal \brief Compute the (5,5)-Pad&eacute; approximant to
* the exponential.
*
* After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
* approximant of \f$ \exp(M) \f$ around \f$ M = 0 \f$.
*
* \param M Argument of matrix exponential
* \param Id Identity matrix of same size as M
* \param tmp Temporary storage, to be provided by the caller
* \param M2 Temporary storage, to be provided by the caller
* \param U Even-degree terms in numerator of Pad&eacute; approximant
* \param V Odd-degree terms in numerator of Pad&eacute; approximant
*/
template <typename MatrixType>
EIGEN_STRONG_INLINE void pade5(const MatrixType &M, const MatrixType& Id, MatrixType& tmp,
MatrixType& M2, MatrixType& U, MatrixType& V)
{
typedef typename ei_traits<MatrixType>::Scalar Scalar;
const Scalar b[] = {30240., 15120., 3360., 420., 30., 1.};
M2 = (M * M).lazy();
MatrixType M4 = (M2 * M2).lazy();
tmp = b[5]*M4 + b[3]*M2 + b[1]*Id;
U = (M * tmp).lazy();
V = b[4]*M4 + b[2]*M2 + b[0]*Id;
}
/** \internal \brief Compute the (7,7)-Pad&eacute; approximant to
* the exponential.
*
* After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
* approximant of \f$ \exp(M) \f$ around \f$ M = 0 \f$.
*
* \param M Argument of matrix exponential
* \param Id Identity matrix of same size as M
* \param tmp Temporary storage, to be provided by the caller
* \param M2 Temporary storage, to be provided by the caller
* \param U Even-degree terms in numerator of Pad&eacute; approximant
* \param V Odd-degree terms in numerator of Pad&eacute; approximant
*/
template <typename MatrixType>
EIGEN_STRONG_INLINE void pade7(const MatrixType &M, const MatrixType& Id, MatrixType& tmp,
MatrixType& M2, MatrixType& U, MatrixType& V)
{
typedef typename ei_traits<MatrixType>::Scalar Scalar;
const Scalar b[] = {17297280., 8648640., 1995840., 277200., 25200., 1512., 56., 1.};
M2 = (M * M).lazy();
MatrixType M4 = (M2 * M2).lazy();
MatrixType M6 = (M4 * M2).lazy();
tmp = b[7]*M6 + b[5]*M4 + b[3]*M2 + b[1]*Id;
U = (M * tmp).lazy();
V = b[6]*M6 + b[4]*M4 + b[2]*M2 + b[0]*Id;
}
/** \internal \brief Compute the (9,9)-Pad&eacute; approximant to
* the exponential.
*
* After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
* approximant of \f$ \exp(M) \f$ around \f$ M = 0 \f$.
*
* \param M Argument of matrix exponential
* \param Id Identity matrix of same size as M
* \param tmp Temporary storage, to be provided by the caller
* \param M2 Temporary storage, to be provided by the caller
* \param U Even-degree terms in numerator of Pad&eacute; approximant
* \param V Odd-degree terms in numerator of Pad&eacute; approximant
*/
template <typename MatrixType>
EIGEN_STRONG_INLINE void pade9(const MatrixType &M, const MatrixType& Id, MatrixType& tmp,
MatrixType& M2, MatrixType& U, MatrixType& V)
{
typedef typename ei_traits<MatrixType>::Scalar Scalar;
const Scalar b[] = {17643225600., 8821612800., 2075673600., 302702400., 30270240.,
typedef typename NumTraits<typename ei_traits<MatrixType>::Scalar>::Real RealScalar;
/** \brief Pointer to matrix whose exponential is to be computed. */
const MatrixType* m_M;
/** \brief Even-degree terms in numerator of Pad&eacute; approximant. */
MatrixType m_U;
/** \brief Odd-degree terms in numerator of Pad&eacute; approximant. */
MatrixType m_V;
/** \brief Used for temporary storage. */
MatrixType m_tmp1;
/** \brief Used for temporary storage. */
MatrixType m_tmp2;
/** \brief Identity matrix of the same size as \c m_M. */
MatrixType m_Id;
/** \brief Number of squarings required in the last step. */
int m_squarings;
/** \brief L1 norm of m_M. */
float m_l1norm;
};
template <typename MatrixType>
MatrixExponential<MatrixType>::MatrixExponential(const MatrixType &M, MatrixType *result) :
m_M(&M),
m_U(M.rows(),M.cols()),
m_V(M.rows(),M.cols()),
m_tmp1(M.rows(),M.cols()),
m_tmp2(M.rows(),M.cols()),
m_Id(MatrixType::Identity(M.rows(), M.cols())),
m_squarings(0),
m_l1norm(static_cast<float>(M.cwise().abs().colwise().sum().maxCoeff()))
{
computeUV(RealScalar());
m_tmp1 = m_U + m_V; // numerator of Pade approximant
m_tmp2 = -m_U + m_V; // denominator of Pade approximant
m_tmp2.partialLu().solve(m_tmp1, result);
for (int i=0; i<m_squarings; i++)
*result *= *result; // undo scaling by repeated squaring
}
template <typename MatrixType>
EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade3(const MatrixType &A)
{
const Scalar b[] = {120., 60., 12., 1.};
m_tmp1.noalias() = A * A;
m_tmp2 = b[3]*m_tmp1 + b[1]*m_Id;
m_U.noalias() = A * m_tmp2;
m_V = b[2]*m_tmp1 + b[0]*m_Id;
}
template <typename MatrixType>
EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade5(const MatrixType &A)
{
const Scalar b[] = {30240., 15120., 3360., 420., 30., 1.};
MatrixType A2 = A * A;
m_tmp1.noalias() = A2 * A2;
m_tmp2 = b[5]*m_tmp1 + b[3]*A2 + b[1]*m_Id;
m_U.noalias() = A * m_tmp2;
m_V = b[4]*m_tmp1 + b[2]*A2 + b[0]*m_Id;
}
template <typename MatrixType>
EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade7(const MatrixType &A)
{
const Scalar b[] = {17297280., 8648640., 1995840., 277200., 25200., 1512., 56., 1.};
MatrixType A2 = A * A;
MatrixType A4 = A2 * A2;
m_tmp1.noalias() = A4 * A2;
m_tmp2 = b[7]*m_tmp1 + b[5]*A4 + b[3]*A2 + b[1]*m_Id;
m_U.noalias() = A * m_tmp2;
m_V = b[6]*m_tmp1 + b[4]*A4 + b[2]*A2 + b[0]*m_Id;
}
template <typename MatrixType>
EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade9(const MatrixType &A)
{
const Scalar b[] = {17643225600., 8821612800., 2075673600., 302702400., 30270240.,
2162160., 110880., 3960., 90., 1.};
M2 = (M * M).lazy();
MatrixType M4 = (M2 * M2).lazy();
MatrixType M6 = (M4 * M2).lazy();
MatrixType M8 = (M6 * M2).lazy();
tmp = b[9]*M8 + b[7]*M6 + b[5]*M4 + b[3]*M2 + b[1]*Id;
U = (M * tmp).lazy();
V = b[8]*M8 + b[6]*M6 + b[4]*M4 + b[2]*M2 + b[0]*Id;
}
/** \internal \brief Compute the (13,13)-Pad&eacute; approximant to
* the exponential.
*
* After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
* approximant of \f$ \exp(M) \f$ around \f$ M = 0 \f$.
*
* \param M Argument of matrix exponential
* \param Id Identity matrix of same size as M
* \param tmp Temporary storage, to be provided by the caller
* \param M2 Temporary storage, to be provided by the caller
* \param U Even-degree terms in numerator of Pad&eacute; approximant
* \param V Odd-degree terms in numerator of Pad&eacute; approximant
*/
template <typename MatrixType>
EIGEN_STRONG_INLINE void pade13(const MatrixType &M, const MatrixType& Id, MatrixType& tmp,
MatrixType& M2, MatrixType& U, MatrixType& V)
{
typedef typename ei_traits<MatrixType>::Scalar Scalar;
const Scalar b[] = {64764752532480000., 32382376266240000., 7771770303897600.,
MatrixType A2 = A * A;
MatrixType A4 = A2 * A2;
MatrixType A6 = A4 * A2;
m_tmp1.noalias() = A6 * A2;
m_tmp2 = b[9]*m_tmp1 + b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*m_Id;
m_U.noalias() = A * m_tmp2;
m_V = b[8]*m_tmp1 + b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*m_Id;
}
template <typename MatrixType>
EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade13(const MatrixType &A)
{
const Scalar b[] = {64764752532480000., 32382376266240000., 7771770303897600.,
1187353796428800., 129060195264000., 10559470521600., 670442572800.,
33522128640., 1323241920., 40840800., 960960., 16380., 182., 1.};
M2 = (M * M).lazy();
MatrixType M4 = (M2 * M2).lazy();
MatrixType M6 = (M4 * M2).lazy();
V = b[13]*M6 + b[11]*M4 + b[9]*M2;
tmp = (M6 * V).lazy();
tmp += b[7]*M6 + b[5]*M4 + b[3]*M2 + b[1]*Id;
U = (M * tmp).lazy();
tmp = b[12]*M6 + b[10]*M4 + b[8]*M2;
V = (M6 * tmp).lazy();
V += b[6]*M6 + b[4]*M4 + b[2]*M2 + b[0]*Id;
}
/** \internal \brief Helper class for computing Pad&eacute;
* approximants to the exponential.
*/
template <typename MatrixType, typename RealScalar = typename NumTraits<typename ei_traits<MatrixType>::Scalar>::Real>
struct computeUV_selector
{
/** \internal \brief Compute Pad&eacute; approximant to the exponential.
*
* Computes \p U, \p V and \p squarings such that \f$ (V+U)(V-U)^{-1} \f$
* is a Pad&eacute; of \f$ \exp(2^{-\mbox{squarings}}M) \f$
* around \f$ M = 0 \f$. The degree of the Pad&eacute;
* approximant and the value of squarings are chosen such that
* the approximation error is no more than the round-off error.
*
* \param M Argument of matrix exponential
* \param Id Identity matrix of same size as M
* \param tmp1 Temporary storage, to be provided by the caller
* \param tmp2 Temporary storage, to be provided by the caller
* \param U Even-degree terms in numerator of Pad&eacute; approximant
* \param V Odd-degree terms in numerator of Pad&eacute; approximant
* \param l1norm L<sub>1</sub> norm of M
* \param squarings Pointer to integer containing number of times
* that the result needs to be squared to find the
* matrix exponential
*/
static void run(const MatrixType &M, const MatrixType& Id, MatrixType& tmp1, MatrixType& tmp2,
MatrixType& U, MatrixType& V, float l1norm, int* squarings);
};
template <typename MatrixType>
struct computeUV_selector<MatrixType, float>
{
static void run(const MatrixType &M, const MatrixType& Id, MatrixType& tmp1, MatrixType& tmp2,
MatrixType& U, MatrixType& V, float l1norm, int* squarings)
{
*squarings = 0;
if (l1norm < 4.258730016922831e-001) {
pade3(M, Id, tmp1, tmp2, U, V);
} else if (l1norm < 1.880152677804762e+000) {
pade5(M, Id, tmp1, tmp2, U, V);
} else {
const float maxnorm = 3.925724783138660;
*squarings = std::max(0, (int)ceil(log2(l1norm / maxnorm)));
MatrixType A = M / std::pow(typename ei_traits<MatrixType>::Scalar(2), *squarings);
pade7(A, Id, tmp1, tmp2, U, V);
}
}
};
template <typename MatrixType>
struct computeUV_selector<MatrixType, double>
{
static void run(const MatrixType &M, const MatrixType& Id, MatrixType& tmp1, MatrixType& tmp2,
MatrixType& U, MatrixType& V, float l1norm, int* squarings)
{
*squarings = 0;
if (l1norm < 1.495585217958292e-002) {
pade3(M, Id, tmp1, tmp2, U, V);
} else if (l1norm < 2.539398330063230e-001) {
pade5(M, Id, tmp1, tmp2, U, V);
} else if (l1norm < 9.504178996162932e-001) {
pade7(M, Id, tmp1, tmp2, U, V);
} else if (l1norm < 2.097847961257068e+000) {
pade9(M, Id, tmp1, tmp2, U, V);
} else {
const double maxnorm = 5.371920351148152;
*squarings = std::max(0, (int)ceil(log2(l1norm / maxnorm)));
MatrixType A = M / std::pow(typename ei_traits<MatrixType>::Scalar(2), *squarings);
pade13(A, Id, tmp1, tmp2, U, V);
}
}
};
/** \internal \brief Compute the matrix exponential.
*
* \param M matrix whose exponential is to be computed.
* \param result pointer to the matrix in which to store the result.
*/
template <typename MatrixType>
void compute(const MatrixType &M, MatrixType* result)
{
MatrixType num, den, U, V;
MatrixType Id = MatrixType::Identity(M.rows(), M.cols());
float l1norm = M.cwise().abs().colwise().sum().maxCoeff();
int squarings;
computeUV_selector<MatrixType>::run(M, Id, num, den, U, V, l1norm, &squarings);
num = U + V; // numerator of Pade approximant
den = -U + V; // denominator of Pade approximant
den.partialLu().solve(num, result);
for (int i=0; i<squarings; i++)
*result *= *result; // undo scaling by repeated squaring
}
MatrixType A2 = A * A;
MatrixType A4 = A2 * A2;
m_tmp1.noalias() = A4 * A2;
m_V = b[13]*m_tmp1 + b[11]*A4 + b[9]*A2; // used for temporary storage
m_tmp2.noalias() = m_tmp1 * m_V;
m_tmp2 += b[7]*m_tmp1 + b[5]*A4 + b[3]*A2 + b[1]*m_Id;
m_U.noalias() = A * m_tmp2;
m_tmp2 = b[12]*m_tmp1 + b[10]*A4 + b[8]*A2;
m_V.noalias() = m_tmp1 * m_tmp2;
m_V += b[6]*m_tmp1 + b[4]*A4 + b[2]*A2 + b[0]*m_Id;
}
} // end of namespace MatrixExponentialInternal
template <typename MatrixType>
void MatrixExponential<MatrixType>::computeUV(float)
{
if (m_l1norm < 4.258730016922831e-001) {
pade3(*m_M);
} else if (m_l1norm < 1.880152677804762e+000) {
pade5(*m_M);
} else {
const float maxnorm = 3.925724783138660f;
m_squarings = std::max(0, (int)ceil(log2(m_l1norm / maxnorm)));
MatrixType A = *m_M / std::pow(Scalar(2), m_squarings);
pade7(A);
}
}
template <typename MatrixType>
void MatrixExponential<MatrixType>::computeUV(double)
{
if (m_l1norm < 1.495585217958292e-002) {
pade3(*m_M);
} else if (m_l1norm < 2.539398330063230e-001) {
pade5(*m_M);
} else if (m_l1norm < 9.504178996162932e-001) {
pade7(*m_M);
} else if (m_l1norm < 2.097847961257068e+000) {
pade9(*m_M);
} else {
const double maxnorm = 5.371920351148152;
m_squarings = std::max(0, (int)ceil(log2(m_l1norm / maxnorm)));
MatrixType A = *m_M / std::pow(Scalar(2), m_squarings);
pade13(A);
}
}
template <typename Derived>
EIGEN_STRONG_INLINE void ei_matrix_exponential(const MatrixBase<Derived> &M,
typename MatrixBase<Derived>::PlainMatrixType* result)
{
ei_assert(M.rows() == M.cols());
MatrixExponentialInternal::compute(M.eval(), result);
MatrixExponential<typename MatrixBase<Derived>::PlainMatrixType>(M, result);
}
#endif // EIGEN_MATRIX_EXPONENTIAL

1
unsupported/doc/Doxyfile.in
View File

@ -211,7 +211,6 @@ ALIASES = "only_for_vectors=This is only for vectors (either row-
"svd_module=This is defined in the %SVD module. \code #include <Eigen/SVD> \endcode" \
"geometry_module=This is defined in the %Geometry module. \code #include <Eigen/Geometry> \endcode" \
"leastsquares_module=This is defined in the %LeastSquares module. \code #include <Eigen/LeastSquares> \endcode" \
"addexample=\anchor" \
"label=\bug" \
"redstar=<a href='#warningarraymodule' style='color:red;text-decoration: none;'>*</a>" \
"nonstableyet=\warning This is not considered to be part of the stable public API yet. Changes may happen in future releases. See \ref Experimental \"Experimental parts of Eigen\""

18
unsupported/test/matrixExponential.cpp
View File

@ -43,10 +43,10 @@ void test2dRotation(double tol)
A << 0, 1, -1, 0;
for (int i=0; i<=20; i++)
{
angle = pow(10, i / 5. - 2);
angle = static_cast<T>(pow(10, i / 5. - 2));
B << cos(angle), sin(angle), -sin(angle), cos(angle);
ei_matrix_exponential(angle*A, &C);
VERIFY(C.isApprox(B, tol));
VERIFY(C.isApprox(B, static_cast<T>(tol)));
}
}
@ -59,13 +59,13 @@ void test2dHyperbolicRotation(double tol)
for (int i=0; i<=20; i++)
{
angle = (i-10) / 2.0;
angle = static_cast<T>((i-10) / 2.0);
ch = std::cosh(angle);
sh = std::sinh(angle);
A << 0, angle*imagUnit, -angle*imagUnit, 0;
B << ch, sh*imagUnit, -sh*imagUnit, ch;
ei_matrix_exponential(A, &C);
VERIFY(C.isApprox(B, tol));
VERIFY(C.isApprox(B, static_cast<T>(tol)));
}
}
@ -77,13 +77,13 @@ void testPascal(double tol)
Matrix<T,Dynamic,Dynamic> A(size,size), B(size,size), C(size,size);
A.setZero();
for (int i=0; i<size-1; i++)
A(i+1,i) = i+1;
A(i+1,i) = static_cast<T>(i+1);
B.setZero();
for (int i=0; i<size; i++)
for (int j=0; j<=i; j++)
B(i,j) = binom(i,j);
B(i,j) = static_cast<T>(binom(i,j));
ei_matrix_exponential(A, &C);
VERIFY(C.isApprox(B, tol));
VERIFY(C.isApprox(B, static_cast<T>(tol)));
}
}
@ -98,11 +98,13 @@ void randomTest(const MatrixType& m, double tol)
MatrixType m1(rows, cols), m2(rows, cols), m3(rows, cols),
identity = MatrixType::Identity(rows, rows);
typedef typename NumTraits<typename ei_traits<MatrixType>::Scalar>::Real RealScalar;
for(int i = 0; i < g_repeat; i++) {
m1 = MatrixType::Random(rows, cols);
ei_matrix_exponential(m1, &m2);
ei_matrix_exponential(-m1, &m3);
VERIFY(identity.isApprox(m2 * m3, tol));
VERIFY(identity.isApprox(m2 * m3, static_cast<RealScalar>(tol)));
}
}