diff options
author | Gael Guennebaud <g.gael@free.fr> | 2015-02-09 11:41:25 +0100 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2015-02-09 11:41:25 +0100 |
commit | 87629cd6395433966977e868ec091c3fc754956c (patch) | |
tree | 7745c6ffdd622583e077ba4496ba924ff84acb9d | |
parent | bde98df03f98f325a8a288324b6ae88fefbbd797 (diff) |
bug #897: makes iterative sparse solvers use a Ref<SparseMatrix> instead of a SparseMatrix pointer. This fixes usage of iterative solvers with a Map<SparseMatrix>.
-rw-r--r-- | Eigen/src/IterativeLinearSolvers/BiCGSTAB.h | 2 | ||||
-rw-r--r-- | Eigen/src/IterativeLinearSolvers/ConjugateGradient.h | 2 | ||||
-rw-r--r-- | Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h | 53 | ||||
-rw-r--r-- | test/sparse_solver.h | 5 |
4 files changed, 41 insertions, 21 deletions
diff --git a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h index 224fe913f..a50680133 100644 --- a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h +++ b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h @@ -193,7 +193,7 @@ public: m_error = Base::m_tolerance; typename Dest::ColXpr xj(x,j); - if(!internal::bicgstab(*mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_error)) + if(!internal::bicgstab(mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_error)) failed = true; } m_info = failed ? NumericalIssue diff --git a/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h b/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h index b5ef6d60f..3e024bda1 100644 --- a/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h +++ b/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h @@ -206,7 +206,7 @@ public: m_error = Base::m_tolerance; typename Dest::ColXpr xj(x,j); - internal::conjugate_gradient(mp_matrix->template selfadjointView<UpLo>(), b.col(j), xj, + internal::conjugate_gradient(mp_matrix.template selfadjointView<UpLo>(), b.col(j), xj, Base::m_preconditioner, m_iterations, m_error); } diff --git a/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h b/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h index f33c868bb..6f48075b4 100644 --- a/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h +++ b/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h @@ -37,7 +37,7 @@ public: /** Default constructor. */ IterativeSolverBase() - : mp_matrix(0) + : m_dummy(0,0), mp_matrix(m_dummy) { init(); } @@ -52,10 +52,11 @@ public: * this class becomes invalid. Call compute() to update it with the new * matrix A, or modify a copy of A. */ - explicit IterativeSolverBase(const MatrixType& A) + template<typename SparseMatrixDerived> + explicit IterativeSolverBase(const SparseMatrixBase<SparseMatrixDerived>& A) { init(); - compute(A); + compute(A.derived()); } ~IterativeSolverBase() {} @@ -65,9 +66,11 @@ public: * Currently, this function mostly calls analyzePattern on the preconditioner. In the future * we might, for instance, implement column reordering for faster matrix vector products. */ - Derived& analyzePattern(const MatrixType& A) + template<typename SparseMatrixDerived> + Derived& analyzePattern(const SparseMatrixBase<SparseMatrixDerived>& A) { - m_preconditioner.analyzePattern(A); + grab(A); + m_preconditioner.analyzePattern(mp_matrix); m_isInitialized = true; m_analysisIsOk = true; m_info = Success; @@ -83,11 +86,12 @@ public: * this class becomes invalid. Call compute() to update it with the new * matrix A, or modify a copy of A. */ - Derived& factorize(const MatrixType& A) + template<typename SparseMatrixDerived> + Derived& factorize(const SparseMatrixBase<SparseMatrixDerived>& A) { eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); - mp_matrix = &A; - m_preconditioner.factorize(A); + grab(A); + m_preconditioner.factorize(mp_matrix); m_factorizationIsOk = true; m_info = Success; return derived(); @@ -103,10 +107,11 @@ public: * this class becomes invalid. Call compute() to update it with the new * matrix A, or modify a copy of A. */ - Derived& compute(const MatrixType& A) + template<typename SparseMatrixDerived> + Derived& compute(const SparseMatrixBase<SparseMatrixDerived>& A) { - mp_matrix = &A; - m_preconditioner.compute(A); + grab(A); + m_preconditioner.compute(mp_matrix); m_isInitialized = true; m_analysisIsOk = true; m_factorizationIsOk = true; @@ -115,9 +120,9 @@ public: } /** \internal */ - Index rows() const { return mp_matrix ? mp_matrix->rows() : 0; } + Index rows() const { return mp_matrix.rows(); } /** \internal */ - Index cols() const { return mp_matrix ? mp_matrix->cols() : 0; } + Index cols() const { return mp_matrix.cols(); } /** \returns the tolerance threshold used by the stopping criteria */ RealScalar tolerance() const { return m_tolerance; } @@ -135,13 +140,18 @@ public: /** \returns a read-only reference to the preconditioner. */ const Preconditioner& preconditioner() const { return m_preconditioner; } - /** \returns the max number of iterations */ + /** \returns the max number of iterations. + * It is either the value setted by setMaxIterations or, by default, + * twice the number of columns of the matrix. + */ int maxIterations() const { - return (mp_matrix && m_maxIterations<0) ? mp_matrix->cols() : m_maxIterations; + return (m_maxIterations<0) ? 2*mp_matrix.cols() : m_maxIterations; } - /** Sets the max number of iterations */ + /** Sets the max number of iterations. + * Default is twice the number of columns of the matrix. + */ Derived& setMaxIterations(int maxIters) { m_maxIterations = maxIters; @@ -210,7 +220,16 @@ protected: m_maxIterations = -1; m_tolerance = NumTraits<Scalar>::epsilon(); } - const MatrixType* mp_matrix; + + template<typename SparseMatrixDerived> + void grab(const SparseMatrixBase<SparseMatrixDerived> &A) + { + mp_matrix.~Ref<const MatrixType>(); + ::new (&mp_matrix) Ref<const MatrixType>(A); + } + + MatrixType m_dummy; + Ref<const MatrixType> mp_matrix; Preconditioner m_preconditioner; int m_maxIterations; diff --git a/test/sparse_solver.h b/test/sparse_solver.h index ee350d561..887ff02a9 100644 --- a/test/sparse_solver.h +++ b/test/sparse_solver.h @@ -32,7 +32,7 @@ void check_sparse_solving(Solver& solver, const typename Solver::MatrixType& A, x = solver.solve(b); if (solver.info() != Success) { - std::cerr << "sparse solver testing: solving failed\n"; + std::cerr << "sparse solver testing: solving failed (" << typeid(Solver).name() << ")\n"; return; } VERIFY(oldb.isApprox(b) && "sparse solver testing: the rhs should not be modified!"); @@ -75,7 +75,8 @@ void check_sparse_solving(Solver& solver, const typename Solver::MatrixType& A, xm = solver.solve(bm); if (solver.info() != Success) { - std::cerr << "sparse solver testing: solving failed\n"; + std::cerr << "sparse solver testing: solving with a Map failed\n"; + exit(0); return; } VERIFY(oldb.isApprox(bm) && "sparse solver testing: the rhs should not be modified!"); |