diff options
author | Desire NUENTSA <desire.nuentsa_wakam@inria.fr> | 2012-08-06 14:55:02 +0200 |
---|---|---|
committer | Desire NUENTSA <desire.nuentsa_wakam@inria.fr> | 2012-08-06 14:55:02 +0200 |
commit | 4d3b7e2a1351d60b9ee26d0fe3442cd5b3a1f8a9 (patch) | |
tree | 8ad3484157af2e83be037c40634ba58d29a00720 | |
parent | a51806993b1a437af308db9c6893cab71e7ca814 (diff) |
Add support for Metis fill-reducing ordering ; it is generally more efficient than COLAMD ordering
-rw-r--r-- | Eigen/MetisSupport | 26 | ||||
-rw-r--r-- | Eigen/src/MetisSupport/CMakeLists.txt | 6 | ||||
-rw-r--r-- | Eigen/src/MetisSupport/MetisSupport.h | 138 | ||||
-rw-r--r-- | bench/spbench/CMakeLists.txt | 6 | ||||
-rw-r--r-- | bench/spbench/test_sparseLU.cpp | 8 | ||||
-rw-r--r-- | cmake/FindMetis.cmake | 3 |
6 files changed, 186 insertions, 1 deletions
diff --git a/Eigen/MetisSupport b/Eigen/MetisSupport new file mode 100644 index 000000000..a44086ad9 --- /dev/null +++ b/Eigen/MetisSupport @@ -0,0 +1,26 @@ +#ifndef EIGEN_METISSUPPORT_MODULE_H +#define EIGEN_METISSUPPORT_MODULE_H + +#include "SparseCore" + +#include "src/Core/util/DisableStupidWarnings.h" + +extern "C" { +#include <metis.h> +} + + +/** \ingroup Support_modules + * \defgroup MetisSupport_Module MetisSupport module + * + * \code + * #include <Eigen/MetisSupport> + * \endcode + */ + + +#include "src/MetisSupport/MetisSupport.h" + +#include "src/Core/util/ReenableStupidWarnings.h" + +#endif // EIGEN_METISSUPPORT_MODULE_H diff --git a/Eigen/src/MetisSupport/CMakeLists.txt b/Eigen/src/MetisSupport/CMakeLists.txt new file mode 100644 index 000000000..2bad31416 --- /dev/null +++ b/Eigen/src/MetisSupport/CMakeLists.txt @@ -0,0 +1,6 @@ +FILE(GLOB Eigen_MetisSupport_SRCS "*.h") + +INSTALL(FILES + ${Eigen_MetisSupport_SRCS} + DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/MetisSupport COMPONENT Devel + ) diff --git a/Eigen/src/MetisSupport/MetisSupport.h b/Eigen/src/MetisSupport/MetisSupport.h new file mode 100644 index 000000000..a762d96f6 --- /dev/null +++ b/Eigen/src/MetisSupport/MetisSupport.h @@ -0,0 +1,138 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. +#ifndef METIS_SUPPORT_H +#define METIS_SUPPORT_H + +namespace Eigen { +/** + * Get the fill-reducing ordering from the METIS package + * + * If A is the original matrix and Ap is the permuted matrix, + * the fill-reducing permutation is defined as follows : + * Row (column) i of A is the matperm(i) row (column) of Ap. + * WARNING: As computed by METIS, this corresponds to the vector iperm (instead of perm) + */ +template <typename Index> +class MetisOrdering +{ +public: + typedef PermutationMatrix<Dynamic,Dynamic,Index> PermutationType; + typedef Matrix<Index,Dynamic,1> IndexVector; + + template <typename MatrixType> + void get_symmetrized_graph(const MatrixType& A) + { + Index m = A.cols(); + + // Get the transpose of the input matrix + MatrixType At = A.transpose(); + // Get the number of nonzeros elements in each row/col of At+A + Index TotNz = 0; + IndexVector visited(m); + visited.setConstant(-1); + for (int j = 0; j < m; j++) + { + // Compute the union structure of of A(j,:) and At(j,:) + visited(j) = j; // Do not include the diagonal element + // Get the nonzeros in row/column j of A + for (typename MatrixType::InnerIterator it(A, j); it; ++it) + { + Index idx = it.index(); // Get the row index (for column major) or column index (for row major) + if (visited(idx) != j ) + { + visited(idx) = j; + ++TotNz; + } + } + //Get the nonzeros in row/column j of At + for (typename MatrixType::InnerIterator it(At, j); it; ++it) + { + Index idx = it.index(); + if(visited(idx) != j) + { + visited(idx) = j; + ++TotNz; + } + } + } + // Reserve place for A + At + m_indexPtr.resize(m+1); + m_innerIndices.resize(TotNz); + + // Now compute the real adjacency list of each column/row + visited.setConstant(-1); + Index CurNz = 0; + for (int j = 0; j < m; j++) + { + m_indexPtr(j) = CurNz; + + visited(j) = j; // Do not include the diagonal element + // Add the pattern of row/column j of A to A+At + for (typename MatrixType::InnerIterator it(A,j); it; ++it) + { + Index idx = it.index(); // Get the row index (for column major) or column index (for row major) + if (visited(idx) != j ) + { + visited(idx) = j; + m_innerIndices(CurNz) = idx; + CurNz++; + } + } + //Add the pattern of row/column j of At to A+At + for (typename MatrixType::InnerIterator it(At, j); it; ++it) + { + Index idx = it.index(); + if(visited(idx) != j) + { + visited(idx) = j; + m_innerIndices(CurNz) = idx; + ++CurNz; + } + } + } + m_indexPtr(m) = CurNz; + } + + template <typename MatrixType> + void operator() (const MatrixType& A, PermutationType& matperm) + { + Index m = A.cols(); + IndexVector perm(m),iperm(m); + // First, symmetrize the matrix graph. + get_symmetrized_graph(A); + int output_error; + + // Call the fill-reducing routine from METIS + output_error = METIS_NodeND(&m, m_indexPtr.data(), m_innerIndices.data(), NULL, NULL, perm.data(), iperm.data()); + + if(output_error != METIS_OK) + { + //FIXME The ordering interface should define a class of possible errors + std::cerr << "ERROR WHILE CALLING THE METIS PACKAGE \n"; + return; + } + + // Get the fill-reducing permutation + //NOTE: If Ap is the permuted matrix then perm and iperm vectors are defined as follows + // Row (column) i of Ap is the perm(i) row(column) of A, and row (column) i of A is the iperm(i) row(column) of Ap + + // To be consistent with the use of the permutation in SparseLU module, we thus keep the iperm vector + matperm.resize(m); + for (int j = 0; j < m; j++) + matperm.indices()(j) = iperm(j); + + } + + protected: + IndexVector m_indexPtr; // Pointer to the adjacenccy list of each row/column + IndexVector m_innerIndices; // Adjacency list +}; + +}// end namespace eigen +#endif
\ No newline at end of file diff --git a/bench/spbench/CMakeLists.txt b/bench/spbench/CMakeLists.txt index a093cc5d9..2eb0befa9 100644 --- a/bench/spbench/CMakeLists.txt +++ b/bench/spbench/CMakeLists.txt @@ -66,5 +66,11 @@ target_link_libraries (spbenchsolver ${SPARSE_LIBS}) add_executable(spsolver sp_solver.cpp) target_link_libraries (spsolver ${SPARSE_LIBS}) +if(METIS_FOUND) + include_directories(${METIS_INCLUDES}) + set (SPARSE_LIBS ${SPARSE_LIBS} ${METIS_LIBRARIES}) + add_definitions("-DEIGEN_METIS_SUPPORT") +endif(METIS_FOUND) + add_executable(test_sparseLU test_sparseLU.cpp) target_link_libraries (test_sparseLU ${SPARSE_LIBS}) diff --git a/bench/spbench/test_sparseLU.cpp b/bench/spbench/test_sparseLU.cpp index 59f8252d0..8c78b0c9b 100644 --- a/bench/spbench/test_sparseLU.cpp +++ b/bench/spbench/test_sparseLU.cpp @@ -7,6 +7,9 @@ #include <unsupported/Eigen/SparseExtra> #include <Eigen/SparseLU> #include <bench/BenchTimer.h> +#ifdef EIGEN_METIS_SUPPORT +#include <Eigen/MetisSupport> +#endif using namespace std; using namespace Eigen; @@ -21,7 +24,12 @@ int main(int argc, char **args) typedef Matrix<scalar, Dynamic, 1> DenseRhs; Matrix<scalar, Dynamic, 1> b, x, tmp; // SparseLU<SparseMatrix<scalar, ColMajor>, AMDOrdering<int> > solver; +#ifdef EIGEN_METIS_SUPPORT + SparseLU<SparseMatrix<scalar, ColMajor>, MetisOrdering<int> > solver; +#else SparseLU<SparseMatrix<scalar, ColMajor>, COLAMDOrdering<int> > solver; +#endif + ifstream matrix_file; string line; int n; diff --git a/cmake/FindMetis.cmake b/cmake/FindMetis.cmake index e4d6ef258..627c3e9ae 100644 --- a/cmake/FindMetis.cmake +++ b/cmake/FindMetis.cmake @@ -12,10 +12,11 @@ find_path(METIS_INCLUDES ${INCLUDE_INSTALL_DIR} PATH_SUFFIXES metis + include ) -find_library(METIS_LIBRARIES metis PATHS $ENV{METISDIR} ${LIB_INSTALL_DIR}) +find_library(METIS_LIBRARIES metis PATHS $ENV{METISDIR} ${LIB_INSTALL_DIR} PATH_SUFFIXES lib) include(FindPackageHandleStandardArgs) find_package_handle_standard_args(METIS DEFAULT_MSG |