1 2// This file is part of Eigen, a lightweight C++ template library 3// for linear algebra. 4// 5// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@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 14namespace Eigen { 15 16#include "Eigen_Colamd.h" 17 18namespace 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 */ 25template<typename MatrixType> 26void 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 */ 50template <typename Index> 51class 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 */ 93template <typename Index> 94class 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 */ 114template<typename Index> 115class 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