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>
|
---|
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
|
---|