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 
12 namespace 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  */
21 template <typename Index>
22 class MetisOrdering
23 {
24 public:
25   typedef PermutationMatrix<Dynamic,Dynamic,Index> PermutationType;
26   typedef Matrix<Index,Dynamic,1> IndexVector;
27 
28   template <typename MatrixType>
get_symmetrized_graph(const MatrixType & A)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>
operator()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