diff options
60 files changed, 1014 insertions, 642 deletions
diff --git a/Eigen/Core b/Eigen/Core index 7e068cbbd..97aa581e2 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -19,6 +19,12 @@ // defined e.g. EIGEN_DONT_ALIGN) so it needs to be done before we do anything with vectorization. #include "src/Core/util/Macros.h" +// Disable the ipa-cp-clone optimization flag with MinGW 6.x or newer (enabled by default with -O3) +// See http://eigen.tuxfamily.org/bz/show_bug.cgi?id=556 for details. +#if defined(__MINGW32__) && EIGEN_GNUC_AT_LEAST(4,6) + #pragma GCC optimize ("-fno-ipa-cp-clone") +#endif + #include <complex> // this include file manages BLAS and MKL related macros 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); diff --git a/cmake/EigenConfigureTesting.cmake b/cmake/EigenConfigureTesting.cmake index 6216a9009..228d29e97 100644 --- a/cmake/EigenConfigureTesting.cmake +++ b/cmake/EigenConfigureTesting.cmake @@ -24,17 +24,19 @@ set(CMAKE_MAKE_PROGRAM "@EIGEN_MAKECOMMAND_PLACEHOLDER@") # This call activates testing and generates the DartConfiguration.tcl include(CTest) +set(EIGEN_TEST_BUILD_FLAGS " " CACHE STRING "Options passed to the build command of unit tests") + # overwrite default DartConfiguration.tcl # The worarounds are different for each version of the MSVC IDE if(MSVC_IDE) if(CMAKE_MAKE_PROGRAM_SAVE MATCHES "devenv") # devenv - set(EIGEN_MAKECOMMAND_PLACEHOLDER "${CMAKE_MAKE_PROGRAM_SAVE} Eigen.sln /build \"Release\" /project buildtests \n# ") + set(EIGEN_MAKECOMMAND_PLACEHOLDER "${CMAKE_MAKE_PROGRAM_SAVE} Eigen.sln /build \"Release\" /project buildtests ${EIGEN_TEST_BUILD_FLAGS} \n# ") else() # msbuild - set(EIGEN_MAKECOMMAND_PLACEHOLDER "${CMAKE_MAKE_PROGRAM_SAVE} buildtests.vcxproj /p:Configuration=\${CTEST_CONFIGURATION_TYPE} \n# ") + set(EIGEN_MAKECOMMAND_PLACEHOLDER "${CMAKE_MAKE_PROGRAM_SAVE} buildtests.vcxproj /p:Configuration=\${CTEST_CONFIGURATION_TYPE} ${EIGEN_TEST_BUILD_FLAGS}\n# ") endif() else() # for make and nmake - set(EIGEN_MAKECOMMAND_PLACEHOLDER "${CMAKE_MAKE_PROGRAM_SAVE} buildtests") + set(EIGEN_MAKECOMMAND_PLACEHOLDER "${CMAKE_MAKE_PROGRAM_SAVE} buildtests ${EIGEN_TEST_BUILD_FLAGS}") endif() # copy ctest properties, which currently diff --git a/doc/FunctionsTakingEigenTypes.dox b/doc/FunctionsTakingEigenTypes.dox index ad315edf0..152dda47d 100644 --- a/doc/FunctionsTakingEigenTypes.dox +++ b/doc/FunctionsTakingEigenTypes.dox @@ -1,18 +1,18 @@ namespace Eigen { -/** \page TopicFunctionTakingEigenTypes Writing Functions Taking Eigen Types as Parameters +/** \page TopicFunctionTakingEigenTypes Writing Functions Taking %Eigen Types as Parameters -Eigen's use of expression templates results in potentially every expression being of a different type. If you pass such an expression to a function taking a parameter of type Matrix, your expression will implicitly be evaluated into a temporary Matrix, which will then be passed to the function. This means that you lose the benefit of expression templates. Concretely, this has two drawbacks: +%Eigen's use of expression templates results in potentially every expression being of a different type. If you pass such an expression to a function taking a parameter of type Matrix, your expression will implicitly be evaluated into a temporary Matrix, which will then be passed to the function. This means that you lose the benefit of expression templates. Concretely, this has two drawbacks: \li The evaluation into a temporary may be useless and inefficient; \li This only allows the function to read from the expression, not to write to it. -Fortunately, all this myriad of expression types have in common that they all inherit a few common, templated base classes. By letting your function take templated parameters of these base types, you can let them play nicely with Eigen's expression templates. +Fortunately, all this myriad of expression types have in common that they all inherit a few common, templated base classes. By letting your function take templated parameters of these base types, you can let them play nicely with %Eigen's expression templates. \eigenAutoToc \section TopicFirstExamples Some First Examples -This section will provide simple examples for different types of objects Eigen is offering. Before starting with the actual examples, we need to recapitulate which base objects we can work with (see also \ref TopicClassHierarchy). +This section will provide simple examples for different types of objects %Eigen is offering. Before starting with the actual examples, we need to recapitulate which base objects we can work with (see also \ref TopicClassHierarchy). \li MatrixBase: The common base class for all dense matrix expressions (as opposed to array expressions, as opposed to sparse and special matrix classes). Use it in functions that are meant to work only on dense matrices. \li ArrayBase: The common base class for all dense array expressions (as opposed to matrix expressions, etc). Use it in functions that are meant to work only on arrays. @@ -20,7 +20,7 @@ This section will provide simple examples for different types of objects Eigen i \li EigenBase: The base class unifying all types of objects that can be evaluated into dense matrices or arrays, for example special matrix classes such as diagonal matrices, permutation matrices, etc. It can be used in functions that are meant to work on any such general type. <b> %EigenBase Example </b><br/><br/> -Prints the dimensions of the most generic object present in Eigen. It coulde be any matrix expressions, any dense or sparse matrix and any array. +Prints the dimensions of the most generic object present in %Eigen. It could be any matrix expressions, any dense or sparse matrix and any array. <table class="example"> <tr><th>Example:</th><th>Output:</th></tr> <tr><td> @@ -76,9 +76,42 @@ where the first argument \c v1 is a vector and the second argument \c 2*v2 is an These examples are just intended to give the reader a first impression of how functions can be written which take a plain and constant Matrix or Array argument. They are also intended to give the reader an idea about the most common base classes being the optimal candidates for functions. In the next section we will look in more detail at an example and the different ways it can be implemented, while discussing each implementation's problems and advantages. For the discussion below, Matrix and Array as well as MatrixBase and ArrayBase can be exchanged and all arguments still hold. + +\section TopicUsingRefClass How to write generic, but non-templated function? + +In all the previous examples, the functions had to be template functions. This approach allows to write very generic code, but it is often desirable to write non templated function and still keep some level of genericity to avoid stupid copies of the arguments. The typical example is to write functions accepting both a MatrixXf or a block of a MatrixXf. This exactly the purpose of the Ref class. Here is a simple example: + +<table class="example"> +<tr><th>Example:</th><th>Output:</th></tr> +<tr><td> +\include function_taking_ref.cpp +</td> +<td> +\verbinclude function_taking_ref.out +</td></tr></table> +In the first two calls to inv_cond, no copy occur because the memory layout of the arguments matches the memory layout accepted by Ref<MatrixXf>. However, in the last call, we have a generic expression that will be automatically evaluated into a temporary MatrixXf by the Ref<> object. + +A Ref object can also be writable. Here is an example of a function computing the covariance matrix of two input matrices where each row is an observation: +\code +void cov(const Ref<const MatrixXf> x, const Ref<const MatrixXf> y, Ref<MatrixXf> C) +{ + const float num_observations = static_cast<float>(x.rows()); + const RowVectorXf x_mean = x.colwise().sum() / num_observations; + const RowVectorXf y_mean = y.colwise().sum() / num_observations; + C = (x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations; +} +\endcode +and here are two examples calling cov without any copy: +\code +MatrixXf m1, m2, m3 +cov(m1, m2, m3); +cov(m1.leftCols<3>(), m2.leftCols<3>(), m3.topLeftCorner<3,3>()); +\endcode +The Ref<> class has two other optional template arguments allowing to control the kind of memory layout that can be accepted without any copy. See the class Ref documentation for the details. + \section TopicPlainFunctionsWorking In which cases do functions taking plain Matrix or Array arguments work? -Let's assume one wants to write a function computing the covariance matrix of two input matrices where each row is an observation. The implementation of this function might look like this +Without using template functions, and without the Ref class, a naive implementation of the previous cov function might look like this \code MatrixXf cov(const MatrixXf& x, const MatrixXf& y) { @@ -88,7 +121,7 @@ MatrixXf cov(const MatrixXf& x, const MatrixXf& y) return (x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations; } \endcode -and contrary to what one might think at first, this implementation is fine unless you require a genric implementation that works with double matrices too and unless you do not care about temporary objects. Why is that the case? Where are temporaries involved? How can code as given below compile? +and contrary to what one might think at first, this implementation is fine unless you require a generic implementation that works with double matrices too and unless you do not care about temporary objects. Why is that the case? Where are temporaries involved? How can code as given below compile? \code MatrixXf x,y,z; MatrixXf C = cov(x,y+z); @@ -97,6 +130,7 @@ In this special case, the example is fine and will be working because both param \b Note: Functions taking \e const references to Matrix (or Array) can process expressions at the cost of temporaries. + \section TopicPlainFunctionsFailing In which cases do functions taking a plain Matrix or Array argument fail? Here, we consider a slightly modified version of the function given above. This time, we do not want to return the result but pass an additional non-const paramter which allows us to store the result. A first naive implementation might look as follows. @@ -149,7 +183,7 @@ MatrixXf y = MatrixXf::Random(100,3); MatrixXf C; cov(x, y, C); \endcode -This is not the case anymore, when we are using an implementation taking MatrixBase as a parameter. In general, Eigen supports automatic resizing but it is not possible to do so on expressions. Why should resizing of a matrix Block be allowed? It is a reference to a sub-matrix and we definitely don't want to resize that. So how can we incorporate resizing if we cannot resize on MatrixBase? The solution is to resize the derived object as in this implementation. +This is not the case anymore, when we are using an implementation taking MatrixBase as a parameter. In general, %Eigen supports automatic resizing but it is not possible to do so on expressions. Why should resizing of a matrix Block be allowed? It is a reference to a sub-matrix and we definitely don't want to resize that. So how can we incorporate resizing if we cannot resize on MatrixBase? The solution is to resize the derived object as in this implementation. \code template <typename Derived, typename OtherDerived> void cov(const MatrixBase<Derived>& x, const MatrixBase<Derived>& y, MatrixBase<OtherDerived> const & C_) diff --git a/doc/TutorialSparse.dox b/doc/TutorialSparse.dox index 41bae6b5c..dbfb4a9eb 100644 --- a/doc/TutorialSparse.dox +++ b/doc/TutorialSparse.dox @@ -116,7 +116,7 @@ Describing the \a buildProblem and \a save functions is out of the scope of this The SparseMatrix and SparseVector classes take three template arguments: * the scalar type (e.g., double) - * the storage order (ColMajor or RowMajor, the default is RowMajor) + * the storage order (ColMajor or RowMajor, the default is ColMajor) * the inner index type (default is \c int). As for dense Matrix objects, constructors takes the size of the object. diff --git a/doc/UsingIntelMKL.dox b/doc/UsingIntelMKL.dox index 7c5981460..4b624a156 100644 --- a/doc/UsingIntelMKL.dox +++ b/doc/UsingIntelMKL.dox @@ -40,7 +40,7 @@ Since Eigen version 3.1 and later, users can benefit from built-in Intel MKL opt <a href="http://eigen.tuxfamily.org/Counter/redirect_to_mkl.php"> Intel MKL </a> provides highly optimized multi-threaded mathematical routines for x86-compatible architectures. Intel MKL is available on Linux, Mac and Windows for both Intel64 and IA32 architectures. -\warning Be aware that Intel® MKL is a proprietary software. It is the responsibility of the users to buy MKL licenses for their products. Moreover, the license of the user product has to allow linking to proprietary software that excludes any unmodified versions of the GPL. As a consequence, this also means that Eigen has to be used through the LGPL3+ license. +\warning Be aware that Intel® MKL is a proprietary software. It is the responsibility of the users to buy MKL licenses for their products. Moreover, the license of the user product has to allow linking to proprietary software that excludes any unmodified versions of the GPL. Using Intel MKL through Eigen is easy: -# define the \c EIGEN_USE_MKL_ALL macro before including any Eigen's header diff --git a/doc/examples/function_taking_ref.cpp b/doc/examples/function_taking_ref.cpp new file mode 100644 index 000000000..162a202e4 --- /dev/null +++ b/doc/examples/function_taking_ref.cpp @@ -0,0 +1,19 @@ +#include <iostream> +#include <Eigen/SVD> +using namespace Eigen; +using namespace std; + +float inv_cond(const Ref<const MatrixXf>& a) +{ + const VectorXf sing_vals = a.jacobiSvd().singularValues(); + return sing_vals(sing_vals.size()-1) / sing_vals(0); +} + +int main() +{ + Matrix4f m = Matrix4f::Random(); + cout << "matrix m:" << endl << m << endl << endl; + cout << "inv_cond(m): " << inv_cond(m) << endl; + cout << "inv_cond(m(1:3,1:3)): " << inv_cond(m.topLeftCorner(3,3)) << endl; + cout << "inv_cond(m+I): " << inv_cond(m+Matrix4f::Identity()) << endl; +} diff --git a/lapack/CMakeLists.txt b/lapack/CMakeLists.txt index 296e05838..7e7444326 100644 --- a/lapack/CMakeLists.txt +++ b/lapack/CMakeLists.txt @@ -29,345 +29,65 @@ set(EigenLapack_SRCS ${EigenLapack_SRCS} dlapy2.f dlapy3.f slapy2.f slapy3.f clacgv.f zlacgv.f slamch.f dlamch.f + second_NONE.f dsecnd_NONE.f ) -get_filename_component(eigen_full_path_to_reference_to_reference_lapack "./reference/" ABSOLUTE) -if(EXISTS ${eigen_full_path_to_reference_to_reference_lapack}) -set(EigenLapack_SRCS ${EigenLapack_SRCS} -# reference/dpotrf.f reference/zpotrf.f reference/cpotrf.f reference/spotrf.f -# reference/dpotrs.f reference/spotrs.f reference/zpotrs.f reference/cpotrs.f -# reference/dgetrf.f reference/cgetrf.f reference/sgetrf.f reference/zgetrf.f -# reference/cgetrs.f reference/dgetrs.f reference/sgetrs.f reference/zgetrs.f -# reference/dsyev.f reference/ssyev.f -reference/dlamch.f reference/ilaver.f reference/lsame.f reference/slamch.f reference/second_NONE.f reference/dsecnd_NONE.f -reference/cbdsqr.f reference/ctbrfs.f reference/dorml2.f reference/sla_porfsx_extended.f reference/zggglm.f -reference/cgbbrd.f reference/ctbtrs.f reference/dormlq.f reference/sla_porpvgrw.f reference/zgghrd.f -reference/cgbcon.f reference/ctfsm.f reference/dormql.f reference/slapy2.f reference/zgglse.f -reference/cgbequb.f reference/ctftri.f reference/dormqr.f reference/slapy3.f reference/zggqrf.f -reference/cgbequ.f reference/ctfttp.f reference/dormr2.f reference/slaqgb.f reference/zggrqf.f -reference/cgbrfs.f reference/ctfttr.f reference/dormr3.f reference/slaqge.f reference/zggsvd.f -reference/cgbrfsx.f reference/ctgevc.f reference/dormrq.f reference/slaqp2.f reference/zggsvp.f -reference/cgbsv.f reference/ctgex2.f reference/dormrz.f reference/slaqps.f reference/zgtcon.f -reference/cgbsvx.f reference/ctgexc.f reference/dormtr.f reference/slaqr0.f reference/zgtrfs.f -reference/cgbsvxx.f reference/ctgsen.f reference/dpbcon.f reference/slaqr1.f reference/zgtsv.f -reference/cgbtf2.f reference/ctgsja.f reference/dpbequ.f reference/slaqr2.f reference/zgtsvx.f -reference/cgbtrf.f reference/ctgsna.f reference/dpbrfs.f reference/slaqr3.f reference/zgttrf.f -reference/cgbtrs.f reference/ctgsy2.f reference/dpbstf.f reference/slaqr4.f reference/zgttrs.f -reference/cgebak.f reference/ctgsyl.f reference/dpbsv.f reference/slaqr5.f reference/zgtts2.f -reference/cgebal.f reference/ctpcon.f reference/dpbsvx.f reference/slaqsb.f reference/zhbevd.f -reference/cgebd2.f reference/ctprfs.f reference/dpbtf2.f reference/slaqsp.f reference/zhbev.f -reference/cgebrd.f reference/ctptri.f reference/dpbtrf.f reference/slaqsy.f reference/zhbevx.f -reference/cgecon.f reference/ctptrs.f reference/dpbtrs.f reference/slaqtr.f reference/zhbgst.f -reference/cgeequb.f reference/ctpttf.f reference/dpftrf.f reference/slar1v.f reference/zhbgvd.f -reference/cgeequ.f reference/ctpttr.f reference/dpftri.f reference/slar2v.f reference/zhbgv.f -reference/cgees.f reference/ctrcon.f reference/dpftrs.f reference/slarfb.f reference/zhbgvx.f -reference/cgeesx.f reference/ctrevc.f reference/dpocon.f reference/slarf.f reference/zhbtrd.f -reference/cgeev.f reference/ctrexc.f reference/dpoequb.f reference/slarfg.f reference/zhecon.f -reference/cgeevx.f reference/ctrrfs.f reference/dpoequ.f reference/slarfp.f reference/zheequb.f -reference/cgegs.f reference/ctrsen.f reference/dporfs.f reference/slarft.f reference/zheevd.f -reference/cgegv.f reference/ctrsna.f reference/dporfsx.f reference/slarfx.f reference/zheev.f -reference/cgehd2.f reference/ctrsyl.f reference/dposv.f reference/slargv.f reference/zheevr.f -reference/cgehrd.f reference/ctrti2.f reference/dposvx.f reference/slarnv.f reference/zheevx.f -reference/cgelq2.f reference/ctrtri.f reference/dposvxx.f reference/sla_rpvgrw.f reference/zhegs2.f -reference/cgelqf.f reference/ctrtrs.f reference/dpotf2.f reference/slarra.f reference/zhegst.f -reference/cgelsd.f reference/ctrttf.f -reference/slarrb.f reference/zhegvd.f -reference/cgels.f reference/ctrttp.f reference/dpotri.f reference/slarrc.f reference/zhegv.f -reference/cgelss.f reference/ctzrqf.f reference/slarrd.f reference/zhegvx.f -reference/cgelsx.f reference/ctzrzf.f reference/dppcon.f reference/slarre.f reference/zherfs.f -reference/cgelsy.f reference/cung2l.f reference/dppequ.f reference/slarrf.f reference/zherfsx.f -reference/cgeql2.f reference/cung2r.f reference/dpprfs.f reference/slarrj.f reference/zhesv.f -reference/cgeqlf.f reference/cungbr.f reference/dppsv.f reference/slarrk.f reference/zhesvx.f -reference/cgeqp3.f reference/cunghr.f reference/dppsvx.f reference/slarrr.f reference/zhesvxx.f -reference/cgeqpf.f reference/cungl2.f reference/dpptrf.f reference/slarrv.f reference/zhetd2.f -reference/cgeqr2.f reference/cunglq.f reference/dpptri.f reference/slarscl2.f reference/zhetf2.f -reference/cgeqrf.f reference/cungql.f reference/dpptrs.f reference/slartg.f reference/zhetrd.f -reference/cgerfs.f reference/cungqr.f reference/dpstf2.f reference/slartv.f reference/zhetrf.f -reference/cgerfsx.f reference/cungr2.f reference/dpstrf.f reference/slaruv.f reference/zhetri.f -reference/cgerq2.f reference/cungrq.f reference/dptcon.f reference/slarzb.f reference/zhetrs.f -reference/cgerqf.f reference/cungtr.f reference/dpteqr.f reference/slarz.f reference/zhfrk.f -reference/cgesc2.f reference/cunm2l.f reference/dptrfs.f reference/slarzt.f reference/zhgeqz.f -reference/cgesdd.f reference/cunm2r.f reference/dptsv.f reference/slas2.f reference/zhpcon.f -reference/cgesvd.f reference/cunmbr.f reference/dptsvx.f reference/slascl2.f reference/zhpevd.f -reference/cgesv.f reference/cunmhr.f reference/dpttrf.f reference/slascl.f reference/zhpev.f -reference/cgesvx.f reference/cunml2.f reference/dpttrs.f reference/slasd0.f reference/zhpevx.f -reference/cgesvxx.f reference/cunmlq.f reference/dptts2.f reference/slasd1.f reference/zhpgst.f -reference/cgetc2.f reference/cunmql.f reference/drscl.f reference/slasd2.f reference/zhpgvd.f -reference/cgetf2.f reference/cunmqr.f reference/dsbevd.f reference/slasd3.f reference/zhpgv.f -reference/cunmr2.f reference/dsbev.f reference/slasd4.f reference/zhpgvx.f -reference/cgetri.f reference/cunmr3.f reference/dsbevx.f reference/slasd5.f reference/zhprfs.f - reference/cunmrq.f reference/dsbgst.f reference/slasd6.f reference/zhpsv.f -reference/cggbak.f reference/cunmrz.f reference/dsbgvd.f reference/slasd7.f reference/zhpsvx.f -reference/cggbal.f reference/cunmtr.f reference/dsbgv.f reference/slasd8.f reference/zhptrd.f -reference/cgges.f reference/cupgtr.f reference/dsbgvx.f reference/slasda.f reference/zhptrf.f -reference/cggesx.f reference/cupmtr.f reference/dsbtrd.f reference/slasdq.f reference/zhptri.f -reference/cggev.f reference/dbdsdc.f reference/dsfrk.f reference/slasdt.f reference/zhptrs.f -reference/cggevx.f reference/dbdsqr.f reference/dsgesv.f reference/slaset.f reference/zhsein.f -reference/cggglm.f reference/ddisna.f reference/dspcon.f reference/slasq1.f reference/zhseqr.f -reference/cgghrd.f reference/dgbbrd.f reference/dspevd.f reference/slasq2.f reference/zlabrd.f -reference/cgglse.f reference/dgbcon.f reference/dspev.f reference/slasq3.f reference/zlacgv.f -reference/cggqrf.f reference/dgbequb.f reference/dspevx.f reference/slasq4.f reference/zlacn2.f -reference/cggrqf.f reference/dgbequ.f reference/dspgst.f reference/slasq5.f reference/zlacon.f -reference/cggsvd.f reference/dgbrfs.f reference/dspgvd.f reference/slasq6.f reference/zlacp2.f -reference/cggsvp.f reference/dgbrfsx.f reference/dspgv.f reference/slasr.f reference/zlacpy.f -reference/cgtcon.f reference/dgbsv.f reference/dspgvx.f reference/slasrt.f reference/zlacrm.f -reference/cgtrfs.f reference/dgbsvx.f reference/dsposv.f reference/slassq.f reference/zlacrt.f -reference/cgtsv.f reference/dgbsvxx.f reference/dsprfs.f reference/slasv2.f reference/zladiv.f -reference/cgtsvx.f reference/dgbtf2.f reference/dspsv.f reference/slaswp.f reference/zlaed0.f -reference/cgttrf.f reference/dgbtrf.f reference/dspsvx.f reference/slasy2.f reference/zlaed7.f -reference/cgttrs.f reference/dgbtrs.f reference/dsptrd.f reference/sla_syamv.f reference/zlaed8.f -reference/cgtts2.f reference/dgebak.f reference/dsptrf.f reference/slasyf.f reference/zlaein.f -reference/chbevd.f reference/dgebal.f reference/dsptri.f reference/sla_syrcond.f reference/zlaesy.f -reference/chbev.f reference/dgebd2.f reference/dsptrs.f reference/sla_syrfsx_extended.f reference/zlaev2.f -reference/chbevx.f reference/dgebrd.f reference/dstebz.f reference/sla_syrpvgrw.f reference/zlag2c.f -reference/chbgst.f reference/dgecon.f reference/dstedc.f reference/slatbs.f reference/zla_gbamv.f -reference/chbgvd.f reference/dgeequb.f reference/dstegr.f reference/slatdf.f reference/zla_gbrcond_c.f -reference/chbgv.f reference/dgeequ.f reference/dstein.f reference/slatps.f reference/zla_gbrcond_x.f -reference/chbgvx.f reference/dgees.f reference/dstemr.f reference/slatrd.f reference/zla_gbrfsx_extended.f -reference/chbtrd.f reference/dgeesx.f reference/dsteqr.f reference/slatrs.f reference/zla_gbrpvgrw.f -reference/checon.f reference/dgeev.f reference/dsterf.f reference/slatrz.f reference/zla_geamv.f -reference/cheequb.f reference/dgeevx.f reference/dstevd.f reference/slatzm.f reference/zla_gercond_c.f -reference/cheevd.f reference/dgegs.f reference/dstev.f reference/slauu2.f reference/zla_gercond_x.f -reference/cheev.f reference/dgegv.f reference/dstevr.f reference/slauum.f reference/zla_gerfsx_extended.f -reference/cheevr.f reference/dgehd2.f reference/dstevx.f reference/sla_wwaddw.f reference/zlags2.f -reference/cheevx.f reference/dgehrd.f reference/dsycon.f reference/sopgtr.f reference/zlagtm.f -reference/chegs2.f reference/dgejsv.f reference/dsyequb.f reference/sopmtr.f reference/zla_heamv.f -reference/chegst.f reference/dgelq2.f reference/dsyevd.f reference/sorg2l.f reference/zlahef.f -reference/chegvd.f reference/dgelqf.f reference/sorg2r.f reference/zla_hercond_c.f -reference/chegv.f reference/dgelsd.f reference/dsyevr.f reference/sorgbr.f reference/zla_hercond_x.f -reference/chegvx.f reference/dgels.f reference/dsyevx.f reference/sorghr.f reference/zla_herfsx_extended.f -reference/cherfs.f reference/dgelss.f reference/dsygs2.f reference/sorgl2.f reference/zla_herpvgrw.f -reference/cherfsx.f reference/dgelsx.f reference/dsygst.f reference/sorglq.f reference/zlahqr.f -reference/chesv.f reference/dgelsy.f reference/dsygvd.f reference/sorgql.f reference/zlahr2.f -reference/chesvx.f reference/dgeql2.f reference/dsygv.f reference/sorgqr.f reference/zlahrd.f -reference/chesvxx.f reference/dgeqlf.f reference/dsygvx.f reference/sorgr2.f reference/zlaic1.f -reference/chetd2.f reference/dgeqp3.f reference/dsyrfs.f reference/sorgrq.f reference/zla_lin_berr.f -reference/chetf2.f reference/dgeqpf.f reference/dsyrfsx.f reference/sorgtr.f reference/zlals0.f -reference/chetrd.f reference/dgeqr2.f reference/dsysv.f reference/sorm2l.f reference/zlalsa.f -reference/chetrf.f reference/dgeqrf.f reference/dsysvx.f reference/sorm2r.f reference/zlalsd.f -reference/chetri.f reference/dgerfs.f reference/dsysvxx.f reference/sormbr.f reference/zlangb.f -reference/chetrs.f reference/dgerfsx.f reference/dsytd2.f reference/sormhr.f reference/zlange.f -reference/chfrk.f reference/dgerq2.f reference/dsytf2.f reference/sorml2.f reference/zlangt.f -reference/chgeqz.f reference/dgerqf.f reference/dsytrd.f reference/sormlq.f reference/zlanhb.f -reference/chla_transtype.f reference/dgesc2.f reference/dsytrf.f reference/sormql.f reference/zlanhe.f -reference/chpcon.f reference/dgesdd.f reference/dsytri.f reference/sormqr.f reference/zlanhf.f -reference/chpevd.f reference/dgesvd.f reference/dsytrs.f reference/sormr2.f reference/zlanhp.f -reference/chpev.f reference/dgesv.f reference/dtbcon.f reference/sormr3.f reference/zlanhs.f -reference/chpevx.f reference/dgesvj.f reference/dtbrfs.f reference/sormrq.f reference/zlanht.f -reference/chpgst.f reference/dgesvx.f reference/dtbtrs.f reference/sormrz.f reference/zlansb.f -reference/chpgvd.f reference/dgesvxx.f reference/dtfsm.f reference/sormtr.f reference/zlansp.f -reference/chpgv.f reference/dgetc2.f reference/dtftri.f reference/spbcon.f reference/zlansy.f -reference/chpgvx.f reference/dgetf2.f reference/dtfttp.f reference/spbequ.f reference/zlantb.f -reference/chprfs.f -reference/dtfttr.f reference/spbrfs.f reference/zlantp.f -reference/chpsv.f reference/dgetri.f reference/dtgevc.f reference/spbstf.f reference/zlantr.f -reference/chpsvx.f reference/dtgex2.f reference/spbsv.f reference/zlapll.f -reference/chptrd.f reference/dggbak.f reference/dtgexc.f reference/spbsvx.f reference/zlapmt.f -reference/chptrf.f reference/dggbal.f reference/dtgsen.f reference/spbtf2.f reference/zla_porcond_c.f -reference/chptri.f reference/dgges.f reference/dtgsja.f reference/spbtrf.f reference/zla_porcond_x.f -reference/chptrs.f reference/dggesx.f reference/dtgsna.f reference/spbtrs.f reference/zla_porfsx_extended.f -reference/chsein.f reference/dggev.f reference/dtgsy2.f reference/spftrf.f reference/zla_porpvgrw.f -reference/chseqr.f reference/dggevx.f reference/dtgsyl.f reference/spftri.f reference/zlaqgb.f -reference/clabrd.f reference/dggglm.f reference/dtpcon.f reference/spftrs.f reference/zlaqge.f -reference/clacgv.f reference/dgghrd.f reference/dtprfs.f reference/spocon.f reference/zlaqhb.f -reference/clacn2.f reference/dgglse.f reference/dtptri.f reference/spoequb.f reference/zlaqhe.f -reference/clacon.f reference/dggqrf.f reference/dtptrs.f reference/spoequ.f reference/zlaqhp.f -reference/clacp2.f reference/dggrqf.f reference/dtpttf.f reference/sporfs.f reference/zlaqp2.f -reference/clacpy.f reference/dggsvd.f reference/dtpttr.f reference/sporfsx.f reference/zlaqps.f -reference/clacrm.f reference/dggsvp.f reference/dtrcon.f reference/sposv.f reference/zlaqr0.f -reference/clacrt.f reference/dgsvj0.f reference/dtrevc.f reference/sposvx.f reference/zlaqr1.f -reference/cladiv.f reference/dgsvj1.f reference/dtrexc.f reference/sposvxx.f reference/zlaqr2.f -reference/claed0.f reference/dgtcon.f reference/dtrrfs.f reference/spotf2.f reference/zlaqr3.f -reference/claed7.f reference/dgtrfs.f reference/dtrsen.f -reference/zlaqr4.f -reference/claed8.f reference/dgtsv.f reference/dtrsna.f reference/spotri.f reference/zlaqr5.f -reference/claein.f reference/dgtsvx.f reference/dtrsyl.f reference/zlaqsb.f -reference/claesy.f reference/dgttrf.f reference/dtrti2.f reference/sppcon.f reference/zlaqsp.f -reference/claev2.f reference/dgttrs.f reference/dtrtri.f reference/sppequ.f reference/zlaqsy.f -reference/clag2z.f reference/dgtts2.f reference/dtrtrs.f reference/spprfs.f reference/zlar1v.f -reference/cla_gbamv.f reference/dhgeqz.f reference/dtrttf.f reference/sppsv.f reference/zlar2v.f -reference/cla_gbrcond_c.f reference/dhsein.f reference/dtrttp.f reference/sppsvx.f reference/zlarcm.f -reference/cla_gbrcond_x.f reference/dhseqr.f reference/dtzrqf.f reference/spptrf.f reference/zlarfb.f -reference/cla_gbrfsx_extended.f reference/disnan.f reference/dtzrzf.f reference/spptri.f reference/zlarf.f -reference/cla_gbrpvgrw.f reference/dlabad.f reference/dzsum1.f reference/spptrs.f reference/zlarfg.f -reference/cla_geamv.f reference/dlabrd.f reference/icmax1.f reference/spstf2.f reference/zlarfp.f -reference/cla_gercond_c.f reference/dlacn2.f reference/ieeeck.f reference/spstrf.f reference/zlarft.f -reference/cla_gercond_x.f reference/dlacon.f reference/ilaclc.f reference/sptcon.f reference/zlarfx.f -reference/cla_gerfsx_extended.f reference/dlacpy.f reference/ilaclr.f reference/spteqr.f reference/zlargv.f -reference/clags2.f reference/dladiv.f reference/iladiag.f reference/sptrfs.f reference/zlarnv.f -reference/clagtm.f reference/dlae2.f reference/iladlc.f reference/sptsv.f reference/zla_rpvgrw.f -reference/cla_heamv.f reference/dlaebz.f reference/iladlr.f reference/sptsvx.f reference/zlarrv.f -reference/clahef.f reference/dlaed0.f reference/ilaenv.f reference/spttrf.f reference/zlarscl2.f -reference/cla_hercond_c.f reference/dlaed1.f reference/ilaprec.f reference/spttrs.f reference/zlartg.f -reference/cla_hercond_x.f reference/dlaed2.f reference/ilaslc.f reference/sptts2.f reference/zlartv.f -reference/cla_herfsx_extended.f reference/dlaed3.f reference/ilaslr.f reference/srscl.f reference/zlarzb.f -reference/cla_herpvgrw.f reference/dlaed4.f reference/ilatrans.f reference/ssbevd.f reference/zlarz.f -reference/clahqr.f reference/dlaed5.f reference/ilauplo.f reference/ssbev.f reference/zlarzt.f -reference/clahr2.f reference/dlaed6.f reference/ilaver.f reference/ssbevx.f reference/zlascl2.f -reference/clahrd.f reference/dlaed7.f reference/ilazlc.f reference/ssbgst.f reference/zlascl.f -reference/claic1.f reference/dlaed8.f reference/ilazlr.f reference/ssbgvd.f reference/zlaset.f -reference/cla_lin_berr.f reference/dlaed9.f reference/iparmq.f reference/ssbgv.f reference/zlasr.f -reference/clals0.f reference/dlaeda.f reference/izmax1.f reference/ssbgvx.f reference/zlassq.f -reference/clalsa.f reference/dlaein.f reference/lsamen.f reference/ssbtrd.f reference/zlaswp.f -reference/clalsd.f reference/dlaev2.f reference/sbdsdc.f reference/ssfrk.f reference/zla_syamv.f -reference/clangb.f reference/dlaexc.f reference/sbdsqr.f reference/sspcon.f reference/zlasyf.f -reference/clange.f reference/dlag2.f reference/scsum1.f reference/sspevd.f reference/zla_syrcond_c.f -reference/clangt.f reference/dlag2s.f reference/sdisna.f reference/sspev.f reference/zla_syrcond_x.f -reference/clanhb.f reference/dla_gbamv.f reference/sgbbrd.f reference/sspevx.f reference/zla_syrfsx_extended.f -reference/clanhe.f reference/dla_gbrcond.f reference/sgbcon.f reference/sspgst.f reference/zla_syrpvgrw.f -reference/clanhf.f reference/dla_gbrfsx_extended.f reference/sgbequb.f reference/sspgvd.f reference/zlat2c.f -reference/clanhp.f reference/dla_gbrpvgrw.f reference/sgbequ.f reference/sspgv.f reference/zlatbs.f -reference/clanhs.f reference/dla_geamv.f reference/sgbrfs.f reference/sspgvx.f reference/zlatdf.f -reference/clanht.f reference/dla_gercond.f reference/sgbrfsx.f reference/ssprfs.f reference/zlatps.f -reference/clansb.f reference/dla_gerfsx_extended.f reference/sgbsv.f reference/sspsv.f reference/zlatrd.f -reference/clansp.f reference/dlags2.f reference/sgbsvx.f reference/sspsvx.f reference/zlatrs.f -reference/clansy.f reference/dlagtf.f reference/sgbsvxx.f reference/ssptrd.f reference/zlatrz.f -reference/clantb.f reference/dlagtm.f reference/sgbtf2.f reference/ssptrf.f reference/zlatzm.f -reference/clantp.f reference/dlagts.f reference/sgbtrf.f reference/ssptri.f reference/zlauu2.f -reference/clantr.f reference/dlagv2.f reference/sgbtrs.f reference/ssptrs.f reference/zlauum.f -reference/clapll.f reference/dlahqr.f reference/sgebak.f reference/sstebz.f reference/zla_wwaddw.f -reference/clapmt.f reference/dlahr2.f reference/sgebal.f reference/sstedc.f reference/zpbcon.f -reference/cla_porcond_c.f reference/dlahrd.f reference/sgebd2.f reference/sstegr.f reference/zpbequ.f -reference/cla_porcond_x.f reference/dlaic1.f reference/sgebrd.f reference/sstein.f reference/zpbrfs.f -reference/cla_porfsx_extended.f reference/dlaisnan.f reference/sgecon.f reference/sstemr.f reference/zpbstf.f -reference/cla_porpvgrw.f reference/dla_lin_berr.f reference/sgeequb.f reference/ssteqr.f reference/zpbsv.f -reference/claqgb.f reference/dlaln2.f reference/sgeequ.f reference/ssterf.f reference/zpbsvx.f -reference/claqge.f reference/dlals0.f reference/sgees.f reference/sstevd.f reference/zpbtf2.f -reference/claqhb.f reference/dlalsa.f reference/sgeesx.f reference/sstev.f reference/zpbtrf.f -reference/claqhe.f reference/dlalsd.f reference/sgeev.f reference/sstevr.f reference/zpbtrs.f -reference/claqhp.f reference/dlamrg.f reference/sgeevx.f reference/sstevx.f reference/zpftrf.f -reference/claqp2.f reference/dlaneg.f reference/sgegs.f reference/ssycon.f reference/zpftri.f -reference/claqps.f reference/dlangb.f reference/sgegv.f reference/ssyequb.f reference/zpftrs.f -reference/claqr0.f reference/dlange.f reference/sgehd2.f reference/ssyevd.f reference/zpocon.f -reference/claqr1.f reference/dlangt.f reference/sgehrd.f reference/zpoequb.f -reference/claqr2.f reference/dlanhs.f reference/sgejsv.f reference/ssyevr.f reference/zpoequ.f -reference/claqr3.f reference/dlansb.f reference/sgelq2.f reference/ssyevx.f reference/zporfs.f -reference/claqr4.f reference/dlansf.f reference/sgelqf.f reference/ssygs2.f reference/zporfsx.f -reference/claqr5.f reference/dlansp.f reference/sgelsd.f reference/ssygst.f reference/zposv.f -reference/claqsb.f reference/dlanst.f reference/sgels.f reference/ssygvd.f reference/zposvx.f -reference/claqsp.f reference/dlansy.f reference/sgelss.f reference/ssygv.f reference/zposvxx.f -reference/claqsy.f reference/dlantb.f reference/sgelsx.f reference/ssygvx.f reference/zpotf2.f -reference/clar1v.f reference/dlantp.f reference/sgelsy.f reference/ssyrfs.f -reference/clar2v.f reference/dlantr.f reference/sgeql2.f reference/ssyrfsx.f reference/zpotri.f -reference/clarcm.f reference/dlanv2.f reference/sgeqlf.f reference/ssysv.f -reference/clarfb.f reference/dlapll.f reference/sgeqp3.f reference/ssysvx.f reference/zppcon.f -reference/clarf.f reference/dlapmt.f reference/sgeqpf.f reference/ssysvxx.f reference/zppequ.f -reference/clarfg.f reference/dla_porcond.f reference/sgeqr2.f reference/ssytd2.f reference/zpprfs.f -reference/clarfp.f reference/dla_porfsx_extended.f reference/sgeqrf.f reference/ssytf2.f reference/zppsv.f -reference/clarft.f reference/dla_porpvgrw.f reference/sgerfs.f reference/ssytrd.f reference/zppsvx.f -reference/clarfx.f reference/dlapy2.f reference/sgerfsx.f reference/ssytrf.f reference/zpptrf.f -reference/clargv.f reference/dlapy3.f reference/sgerq2.f reference/ssytri.f reference/zpptri.f -reference/clarnv.f reference/dlaqgb.f reference/sgerqf.f reference/ssytrs.f reference/zpptrs.f -reference/cla_rpvgrw.f reference/dlaqge.f reference/sgesc2.f reference/stbcon.f reference/zpstf2.f -reference/clarrv.f reference/dlaqp2.f reference/sgesdd.f reference/stbrfs.f reference/zpstrf.f -reference/clarscl2.f reference/dlaqps.f reference/sgesvd.f reference/stbtrs.f reference/zptcon.f -reference/clartg.f reference/dlaqr0.f reference/sgesv.f reference/stfsm.f reference/zpteqr.f -reference/clartv.f reference/dlaqr1.f reference/sgesvj.f reference/stftri.f reference/zptrfs.f -reference/clarzb.f reference/dlaqr2.f reference/sgesvx.f reference/stfttp.f reference/zptsv.f -reference/clarz.f reference/dlaqr3.f reference/sgesvxx.f reference/stfttr.f reference/zptsvx.f -reference/clarzt.f reference/dlaqr4.f reference/sgetc2.f reference/stgevc.f reference/zpttrf.f -reference/clascl2.f reference/dlaqr5.f reference/sgetf2.f reference/stgex2.f reference/zpttrs.f -reference/clascl.f reference/dlaqsb.f -reference/stgexc.f reference/zptts2.f -reference/claset.f reference/dlaqsp.f reference/sgetri.f reference/stgsen.f reference/zrot.f -reference/clasr.f reference/dlaqsy.f reference/stgsja.f reference/zspcon.f -reference/classq.f reference/dlaqtr.f reference/sggbak.f reference/stgsna.f reference/zspmv.f -reference/claswp.f reference/dlar1v.f reference/sggbal.f reference/stgsy2.f reference/zspr.f -reference/cla_syamv.f reference/dlar2v.f reference/sgges.f reference/stgsyl.f reference/zsprfs.f -reference/clasyf.f reference/dlarfb.f reference/sggesx.f reference/stpcon.f reference/zspsv.f -reference/cla_syrcond_c.f reference/dlarf.f reference/sggev.f reference/stprfs.f reference/zspsvx.f -reference/cla_syrcond_x.f reference/dlarfg.f reference/sggevx.f reference/stptri.f reference/zsptrf.f -reference/cla_syrfsx_extended.f reference/dlarfp.f reference/sggglm.f reference/stptrs.f reference/zsptri.f -reference/cla_syrpvgrw.f reference/dlarft.f reference/sgghrd.f reference/stpttf.f reference/zsptrs.f -reference/clatbs.f reference/dlarfx.f reference/sgglse.f reference/stpttr.f reference/zstedc.f -reference/clatdf.f reference/dlargv.f reference/sggqrf.f reference/strcon.f reference/zstegr.f -reference/clatps.f reference/dlarnv.f reference/sggrqf.f reference/strevc.f reference/zstein.f -reference/clatrd.f reference/dla_rpvgrw.f reference/sggsvd.f reference/strexc.f reference/zstemr.f -reference/clatrs.f reference/dlarra.f reference/sggsvp.f reference/strrfs.f reference/zsteqr.f -reference/clatrz.f reference/dlarrb.f reference/sgsvj0.f reference/strsen.f reference/zsycon.f -reference/clatzm.f reference/dlarrc.f reference/sgsvj1.f reference/strsna.f reference/zsyequb.f -reference/clauu2.f reference/dlarrd.f reference/sgtcon.f reference/strsyl.f reference/zsymv.f -reference/clauum.f reference/dlarre.f reference/sgtrfs.f reference/strti2.f reference/zsyr.f -reference/cla_wwaddw.f reference/dlarrf.f reference/sgtsv.f reference/strtri.f reference/zsyrfs.f -reference/cpbcon.f reference/dlarrj.f reference/sgtsvx.f reference/strtrs.f reference/zsyrfsx.f -reference/cpbequ.f reference/dlarrk.f reference/sgttrf.f reference/strttf.f reference/zsysv.f -reference/cpbrfs.f reference/dlarrr.f reference/sgttrs.f reference/strttp.f reference/zsysvx.f -reference/cpbstf.f reference/dlarrv.f reference/sgtts2.f reference/stzrqf.f reference/zsysvxx.f -reference/cpbsv.f reference/dlarscl2.f reference/shgeqz.f reference/stzrzf.f reference/zsytf2.f -reference/cpbsvx.f reference/dlartg.f reference/shsein.f reference/xerbla_array.f reference/zsytrf.f -reference/cpbtf2.f reference/dlartv.f reference/shseqr.f reference/xerbla.f reference/zsytri.f -reference/cpbtrf.f reference/dlaruv.f reference/sisnan.f reference/zbdsqr.f reference/zsytrs.f -reference/cpbtrs.f reference/dlarzb.f reference/slabad.f reference/zcgesv.f reference/ztbcon.f -reference/cpftrf.f reference/dlarz.f reference/slabrd.f reference/zcposv.f reference/ztbrfs.f -reference/cpftri.f reference/dlarzt.f reference/slacn2.f reference/zdrscl.f reference/ztbtrs.f -reference/cpftrs.f reference/dlas2.f reference/slacon.f reference/zgbbrd.f reference/ztfsm.f -reference/cpocon.f reference/dlascl2.f reference/slacpy.f reference/zgbcon.f reference/ztftri.f -reference/cpoequb.f reference/dlascl.f reference/sladiv.f reference/zgbequb.f reference/ztfttp.f -reference/cpoequ.f reference/dlasd0.f reference/slae2.f reference/zgbequ.f reference/ztfttr.f -reference/cporfs.f reference/dlasd1.f reference/slaebz.f reference/zgbrfs.f reference/ztgevc.f -reference/cporfsx.f reference/dlasd2.f reference/slaed0.f reference/zgbrfsx.f reference/ztgex2.f -reference/cposv.f reference/dlasd3.f reference/slaed1.f reference/zgbsv.f reference/ztgexc.f -reference/cposvx.f reference/dlasd4.f reference/slaed2.f reference/zgbsvx.f reference/ztgsen.f -reference/cposvxx.f reference/dlasd5.f reference/slaed3.f reference/zgbsvxx.f reference/ztgsja.f -reference/cpotf2.f reference/dlasd6.f reference/slaed4.f reference/zgbtf2.f reference/ztgsna.f -reference/dlasd7.f reference/slaed5.f reference/zgbtrf.f reference/ztgsy2.f -reference/cpotri.f reference/dlasd8.f reference/slaed6.f reference/zgbtrs.f reference/ztgsyl.f - reference/dlasda.f reference/slaed7.f reference/zgebak.f reference/ztpcon.f -reference/cppcon.f reference/dlasdq.f reference/slaed8.f reference/zgebal.f reference/ztprfs.f -reference/cppequ.f reference/dlasdt.f reference/slaed9.f reference/zgebd2.f reference/ztptri.f -reference/cpprfs.f reference/dlaset.f reference/slaeda.f reference/zgebrd.f reference/ztptrs.f -reference/cppsv.f reference/dlasq1.f reference/slaein.f reference/zgecon.f reference/ztpttf.f -reference/cppsvx.f reference/dlasq2.f reference/slaev2.f reference/zgeequb.f reference/ztpttr.f -reference/cpptrf.f reference/dlasq3.f reference/slaexc.f reference/zgeequ.f reference/ztrcon.f -reference/cpptri.f reference/dlasq4.f reference/slag2d.f reference/zgees.f reference/ztrevc.f -reference/cpptrs.f reference/dlasq5.f reference/slag2.f reference/zgeesx.f reference/ztrexc.f -reference/cpstf2.f reference/dlasq6.f reference/sla_gbamv.f reference/zgeev.f reference/ztrrfs.f -reference/cpstrf.f reference/dlasr.f reference/sla_gbrcond.f reference/zgeevx.f reference/ztrsen.f -reference/cptcon.f reference/dlasrt.f reference/sla_gbrfsx_extended.f reference/zgegs.f reference/ztrsna.f -reference/cpteqr.f reference/dlassq.f reference/sla_gbrpvgrw.f reference/zgegv.f reference/ztrsyl.f -reference/cptrfs.f reference/dlasv2.f reference/sla_geamv.f reference/zgehd2.f reference/ztrti2.f -reference/cptsv.f reference/dlaswp.f reference/sla_gercond.f reference/zgehrd.f reference/ztrtri.f -reference/cptsvx.f reference/dlasy2.f reference/sla_gerfsx_extended.f reference/zgelq2.f reference/ztrtrs.f -reference/cpttrf.f reference/dla_syamv.f reference/slags2.f reference/zgelqf.f reference/ztrttf.f -reference/cpttrs.f reference/dlasyf.f reference/slagtf.f reference/zgelsd.f reference/ztrttp.f -reference/cptts2.f reference/dla_syrcond.f reference/slagtm.f reference/zgels.f reference/ztzrqf.f -reference/crot.f reference/dla_syrfsx_extended.f reference/slagts.f reference/zgelss.f reference/ztzrzf.f -reference/cspcon.f reference/dla_syrpvgrw.f reference/slagv2.f reference/zgelsx.f reference/zung2l.f -reference/cspmv.f reference/dlat2s.f reference/slahqr.f reference/zgelsy.f reference/zung2r.f -reference/cspr.f reference/dlatbs.f reference/slahr2.f reference/zgeql2.f reference/zungbr.f -reference/csprfs.f reference/dlatdf.f reference/slahrd.f reference/zgeqlf.f reference/zunghr.f -reference/cspsv.f reference/dlatps.f reference/slaic1.f reference/zgeqp3.f reference/zungl2.f -reference/cspsvx.f reference/dlatrd.f reference/slaisnan.f reference/zgeqpf.f reference/zunglq.f -reference/csptrf.f reference/dlatrs.f reference/sla_lin_berr.f reference/zgeqr2.f reference/zungql.f -reference/csptri.f reference/dlatrz.f reference/slaln2.f reference/zgeqrf.f reference/zungqr.f -reference/csptrs.f reference/dlatzm.f reference/slals0.f reference/zgerfs.f reference/zungr2.f -reference/csrscl.f reference/dlauu2.f reference/slalsa.f reference/zgerfsx.f reference/zungrq.f -reference/cstedc.f reference/dlauum.f reference/slalsd.f reference/zgerq2.f reference/zungtr.f -reference/cstegr.f reference/dla_wwaddw.f reference/slamrg.f reference/zgerqf.f reference/zunm2l.f -reference/cstein.f reference/dopgtr.f reference/slaneg.f reference/zgesc2.f reference/zunm2r.f -reference/cstemr.f reference/dopmtr.f reference/slangb.f reference/zgesdd.f reference/zunmbr.f -reference/csteqr.f reference/dorg2l.f reference/slange.f reference/zgesvd.f reference/zunmhr.f -reference/csycon.f reference/dorg2r.f reference/slangt.f reference/zgesv.f reference/zunml2.f -reference/csyequb.f reference/dorgbr.f reference/slanhs.f reference/zgesvx.f reference/zunmlq.f -reference/csymv.f reference/dorghr.f reference/slansb.f reference/zgesvxx.f reference/zunmql.f -reference/csyr.f reference/dorgl2.f reference/slansf.f reference/zgetc2.f reference/zunmqr.f -reference/csyrfs.f reference/dorglq.f reference/slansp.f reference/zgetf2.f reference/zunmr2.f -reference/csyrfsx.f reference/dorgql.f reference/slanst.f -reference/zunmr3.f -reference/csysv.f reference/dorgqr.f reference/slansy.f reference/zgetri.f reference/zunmrq.f -reference/csysvx.f reference/dorgr2.f reference/slantb.f reference/zunmrz.f -reference/csysvxx.f reference/dorgrq.f reference/slantp.f reference/zggbak.f reference/zunmtr.f -reference/csytf2.f reference/dorgtr.f reference/slantr.f reference/zggbal.f reference/zupgtr.f -reference/csytrf.f reference/dorm2l.f reference/slanv2.f reference/zgges.f reference/zupmtr.f -reference/csytri.f reference/dorm2r.f reference/slapll.f reference/zggesx.f -reference/csytrs.f reference/dormbr.f reference/slapmt.f reference/zggev.f -reference/ctbcon.f reference/dormhr.f reference/sla_porcond.f reference/zggevx.f -) -endif() +option(EIGEN_ENABLE_LAPACK_TESTS OFF "Enbale the Lapack unit tests") + +if(EIGEN_ENABLE_LAPACK_TESTS) + + get_filename_component(eigen_full_path_to_reference_lapack "./reference/" ABSOLUTE) + if(NOT EXISTS ${eigen_full_path_to_reference_lapack}) + # Download lapack and install sources and testing at the right place + message(STATUS "Download lapack_addons_3.4.1.tgz...") + + file(DOWNLOAD "http://downloads.tuxfamily.org/eigen/lapack_addons_3.4.1.tgz" + "${CMAKE_CURRENT_SOURCE_DIR}/lapack_addons_3.4.1.tgz" + INACTIVITY_TIMEOUT 15 + TIMEOUT 240 + STATUS download_status + EXPECTED_MD5 5758ce55afcf79da98de8b9de1615ad5 + SHOW_PROGRESS) + + message(STATUS ${download_status}) + list(GET download_status 0 download_status_num) + set(download_status_num 0) + if(download_status_num EQUAL 0) + message(STATUS "Setup lapack reference and lapack unit tests") + execute_process(COMMAND tar xzf "lapack_addons_3.4.1.tgz" WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}) + else() + message(STATUS "Download of lapack_addons_3.4.1.tgz failed, LAPACK unit tests wont be enabled") + set(EIGEN_ENABLE_LAPACK_TESTS false) + endif() + + endif() + + get_filename_component(eigen_full_path_to_reference_lapack "./reference/" ABSOLUTE) + if(EXISTS ${eigen_full_path_to_reference_lapack}) + set(EigenLapack_funcfilenames + ssyev.f dsyev.f csyev.f zsyev.f + spotrf.f dpotrf.f cpotrf.f zpotrf.f + spotrs.f dpotrs.f cpotrs.f zpotrs.f + sgetrf.f dgetrf.f cgetrf.f zgetrf.f + sgetrs.f dgetrs.f cgetrs.f zgetrs.f) + + FILE(GLOB ReferenceLapack_SRCS0 RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} "reference/*.f") + foreach(filename1 IN LISTS ReferenceLapack_SRCS0) + string(REPLACE "reference/" "" filename ${filename1}) + list(FIND EigenLapack_SRCS ${filename} id1) + list(FIND EigenLapack_funcfilenames ${filename} id2) + if((id1 EQUAL -1) AND (id2 EQUAL -1)) + set(ReferenceLapack_SRCS ${ReferenceLapack_SRCS} reference/${filename}) + endif() + endforeach() + endif() + + +endif(EIGEN_ENABLE_LAPACK_TESTS) endif(EIGEN_Fortran_COMPILER_WORKS) -add_library(eigen_lapack_static ${EigenLapack_SRCS}) +add_library(eigen_lapack_static ${EigenLapack_SRCS} ${ReferenceLapack_SRCS}) add_library(eigen_lapack SHARED ${EigenLapack_SRCS}) target_link_libraries(eigen_lapack eigen_blas) @@ -384,5 +104,343 @@ install(TARGETS eigen_lapack eigen_lapack_static LIBRARY DESTINATION lib ARCHIVE DESTINATION lib) -# add_subdirectory(testing) + + +get_filename_component(eigen_full_path_to_testing_lapack "./testing/" ABSOLUTE) +if(EXISTS ${eigen_full_path_to_testing_lapack}) + + # The following comes from lapack/TESTING/CMakeLists.txt + # Get Python + find_package(PythonInterp) + message(STATUS "Looking for Python found - ${PYTHONINTERP_FOUND}") + if (PYTHONINTERP_FOUND) + message(STATUS "Using Python version ${PYTHON_VERSION_STRING}") + endif() + + set(LAPACK_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}) + set(LAPACK_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}) + set(BUILD_SINGLE true) + set(BUILD_DOUBLE true) + set(BUILD_COMPLEX true) + set(BUILD_COMPLEX16E true) + + if(MSVC_VERSION) +# string(REPLACE "/STACK:10000000" "/STACK:900000000000000000" +# CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS}") + string(REGEX REPLACE "(.*)/STACK:(.*) (.*)" "\\1/STACK:900000000000000000 \\3" + CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS}") + endif() + add_subdirectory(testing/MATGEN) + add_subdirectory(testing/LIN) + add_subdirectory(testing/EIG) + macro(add_lapack_test output input target) + set(TEST_INPUT "${LAPACK_SOURCE_DIR}/testing/${input}") + set(TEST_OUTPUT "${LAPACK_BINARY_DIR}/testing/${output}") + get_target_property(TEST_LOC ${target} LOCATION) + string(REPLACE "." "_" input_name ${input}) + set(testName "${target}_${input_name}") + if(EXISTS "${TEST_INPUT}") + add_test(LAPACK-${testName} "${CMAKE_COMMAND}" + -DTEST=${TEST_LOC} + -DINPUT=${TEST_INPUT} + -DOUTPUT=${TEST_OUTPUT} + -DINTDIR=${CMAKE_CFG_INTDIR} + -P "${LAPACK_SOURCE_DIR}/testing/runtest.cmake") + endif() + endmacro(add_lapack_test) + + if (BUILD_SINGLE) + add_lapack_test(stest.out stest.in xlintsts) + # + # ======== SINGLE RFP LIN TESTS ======================== + add_lapack_test(stest_rfp.out stest_rfp.in xlintstrfs) + # + # + # ======== SINGLE EIG TESTS =========================== + # + + add_lapack_test(snep.out nep.in xeigtsts) + + + add_lapack_test(ssep.out sep.in xeigtsts) + + + add_lapack_test(ssvd.out svd.in xeigtsts) + + + add_lapack_test(sec.out sec.in xeigtsts) + + + add_lapack_test(sed.out sed.in xeigtsts) + + + add_lapack_test(sgg.out sgg.in xeigtsts) + + + add_lapack_test(sgd.out sgd.in xeigtsts) + + + add_lapack_test(ssb.out ssb.in xeigtsts) + + + add_lapack_test(ssg.out ssg.in xeigtsts) + + + add_lapack_test(sbal.out sbal.in xeigtsts) + + + add_lapack_test(sbak.out sbak.in xeigtsts) + + + add_lapack_test(sgbal.out sgbal.in xeigtsts) + + + add_lapack_test(sgbak.out sgbak.in xeigtsts) + + + add_lapack_test(sbb.out sbb.in xeigtsts) + + + add_lapack_test(sglm.out glm.in xeigtsts) + + + add_lapack_test(sgqr.out gqr.in xeigtsts) + + + add_lapack_test(sgsv.out gsv.in xeigtsts) + + + add_lapack_test(scsd.out csd.in xeigtsts) + + + add_lapack_test(slse.out lse.in xeigtsts) + endif() + + if (BUILD_DOUBLE) + # + # ======== DOUBLE LIN TESTS =========================== + add_lapack_test(dtest.out dtest.in xlintstd) + # + # ======== DOUBLE RFP LIN TESTS ======================== + add_lapack_test(dtest_rfp.out dtest_rfp.in xlintstrfd) + # + # ======== DOUBLE EIG TESTS =========================== + + add_lapack_test(dnep.out nep.in xeigtstd) + + + add_lapack_test(dsep.out sep.in xeigtstd) + + + add_lapack_test(dsvd.out svd.in xeigtstd) + + + add_lapack_test(dec.out dec.in xeigtstd) + + + add_lapack_test(ded.out ded.in xeigtstd) + + + add_lapack_test(dgg.out dgg.in xeigtstd) + + + add_lapack_test(dgd.out dgd.in xeigtstd) + + + add_lapack_test(dsb.out dsb.in xeigtstd) + + + add_lapack_test(dsg.out dsg.in xeigtstd) + + + add_lapack_test(dbal.out dbal.in xeigtstd) + + + add_lapack_test(dbak.out dbak.in xeigtstd) + + + add_lapack_test(dgbal.out dgbal.in xeigtstd) + + + add_lapack_test(dgbak.out dgbak.in xeigtstd) + + + add_lapack_test(dbb.out dbb.in xeigtstd) + + + add_lapack_test(dglm.out glm.in xeigtstd) + + + add_lapack_test(dgqr.out gqr.in xeigtstd) + + + add_lapack_test(dgsv.out gsv.in xeigtstd) + + + add_lapack_test(dcsd.out csd.in xeigtstd) + + + add_lapack_test(dlse.out lse.in xeigtstd) + endif() + + if (BUILD_COMPLEX) + add_lapack_test(ctest.out ctest.in xlintstc) + # + # ======== COMPLEX RFP LIN TESTS ======================== + add_lapack_test(ctest_rfp.out ctest_rfp.in xlintstrfc) + # + # ======== COMPLEX EIG TESTS =========================== + + add_lapack_test(cnep.out nep.in xeigtstc) + + + add_lapack_test(csep.out sep.in xeigtstc) + + + add_lapack_test(csvd.out svd.in xeigtstc) + + + add_lapack_test(cec.out cec.in xeigtstc) + + + add_lapack_test(ced.out ced.in xeigtstc) + + + add_lapack_test(cgg.out cgg.in xeigtstc) + + + add_lapack_test(cgd.out cgd.in xeigtstc) + + + add_lapack_test(csb.out csb.in xeigtstc) + + + add_lapack_test(csg.out csg.in xeigtstc) + + + add_lapack_test(cbal.out cbal.in xeigtstc) + + + add_lapack_test(cbak.out cbak.in xeigtstc) + + + add_lapack_test(cgbal.out cgbal.in xeigtstc) + + + add_lapack_test(cgbak.out cgbak.in xeigtstc) + + + add_lapack_test(cbb.out cbb.in xeigtstc) + + + add_lapack_test(cglm.out glm.in xeigtstc) + + + add_lapack_test(cgqr.out gqr.in xeigtstc) + + + add_lapack_test(cgsv.out gsv.in xeigtstc) + + + add_lapack_test(ccsd.out csd.in xeigtstc) + + + add_lapack_test(clse.out lse.in xeigtstc) + endif() + + if (BUILD_COMPLEX16) + # + # ======== COMPLEX16 LIN TESTS ======================== + add_lapack_test(ztest.out ztest.in xlintstz) + # + # ======== COMPLEX16 RFP LIN TESTS ======================== + add_lapack_test(ztest_rfp.out ztest_rfp.in xlintstrfz) + # + # ======== COMPLEX16 EIG TESTS =========================== + + add_lapack_test(znep.out nep.in xeigtstz) + + + add_lapack_test(zsep.out sep.in xeigtstz) + + + add_lapack_test(zsvd.out svd.in xeigtstz) + + + add_lapack_test(zec.out zec.in xeigtstz) + + + add_lapack_test(zed.out zed.in xeigtstz) + + + add_lapack_test(zgg.out zgg.in xeigtstz) + + + add_lapack_test(zgd.out zgd.in xeigtstz) + + + add_lapack_test(zsb.out zsb.in xeigtstz) + + + add_lapack_test(zsg.out zsg.in xeigtstz) + + + add_lapack_test(zbal.out zbal.in xeigtstz) + + + add_lapack_test(zbak.out zbak.in xeigtstz) + + + add_lapack_test(zgbal.out zgbal.in xeigtstz) + + + add_lapack_test(zgbak.out zgbak.in xeigtstz) + + + add_lapack_test(zbb.out zbb.in xeigtstz) + + + add_lapack_test(zglm.out glm.in xeigtstz) + + + add_lapack_test(zgqr.out gqr.in xeigtstz) + + + add_lapack_test(zgsv.out gsv.in xeigtstz) + + + add_lapack_test(zcsd.out csd.in xeigtstz) + + + add_lapack_test(zlse.out lse.in xeigtstz) + endif() + + + if (BUILD_SIMPLE) + if (BUILD_DOUBLE) + # + # ======== SINGLE-DOUBLE PROTO LIN TESTS ============== + add_lapack_test(dstest.out dstest.in xlintstds) + endif() + endif() + + + if (BUILD_COMPLEX) + if (BUILD_COMPLEX16) + # + # ======== COMPLEX-COMPLEX16 LIN TESTS ======================== + add_lapack_test(zctest.out zctest.in xlintstzc) + endif() + endif() + + # ============================================================================== + + execute_process(COMMAND ${CMAKE_COMMAND} -E copy ${LAPACK_SOURCE_DIR}/testing/lapack_testing.py ${LAPACK_BINARY_DIR}) + add_test( + NAME LAPACK_Test_Summary + WORKING_DIRECTORY ${LAPACK_BINARY_DIR} + COMMAND ${PYTHON_EXECUTABLE} "lapack_testing.py" + ) + +endif() diff --git a/lapack/dsecnd_NONE.f b/lapack/dsecnd_NONE.f new file mode 100644 index 000000000..61a8dff13 --- /dev/null +++ b/lapack/dsecnd_NONE.f @@ -0,0 +1,52 @@ +*> \brief \b DSECND returns nothing +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* DOUBLE PRECISION FUNCTION DSECND( ) +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DSECND returns nothing instead of returning the user time for a process in seconds. +*> If you are using that routine, it means that neither EXTERNAL ETIME, +*> EXTERNAL ETIME_, INTERNAL ETIME, INTERNAL CPU_TIME is available on +*> your machine. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup auxOTHERauxiliary +* +* ===================================================================== + DOUBLE PRECISION FUNCTION DSECND( ) +* +* -- LAPACK auxiliary routine (version 3.4.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* ===================================================================== +* + DSECND = 0.0D+0 + RETURN +* +* End of DSECND +* + END diff --git a/lapack/second_NONE.f b/lapack/second_NONE.f new file mode 100644 index 000000000..d3e6d3319 --- /dev/null +++ b/lapack/second_NONE.f @@ -0,0 +1,52 @@ +*> \brief \b SECOND returns nothing +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* REAL FUNCTION SECOND( ) +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> SECOND returns nothing instead of returning the user time for a process in seconds. +*> If you are using that routine, it means that neither EXTERNAL ETIME, +*> EXTERNAL ETIME_, INTERNAL ETIME, INTERNAL CPU_TIME is available on +*> your machine. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup auxOTHERauxiliary +* +* ===================================================================== + REAL FUNCTION SECOND( ) +* +* -- LAPACK auxiliary routine (version 3.4.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* ===================================================================== +* + SECOND = 0.0E+0 + RETURN +* +* End of SECOND +* + END diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index be9617d85..580e2ef5e 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -153,7 +153,7 @@ ei_add_test(diagonal) ei_add_test(miscmatrices) ei_add_test(commainitializer) ei_add_test(smallvectors) -ei_add_test(map) +ei_add_test(mapped_matrix) ei_add_test(mapstride) ei_add_test(mapstaticmethods) ei_add_test(array) diff --git a/test/adjoint.cpp b/test/adjoint.cpp index b63e843c6..ea36f7841 100644 --- a/test/adjoint.cpp +++ b/test/adjoint.cpp @@ -28,6 +28,7 @@ template<> struct adjoint_specific<false> { template<typename Vec, typename Mat, typename Scalar> static void run(const Vec& v1, const Vec& v2, Vec& v3, const Mat& square, Scalar s1, Scalar s2) { typedef typename NumTraits<Scalar>::Real RealScalar; + using std::abs; RealScalar ref = NumTraits<Scalar>::IsInteger ? RealScalar(0) : (std::max)((s1 * v1 + s2 * v2).norm(),v3.norm()); VERIFY(test_isApproxWithRef((s1 * v1 + s2 * v2).dot(v3), numext::conj(s1) * v1.dot(v3) + numext::conj(s2) * v2.dot(v3), ref)); @@ -44,7 +45,7 @@ template<> struct adjoint_specific<false> { // check compatibility of dot and adjoint ref = NumTraits<Scalar>::IsInteger ? 0 : (std::max)((std::max)(v1.norm(),v2.norm()),(std::max)((square * v2).norm(),(square.adjoint() * v1).norm())); - VERIFY(test_isApproxWithRef(v1.dot(square * v2), (square.adjoint() * v1).dot(v2), ref)); + VERIFY(internal::isMuchSmallerThan(abs(v1.dot(square * v2) - (square.adjoint() * v1).dot(v2)), ref, test_precision<Scalar>())); // check that Random().normalized() works: tricky as the random xpr must be evaluated by // normalized() in order to produce a consistent result. diff --git a/test/array.cpp b/test/array.cpp index 631a29432..8960e49f8 100644 --- a/test/array.cpp +++ b/test/array.cpp @@ -152,6 +152,7 @@ template<typename ArrayType> void comparisons(const ArrayType& m) template<typename ArrayType> void array_real(const ArrayType& m) { using std::abs; + using std::sqrt; typedef typename ArrayType::Index Index; typedef typename ArrayType::Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; @@ -189,7 +190,7 @@ template<typename ArrayType> void array_real(const ArrayType& m) if(!NumTraits<Scalar>::IsComplex) VERIFY_IS_APPROX(numext::real(m1), m1); - VERIFY((m1.abs().log() == log(abs(m1))).all()); + VERIFY_IS_APPROX(m1.abs().log() , log(abs(m1))); // VERIFY_IS_APPROX(m1.exp(), std::exp(m1)); VERIFY_IS_APPROX(m1.exp() * m2.exp(), exp(m1+m2)); diff --git a/test/cholesky.cpp b/test/cholesky.cpp index 38862924a..378525a83 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -328,4 +328,5 @@ void test_cholesky() CALL_SUBTEST_9( LDLT<MatrixXf>(10) ); TEST_SET_BUT_UNUSED_VARIABLE(s) + TEST_SET_BUT_UNUSED_VARIABLE(nb_temporaries) } diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp index 5c6ecd875..06a6a8654 100644 --- a/test/eigensolver_selfadjoint.cpp +++ b/test/eigensolver_selfadjoint.cpp @@ -99,6 +99,15 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m) // FIXME tridiag.matrixQ().adjoint() does not work VERIFY_IS_APPROX(MatrixType(symmA.template selfadjointView<Lower>()), tridiag.matrixQ() * tridiag.matrixT().eval() * MatrixType(tridiag.matrixQ()).adjoint()); + // Test computation of eigenvalues from tridiagonal matrix + if(rows > 1) + { + SelfAdjointEigenSolver<MatrixType> eiSymmTridiag; + eiSymmTridiag.computeFromTridiagonal(tridiag.matrixT().diagonal(), tridiag.matrixT().diagonal(-1), ComputeEigenvectors); + VERIFY_IS_APPROX(eiSymm.eigenvalues(), eiSymmTridiag.eigenvalues()); + VERIFY_IS_APPROX(tridiag.matrixT(), eiSymmTridiag.eigenvectors().real() * eiSymmTridiag.eigenvalues().asDiagonal() * eiSymmTridiag.eigenvectors().real().transpose()); + } + if (rows > 1) { // Test matrix with NaN diff --git a/test/geo_alignedbox.cpp b/test/geo_alignedbox.cpp index e9fbfddf1..8e36adbe3 100644 --- a/test/geo_alignedbox.cpp +++ b/test/geo_alignedbox.cpp @@ -15,6 +15,10 @@ #include<iostream> using namespace std; +template<typename T> EIGEN_DONT_INLINE +void kill_extra_precision(T& x) { eigen_assert(&x != 0); } + + template<typename BoxType> void alignedbox(const BoxType& _box) { /* this test covers the following files: @@ -36,6 +40,10 @@ template<typename BoxType> void alignedbox(const BoxType& _box) BoxType b0(dim); BoxType b1(VectorType::Random(dim),VectorType::Random(dim)); BoxType b2; + + kill_extra_precision(b1); + kill_extra_precision(p0); + kill_extra_precision(p1); b0.extend(p0); b0.extend(p1); diff --git a/test/geo_eulerangles.cpp b/test/geo_eulerangles.cpp index 26456beee..5445cd81a 100644 --- a/test/geo_eulerangles.cpp +++ b/test/geo_eulerangles.cpp @@ -17,13 +17,16 @@ template<typename Scalar> void check_all_var(const Matrix<Scalar,3,1>& ea) typedef Matrix<Scalar,3,3> Matrix3; typedef Matrix<Scalar,3,1> Vector3; typedef AngleAxis<Scalar> AngleAxisx; + using std::abs; #define VERIFY_EULER(I,J,K, X,Y,Z) { \ Matrix3 m(AngleAxisx(ea[0], Vector3::Unit##X()) * AngleAxisx(ea[1], Vector3::Unit##Y()) * AngleAxisx(ea[2], Vector3::Unit##Z())); \ Vector3 eabis = m.eulerAngles(I,J,K); \ Matrix3 mbis(AngleAxisx(eabis[0], Vector3::Unit##X()) * AngleAxisx(eabis[1], Vector3::Unit##Y()) * AngleAxisx(eabis[2], Vector3::Unit##Z())); \ VERIFY_IS_APPROX(m, mbis); \ - if(I!=K || ea[1]!=0) VERIFY_IS_APPROX(ea, eabis); \ + /* If I==K, and ea[1]==0, then there no unique solution. */ \ + /* The remark apply in the case where I!=K, and |ea[1]| is close to pi/2. */ \ + if( (I!=K || ea[1]!=0) && (I==K || !internal::isApprox(abs(ea[1]),Scalar(M_PI/2),test_precision<Scalar>())) ) VERIFY((ea-eabis).norm() <= test_precision<Scalar>()); \ } VERIFY_EULER(0,1,2, X,Y,Z); VERIFY_EULER(0,1,0, X,Y,X); @@ -41,7 +44,7 @@ template<typename Scalar> void check_all_var(const Matrix<Scalar,3,1>& ea) VERIFY_EULER(2,1,2, Z,Y,Z); } -template<typename Scalar> void eulerangles(void) +template<typename Scalar> void eulerangles() { typedef Matrix<Scalar,3,3> Matrix3; typedef Matrix<Scalar,3,1> Vector3; @@ -60,13 +63,13 @@ template<typename Scalar> void eulerangles(void) ea = m.eulerAngles(0,1,0); check_all_var(ea); - ea = (Array3::Random() + Array3(1,1,0))*M_PI*Array3(0.5,0.5,1); + ea = (Array3::Random() + Array3(1,1,0))*Scalar(M_PI)*Array3(0.5,0.5,1); check_all_var(ea); - ea[2] = ea[0] = internal::random<Scalar>(0,M_PI); + ea[2] = ea[0] = internal::random<Scalar>(0,Scalar(M_PI)); check_all_var(ea); - ea[0] = ea[1] = internal::random<Scalar>(0,M_PI); + ea[0] = ea[1] = internal::random<Scalar>(0,Scalar(M_PI)); check_all_var(ea); ea[1] = 0; diff --git a/test/geo_quaternion.cpp b/test/geo_quaternion.cpp index 06b3af7c1..1694b32c7 100644 --- a/test/geo_quaternion.cpp +++ b/test/geo_quaternion.cpp @@ -31,14 +31,14 @@ template<typename QuatType> void check_slerp(const QuatType& q0, const QuatType& Scalar theta_tot = AA(q1*q0.inverse()).angle(); if(theta_tot>M_PI) - theta_tot = 2.*M_PI-theta_tot; - for(Scalar t=0; t<=1.001; t+=0.1) + theta_tot = Scalar(2.*M_PI)-theta_tot; + for(Scalar t=0; t<=Scalar(1.001); t+=Scalar(0.1)) { QuatType q = q0.slerp(t,q1); Scalar theta = AA(q*q0.inverse()).angle(); VERIFY(abs(q.norm() - 1) < largeEps); if(theta_tot==0) VERIFY(theta_tot==0); - else VERIFY(abs(theta/theta_tot - t) < largeEps); + else VERIFY(abs(theta - t * theta_tot) < largeEps); } } @@ -156,7 +156,7 @@ template<typename Scalar, int Options> void quaternion(void) check_slerp(q1,q2); q1 = AngleAxisx(b, v1.normalized()); - q2 = AngleAxisx(b+M_PI, v1.normalized()); + q2 = AngleAxisx(b+Scalar(M_PI), v1.normalized()); check_slerp(q1,q2); q1 = AngleAxisx(b, v1.normalized()); diff --git a/test/jacobi.cpp b/test/jacobi.cpp index 0e6f1de49..7ccd4124b 100644 --- a/test/jacobi.cpp +++ b/test/jacobi.cpp @@ -74,7 +74,8 @@ void test_jacobi() // complex<float> is really important to test as it is the only way to cover conjugation issues in certain unaligned paths CALL_SUBTEST_6(( jacobi<MatrixXcf, float>(MatrixXcf(r,c)) )); CALL_SUBTEST_6(( jacobi<MatrixXcf, std::complex<float> >(MatrixXcf(r,c)) )); - (void) r; - (void) c; + + TEST_SET_BUT_UNUSED_VARIABLE(r); + TEST_SET_BUT_UNUSED_VARIABLE(c); } } diff --git a/test/main.h b/test/main.h index 0f55fef71..14f0d2f78 100644 --- a/test/main.h +++ b/test/main.h @@ -270,6 +270,7 @@ inline bool test_isApprox(const Type1& a, const Type2& b) // The idea behind this function is to compare the two scalars a and b where // the scalar ref is a hint about the expected order of magnitude of a and b. +// WARNING: the scalar a and b must be positive // Therefore, if for some reason a and b are very small compared to ref, // we won't issue a false negative. // This test could be: abs(a-b) <= eps * ref diff --git a/test/map.cpp b/test/mapped_matrix.cpp index 2b52e4f38..de9dbbde3 100644 --- a/test/map.cpp +++ b/test/mapped_matrix.cpp @@ -114,7 +114,7 @@ template<typename PlainObjectType> void check_const_correctness(const PlainObjec VERIFY( !(Map<ConstPlainObjectType, Aligned>::Flags & LvalueBit) ); } -void test_map() +void test_mapped_matrix() { for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1( map_class_vector(Matrix<float, 1, 1>()) ); diff --git a/test/mapstride.cpp b/test/mapstride.cpp index fe35b9d23..b1dc9de2a 100644 --- a/test/mapstride.cpp +++ b/test/mapstride.cpp @@ -116,7 +116,7 @@ template<int Alignment,typename MatrixType> void map_class_matrix(const MatrixTy void test_mapstride() { for(int i = 0; i < g_repeat; i++) { - EIGEN_UNUSED int maxn = 30; + int maxn = 30; CALL_SUBTEST_1( map_class_vector<Aligned>(Matrix<float, 1, 1>()) ); CALL_SUBTEST_1( map_class_vector<Unaligned>(Matrix<float, 1, 1>()) ); CALL_SUBTEST_2( map_class_vector<Aligned>(Vector4d()) ); @@ -142,5 +142,7 @@ void test_mapstride() CALL_SUBTEST_5( map_class_matrix<Unaligned>(MatrixXi(internal::random<int>(1,maxn),internal::random<int>(1,maxn))) ); CALL_SUBTEST_6( map_class_matrix<Aligned>(MatrixXcd(internal::random<int>(1,maxn),internal::random<int>(1,maxn))) ); CALL_SUBTEST_6( map_class_matrix<Unaligned>(MatrixXcd(internal::random<int>(1,maxn),internal::random<int>(1,maxn))) ); + + TEST_SET_BUT_UNUSED_VARIABLE(maxn); } } diff --git a/test/redux.cpp b/test/redux.cpp index bb65f9461..0d176e500 100644 --- a/test/redux.cpp +++ b/test/redux.cpp @@ -22,7 +22,7 @@ template<typename MatrixType> void matrixRedux(const MatrixType& m) // The entries of m1 are uniformly distributed in [0,1], so m1.prod() is very small. This may lead to test // failures if we underflow into denormals. Thus, we scale so that entires are close to 1. - MatrixType m1_for_prod = MatrixType::Ones(rows, cols) + Scalar(0.2) * m1; + MatrixType m1_for_prod = MatrixType::Ones(rows, cols) + RealScalar(0.2) * m1; VERIFY_IS_MUCH_SMALLER_THAN(MatrixType::Zero(rows, cols).sum(), Scalar(1)); VERIFY_IS_APPROX(MatrixType::Ones(rows, cols).sum(), Scalar(float(rows*cols))); // the float() here to shut up excessive MSVC warning about int->complex conversion being lossy diff --git a/test/sizeof.cpp b/test/sizeof.cpp index 68463c9b6..d9ad35620 100644 --- a/test/sizeof.cpp +++ b/test/sizeof.cpp @@ -13,7 +13,7 @@ template<typename MatrixType> void verifySizeOf(const MatrixType&) { typedef typename MatrixType::Scalar Scalar; if (MatrixType::RowsAtCompileTime!=Dynamic && MatrixType::ColsAtCompileTime!=Dynamic) - VERIFY(sizeof(MatrixType)==sizeof(Scalar)*size_t(MatrixType::SizeAtCompileTime)); + VERIFY(std::ptrdiff_t(sizeof(MatrixType))==std::ptrdiff_t(sizeof(Scalar))*std::ptrdiff_t(MatrixType::SizeAtCompileTime)); else VERIFY(sizeof(MatrixType)==sizeof(Scalar*) + 2 * sizeof(typename MatrixType::Index)); } diff --git a/test/sparse.h b/test/sparse.h index 7e2b98494..a09c65e5f 100644 --- a/test/sparse.h +++ b/test/sparse.h @@ -58,18 +58,18 @@ initSparse(double density, Matrix<Scalar,Dynamic,Dynamic,Opt1>& refMat, SparseMatrix<Scalar,Opt2,Index>& sparseMat, int flags = 0, - std::vector<Vector2i>* zeroCoords = 0, - std::vector<Vector2i>* nonzeroCoords = 0) + std::vector<Matrix<Index,2,1> >* zeroCoords = 0, + std::vector<Matrix<Index,2,1> >* nonzeroCoords = 0) { enum { IsRowMajor = SparseMatrix<Scalar,Opt2,Index>::IsRowMajor }; sparseMat.setZero(); //sparseMat.reserve(int(refMat.rows()*refMat.cols()*density)); sparseMat.reserve(VectorXi::Constant(IsRowMajor ? refMat.rows() : refMat.cols(), int((1.5*density)*(IsRowMajor?refMat.cols():refMat.rows())))); - for(int j=0; j<sparseMat.outerSize(); j++) + for(Index j=0; j<sparseMat.outerSize(); j++) { //sparseMat.startVec(j); - for(int i=0; i<sparseMat.innerSize(); i++) + for(Index i=0; i<sparseMat.innerSize(); i++) { int ai(i), aj(j); if(IsRowMajor) @@ -93,11 +93,11 @@ initSparse(double density, //sparseMat.insertBackByOuterInner(j,i) = v; sparseMat.insertByOuterInner(j,i) = v; if (nonzeroCoords) - nonzeroCoords->push_back(Vector2i(ai,aj)); + nonzeroCoords->push_back(Matrix<Index,2,1> (ai,aj)); } else if (zeroCoords) { - zeroCoords->push_back(Vector2i(ai,aj)); + zeroCoords->push_back(Matrix<Index,2,1> (ai,aj)); } refMat(ai,aj) = v; } @@ -110,8 +110,8 @@ initSparse(double density, Matrix<Scalar,Dynamic,Dynamic, Opt1>& refMat, DynamicSparseMatrix<Scalar, Opt2, Index>& sparseMat, int flags = 0, - std::vector<Vector2i>* zeroCoords = 0, - std::vector<Vector2i>* nonzeroCoords = 0) + std::vector<Matrix<Index,2,1> >* zeroCoords = 0, + std::vector<Matrix<Index,2,1> >* nonzeroCoords = 0) { enum { IsRowMajor = DynamicSparseMatrix<Scalar,Opt2,Index>::IsRowMajor }; sparseMat.setZero(); @@ -142,11 +142,11 @@ initSparse(double density, { sparseMat.insertBackByOuterInner(j,i) = v; if (nonzeroCoords) - nonzeroCoords->push_back(Vector2i(ai,aj)); + nonzeroCoords->push_back(Matrix<Index,2,1> (ai,aj)); } else if (zeroCoords) { - zeroCoords->push_back(Vector2i(ai,aj)); + zeroCoords->push_back(Matrix<Index,2,1> (ai,aj)); } refMat(ai,aj) = v; } diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp index 8fc1904b1..498ecfe29 100644 --- a/test/sparse_basic.cpp +++ b/test/sparse_basic.cpp @@ -14,7 +14,8 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& ref) { typedef typename SparseMatrixType::Index Index; - + typedef Matrix<Index,2,1> Vector2; + const Index rows = ref.rows(); const Index cols = ref.cols(); typedef typename SparseMatrixType::Scalar Scalar; @@ -31,8 +32,8 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re DenseMatrix refMat = DenseMatrix::Zero(rows, cols); DenseVector vec1 = DenseVector::Random(rows); - std::vector<Vector2i> zeroCoords; - std::vector<Vector2i> nonzeroCoords; + std::vector<Vector2> zeroCoords; + std::vector<Vector2> nonzeroCoords; initSparse<Scalar>(density, refMat, m, 0, &zeroCoords, &nonzeroCoords); if (zeroCoords.size()==0 || nonzeroCoords.size()==0) @@ -104,11 +105,11 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re SparseMatrixType m2(rows,cols); if(internal::random<int>()%2) m2.reserve(VectorXi::Constant(m2.outerSize(), 2)); - for (int j=0; j<cols; ++j) + for (Index j=0; j<cols; ++j) { - for (int k=0; k<rows/2; ++k) + for (Index k=0; k<rows/2; ++k) { - int i = internal::random<int>(0,rows-1); + Index i = internal::random<Index>(0,rows-1); if (m1.coeff(i,j)==Scalar(0)) m2.insert(i,j) = m1(i,j) = internal::random<Scalar>(); } @@ -126,8 +127,8 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re m2.reserve(VectorXi::Constant(m2.outerSize(), 2)); for (int k=0; k<rows*cols; ++k) { - int i = internal::random<int>(0,rows-1); - int j = internal::random<int>(0,cols-1); + Index i = internal::random<Index>(0,rows-1); + Index j = internal::random<Index>(0,cols-1); if ((m1.coeff(i,j)==Scalar(0)) && (internal::random<int>()%2)) m2.insert(i,j) = m1(i,j) = internal::random<Scalar>(); else @@ -150,8 +151,8 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re m2.reserve(r); for (int k=0; k<rows*cols; ++k) { - int i = internal::random<int>(0,rows-1); - int j = internal::random<int>(0,cols-1); + Index i = internal::random<Index>(0,rows-1); + Index j = internal::random<Index>(0,cols-1); if (m1.coeff(i,j)==Scalar(0)) m2.insert(i,j) = m1(i,j) = internal::random<Scalar>(); if(mode==3) @@ -167,8 +168,8 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); SparseMatrixType m2(rows, rows); initSparse<Scalar>(density, refMat2, m2); - int j0 = internal::random<int>(0,rows-1); - int j1 = internal::random<int>(0,rows-1); + Index j0 = internal::random<Index>(0,rows-1); + Index j1 = internal::random<Index>(0,rows-1); if(SparseMatrixType::IsRowMajor) VERIFY_IS_APPROX(m2.innerVector(j0), refMat2.row(j0)); else @@ -181,17 +182,17 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re SparseMatrixType m3(rows,rows); m3.reserve(VectorXi::Constant(rows,rows/2)); - for(int j=0; j<rows; ++j) - for(int k=0; k<j; ++k) + for(Index j=0; j<rows; ++j) + for(Index k=0; k<j; ++k) m3.insertByOuterInner(j,k) = k+1; - for(int j=0; j<rows; ++j) + for(Index j=0; j<rows; ++j) { VERIFY(j==numext::real(m3.innerVector(j).nonZeros())); if(j>0) VERIFY(j==numext::real(m3.innerVector(j).lastCoeff())); } m3.makeCompressed(); - for(int j=0; j<rows; ++j) + for(Index j=0; j<rows; ++j) { VERIFY(j==numext::real(m3.innerVector(j).nonZeros())); if(j>0) @@ -210,9 +211,9 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re initSparse<Scalar>(density, refMat2, m2); if(internal::random<float>(0,1)>0.5) m2.makeCompressed(); - int j0 = internal::random<int>(0,rows-2); - int j1 = internal::random<int>(0,rows-2); - int n0 = internal::random<int>(1,rows-(std::max)(j0,j1)); + Index j0 = internal::random<Index>(0,rows-2); + Index j1 = internal::random<Index>(0,rows-2); + Index n0 = internal::random<Index>(1,rows-(std::max)(j0,j1)); if(SparseMatrixType::IsRowMajor) VERIFY_IS_APPROX(m2.innerVectors(j0,n0), refMat2.block(j0,0,n0,cols)); else @@ -300,9 +301,9 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); SparseMatrixType m2(rows, rows); initSparse<Scalar>(density, refMat2, m2); - int j0 = internal::random<int>(0,rows-2); - int j1 = internal::random<int>(0,rows-2); - int n0 = internal::random<int>(1,rows-(std::max)(j0,j1)); + Index j0 = internal::random<Index>(0,rows-2); + Index j1 = internal::random<Index>(0,rows-2); + Index n0 = internal::random<Index>(1,rows-(std::max)(j0,j1)); if(SparseMatrixType::IsRowMajor) VERIFY_IS_APPROX(m2.block(j0,0,n0,cols), refMat2.block(j0,0,n0,cols)); else @@ -315,7 +316,7 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re VERIFY_IS_APPROX(m2.block(0,j0,rows,n0)+m2.block(0,j1,rows,n0), refMat2.block(0,j0,rows,n0)+refMat2.block(0,j1,rows,n0)); - int i = internal::random<int>(0,m2.outerSize()-1); + Index i = internal::random<Index>(0,m2.outerSize()-1); if(SparseMatrixType::IsRowMajor) { m2.innerVector(i) = m2.innerVector(i) * s1; refMat2.row(i) = refMat2.row(i) * s1; @@ -334,10 +335,10 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re refM2.setZero(); int countFalseNonZero = 0; int countTrueNonZero = 0; - for (int j=0; j<m2.outerSize(); ++j) + for (Index j=0; j<m2.outerSize(); ++j) { m2.startVec(j); - for (int i=0; i<m2.innerSize(); ++i) + for (Index i=0; i<m2.innerSize(); ++i) { float x = internal::random<float>(0,1); if (x<0.1) @@ -378,8 +379,8 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re refMat.setZero(); for(int i=0;i<ntriplets;++i) { - int r = internal::random<int>(0,rows-1); - int c = internal::random<int>(0,cols-1); + Index r = internal::random<Index>(0,rows-1); + Index c = internal::random<Index>(0,cols-1); Scalar v = internal::random<Scalar>(); triplets.push_back(TripletType(r,c,v)); refMat(r,c) += v; @@ -448,16 +449,16 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re // test conservative resize { - std::vector< std::pair<int,int> > inc; - inc.push_back(std::pair<int,int>(-3,-2)); - inc.push_back(std::pair<int,int>(0,0)); - inc.push_back(std::pair<int,int>(3,2)); - inc.push_back(std::pair<int,int>(3,0)); - inc.push_back(std::pair<int,int>(0,3)); + std::vector< std::pair<Index,Index> > inc; + inc.push_back(std::pair<Index,Index>(-3,-2)); + inc.push_back(std::pair<Index,Index>(0,0)); + inc.push_back(std::pair<Index,Index>(3,2)); + inc.push_back(std::pair<Index,Index>(3,0)); + inc.push_back(std::pair<Index,Index>(0,3)); for(size_t i = 0; i< inc.size(); i++) { - int incRows = inc[i].first; - int incCols = inc[i].second; + Index incRows = inc[i].first; + Index incCols = inc[i].second; SparseMatrixType m1(rows, cols); DenseMatrix refMat1 = DenseMatrix::Zero(rows, cols); initSparse<Scalar>(density, refMat1, m1); @@ -471,9 +472,9 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re // Insert new values if (incRows > 0) - m1.insert(refMat1.rows()-1, 0) = refMat1(refMat1.rows()-1, 0) = 1; + m1.insert(m1.rows()-1, 0) = refMat1(refMat1.rows()-1, 0) = 1; if (incCols > 0) - m1.insert(0, refMat1.cols()-1) = refMat1(0, refMat1.cols()-1) = 1; + m1.insert(0, m1.cols()-1) = refMat1(0, refMat1.cols()-1) = 1; VERIFY_IS_APPROX(m1, refMat1); @@ -502,7 +503,7 @@ void test_sparse_basic() CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double,ColMajor,long int>(s, s)) )); CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double,RowMajor,long int>(s, s)) )); - CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double,ColMajor,short int>(s, s)) )); - CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double,RowMajor,short int>(s, s)) )); + CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double,ColMajor,short int>(short(s), short(s))) )); + CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double,RowMajor,short int>(short(s), short(s))) )); } } diff --git a/test/sparse_product.cpp b/test/sparse_product.cpp index 51eed428b..664e33887 100644 --- a/test/sparse_product.cpp +++ b/test/sparse_product.cpp @@ -161,6 +161,8 @@ template<typename SparseMatrixType> void sparse_product() VERIFY_IS_APPROX(m3=v2.asDiagonal()*m2, refM3=v2.asDiagonal()*refM2); VERIFY_IS_APPROX(m3=v1.asDiagonal()*m2.transpose(), refM3=v1.asDiagonal()*refM2.transpose()); + VERIFY_IS_APPROX(m3=v2.asDiagonal()*m2*v1.asDiagonal(), refM3=v2.asDiagonal()*refM2*v1.asDiagonal()); + // evaluate to a dense matrix to check the .row() and .col() iterator functions VERIFY_IS_APPROX(d3=m2*d1, refM3=refM2*d1); VERIFY_IS_APPROX(d3=m2.transpose()*d2, refM3=refM2.transpose()*d2); diff --git a/test/sparselu.cpp b/test/sparselu.cpp index 6a9eac065..37980defc 100644 --- a/test/sparselu.cpp +++ b/test/sparselu.cpp @@ -26,7 +26,7 @@ // SparseLU solve does not accept column major matrices for the destination. // However, as expected, the generic check_sparse_square_solving routines produces row-major // rhs and destination matrices when compiled with EIGEN_DEFAULT_TO_ROW_MAJOR -// + #ifdef EIGEN_DEFAULT_TO_ROW_MAJOR #undef EIGEN_DEFAULT_TO_ROW_MAJOR #endif @@ -37,7 +37,7 @@ template<typename T> void test_sparselu_T() { - SparseLU<SparseMatrix<T, ColMajor>, COLAMDOrdering<int> > sparselu_colamd; + SparseLU<SparseMatrix<T, ColMajor> /*, COLAMDOrdering<int>*/ > sparselu_colamd; // COLAMDOrdering is the default SparseLU<SparseMatrix<T, ColMajor>, AMDOrdering<int> > sparselu_amd; SparseLU<SparseMatrix<T, ColMajor, long int>, NaturalOrdering<long int> > sparselu_natural; diff --git a/test/special_numbers.cpp b/test/special_numbers.cpp index b5c83af8b..2f1b704be 100644 --- a/test/special_numbers.cpp +++ b/test/special_numbers.cpp @@ -15,8 +15,8 @@ template<typename Scalar> void special_numbers() int rows = internal::random<int>(1,300); int cols = internal::random<int>(1,300); - Scalar nan = Scalar(0)/Scalar(0); - Scalar inf = Scalar(1)/Scalar(0); + Scalar nan = std::numeric_limits<Scalar>::quiet_NaN(); + Scalar inf = std::numeric_limits<Scalar>::infinity(); Scalar s1 = internal::random<Scalar>(); MatType m1 = MatType::Random(rows,cols), @@ -33,7 +33,7 @@ template<typename Scalar> void special_numbers() mboth = mnan + minf; VERIFY(!m1.hasNaN()); - VERIFY(m1.isFinite()); + VERIFY(m1.allFinite()); VERIFY(mnan.hasNaN()); VERIFY((s1*mnan).hasNaN()); @@ -42,11 +42,11 @@ template<typename Scalar> void special_numbers() VERIFY(mboth.hasNaN()); VERIFY(mboth.array().hasNaN()); - VERIFY(!mnan.isFinite()); - VERIFY(!minf.isFinite()); - VERIFY(!(minf-mboth).isFinite()); - VERIFY(!mboth.isFinite()); - VERIFY(!mboth.array().isFinite()); + VERIFY(!mnan.allFinite()); + VERIFY(!minf.allFinite()); + VERIFY(!(minf-mboth).allFinite()); + VERIFY(!mboth.allFinite()); + VERIFY(!mboth.array().allFinite()); } void test_special_numbers() diff --git a/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h b/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h index ad47d3d84..51dd1d3c4 100644 --- a/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h +++ b/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h @@ -107,7 +107,7 @@ void lmpar2(const QRSolver &qr, const VectorType &diag, const VectorType &qtb, * http://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm */ template<typename _FunctorType> -class LevenbergMarquardt +class LevenbergMarquardt : internal::no_assignment_operator { public: typedef _FunctorType FunctorType; diff --git a/unsupported/Eigen/src/SVD/BDCSVD.h b/unsupported/Eigen/src/SVD/BDCSVD.h index 87ec33a05..11d4882e4 100644 --- a/unsupported/Eigen/src/SVD/BDCSVD.h +++ b/unsupported/Eigen/src/SVD/BDCSVD.h @@ -37,10 +37,14 @@ namespace Eigen { template<typename _MatrixType> class BDCSVD : public SVDBase<_MatrixType> { + typedef SVDBase<_MatrixType> Base; public: + using Base::rows; + using Base::cols; + typedef _MatrixType MatrixType; - typedef typename SVDBase<_MatrixType>::MatrixType::Scalar Scalar; + typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; typedef typename MatrixType::Index Index; enum { diff --git a/unsupported/Eigen/src/SVD/doneInBDCSVD.txt b/unsupported/Eigen/src/SVD/doneInBDCSVD.txt index 6fb6b84a6..8563ddab8 100644 --- a/unsupported/Eigen/src/SVD/doneInBDCSVD.txt +++ b/unsupported/Eigen/src/SVD/doneInBDCSVD.txt @@ -1,7 +1,7 @@ This unsupported package is about a divide and conquer algorithm to compute SVD. The implementation follows as closely as possible the following reference paper : -www.cs.yale.edu/publications/techreports/tr933.pdf +http://www.cs.yale.edu/publications/techreports/tr933.pdf The code documentation uses the same names for variables as the reference paper. The code, deflation included, is working but there are a few things that could be optimised as explained in the TODOBdsvd. diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt index 78b9610d4..a94a3b5e5 100644 --- a/unsupported/test/CMakeLists.txt +++ b/unsupported/test/CMakeLists.txt @@ -90,3 +90,4 @@ ei_add_test(splines) ei_add_test(gmres) ei_add_test(minres) ei_add_test(levenberg_marquardt) +ei_add_test(bdcsvd) |