mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-30 17:40:05 +08:00
Code Issues Packages Projects Releases Wiki Activity

bug #963: make IncompleteLUT compatible with non-default storage index types.

Browse Source
This commit is contained in:
Gael Guennebaud 2015-03-09 13:55:20 +01:00
parent cf9940e17b
commit 224a1fe4c6
2 changed files with 28 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

41
Eigen/src/IterativeLinearSolvers/IncompleteLUT.h
View File

@ -93,21 +93,23 @@ Index QuickSplit(VectorV &row, VectorI &ind, Index ncut)
* alternatively, on GMANE:
* http://comments.gmane.org/gmane.comp.lib.eigen/3302
*/
template <typename _Scalar>
class IncompleteLUT : public SparseSolverBase<IncompleteLUT<_Scalar> >
template <typename _Scalar, typename _StorageIndex = int>
class IncompleteLUT : public SparseSolverBase<IncompleteLUT<_Scalar, _StorageIndex> >
{
protected:
typedef SparseSolverBase<IncompleteLUT<_Scalar> > Base;
typedef SparseSolverBase<IncompleteLUT> Base;
using Base::m_isInitialized;
public:
typedef _Scalar Scalar;
typedef _StorageIndex StorageIndex;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef Matrix<Scalar,Dynamic,1> Vector;
typedef SparseMatrix<Scalar,RowMajor> FactorType;
typedef SparseMatrix<Scalar,ColMajor> PermutType;
typedef typename FactorType::StorageIndex StorageIndex;
typedef Matrix<StorageIndex,Dynamic,1> VectorI;
typedef SparseMatrix<Scalar,RowMajor,StorageIndex> FactorType;
public:
// this typedef is only to export the scalar type and compile-time dimensions to solve_retval
typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
IncompleteLUT()
@ -151,7 +153,7 @@ class IncompleteLUT : public SparseSolverBase<IncompleteLUT<_Scalar> >
*
**/
template<typename MatrixType>
IncompleteLUT<Scalar>& compute(const MatrixType& amat)
IncompleteLUT& compute(const MatrixType& amat)
{
analyzePattern(amat);
factorize(amat);
@ -197,8 +199,8 @@ protected:
* Set control parameter droptol
* \param droptol Drop any element whose magnitude is less than this tolerance
**/
template<typename Scalar>
void IncompleteLUT<Scalar>::setDroptol(const RealScalar& droptol)
template<typename Scalar, typename StorageIndex>
void IncompleteLUT<Scalar,StorageIndex>::setDroptol(const RealScalar& droptol)
{
this->m_droptol = droptol;
}
@ -207,15 +209,15 @@ void IncompleteLUT<Scalar>::setDroptol(const RealScalar& droptol)
* Set control parameter fillfactor
* \param fillfactor This is used to compute the number @p fill_in of largest elements to keep on each row.
**/
template<typename Scalar>
void IncompleteLUT<Scalar>::setFillfactor(int fillfactor)
template<typename Scalar, typename StorageIndex>
void IncompleteLUT<Scalar,StorageIndex>::setFillfactor(int fillfactor)
{
this->m_fillfactor = fillfactor;
}
template <typename Scalar>
template <typename Scalar, typename StorageIndex>
template<typename _MatrixType>
void IncompleteLUT<Scalar>::analyzePattern(const _MatrixType& amat)
void IncompleteLUT<Scalar,StorageIndex>::analyzePattern(const _MatrixType& amat)
{
// Compute the Fill-reducing permutation
SparseMatrix<Scalar,ColMajor, StorageIndex> mat1 = amat;
@ -232,9 +234,9 @@ void IncompleteLUT<Scalar>::analyzePattern(const _MatrixType& amat)
m_analysisIsOk = true;
}
template <typename Scalar>
template <typename Scalar, typename StorageIndex>
template<typename _MatrixType>
void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat)
void IncompleteLUT<Scalar,StorageIndex>::factorize(const _MatrixType& amat)
{
using std::sqrt;
using std::swap;
@ -246,8 +248,8 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat)
m_lu.resize(n,n);
// Declare Working vectors and variables
Vector u(n) ; // real values of the row -- maximum size is n --
VectorXi ju(n); // column position of the values in u -- maximum size is n
VectorXi jr(n); // Indicate the position of the nonzero elements in the vector u -- A zero location is indicated by -1
VectorI ju(n); // column position of the values in u -- maximum size is n
VectorI jr(n); // Indicate the position of the nonzero elements in the vector u -- A zero location is indicated by -1
// Apply the fill-reducing permutation
eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
@ -398,7 +400,7 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat)
sizel = len;
len = (std::min)(sizel, nnzL);
typename Vector::SegmentReturnType ul(u.segment(0, sizel));
typename VectorXi::SegmentReturnType jul(ju.segment(0, sizel));
typename VectorI::SegmentReturnType jul(ju.segment(0, sizel));
internal::QuickSplit(ul, jul, len);
// store the largest m_fill elements of the L part
@ -427,14 +429,13 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat)
sizeu = len + 1; // +1 to take into account the diagonal element
len = (std::min)(sizeu, nnzU);
typename Vector::SegmentReturnType uu(u.segment(ii+1, sizeu-1));
typename VectorXi::SegmentReturnType juu(ju.segment(ii+1, sizeu-1));
typename VectorI::SegmentReturnType juu(ju.segment(ii+1, sizeu-1));
internal::QuickSplit(uu, juu, len);
// store the largest elements of the U part
for(Index k = ii + 1; k < ii + len; k++)
m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k);
}
m_lu.finalize();
m_lu.makeCompressed();

13
test/bicgstab.cpp
View File

@ -10,11 +10,11 @@
#include "sparse_solver.h"
#include <Eigen/IterativeLinearSolvers>
template<typename T> void test_bicgstab_T()
template<typename T, typename I> void test_bicgstab_T()
{
BiCGSTAB<SparseMatrix<T>, DiagonalPreconditioner<T> > bicgstab_colmajor_diag;
BiCGSTAB<SparseMatrix<T>, IdentityPreconditioner > bicgstab_colmajor_I;
BiCGSTAB<SparseMatrix<T>, IncompleteLUT<T> > bicgstab_colmajor_ilut;
BiCGSTAB<SparseMatrix<T,0,I>, DiagonalPreconditioner<T> > bicgstab_colmajor_diag;
BiCGSTAB<SparseMatrix<T,0,I>, IdentityPreconditioner > bicgstab_colmajor_I;
BiCGSTAB<SparseMatrix<T,0,I>, IncompleteLUT<T,I> > bicgstab_colmajor_ilut;
//BiCGSTAB<SparseMatrix<T>, SSORPreconditioner<T> > bicgstab_colmajor_ssor;
CALL_SUBTEST( check_sparse_square_solving(bicgstab_colmajor_diag) );
@ -25,6 +25,7 @@ template<typename T> void test_bicgstab_T()
void test_bicgstab()
{
CALL_SUBTEST_1(test_bicgstab_T<double>());
CALL_SUBTEST_2(test_bicgstab_T<std::complex<double> >());
CALL_SUBTEST_1((test_bicgstab_T<double,int>()) );
CALL_SUBTEST_1((test_bicgstab_T<double,long int>()));
CALL_SUBTEST_2((test_bicgstab_T<std::complex<double>, int>()));
}