src/OrderingMethods/Ordering.h

0001
0002 // This file is part of Eigen, a lightweight C++ template library
0003 // for linear algebra.
0004 //
0005 // Copyright (C) 2012  Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
0006 //
0007 // This Source Code Form is subject to the terms of the Mozilla
0008 // Public License v. 2.0. If a copy of the MPL was not distributed
0009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0010
0011 #ifndef EIGEN_ORDERING_H
0012 #define EIGEN_ORDERING_H
0013
0014 namespace Eigen {
0015
0016 #include "Eigen_Colamd.h"
0017
0018 namespace internal {
0019
0020 /** \internal
0021   * \ingroup OrderingMethods_Module
0022   * \param[in] A the input non-symmetric matrix
0023   * \param[out] symmat the symmetric pattern A^T+A from the input matrix \a A.
0024   * FIXME: The values should not be considered here
0025   */
0026 template<typename MatrixType>
0027 void ordering_helper_at_plus_a(const MatrixType& A, MatrixType& symmat)
0028 {
0029   MatrixType C;
0030   C = A.transpose(); // NOTE: Could be  costly
0031   for (int i = 0; i < C.rows(); i++)
0032   {
0033       for (typename MatrixType::InnerIterator it(C, i); it; ++it)
0034         it.valueRef() = typename MatrixType::Scalar(0);
0035   }
0036   symmat = C + A;
0037 }
0038
0039 }
0040
0041 /** \ingroup OrderingMethods_Module
0042   * \class AMDOrdering
0043   *
0044   * Functor computing the \em approximate \em minimum \em degree ordering
0045   * If the matrix is not structurally symmetric, an ordering of A^T+A is computed
0046   * \tparam  StorageIndex The type of indices of the matrix
0047   * \sa COLAMDOrdering
0048   */
0049 template <typename StorageIndex>
0050 class AMDOrdering
0051 {
0052   public:
0053     typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType;
0054
0055     /** Compute the permutation vector from a sparse matrix
0056      * This routine is much faster if the input matrix is column-major
0057      */
0058     template <typename MatrixType>
0059     void operator()(const MatrixType& mat, PermutationType& perm)
0060     {
0061       // Compute the symmetric pattern
0062       SparseMatrix<typename MatrixType::Scalar, ColMajor, StorageIndex> symm;
0063       internal::ordering_helper_at_plus_a(mat,symm);
0064
0065       // Call the AMD routine
0066       //m_mat.prune(keep_diag());
0067       internal::minimum_degree_ordering(symm, perm);
0068     }
0069
0070     /** Compute the permutation with a selfadjoint matrix */
0071     template <typename SrcType, unsigned int SrcUpLo>
0072     void operator()(const SparseSelfAdjointView<SrcType, SrcUpLo>& mat, PermutationType& perm)
0073     {
0074       SparseMatrix<typename SrcType::Scalar, ColMajor, StorageIndex> C; C = mat;
0075
0076       // Call the AMD routine
0077       // m_mat.prune(keep_diag()); //Remove the diagonal elements
0078       internal::minimum_degree_ordering(C, perm);
0079     }
0080 };
0081
0082 /** \ingroup OrderingMethods_Module
0083   * \class NaturalOrdering
0084   *
0085   * Functor computing the natural ordering (identity)
0086   *
0087   * \note Returns an empty permutation matrix
0088   * \tparam  StorageIndex The type of indices of the matrix
0089   */
0090 template <typename StorageIndex>
0091 class NaturalOrdering
0092 {
0093   public:
0094     typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType;
0095
0096     /** Compute the permutation vector from a column-major sparse matrix */
0097     template <typename MatrixType>
0098     void operator()(const MatrixType& /*mat*/, PermutationType& perm)
0099     {
0100       perm.resize(0);
0101     }
0102
0103 };
0104
0105 /** \ingroup OrderingMethods_Module
0106   * \class COLAMDOrdering
0107   *
0108   * \tparam  StorageIndex The type of indices of the matrix
0109   *
0110   * Functor computing the \em column \em approximate \em minimum \em degree ordering
0111   * The matrix should be in column-major and \b compressed format (see SparseMatrix::makeCompressed()).
0112   */
0113 template<typename StorageIndex>
0114 class COLAMDOrdering
0115 {
0116   public:
0117     typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType;
0118     typedef Matrix<StorageIndex, Dynamic, 1> IndexVector;
0119
0120     /** Compute the permutation vector \a perm form the sparse matrix \a mat
0121       * \warning The input sparse matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()).
0122       */
0123     template <typename MatrixType>
0124     void operator() (const MatrixType& mat, PermutationType& perm)
0125     {
0126       eigen_assert(mat.isCompressed() && "COLAMDOrdering requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to COLAMDOrdering");
0127
0128       StorageIndex m = StorageIndex(mat.rows());
0129       StorageIndex n = StorageIndex(mat.cols());
0130       StorageIndex nnz = StorageIndex(mat.nonZeros());
0131       // Get the recommended value of Alen to be used by colamd
0132       StorageIndex Alen = internal::Colamd::recommended(nnz, m, n);
0133       // Set the default parameters
0134       double knobs [internal::Colamd::NKnobs];
0135       StorageIndex stats [internal::Colamd::NStats];
0136       internal::Colamd::set_defaults(knobs);
0137
0138       IndexVector p(n+1), A(Alen);
0139       for(StorageIndex i=0; i <= n; i++)   p(i) = mat.outerIndexPtr()[i];
0140       for(StorageIndex i=0; i < nnz; i++)  A(i) = mat.innerIndexPtr()[i];
0141       // Call Colamd routine to compute the ordering
0142       StorageIndex info = internal::Colamd::compute_ordering(m, n, Alen, A.data(), p.data(), knobs, stats);
0143       EIGEN_UNUSED_VARIABLE(info);
0144       eigen_assert( info && "COLAMD failed " );
0145
0146       perm.resize(n);
0147       for (StorageIndex i = 0; i < n; i++) perm.indices()(p(i)) = i;
0148     }
0149 };
0150
0151 } // end namespace Eigen
0152
0153 #endif