From ce30d50e3ed9723ed3ecd38e7c99661730c12813 Mon Sep 17 00:00:00 2001 From: "Desire NUENTSA W." Date: Fri, 27 Jul 2012 16:38:20 +0200 Subject: Improve the permutation --- Eigen/src/SparseCore/SparseMatrix.h | 12 ++++++++++++ Eigen/src/SparseLU/SparseLU.h | 26 +++++++++++++++++++++----- Eigen/src/SparseLU/SparseLU_snode_dfs.h | 13 ++++++------- 3 files changed, 39 insertions(+), 12 deletions(-) (limited to 'Eigen') diff --git a/Eigen/src/SparseCore/SparseMatrix.h b/Eigen/src/SparseCore/SparseMatrix.h index 214f130f5..52a9dab70 100644 --- a/Eigen/src/SparseCore/SparseMatrix.h +++ b/Eigen/src/SparseCore/SparseMatrix.h @@ -477,6 +477,18 @@ class SparseMatrix m_data.squeeze(); } + /** Turns the matrix into the uncompressed mode */ + void Uncompress() + { + if(m_innerNonZeros != 0) + return; + m_innerNonZeros = new Index[m_outerSize]; + for (int i = 0; i < m_outerSize; i++) + { + m_innerNonZeros[i] = m_outerIndex[i+1] - m_outerIndex[i]; + } + } + /** Suppresses all nonzeros which are \b much \b smaller \b than \a reference under the tolerence \a epsilon */ void prune(Scalar reference, RealScalar epsilon = NumTraits::dummy_precision()) { diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index 474dfdedc..70898958b 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h @@ -346,8 +346,17 @@ void SparseLU::analyzePattern(const MatrixType& mat) // Apply the permutation to the column of the input matrix - m_mat = mat * m_perm_c.inverse(); //FIXME It should be less expensive here to permute only the structural pattern of the matrix - +// m_mat = mat * m_perm_c.inverse(); //FIXME It should be less expensive here to permute only the structural pattern of the matrix + + //First copy the whole input matrix. + m_mat = mat; + m_mat.Uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers. FIXME : This vector is filled but not subsequently used. + //Then, permute only the column pointers + for (int i = 0; i < mat.cols(); i++) + { + m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = mat.outerIndexPtr()[i]; + m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = mat.outerIndexPtr()[i+1] - mat.outerIndexPtr()[i]; + } // Compute the column elimination tree of the permuted matrix if (m_etree.size() == 0) m_etree.resize(m_mat.cols()); @@ -424,8 +433,15 @@ void SparseLU::factorize(const MatrixType& matrix) // Apply the column permutation computed in analyzepattern() - m_mat = matrix * m_perm_c.inverse(); - m_mat.makeCompressed(); + // m_mat = matrix * m_perm_c.inverse(); + m_mat = matrix; + m_mat.Uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers. + //Then, permute only the column pointers + for (int i = 0; i < matrix.cols(); i++) + { + m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = matrix.outerIndexPtr()[i]; + m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = matrix.outerIndexPtr()[i+1] - matrix.outerIndexPtr()[i]; + } int m = m_mat.rows(); int n = m_mat.cols(); @@ -504,7 +520,7 @@ void SparseLU::factorize(const MatrixType& matrix) // Factorize the relaxed supernode(jcol:kcol) // First, determine the union of the row structure of the snode - info = LU_snode_dfs(jcol, kcol, m_mat.innerIndexPtr(), m_mat.outerIndexPtr(), xprune, marker, m_glu); + info = LU_snode_dfs(jcol, kcol, m_mat, xprune, marker, m_glu); if ( info ) { std::cerr << "MEMORY ALLOCATION FAILED IN SNODE_DFS() \n"; diff --git a/Eigen/src/SparseLU/SparseLU_snode_dfs.h b/Eigen/src/SparseLU/SparseLU_snode_dfs.h index 6b2817262..150d9d0ef 100644 --- a/Eigen/src/SparseLU/SparseLU_snode_dfs.h +++ b/Eigen/src/SparseLU/SparseLU_snode_dfs.h @@ -57,8 +57,8 @@ * \param marker (in/out) working vector * \return 0 on success, > 0 size of the memory when memory allocation failed */ - template - int LU_snode_dfs(const int jcol, const int kcol, const typename IndexVector::Scalar* asub, const typename IndexVector::Scalar* colptr, IndexVector& xprune, IndexVector& marker, LU_GlobalLU_t& glu) + template + int LU_snode_dfs(const int jcol, const int kcol,const MatrixType& mat, IndexVector& xprune, IndexVector& marker, LU_GlobalLU_t& glu) { typedef typename IndexVector::Scalar Index; IndexVector& xsup = glu.xsup; @@ -69,14 +69,13 @@ int mem; Index nsuper = ++supno(jcol); // Next available supernode number int nextl = xlsub(jcol); //Index of the starting location of the jcol-th column in lsub - int i,k; int krow,kmark; - for (i = jcol; i <=kcol; i++) + for (int i = jcol; i <=kcol; i++) { // For each nonzero in A(*,i) - for (k = colptr[i]; k < colptr[i+1]; k++) + for (typename MatrixType::InnerIterator it(mat, i); it; ++it) { - krow = asub[k]; + krow = it.row(); kmark = marker(krow); if ( kmark != kcol ) { @@ -105,7 +104,7 @@ Index ifrom, ito = nextl; for (ifrom = xlsub(jcol); ifrom < nextl;) lsub(ito++) = lsub(ifrom++); - for (i = jcol+1; i <=kcol; i++) xlsub(i) = nextl; + for (int i = jcol+1; i <=kcol; i++) xlsub(i) = nextl; nextl = ito; } xsup(nsuper+1) = kcol + 1; // Start of next available supernode -- cgit v1.2.3