This commit is contained in:
Gael Guennebaud 2009-05-20 00:05:28 +02:00
commit 6ecd02d7ec
14 changed files with 160 additions and 77 deletions

View File

@ -37,6 +37,9 @@ if(NOT WIN32)
option(EIGEN_BUILD_LIB "Build the binary shared library" OFF) option(EIGEN_BUILD_LIB "Build the binary shared library" OFF)
endif(NOT WIN32) endif(NOT WIN32)
option(EIGEN_BUILD_BTL "Build benchmark suite" OFF) option(EIGEN_BUILD_BTL "Build benchmark suite" OFF)
if(NOT WIN32)
option(EIGEN_BUILD_PKGCONFIG "Build pkg-config .pc file for Eigen" ON)
endif(NOT WIN32)
if(EIGEN_BUILD_LIB) if(EIGEN_BUILD_LIB)
option(EIGEN_TEST_LIB "Build the unit tests using the library (disable -pedantic)" OFF) option(EIGEN_TEST_LIB "Build the unit tests using the library (disable -pedantic)" OFF)
@ -108,6 +111,13 @@ set(INCLUDE_INSTALL_DIR
"The directory where we install the header files" "The directory where we install the header files"
FORCE) FORCE)
if(EIGEN_BUILD_PKGCONFIG)
configure_file(eigen2.pc.in eigen2.pc)
install(FILES eigen2.pc
DESTINATION lib/pkgconfig
)
endif(EIGEN_BUILD_PKGCONFIG)
add_subdirectory(Eigen) add_subdirectory(Eigen)
add_subdirectory(doc) add_subdirectory(doc)

View File

@ -3,11 +3,11 @@
## project to incorporate the testing dashboard. ## project to incorporate the testing dashboard.
## # The following are required to uses Dart and the Cdash dashboard ## # The following are required to uses Dart and the Cdash dashboard
## ENABLE_TESTING() ## ENABLE_TESTING()
## INCLUDE(Dart) ## INCLUDE(CTest)
set(CTEST_PROJECT_NAME "Eigen") set(CTEST_PROJECT_NAME "Eigen")
set(CTEST_NIGHTLY_START_TIME "05:00:00 UTC") set(CTEST_NIGHTLY_START_TIME "05:00:00 UTC")
set(CTEST_DROP_METHOD "http") set(CTEST_DROP_METHOD "http")
set(CTEST_DROP_SITE "www.cdash.org") set(CTEST_DROP_SITE "my.cdash.org")
set(CTEST_DROP_LOCATION "/CDashPublic/submit.php?project=Eigen") set(CTEST_DROP_LOCATION "/submit.php?project=Eigen")
set(CTEST_DROP_SITE_CDASH TRUE) set(CTEST_DROP_SITE_CDASH TRUE)

View File

@ -225,8 +225,8 @@ void PartialLU<MatrixType>::solve(
/* The decomposition PA = LU can be rewritten as A = P^{-1} L U. /* The decomposition PA = LU can be rewritten as A = P^{-1} L U.
* So we proceed as follows: * So we proceed as follows:
* Step 1: compute c = Pb. * Step 1: compute c = Pb.
* Step 2: replace c by the solution x to Lx = c. Exists because L is invertible. * Step 2: replace c by the solution x to Lx = c.
* Step 3: replace c by the solution x to Ux = c. Check if a solution really exists. * Step 3: replace c by the solution x to Ux = c.
*/ */
const int size = m_lu.rows(); const int size = m_lu.rows();

View File

@ -67,10 +67,10 @@ class DynamicSparseMatrix
// EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(DynamicSparseMatrix, +=) // EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(DynamicSparseMatrix, +=)
// EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(DynamicSparseMatrix, -=) // EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(DynamicSparseMatrix, -=)
typedef MappedSparseMatrix<Scalar,Flags> Map; typedef MappedSparseMatrix<Scalar,Flags> Map;
using Base::IsRowMajor;
protected: protected:
enum { IsRowMajor = Base::IsRowMajor };
typedef DynamicSparseMatrix<Scalar,(Flags&~RowMajorBit)|(IsRowMajor?RowMajorBit:0)> TransposedSparseMatrix; typedef DynamicSparseMatrix<Scalar,(Flags&~RowMajorBit)|(IsRowMajor?RowMajorBit:0)> TransposedSparseMatrix;
int m_innerSize; int m_innerSize;

View File

@ -43,16 +43,21 @@ template<typename MatrixType, int Size>
class SparseInnerVectorSet : ei_no_assignment_operator, class SparseInnerVectorSet : ei_no_assignment_operator,
public SparseMatrixBase<SparseInnerVectorSet<MatrixType, Size> > public SparseMatrixBase<SparseInnerVectorSet<MatrixType, Size> >
{ {
enum { IsRowMajor = ei_traits<SparseInnerVectorSet>::IsRowMajor };
public: public:
enum { IsRowMajor = ei_traits<SparseInnerVectorSet>::IsRowMajor };
EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseInnerVectorSet) EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseInnerVectorSet)
class InnerIterator: public MatrixType::InnerIterator class InnerIterator: public MatrixType::InnerIterator
{ {
public: public:
inline InnerIterator(const SparseInnerVectorSet& xpr, int outer) inline InnerIterator(const SparseInnerVectorSet& xpr, int outer)
: MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer) : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer)
{} {}
inline int row() const { return IsRowMajor ? m_outer : this->index(); }
inline int col() const { return IsRowMajor ? this->index() : m_outer; }
protected:
int m_outer;
}; };
inline SparseInnerVectorSet(const MatrixType& matrix, int outerStart, int outerSize) inline SparseInnerVectorSet(const MatrixType& matrix, int outerStart, int outerSize)
@ -100,16 +105,21 @@ class SparseInnerVectorSet<DynamicSparseMatrix<_Scalar, _Options>, Size>
: public SparseMatrixBase<SparseInnerVectorSet<DynamicSparseMatrix<_Scalar, _Options>, Size> > : public SparseMatrixBase<SparseInnerVectorSet<DynamicSparseMatrix<_Scalar, _Options>, Size> >
{ {
typedef DynamicSparseMatrix<_Scalar, _Options> MatrixType; typedef DynamicSparseMatrix<_Scalar, _Options> MatrixType;
enum { IsRowMajor = ei_traits<SparseInnerVectorSet>::IsRowMajor };
public: public:
enum { IsRowMajor = ei_traits<SparseInnerVectorSet>::IsRowMajor };
EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseInnerVectorSet) EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseInnerVectorSet)
class InnerIterator: public MatrixType::InnerIterator class InnerIterator: public MatrixType::InnerIterator
{ {
public: public:
inline InnerIterator(const SparseInnerVectorSet& xpr, int outer) inline InnerIterator(const SparseInnerVectorSet& xpr, int outer)
: MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer) : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer)
{} {}
inline int row() const { return IsRowMajor ? m_outer : this->index(); }
inline int col() const { return IsRowMajor ? this->index() : m_outer; }
protected:
int m_outer;
}; };
inline SparseInnerVectorSet(const MatrixType& matrix, int outerStart, int outerSize) inline SparseInnerVectorSet(const MatrixType& matrix, int outerStart, int outerSize)
@ -193,16 +203,21 @@ class SparseInnerVectorSet<SparseMatrix<_Scalar, _Options>, Size>
: public SparseMatrixBase<SparseInnerVectorSet<SparseMatrix<_Scalar, _Options>, Size> > : public SparseMatrixBase<SparseInnerVectorSet<SparseMatrix<_Scalar, _Options>, Size> >
{ {
typedef SparseMatrix<_Scalar, _Options> MatrixType; typedef SparseMatrix<_Scalar, _Options> MatrixType;
enum { IsRowMajor = ei_traits<SparseInnerVectorSet>::IsRowMajor };
public: public:
enum { IsRowMajor = ei_traits<SparseInnerVectorSet>::IsRowMajor };
EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseInnerVectorSet) EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseInnerVectorSet)
class InnerIterator: public MatrixType::InnerIterator class InnerIterator: public MatrixType::InnerIterator
{ {
public: public:
inline InnerIterator(const SparseInnerVectorSet& xpr, int outer) inline InnerIterator(const SparseInnerVectorSet& xpr, int outer)
: MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer) : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer)
{} {}
inline int row() const { return IsRowMajor ? m_outer : this->index(); }
inline int col() const { return IsRowMajor ? this->index() : m_outer; }
protected:
int m_outer;
}; };
inline SparseInnerVectorSet(const MatrixType& matrix, int outerStart, int outerSize) inline SparseInnerVectorSet(const MatrixType& matrix, int outerStart, int outerSize)

View File

@ -186,8 +186,8 @@ class ei_sparse_cwise_binary_op_inner_iterator_selector<BinaryOp, Lhs, Rhs, Deri
EIGEN_STRONG_INLINE Scalar value() const { return m_value; } EIGEN_STRONG_INLINE Scalar value() const { return m_value; }
EIGEN_STRONG_INLINE int index() const { return m_id; } EIGEN_STRONG_INLINE int index() const { return m_id; }
EIGEN_STRONG_INLINE int row() const { return m_lhsIter.row(); } EIGEN_STRONG_INLINE int row() const { return Lhs::IsRowMajor ? m_lhsIter.row() : index(); }
EIGEN_STRONG_INLINE int col() const { return m_lhsIter.col(); } EIGEN_STRONG_INLINE int col() const { return Lhs::IsRowMajor ? index() : m_lhsIter.col(); }
EIGEN_STRONG_INLINE operator bool() const { return m_id>=0; } EIGEN_STRONG_INLINE operator bool() const { return m_id>=0; }

View File

@ -25,17 +25,24 @@
#ifndef EIGEN_SPARSEMATRIX_H #ifndef EIGEN_SPARSEMATRIX_H
#define EIGEN_SPARSEMATRIX_H #define EIGEN_SPARSEMATRIX_H
/** \class SparseMatrix /** \ingroup Sparse_Module
* *
* \brief Sparse matrix * \class SparseMatrix
*
* \brief The main sparse matrix class
*
* This class implements a sparse matrix using the very common compressed row/column storage
* scheme.
* *
* \param _Scalar the scalar type, i.e. the type of the coefficients * \param _Scalar the scalar type, i.e. the type of the coefficients
* \param _Options Union of bit flags controlling the storage scheme. Currently the only possibility
* is RowMajor. The default is 0 which means column-major.
* *
* See http://www.netlib.org/linalg/html_templates/node91.html for details on the storage scheme. * See http://www.netlib.org/linalg/html_templates/node91.html for details on the storage scheme.
* *
*/ */
template<typename _Scalar, int _Flags> template<typename _Scalar, int _Options>
struct ei_traits<SparseMatrix<_Scalar, _Flags> > struct ei_traits<SparseMatrix<_Scalar, _Options> >
{ {
typedef _Scalar Scalar; typedef _Scalar Scalar;
enum { enum {
@ -43,17 +50,15 @@ struct ei_traits<SparseMatrix<_Scalar, _Flags> >
ColsAtCompileTime = Dynamic, ColsAtCompileTime = Dynamic,
MaxRowsAtCompileTime = Dynamic, MaxRowsAtCompileTime = Dynamic,
MaxColsAtCompileTime = Dynamic, MaxColsAtCompileTime = Dynamic,
Flags = SparseBit | _Flags, Flags = SparseBit | _Options,
CoeffReadCost = NumTraits<Scalar>::ReadCost, CoeffReadCost = NumTraits<Scalar>::ReadCost,
SupportedAccessPatterns = InnerRandomAccessPattern SupportedAccessPatterns = InnerRandomAccessPattern
}; };
}; };
template<typename _Scalar, int _Options>
template<typename _Scalar, int _Flags>
class SparseMatrix class SparseMatrix
: public SparseMatrixBase<SparseMatrix<_Scalar, _Flags> > : public SparseMatrixBase<SparseMatrix<_Scalar, _Options> >
{ {
public: public:
EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseMatrix) EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseMatrix)
@ -64,10 +69,10 @@ class SparseMatrix
// EIGEN_SPARSE_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(SparseMatrix, /=) // EIGEN_SPARSE_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(SparseMatrix, /=)
typedef MappedSparseMatrix<Scalar,Flags> Map; typedef MappedSparseMatrix<Scalar,Flags> Map;
using Base::IsRowMajor;
protected: protected:
enum { IsRowMajor = Base::IsRowMajor };
typedef SparseMatrix<Scalar,(Flags&~RowMajorBit)|(IsRowMajor?RowMajorBit:0)> TransposedSparseMatrix; typedef SparseMatrix<Scalar,(Flags&~RowMajorBit)|(IsRowMajor?RowMajorBit:0)> TransposedSparseMatrix;
int m_outerSize; int m_outerSize;
@ -508,10 +513,13 @@ class SparseMatrix
{ {
delete[] m_outerIndex; delete[] m_outerIndex;
} }
/** Overloaded for performance */
Scalar sum() const;
}; };
template<typename Scalar, int _Flags> template<typename Scalar, int _Options>
class SparseMatrix<Scalar,_Flags>::InnerIterator class SparseMatrix<Scalar,_Options>::InnerIterator
{ {
public: public:
InnerIterator(const SparseMatrix& mat, int outer) InnerIterator(const SparseMatrix& mat, int outer)

View File

@ -25,6 +25,17 @@
#ifndef EIGEN_SPARSEMATRIXBASE_H #ifndef EIGEN_SPARSEMATRIXBASE_H
#define EIGEN_SPARSEMATRIXBASE_H #define EIGEN_SPARSEMATRIXBASE_H
/** \ingroup Sparse_Module
*
* \class SparseMatrixBase
*
* \brief Base class of any sparse matrices or sparse expressions
*
* \param Derived
*
*
*
*/
template<typename Derived> class SparseMatrixBase template<typename Derived> class SparseMatrixBase
{ {
public: public:
@ -432,10 +443,16 @@ template<typename Derived> class SparseMatrixBase
for (int j=0; j<outerSize(); ++j) for (int j=0; j<outerSize(); ++j)
{ {
for (typename Derived::InnerIterator i(derived(),j); i; ++i) for (typename Derived::InnerIterator i(derived(),j); i; ++i)
{
if(IsRowMajor) if(IsRowMajor)
res.coeffRef(j,i.index()) = i.value(); std::cerr << i.row() << "," << i.col() << " == " << j << "," << i.index() << "\n";
else else
res.coeffRef(i.index(),j) = i.value(); std::cerr << i.row() << "," << i.col() << " == " << i.index() << "," << j << "\n";
// if(IsRowMajor)
res.coeffRef(i.row(),i.col()) = i.value();
// else
// res.coeffRef(i.index(),j) = i.value();
}
} }
return res; return res;
} }

View File

@ -37,4 +37,20 @@ SparseMatrixBase<Derived>::sum() const
return res; return res;
} }
template<typename _Scalar, int _Options>
typename ei_traits<SparseMatrix<_Scalar,_Options> >::Scalar
SparseMatrix<_Scalar,_Options>::sum() const
{
ei_assert(rows()>0 && cols()>0 && "you are using a non initialized matrix");
return Matrix<Scalar,1,Dynamic>::Map(m_data.value(0), m_data.size()).sum();
}
template<typename _Scalar, int _Options>
typename ei_traits<SparseVector<_Scalar,_Options> >::Scalar
SparseVector<_Scalar,_Options>::sum() const
{
ei_assert(rows()>0 && cols()>0 && "you are using a non initialized matrix");
return Matrix<Scalar,1,Dynamic>::Map(m_data.value(0), m_data.size()).sum();
}
#endif // EIGEN_SPARSEREDUX_H #endif // EIGEN_SPARSEREDUX_H

View File

@ -66,20 +66,26 @@ template<typename MatrixType> class SparseTranspose
template<typename MatrixType> class SparseTranspose<MatrixType>::InnerIterator : public MatrixType::InnerIterator template<typename MatrixType> class SparseTranspose<MatrixType>::InnerIterator : public MatrixType::InnerIterator
{ {
typedef typename MatrixType::InnerIterator Base;
public: public:
EIGEN_STRONG_INLINE InnerIterator(const SparseTranspose& trans, int outer) EIGEN_STRONG_INLINE InnerIterator(const SparseTranspose& trans, int outer)
: MatrixType::InnerIterator(trans.m_matrix, outer) : Base(trans.m_matrix, outer)
{} {}
inline int row() const { return Base::col(); }
inline int col() const { return Base::row(); }
}; };
template<typename MatrixType> class SparseTranspose<MatrixType>::ReverseInnerIterator : public MatrixType::ReverseInnerIterator template<typename MatrixType> class SparseTranspose<MatrixType>::ReverseInnerIterator : public MatrixType::ReverseInnerIterator
{ {
typedef typename MatrixType::ReverseInnerIterator Base;
public: public:
EIGEN_STRONG_INLINE ReverseInnerIterator(const SparseTranspose& xpr, int outer) EIGEN_STRONG_INLINE ReverseInnerIterator(const SparseTranspose& xpr, int outer)
: MatrixType::ReverseInnerIterator(xpr.m_matrix, outer) : Base(xpr.m_matrix, outer)
{} {}
inline int row() const { return Base::col(); }
inline int col() const { return Base::row(); }
}; };
#endif // EIGEN_SPARSETRANSPOSE_H #endif // EIGEN_SPARSETRANSPOSE_H

View File

@ -34,26 +34,26 @@
* See http://www.netlib.org/linalg/html_templates/node91.html for details on the storage scheme. * See http://www.netlib.org/linalg/html_templates/node91.html for details on the storage scheme.
* *
*/ */
template<typename _Scalar, int _Flags> template<typename _Scalar, int _Options>
struct ei_traits<SparseVector<_Scalar, _Flags> > struct ei_traits<SparseVector<_Scalar, _Options> >
{ {
typedef _Scalar Scalar; typedef _Scalar Scalar;
enum { enum {
IsColVector = _Flags & RowMajorBit ? 0 : 1, IsColVector = _Options & RowMajorBit ? 0 : 1,
RowsAtCompileTime = IsColVector ? Dynamic : 1, RowsAtCompileTime = IsColVector ? Dynamic : 1,
ColsAtCompileTime = IsColVector ? 1 : Dynamic, ColsAtCompileTime = IsColVector ? 1 : Dynamic,
MaxRowsAtCompileTime = RowsAtCompileTime, MaxRowsAtCompileTime = RowsAtCompileTime,
MaxColsAtCompileTime = ColsAtCompileTime, MaxColsAtCompileTime = ColsAtCompileTime,
Flags = SparseBit | _Flags, Flags = SparseBit | _Options,
CoeffReadCost = NumTraits<Scalar>::ReadCost, CoeffReadCost = NumTraits<Scalar>::ReadCost,
SupportedAccessPatterns = InnerRandomAccessPattern SupportedAccessPatterns = InnerRandomAccessPattern
}; };
}; };
template<typename _Scalar, int _Flags> template<typename _Scalar, int _Options>
class SparseVector class SparseVector
: public SparseMatrixBase<SparseVector<_Scalar, _Flags> > : public SparseMatrixBase<SparseVector<_Scalar, _Options> >
{ {
public: public:
EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseVector) EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseVector)
@ -357,10 +357,13 @@ class SparseVector
/** Destructor */ /** Destructor */
inline ~SparseVector() {} inline ~SparseVector() {}
/** Overloaded for performance */
Scalar sum() const;
}; };
template<typename Scalar, int _Flags> template<typename Scalar, int _Options>
class SparseVector<Scalar,_Flags>::InnerIterator class SparseVector<Scalar,_Options>::InnerIterator
{ {
public: public:
InnerIterator(const SparseVector& vec, int outer=0) InnerIterator(const SparseVector& vec, int outer=0)

7
eigen2.pc.in Normal file
View File

@ -0,0 +1,7 @@
Name: Eigen2
Description: A C++ template library for linear algebra: vectors, matrices, and related algorithms
Requires:
Version: ${EIGEN_VERSION_NUMBER}
Libs:
Cflags: -I${INCLUDE_INSTALL_DIR}

View File

@ -39,7 +39,7 @@
# VERSION=opensuse-11.1 # VERSION=opensuse-11.1
# WORK_DIR=/home/gael/Coding/eigen2/cdash # WORK_DIR=/home/gael/Coding/eigen2/cdash
# # get the last version of the script # # get the last version of the script
# svn cat svn://anonsvn.kde.org/home/kde/trunk/kdesupport/eigen2/test/testsuite.cmake > $WORK_DIR/testsuite.cmake # wget http://bitbucket.org/eigen/eigen2/raw/tip/test/testsuite.cmake -o $WORK_DIR/testsuite.cmake
# COMMON="ctest -S $WORK_DIR/testsuite.cmake,EIGEN_WORK_DIR=$WORK_DIR,EIGEN_SITE=$SITE,EIGEN_MODE=$1,EIGEN_BUILD_STRING=$OS_VERSION-$ARCH" # COMMON="ctest -S $WORK_DIR/testsuite.cmake,EIGEN_WORK_DIR=$WORK_DIR,EIGEN_SITE=$SITE,EIGEN_MODE=$1,EIGEN_BUILD_STRING=$OS_VERSION-$ARCH"
# $COMMON-gcc-3.4.6,EIGEN_CXX=g++-3.4 # $COMMON-gcc-3.4.6,EIGEN_CXX=g++-3.4
# $COMMON-gcc-4.0.1,EIGEN_CXX=g++-4.0.1 # $COMMON-gcc-4.0.1,EIGEN_CXX=g++-4.0.1
@ -132,8 +132,8 @@ endif(NOT EIGEN_MODE)
## mandatory variables (the default should be ok in most cases): ## mandatory variables (the default should be ok in most cases):
SET (CTEST_CVS_COMMAND "svn") SET (CTEST_CVS_COMMAND "hg")
SET (CTEST_CVS_CHECKOUT "${CTEST_CVS_COMMAND} co svn://anonsvn.kde.org/home/kde/trunk/kdesupport/eigen2 \"${CTEST_SOURCE_DIRECTORY}\"") SET (CTEST_CVS_CHECKOUT "${CTEST_CVS_COMMAND} clone http://bitbucket.org/eigen/eigen2 \"${CTEST_SOURCE_DIRECTORY}\"")
# which ctest command to use for running the dashboard # which ctest command to use for running the dashboard
SET (CTEST_COMMAND "${EIGEN_CMAKE_DIR}ctest -D ${EIGEN_MODE}") SET (CTEST_COMMAND "${EIGEN_CMAKE_DIR}ctest -D ${EIGEN_MODE}")

View File

@ -72,7 +72,8 @@ void ei_matrix_exponential(const MatrixBase<Derived> &M, typename ei_plain_matri
PlainMatrixType num, den, U, V; PlainMatrixType num, den, U, V;
PlainMatrixType Id = PlainMatrixType::Identity(M.rows(), M.cols()); PlainMatrixType Id = PlainMatrixType::Identity(M.rows(), M.cols());
RealScalar l1norm = M.cwise().abs().colwise().sum().maxCoeff(); typename ei_eval<Derived>::type Meval = M.eval();
RealScalar l1norm = Meval.cwise().abs().colwise().sum().maxCoeff();
int squarings = 0; int squarings = 0;
// Choose degree of Pade approximant, depending on norm of M // Choose degree of Pade approximant, depending on norm of M
@ -81,9 +82,9 @@ void ei_matrix_exponential(const MatrixBase<Derived> &M, typename ei_plain_matri
// Use (3,3)-Pade // Use (3,3)-Pade
const Scalar b[] = {120., 60., 12., 1.}; const Scalar b[] = {120., 60., 12., 1.};
PlainMatrixType M2; PlainMatrixType M2;
M2 = (M * M).lazy(); M2 = (Meval * Meval).lazy();
num = b[3]*M2 + b[1]*Id; num = b[3]*M2 + b[1]*Id;
U = (M * num).lazy(); U = (Meval * num).lazy();
V = b[2]*M2 + b[0]*Id; V = b[2]*M2 + b[0]*Id;
} else if (l1norm < 2.539398330063230e-001) { } else if (l1norm < 2.539398330063230e-001) {
@ -91,10 +92,10 @@ void ei_matrix_exponential(const MatrixBase<Derived> &M, typename ei_plain_matri
// Use (5,5)-Pade // Use (5,5)-Pade
const Scalar b[] = {30240., 15120., 3360., 420., 30., 1.}; const Scalar b[] = {30240., 15120., 3360., 420., 30., 1.};
PlainMatrixType M2, M4; PlainMatrixType M2, M4;
M2 = (M * M).lazy(); M2 = (Meval * Meval).lazy();
M4 = (M2 * M2).lazy(); M4 = (M2 * M2).lazy();
num = b[5]*M4 + b[3]*M2 + b[1]*Id; num = b[5]*M4 + b[3]*M2 + b[1]*Id;
U = (M * num).lazy(); U = (Meval * num).lazy();
V = b[4]*M4 + b[2]*M2 + b[0]*Id; V = b[4]*M4 + b[2]*M2 + b[0]*Id;
} else if (l1norm < 9.504178996162932e-001) { } else if (l1norm < 9.504178996162932e-001) {
@ -102,11 +103,11 @@ void ei_matrix_exponential(const MatrixBase<Derived> &M, typename ei_plain_matri
// Use (7,7)-Pade // Use (7,7)-Pade
const Scalar b[] = {17297280., 8648640., 1995840., 277200., 25200., 1512., 56., 1.}; const Scalar b[] = {17297280., 8648640., 1995840., 277200., 25200., 1512., 56., 1.};
PlainMatrixType M2, M4, M6; PlainMatrixType M2, M4, M6;
M2 = (M * M).lazy(); M2 = (Meval * Meval).lazy();
M4 = (M2 * M2).lazy(); M4 = (M2 * M2).lazy();
M6 = (M4 * M2).lazy(); M6 = (M4 * M2).lazy();
num = b[7]*M6 + b[5]*M4 + b[3]*M2 + b[1]*Id; num = b[7]*M6 + b[5]*M4 + b[3]*M2 + b[1]*Id;
U = (M * num).lazy(); U = (Meval * num).lazy();
V = b[6]*M6 + b[4]*M4 + b[2]*M2 + b[0]*Id; V = b[6]*M6 + b[4]*M4 + b[2]*M2 + b[0]*Id;
} else if (l1norm < 2.097847961257068e+000) { } else if (l1norm < 2.097847961257068e+000) {
@ -115,12 +116,12 @@ void ei_matrix_exponential(const MatrixBase<Derived> &M, typename ei_plain_matri
const Scalar b[] = {17643225600., 8821612800., 2075673600., 302702400., 30270240., const Scalar b[] = {17643225600., 8821612800., 2075673600., 302702400., 30270240.,
2162160., 110880., 3960., 90., 1.}; 2162160., 110880., 3960., 90., 1.};
PlainMatrixType M2, M4, M6, M8; PlainMatrixType M2, M4, M6, M8;
M2 = (M * M).lazy(); M2 = (Meval * Meval).lazy();
M4 = (M2 * M2).lazy(); M4 = (M2 * M2).lazy();
M6 = (M4 * M2).lazy(); M6 = (M4 * M2).lazy();
M8 = (M6 * M2).lazy(); M8 = (M6 * M2).lazy();
num = b[9]*M8 + b[7]*M6 + b[5]*M4 + b[3]*M2 + b[1]*Id; num = b[9]*M8 + b[7]*M6 + b[5]*M4 + b[3]*M2 + b[1]*Id;
U = (M * num).lazy(); U = (Meval * num).lazy();
V = b[8]*M8 + b[6]*M6 + b[4]*M4 + b[2]*M2 + b[0]*Id; V = b[8]*M8 + b[6]*M6 + b[4]*M4 + b[2]*M2 + b[0]*Id;
} else { } else {
@ -135,7 +136,7 @@ void ei_matrix_exponential(const MatrixBase<Derived> &M, typename ei_plain_matri
squarings = std::max(0, (int)ceil(log2(l1norm / maxnorm))); squarings = std::max(0, (int)ceil(log2(l1norm / maxnorm)));
PlainMatrixType A, A2, A4, A6; PlainMatrixType A, A2, A4, A6;
A = M / pow(Scalar(2), squarings); A = Meval / pow(Scalar(2), squarings);
A2 = (A * A).lazy(); A2 = (A * A).lazy();
A4 = (A2 * A2).lazy(); A4 = (A2 * A2).lazy();
A6 = (A4 * A2).lazy(); A6 = (A4 * A2).lazy();