1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
|
// 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
|
