diff options
Diffstat (limited to 'Eigen/src')
28 files changed, 328 insertions, 214 deletions
diff --git a/Eigen/src/Core/BooleanRedux.h b/Eigen/src/Core/BooleanRedux.h index f6afeb034..6e37e031a 100644 --- a/Eigen/src/Core/BooleanRedux.h +++ b/Eigen/src/Core/BooleanRedux.h @@ -131,7 +131,7 @@ inline typename DenseBase<Derived>::Index DenseBase<Derived>::count() const /** \returns true is \c *this contains at least one Not A Number (NaN). * - * \sa isFinite() + * \sa allFinite() */ template<typename Derived> inline bool DenseBase<Derived>::hasNaN() const @@ -144,7 +144,7 @@ inline bool DenseBase<Derived>::hasNaN() const * \sa hasNaN() */ template<typename Derived> -inline bool DenseBase<Derived>::isFinite() const +inline bool DenseBase<Derived>::allFinite() const { return !((derived()-derived()).hasNaN()); } diff --git a/Eigen/src/Core/DenseBase.h b/Eigen/src/Core/DenseBase.h index 4e8b820bb..c5800f6c8 100644 --- a/Eigen/src/Core/DenseBase.h +++ b/Eigen/src/Core/DenseBase.h @@ -281,7 +281,7 @@ template<typename Derived> class DenseBase CommaInitializer<Derived> operator<< (const DenseBase<OtherDerived>& other); Eigen::Transpose<Derived> transpose(); - typedef const Transpose<const Derived> ConstTransposeReturnType; + typedef typename internal::add_const<Transpose<const Derived> >::type ConstTransposeReturnType; ConstTransposeReturnType transpose() const; void transposeInPlace(); #ifndef EIGEN_NO_DEBUG @@ -348,7 +348,7 @@ template<typename Derived> class DenseBase bool isOnes(const RealScalar& prec = NumTraits<Scalar>::dummy_precision()) const; inline bool hasNaN() const; - inline bool isFinite() const; + inline bool allFinite() const; inline Derived& operator*=(const Scalar& other); inline Derived& operator/=(const Scalar& other); diff --git a/Eigen/src/Core/Diagonal.h b/Eigen/src/Core/Diagonal.h index c106f93c4..aab8007b3 100644 --- a/Eigen/src/Core/Diagonal.h +++ b/Eigen/src/Core/Diagonal.h @@ -172,7 +172,7 @@ MatrixBase<Derived>::diagonal() /** This is the const version of diagonal(). */ template<typename Derived> -inline const typename MatrixBase<Derived>::ConstDiagonalReturnType +inline typename MatrixBase<Derived>::ConstDiagonalReturnType MatrixBase<Derived>::diagonal() const { return ConstDiagonalReturnType(derived()); diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index d415bc719..454f9ce52 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h @@ -518,11 +518,10 @@ struct random_default_impl<Scalar, false, true> #else enum { rand_bits = floor_log2<(unsigned int)(RAND_MAX)+1>::value, scalar_bits = sizeof(Scalar) * CHAR_BIT, - shift = EIGEN_PLAIN_ENUM_MAX(0, int(rand_bits) - int(scalar_bits)) + shift = EIGEN_PLAIN_ENUM_MAX(0, int(rand_bits) - int(scalar_bits)), + offset = NumTraits<Scalar>::IsSigned ? (1 << (EIGEN_PLAIN_ENUM_MIN(rand_bits,scalar_bits)-1)) : 0 }; - Scalar x = Scalar(std::rand() >> shift); - Scalar offset = NumTraits<Scalar>::IsSigned ? Scalar(1 << (rand_bits-1)) : Scalar(0); - return x - offset; + return Scalar((std::rand() >> shift) - offset); #endif } }; diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 0628ebd1f..fbed47233 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -212,8 +212,8 @@ template<typename Derived> class MatrixBase typedef Diagonal<Derived> DiagonalReturnType; DiagonalReturnType diagonal(); - typedef const Diagonal<const Derived> ConstDiagonalReturnType; - const ConstDiagonalReturnType diagonal() const; + typedef typename internal::add_const<Diagonal<const Derived> >::type ConstDiagonalReturnType; + ConstDiagonalReturnType diagonal() const; template<int Index> struct DiagonalIndexReturnType { typedef Diagonal<Derived,Index> Type; }; template<int Index> struct ConstDiagonalIndexReturnType { typedef const Diagonal<const Derived,Index> Type; }; diff --git a/Eigen/src/Core/ReturnByValue.h b/Eigen/src/Core/ReturnByValue.h index 613912ffa..d66c24ba0 100644 --- a/Eigen/src/Core/ReturnByValue.h +++ b/Eigen/src/Core/ReturnByValue.h @@ -48,7 +48,7 @@ struct nested<ReturnByValue<Derived>, n, PlainObject> } // end namespace internal template<typename Derived> class ReturnByValue - : public internal::dense_xpr_base< ReturnByValue<Derived> >::type + : internal::no_assignment_operator, public internal::dense_xpr_base< ReturnByValue<Derived> >::type { public: typedef typename internal::traits<Derived>::ReturnType ReturnType; diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h index 798120bf4..f21b3aa65 100644 --- a/Eigen/src/Core/Transpose.h +++ b/Eigen/src/Core/Transpose.h @@ -207,7 +207,7 @@ DenseBase<Derived>::transpose() * * \sa transposeInPlace(), adjoint() */ template<typename Derived> -inline const typename DenseBase<Derived>::ConstTransposeReturnType +inline typename DenseBase<Derived>::ConstTransposeReturnType DenseBase<Derived>::transpose() const { return ConstTransposeReturnType(derived()); diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h index 96737456a..6798a3e0f 100644 --- a/Eigen/src/Core/util/Macros.h +++ b/Eigen/src/Core/util/Macros.h @@ -12,8 +12,8 @@ #define EIGEN_MACROS_H #define EIGEN_WORLD_VERSION 3 -#define EIGEN_MAJOR_VERSION 1 -#define EIGEN_MINOR_VERSION 91 +#define EIGEN_MAJOR_VERSION 2 +#define EIGEN_MINOR_VERSION 90 #define EIGEN_VERSION_AT_LEAST(x,y,z) (EIGEN_WORLD_VERSION>x || (EIGEN_WORLD_VERSION>=x && \ (EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \ diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 3993046a8..4e06809c4 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -20,6 +20,8 @@ class GeneralizedSelfAdjointEigenSolver; namespace internal { template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues; +template<typename MatrixType, typename DiagType, typename SubDiagType> +ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag, const typename MatrixType::Index maxIterations, bool computeEigenvectors, MatrixType& eivec); } /** \eigenvalues_module \ingroup Eigenvalues_Module @@ -98,6 +100,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver */ typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType; typedef Tridiagonalization<MatrixType> TridiagonalizationType; + typedef typename TridiagonalizationType::SubDiagonalType SubDiagonalType; /** \brief Default constructor for fixed-size matrices. * @@ -207,6 +210,19 @@ template<typename _MatrixType> class SelfAdjointEigenSolver */ SelfAdjointEigenSolver& computeDirect(const MatrixType& matrix, int options = ComputeEigenvectors); + /** + *\brief Computes the eigen decomposition from a tridiagonal symmetric matrix + * + * \param[in] diag The vector containing the diagonal of the matrix. + * \param[in] subdiag The subdiagonal of the matrix. + * \returns Reference to \c *this + * + * This function assumes that the matrix has been reduced to tridiagonal form. + * + * \sa compute(const MatrixType&, int) for more information + */ + SelfAdjointEigenSolver& computeFromTridiagonal(const RealVectorType& diag, const SubDiagonalType& subdiag , int options=ComputeEigenvectors); + /** \brief Returns the eigenvectors of given matrix. * * \returns A const reference to the matrix whose columns are the eigenvectors. @@ -415,7 +431,58 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> mat.template triangularView<Lower>() /= scale; m_subdiag.resize(n-1); internal::tridiagonalization_inplace(mat, diag, m_subdiag, computeEigenvectors); + + m_info = internal::computeFromTridiagonal_impl(diag, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec); + // scale back the eigen values + m_eivalues *= scale; + + m_isInitialized = true; + m_eigenvectorsOk = computeEigenvectors; + return *this; +} + +template<typename MatrixType> +SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> +::computeFromTridiagonal(const RealVectorType& diag, const SubDiagonalType& subdiag , int options) +{ + //TODO : Add an option to scale the values beforehand + bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors; + + m_eivalues = diag; + m_subdiag = subdiag; + if (computeEigenvectors) + { + m_eivec.setIdentity(diag.size(), diag.size()); + } + m_info = computeFromTridiagonal_impl(m_eivalues, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec); + + m_isInitialized = true; + m_eigenvectorsOk = computeEigenvectors; + return *this; +} + +namespace internal { +/** + * \internal + * \brief Compute the eigendecomposition from a tridiagonal matrix + * + * \param[in,out] diag : On input, the diagonal of the matrix, on output the eigenvalues + * \param[in] subdiag : The subdiagonal part of the matrix. + * \param[in,out] : On input, the maximum number of iterations, on output, the effective number of iterations. + * \param[out] eivec : The matrix to store the eigenvectors... if needed. allocated on input + * \returns \c Success or \c NoConvergence + */ +template<typename MatrixType, typename DiagType, typename SubDiagType> +ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag, const typename MatrixType::Index maxIterations, bool computeEigenvectors, MatrixType& eivec) +{ + using std::abs; + + ComputationInfo info; + typedef typename MatrixType::Index Index; + typedef typename MatrixType::Scalar Scalar; + + Index n = diag.size(); Index end = n-1; Index start = 0; Index iter = 0; // total number of iterations @@ -423,11 +490,11 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> while (end>0) { for (Index i = start; i<end; ++i) - if (internal::isMuchSmallerThan(abs(m_subdiag[i]),(abs(diag[i])+abs(diag[i+1])))) - m_subdiag[i] = 0; + if (internal::isMuchSmallerThan(abs(subdiag[i]),(abs(diag[i])+abs(diag[i+1])))) + subdiag[i] = 0; // find the largest unreduced block - while (end>0 && m_subdiag[end-1]==0) + while (end>0 && subdiag[end-1]==0) { end--; } @@ -436,48 +503,38 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> // if we spent too many iterations, we give up iter++; - if(iter > m_maxIterations * n) break; + if(iter > maxIterations * n) break; start = end - 1; - while (start>0 && m_subdiag[start-1]!=0) + while (start>0 && subdiag[start-1]!=0) start--; - internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n); + internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), subdiag.data(), start, end, computeEigenvectors ? eivec.data() : (Scalar*)0, n); } - - if (iter <= m_maxIterations * n) - m_info = Success; + if (iter <= maxIterations * n) + info = Success; else - m_info = NoConvergence; + info = NoConvergence; // Sort eigenvalues and corresponding vectors. // TODO make the sort optional ? // TODO use a better sort algorithm !! - if (m_info == Success) + if (info == Success) { for (Index i = 0; i < n-1; ++i) { Index k; - m_eivalues.segment(i,n-i).minCoeff(&k); + diag.segment(i,n-i).minCoeff(&k); if (k > 0) { - std::swap(m_eivalues[i], m_eivalues[k+i]); + std::swap(diag[i], diag[k+i]); if(computeEigenvectors) - m_eivec.col(i).swap(m_eivec.col(k+i)); + eivec.col(i).swap(eivec.col(k+i)); } } } - - // scale back the eigen values - m_eivalues *= scale; - - m_isInitialized = true; - m_eigenvectorsOk = computeEigenvectors; - return *this; + return info; } - - -namespace internal { template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues { diff --git a/Eigen/src/Geometry/Homogeneous.h b/Eigen/src/Geometry/Homogeneous.h index df03feb55..00e71d190 100644 --- a/Eigen/src/Geometry/Homogeneous.h +++ b/Eigen/src/Geometry/Homogeneous.h @@ -59,7 +59,7 @@ template<typename MatrixType,typename Rhs> struct homogeneous_right_product_impl } // end namespace internal template<typename MatrixType,int _Direction> class Homogeneous - : public MatrixBase<Homogeneous<MatrixType,_Direction> > + : internal::no_assignment_operator, public MatrixBase<Homogeneous<MatrixType,_Direction> > { public: diff --git a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h index 50a870aec..b55afc136 100644 --- a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h +++ b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h @@ -150,8 +150,7 @@ class IncompleteLUT : internal::noncopyable { analyzePattern(amat); factorize(amat); - eigen_assert(m_factorizationIsOk == true); - m_isInitialized = true; + m_isInitialized = m_factorizationIsOk; return *this; } diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h index 8a7c9b9e2..8b01f8179 100644 --- a/Eigen/src/QR/ColPivHouseholderQR.h +++ b/Eigen/src/QR/ColPivHouseholderQR.h @@ -525,8 +525,8 @@ struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs> { eigen_assert(rhs().rows() == dec().rows()); - const int cols = dec().cols(), - nonzero_pivots = dec().nonzeroPivots(); + const Index cols = dec().cols(), + nonzero_pivots = dec().nonzeroPivots(); if(nonzero_pivots == 0) { diff --git a/Eigen/src/SparseCore/SparseBlock.h b/Eigen/src/SparseCore/SparseBlock.h index e025e4d40..0b3e193db 100644 --- a/Eigen/src/SparseCore/SparseBlock.h +++ b/Eigen/src/SparseCore/SparseBlock.h @@ -27,6 +27,7 @@ public: class InnerIterator: public XprType::InnerIterator { + typedef typename BlockImpl::Index Index; public: inline InnerIterator(const BlockType& xpr, Index outer) : XprType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) @@ -38,6 +39,7 @@ public: }; class ReverseInnerIterator: public XprType::ReverseInnerIterator { + typedef typename BlockImpl::Index Index; public: inline ReverseInnerIterator(const BlockType& xpr, Index outer) : XprType::ReverseInnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) diff --git a/Eigen/src/SparseCore/SparseCwiseBinaryOp.h b/Eigen/src/SparseCore/SparseCwiseBinaryOp.h index 64b8c8547..ec86ca933 100644 --- a/Eigen/src/SparseCore/SparseCwiseBinaryOp.h +++ b/Eigen/src/SparseCore/SparseCwiseBinaryOp.h @@ -73,7 +73,7 @@ class CwiseBinaryOpImpl<BinaryOp,Lhs,Rhs,Sparse>::InnerIterator typedef internal::sparse_cwise_binary_op_inner_iterator_selector< BinaryOp,Lhs,Rhs, InnerIterator> Base; - EIGEN_STRONG_INLINE InnerIterator(const CwiseBinaryOpImpl& binOp, typename CwiseBinaryOpImpl::Index outer) + EIGEN_STRONG_INLINE InnerIterator(const CwiseBinaryOpImpl& binOp, Index outer) : Base(binOp.derived(),outer) {} }; diff --git a/Eigen/src/SparseCore/SparseDenseProduct.h b/Eigen/src/SparseCore/SparseDenseProduct.h index 8c608a622..30975c29c 100644 --- a/Eigen/src/SparseCore/SparseDenseProduct.h +++ b/Eigen/src/SparseCore/SparseDenseProduct.h @@ -111,6 +111,7 @@ template<typename Lhs, typename Rhs, bool Transpose> class SparseDenseOuterProduct<Lhs,Rhs,Transpose>::InnerIterator : public _LhsNested::InnerIterator { typedef typename _LhsNested::InnerIterator Base; + typedef typename SparseDenseOuterProduct::Index Index; public: EIGEN_STRONG_INLINE InnerIterator(const SparseDenseOuterProduct& prod, Index outer) : Base(prod.lhs(), 0), m_outer(outer), m_factor(prod.rhs().coeff(outer)) diff --git a/Eigen/src/SparseCore/SparseDiagonalProduct.h b/Eigen/src/SparseCore/SparseDiagonalProduct.h index 3e314bcfc..1bb590e64 100644 --- a/Eigen/src/SparseCore/SparseDiagonalProduct.h +++ b/Eigen/src/SparseCore/SparseDiagonalProduct.h @@ -78,7 +78,11 @@ class SparseDiagonalProduct EIGEN_SPARSE_PUBLIC_INTERFACE(SparseDiagonalProduct) typedef internal::sparse_diagonal_product_inner_iterator_selector - <_LhsNested,_RhsNested,SparseDiagonalProduct,LhsMode,RhsMode> InnerIterator; + <_LhsNested,_RhsNested,SparseDiagonalProduct,LhsMode,RhsMode> InnerIterator; + + // We do not want ReverseInnerIterator for diagonal-sparse products, + // but this dummy declaration is needed to make diag * sparse * diag compile. + class ReverseInnerIterator; EIGEN_STRONG_INLINE SparseDiagonalProduct(const Lhs& lhs, const Rhs& rhs) : m_lhs(lhs), m_rhs(rhs) diff --git a/Eigen/src/SparseCore/SparseMatrix.h b/Eigen/src/SparseCore/SparseMatrix.h index 2386dfecc..adceafe18 100644 --- a/Eigen/src/SparseCore/SparseMatrix.h +++ b/Eigen/src/SparseCore/SparseMatrix.h @@ -531,59 +531,63 @@ class SparseMatrix */ void conservativeResize(Index rows, Index cols) { - // No change - if (this->rows() == rows && this->cols() == cols) return; + // No change + if (this->rows() == rows && this->cols() == cols) return; + + // If one dimension is null, then there is nothing to be preserved + if(rows==0 || cols==0) return resize(rows,cols); - Index innerChange = IsRowMajor ? cols - this->cols() : rows - this->rows(); - Index outerChange = IsRowMajor ? rows - this->rows() : cols - this->cols(); - Index newInnerSize = IsRowMajor ? cols : rows; + Index innerChange = IsRowMajor ? cols - this->cols() : rows - this->rows(); + Index outerChange = IsRowMajor ? rows - this->rows() : cols - this->cols(); + Index newInnerSize = IsRowMajor ? cols : rows; - // Deals with inner non zeros - if (m_innerNonZeros) - { - // Resize m_innerNonZeros - Index *newInnerNonZeros = static_cast<Index*>(std::realloc(m_innerNonZeros, (m_outerSize + outerChange) * sizeof(Index))); - if (!newInnerNonZeros) internal::throw_std_bad_alloc(); - m_innerNonZeros = newInnerNonZeros; - - for(Index i=m_outerSize; i<m_outerSize+outerChange; i++) - m_innerNonZeros[i] = 0; - } - else if (innerChange < 0) - { - // Inner size decreased: allocate a new m_innerNonZeros - m_innerNonZeros = static_cast<Index*>(std::malloc((m_outerSize+outerChange+1) * sizeof(Index))); - if (!m_innerNonZeros) internal::throw_std_bad_alloc(); - for(Index i = 0; i < m_outerSize; i++) - m_innerNonZeros[i] = m_outerIndex[i+1] - m_outerIndex[i]; - } + // Deals with inner non zeros + if (m_innerNonZeros) + { + // Resize m_innerNonZeros + Index *newInnerNonZeros = static_cast<Index*>(std::realloc(m_innerNonZeros, (m_outerSize + outerChange) * sizeof(Index))); + if (!newInnerNonZeros) internal::throw_std_bad_alloc(); + m_innerNonZeros = newInnerNonZeros; - // Change the m_innerNonZeros in case of a decrease of inner size - if (m_innerNonZeros && innerChange < 0) { - for(Index i = 0; i < m_outerSize + (std::min)(outerChange, Index(0)); i++) - { - Index &n = m_innerNonZeros[i]; - Index start = m_outerIndex[i]; - while (n > 0 && m_data.index(start+n-1) >= newInnerSize) --n; - } + for(Index i=m_outerSize; i<m_outerSize+outerChange; i++) + m_innerNonZeros[i] = 0; + } + else if (innerChange < 0) + { + // Inner size decreased: allocate a new m_innerNonZeros + m_innerNonZeros = static_cast<Index*>(std::malloc((m_outerSize+outerChange+1) * sizeof(Index))); + if (!m_innerNonZeros) internal::throw_std_bad_alloc(); + for(Index i = 0; i < m_outerSize; i++) + m_innerNonZeros[i] = m_outerIndex[i+1] - m_outerIndex[i]; + } + + // Change the m_innerNonZeros in case of a decrease of inner size + if (m_innerNonZeros && innerChange < 0) + { + for(Index i = 0; i < m_outerSize + (std::min)(outerChange, Index(0)); i++) + { + Index &n = m_innerNonZeros[i]; + Index start = m_outerIndex[i]; + while (n > 0 && m_data.index(start+n-1) >= newInnerSize) --n; } - - m_innerSize = newInnerSize; + } + + m_innerSize = newInnerSize; - // Re-allocate outer index structure if necessary - if (outerChange == 0) - return; - - Index *newOuterIndex = static_cast<Index*>(std::realloc(m_outerIndex, (m_outerSize + outerChange + 1) * sizeof(Index))); - if (!newOuterIndex) internal::throw_std_bad_alloc(); - m_outerIndex = newOuterIndex; - if (outerChange > 0) { - Index last = m_outerSize == 0 ? 0 : m_outerIndex[m_outerSize]; - for(Index i=m_outerSize; i<m_outerSize+outerChange+1; i++) - m_outerIndex[i] = last; - } - m_outerSize += outerChange; - + // Re-allocate outer index structure if necessary + if (outerChange == 0) + return; + + Index *newOuterIndex = static_cast<Index*>(std::realloc(m_outerIndex, (m_outerSize + outerChange + 1) * sizeof(Index))); + if (!newOuterIndex) internal::throw_std_bad_alloc(); + m_outerIndex = newOuterIndex; + if (outerChange > 0) + { + Index last = m_outerSize == 0 ? 0 : m_outerIndex[m_outerSize]; + for(Index i=m_outerSize; i<m_outerSize+outerChange+1; i++) + m_outerIndex[i] = last; + } + m_outerSize += outerChange; } /** Resizes the matrix to a \a rows x \a cols matrix and initializes it to zero. diff --git a/Eigen/src/SparseCore/SparseMatrixBase.h b/Eigen/src/SparseCore/SparseMatrixBase.h index 90fee01bc..706f699b8 100644 --- a/Eigen/src/SparseCore/SparseMatrixBase.h +++ b/Eigen/src/SparseCore/SparseMatrixBase.h @@ -89,6 +89,9 @@ template<typename Derived> class SparseMatrixBase : public EigenBase<Derived> */ IsRowMajor = Flags&RowMajorBit ? 1 : 0, + + InnerSizeAtCompileTime = int(IsVectorAtCompileTime) ? int(SizeAtCompileTime) + : int(IsRowMajor) ? int(ColsAtCompileTime) : int(RowsAtCompileTime), #ifndef EIGEN_PARSED_BY_DOXYGEN _HasDirectAccess = (int(Flags)&DirectAccessBit) ? 1 : 0 // workaround sunCC @@ -102,7 +105,7 @@ template<typename Derived> class SparseMatrixBase : public EigenBase<Derived> >::type AdjointReturnType; - typedef SparseMatrix<Scalar, Flags&RowMajorBit ? RowMajor : ColMajor> PlainObject; + typedef SparseMatrix<Scalar, Flags&RowMajorBit ? RowMajor : ColMajor, Index> PlainObject; #ifndef EIGEN_PARSED_BY_DOXYGEN diff --git a/Eigen/src/SparseCore/SparseSelfAdjointView.h b/Eigen/src/SparseCore/SparseSelfAdjointView.h index 60fcf3f40..0eda96bc4 100644 --- a/Eigen/src/SparseCore/SparseSelfAdjointView.h +++ b/Eigen/src/SparseCore/SparseSelfAdjointView.h @@ -75,22 +75,22 @@ template<typename MatrixType, unsigned int UpLo> class SparseSelfAdjointView * Indeed, the SparseSelfadjointView operand is first copied into a temporary SparseMatrix before computing the product. */ template<typename OtherDerived> - SparseSparseProduct<SparseMatrix<Scalar, (internal::traits<OtherDerived>::Flags&RowMajorBit) ? RowMajor : ColMajor,Index>, OtherDerived> + SparseSparseProduct<typename OtherDerived::PlainObject, OtherDerived> operator*(const SparseMatrixBase<OtherDerived>& rhs) const { - return SparseSparseProduct<SparseMatrix<Scalar, (internal::traits<OtherDerived>::Flags&RowMajorBit) ? RowMajor : ColMajor, Index>, OtherDerived>(*this, rhs.derived()); + return SparseSparseProduct<typename OtherDerived::PlainObject, OtherDerived>(*this, rhs.derived()); } - + /** \returns an expression of the matrix product between a sparse matrix \a lhs and a sparse self-adjoint matrix \a rhs. * * Note that there is no algorithmic advantage of performing such a product compared to a general sparse-sparse matrix product. * Indeed, the SparseSelfadjointView operand is first copied into a temporary SparseMatrix before computing the product. */ - template<typename OtherDerived> friend - SparseSparseProduct<OtherDerived, SparseMatrix<Scalar, (internal::traits<OtherDerived>::Flags&RowMajorBit) ? RowMajor : ColMajor,Index> > + template<typename OtherDerived> friend + SparseSparseProduct<OtherDerived, typename OtherDerived::PlainObject > operator*(const SparseMatrixBase<OtherDerived>& lhs, const SparseSelfAdjointView& rhs) { - return SparseSparseProduct< OtherDerived, SparseMatrix<Scalar, (internal::traits<OtherDerived>::Flags&RowMajorBit) ? RowMajor : ColMajor, Index> >(lhs.derived(), rhs.derived()); + return SparseSparseProduct<OtherDerived, typename OtherDerived::PlainObject>(lhs.derived(), rhs); } /** Efficient sparse self-adjoint matrix times dense vector/matrix product */ diff --git a/Eigen/src/SparseCore/SparseTranspose.h b/Eigen/src/SparseCore/SparseTranspose.h index c78c20a2f..7c300ee8d 100644 --- a/Eigen/src/SparseCore/SparseTranspose.h +++ b/Eigen/src/SparseCore/SparseTranspose.h @@ -34,26 +34,28 @@ template<typename MatrixType> class TransposeImpl<MatrixType,Sparse>::InnerItera : public _MatrixTypeNested::InnerIterator { typedef typename _MatrixTypeNested::InnerIterator Base; + typedef typename TransposeImpl::Index Index; public: EIGEN_STRONG_INLINE InnerIterator(const TransposeImpl& trans, typename TransposeImpl<MatrixType,Sparse>::Index outer) : Base(trans.derived().nestedExpression(), outer) {} - inline typename TransposeImpl<MatrixType,Sparse>::Index row() const { return Base::col(); } - inline typename TransposeImpl<MatrixType,Sparse>::Index col() const { return Base::row(); } + Index row() const { return Base::col(); } + Index col() const { return Base::row(); } }; template<typename MatrixType> class TransposeImpl<MatrixType,Sparse>::ReverseInnerIterator : public _MatrixTypeNested::ReverseInnerIterator { typedef typename _MatrixTypeNested::ReverseInnerIterator Base; + typedef typename TransposeImpl::Index Index; public: EIGEN_STRONG_INLINE ReverseInnerIterator(const TransposeImpl& xpr, typename TransposeImpl<MatrixType,Sparse>::Index outer) : Base(xpr.derived().nestedExpression(), outer) {} - inline typename TransposeImpl<MatrixType,Sparse>::Index row() const { return Base::col(); } - inline typename TransposeImpl<MatrixType,Sparse>::Index col() const { return Base::row(); } + Index row() const { return Base::col(); } + Index col() const { return Base::row(); } }; } // end namespace Eigen diff --git a/Eigen/src/SparseCore/SparseTriangularView.h b/Eigen/src/SparseCore/SparseTriangularView.h index 88a345b22..333127b78 100644 --- a/Eigen/src/SparseCore/SparseTriangularView.h +++ b/Eigen/src/SparseCore/SparseTriangularView.h @@ -66,6 +66,7 @@ template<typename MatrixType, int Mode> class SparseTriangularView<MatrixType,Mode>::InnerIterator : public MatrixTypeNestedCleaned::InnerIterator { typedef typename MatrixTypeNestedCleaned::InnerIterator Base; + typedef typename SparseTriangularView::Index Index; public: EIGEN_STRONG_INLINE InnerIterator(const SparseTriangularView& view, Index outer) @@ -135,6 +136,7 @@ template<typename MatrixType, int Mode> class SparseTriangularView<MatrixType,Mode>::ReverseInnerIterator : public MatrixTypeNestedCleaned::ReverseInnerIterator { typedef typename MatrixTypeNestedCleaned::ReverseInnerIterator Base; + typedef typename SparseTriangularView::Index Index; public: EIGEN_STRONG_INLINE ReverseInnerIterator(const SparseTriangularView& view, Index outer) diff --git a/Eigen/src/SparseCore/SparseVector.h b/Eigen/src/SparseCore/SparseVector.h index b05d409c3..7e15c814b 100644 --- a/Eigen/src/SparseCore/SparseVector.h +++ b/Eigen/src/SparseCore/SparseVector.h @@ -45,6 +45,20 @@ struct traits<SparseVector<_Scalar, _Options, _Index> > SupportedAccessPatterns = InnerRandomAccessPattern }; }; + +// Sparse-Vector-Assignment kinds: +enum { + SVA_RuntimeSwitch, + SVA_Inner, + SVA_Outer +}; + +template< typename Dest, typename Src, + int AssignmentKind = !bool(Src::IsVectorAtCompileTime) ? SVA_RuntimeSwitch + : Src::InnerSizeAtCompileTime==1 ? SVA_Outer + : SVA_Inner> +struct sparse_vector_assign_selector; + } template<typename _Scalar, int _Options, typename _Index> @@ -241,11 +255,10 @@ class SparseVector template<typename OtherDerived> inline SparseVector& operator=(const SparseMatrixBase<OtherDerived>& other) { - if ( (bool(OtherDerived::IsVectorAtCompileTime) && int(RowsAtCompileTime)!=int(OtherDerived::RowsAtCompileTime)) - || ((!bool(OtherDerived::IsVectorAtCompileTime)) && ( bool(IsColVector) ? other.cols()>1 : other.rows()>1 ))) - return assign(other.transpose()); - else - return assign(other); + SparseVector tmp(other.size()); + internal::sparse_vector_assign_selector<SparseVector,OtherDerived>::run(tmp,other.derived()); + this->swap(tmp); + return *this; } #ifndef EIGEN_PARSED_BY_DOXYGEN @@ -327,9 +340,6 @@ protected: EIGEN_STATIC_ASSERT((_Options&(ColMajor|RowMajor))==Options,INVALID_MATRIX_TEMPLATE_PARAMETERS); } - template<typename OtherDerived> - EIGEN_DONT_INLINE SparseVector& assign(const SparseMatrixBase<OtherDerived>& _other); - Storage m_data; Index m_size; }; @@ -398,33 +408,40 @@ class SparseVector<Scalar,_Options,_Index>::ReverseInnerIterator const Index m_start; }; -template<typename Scalar, int _Options, typename _Index> -template<typename OtherDerived> -EIGEN_DONT_INLINE SparseVector<Scalar,_Options,_Index>& SparseVector<Scalar,_Options,_Index>::assign(const SparseMatrixBase<OtherDerived>& _other) -{ - const OtherDerived& other(_other.derived()); - const bool needToTranspose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit); - if(needToTranspose) - { - Index size = other.size(); - Index nnz = other.nonZeros(); - resize(size); - reserve(nnz); - for(Index i=0; i<size; ++i) +namespace internal { + +template< typename Dest, typename Src> +struct sparse_vector_assign_selector<Dest,Src,SVA_Inner> { + static void run(Dest& dst, const Src& src) { + eigen_internal_assert(src.innerSize()==src.size()); + for(typename Src::InnerIterator it(src, 0); it; ++it) + dst.insert(it.index()) = it.value(); + } +}; + +template< typename Dest, typename Src> +struct sparse_vector_assign_selector<Dest,Src,SVA_Outer> { + static void run(Dest& dst, const Src& src) { + eigen_internal_assert(src.outerSize()==src.size()); + for(typename Dest::Index i=0; i<src.size(); ++i) { - typename OtherDerived::InnerIterator it(other, i); + typename Src::InnerIterator it(src, i); if(it) - insert(i) = it.value(); + dst.insert(i) = it.value(); } - return *this; } - else - { - // there is no special optimization - return Base::operator=(other); +}; + +template< typename Dest, typename Src> +struct sparse_vector_assign_selector<Dest,Src,SVA_RuntimeSwitch> { + static void run(Dest& dst, const Src& src) { + if(src.outerSize()==1) sparse_vector_assign_selector<Dest,Src,SVA_Inner>::run(dst, src); + else sparse_vector_assign_selector<Dest,Src,SVA_Outer>::run(dst, src); } +}; + } - + } // end namespace Eigen #endif // EIGEN_SPARSEVECTOR_H diff --git a/Eigen/src/SparseCore/SparseView.h b/Eigen/src/SparseCore/SparseView.h index 67eb93245..fd8450463 100644 --- a/Eigen/src/SparseCore/SparseView.h +++ b/Eigen/src/SparseCore/SparseView.h @@ -56,6 +56,7 @@ protected: template<typename MatrixType> class SparseView<MatrixType>::InnerIterator : public _MatrixTypeNested::InnerIterator { + typedef typename SparseView::Index Index; public: typedef typename _MatrixTypeNested::InnerIterator IterBase; InnerIterator(const SparseView& view, Index outer) : diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index ee79c7762..dd9eab2c2 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h @@ -14,9 +14,10 @@ namespace Eigen { -template <typename _MatrixType, typename _OrderingType> class SparseLU; +template <typename _MatrixType, typename _OrderingType = COLAMDOrdering<typename _MatrixType::Index> > class SparseLU; template <typename MappedSparseMatrixType> struct SparseLUMatrixLReturnType; template <typename MatrixLType, typename MatrixUType> struct SparseLUMatrixUReturnType; + /** \ingroup SparseLU_Module * \class SparseLU * @@ -62,7 +63,7 @@ template <typename MatrixLType, typename MatrixUType> struct SparseLUMatrixURetu * "unsupported/Eigen/src/IterativeSolvers/Scaling.h" * * \tparam _MatrixType The type of the sparse matrix. It must be a column-major SparseMatrix<> - * \tparam _OrderingType The ordering method to use, either AMD, COLAMD or METIS + * \tparam _OrderingType The ordering method to use, either AMD, COLAMD or METIS. Default is COLMAD * * * \sa \ref TutorialSparseDirectSolvers @@ -105,9 +106,9 @@ class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typ void simplicialfactorize(const MatrixType& matrix); /** - * Compute the symbolic and numeric factorization of the input sparse matrix. - * The input matrix should be in column-major storage. - */ + * Compute the symbolic and numeric factorization of the input sparse matrix. + * The input matrix should be in column-major storage. + */ void compute (const MatrixType& matrix) { // Analyze @@ -125,38 +126,38 @@ class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typ } /** \returns an expression of the matrix L, internally stored as supernodes - * The only operation available with this expression is the triangular solve - * \code - * y = b; matrixL().solveInPlace(y); - * \endcode - */ + * The only operation available with this expression is the triangular solve + * \code + * y = b; matrixL().solveInPlace(y); + * \endcode + */ SparseLUMatrixLReturnType<SCMatrix> matrixL() const { return SparseLUMatrixLReturnType<SCMatrix>(m_Lstore); } /** \returns an expression of the matrix U, - * The only operation available with this expression is the triangular solve - * \code - * y = b; matrixU().solveInPlace(y); - * \endcode - */ + * The only operation available with this expression is the triangular solve + * \code + * y = b; matrixU().solveInPlace(y); + * \endcode + */ SparseLUMatrixUReturnType<SCMatrix,MappedSparseMatrix<Scalar,ColMajor,Index> > matrixU() const { return SparseLUMatrixUReturnType<SCMatrix, MappedSparseMatrix<Scalar,ColMajor,Index> >(m_Lstore, m_Ustore); } /** - * \returns a reference to the row matrix permutation \f$ P_r \f$ such that \f$P_r A P_c^T = L U\f$ - * \sa colsPermutation() - */ + * \returns a reference to the row matrix permutation \f$ P_r \f$ such that \f$P_r A P_c^T = L U\f$ + * \sa colsPermutation() + */ inline const PermutationType& rowsPermutation() const { return m_perm_r; } /** - * \returns a reference to the column matrix permutation\f$ P_c^T \f$ such that \f$P_r A P_c^T = L U\f$ - * \sa rowsPermutation() - */ + * \returns a reference to the column matrix permutation\f$ P_c^T \f$ such that \f$P_r A P_c^T = L U\f$ + * \sa rowsPermutation() + */ inline const PermutationType& colsPermutation() const { return m_perm_c; @@ -182,7 +183,7 @@ class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typ return internal::solve_retval<SparseLU, Rhs>(*this, B.derived()); } - /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. + /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. * * \sa compute() */ @@ -195,7 +196,7 @@ class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typ return internal::sparse_solve_retval<SparseLU, Rhs>(*this, B.derived()); } - /** \brief Reports whether previous computation was successful. + /** \brief Reports whether previous computation was successful. * * \returns \c Success if computation was succesful, * \c NumericalIssue if the LU factorization reports a problem, zero diagonal for instance @@ -208,9 +209,10 @@ class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typ eigen_assert(m_isInitialized && "Decomposition is not initialized."); return m_info; } + /** - * \returns A string describing the type of error - */ + * \returns A string describing the type of error + */ std::string lastErrorMessage() const { return m_lastError; @@ -240,6 +242,7 @@ class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typ return true; } + /** * \returns the absolute value of the determinant of the matrix of which * *this is the QR decomposition. @@ -249,7 +252,7 @@ class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typ * One way to work around that is to use logAbsDeterminant() instead. * * \sa logAbsDeterminant(), signDeterminant() - */ + */ Scalar absDeterminant() { eigen_assert(m_factorizationIsOk && "The matrix should be factorized first."); @@ -276,8 +279,8 @@ class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typ * of which **this is the QR decomposition * * \note This method is useful to work around the risk of overflow/underflow that's - * inherent to the determinant computation - *a + * inherent to the determinant computation. + * * \sa absDeterminant(), signDeterminant() */ Scalar logAbsDeterminant() const @@ -353,15 +356,15 @@ class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typ // Functions needed by the anaysis phase /** - * Compute the column permutation to minimize the fill-in - * - * - Apply this permutation to the input matrix - - * - * - Compute the column elimination tree on the permuted matrix - * - * - Postorder the elimination tree and the column permutation - * - */ + * Compute the column permutation to minimize the fill-in + * + * - Apply this permutation to the input matrix - + * + * - Compute the column elimination tree on the permuted matrix + * + * - Postorder the elimination tree and the column permutation + * + */ template <typename MatrixType, typename OrderingType> void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat) { @@ -377,11 +380,20 @@ void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat) if (m_perm_c.size()) { m_mat.uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers. FIXME : This vector is filled but not subsequently used. //Then, permute only the column pointers + const Index * outerIndexPtr; + if (mat.isCompressed()) outerIndexPtr = mat.outerIndexPtr(); + else + { + Index *outerIndexPtr_t = new Index[mat.cols()+1]; + for(Index i = 0; i <= mat.cols(); i++) outerIndexPtr_t[i] = m_mat.outerIndexPtr()[i]; + outerIndexPtr = outerIndexPtr_t; + } for (Index i = 0; i < mat.cols(); i++) { - m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = mat.outerIndexPtr()[i]; - m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = mat.outerIndexPtr()[i+1] - mat.outerIndexPtr()[i]; + m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i]; + m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i]; } + if(!mat.isCompressed()) delete[] outerIndexPtr; } // Compute the column elimination tree of the permuted matrix IndexVector firstRowElt; @@ -419,23 +431,23 @@ void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat) /** - * - Numerical factorization - * - Interleaved with the symbolic factorization - * On exit, info is - * - * = 0: successful factorization - * - * > 0: if info = i, and i is - * - * <= A->ncol: U(i,i) is exactly zero. The factorization has - * been completed, but the factor U is exactly singular, - * and division by zero will occur if it is used to solve a - * system of equations. - * - * > A->ncol: number of bytes allocated when memory allocation - * failure occurred, plus A->ncol. If lwork = -1, it is - * the estimated amount of space needed, plus A->ncol. - */ + * - Numerical factorization + * - Interleaved with the symbolic factorization + * On exit, info is + * + * = 0: successful factorization + * + * > 0: if info = i, and i is + * + * <= A->ncol: U(i,i) is exactly zero. The factorization has + * been completed, but the factor U is exactly singular, + * and division by zero will occur if it is used to solve a + * system of equations. + * + * > A->ncol: number of bytes allocated when memory allocation + * failure occurred, plus A->ncol. If lwork = -1, it is + * the estimated amount of space needed, plus A->ncol. + */ template <typename MatrixType, typename OrderingType> void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix) { @@ -453,11 +465,20 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix) { m_mat.uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers. //Then, permute only the column pointers + const Index * outerIndexPtr; + if (matrix.isCompressed()) outerIndexPtr = matrix.outerIndexPtr(); + else + { + Index* outerIndexPtr_t = new Index[matrix.cols()+1]; + for(Index i = 0; i <= matrix.cols(); i++) outerIndexPtr_t[i] = m_mat.outerIndexPtr()[i]; + outerIndexPtr = outerIndexPtr_t; + } for (Index i = 0; i < matrix.cols(); i++) { - m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = matrix.outerIndexPtr()[i]; - m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = matrix.outerIndexPtr()[i+1] - matrix.outerIndexPtr()[i]; + m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i]; + m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i]; } + if(!matrix.isCompressed()) delete[] outerIndexPtr; } else { //FIXME This should not be needed if the empty permutation is handled transparently @@ -511,7 +532,7 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix) m_perm_r.resize(m); m_perm_r.indices().setConstant(-1); marker.setConstant(-1); - m_detPermR = 1.0; // Record the determinant of the row permutation + m_detPermR = 1; // Record the determinant of the row permutation m_glu.supno(0) = emptyIdxLU; m_glu.xsup.setConstant(0); m_glu.xsup(0) = m_glu.xlsub(0) = m_glu.xusub(0) = m_glu.xlusup(0) = Index(0); @@ -630,7 +651,7 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix) } template<typename MappedSupernodalType> -struct SparseLUMatrixLReturnType +struct SparseLUMatrixLReturnType : internal::no_assignment_operator { typedef typename MappedSupernodalType::Index Index; typedef typename MappedSupernodalType::Scalar Scalar; @@ -647,7 +668,7 @@ struct SparseLUMatrixLReturnType }; template<typename MatrixLType, typename MatrixUType> -struct SparseLUMatrixUReturnType +struct SparseLUMatrixUReturnType : internal::no_assignment_operator { typedef typename MatrixLType::Index Index; typedef typename MatrixLType::Scalar Scalar; @@ -700,6 +721,7 @@ struct SparseLUMatrixUReturnType const MatrixLType& m_mapL; const MatrixUType& m_mapU; }; + namespace internal { template<typename _MatrixType, typename Derived, typename Rhs> diff --git a/Eigen/src/SparseLU/SparseLU_Memory.h b/Eigen/src/SparseLU/SparseLU_Memory.h index 6d9570d19..a5158025c 100644 --- a/Eigen/src/SparseLU/SparseLU_Memory.h +++ b/Eigen/src/SparseLU/SparseLU_Memory.h @@ -70,7 +70,7 @@ Index SparseLUImpl<Scalar,Index>::expand(VectorType& vec, Index& length, Index if(num_expansions == 0 || keep_prev) new_len = length ; // First time allocate requested else - new_len = alpha * length ; + new_len = Index(alpha * length); VectorType old_vec; // Temporary vector to hold the previous values if (nbElts > 0 ) @@ -100,7 +100,7 @@ Index SparseLUImpl<Scalar,Index>::expand(VectorType& vec, Index& length, Index do { alpha = (alpha + 1)/2; - new_len = alpha * length ; + new_len = Index(alpha * length); try { vec.resize(new_len); @@ -141,7 +141,7 @@ Index SparseLUImpl<Scalar,Index>::memInit(Index m, Index n, Index annz, Index lw Index& num_expansions = glu.num_expansions; //No memory expansions so far num_expansions = 0; glu.nzumax = glu.nzlumax = (std::max)(fillratio * annz, m*n); // estimated number of nonzeros in U - glu.nzlmax = (std::max)(1., fillratio/4.) * annz; // estimated nnz in L factor + glu.nzlmax = (std::max)(Index(4), fillratio) * annz / 4; // estimated nnz in L factor // Return the estimated size to the user if necessary Index tempSpace; diff --git a/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h b/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h index 3836d1096..ad6f2183f 100644 --- a/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h +++ b/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h @@ -216,13 +216,13 @@ class MappedSuperNodalMatrix<Scalar,Index>::InnerIterator protected: const MappedSuperNodalMatrix& m_matrix; // Supernodal lower triangular matrix - const Index m_outer; // Current column - const Index m_supno; // Current SuperNode number - Index m_idval; //Index to browse the values in the current column - const Index m_startidval; // Start of the column value - const Index m_endidval; // End of the column value - Index m_idrow; //Index to browse the row indices - Index m_endidrow; // End index of row indices of the current column + const Index m_outer; // Current column + const Index m_supno; // Current SuperNode number + Index m_idval; // Index to browse the values in the current column + const Index m_startidval; // Start of the column value + const Index m_endidval; // End of the column value + Index m_idrow; // Index to browse the row indices + Index m_endidrow; // End index of row indices of the current column }; /** @@ -235,17 +235,17 @@ void MappedSuperNodalMatrix<Scalar,Index>::solveInPlace( MatrixBase<Dest>&X) con { Index n = X.rows(); Index nrhs = X.cols(); - const Scalar * Lval = valuePtr(); // Nonzero values - Matrix<Scalar,Dynamic,Dynamic> work(n, nrhs); // working vector + const Scalar * Lval = valuePtr(); // Nonzero values + Matrix<Scalar,Dynamic,Dynamic> work(n, nrhs); // working vector work.setZero(); for (Index k = 0; k <= nsuper(); k ++) { - Index fsupc = supToCol()[k]; // First column of the current supernode - Index istart = rowIndexPtr()[fsupc]; // Pointer index to the subscript of the current column + Index fsupc = supToCol()[k]; // First column of the current supernode + Index istart = rowIndexPtr()[fsupc]; // Pointer index to the subscript of the current column Index nsupr = rowIndexPtr()[fsupc+1] - istart; // Number of rows in the current supernode - Index nsupc = supToCol()[k+1] - fsupc; // Number of columns in the current supernode - Index nrow = nsupr - nsupc; // Number of rows in the non-diagonal part of the supernode - Index irow; //Current index row + Index nsupc = supToCol()[k+1] - fsupc; // Number of columns in the current supernode + Index nrow = nsupr - nsupc; // Number of rows in the non-diagonal part of the supernode + Index irow; //Current index row if (nsupc == 1 ) { @@ -294,4 +294,5 @@ void MappedSuperNodalMatrix<Scalar,Index>::solveInPlace( MatrixBase<Dest>&X) con } // end namespace internal } // end namespace Eigen + #endif // EIGEN_SPARSELU_MATRIX_H diff --git a/Eigen/src/SparseLU/SparseLU_column_dfs.h b/Eigen/src/SparseLU/SparseLU_column_dfs.h index bc4cfbf37..4c04b0e44 100644 --- a/Eigen/src/SparseLU/SparseLU_column_dfs.h +++ b/Eigen/src/SparseLU/SparseLU_column_dfs.h @@ -36,7 +36,7 @@ namespace Eigen { namespace internal { template<typename IndexVector, typename ScalarVector> -struct column_dfs_traits +struct column_dfs_traits : no_assignment_operator { typedef typename ScalarVector::Scalar Scalar; typedef typename IndexVector::Scalar Index; diff --git a/Eigen/src/SparseLU/SparseLU_pruneL.h b/Eigen/src/SparseLU/SparseLU_pruneL.h index 5a855f82f..66460d168 100644 --- a/Eigen/src/SparseLU/SparseLU_pruneL.h +++ b/Eigen/src/SparseLU/SparseLU_pruneL.h @@ -56,7 +56,7 @@ void SparseLUImpl<Scalar,Index>::pruneL(const Index jcol, const IndexVector& per Index jsupno = glu.supno(jcol); Index i,irep,irep1; bool movnum, do_prune = false; - Index kmin, kmax, minloc, maxloc,krow; + Index kmin = 0, kmax = 0, minloc, maxloc,krow; for (i = 0; i < nseg; i++) { irep = segrep(i); |