From d8df318d77b8a9bd9d6274f25145639603c2e8d4 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 2 Sep 2008 19:55:26 +0000 Subject: [PATCH] resurrected sparse triangular solver --- Eigen/src/Core/Flagged.h | 1 + Eigen/src/Core/SolveTriangular.h | 6 +- Eigen/src/Core/util/Constants.h | 10 +--- Eigen/src/Core/util/XprHelper.h | 2 +- Eigen/src/Geometry/Quaternion.h | 2 +- Eigen/src/Sparse/SparseMatrix.h | 6 ++ Eigen/src/Sparse/SparseMatrixBase.h | 3 - Eigen/src/Sparse/SparseUtil.h | 2 +- Eigen/src/Sparse/TriangularSolver.h | 88 ++++++++++++----------------- test/sparse.cpp | 83 +++++++++++++++++++++------ 10 files changed, 117 insertions(+), 86 deletions(-) diff --git a/Eigen/src/Core/Flagged.h b/Eigen/src/Core/Flagged.h index a6abe617b..387854b20 100644 --- a/Eigen/src/Core/Flagged.h +++ b/Eigen/src/Core/Flagged.h @@ -62,6 +62,7 @@ template clas EIGEN_GENERIC_PUBLIC_INTERFACE(Flagged) typedef typename ei_meta_if::ret, ExpressionType, const ExpressionType&>::ret ExpressionTypeNested; + typedef typename ExpressionType::InnerIterator InnerIterator; inline Flagged(const ExpressionType& matrix) : m_matrix(matrix) {} diff --git a/Eigen/src/Core/SolveTriangular.h b/Eigen/src/Core/SolveTriangular.h index 0f7aa3b3e..aaebc5989 100755 --- a/Eigen/src/Core/SolveTriangular.h +++ b/Eigen/src/Core/SolveTriangular.h @@ -35,7 +35,7 @@ template::value ? -1 // this is to solve ambiguous specializations - : int(Lhs::Flags) & RowMajorBit ? RowMajor : ColMajor + : int(Lhs::Flags) & (RowMajorBit|SparseBit) > struct ei_solve_triangular_selector; @@ -51,7 +51,7 @@ struct ei_solve_triangular_selector,Rhs,UpLo,StorageOrder> // forward substitution, row-major template -struct ei_solve_triangular_selector +struct ei_solve_triangular_selector { typedef typename Rhs::Scalar Scalar; static void run(const Lhs& lhs, Rhs& other) @@ -138,7 +138,7 @@ struct ei_solve_triangular_selector // - inv(Upper, ColMajor) * Column vector // - inv(Upper,UnitDiag,ColMajor) * Column vector template -struct ei_solve_triangular_selector +struct ei_solve_triangular_selector { typedef typename Rhs::Scalar Scalar; typedef typename ei_packet_traits::type Packet; diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index d08af60e0..e852fffa0 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -199,20 +199,16 @@ enum { NoUnrolling }; -enum { - Dense = 0, - Sparse = SparseBit -}; - enum { ColMajor = 0, RowMajor = RowMajorBit }; enum { - NoDirectAccess = 0, + IsDense = 0, + NoDirectAccess = 0, HasDirectAccess = DirectAccessBit, - IsSparse = SparseBit + IsSparse = SparseBit }; const int FullyCoherentAccessPattern = 0x1; diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h index 225100118..00f1a39ea 100644 --- a/Eigen/src/Core/util/XprHelper.h +++ b/Eigen/src/Core/util/XprHelper.h @@ -110,7 +110,7 @@ template struct ei_size_at_compile_time template::Flags&SparseBit> class ei_eval; -template struct ei_eval +template struct ei_eval { typedef Matrix::Scalar, ei_traits::RowsAtCompileTime, diff --git a/Eigen/src/Geometry/Quaternion.h b/Eigen/src/Geometry/Quaternion.h index a2cc785c7..970280d33 100644 --- a/Eigen/src/Geometry/Quaternion.h +++ b/Eigen/src/Geometry/Quaternion.h @@ -375,7 +375,7 @@ Quaternion Quaternion::slerp(Scalar t, const Quaternion& other) Scalar theta = std::acos(ei_abs(d)); Scalar sinTheta = ei_sin(theta); - Scalar scale0 = ei_sin( ( 1 - t ) * theta) / sinTheta; + Scalar scale0 = ei_sin( ( Scalar(1) - t ) * theta) / sinTheta; Scalar scale1 = ei_sin( ( t * theta) ) / sinTheta; if (d<0) scale1 = -scale1; diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h index d62718491..b4c4fe563 100644 --- a/Eigen/src/Sparse/SparseMatrix.h +++ b/Eigen/src/Sparse/SparseMatrix.h @@ -261,6 +261,12 @@ class SparseMatrix::InnerIterator : m_matrix(mat), m_id(mat.m_outerIndex[outer]), m_start(m_id), m_end(mat.m_outerIndex[outer+1]) {} + template + InnerIterator(const Flagged& mat, int outer) + : m_matrix(mat._expression()), m_id(m_matrix.m_outerIndex[outer]), + m_start(m_id), m_end(m_matrix.m_outerIndex[outer+1]) + {} + InnerIterator& operator++() { m_id++; return *this; } Scalar value() const { return m_matrix.m_data.value(m_id); } diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h index 83dcefe06..1dcf83c24 100644 --- a/Eigen/src/Sparse/SparseMatrixBase.h +++ b/Eigen/src/Sparse/SparseMatrixBase.h @@ -160,9 +160,6 @@ class SparseMatrixBase : public MatrixBase return s; } - template - OtherDerived solveTriangular(const MatrixBase& other) const; - protected: bool m_isRValue; diff --git a/Eigen/src/Sparse/SparseUtil.h b/Eigen/src/Sparse/SparseUtil.h index bcf87ad90..1a860b31d 100644 --- a/Eigen/src/Sparse/SparseUtil.h +++ b/Eigen/src/Sparse/SparseUtil.h @@ -48,7 +48,7 @@ template struct ei_support_access_patter }; }; -template class ei_eval +template class ei_eval { typedef typename ei_traits::Scalar _Scalar; enum { diff --git a/Eigen/src/Sparse/TriangularSolver.h b/Eigen/src/Sparse/TriangularSolver.h index 45679fe86..5fb8396f2 100644 --- a/Eigen/src/Sparse/TriangularSolver.h +++ b/Eigen/src/Sparse/TriangularSolver.h @@ -25,42 +25,42 @@ #ifndef EIGEN_SPARSETRIANGULARSOLVER_H #define EIGEN_SPARSETRIANGULARSOLVER_H -template -struct ei_sparse_trisolve_selector; +// template +// struct ei_sparse_trisolve_selector; // forward substitution, row-major template -struct ei_sparse_trisolve_selector +struct ei_solve_triangular_selector { typedef typename Rhs::Scalar Scalar; - static void run(const Lhs& lhs, const Rhs& rhs, Rhs& res) + static void run(const Lhs& lhs, Rhs& other) { - for(int col=0 ; col // backward substitution, row-major template -struct ei_sparse_trisolve_selector +struct ei_solve_triangular_selector { typedef typename Rhs::Scalar Scalar; - static void run(const Lhs& lhs, const Rhs& rhs, Rhs& res) + static void run(const Lhs& lhs, Rhs& other) { - for(int col=0 ; col=0 ; --i) { - Scalar tmp = rhs.coeff(i,col); + Scalar tmp = other.coeff(i,col); typename Lhs::InnerIterator it(lhs, i); for(++it; it; ++it) { - tmp -= it.value() * res.coeff(it.index(),col); + tmp -= it.value() * other.coeff(it.index(),col); } if (Lhs::Flags & UnitDiagBit) - res.coeffRef(i,col) = tmp; + other.coeffRef(i,col) = tmp; else { typename Lhs::InnerIterator it(lhs, i); ei_assert(it.index() == i); - res.coeffRef(i,col) = tmp/it.value(); + other.coeffRef(i,col) = tmp/it.value(); } } } @@ -100,26 +100,25 @@ struct ei_sparse_trisolve_selector // forward substitution, col-major template -struct ei_sparse_trisolve_selector +struct ei_solve_triangular_selector { typedef typename Rhs::Scalar Scalar; - static void run(const Lhs& lhs, const Rhs& rhs, Rhs& res) + static void run(const Lhs& lhs, Rhs& other) { - // NOTE we could avoid this copy using an in-place API - res = rhs; - for(int col=0 ; col // backward substitution, col-major template -struct ei_sparse_trisolve_selector +struct ei_solve_triangular_selector { typedef typename Rhs::Scalar Scalar; - static void run(const Lhs& lhs, const Rhs& rhs, Rhs& res) + static void run(const Lhs& lhs, Rhs& other) { - // NOTE we could avoid this copy using an in-place API - res = rhs; - for(int col=0 ; col=0; --i) { @@ -142,28 +139,15 @@ struct ei_sparse_trisolve_selector { // FIXME lhs.coeff(i,i) might not be always efficient while it must simply be the // last element of the column ! - res.coeffRef(i,col) /= lhs.coeff(i,i); + other.coeffRef(i,col) /= lhs.coeff(i,i); } - Scalar tmp = res.coeffRef(i,col); + Scalar tmp = other.coeffRef(i,col); typename Lhs::InnerIterator it(lhs, i); for(; it && it.index() -template -OtherDerived SparseMatrixBase::solveTriangular(const MatrixBase& other) const -{ - ei_assert(derived().cols() == other.rows()); - ei_assert(!(Flags & ZeroDiagBit)); - ei_assert(Flags & (UpperTriangularBit|LowerTriangularBit)); - - OtherDerived res(other.rows(), other.cols()); - ei_sparse_trisolve_selector::run(derived(), other.derived(), res); - return res; -} - #endif // EIGEN_SPARSETRIANGULARSOLVER_H diff --git a/test/sparse.cpp b/test/sparse.cpp index 7decb4111..3095297ce 100644 --- a/test/sparse.cpp +++ b/test/sparse.cpp @@ -25,36 +25,63 @@ #include "main.h" #include -template void sparse(int rows, int cols) +enum { + ForceNonZeroDiag = 1, + MakeLowerTriangular = 2, + MakeUpperTriangular = 4 +}; + +template void +initSparse(double density, + Matrix& refMat, + SparseMatrix& sparseMat, + int flags = 0, + std::vector* zeroCoords = 0, + std::vector* nonzeroCoords = 0) { - double density = std::max(8./(rows*cols), 0.01); - typedef Matrix DenseMatrix; - Scalar eps = 1e-6; - - SparseMatrix m(rows, cols); - DenseMatrix refMat = DenseMatrix::Zero(rows, cols); - - std::vector zeroCoords; - std::vector nonzeroCoords; - m.startFill(rows*cols*density); - for(int j=0; j(0,1) < density) ? ei_random() : 0; + if ((flags&ForceNonZeroDiag) && (i==j)) + while (ei_abs(v)<1e-2) + v = ei_random(); + if ((flags & MakeLowerTriangular) && j>i) + v = 0; + else if ((flags & MakeUpperTriangular) && jpush_back(Vector2i(i,j)); } - else + else if (zeroCoords) { - zeroCoords.push_back(Vector2i(i,j)); + zeroCoords->push_back(Vector2i(i,j)); } refMat(i,j) = v; } } - m.endFill(); + sparseMat.endFill(); +} + +template void sparse(int rows, int cols) +{ + double density = std::max(8./(rows*cols), 0.01); + typedef Matrix DenseMatrix; + typedef Matrix DenseVector; + Scalar eps = 1e-6; + + SparseMatrix m(rows, cols); + DenseMatrix refMat = DenseMatrix::Zero(rows, cols); + DenseVector vec1 = DenseVector::Random(rows); + + std::vector zeroCoords; + std::vector nonzeroCoords; + initSparse(density, refMat, m, 0, &zeroCoords, &nonzeroCoords); VERIFY(zeroCoords.size()>0 && "re-run the test"); VERIFY(nonzeroCoords.size()>0 && "re-run the test"); @@ -145,10 +172,30 @@ template void sparse(int rows, int cols) } VERIFY_IS_APPROX(m, refMat); + // test triangular solver + { + DenseVector vec2 = vec1, vec3 = vec1; + SparseMatrix m2(rows, cols); + DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols); + + // lower + initSparse(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); + VERIFY_IS_APPROX(refMat2.template marked().solveTriangular(vec2), + m2.template marked().solveTriangular(vec3)); + + // upper + initSparse(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords); + VERIFY_IS_APPROX(refMat2.template marked().solveTriangular(vec2), + m2.template marked().solveTriangular(vec3)); + + // TODO test row major + } + } void test_sparse() { sparse(8, 8); sparse(16, 16); + sparse(33, 33); }