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>
|
---|
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
|
---|