[Rcpp-commits] r4052 - in pkg/RcppEigen: . inst/include/Eigen/src/Core inst/include/Eigen/src/Core/products inst/include/Eigen/src/Core/util inst/include/Eigen/src/Eigenvalues inst/include/Eigen/src/Geometry inst/include/Eigen/src/PardisoSupport inst/include/Eigen/src/QR inst/include/Eigen/src/SVD inst/include/Eigen/src/SparseCore inst/include/Eigen/src/SuperLUSupport inst/include/Eigen/src/plugins inst/include/unsupported/Eigen/src/MatrixFunctions
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Nov 29 19:00:38 CET 2012
Author: dmbates
Date: 2012-11-29 19:00:37 +0100 (Thu, 29 Nov 2012)
New Revision: 4052
Removed:
pkg/RcppEigen/inst/include/Eigen/src/Eigenvalues/EigenvaluesCommon.h
Modified:
pkg/RcppEigen/ChangeLog
pkg/RcppEigen/DESCRIPTION
pkg/RcppEigen/inst/include/Eigen/src/Core/ArrayWrapper.h
pkg/RcppEigen/inst/include/Eigen/src/Core/CommaInitializer.h
pkg/RcppEigen/inst/include/Eigen/src/Core/DiagonalMatrix.h
pkg/RcppEigen/inst/include/Eigen/src/Core/Functors.h
pkg/RcppEigen/inst/include/Eigen/src/Core/MatrixBase.h
pkg/RcppEigen/inst/include/Eigen/src/Core/StableNorm.h
pkg/RcppEigen/inst/include/Eigen/src/Core/TriangularMatrix.h
pkg/RcppEigen/inst/include/Eigen/src/Core/products/GeneralMatrixVector.h
pkg/RcppEigen/inst/include/Eigen/src/Core/products/TriangularMatrixMatrix_MKL.h
pkg/RcppEigen/inst/include/Eigen/src/Core/products/TriangularMatrixVector_MKL.h
pkg/RcppEigen/inst/include/Eigen/src/Core/util/Macros.h
pkg/RcppEigen/inst/include/Eigen/src/Core/util/Memory.h
pkg/RcppEigen/inst/include/Eigen/src/Core/util/XprHelper.h
pkg/RcppEigen/inst/include/Eigen/src/Eigenvalues/ComplexSchur.h
pkg/RcppEigen/inst/include/Eigen/src/Eigenvalues/ComplexSchur_MKL.h
pkg/RcppEigen/inst/include/Eigen/src/Eigenvalues/RealSchur.h
pkg/RcppEigen/inst/include/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
pkg/RcppEigen/inst/include/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h
pkg/RcppEigen/inst/include/Eigen/src/Geometry/AlignedBox.h
pkg/RcppEigen/inst/include/Eigen/src/PardisoSupport/PardisoSupport.h
pkg/RcppEigen/inst/include/Eigen/src/QR/ColPivHouseholderQR_MKL.h
pkg/RcppEigen/inst/include/Eigen/src/SVD/JacobiSVD_MKL.h
pkg/RcppEigen/inst/include/Eigen/src/SparseCore/SparseDiagonalProduct.h
pkg/RcppEigen/inst/include/Eigen/src/SparseCore/SparseMatrix.h
pkg/RcppEigen/inst/include/Eigen/src/SparseCore/SparseUtil.h
pkg/RcppEigen/inst/include/Eigen/src/SuperLUSupport/SuperLUSupport.h
pkg/RcppEigen/inst/include/Eigen/src/plugins/ArrayCwiseBinaryOps.h
pkg/RcppEigen/inst/include/Eigen/src/plugins/ArrayCwiseUnaryOps.h
pkg/RcppEigen/inst/include/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h
Log:
Upgrade to version 3.1.2 of Eigen, change RcppEigen version number to 0.3.1.2 to reflect this.
Modified: pkg/RcppEigen/ChangeLog
===================================================================
--- pkg/RcppEigen/ChangeLog 2012-11-29 15:32:15 UTC (rev 4051)
+++ pkg/RcppEigen/ChangeLog 2012-11-29 18:00:37 UTC (rev 4052)
@@ -1,3 +1,9 @@
+2012-11-29 Douglas Bates <bates at stat.wisc.edu>
+
+ * DESCRIPTION: New version
+ * inst/include/Eigen, inst/include/unsupported: Update to version
+ 3.1.2 of Eigen.
+
2012-11-29 Romain Francois <romain at r-enthusiasts.com>
* include/RcppEigenWrap.h: compatibility issue with Rcpp 0.10.1
Modified: pkg/RcppEigen/DESCRIPTION
===================================================================
--- pkg/RcppEigen/DESCRIPTION 2012-11-29 15:32:15 UTC (rev 4051)
+++ pkg/RcppEigen/DESCRIPTION 2012-11-29 18:00:37 UTC (rev 4052)
@@ -1,8 +1,8 @@
Package: RcppEigen
Type: Package
Title: Rcpp integration for the Eigen templated linear algebra library.
-Version: 0.3.2
-Date: 2012-08-07
+Version: 0.3.1.2
+Date: 2012-11-29
Author: Douglas Bates, Romain Francois and Dirk Eddelbuettel
Maintainer: Douglas Bates <bates at stat.wisc.edu>
Description: R and Eigen integration using Rcpp.
Modified: pkg/RcppEigen/inst/include/Eigen/src/Core/ArrayWrapper.h
===================================================================
--- pkg/RcppEigen/inst/include/Eigen/src/Core/ArrayWrapper.h 2012-11-29 15:32:15 UTC (rev 4051)
+++ pkg/RcppEigen/inst/include/Eigen/src/Core/ArrayWrapper.h 2012-11-29 18:00:37 UTC (rev 4052)
@@ -121,6 +121,13 @@
return m_expression;
}
+ /** Forwards the resizing request to the nested expression
+ * \sa DenseBase::resize(Index) */
+ void resize(Index newSize) { m_expression.const_cast_derived().resize(newSize); }
+ /** Forwards the resizing request to the nested expression
+ * \sa DenseBase::resize(Index,Index)*/
+ void resize(Index nbRows, Index nbCols) { m_expression.const_cast_derived().resize(nbRows,nbCols); }
+
protected:
NestedExpressionType m_expression;
};
@@ -231,6 +238,13 @@
return m_expression;
}
+ /** Forwards the resizing request to the nested expression
+ * \sa DenseBase::resize(Index) */
+ void resize(Index newSize) { m_expression.const_cast_derived().resize(newSize); }
+ /** Forwards the resizing request to the nested expression
+ * \sa DenseBase::resize(Index,Index)*/
+ void resize(Index nbRows, Index nbCols) { m_expression.const_cast_derived().resize(nbRows,nbCols); }
+
protected:
NestedExpressionType m_expression;
};
Modified: pkg/RcppEigen/inst/include/Eigen/src/Core/CommaInitializer.h
===================================================================
--- pkg/RcppEigen/inst/include/Eigen/src/Core/CommaInitializer.h 2012-11-29 15:32:15 UTC (rev 4051)
+++ pkg/RcppEigen/inst/include/Eigen/src/Core/CommaInitializer.h 2012-11-29 18:00:37 UTC (rev 4052)
@@ -65,6 +65,8 @@
template<typename OtherDerived>
CommaInitializer& operator,(const DenseBase<OtherDerived>& other)
{
+ if(other.cols()==0 || other.rows()==0)
+ return *this;
if (m_col==m_xpr.cols())
{
m_row+=m_currentBlockRows;
Modified: pkg/RcppEigen/inst/include/Eigen/src/Core/DiagonalMatrix.h
===================================================================
--- pkg/RcppEigen/inst/include/Eigen/src/Core/DiagonalMatrix.h 2012-11-29 15:32:15 UTC (rev 4051)
+++ pkg/RcppEigen/inst/include/Eigen/src/Core/DiagonalMatrix.h 2012-11-29 18:00:37 UTC (rev 4052)
@@ -20,6 +20,7 @@
public:
typedef typename internal::traits<Derived>::DiagonalVectorType DiagonalVectorType;
typedef typename DiagonalVectorType::Scalar Scalar;
+ typedef typename DiagonalVectorType::RealScalar RealScalar;
typedef typename internal::traits<Derived>::StorageKind StorageKind;
typedef typename internal::traits<Derived>::Index Index;
@@ -65,6 +66,17 @@
return diagonal().cwiseInverse();
}
+ inline const DiagonalWrapper<const CwiseUnaryOp<internal::scalar_multiple_op<Scalar>, const DiagonalVectorType> >
+ operator*(const Scalar& scalar) const
+ {
+ return diagonal() * scalar;
+ }
+ friend inline const DiagonalWrapper<const CwiseUnaryOp<internal::scalar_multiple_op<Scalar>, const DiagonalVectorType> >
+ operator*(const Scalar& scalar, const DiagonalBase& other)
+ {
+ return other.diagonal() * scalar;
+ }
+
#ifdef EIGEN2_SUPPORT
template<typename OtherDerived>
bool isApprox(const DiagonalBase<OtherDerived>& other, typename NumTraits<Scalar>::Real precision = NumTraits<Scalar>::dummy_precision()) const
Modified: pkg/RcppEigen/inst/include/Eigen/src/Core/Functors.h
===================================================================
--- pkg/RcppEigen/inst/include/Eigen/src/Core/Functors.h 2012-11-29 15:32:15 UTC (rev 4051)
+++ pkg/RcppEigen/inst/include/Eigen/src/Core/Functors.h 2012-11-29 18:00:37 UTC (rev 4052)
@@ -447,7 +447,7 @@
* indeed it seems better to declare m_other as a Packet and do the pset1() once
* in the constructor. However, in practice:
* - GCC does not like m_other as a Packet and generate a load every time it needs it
- * - on the other hand GCC is able to moves the pset1() away the loop :)
+ * - on the other hand GCC is able to moves the pset1() outside the loop :)
* - simpler code ;)
* (ICC and gcc 4.4 seems to perform well in both cases, the issue is visible with y = a*x + b*y)
*/
@@ -478,33 +478,6 @@
struct functor_traits<scalar_multiple2_op<Scalar1,Scalar2> >
{ enum { Cost = NumTraits<Scalar1>::MulCost, PacketAccess = false }; };
-template<typename Scalar, bool IsInteger>
-struct scalar_quotient1_impl {
- typedef typename packet_traits<Scalar>::type Packet;
- // FIXME default copy constructors seems bugged with std::complex<>
- EIGEN_STRONG_INLINE scalar_quotient1_impl(const scalar_quotient1_impl& other) : m_other(other.m_other) { }
- EIGEN_STRONG_INLINE scalar_quotient1_impl(const Scalar& other) : m_other(static_cast<Scalar>(1) / other) {}
- EIGEN_STRONG_INLINE Scalar operator() (const Scalar& a) const { return a * m_other; }
- EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a) const
- { return internal::pmul(a, pset1<Packet>(m_other)); }
- const Scalar m_other;
-};
-template<typename Scalar>
-struct functor_traits<scalar_quotient1_impl<Scalar,false> >
-{ enum { Cost = NumTraits<Scalar>::MulCost, PacketAccess = packet_traits<Scalar>::HasMul }; };
-
-template<typename Scalar>
-struct scalar_quotient1_impl<Scalar,true> {
- // FIXME default copy constructors seems bugged with std::complex<>
- EIGEN_STRONG_INLINE scalar_quotient1_impl(const scalar_quotient1_impl& other) : m_other(other.m_other) { }
- EIGEN_STRONG_INLINE scalar_quotient1_impl(const Scalar& other) : m_other(other) {}
- EIGEN_STRONG_INLINE Scalar operator() (const Scalar& a) const { return a / m_other; }
- typename add_const_on_value_type<typename NumTraits<Scalar>::Nested>::type m_other;
-};
-template<typename Scalar>
-struct functor_traits<scalar_quotient1_impl<Scalar,true> >
-{ enum { Cost = 2 * NumTraits<Scalar>::MulCost, PacketAccess = false }; };
-
/** \internal
* \brief Template functor to divide a scalar by a fixed other one
*
@@ -514,14 +487,19 @@
* \sa class CwiseUnaryOp, MatrixBase::operator/
*/
template<typename Scalar>
-struct scalar_quotient1_op : scalar_quotient1_impl<Scalar, NumTraits<Scalar>::IsInteger > {
- EIGEN_STRONG_INLINE scalar_quotient1_op(const Scalar& other)
- : scalar_quotient1_impl<Scalar, NumTraits<Scalar>::IsInteger >(other) {}
+struct scalar_quotient1_op {
+ typedef typename packet_traits<Scalar>::type Packet;
+ // FIXME default copy constructors seems bugged with std::complex<>
+ EIGEN_STRONG_INLINE scalar_quotient1_op(const scalar_quotient1_op& other) : m_other(other.m_other) { }
+ EIGEN_STRONG_INLINE scalar_quotient1_op(const Scalar& other) : m_other(other) {}
+ EIGEN_STRONG_INLINE Scalar operator() (const Scalar& a) const { return a / m_other; }
+ EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a) const
+ { return internal::pdiv(a, pset1<Packet>(m_other)); }
+ typename add_const_on_value_type<typename NumTraits<Scalar>::Nested>::type m_other;
};
template<typename Scalar>
struct functor_traits<scalar_quotient1_op<Scalar> >
-: functor_traits<scalar_quotient1_impl<Scalar, NumTraits<Scalar>::IsInteger> >
-{};
+{ enum { Cost = 2 * NumTraits<Scalar>::MulCost, PacketAccess = packet_traits<Scalar>::HasDiv }; };
// nullary functors
Modified: pkg/RcppEigen/inst/include/Eigen/src/Core/MatrixBase.h
===================================================================
--- pkg/RcppEigen/inst/include/Eigen/src/Core/MatrixBase.h 2012-11-29 15:32:15 UTC (rev 4051)
+++ pkg/RcppEigen/inst/include/Eigen/src/Core/MatrixBase.h 2012-11-29 18:00:37 UTC (rev 4052)
@@ -237,7 +237,7 @@
// huuuge hack. make Eigen2's matrix.part<Diagonal>() work in eigen3. Problem: Diagonal is now a class template instead
// of an integer constant. Solution: overload the part() method template wrt template parameters list.
- template<template<typename T, int n> class U>
+ template<template<typename T, int N> class U>
const DiagonalWrapper<ConstDiagonalReturnType> part() const
{ return diagonal().asDiagonal(); }
#endif // EIGEN2_SUPPORT
Modified: pkg/RcppEigen/inst/include/Eigen/src/Core/StableNorm.h
===================================================================
--- pkg/RcppEigen/inst/include/Eigen/src/Core/StableNorm.h 2012-11-29 15:32:15 UTC (rev 4051)
+++ pkg/RcppEigen/inst/include/Eigen/src/Core/StableNorm.h 2012-11-29 18:00:37 UTC (rev 4052)
@@ -131,7 +131,6 @@
abig = internal::sqrt(abig);
if(abig > overfl)
{
- eigen_assert(false && "overflow");
return rbig;
}
if(amed > RealScalar(0))
Modified: pkg/RcppEigen/inst/include/Eigen/src/Core/TriangularMatrix.h
===================================================================
--- pkg/RcppEigen/inst/include/Eigen/src/Core/TriangularMatrix.h 2012-11-29 15:32:15 UTC (rev 4051)
+++ pkg/RcppEigen/inst/include/Eigen/src/Core/TriangularMatrix.h 2012-11-29 18:00:37 UTC (rev 4052)
@@ -511,6 +511,7 @@
struct triangular_assignment_selector<Derived1, Derived2, StrictlyUpper, Dynamic, ClearOpposite>
{
typedef typename Derived1::Index Index;
+ typedef typename Derived1::Scalar Scalar;
static inline void run(Derived1 &dst, const Derived2 &src)
{
for(Index j = 0; j < dst.cols(); ++j)
@@ -520,7 +521,7 @@
dst.copyCoeff(i, j, src);
if (ClearOpposite)
for(Index i = maxi; i < dst.rows(); ++i)
- dst.coeffRef(i, j) = 0;
+ dst.coeffRef(i, j) = Scalar(0);
}
}
};
Modified: pkg/RcppEigen/inst/include/Eigen/src/Core/products/GeneralMatrixVector.h
===================================================================
--- pkg/RcppEigen/inst/include/Eigen/src/Core/products/GeneralMatrixVector.h 2012-11-29 15:32:15 UTC (rev 4051)
+++ pkg/RcppEigen/inst/include/Eigen/src/Core/products/GeneralMatrixVector.h 2012-11-29 18:00:37 UTC (rev 4052)
@@ -81,14 +81,13 @@
const Index peels = 2;
const Index LhsPacketAlignedMask = LhsPacketSize-1;
const Index ResPacketAlignedMask = ResPacketSize-1;
- const Index PeelAlignedMask = ResPacketSize*peels-1;
const Index size = rows;
// How many coeffs of the result do we have to skip to be aligned.
// Here we assume data are at least aligned on the base scalar type.
Index alignedStart = internal::first_aligned(res,size);
Index alignedSize = ResPacketSize>1 ? alignedStart + ((size-alignedStart) & ~ResPacketAlignedMask) : 0;
- const Index peeledSize = peels>1 ? alignedStart + ((alignedSize-alignedStart) & ~PeelAlignedMask) : alignedStart;
+ const Index peeledSize = alignedSize - RhsPacketSize*peels - RhsPacketSize + 1;
const Index alignmentStep = LhsPacketSize>1 ? (LhsPacketSize - lhsStride % LhsPacketSize) & LhsPacketAlignedMask : 0;
Index alignmentPattern = alignmentStep==0 ? AllAligned
@@ -177,6 +176,8 @@
_EIGEN_ACCUMULATE_PACKETS(d,du,d);
break;
case FirstAligned:
+ {
+ Index j = alignedStart;
if(peels>1)
{
LhsPacket A00, A01, A02, A03, A10, A11, A12, A13;
@@ -186,7 +187,7 @@
A02 = pload<LhsPacket>(&lhs2[alignedStart-2]);
A03 = pload<LhsPacket>(&lhs3[alignedStart-3]);
- for (Index j = alignedStart; j<peeledSize; j+=peels*ResPacketSize)
+ for (; j<peeledSize; j+=peels*ResPacketSize)
{
A11 = pload<LhsPacket>(&lhs1[j-1+LhsPacketSize]); palign<1>(A01,A11);
A12 = pload<LhsPacket>(&lhs2[j-2+LhsPacketSize]); palign<2>(A02,A12);
@@ -210,9 +211,10 @@
pstore(&res[j+ResPacketSize],T1);
}
}
- for (Index j = peeledSize; j<alignedSize; j+=ResPacketSize)
+ for (; j<alignedSize; j+=ResPacketSize)
_EIGEN_ACCUMULATE_PACKETS(d,du,du);
break;
+ }
default:
for (Index j = alignedStart; j<alignedSize; j+=ResPacketSize)
_EIGEN_ACCUMULATE_PACKETS(du,du,du);
@@ -332,7 +334,6 @@
const Index peels = 2;
const Index RhsPacketAlignedMask = RhsPacketSize-1;
const Index LhsPacketAlignedMask = LhsPacketSize-1;
- const Index PeelAlignedMask = RhsPacketSize*peels-1;
const Index depth = cols;
// How many coeffs of the result do we have to skip to be aligned.
@@ -340,7 +341,7 @@
// if that's not the case then vectorization is discarded, see below.
Index alignedStart = internal::first_aligned(rhs, depth);
Index alignedSize = RhsPacketSize>1 ? alignedStart + ((depth-alignedStart) & ~RhsPacketAlignedMask) : 0;
- const Index peeledSize = peels>1 ? alignedStart + ((alignedSize-alignedStart) & ~PeelAlignedMask) : alignedStart;
+ const Index peeledSize = alignedSize - RhsPacketSize*peels - RhsPacketSize + 1;
const Index alignmentStep = LhsPacketSize>1 ? (LhsPacketSize - lhsStride % LhsPacketSize) & LhsPacketAlignedMask : 0;
Index alignmentPattern = alignmentStep==0 ? AllAligned
@@ -430,10 +431,12 @@
_EIGEN_ACCUMULATE_PACKETS(d,du,d);
break;
case FirstAligned:
+ {
+ Index j = alignedStart;
if (peels>1)
{
/* Here we proccess 4 rows with with two peeled iterations to hide
- * tghe overhead of unaligned loads. Moreover unaligned loads are handled
+ * the overhead of unaligned loads. Moreover unaligned loads are handled
* using special shift/move operations between the two aligned packets
* overlaping the desired unaligned packet. This is *much* more efficient
* than basic unaligned loads.
@@ -443,7 +446,7 @@
A02 = pload<LhsPacket>(&lhs2[alignedStart-2]);
A03 = pload<LhsPacket>(&lhs3[alignedStart-3]);
- for (Index j = alignedStart; j<peeledSize; j+=peels*RhsPacketSize)
+ for (; j<peeledSize; j+=peels*RhsPacketSize)
{
RhsPacket b = pload<RhsPacket>(&rhs[j]);
A11 = pload<LhsPacket>(&lhs1[j-1+LhsPacketSize]); palign<1>(A01,A11);
@@ -465,9 +468,10 @@
ptmp3 = pcj.pmadd(A13, b, ptmp3);
}
}
- for (Index j = peeledSize; j<alignedSize; j+=RhsPacketSize)
+ for (; j<alignedSize; j+=RhsPacketSize)
_EIGEN_ACCUMULATE_PACKETS(d,du,du);
break;
+ }
default:
for (Index j = alignedStart; j<alignedSize; j+=RhsPacketSize)
_EIGEN_ACCUMULATE_PACKETS(du,du,du);
Modified: pkg/RcppEigen/inst/include/Eigen/src/Core/products/TriangularMatrixMatrix_MKL.h
===================================================================
--- pkg/RcppEigen/inst/include/Eigen/src/Core/products/TriangularMatrixMatrix_MKL.h 2012-11-29 15:32:15 UTC (rev 4051)
+++ pkg/RcppEigen/inst/include/Eigen/src/Core/products/TriangularMatrixMatrix_MKL.h 2012-11-29 18:00:37 UTC (rev 4052)
@@ -57,11 +57,11 @@
struct product_triangular_matrix_matrix<Scalar,Index, Mode, LhsIsTriangular, \
LhsStorageOrder,ConjugateLhs, RhsStorageOrder,ConjugateRhs,ColMajor,Specialized> { \
static inline void run(Index _rows, Index _cols, Index _depth, const Scalar* _lhs, Index lhsStride,\
- const Scalar* _rhs, Index rhsStride, Scalar* res, Index resStride, Scalar alpha) { \
+ const Scalar* _rhs, Index rhsStride, Scalar* res, Index resStride, Scalar alpha, level3_blocking<Scalar,Scalar>& blocking) { \
product_triangular_matrix_matrix_trmm<Scalar,Index,Mode, \
LhsIsTriangular,LhsStorageOrder,ConjugateLhs, \
RhsStorageOrder, ConjugateRhs, ColMajor>::run( \
- _rows, _cols, _depth, _lhs, lhsStride, _rhs, rhsStride, res, resStride, alpha); \
+ _rows, _cols, _depth, _lhs, lhsStride, _rhs, rhsStride, res, resStride, alpha, blocking); \
} \
};
@@ -96,7 +96,7 @@
const EIGTYPE* _lhs, Index lhsStride, \
const EIGTYPE* _rhs, Index rhsStride, \
EIGTYPE* res, Index resStride, \
- EIGTYPE alpha) \
+ EIGTYPE alpha, level3_blocking<EIGTYPE,EIGTYPE>& blocking) \
{ \
Index diagSize = (std::min)(_rows,_depth); \
Index rows = IsLower ? _rows : diagSize; \
@@ -115,16 +115,16 @@
/* Most likely no benefit to call TRMM or GEMM from MKL*/ \
product_triangular_matrix_matrix<EIGTYPE,Index,Mode,true, \
LhsStorageOrder,ConjugateLhs, RhsStorageOrder, ConjugateRhs, ColMajor, BuiltIn>::run( \
- _rows, _cols, _depth, _lhs, lhsStride, _rhs, rhsStride, res, resStride, alpha); \
+ _rows, _cols, _depth, _lhs, lhsStride, _rhs, rhsStride, res, resStride, alpha, blocking); \
/*std::cout << "TRMM_L: A is not square! Go to Eigen TRMM implementation!\n";*/ \
} else { \
/* Make sense to call GEMM */ \
Map<const MatrixLhs, 0, OuterStride<> > lhsMap(_lhs,rows,depth,OuterStride<>(lhsStride)); \
MatrixLhs aa_tmp=lhsMap.template triangularView<Mode>(); \
MKL_INT aStride = aa_tmp.outerStride(); \
- gemm_blocking_space<ColMajor,EIGTYPE,EIGTYPE,Dynamic,Dynamic,Dynamic> blocking(_rows,_cols,_depth); \
+ gemm_blocking_space<ColMajor,EIGTYPE,EIGTYPE,Dynamic,Dynamic,Dynamic> gemm_blocking(_rows,_cols,_depth); \
general_matrix_matrix_product<Index,EIGTYPE,LhsStorageOrder,ConjugateLhs,EIGTYPE,RhsStorageOrder,ConjugateRhs,ColMajor>::run( \
- rows, cols, depth, aa_tmp.data(), aStride, _rhs, rhsStride, res, resStride, alpha, blocking, 0); \
+ rows, cols, depth, aa_tmp.data(), aStride, _rhs, rhsStride, res, resStride, alpha, gemm_blocking, 0); \
\
/*std::cout << "TRMM_L: A is not square! Go to MKL GEMM implementation! " << nthr<<" \n";*/ \
} \
@@ -210,7 +210,7 @@
const EIGTYPE* _lhs, Index lhsStride, \
const EIGTYPE* _rhs, Index rhsStride, \
EIGTYPE* res, Index resStride, \
- EIGTYPE alpha) \
+ EIGTYPE alpha, level3_blocking<EIGTYPE,EIGTYPE>& blocking) \
{ \
Index diagSize = (std::min)(_cols,_depth); \
Index rows = _rows; \
@@ -229,16 +229,16 @@
/* Most likely no benefit to call TRMM or GEMM from MKL*/ \
product_triangular_matrix_matrix<EIGTYPE,Index,Mode,false, \
LhsStorageOrder,ConjugateLhs, RhsStorageOrder, ConjugateRhs, ColMajor, BuiltIn>::run( \
- _rows, _cols, _depth, _lhs, lhsStride, _rhs, rhsStride, res, resStride, alpha); \
+ _rows, _cols, _depth, _lhs, lhsStride, _rhs, rhsStride, res, resStride, alpha, blocking); \
/*std::cout << "TRMM_R: A is not square! Go to Eigen TRMM implementation!\n";*/ \
} else { \
/* Make sense to call GEMM */ \
Map<const MatrixRhs, 0, OuterStride<> > rhsMap(_rhs,depth,cols, OuterStride<>(rhsStride)); \
MatrixRhs aa_tmp=rhsMap.template triangularView<Mode>(); \
MKL_INT aStride = aa_tmp.outerStride(); \
- gemm_blocking_space<ColMajor,EIGTYPE,EIGTYPE,Dynamic,Dynamic,Dynamic> blocking(_rows,_cols,_depth); \
+ gemm_blocking_space<ColMajor,EIGTYPE,EIGTYPE,Dynamic,Dynamic,Dynamic> gemm_blocking(_rows,_cols,_depth); \
general_matrix_matrix_product<Index,EIGTYPE,LhsStorageOrder,ConjugateLhs,EIGTYPE,RhsStorageOrder,ConjugateRhs,ColMajor>::run( \
- rows, cols, depth, _lhs, lhsStride, aa_tmp.data(), aStride, res, resStride, alpha, blocking, 0); \
+ rows, cols, depth, _lhs, lhsStride, aa_tmp.data(), aStride, res, resStride, alpha, gemm_blocking, 0); \
\
/*std::cout << "TRMM_R: A is not square! Go to MKL GEMM implementation! " << nthr<<" \n";*/ \
} \
Modified: pkg/RcppEigen/inst/include/Eigen/src/Core/products/TriangularMatrixVector_MKL.h
===================================================================
--- pkg/RcppEigen/inst/include/Eigen/src/Core/products/TriangularMatrixVector_MKL.h 2012-11-29 15:32:15 UTC (rev 4051)
+++ pkg/RcppEigen/inst/include/Eigen/src/Core/products/TriangularMatrixVector_MKL.h 2012-11-29 18:00:37 UTC (rev 4052)
@@ -82,11 +82,11 @@
LowUp = IsLower ? Lower : Upper \
}; \
static EIGEN_DONT_INLINE void run(Index _rows, Index _cols, const EIGTYPE* _lhs, Index lhsStride, \
- const EIGTYPE* _rhs, Index rhsIncr, EIGTYPE* _res, Index resIncr, EIGTYPE alpha, level3_blocking<EIGTYPE,EIGTYPE>& blocking) \
+ const EIGTYPE* _rhs, Index rhsIncr, EIGTYPE* _res, Index resIncr, EIGTYPE alpha) \
{ \
if (ConjLhs || IsZeroDiag) { \
triangular_matrix_vector_product<Index,Mode,EIGTYPE,ConjLhs,EIGTYPE,ConjRhs,ColMajor,BuiltIn>::run( \
- _rows, _cols, _lhs, lhsStride, _rhs, rhsIncr, _res, resIncr, alpha, blocking); \
+ _rows, _cols, _lhs, lhsStride, _rhs, rhsIncr, _res, resIncr, alpha); \
return; \
}\
Index size = (std::min)(_rows,_cols); \
@@ -167,11 +167,11 @@
LowUp = IsLower ? Lower : Upper \
}; \
static EIGEN_DONT_INLINE void run(Index _rows, Index _cols, const EIGTYPE* _lhs, Index lhsStride, \
- const EIGTYPE* _rhs, Index rhsIncr, EIGTYPE* _res, Index resIncr, EIGTYPE alpha, level3_blocking<EIGTYPE,EIGTYPE>& blocking) \
+ const EIGTYPE* _rhs, Index rhsIncr, EIGTYPE* _res, Index resIncr, EIGTYPE alpha) \
{ \
if (IsZeroDiag) { \
triangular_matrix_vector_product<Index,Mode,EIGTYPE,ConjLhs,EIGTYPE,ConjRhs,RowMajor,BuiltIn>::run( \
- _rows, _cols, _lhs, lhsStride, _rhs, rhsIncr, _res, resIncr, alpha, blocking); \
+ _rows, _cols, _lhs, lhsStride, _rhs, rhsIncr, _res, resIncr, alpha); \
return; \
}\
Index size = (std::min)(_rows,_cols); \
Modified: pkg/RcppEigen/inst/include/Eigen/src/Core/util/Macros.h
===================================================================
--- pkg/RcppEigen/inst/include/Eigen/src/Core/util/Macros.h 2012-11-29 15:32:15 UTC (rev 4051)
+++ pkg/RcppEigen/inst/include/Eigen/src/Core/util/Macros.h 2012-11-29 18:00:37 UTC (rev 4052)
@@ -13,7 +13,7 @@
#define EIGEN_WORLD_VERSION 3
#define EIGEN_MAJOR_VERSION 1
-#define EIGEN_MINOR_VERSION 1
+#define EIGEN_MINOR_VERSION 2
#define EIGEN_VERSION_AT_LEAST(x,y,z) (EIGEN_WORLD_VERSION>x || (EIGEN_WORLD_VERSION>=x && \
(EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \
Modified: pkg/RcppEigen/inst/include/Eigen/src/Core/util/Memory.h
===================================================================
--- pkg/RcppEigen/inst/include/Eigen/src/Core/util/Memory.h 2012-11-29 15:32:15 UTC (rev 4051)
+++ pkg/RcppEigen/inst/include/Eigen/src/Core/util/Memory.h 2012-11-29 18:00:37 UTC (rev 4052)
@@ -204,7 +204,7 @@
if(posix_memalign(&result, 16, size)) result = 0;
#elif EIGEN_HAS_MM_MALLOC
result = _mm_malloc(size, 16);
- #elif (defined _MSC_VER)
+#elif defined(_MSC_VER) && (!defined(_WIN32_WCE))
result = _aligned_malloc(size, 16);
#else
result = handmade_aligned_malloc(size);
@@ -745,7 +745,7 @@
__asm__ __volatile__ ("cpuid": "=a" (abcd[0]), "=b" (abcd[1]), "=c" (abcd[2]), "=d" (abcd[3]) : "a" (func), "c" (id) );
# endif
# elif defined(_MSC_VER)
-# if (_MSC_VER > 1500)
+# if (_MSC_VER > 1500) && ( defined(_M_IX86) || defined(_M_X64) )
# define EIGEN_CPUID(abcd,func,id) __cpuidex((int*)abcd,func,id)
# endif
# endif
Modified: pkg/RcppEigen/inst/include/Eigen/src/Core/util/XprHelper.h
===================================================================
--- pkg/RcppEigen/inst/include/Eigen/src/Core/util/XprHelper.h 2012-11-29 15:32:15 UTC (rev 4051)
+++ pkg/RcppEigen/inst/include/Eigen/src/Core/util/XprHelper.h 2012-11-29 18:00:37 UTC (rev 4052)
@@ -301,9 +301,9 @@
// it's important that this value can still be squared without integer overflowing.
DynamicAsInteger = 10000,
ScalarReadCost = NumTraits<typename traits<T>::Scalar>::ReadCost,
- ScalarReadCostAsInteger = ScalarReadCost == Dynamic ? DynamicAsInteger : ScalarReadCost,
+ ScalarReadCostAsInteger = ScalarReadCost == Dynamic ? int(DynamicAsInteger) : int(ScalarReadCost),
CoeffReadCost = traits<T>::CoeffReadCost,
- CoeffReadCostAsInteger = CoeffReadCost == Dynamic ? DynamicAsInteger : CoeffReadCost,
+ CoeffReadCostAsInteger = CoeffReadCost == Dynamic ? int(DynamicAsInteger) : int(CoeffReadCost),
NAsInteger = n == Dynamic ? int(DynamicAsInteger) : n,
CostEvalAsInteger = (NAsInteger+1) * ScalarReadCostAsInteger + CoeffReadCostAsInteger,
CostNoEvalAsInteger = NAsInteger * CoeffReadCostAsInteger
Modified: pkg/RcppEigen/inst/include/Eigen/src/Eigenvalues/ComplexSchur.h
===================================================================
--- pkg/RcppEigen/inst/include/Eigen/src/Eigenvalues/ComplexSchur.h 2012-11-29 15:32:15 UTC (rev 4051)
+++ pkg/RcppEigen/inst/include/Eigen/src/Eigenvalues/ComplexSchur.h 2012-11-29 18:00:37 UTC (rev 4052)
@@ -336,6 +336,7 @@
Index iu = m_matT.cols() - 1;
Index il;
Index iter = 0; // number of iterations we are working on the (iu,iu) element
+ Index totalIter = 0; // number of iterations for whole matrix
while(true)
{
@@ -350,9 +351,10 @@
// if iu is zero then we are done; the whole matrix is triangularized
if(iu==0) break;
- // if we spent too many iterations on the current element, we give up
+ // if we spent too many iterations, we give up
iter++;
- if(iter > m_maxIterations * m_matT.cols()) break;
+ totalIter++;
+ if(totalIter > m_maxIterations * m_matT.cols()) break;
// find il, the top row of the active submatrix
il = iu-1;
@@ -382,7 +384,7 @@
}
}
- if(iter <= m_maxIterations * m_matT.cols())
+ if(totalIter <= m_maxIterations * m_matT.cols())
m_info = Success;
else
m_info = NoConvergence;
Modified: pkg/RcppEigen/inst/include/Eigen/src/Eigenvalues/ComplexSchur_MKL.h
===================================================================
--- pkg/RcppEigen/inst/include/Eigen/src/Eigenvalues/ComplexSchur_MKL.h 2012-11-29 15:32:15 UTC (rev 4051)
+++ pkg/RcppEigen/inst/include/Eigen/src/Eigenvalues/ComplexSchur_MKL.h 2012-11-29 18:00:37 UTC (rev 4052)
@@ -40,7 +40,7 @@
/** \internal Specialization for the data types supported by MKL */
#define EIGEN_MKL_SCHUR_COMPLEX(EIGTYPE, MKLTYPE, MKLPREFIX, MKLPREFIX_U, EIGCOLROW, MKLCOLROW) \
-template<> inline\
+template<> inline \
ComplexSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >& \
ComplexSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(const Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW>& matrix, bool computeU) \
{ \
Deleted: pkg/RcppEigen/inst/include/Eigen/src/Eigenvalues/EigenvaluesCommon.h
===================================================================
--- pkg/RcppEigen/inst/include/Eigen/src/Eigenvalues/EigenvaluesCommon.h 2012-11-29 15:32:15 UTC (rev 4051)
+++ pkg/RcppEigen/inst/include/Eigen/src/Eigenvalues/EigenvaluesCommon.h 2012-11-29 18:00:37 UTC (rev 4052)
@@ -1,31 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2010 Jitse Niesen <jitse at maths.leeds.ac.uk>
-//
-// 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_EIGENVALUES_COMMON_H
-#define EIGEN_EIGENVALUES_COMMON_H
-
-
-
-#endif // EIGEN_EIGENVALUES_COMMON_H
-
Modified: pkg/RcppEigen/inst/include/Eigen/src/Eigenvalues/RealSchur.h
===================================================================
--- pkg/RcppEigen/inst/include/Eigen/src/Eigenvalues/RealSchur.h 2012-11-29 15:32:15 UTC (rev 4051)
+++ pkg/RcppEigen/inst/include/Eigen/src/Eigenvalues/RealSchur.h 2012-11-29 18:00:37 UTC (rev 4052)
@@ -220,8 +220,9 @@
// Rows il,...,iu is the part we are working on (the active window).
// Rows iu+1,...,end are already brought in triangular form.
Index iu = m_matT.cols() - 1;
- Index iter = 0; // iteration count
- Scalar exshift(0); // sum of exceptional shifts
+ Index iter = 0; // iteration count for current eigenvalue
+ Index totalIter = 0; // iteration count for whole matrix
+ Scalar exshift(0); // sum of exceptional shifts
Scalar norm = computeNormOfT();
if(norm!=0)
@@ -251,14 +252,15 @@
Vector3s firstHouseholderVector(0,0,0), shiftInfo;
computeShift(iu, iter, exshift, shiftInfo);
iter = iter + 1;
- if (iter > m_maxIterations * m_matT.cols()) break;
+ totalIter = totalIter + 1;
+ if (totalIter > m_maxIterations * matrix.cols()) break;
Index im;
initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
}
}
}
- if(iter <= m_maxIterations * m_matT.cols())
+ if(totalIter <= m_maxIterations * matrix.cols())
m_info = Success;
else
m_info = NoConvergence;
Modified: pkg/RcppEigen/inst/include/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
===================================================================
--- pkg/RcppEigen/inst/include/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h 2012-11-29 15:32:15 UTC (rev 4051)
+++ pkg/RcppEigen/inst/include/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h 2012-11-29 18:00:37 UTC (rev 4052)
@@ -743,7 +743,16 @@
// RealScalar e2 = abs2(subdiag[end-1]);
// RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2));
// This explain the following, somewhat more complicated, version:
- RealScalar mu = diag[end] - (e / (td + (td>0 ? 1 : -1))) * (e / hypot(td,e));
+ RealScalar mu = diag[end];
+ if(td==0)
+ mu -= abs(e);
+ else
+ {
+ RealScalar e2 = abs2(subdiag[end-1]);
+ RealScalar h = hypot(td,e);
+ if(e2==0) mu -= (e / (td + (td>0 ? 1 : -1))) * (e / h);
+ else mu -= e2 / (td + (td>0 ? h : -h));
+ }
RealScalar x = diag[start] - mu;
RealScalar z = subdiag[start];
Modified: pkg/RcppEigen/inst/include/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h
===================================================================
--- pkg/RcppEigen/inst/include/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h 2012-11-29 15:32:15 UTC (rev 4051)
+++ pkg/RcppEigen/inst/include/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h 2012-11-29 18:00:37 UTC (rev 4052)
@@ -40,7 +40,7 @@
/** \internal Specialization for the data types supported by MKL */
#define EIGEN_MKL_EIG_SELFADJ(EIGTYPE, MKLTYPE, MKLRTYPE, MKLNAME, EIGCOLROW, MKLCOLROW ) \
-template<> inline\
+template<> inline \
SelfAdjointEigenSolver<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >& \
SelfAdjointEigenSolver<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(const Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW>& matrix, int options) \
{ \
Modified: pkg/RcppEigen/inst/include/Eigen/src/Geometry/AlignedBox.h
===================================================================
--- pkg/RcppEigen/inst/include/Eigen/src/Geometry/AlignedBox.h 2012-11-29 15:32:15 UTC (rev 4051)
+++ pkg/RcppEigen/inst/include/Eigen/src/Geometry/AlignedBox.h 2012-11-29 18:00:37 UTC (rev 4052)
@@ -79,7 +79,7 @@
~AlignedBox() {}
/** \returns the dimension in which the box holds */
- inline Index dim() const { return AmbientDimAtCompileTime==Dynamic ? m_min.size()-1 : Index(AmbientDimAtCompileTime); }
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rcpp -r 4052
More information about the Rcpp-commits
mailing list