Home | History | Annotate | Download | only in MetisSupport 
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2012 Dsir Nuentsa-Wakam <desire.nuentsa_wakam (at) 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 StorageIndex>
     22 class MetisOrdering
     23 {
     24 public:
     25   typedef PermutationMatrix<Dynamic,Dynamic,StorageIndex> PermutationType;
     26   typedef Matrix<StorageIndex,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 (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
    138