MetisSupport.h
Go to the documentation of this file.
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 {
21 template <typename Index>
23 {
24 public:
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.
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
IndexVector m_indexPtr
Definition: MetisSupport.h:132
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
void get_symmetrized_graph(const MatrixType &A)
Definition: MetisSupport.h:29
EIGEN_STRONG_INLINE void resize(Index nbRows, Index nbCols)
PermutationMatrix< Dynamic, Dynamic, Index > PermutationType
Definition: MetisSupport.h:25
Derived & setConstant(Index size, const Scalar &value)
EIGEN_STRONG_INLINE const Scalar * data() const
Matrix< Index, Dynamic, 1 > IndexVector
Definition: MetisSupport.h:26
const IndicesType & indices() const
IndexVector m_innerIndices
Definition: MetisSupport.h:133
void operator()(const MatrixType &A, PermutationType &matperm)
Definition: MetisSupport.h:103
#define eigen_assert(x)
void resize(Index newSize)


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:34:54