aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Eigen/src/OrderingMethods/Ordering.h26
-rw-r--r--Eigen/src/SparseLU/SparseLU.h35
-rw-r--r--Eigen/src/SparseLU/SparseLU_Coletree.h9
-rw-r--r--Eigen/src/SparseLU/SparseLU_relax_snode.h5
-rw-r--r--Eigen/src/SparseLU/SparseLU_snode_dfs.h8
-rw-r--r--Eigen/src/SuperLUSupport/SuperLUSupport.h3
-rw-r--r--bench/spbench/test_sparseLU.cpp4
7 files changed, 67 insertions, 23 deletions
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<Dynamic, Dynamic, Index> 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 <typename MatrixType>
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 <typename SrcType, unsigned int SrcUpLo>
void operator()(const SparseSelfAdjointView<SrcType, SrcUpLo>& 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 <typename Index>
+class NaturalOrdering
+{
+ public:
+ typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
+
+ /** Compute the permutation vector from a column-major sparse matrix */
+ template <typename MatrixType>
+ 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<MatrixType, OrderingType>::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<MatrixType, OrderingType>::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<MatrixType, OrderingType>::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<MatrixType, OrderingType>::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<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
else
LU_relax_snode<IndexVector>(n, m_etree, m_relax, marker, relax_end);
+ //DEBUG
+// std::cout<< "relax_end " <<relax_end.transpose() << std::endl;
+
m_perm_r.resize(m);
m_perm_r.indices().setConstant(-1); //FIXME
marker.setConstant(-1);
diff --git a/Eigen/src/SparseLU/SparseLU_Coletree.h b/Eigen/src/SparseLU/SparseLU_Coletree.h
index 585b02fdf..1329d383f 100644
--- a/Eigen/src/SparseLU/SparseLU_Coletree.h
+++ b/Eigen/src/SparseLU/SparseLU_Coletree.h
@@ -75,7 +75,6 @@ int LU_sp_coletree(const MatrixType& mat, IndexVector& parent)
IndexVector pp(nc); // disjoint sets
pp.setZero(); // Initialize disjoint sets
IndexVector firstcol(nr); // First nonzero column in each row
- firstcol.setZero();
//Compute first nonzero column in each row
int row,col;
@@ -95,7 +94,9 @@ int LU_sp_coletree(const MatrixType& mat, IndexVector& parent)
int rset, cset, rroot;
for (col = 0; col < nc; col++)
{
- cset = pp(col) = col; // Initially, each element is in its own set //FIXME
+// cset = pp(col) = col; // Initially, each element is in its own set //FIXME
+ pp(col) = col;
+ cset = col;
root(cset) = col;
parent(col) = nc;
for (typename MatrixType::InnerIterator it(mat, col); it; ++it)
@@ -107,7 +108,9 @@ int LU_sp_coletree(const MatrixType& mat, IndexVector& parent)
if (rroot != col)
{
parent(rroot) = col;
- cset = pp(cset) = rset; // Get the union of cset and rset //FIXME
+// cset = pp(cset) = rset; // Get the union of cset and rset //FIXME
+ pp(cset) = rset;
+ cset = rset;
root(cset) = col;
}
}
diff --git a/Eigen/src/SparseLU/SparseLU_relax_snode.h b/Eigen/src/SparseLU/SparseLU_relax_snode.h
index 0006dde33..5123e94bf 100644
--- a/Eigen/src/SparseLU/SparseLU_relax_snode.h
+++ b/Eigen/src/SparseLU/SparseLU_relax_snode.h
@@ -57,7 +57,7 @@ void LU_relax_snode (const int n, IndexVector& et, const int relax_columns, Inde
{
// compute the number of descendants of each node in the etree
- register int j, parent;
+ int j, parent;
relax_end.setConstant(IND_EMPTY);
descendants.setZero();
for (j = 0; j < n; j++)
@@ -66,9 +66,8 @@ void LU_relax_snode (const int n, IndexVector& et, const int relax_columns, Inde
if (parent != n) // not the dummy root
descendants(parent) += descendants(j) + 1;
}
-
// Identify the relaxed supernodes by postorder traversal of the etree
- register int snode_start; // beginning of a snode
+ int snode_start; // beginning of a snode
for (j = 0; j < n; )
{
parent = et(j);
diff --git a/Eigen/src/SparseLU/SparseLU_snode_dfs.h b/Eigen/src/SparseLU/SparseLU_snode_dfs.h
index 3e7033c67..6b2817262 100644
--- a/Eigen/src/SparseLU/SparseLU_snode_dfs.h
+++ b/Eigen/src/SparseLU/SparseLU_snode_dfs.h
@@ -68,8 +68,8 @@
Index& nzlmax = glu.nzlmax;
int mem;
Index nsuper = ++supno(jcol); // Next available supernode number
- register int nextl = xlsub(jcol); //Index of the starting location of the jcol-th column in lsub
- register int i,k;
+ 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++)
{
@@ -86,7 +86,7 @@
if( nextl >= nzlmax )
{
mem = LUMemXpand<IndexVector>(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<IndexVector>(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<MatrixType>::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<double, Dynamic, Dynamic> DenseMatrix;
typedef Matrix<double, Dynamic, 1> DenseRhs;
VectorXd b, x, tmp;
- SparseLU<SparseMatrix<double, ColMajor>, AMDOrdering<double, int> > solver;
+ SparseLU<SparseMatrix<double, ColMajor>, AMDOrdering<int> > 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);
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-29 15:52:20 +0000