aboutsummaryrefslogtreecommitdiffhomepage
path: root/third_party/eigen3/Eigen/src/SparseCore/TriangularSolver.h
diff options
context:
space:
mode:
Diffstat (limited to 'third_party/eigen3/Eigen/src/SparseCore/TriangularSolver.h')
-rw-r--r--third_party/eigen3/Eigen/src/SparseCore/TriangularSolver.h334
1 files changed, 0 insertions, 334 deletions
diff --git a/third_party/eigen3/Eigen/src/SparseCore/TriangularSolver.h b/third_party/eigen3/Eigen/src/SparseCore/TriangularSolver.h
deleted file mode 100644
index cb8ad82b4f..0000000000
--- a/third_party/eigen3/Eigen/src/SparseCore/TriangularSolver.h
+++ /dev/null
@@ -1,334 +0,0 @@
-// 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
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-30 03:03:07 +0000