diff options
Diffstat (limited to 'third_party/eigen3/Eigen/src/SparseCore/TriangularSolver.h')
| -rw-r--r-- | third_party/eigen3/Eigen/src/SparseCore/TriangularSolver.h | 334 |
1 files changed, 334 insertions, 0 deletions
diff --git a/third_party/eigen3/Eigen/src/SparseCore/TriangularSolver.h b/third_party/eigen3/Eigen/src/SparseCore/TriangularSolver.h new file mode 100644 index 0000000000..cb8ad82b4f --- /dev/null +++ b/third_party/eigen3/Eigen/src/SparseCore/TriangularSolver.h @@ -0,0 +1,334 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_SPARSETRIANGULARSOLVER_H +#define EIGEN_SPARSETRIANGULARSOLVER_H + +namespace Eigen { + +namespace internal { + +template<typename Lhs, typename Rhs, int Mode, + int UpLo = (Mode & Lower) + ? Lower + : (Mode & Upper) + ? Upper + : -1, + int StorageOrder = int(traits<Lhs>::Flags) & RowMajorBit> +struct sparse_solve_triangular_selector; + +// forward substitution, row-major +template<typename Lhs, typename Rhs, int Mode> +struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Lower,RowMajor> +{ + typedef typename Rhs::Scalar Scalar; + static void run(const Lhs& lhs, Rhs& other) + { + for(int col=0 ; col<other.cols() ; ++col) + { + for(int i=0; i<lhs.rows(); ++i) + { + Scalar tmp = other.coeff(i,col); + Scalar lastVal(0); + int lastIndex = 0; + for(typename Lhs::InnerIterator it(lhs, i); it; ++it) + { + lastVal = it.value(); + lastIndex = it.index(); + if(lastIndex==i) + break; + tmp -= lastVal * other.coeff(lastIndex,col); + } + if (Mode & UnitDiag) + other.coeffRef(i,col) = tmp; + else + { + eigen_assert(lastIndex==i); + other.coeffRef(i,col) = tmp/lastVal; + } + } + } + } +}; + +// backward substitution, row-major +template<typename Lhs, typename Rhs, int Mode> +struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,RowMajor> +{ + typedef typename Rhs::Scalar Scalar; + static void run(const Lhs& lhs, Rhs& other) + { + for(int col=0 ; col<other.cols() ; ++col) + { + for(int i=lhs.rows()-1 ; i>=0 ; --i) + { + Scalar tmp = other.coeff(i,col); + Scalar l_ii = 0; + typename Lhs::InnerIterator it(lhs, i); + while(it && it.index()<i) + ++it; + if(!(Mode & UnitDiag)) + { + eigen_assert(it && it.index()==i); + l_ii = it.value(); + ++it; + } + else if (it && it.index() == i) + ++it; + for(; it; ++it) + { + tmp -= it.value() * other.coeff(it.index(),col); + } + + if (Mode & UnitDiag) + other.coeffRef(i,col) = tmp; + else + other.coeffRef(i,col) = tmp/l_ii; + } + } + } +}; + +// forward substitution, col-major +template<typename Lhs, typename Rhs, int Mode> +struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Lower,ColMajor> +{ + typedef typename Rhs::Scalar Scalar; + static void run(const Lhs& lhs, Rhs& other) + { + for(int col=0 ; col<other.cols() ; ++col) + { + for(int i=0; i<lhs.cols(); ++i) + { + Scalar& tmp = other.coeffRef(i,col); + if (tmp!=Scalar(0)) // optimization when other is actually sparse + { + typename Lhs::InnerIterator it(lhs, i); + while(it && it.index()<i) + ++it; + if(!(Mode & UnitDiag)) + { + eigen_assert(it && it.index()==i); + tmp /= it.value(); + } + if (it && it.index()==i) + ++it; + for(; it; ++it) + other.coeffRef(it.index(), col) -= tmp * it.value(); + } + } + } + } +}; + +// backward substitution, col-major +template<typename Lhs, typename Rhs, int Mode> +struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,ColMajor> +{ + typedef typename Rhs::Scalar Scalar; + static void run(const Lhs& lhs, Rhs& other) + { + for(int col=0 ; col<other.cols() ; ++col) + { + for(int i=lhs.cols()-1; i>=0; --i) + { + Scalar& tmp = other.coeffRef(i,col); + if (tmp!=Scalar(0)) // optimization when other is actually sparse + { + if(!(Mode & UnitDiag)) + { + // TODO replace this by a binary search. make sure the binary search is safe for partially sorted elements + typename Lhs::ReverseInnerIterator it(lhs, i); + while(it && it.index()!=i) + --it; + eigen_assert(it && it.index()==i); + other.coeffRef(i,col) /= it.value(); + } + typename Lhs::InnerIterator it(lhs, i); + for(; it && it.index()<i; ++it) + other.coeffRef(it.index(), col) -= tmp * it.value(); + } + } + } + } +}; + +} // end namespace internal + +template<typename ExpressionType,int Mode> +template<typename OtherDerived> +void SparseTriangularView<ExpressionType,Mode>::solveInPlace(MatrixBase<OtherDerived>& other) const +{ + eigen_assert(m_matrix.cols() == m_matrix.rows() && m_matrix.cols() == other.rows()); + eigen_assert((!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower))); + + enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit }; + + typedef typename internal::conditional<copy, + typename internal::plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::type OtherCopy; + OtherCopy otherCopy(other.derived()); + + internal::sparse_solve_triangular_selector<ExpressionType, typename internal::remove_reference<OtherCopy>::type, Mode>::run(m_matrix, otherCopy); + + if (copy) + other = otherCopy; +} + +template<typename ExpressionType,int Mode> +template<typename OtherDerived> +typename internal::plain_matrix_type_column_major<OtherDerived>::type +SparseTriangularView<ExpressionType,Mode>::solve(const MatrixBase<OtherDerived>& other) const +{ + typename internal::plain_matrix_type_column_major<OtherDerived>::type res(other); + solveInPlace(res); + return res; +} + +// pure sparse path + +namespace internal { + +template<typename Lhs, typename Rhs, int Mode, + int UpLo = (Mode & Lower) + ? Lower + : (Mode & Upper) + ? Upper + : -1, + int StorageOrder = int(Lhs::Flags) & (RowMajorBit)> +struct sparse_solve_triangular_sparse_selector; + +// forward substitution, col-major +template<typename Lhs, typename Rhs, int Mode, int UpLo> +struct sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor> +{ + typedef typename Rhs::Scalar Scalar; + typedef typename promote_index_type<typename traits<Lhs>::Index, + typename traits<Rhs>::Index>::type Index; + static void run(const Lhs& lhs, Rhs& other) + { + const bool IsLower = (UpLo==Lower); + AmbiVector<Scalar,Index> tempVector(other.rows()*2); + tempVector.setBounds(0,other.rows()); + + Rhs res(other.rows(), other.cols()); + res.reserve(other.nonZeros()); + + for(int col=0 ; col<other.cols() ; ++col) + { + // FIXME estimate number of non zeros + tempVector.init(.99/*float(other.col(col).nonZeros())/float(other.rows())*/); + tempVector.setZero(); + tempVector.restart(); + for (typename Rhs::InnerIterator rhsIt(other, col); rhsIt; ++rhsIt) + { + tempVector.coeffRef(rhsIt.index()) = rhsIt.value(); + } + + for(int i=IsLower?0:lhs.cols()-1; + IsLower?i<lhs.cols():i>=0; + i+=IsLower?1:-1) + { + tempVector.restart(); + Scalar& ci = tempVector.coeffRef(i); + if (ci!=Scalar(0)) + { + // find + typename Lhs::InnerIterator it(lhs, i); + if(!(Mode & UnitDiag)) + { + if (IsLower) + { + eigen_assert(it.index()==i); + ci /= it.value(); + } + else + ci /= lhs.coeff(i,i); + } + tempVector.restart(); + if (IsLower) + { + if (it.index()==i) + ++it; + for(; it; ++it) + tempVector.coeffRef(it.index()) -= ci * it.value(); + } + else + { + for(; it && it.index()<i; ++it) + tempVector.coeffRef(it.index()) -= ci * it.value(); + } + } + } + + + int count = 0; + // FIXME compute a reference value to filter zeros + for (typename AmbiVector<Scalar,Index>::Iterator it(tempVector/*,1e-12*/); it; ++it) + { + ++ count; +// std::cerr << "fill " << it.index() << ", " << col << "\n"; +// std::cout << it.value() << " "; + // FIXME use insertBack + res.insert(it.index(), col) = it.value(); + } +// std::cout << "tempVector.nonZeros() == " << int(count) << " / " << (other.rows()) << "\n"; + } + res.finalize(); + other = res.markAsRValue(); + } +}; + +} // end namespace internal + +template<typename ExpressionType,int Mode> +template<typename OtherDerived> +void SparseTriangularView<ExpressionType,Mode>::solveInPlace(SparseMatrixBase<OtherDerived>& other) const +{ + eigen_assert(m_matrix.cols() == m_matrix.rows() && m_matrix.cols() == other.rows()); + eigen_assert( (!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower))); + +// enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit }; + +// typedef typename internal::conditional<copy, +// typename internal::plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::type OtherCopy; +// OtherCopy otherCopy(other.derived()); + + internal::sparse_solve_triangular_sparse_selector<ExpressionType, OtherDerived, Mode>::run(m_matrix, other.derived()); + +// if (copy) +// other = otherCopy; +} + +#ifdef EIGEN2_SUPPORT + +// deprecated stuff: + +/** \deprecated */ +template<typename Derived> +template<typename OtherDerived> +void SparseMatrixBase<Derived>::solveTriangularInPlace(MatrixBase<OtherDerived>& other) const +{ + this->template triangular<Flags&(Upper|Lower)>().solveInPlace(other); +} + +/** \deprecated */ +template<typename Derived> +template<typename OtherDerived> +typename internal::plain_matrix_type_column_major<OtherDerived>::type +SparseMatrixBase<Derived>::solveTriangular(const MatrixBase<OtherDerived>& other) const +{ + typename internal::plain_matrix_type_column_major<OtherDerived>::type res(other); + derived().solveTriangularInPlace(res); + return res; +} +#endif // EIGEN2_SUPPORT + +} // end namespace Eigen + +#endif // EIGEN_SPARSETRIANGULARSOLVER_H |