1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
5//
6// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9#ifndef METIS_SUPPORT_H
10#define METIS_SUPPORT_H
11
12namespace Eigen {
13/**
14 * Get the fill-reducing ordering from the METIS package
15 *
16 * If A is the original matrix and Ap is the permuted matrix,
17 * the fill-reducing permutation is defined as follows :
18 * Row (column) i of A is the matperm(i) row (column) of Ap.
19 * WARNING: As computed by METIS, this corresponds to the vector iperm (instead of perm)
20 */
21template <typename Index>
22class MetisOrdering
23{
24public:
25  typedef PermutationMatrix<Dynamic,Dynamic,Index> PermutationType;
26  typedef Matrix<Index,Dynamic,1> IndexVector;
27
28  template <typename MatrixType>
29  void get_symmetrized_graph(const MatrixType& A)
30  {
31    Index m = A.cols();
32    eigen_assert((A.rows() == A.cols()) && "ONLY FOR SQUARED MATRICES");
33    // Get the transpose of the input matrix
34    MatrixType At = A.transpose();
35    // Get the number of nonzeros elements in each row/col of At+A
36    Index TotNz = 0;
37    IndexVector visited(m);
38    visited.setConstant(-1);
39    for (int j = 0; j < m; j++)
40    {
41      // Compute the union structure of of A(j,:) and At(j,:)
42      visited(j) = j; // Do not include the diagonal element
43      // Get the nonzeros in row/column j of A
44      for (typename MatrixType::InnerIterator it(A, j); it; ++it)
45      {
46        Index idx = it.index(); // Get the row index (for column major) or column index (for row major)
47        if (visited(idx) != j )
48        {
49          visited(idx) = j;
50          ++TotNz;
51        }
52      }
53      //Get the nonzeros in row/column j of At
54      for (typename MatrixType::InnerIterator it(At, j); it; ++it)
55      {
56        Index idx = it.index();
57        if(visited(idx) != j)
58        {
59          visited(idx) = j;
60          ++TotNz;
61        }
62      }
63    }
64    // Reserve place for A + At
65    m_indexPtr.resize(m+1);
66    m_innerIndices.resize(TotNz);
67
68    // Now compute the real adjacency list of each column/row
69    visited.setConstant(-1);
70    Index CurNz = 0;
71    for (int j = 0; j < m; j++)
72    {
73      m_indexPtr(j) = CurNz;
74
75      visited(j) = j; // Do not include the diagonal element
76      // Add the pattern of row/column j of A to A+At
77      for (typename MatrixType::InnerIterator it(A,j); it; ++it)
78      {
79        Index idx = it.index(); // Get the row index (for column major) or column index (for row major)
80        if (visited(idx) != j )
81        {
82          visited(idx) = j;
83          m_innerIndices(CurNz) = idx;
84          CurNz++;
85        }
86      }
87      //Add the pattern of row/column j of At to A+At
88      for (typename MatrixType::InnerIterator it(At, j); it; ++it)
89      {
90        Index idx = it.index();
91        if(visited(idx) != j)
92        {
93          visited(idx) = j;
94          m_innerIndices(CurNz) = idx;
95          ++CurNz;
96        }
97      }
98    }
99    m_indexPtr(m) = CurNz;
100  }
101
102  template <typename MatrixType>
103  void operator() (const MatrixType& A, PermutationType& matperm)
104  {
105     Index m = A.cols();
106     IndexVector perm(m),iperm(m);
107    // First, symmetrize the matrix graph.
108     get_symmetrized_graph(A);
109     int output_error;
110
111     // Call the fill-reducing routine from METIS
112     output_error = METIS_NodeND(&m, m_indexPtr.data(), m_innerIndices.data(), NULL, NULL, perm.data(), iperm.data());
113
114    if(output_error != METIS_OK)
115    {
116      //FIXME The ordering interface should define a class of possible errors
117     std::cerr << "ERROR WHILE CALLING THE METIS PACKAGE \n";
118     return;
119    }
120
121    // Get the fill-reducing permutation
122    //NOTE:  If Ap is the permuted matrix then perm and iperm vectors are defined as follows
123    // 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
124
125     matperm.resize(m);
126     for (int j = 0; j < m; j++)
127       matperm.indices()(iperm(j)) = j;
128
129  }
130
131  protected:
132    IndexVector m_indexPtr; // Pointer to the adjacenccy list of each row/column
133    IndexVector m_innerIndices; // Adjacency list
134};
135
136}// end namespace eigen
137#endif
138