1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
|
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_ORDERING_H
#define EIGEN_ORDERING_H
namespace Eigen {
#include "Eigen_Colamd.h"
namespace internal {
/** \internal
* \ingroup OrderingMethods_Module
* \param[in] A the input non-symmetric matrix
* \param[out] symmat the symmetric pattern A^T+A from the input matrix \a A.
* FIXME: The values should not be considered here
*/
template<typename MatrixType>
void ordering_helper_at_plus_a(const MatrixType& A, MatrixType& symmat)
{
MatrixType C;
C = A.transpose(); // NOTE: Could be costly
for (int i = 0; i < C.rows(); i++)
{
for (typename MatrixType::InnerIterator it(C, i); it; ++it)
it.valueRef() = 0.0;
}
symmat = C + A;
}
}
#ifndef EIGEN_MPL2_ONLY
/** \ingroup OrderingMethods_Module
* \class AMDOrdering
*
* Functor computing the \em approximate \em minimum \em degree ordering
* If the matrix is not structurally symmetric, an ordering of A^T+A is computed
* \tparam StorageIndex The type of indices of the matrix
* \sa COLAMDOrdering
*/
template <typename StorageIndex>
class AMDOrdering
{
public:
typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType;
/** Compute the permutation vector from a sparse matrix
* This routine is much faster if the input matrix is column-major
*/
template <typename MatrixType>
void operator()(const MatrixType& mat, PermutationType& perm)
{
// Compute the symmetric pattern
SparseMatrix<typename MatrixType::Scalar, ColMajor, StorageIndex> symm;
internal::ordering_helper_at_plus_a(mat,symm);
// Call the AMD routine
//m_mat.prune(keep_diag());
internal::minimum_degree_ordering(symm, perm);
}
/** Compute the permutation with a selfadjoint matrix */
template <typename SrcType, unsigned int SrcUpLo>
void operator()(const SparseSelfAdjointView<SrcType, SrcUpLo>& mat, PermutationType& perm)
{
SparseMatrix<typename SrcType::Scalar, ColMajor, StorageIndex> C; C = mat;
// Call the AMD routine
// m_mat.prune(keep_diag()); //Remove the diagonal elements
internal::minimum_degree_ordering(C, perm);
}
};
#endif // EIGEN_MPL2_ONLY
/** \ingroup OrderingMethods_Module
* \class NaturalOrdering
*
* Functor computing the natural ordering (identity)
*
* \note Returns an empty permutation matrix
* \tparam StorageIndex The type of indices of the matrix
*/
template <typename StorageIndex>
class NaturalOrdering
{
public:
typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType;
/** Compute the permutation vector from a column-major sparse matrix */
template <typename MatrixType>
void operator()(const MatrixType& /*mat*/, PermutationType& perm)
{
perm.resize(0);
}
};
/** \ingroup OrderingMethods_Module
* \class COLAMDOrdering
*
* \tparam StorageIndex The type of indices of the matrix
*
* Functor computing the \em column \em approximate \em minimum \em degree ordering
* The matrix should be in column-major and \b compressed format (see SparseMatrix::makeCompressed()).
*/
template<typename StorageIndex>
class COLAMDOrdering
{
public:
typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType;
typedef Matrix<StorageIndex, Dynamic, 1> IndexVector;
/** Compute the permutation vector \a perm form the sparse matrix \a mat
* \warning The input sparse matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()).
*/
template <typename MatrixType>
void operator() (const MatrixType& mat, PermutationType& perm)
{
eigen_assert(mat.isCompressed() && "COLAMDOrdering requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to COLAMDOrdering");
StorageIndex m = StorageIndex(mat.rows());
StorageIndex n = StorageIndex(mat.cols());
StorageIndex nnz = StorageIndex(mat.nonZeros());
// Get the recommended value of Alen to be used by colamd
StorageIndex Alen = internal::colamd_recommended(nnz, m, n);
// Set the default parameters
double knobs [COLAMD_KNOBS];
StorageIndex stats [COLAMD_STATS];
internal::colamd_set_defaults(knobs);
IndexVector p(n+1), A(Alen);
for(StorageIndex i=0; i <= n; i++) p(i) = mat.outerIndexPtr()[i];
for(StorageIndex i=0; i < nnz; i++) A(i) = mat.innerIndexPtr()[i];
// Call Colamd routine to compute the ordering
StorageIndex info = internal::colamd(m, n, Alen, A.data(), p.data(), knobs, stats);
EIGEN_UNUSED_VARIABLE(info);
eigen_assert( info && "COLAMD failed " );
perm.resize(n);
for (StorageIndex i = 0; i < n; i++) perm.indices()(p(i)) = i;
}
};
} // end namespace Eigen
#endif
|
