mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-06-04 02:33:59 +08:00
335 lines
9.9 KiB
C++
335 lines
9.9 KiB
C++
// 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
|