MetisSupport.h
Go to the documentation of this file.
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 Index>
00022 class MetisOrdering
00023 {
00024 public:
00025   typedef PermutationMatrix<Dynamic,Dynamic,Index> PermutationType;
00026   typedef Matrix<Index,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 (int 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     Index CurNz = 0; 
00071     for (int 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         Index 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         Index 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      Index m = A.cols();
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


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Thu Aug 27 2015 11:59:20