Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.0
 
Loading...
Searching...
No Matches
MetisSupport.h
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 {
21template <typename StorageIndex>
23{
24public:
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 (StorageIndex 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 StorageIndex CurNz = 0;
71 for (StorageIndex 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 StorageIndex 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 StorageIndex 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 StorageIndex m = internal::convert_index<StorageIndex>(A.cols()); // must be StorageIndex, because it is passed by address to METIS
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
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
Definition: MetisSupport.h:23
void resize(Index newSize)
Definition: PermutationMatrix.h:125
Permutation matrix.
Definition: PermutationMatrix.h:298
const IndicesType & indices() const
Definition: PermutationMatrix.h:360
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:271
Derived & setConstant(Index size, const Scalar &val)
Definition: CwiseNullaryOp.h:361
const Scalar * data() const
Definition: PlainObjectBase.h:247
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74