|
|
|
|
@@ -32,7 +32,7 @@ Index QuickSplit(VectorV &row, VectorI &ind, Index ncut)
|
|
|
|
|
using std::swap;
|
|
|
|
|
using std::abs;
|
|
|
|
|
Index mid;
|
|
|
|
|
Index n = row.size(); /* length of the vector */
|
|
|
|
|
Index n = convert_index<Index>(row.size()); /* length of the vector */
|
|
|
|
|
Index first, last ;
|
|
|
|
|
|
|
|
|
|
ncut--; /* to fit the zero-based indices */
|
|
|
|
|
@@ -105,7 +105,7 @@ class IncompleteLUT : public SparseSolverBase<IncompleteLUT<_Scalar> >
|
|
|
|
|
typedef Matrix<Scalar,Dynamic,1> Vector;
|
|
|
|
|
typedef SparseMatrix<Scalar,RowMajor> FactorType;
|
|
|
|
|
typedef SparseMatrix<Scalar,ColMajor> PermutType;
|
|
|
|
|
typedef typename FactorType::Index Index;
|
|
|
|
|
typedef typename FactorType::StorageIndex StorageIndex;
|
|
|
|
|
|
|
|
|
|
public:
|
|
|
|
|
typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
|
|
|
|
|
@@ -124,9 +124,9 @@ class IncompleteLUT : public SparseSolverBase<IncompleteLUT<_Scalar> >
|
|
|
|
|
compute(mat);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
Index rows() const { return m_lu.rows(); }
|
|
|
|
|
StorageIndex rows() const { return m_lu.rows(); }
|
|
|
|
|
|
|
|
|
|
Index cols() const { return m_lu.cols(); }
|
|
|
|
|
StorageIndex cols() const { return m_lu.cols(); }
|
|
|
|
|
|
|
|
|
|
/** \brief Reports whether previous computation was successful.
|
|
|
|
|
*
|
|
|
|
|
@@ -189,8 +189,8 @@ protected:
|
|
|
|
|
bool m_analysisIsOk;
|
|
|
|
|
bool m_factorizationIsOk;
|
|
|
|
|
ComputationInfo m_info;
|
|
|
|
|
PermutationMatrix<Dynamic,Dynamic,Index> m_P; // Fill-reducing permutation
|
|
|
|
|
PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv; // Inverse permutation
|
|
|
|
|
PermutationMatrix<Dynamic,Dynamic,StorageIndex> m_P; // Fill-reducing permutation
|
|
|
|
|
PermutationMatrix<Dynamic,Dynamic,StorageIndex> m_Pinv; // Inverse permutation
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
/**
|
|
|
|
|
@@ -218,14 +218,14 @@ template<typename _MatrixType>
|
|
|
|
|
void IncompleteLUT<Scalar>::analyzePattern(const _MatrixType& amat)
|
|
|
|
|
{
|
|
|
|
|
// Compute the Fill-reducing permutation
|
|
|
|
|
SparseMatrix<Scalar,ColMajor, Index> mat1 = amat;
|
|
|
|
|
SparseMatrix<Scalar,ColMajor, Index> mat2 = amat.transpose();
|
|
|
|
|
SparseMatrix<Scalar,ColMajor, StorageIndex> mat1 = amat;
|
|
|
|
|
SparseMatrix<Scalar,ColMajor, StorageIndex> mat2 = amat.transpose();
|
|
|
|
|
// Symmetrize the pattern
|
|
|
|
|
// FIXME for a matrix with nearly symmetric pattern, mat2+mat1 is the appropriate choice.
|
|
|
|
|
// on the other hand for a really non-symmetric pattern, mat2*mat1 should be prefered...
|
|
|
|
|
SparseMatrix<Scalar,ColMajor, Index> AtA = mat2 + mat1;
|
|
|
|
|
SparseMatrix<Scalar,ColMajor, StorageIndex> AtA = mat2 + mat1;
|
|
|
|
|
AtA.prune(keep_diag());
|
|
|
|
|
internal::minimum_degree_ordering<Scalar, Index>(AtA, m_P); // Then compute the AMD ordering...
|
|
|
|
|
internal::minimum_degree_ordering<Scalar, StorageIndex>(AtA, m_P); // Then compute the AMD ordering...
|
|
|
|
|
|
|
|
|
|
m_Pinv = m_P.inverse(); // ... and the inverse permutation
|
|
|
|
|
|
|
|
|
|
@@ -241,7 +241,7 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat)
|
|
|
|
|
using std::abs;
|
|
|
|
|
|
|
|
|
|
eigen_assert((amat.rows() == amat.cols()) && "The factorization should be done on a square matrix");
|
|
|
|
|
Index n = amat.cols(); // Size of the matrix
|
|
|
|
|
StorageIndex n = amat.cols(); // Size of the matrix
|
|
|
|
|
m_lu.resize(n,n);
|
|
|
|
|
// Declare Working vectors and variables
|
|
|
|
|
Vector u(n) ; // real values of the row -- maximum size is n --
|
|
|
|
|
@@ -250,7 +250,7 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat)
|
|
|
|
|
|
|
|
|
|
// Apply the fill-reducing permutation
|
|
|
|
|
eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
|
|
|
|
|
SparseMatrix<Scalar,RowMajor, Index> mat;
|
|
|
|
|
SparseMatrix<Scalar,RowMajor, StorageIndex> mat;
|
|
|
|
|
mat = amat.twistedBy(m_Pinv);
|
|
|
|
|
|
|
|
|
|
// Initialization
|
|
|
|
|
@@ -259,21 +259,21 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat)
|
|
|
|
|
u.fill(0);
|
|
|
|
|
|
|
|
|
|
// number of largest elements to keep in each row:
|
|
|
|
|
Index fill_in = static_cast<Index> (amat.nonZeros()*m_fillfactor)/n+1;
|
|
|
|
|
StorageIndex fill_in = static_cast<StorageIndex> (amat.nonZeros()*m_fillfactor)/n+1;
|
|
|
|
|
if (fill_in > n) fill_in = n;
|
|
|
|
|
|
|
|
|
|
// number of largest nonzero elements to keep in the L and the U part of the current row:
|
|
|
|
|
Index nnzL = fill_in/2;
|
|
|
|
|
Index nnzU = nnzL;
|
|
|
|
|
StorageIndex nnzL = fill_in/2;
|
|
|
|
|
StorageIndex nnzU = nnzL;
|
|
|
|
|
m_lu.reserve(n * (nnzL + nnzU + 1));
|
|
|
|
|
|
|
|
|
|
// global loop over the rows of the sparse matrix
|
|
|
|
|
for (Index ii = 0; ii < n; ii++)
|
|
|
|
|
for (StorageIndex ii = 0; ii < n; ii++)
|
|
|
|
|
{
|
|
|
|
|
// 1 - copy the lower and the upper part of the row i of mat in the working vector u
|
|
|
|
|
|
|
|
|
|
Index sizeu = 1; // number of nonzero elements in the upper part of the current row
|
|
|
|
|
Index sizel = 0; // number of nonzero elements in the lower part of the current row
|
|
|
|
|
StorageIndex sizeu = 1; // number of nonzero elements in the upper part of the current row
|
|
|
|
|
StorageIndex sizel = 0; // number of nonzero elements in the lower part of the current row
|
|
|
|
|
ju(ii) = ii;
|
|
|
|
|
u(ii) = 0;
|
|
|
|
|
jr(ii) = ii;
|
|
|
|
|
@@ -282,7 +282,7 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat)
|
|
|
|
|
typename FactorType::InnerIterator j_it(mat, ii); // Iterate through the current row ii
|
|
|
|
|
for (; j_it; ++j_it)
|
|
|
|
|
{
|
|
|
|
|
Index k = j_it.index();
|
|
|
|
|
StorageIndex k = j_it.index();
|
|
|
|
|
if (k < ii)
|
|
|
|
|
{
|
|
|
|
|
// copy the lower part
|
|
|
|
|
@@ -298,7 +298,7 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat)
|
|
|
|
|
else
|
|
|
|
|
{
|
|
|
|
|
// copy the upper part
|
|
|
|
|
Index jpos = ii + sizeu;
|
|
|
|
|
StorageIndex jpos = ii + sizeu;
|
|
|
|
|
ju(jpos) = k;
|
|
|
|
|
u(jpos) = j_it.value();
|
|
|
|
|
jr(k) = jpos;
|
|
|
|
|
@@ -317,19 +317,19 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat)
|
|
|
|
|
rownorm = sqrt(rownorm);
|
|
|
|
|
|
|
|
|
|
// 3 - eliminate the previous nonzero rows
|
|
|
|
|
Index jj = 0;
|
|
|
|
|
Index len = 0;
|
|
|
|
|
StorageIndex jj = 0;
|
|
|
|
|
StorageIndex len = 0;
|
|
|
|
|
while (jj < sizel)
|
|
|
|
|
{
|
|
|
|
|
// In order to eliminate in the correct order,
|
|
|
|
|
// we must select first the smallest column index among ju(jj:sizel)
|
|
|
|
|
Index k;
|
|
|
|
|
Index minrow = ju.segment(jj,sizel-jj).minCoeff(&k); // k is relative to the segment
|
|
|
|
|
StorageIndex k;
|
|
|
|
|
StorageIndex minrow = ju.segment(jj,sizel-jj).minCoeff(&k); // k is relative to the segment
|
|
|
|
|
k += jj;
|
|
|
|
|
if (minrow != ju(jj))
|
|
|
|
|
{
|
|
|
|
|
// swap the two locations
|
|
|
|
|
Index j = ju(jj);
|
|
|
|
|
StorageIndex j = ju(jj);
|
|
|
|
|
swap(ju(jj), ju(k));
|
|
|
|
|
jr(minrow) = jj; jr(j) = k;
|
|
|
|
|
swap(u(jj), u(k));
|
|
|
|
|
@@ -355,11 +355,11 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat)
|
|
|
|
|
for (; ki_it; ++ki_it)
|
|
|
|
|
{
|
|
|
|
|
Scalar prod = fact * ki_it.value();
|
|
|
|
|
Index j = ki_it.index();
|
|
|
|
|
Index jpos = jr(j);
|
|
|
|
|
StorageIndex j = ki_it.index();
|
|
|
|
|
StorageIndex jpos = jr(j);
|
|
|
|
|
if (jpos == -1) // fill-in element
|
|
|
|
|
{
|
|
|
|
|
Index newpos;
|
|
|
|
|
StorageIndex newpos;
|
|
|
|
|
if (j >= ii) // dealing with the upper part
|
|
|
|
|
{
|
|
|
|
|
newpos = ii + sizeu;
|
|
|
|
|
@@ -388,7 +388,7 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat)
|
|
|
|
|
} // end of the elimination on the row ii
|
|
|
|
|
|
|
|
|
|
// reset the upper part of the pointer jr to zero
|
|
|
|
|
for(Index k = 0; k <sizeu; k++) jr(ju(ii+k)) = -1;
|
|
|
|
|
for(StorageIndex k = 0; k <sizeu; k++) jr(ju(ii+k)) = -1;
|
|
|
|
|
|
|
|
|
|
// 4 - partially sort and insert the elements in the m_lu matrix
|
|
|
|
|
|
|
|
|
|
@@ -401,7 +401,7 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat)
|
|
|
|
|
|
|
|
|
|
// store the largest m_fill elements of the L part
|
|
|
|
|
m_lu.startVec(ii);
|
|
|
|
|
for(Index k = 0; k < len; k++)
|
|
|
|
|
for(StorageIndex k = 0; k < len; k++)
|
|
|
|
|
m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k);
|
|
|
|
|
|
|
|
|
|
// store the diagonal element
|
|
|
|
|
@@ -413,7 +413,7 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat)
|
|
|
|
|
// sort the U-part of the row
|
|
|
|
|
// apply the dropping rule first
|
|
|
|
|
len = 0;
|
|
|
|
|
for(Index k = 1; k < sizeu; k++)
|
|
|
|
|
for(StorageIndex k = 1; k < sizeu; k++)
|
|
|
|
|
{
|
|
|
|
|
if(abs(u(ii+k)) > m_droptol * rownorm )
|
|
|
|
|
{
|
|
|
|
|
@@ -429,7 +429,7 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat)
|
|
|
|
|
internal::QuickSplit(uu, juu, len);
|
|
|
|
|
|
|
|
|
|
// store the largest elements of the U part
|
|
|
|
|
for(Index k = ii + 1; k < ii + len; k++)
|
|
|
|
|
for(StorageIndex k = ii + 1; k < ii + len; k++)
|
|
|
|
|
m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|