Merged in benoitsteiner/opencl (pull request PR-246) Improved support for OpenCL
diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index ddf4875..87ca8d4 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h
@@ -351,7 +351,7 @@ Index ret; if((ret=unblocked(A11))>=0) return k+ret; if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21); - if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,-1); // bottleneck + if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,typename NumTraits<RealScalar>::Literal(-1)); // bottleneck } return -1; }
diff --git a/Eigen/src/CholmodSupport/CholmodSupport.h b/Eigen/src/CholmodSupport/CholmodSupport.h index 4a4aa21..5719720 100644 --- a/Eigen/src/CholmodSupport/CholmodSupport.h +++ b/Eigen/src/CholmodSupport/CholmodSupport.h
@@ -55,7 +55,7 @@ * Note that the data are shared. */ template<typename _Scalar, int _Options, typename _StorageIndex> -cholmod_sparse viewAsCholmod(SparseMatrix<_Scalar,_Options,_StorageIndex>& mat) +cholmod_sparse viewAsCholmod(Ref<SparseMatrix<_Scalar,_Options,_StorageIndex> > mat) { cholmod_sparse res; res.nzmax = mat.nonZeros(); @@ -104,7 +104,14 @@ template<typename _Scalar, int _Options, typename _Index> const cholmod_sparse viewAsCholmod(const SparseMatrix<_Scalar,_Options,_Index>& mat) { - cholmod_sparse res = viewAsCholmod(mat.const_cast_derived()); + cholmod_sparse res = viewAsCholmod(Ref<SparseMatrix<_Scalar,_Options,_Index> >(mat.const_cast_derived())); + return res; +} + +template<typename _Scalar, int _Options, typename _Index> +const cholmod_sparse viewAsCholmod(const SparseVector<_Scalar,_Options,_Index>& mat) +{ + cholmod_sparse res = viewAsCholmod(Ref<SparseMatrix<_Scalar,_Options,_Index> >(mat.const_cast_derived())); return res; } @@ -113,7 +120,7 @@ template<typename _Scalar, int _Options, typename _Index, unsigned int UpLo> cholmod_sparse viewAsCholmod(const SparseSelfAdjointView<const SparseMatrix<_Scalar,_Options,_Index>, UpLo>& mat) { - cholmod_sparse res = viewAsCholmod(mat.matrix().const_cast_derived()); + cholmod_sparse res = viewAsCholmod(Ref<SparseMatrix<_Scalar,_Options,_Index> >(mat.matrix().const_cast_derived())); if(UpLo==Upper) res.stype = 1; if(UpLo==Lower) res.stype = -1; @@ -298,8 +305,8 @@ } /** \internal */ - template<typename RhsScalar, int RhsOptions, typename RhsIndex, typename DestScalar, int DestOptions, typename DestIndex> - void _solve_impl(const SparseMatrix<RhsScalar,RhsOptions,RhsIndex> &b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const + template<typename RhsDerived, typename DestDerived> + void _solve_impl(const SparseMatrixBase<RhsDerived> &b, SparseMatrixBase<DestDerived> &dest) const { eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); const Index size = m_cholmodFactor->n; @@ -307,7 +314,8 @@ eigen_assert(size==b.rows()); // note: cs stands for Cholmod Sparse - cholmod_sparse b_cs = viewAsCholmod(b); + Ref<SparseMatrix<typename RhsDerived::Scalar,ColMajor,typename RhsDerived::StorageIndex> > b_ref(b.const_cast_derived()); + cholmod_sparse b_cs = viewAsCholmod(b_ref); cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod); if(!x_cs) { @@ -315,7 +323,7 @@ return; } // TODO optimize this copy by swapping when possible (be careful with alignment, etc.) - dest = viewAsEigen<DestScalar,DestOptions,DestIndex>(*x_cs); + dest.derived() = viewAsEigen<typename DestDerived::Scalar,ColMajor,typename DestDerived::StorageIndex>(*x_cs); cholmod_free_sparse(&x_cs, &m_cholmod); } #endif // EIGEN_PARSED_BY_DOXYGEN @@ -570,7 +578,7 @@ * This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed. * * \warning Only double precision real and complex scalar types are supported by Cholmod. - * + * * \sa \ref TutorialSparseSolverConcept */ template<typename _MatrixType, int _UpLo = Lower>
diff --git a/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h b/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h index 3d62fef..7c2326e 100644 --- a/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h +++ b/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h
@@ -330,25 +330,27 @@ } /** \internal */ - template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex> - void _solve_impl(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const + template<typename Rhs, typename DestDerived> + void _solve_impl(const Rhs& b, SparseMatrixBase<DestDerived> &aDest) const { eigen_assert(rows()==b.rows()); Index rhsCols = b.cols(); Index size = b.rows(); + DestDerived& dest(aDest.derived()); + typedef typename DestDerived::Scalar DestScalar; Eigen::Matrix<DestScalar,Dynamic,1> tb(size); Eigen::Matrix<DestScalar,Dynamic,1> tx(cols()); // We do not directly fill dest because sparse expressions have to be free of aliasing issue. // For non square least-square problems, b and dest might not have the same size whereas they might alias each-other. - SparseMatrix<DestScalar,DestOptions,DestIndex> tmp(cols(),rhsCols); + typename DestDerived::PlainObject tmp(cols(),rhsCols); for(Index k=0; k<rhsCols; ++k) { tb = b.col(k); tx = derived().solve(tb); tmp.col(k) = tx.sparseView(0); } - tmp.swap(dest); + dest.swap(tmp); } protected:
diff --git a/Eigen/src/SparseCore/SparseMatrix.h b/Eigen/src/SparseCore/SparseMatrix.h index 64ca5fc..ec650b2 100644 --- a/Eigen/src/SparseCore/SparseMatrix.h +++ b/Eigen/src/SparseCore/SparseMatrix.h
@@ -94,6 +94,7 @@ { typedef SparseCompressedBase<SparseMatrix> Base; using Base::convert_index; + friend class SparseVector<_Scalar,0,_Index>; public: using Base::isCompressed; using Base::nonZeros;
diff --git a/Eigen/src/SparseCore/SparseSolverBase.h b/Eigen/src/SparseCore/SparseSolverBase.h index 1cb7080..b4c9a42 100644 --- a/Eigen/src/SparseCore/SparseSolverBase.h +++ b/Eigen/src/SparseCore/SparseSolverBase.h
@@ -19,7 +19,8 @@ * The rhs is decomposed into small vertical panels which are solved through dense temporaries. */ template<typename Decomposition, typename Rhs, typename Dest> -void solve_sparse_through_dense_panels(const Decomposition &dec, const Rhs& rhs, Dest &dest) +typename enable_if<Rhs::ColsAtCompileTime!=1 && Dest::ColsAtCompileTime!=1>::type +solve_sparse_through_dense_panels(const Decomposition &dec, const Rhs& rhs, Dest &dest) { EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); typedef typename Dest::Scalar DestScalar; @@ -40,6 +41,19 @@ } } +// Overload for vector as rhs +template<typename Decomposition, typename Rhs, typename Dest> +typename enable_if<Rhs::ColsAtCompileTime==1 || Dest::ColsAtCompileTime==1>::type +solve_sparse_through_dense_panels(const Decomposition &dec, const Rhs& rhs, Dest &dest) +{ + typedef typename Dest::Scalar DestScalar; + Index size = rhs.rows(); + Eigen::Matrix<DestScalar,Dynamic,1> rhs_dense(rhs); + Eigen::Matrix<DestScalar,Dynamic,1> dest_dense(size); + dest_dense = dec.solve(rhs_dense); + dest = dest_dense.sparseView(); +} + } // end namespace internal /** \class SparseSolverBase
diff --git a/Eigen/src/SparseCore/SparseVector.h b/Eigen/src/SparseCore/SparseVector.h index 00ee6ec..19b0fbc 100644 --- a/Eigen/src/SparseCore/SparseVector.h +++ b/Eigen/src/SparseCore/SparseVector.h
@@ -290,6 +290,14 @@ m_data.swap(other.m_data); } + template<int OtherOptions> + inline void swap(SparseMatrix<Scalar,OtherOptions,StorageIndex>& other) + { + eigen_assert(other.outerSize()==1); + std::swap(m_size, other.m_innerSize); + m_data.swap(other.m_data); + } + inline SparseVector& operator=(const SparseVector& other) { if (other.isRValue()) @@ -403,6 +411,7 @@ : evaluator_base<SparseVector<_Scalar,_Options,_Index> > { typedef SparseVector<_Scalar,_Options,_Index> SparseVectorType; + typedef evaluator_base<SparseVectorType> Base; typedef typename SparseVectorType::InnerIterator InnerIterator; typedef typename SparseVectorType::ReverseInnerIterator ReverseInnerIterator; @@ -410,20 +419,22 @@ CoeffReadCost = NumTraits<_Scalar>::ReadCost, Flags = SparseVectorType::Flags }; + + evaluator() : Base() {} - explicit evaluator(const SparseVectorType &mat) : m_matrix(mat) + explicit evaluator(const SparseVectorType &mat) : m_matrix(&mat) { EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); } inline Index nonZerosEstimate() const { - return m_matrix.nonZeros(); + return m_matrix->nonZeros(); } - operator SparseVectorType&() { return m_matrix.const_cast_derived(); } - operator const SparseVectorType&() const { return m_matrix; } + operator SparseVectorType&() { return m_matrix->const_cast_derived(); } + operator const SparseVectorType&() const { return *m_matrix; } - const SparseVectorType &m_matrix; + const SparseVectorType *m_matrix; }; template< typename Dest, typename Src>
diff --git a/test/sparse_solver.h b/test/sparse_solver.h index fd6199f..5145bc3 100644 --- a/test/sparse_solver.h +++ b/test/sparse_solver.h
@@ -272,6 +272,7 @@ typedef typename Mat::Scalar Scalar; typedef typename Mat::StorageIndex StorageIndex; typedef SparseMatrix<Scalar,ColMajor, StorageIndex> SpMat; + typedef SparseVector<Scalar, 0, StorageIndex> SpVec; typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; typedef Matrix<Scalar,Dynamic,1> DenseVector; @@ -288,6 +289,8 @@ DenseVector b = DenseVector::Random(size); DenseMatrix dB(size,rhsCols); initSparse<Scalar>(density, dB, B, ForceNonZeroDiag); + SpVec c = B.col(0); + DenseVector dc = dB.col(0); CALL_SUBTEST( check_sparse_solving(solver, A, b, dA, b) ); CALL_SUBTEST( check_sparse_solving(solver, halfA, b, dA, b) ); @@ -295,6 +298,8 @@ CALL_SUBTEST( check_sparse_solving(solver, halfA, dB, dA, dB) ); CALL_SUBTEST( check_sparse_solving(solver, A, B, dA, dB) ); CALL_SUBTEST( check_sparse_solving(solver, halfA, B, dA, dB) ); + CALL_SUBTEST( check_sparse_solving(solver, A, c, dA, dc) ); + CALL_SUBTEST( check_sparse_solving(solver, halfA, c, dA, dc) ); // check only once if(i==0) @@ -396,6 +401,7 @@ typedef typename Solver::MatrixType Mat; typedef typename Mat::Scalar Scalar; typedef SparseMatrix<Scalar,ColMajor, typename Mat::StorageIndex> SpMat; + typedef SparseVector<Scalar, 0, typename Mat::StorageIndex> SpVec; typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; typedef Matrix<Scalar,Dynamic,1> DenseVector; @@ -413,9 +419,12 @@ double density = (std::max)(8./(size*rhsCols), 0.1); initSparse<Scalar>(density, dB, B, ForceNonZeroDiag); B.makeCompressed(); + SpVec c = B.col(0); + DenseVector dc = dB.col(0); CALL_SUBTEST(check_sparse_solving(solver, A, b, dA, b)); CALL_SUBTEST(check_sparse_solving(solver, A, dB, dA, dB)); CALL_SUBTEST(check_sparse_solving(solver, A, B, dA, dB)); + CALL_SUBTEST(check_sparse_solving(solver, A, c, dA, dc)); // check only once if(i==0)
diff --git a/unsupported/test/NonLinearOptimization.cpp b/unsupported/test/NonLinearOptimization.cpp index 6a5ed05..1d682dd 100644 --- a/unsupported/test/NonLinearOptimization.cpp +++ b/unsupported/test/NonLinearOptimization.cpp
@@ -12,8 +12,6 @@ // It is intended to be done for this test only. #include <Eigen/src/Core/util/DisableStupidWarnings.h> -using std::sqrt; - // tolerance for chekcing number of iterations #define LM_EVAL_COUNT_TOL 4/3
diff --git a/unsupported/test/mpreal_support.cpp b/unsupported/test/mpreal_support.cpp index ffa5691..685e7ea 100644 --- a/unsupported/test/mpreal_support.cpp +++ b/unsupported/test/mpreal_support.cpp
@@ -12,6 +12,7 @@ // set precision to 256 bits (double has only 53 bits) mpreal::set_default_prec(256); typedef Matrix<mpreal,Eigen::Dynamic,Eigen::Dynamic> MatrixXmp; + typedef Matrix<std::complex<mpreal>,Eigen::Dynamic,Eigen::Dynamic> MatrixXcmp; std::cerr << "epsilon = " << NumTraits<mpreal>::epsilon() << "\n"; std::cerr << "dummy_precision = " << NumTraits<mpreal>::dummy_precision() << "\n"; @@ -25,6 +26,10 @@ MatrixXmp B = MatrixXmp::Random(s,s); MatrixXmp S = A.adjoint() * A; MatrixXmp X; + MatrixXcmp Ac = MatrixXcmp::Random(s,s); + MatrixXcmp Bc = MatrixXcmp::Random(s,s); + MatrixXcmp Sc = Ac.adjoint() * Ac; + MatrixXcmp Xc; // Basic stuffs VERIFY_IS_APPROX(A.real(), A); @@ -37,6 +42,9 @@ // Cholesky X = S.selfadjointView<Lower>().llt().solve(B); VERIFY_IS_APPROX((S.selfadjointView<Lower>()*X).eval(),B); + + Xc = Sc.selfadjointView<Lower>().llt().solve(Bc); + VERIFY_IS_APPROX((Sc.selfadjointView<Lower>()*Xc).eval(),Bc); // partial LU X = A.lu().solve(B);