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
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>
ordering_helper_at_plus_a(const MatrixType & mat,MatrixType & symmat)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>
operator()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>
operator()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>
operator()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>
operator()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