![]() |
Eigen
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 #ifndef METIS_SUPPORT_H 00010 #define METIS_SUPPORT_H 00011 00012 namespace Eigen { 00021 template <typename StorageIndex> 00022 class MetisOrdering 00023 { 00024 public: 00025 typedef PermutationMatrix<Dynamic,Dynamic,StorageIndex> PermutationType; 00026 typedef Matrix<StorageIndex,Dynamic,1> IndexVector; 00027 00028 template <typename MatrixType> 00029 void get_symmetrized_graph(const MatrixType& A) 00030 { 00031 Index m = A.cols(); 00032 eigen_assert((A.rows() == A.cols()) && "ONLY FOR SQUARED MATRICES"); 00033 // Get the transpose of the input matrix 00034 MatrixType At = A.transpose(); 00035 // Get the number of nonzeros elements in each row/col of At+A 00036 Index TotNz = 0; 00037 IndexVector visited(m); 00038 visited.setConstant(-1); 00039 for (StorageIndex j = 0; j < m; j++) 00040 { 00041 // Compute the union structure of of A(j,:) and At(j,:) 00042 visited(j) = j; // Do not include the diagonal element 00043 // Get the nonzeros in row/column j of A 00044 for (typename MatrixType::InnerIterator it(A, j); it; ++it) 00045 { 00046 Index idx = it.index(); // Get the row index (for column major) or column index (for row major) 00047 if (visited(idx) != j ) 00048 { 00049 visited(idx) = j; 00050 ++TotNz; 00051 } 00052 } 00053 //Get the nonzeros in row/column j of At 00054 for (typename MatrixType::InnerIterator it(At, j); it; ++it) 00055 { 00056 Index idx = it.index(); 00057 if(visited(idx) != j) 00058 { 00059 visited(idx) = j; 00060 ++TotNz; 00061 } 00062 } 00063 } 00064 // Reserve place for A + At 00065 m_indexPtr.resize(m+1); 00066 m_innerIndices.resize(TotNz); 00067 00068 // Now compute the real adjacency list of each column/row 00069 visited.setConstant(-1); 00070 StorageIndex CurNz = 0; 00071 for (StorageIndex j = 0; j < m; j++) 00072 { 00073 m_indexPtr(j) = CurNz; 00074 00075 visited(j) = j; // Do not include the diagonal element 00076 // Add the pattern of row/column j of A to A+At 00077 for (typename MatrixType::InnerIterator it(A,j); it; ++it) 00078 { 00079 StorageIndex idx = it.index(); // Get the row index (for column major) or column index (for row major) 00080 if (visited(idx) != j ) 00081 { 00082 visited(idx) = j; 00083 m_innerIndices(CurNz) = idx; 00084 CurNz++; 00085 } 00086 } 00087 //Add the pattern of row/column j of At to A+At 00088 for (typename MatrixType::InnerIterator it(At, j); it; ++it) 00089 { 00090 StorageIndex idx = it.index(); 00091 if(visited(idx) != j) 00092 { 00093 visited(idx) = j; 00094 m_innerIndices(CurNz) = idx; 00095 ++CurNz; 00096 } 00097 } 00098 } 00099 m_indexPtr(m) = CurNz; 00100 } 00101 00102 template <typename MatrixType> 00103 void operator() (const MatrixType& A, PermutationType& matperm) 00104 { 00105 StorageIndex m = internal::convert_index<StorageIndex>(A.cols()); // must be StorageIndex, because it is passed by address to METIS 00106 IndexVector perm(m),iperm(m); 00107 // First, symmetrize the matrix graph. 00108 get_symmetrized_graph(A); 00109 int output_error; 00110 00111 // Call the fill-reducing routine from METIS 00112 output_error = METIS_NodeND(&m, m_indexPtr.data(), m_innerIndices.data(), NULL, NULL, perm.data(), iperm.data()); 00113 00114 if(output_error != METIS_OK) 00115 { 00116 //FIXME The ordering interface should define a class of possible errors 00117 std::cerr << "ERROR WHILE CALLING THE METIS PACKAGE \n"; 00118 return; 00119 } 00120 00121 // Get the fill-reducing permutation 00122 //NOTE: If Ap is the permuted matrix then perm and iperm vectors are defined as follows 00123 // 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 00124 00125 matperm.resize(m); 00126 for (int j = 0; j < m; j++) 00127 matperm.indices()(iperm(j)) = j; 00128 00129 } 00130 00131 protected: 00132 IndexVector m_indexPtr; // Pointer to the adjacenccy list of each row/column 00133 IndexVector m_innerIndices; // Adjacency list 00134 }; 00135 00136 }// end namespace eigen 00137 #endif