[136] | 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
|
---|