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