1 2 // This file is part of Eigen, a lightweight C++ template library 3 // for linear algebra. 4 // 5 // Copyright (C) 2012 Dsir Nuentsa-Wakam <desire.nuentsa_wakam (at) inria.fr> 6 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 #ifndef EIGEN_ORDERING_H 12 #define EIGEN_ORDERING_H 13 14 namespace Eigen { 15 16 #include "Eigen_Colamd.h" 17 18 namespace internal { 19 20 /** \internal 21 * \ingroup OrderingMethods_Module 22 * \returns the symmetric pattern A^T+A from the input matrix A. 23 * FIXME: The values should not be considered here 24 */ 25 template<typename MatrixType> 26 void ordering_helper_at_plus_a(const MatrixType& mat, MatrixType& symmat) 27 { 28 MatrixType C; 29 C = mat.transpose(); // NOTE: Could be costly 30 for (int i = 0; i < C.rows(); i++) 31 { 32 for (typename MatrixType::InnerIterator it(C, i); it; ++it) 33 it.valueRef() = 0.0; 34 } 35 symmat = C + mat; 36 } 37 38 } 39 40 #ifndef EIGEN_MPL2_ONLY 41 42 /** \ingroup OrderingMethods_Module 43 * \class AMDOrdering 44 * 45 * Functor computing the \em approximate \em minimum \em degree ordering 46 * If the matrix is not structurally symmetric, an ordering of A^T+A is computed 47 * \tparam Index The type of indices of the matrix 48 * \sa COLAMDOrdering 49 */ 50 template <typename Index> 51 class AMDOrdering 52 { 53 public: 54 typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; 55 56 /** Compute the permutation vector from a sparse matrix 57 * This routine is much faster if the input matrix is column-major 58 */ 59 template <typename MatrixType> 60 void operator()(const MatrixType& mat, PermutationType& perm) 61 { 62 // Compute the symmetric pattern 63 SparseMatrix<typename MatrixType::Scalar, ColMajor, Index> symm; 64 internal::ordering_helper_at_plus_a(mat,symm); 65 66 // Call the AMD routine 67 //m_mat.prune(keep_diag()); 68 internal::minimum_degree_ordering(symm, perm); 69 } 70 71 /** Compute the permutation with a selfadjoint matrix */ 72 template <typename SrcType, unsigned int SrcUpLo> 73 void operator()(const SparseSelfAdjointView<SrcType, SrcUpLo>& mat, PermutationType& perm) 74 { 75 SparseMatrix<typename SrcType::Scalar, ColMajor, Index> C; C = mat; 76 77 // Call the AMD routine 78 // m_mat.prune(keep_diag()); //Remove the diagonal elements 79 internal::minimum_degree_ordering(C, perm); 80 } 81 }; 82 83 #endif // EIGEN_MPL2_ONLY 84 85 /** \ingroup OrderingMethods_Module 86 * \class NaturalOrdering 87 * 88 * Functor computing the natural ordering (identity) 89 * 90 * \note Returns an empty permutation matrix 91 * \tparam Index The type of indices of the matrix 92 */ 93 template <typename Index> 94 class NaturalOrdering 95 { 96 public: 97 typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; 98 99 /** Compute the permutation vector from a column-major sparse matrix */ 100 template <typename MatrixType> 101 void operator()(const MatrixType& /*mat*/, PermutationType& perm) 102 { 103 perm.resize(0); 104 } 105 106 }; 107 108 /** \ingroup OrderingMethods_Module 109 * \class COLAMDOrdering 110 * 111 * Functor computing the \em column \em approximate \em minimum \em degree ordering 112 * The matrix should be in column-major and \b compressed format (see SparseMatrix::makeCompressed()). 113 */ 114 template<typename Index> 115 class COLAMDOrdering 116 { 117 public: 118 typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; 119 typedef Matrix<Index, Dynamic, 1> IndexVector; 120 121 /** Compute the permutation vector \a perm form the sparse matrix \a mat 122 * \warning The input sparse matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()). 123 */ 124 template <typename MatrixType> 125 void operator() (const MatrixType& mat, PermutationType& perm) 126 { 127 eigen_assert(mat.isCompressed() && "COLAMDOrdering requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to COLAMDOrdering"); 128 129 Index m = mat.rows(); 130 Index n = mat.cols(); 131 Index nnz = mat.nonZeros(); 132 // Get the recommended value of Alen to be used by colamd 133 Index Alen = internal::colamd_recommended(nnz, m, n); 134 // Set the default parameters 135 double knobs [COLAMD_KNOBS]; 136 Index stats [COLAMD_STATS]; 137 internal::colamd_set_defaults(knobs); 138 139 IndexVector p(n+1), A(Alen); 140 for(Index i=0; i <= n; i++) p(i) = mat.outerIndexPtr()[i]; 141 for(Index i=0; i < nnz; i++) A(i) = mat.innerIndexPtr()[i]; 142 // Call Colamd routine to compute the ordering 143 Index info = internal::colamd(m, n, Alen, A.data(), p.data(), knobs, stats); 144 EIGEN_UNUSED_VARIABLE(info); 145 eigen_assert( info && "COLAMD failed " ); 146 147 perm.resize(n); 148 for (Index i = 0; i < n; i++) perm.indices()(p(i)) = i; 149 } 150 }; 151 152 } // end namespace Eigen 153 154 #endif 155
