src/MetisSupport/MetisSupport.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009 #ifndef METIS_SUPPORT_H
0010 #define METIS_SUPPORT_H
0011
0012 namespace Eigen {
0013 /**
0014  * Get the fill-reducing ordering from the METIS package
0015  *
0016  * If A is the original matrix and Ap is the permuted matrix,
0017  * the fill-reducing permutation is defined as follows :
0018  * Row (column) i of A is the matperm(i) row (column) of Ap.
0019  * WARNING: As computed by METIS, this corresponds to the vector iperm (instead of perm)
0020  */
0021 template <typename StorageIndex>
0022 class MetisOrdering
0023 {
0024 public:
0025   typedef PermutationMatrix<Dynamic,Dynamic,StorageIndex> PermutationType;
0026   typedef Matrix<StorageIndex,Dynamic,1> IndexVector;
0027
0028   template <typename MatrixType>
0029   void get_symmetrized_graph(const MatrixType& A)
0030   {
0031     Index m = A.cols();
0032     eigen_assert((A.rows() == A.cols()) && "ONLY FOR SQUARED MATRICES");
0033     // Get the transpose of the input matrix
0034     MatrixType At = A.transpose();
0035     // Get the number of nonzeros elements in each row/col of At+A
0036     Index TotNz = 0;
0037     IndexVector visited(m);
0038     visited.setConstant(-1);
0039     for (StorageIndex j = 0; j < m; j++)
0040     {
0041       // Compute the union structure of of A(j,:) and At(j,:)
0042       visited(j) = j; // Do not include the diagonal element
0043       // Get the nonzeros in row/column j of A
0044       for (typename MatrixType::InnerIterator it(A, j); it; ++it)
0045       {
0046         Index idx = it.index(); // Get the row index (for column major) or column index (for row major)
0047         if (visited(idx) != j )
0048         {
0049           visited(idx) = j;
0050           ++TotNz;
0051         }
0052       }
0053       //Get the nonzeros in row/column j of At
0054       for (typename MatrixType::InnerIterator it(At, j); it; ++it)
0055       {
0056         Index idx = it.index();
0057         if(visited(idx) != j)
0058         {
0059           visited(idx) = j;
0060           ++TotNz;
0061         }
0062       }
0063     }
0064     // Reserve place for A + At
0065     m_indexPtr.resize(m+1);
0066     m_innerIndices.resize(TotNz);
0067
0068     // Now compute the real adjacency list of each column/row
0069     visited.setConstant(-1);
0070     StorageIndex CurNz = 0;
0071     for (StorageIndex j = 0; j < m; j++)
0072     {
0073       m_indexPtr(j) = CurNz;
0074
0075       visited(j) = j; // Do not include the diagonal element
0076       // Add the pattern of row/column j of A to A+At
0077       for (typename MatrixType::InnerIterator it(A,j); it; ++it)
0078       {
0079         StorageIndex idx = it.index(); // Get the row index (for column major) or column index (for row major)
0080         if (visited(idx) != j )
0081         {
0082           visited(idx) = j;
0083           m_innerIndices(CurNz) = idx;
0084           CurNz++;
0085         }
0086       }
0087       //Add the pattern of row/column j of At to A+At
0088       for (typename MatrixType::InnerIterator it(At, j); it; ++it)
0089       {
0090         StorageIndex idx = it.index();
0091         if(visited(idx) != j)
0092         {
0093           visited(idx) = j;
0094           m_innerIndices(CurNz) = idx;
0095           ++CurNz;
0096         }
0097       }
0098     }
0099     m_indexPtr(m) = CurNz;
0100   }
0101
0102   template <typename MatrixType>
0103   void operator() (const MatrixType& A, PermutationType& matperm)
0104   {
0105      StorageIndex m = internal::convert_index<StorageIndex>(A.cols()); // must be StorageIndex, because it is passed by address to METIS
0106      IndexVector perm(m),iperm(m);
0107     // First, symmetrize the matrix graph.
0108      get_symmetrized_graph(A);
0109      int output_error;
0110
0111      // Call the fill-reducing routine from METIS
0112      output_error = METIS_NodeND(&m, m_indexPtr.data(), m_innerIndices.data(), NULL, NULL, perm.data(), iperm.data());
0113
0114     if(output_error != METIS_OK)
0115     {
0116       //FIXME The ordering interface should define a class of possible errors
0117      std::cerr << "ERROR WHILE CALLING THE METIS PACKAGE \n";
0118      return;
0119     }
0120
0121     // Get the fill-reducing permutation
0122     //NOTE:  If Ap is the permuted matrix then perm and iperm vectors are defined as follows
0123     // 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
0124
0125      matperm.resize(m);
0126      for (int j = 0; j < m; j++)
0127        matperm.indices()(iperm(j)) = j;
0128
0129   }
0130
0131   protected:
0132     IndexVector m_indexPtr; // Pointer to the adjacenccy list of each row/column
0133     IndexVector m_innerIndices; // Adjacency list
0134 };
0135
0136 }// end namespace eigen
0137 #endif