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 { 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> get_symmetrized_graph(const MatrixType & A)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> operator()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