123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137 |
- // 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 StorageIndex>
- class MetisOrdering
- {
- public:
- typedef PermutationMatrix<Dynamic,Dynamic,StorageIndex> PermutationType;
- typedef Matrix<StorageIndex,Dynamic,1> IndexVector;
-
- template <typename MatrixType>
- void get_symmetrized_graph(const MatrixType& A)
- {
- Index m = A.cols();
- eigen_assert((A.rows() == A.cols()) && "ONLY FOR SQUARED MATRICES");
- // 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 (StorageIndex 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);
- StorageIndex CurNz = 0;
- for (StorageIndex 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)
- {
- StorageIndex 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)
- {
- StorageIndex 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)
- {
- StorageIndex m = internal::convert_index<StorageIndex>(A.cols()); // must be StorageIndex, because it is passed by address to METIS
- 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
-
- matperm.resize(m);
- for (int j = 0; j < m; j++)
- matperm.indices()(iperm(j)) = j;
-
- }
-
- protected:
- IndexVector m_indexPtr; // Pointer to the adjacenccy list of each row/column
- IndexVector m_innerIndices; // Adjacency list
- };
- }// end namespace eigen
- #endif
|