aboutsummaryrefslogtreecommitdiff
path: root/Eigen/src/SparseCore/TriangularSolver.h
diff options
context:
space:
mode:
Diffstat (limited to 'Eigen/src/SparseCore/TriangularSolver.h')
-rw-r--r--Eigen/src/SparseCore/TriangularSolver.h117
1 files changed, 49 insertions, 68 deletions
diff --git a/Eigen/src/SparseCore/TriangularSolver.h b/Eigen/src/SparseCore/TriangularSolver.h
index ccc12af79..f9c56ba79 100644
--- a/Eigen/src/SparseCore/TriangularSolver.h
+++ b/Eigen/src/SparseCore/TriangularSolver.h
@@ -28,16 +28,19 @@ template<typename Lhs, typename Rhs, int Mode>
struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Lower,RowMajor>
{
typedef typename Rhs::Scalar Scalar;
+ typedef evaluator<Lhs> LhsEval;
+ typedef typename evaluator<Lhs>::InnerIterator LhsIterator;
static void run(const Lhs& lhs, Rhs& other)
{
- for(int col=0 ; col<other.cols() ; ++col)
+ LhsEval lhsEval(lhs);
+ for(Index col=0 ; col<other.cols() ; ++col)
{
- for(int i=0; i<lhs.rows(); ++i)
+ for(Index 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)
+ Index lastIndex = 0;
+ for(LhsIterator it(lhsEval, i); it; ++it)
{
lastVal = it.value();
lastIndex = it.index();
@@ -62,15 +65,18 @@ template<typename Lhs, typename Rhs, int Mode>
struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,RowMajor>
{
typedef typename Rhs::Scalar Scalar;
+ typedef evaluator<Lhs> LhsEval;
+ typedef typename evaluator<Lhs>::InnerIterator LhsIterator;
static void run(const Lhs& lhs, Rhs& other)
{
- for(int col=0 ; col<other.cols() ; ++col)
+ LhsEval lhsEval(lhs);
+ for(Index col=0 ; col<other.cols() ; ++col)
{
- for(int i=lhs.rows()-1 ; i>=0 ; --i)
+ for(Index i=lhs.rows()-1 ; i>=0 ; --i)
{
Scalar tmp = other.coeff(i,col);
Scalar l_ii(0);
- typename Lhs::InnerIterator it(lhs, i);
+ LhsIterator it(lhsEval, i);
while(it && it.index()<i)
++it;
if(!(Mode & UnitDiag))
@@ -86,10 +92,8 @@ struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,RowMajor>
tmp -= it.value() * other.coeff(it.index(),col);
}
- if (Mode & UnitDiag)
- other.coeffRef(i,col) = tmp;
- else
- other.coeffRef(i,col) = tmp/l_ii;
+ if (Mode & UnitDiag) other.coeffRef(i,col) = tmp;
+ else other.coeffRef(i,col) = tmp/l_ii;
}
}
}
@@ -100,16 +104,19 @@ template<typename Lhs, typename Rhs, int Mode>
struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Lower,ColMajor>
{
typedef typename Rhs::Scalar Scalar;
+ typedef evaluator<Lhs> LhsEval;
+ typedef typename evaluator<Lhs>::InnerIterator LhsIterator;
static void run(const Lhs& lhs, Rhs& other)
{
- for(int col=0 ; col<other.cols() ; ++col)
+ LhsEval lhsEval(lhs);
+ for(Index col=0 ; col<other.cols() ; ++col)
{
- for(int i=0; i<lhs.cols(); ++i)
+ for(Index 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);
+ LhsIterator it(lhsEval, i);
while(it && it.index()<i)
++it;
if(!(Mode & UnitDiag))
@@ -132,11 +139,14 @@ template<typename Lhs, typename Rhs, int Mode>
struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,ColMajor>
{
typedef typename Rhs::Scalar Scalar;
+ typedef evaluator<Lhs> LhsEval;
+ typedef typename evaluator<Lhs>::InnerIterator LhsIterator;
static void run(const Lhs& lhs, Rhs& other)
{
- for(int col=0 ; col<other.cols() ; ++col)
+ LhsEval lhsEval(lhs);
+ for(Index col=0 ; col<other.cols() ; ++col)
{
- for(int i=lhs.cols()-1; i>=0; --i)
+ for(Index i=lhs.cols()-1; i>=0; --i)
{
Scalar& tmp = other.coeffRef(i,col);
if (tmp!=Scalar(0)) // optimization when other is actually sparse
@@ -144,13 +154,13 @@ struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,ColMajor>
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);
+ LhsIterator it(lhsEval, i);
while(it && it.index()!=i)
- --it;
+ ++it;
eigen_assert(it && it.index()==i);
other.coeffRef(i,col) /= it.value();
}
- typename Lhs::InnerIterator it(lhs, i);
+ LhsIterator it(lhsEval, i);
for(; it && it.index()<i; ++it)
other.coeffRef(it.index(), col) -= tmp * it.value();
}
@@ -161,11 +171,13 @@ struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,ColMajor>
} // end namespace internal
-template<typename ExpressionType,int Mode>
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+
+template<typename ExpressionType,unsigned int Mode>
template<typename OtherDerived>
-void SparseTriangularView<ExpressionType,Mode>::solveInPlace(MatrixBase<OtherDerived>& other) const
+void TriangularViewImpl<ExpressionType,Mode,Sparse>::solveInPlace(MatrixBase<OtherDerived>& other) const
{
- eigen_assert(m_matrix.cols() == m_matrix.rows() && m_matrix.cols() == other.rows());
+ eigen_assert(derived().cols() == derived().rows() && derived().cols() == other.rows());
eigen_assert((!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower)));
enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit };
@@ -174,21 +186,12 @@ void SparseTriangularView<ExpressionType,Mode>::solveInPlace(MatrixBase<OtherDer
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);
+ internal::sparse_solve_triangular_selector<ExpressionType, typename internal::remove_reference<OtherCopy>::type, Mode>::run(derived().nestedExpression(), 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;
-}
+#endif
// pure sparse path
@@ -208,18 +211,18 @@ 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;
+ typedef typename promote_index_type<typename traits<Lhs>::StorageIndex,
+ typename traits<Rhs>::StorageIndex>::type StorageIndex;
static void run(const Lhs& lhs, Rhs& other)
{
const bool IsLower = (UpLo==Lower);
- AmbiVector<Scalar,Index> tempVector(other.rows()*2);
+ AmbiVector<Scalar,StorageIndex> 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)
+ for(Index col=0 ; col<other.cols() ; ++col)
{
// FIXME estimate number of non zeros
tempVector.init(.99/*float(other.col(col).nonZeros())/float(other.rows())*/);
@@ -230,7 +233,7 @@ struct sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor>
tempVector.coeffRef(rhsIt.index()) = rhsIt.value();
}
- for(int i=IsLower?0:lhs.cols()-1;
+ for(Index i=IsLower?0:lhs.cols()-1;
IsLower?i<lhs.cols():i>=0;
i+=IsLower?1:-1)
{
@@ -267,9 +270,9 @@ struct sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor>
}
- int count = 0;
+ Index count = 0;
// FIXME compute a reference value to filter zeros
- for (typename AmbiVector<Scalar,Index>::Iterator it(tempVector/*,1e-12*/); it; ++it)
+ for (typename AmbiVector<Scalar,StorageIndex>::Iterator it(tempVector/*,1e-12*/); it; ++it)
{
++ count;
// std::cerr << "fill " << it.index() << ", " << col << "\n";
@@ -286,11 +289,12 @@ struct sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor>
} // end namespace internal
-template<typename ExpressionType,int Mode>
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+template<typename ExpressionType,unsigned int Mode>
template<typename OtherDerived>
-void SparseTriangularView<ExpressionType,Mode>::solveInPlace(SparseMatrixBase<OtherDerived>& other) const
+void TriangularViewImpl<ExpressionType,Mode,Sparse>::solveInPlace(SparseMatrixBase<OtherDerived>& other) const
{
- eigen_assert(m_matrix.cols() == m_matrix.rows() && m_matrix.cols() == other.rows());
+ eigen_assert(derived().cols() == derived().rows() && derived().cols() == other.rows());
eigen_assert( (!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower)));
// enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit };
@@ -299,35 +303,12 @@ void SparseTriangularView<ExpressionType,Mode>::solveInPlace(SparseMatrixBase<Ot
// 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());
+ internal::sparse_solve_triangular_sparse_selector<ExpressionType, OtherDerived, Mode>::run(derived().nestedExpression(), 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
+#endif
} // end namespace Eigen
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-19 17:53:47 +0000