From bccf64d34281066da48cf2da29fd61f7ed703025 Mon Sep 17 00:00:00 2001 From: Desire NUENTSA Date: Mon, 11 Jun 2012 18:52:26 +0200 Subject: Checking Syntax... --- Eigen/src/OrderingMethods/Ordering.h | 214 +++++++++++++++++++++++++++++++++++ 1 file changed, 214 insertions(+) create mode 100644 Eigen/src/OrderingMethods/Ordering.h (limited to 'Eigen/src/OrderingMethods/Ordering.h') diff --git a/Eigen/src/OrderingMethods/Ordering.h b/Eigen/src/OrderingMethods/Ordering.h new file mode 100644 index 000000000..c43c381a4 --- /dev/null +++ b/Eigen/src/OrderingMethods/Ordering.h @@ -0,0 +1,214 @@ + +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Désiré Nuentsa-Wakam +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_ORDERING_H +#define EIGEN_ORDERING_H + +#include +#include + +namespace Eigen { +template +class OrderingBase +{ + public: + typedef typename internal::traits::MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::Index Index; + typedef PermutationMatrix PermutationType; + + public: + OrderingBase():m_isInitialized(false) + { + + } + OrderingBase(const MatrixType& mat):OrderingBase() + { + compute(mat); + } + Derived& compute(const MatrixType& mat) + { + return derived().compute(mat); + } + Derived& derived() + { + return *static_cast(this); + } + const Derived& derived() const + { + return *static_cast(this); + } + /** + * Get the permutation vector + */ + PermutationType& get_perm(const MatrixType& mat) + { + if (m_isInitialized = true) return m_P; + else abort(); // FIXME Should find a smoother way to exit with error code + } + template + void at_plus_a(const MatrixType& mat); + + /** keeps off-diagonal entries; drops diagonal entries */ + struct keep_diag { + inline bool operator() (const Index& row, const Index& col, const Scalar&) const + { + return row!=col; + } + }; + + protected: + void init() + { + m_isInitialized = false; + } + PermutationType m_P; // The computed permutation + mutable bool m_isInitialized; + SparseMatrix m_mat; // Stores the (symmetrized) matrix to permute +} +/** + * Get the symmetric pattern A^T+A from the input matrix A. + * NOTE: The values should not be considered here + */ +template +void OrderingBase::at_plus_a(const MatrixType& mat) +{ + 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; + } + m_mat = C + mat; + +/** + * Get the column approximate minimum degree ordering + * The matrix should be in column-major format + */ +template +class COLAMDOrdering: public OrderingBase< ColamdOrdering > +{ + public: + typedef OrderingBase< ColamdOrdering > Base; + typedef SparseMatrix MatrixType; + + public: + COLAMDOrdering():Base() {} + + COLAMDOrdering(const MatrixType& matrix):Base() + { + compute(matrix); + } + COLAMDOrdering(const MatrixType& mat, PermutationType& perm_c):Base() + { + compute(matrix); + perm_c = this.get_perm(); + } + void compute(const MatrixType& mat) + { + // Test if the matrix is column major... + + int m = mat.rows(); + int n = mat.cols(); + int nnz = mat.nonZeros(); + // Get the recommended value of Alen to be used by colamd + int Alen = colamd_recommended(nnz, m, n); + // Set the default parameters + double knobs[COLAMD_KNOBS]; + colamd_set_defaults(knobs); + + int info; + VectorXi p(n), A(nnz); + for(int i=0; i < n; i++) p(i) = mat.outerIndexPtr()(i); + for(int i=0; i < nnz; i++) A(i) = mat.innerIndexPtr()(i); + // Call Colamd routine to compute the ordering + info = colamd(m, n, Alen, A,p , knobs, stats) + eigen_assert( (info != FALSE)&& "COLAMD failed " ); + + m_P.resize(n); + for (int i = 0; i < n; i++) m_P(p(i)) = i; + m_isInitialized = true; + } + protected: + using Base::m_isInitialized; + using Base m_P; +} + +/** + * Get the approximate minimum degree ordering + * If the matrix is not structurally symmetric, an ordering of A^T+A is computed + * \tparam Scalar The type of the scalar of the matrix for which the ordering is applied + * \tparam Index The type of indices of the matrix + * \tparam _UpLo If the matrix is symmetric, indicates which part to use + */ +template +class AMDordering : public OrderingBase > +{ + public: + enum { UpLo = _UpLo }; + typedef OrderingBase< AMDOrdering > Base; + typedef SparseMatrix MatrixType; + public: + AMDOrdering():Base(){} + AMDOrdering(const MatrixType& mat):Base() + { + compute(matrix); + } + AMDOrdering(const MatrixType& mat, PermutationType& perm_c):Base() + { + compute(matrix); + perm_c = this.get_perm(); + } + /** Compute the permutation vector from a column-major sparse matrix */ + void compute(const MatrixType& mat) + { + // Compute the symmetric pattern + at_plus_a(mat); + + // Call the AMD routine + m_mat.prune(keep_diag()); + internal::minimum_degree_ordering(m_mat, m_P); + if (m_P.size()>0) m_isInitialized = true; + } + /** Compute the permutation with a self adjoint matrix */ + template + void compute(const SparseSelfAdjointView& mat) + { + m_mat = mat; + + // Call the AMD routine + m_mat.prune(keep_diag()); + internal::minimum_degree_ordering(m_mat, m_P); + if (m_P.size()>0) m_isInitialized = true; + } + protected: + using Base::m_isInitialized; + using Base::m_P; + using Base::m_mat; +} + +} // end namespace Eigen +#endif \ No newline at end of file -- cgit v1.2.3