diff --git a/Eigen/src/OrderingMethods/Ordering.h b/Eigen/src/OrderingMethods/Ordering.h index 3751f9bee..670cca9c4 100644 --- a/Eigen/src/OrderingMethods/Ordering.h +++ b/Eigen/src/OrderingMethods/Ordering.h @@ -60,7 +60,9 @@ class AMDOrdering public: typedef PermutationMatrix PermutationType; - /** Compute the permutation vector from a column-major sparse matrix */ + /** Compute the permutation vector from a sparse matrix + * This routine is much faster if the input matrix is column-major + */ template void operator()(const MatrixType& mat, PermutationType& perm) { @@ -73,7 +75,7 @@ class AMDOrdering internal::minimum_degree_ordering(symm, perm); } - /** Compute the permutation with a self adjoint matrix */ + /** Compute the permutation with a selfadjoint matrix */ template void operator()(const SparseSelfAdjointView& mat, PermutationType& perm) { @@ -85,6 +87,26 @@ class AMDOrdering } }; +/** + * Get the natural ordering + * + *NOTE Returns an empty permutation matrix + * \tparam Index The type of indices of the matrix + */ +template +class NaturalOrdering +{ + public: + typedef PermutationMatrix PermutationType; + + /** Compute the permutation vector from a column-major sparse matrix */ + template + void operator()(const MatrixType& mat, PermutationType& perm) + { + perm.resize(0); + } + +}; /** * Get the column approximate minimum degree ordering diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index e4a4c3a7b..74f710563 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h @@ -255,7 +255,7 @@ class SparseLU void initperfvalues() { m_panel_size = 12; - m_relax = 1; + m_relax = 6; m_maxsuper = 100; m_rowblk = 200; m_colblk = 60; @@ -320,26 +320,31 @@ void SparseLU::analyzePattern(const MatrixType& mat) // Compute the fill-reducing ordering // TODO Currently, the only available ordering method is AMD. - OrderingType ord(mat); - m_perm_c = ord.get_perm(); + OrderingType ord; + ord(mat,m_perm_c); //FIXME Check the right semantic behind m_perm_c // that is, column j of mat goes to column m_perm_c(j) of mat * m_perm_c; + //DEBUG : Set the natural ordering + for (int i = 0; i < mat.cols(); i++) + m_perm_c.indices()(i) = i; // Apply the permutation to the column of the input matrix - m_mat = mat * m_perm_c; + m_mat = mat * m_perm_c.inverse(); // Compute the column elimination tree of the permuted matrix if (m_etree.size() == 0) m_etree.resize(m_mat.cols()); + LU_sp_coletree(m_mat, m_etree); - + // In symmetric mode, do not do postorder here if (!m_symmetricmode) { IndexVector post, iwork; // Post order etree LU_TreePostorder(m_mat.cols(), m_etree, post); + // Renumber etree in postorder int m = m_mat.cols(); iwork.resize(m+1); @@ -348,12 +353,15 @@ void SparseLU::analyzePattern(const MatrixType& mat) // Postmultiply A*Pc by post, i.e reorder the matrix according to the postorder of the etree - PermutationType post_perm(m);; + PermutationType post_perm(m); //FIXME Use vector constructor for (int i = 0; i < m; i++) post_perm.indices()(i) = post(i); - //m_mat = m_mat * post_perm; // FIXME This should surely be in factorize() + +// m_mat = m_mat * post_perm.inverse(); // FIXME This should surely be in factorize() + // Composition of the two permutations m_perm_c = m_perm_c * post_perm; + } // end postordering m_analysisIsOk = true; @@ -402,9 +410,14 @@ void SparseLU::factorize(const MatrixType& matrix) // Apply the column permutation computed in analyzepattern() - m_mat = matrix * m_perm_c; + m_mat = matrix * m_perm_c.inverse(); m_mat.makeCompressed(); + // DEBUG ... Watch matrix permutation + const int *asub_in = matrix.innerIndexPtr(); + const int *colptr_in = matrix.outerIndexPtr(); + int * asub = m_mat.innerIndexPtr(); + int * colptr = m_mat.outerIndexPtr(); int m = m_mat.rows(); int n = m_mat.cols(); int nnz = m_mat.nonZeros(); @@ -455,7 +468,8 @@ void SparseLU::factorize(const MatrixType& matrix) // Setup Permutation vectors // Compute the inverse of perm_c - PermutationType iperm_c (m_perm_c.inverse() ); +// PermutationType iperm_c (m_perm_c.inverse() ); + PermutationType iperm_c (m_perm_c); // Identify initial relaxed snodes IndexVector relax_end(n); @@ -464,6 +478,9 @@ void SparseLU::factorize(const MatrixType& matrix) else LU_relax_snode(n, m_etree, m_relax, marker, relax_end); + //DEBUG +// std::cout<< "relax_end " <= nzlmax ) { mem = LUMemXpand(lsub, nzlmax, nextl, LSUB, glu.num_expansions); - if (mem) return mem; + if (mem) return mem; // Memory expansion failed... Return the memory allocated so far } } } @@ -100,7 +100,7 @@ while (new_next > nzlmax) { mem = LUMemXpand(lsub, nzlmax, nextl, LSUB, glu.num_expansions); - if (mem) return mem; + if (mem) return mem; // Memory expansion failed... Return the memory allocated so far } Index ifrom, ito = nextl; for (ifrom = xlsub(jcol); ifrom < nextl;) diff --git a/Eigen/src/SuperLUSupport/SuperLUSupport.h b/Eigen/src/SuperLUSupport/SuperLUSupport.h index 60a3eb09a..9c2e6e17e 100644 --- a/Eigen/src/SuperLUSupport/SuperLUSupport.h +++ b/Eigen/src/SuperLUSupport/SuperLUSupport.h @@ -627,6 +627,9 @@ void SuperLU::factorize(const MatrixType& a) this->initFactorization(a); + //DEBUG + m_sluOptions.ColPerm = NATURAL; + m_sluOptions.Equil = NO; int info = 0; RealScalar recip_pivot_growth, rcond; RealScalar ferr, berr; diff --git a/bench/spbench/test_sparseLU.cpp b/bench/spbench/test_sparseLU.cpp index 4727cc12b..841011f30 100644 --- a/bench/spbench/test_sparseLU.cpp +++ b/bench/spbench/test_sparseLU.cpp @@ -17,7 +17,7 @@ int main(int argc, char **args) typedef Matrix DenseMatrix; typedef Matrix DenseRhs; VectorXd b, x, tmp; - SparseLU, AMDOrdering > solver; + SparseLU, AMDOrdering > solver; ifstream matrix_file; string line; int n; @@ -52,7 +52,7 @@ int main(int argc, char **args) } /* Compute the factorization */ - solver.isSymmetric(true); +// solver.isSymmetric(true); solver.compute(A); solver._solve(b, x);