diff options
author | Sergey Sharybin <sergey.vfx@gmail.com> | 2014-03-21 14:04:53 +0400 |
---|---|---|
committer | Sergey Sharybin <sergey.vfx@gmail.com> | 2014-03-21 14:04:53 +0400 |
commit | 3411146984316c97f56543333a46f47aeb7f9d35 (patch) | |
tree | 5de608a3c18ff2a5459fd2191609dd882ad86213 /extern/Eigen3/Eigen/src/OrderingMethods/Ordering.h | |
parent | 1781928f9d720fa1dc4811afb69935b35aa89878 (diff) |
Update Eigen to version 3.2.1
Main purpose of this is to have SparseLU solver which
we can use now as a replacement to opennl library.
Diffstat (limited to 'extern/Eigen3/Eigen/src/OrderingMethods/Ordering.h')
-rw-r--r-- | extern/Eigen3/Eigen/src/OrderingMethods/Ordering.h | 150 |
1 files changed, 150 insertions, 0 deletions
diff --git a/extern/Eigen3/Eigen/src/OrderingMethods/Ordering.h b/extern/Eigen3/Eigen/src/OrderingMethods/Ordering.h new file mode 100644 index 00000000000..b4da6531a1d --- /dev/null +++ b/extern/Eigen3/Eigen/src/OrderingMethods/Ordering.h @@ -0,0 +1,150 @@ + +// 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 + * \returns the symmetric pattern A^T+A from the input matrix A. + * FIXME: The values should not be considered here + */ +template<typename MatrixType> +void ordering_helper_at_plus_a(const MatrixType& mat, MatrixType& symmat) +{ + MatrixType C; + C = mat.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 + mat; +} + +} + +#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 Index The type of indices of the matrix + * \sa COLAMDOrdering + */ +template <typename Index> +class AMDOrdering +{ + public: + typedef PermutationMatrix<Dynamic, Dynamic, Index> 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, Index> 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, Index> 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 Index The type of indices of the matrix + */ +template <typename Index> +class NaturalOrdering +{ + public: + typedef PermutationMatrix<Dynamic, Dynamic, Index> 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 + * + * Functor computing the \em column \em approximate \em minimum \em degree ordering + * The matrix should be in column-major format + */ +template<typename Index> +class COLAMDOrdering +{ + public: + typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; + typedef Matrix<Index, Dynamic, 1> IndexVector; + + /** Compute the permutation vector form a sparse matrix */ + template <typename MatrixType> + void operator() (const MatrixType& mat, PermutationType& perm) + { + Index m = mat.rows(); + Index n = mat.cols(); + Index nnz = mat.nonZeros(); + // Get the recommended value of Alen to be used by colamd + Index Alen = internal::colamd_recommended(nnz, m, n); + // Set the default parameters + double knobs [COLAMD_KNOBS]; + Index stats [COLAMD_STATS]; + internal::colamd_set_defaults(knobs); + + Index info; + IndexVector p(n+1), A(Alen); + for(Index i=0; i <= n; i++) p(i) = mat.outerIndexPtr()[i]; + for(Index i=0; i < nnz; i++) A(i) = mat.innerIndexPtr()[i]; + // Call Colamd routine to compute the ordering + info = internal::colamd(m, n, Alen, A.data(), p.data(), knobs, stats); + eigen_assert( info && "COLAMD failed " ); + + perm.resize(n); + for (Index i = 0; i < n; i++) perm.indices()(p(i)) = i; + } +}; + +} // end namespace Eigen + +#endif |