17faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 27faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This file is part of Eigen, a lightweight C++ template library 37faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// for linear algebra. 47faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 57faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> 67faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 77faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This Source Code Form is subject to the terms of the Mozilla 87faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Public License v. 2.0. If a copy of the MPL was not distributed 97faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#ifndef EIGEN_ORDERING_H 127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#define EIGEN_ORDERING_H 137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeznamespace Eigen { 157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#include "Eigen_Colamd.h" 177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeznamespace internal { 197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/** \internal 217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \ingroup OrderingMethods_Module 227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \returns the symmetric pattern A^T+A from the input matrix A. 237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * FIXME: The values should not be considered here 247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType> 267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid ordering_helper_at_plus_a(const MatrixType& mat, MatrixType& symmat) 277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixType C; 297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez C = mat.transpose(); // NOTE: Could be costly 307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (int i = 0; i < C.rows(); i++) 317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (typename MatrixType::InnerIterator it(C, i); it; ++it) 337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez it.valueRef() = 0.0; 347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez symmat = C + mat; 367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#ifndef EIGEN_MPL2_ONLY 417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/** \ingroup OrderingMethods_Module 437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \class AMDOrdering 447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Functor computing the \em approximate \em minimum \em degree ordering 467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * If the matrix is not structurally symmetric, an ordering of A^T+A is computed 477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \tparam Index The type of indices of the matrix 487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa COLAMDOrdering 497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate <typename Index> 517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezclass AMDOrdering 527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez public: 547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; 557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** Compute the permutation vector from a sparse matrix 577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This routine is much faster if the input matrix is column-major 587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template <typename MatrixType> 607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void operator()(const MatrixType& mat, PermutationType& perm) 617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Compute the symmetric pattern 637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SparseMatrix<typename MatrixType::Scalar, ColMajor, Index> symm; 647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::ordering_helper_at_plus_a(mat,symm); 657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Call the AMD routine 677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez //m_mat.prune(keep_diag()); 687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::minimum_degree_ordering(symm, perm); 697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** Compute the permutation with a selfadjoint matrix */ 727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template <typename SrcType, unsigned int SrcUpLo> 737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void operator()(const SparseSelfAdjointView<SrcType, SrcUpLo>& mat, PermutationType& perm) 747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SparseMatrix<typename SrcType::Scalar, ColMajor, Index> C; C = mat; 767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Call the AMD routine 787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // m_mat.prune(keep_diag()); //Remove the diagonal elements 797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::minimum_degree_ordering(C, perm); 807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#endif // EIGEN_MPL2_ONLY 847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/** \ingroup OrderingMethods_Module 867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \class NaturalOrdering 877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Functor computing the natural ordering (identity) 897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \note Returns an empty permutation matrix 917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \tparam Index The type of indices of the matrix 927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate <typename Index> 947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezclass NaturalOrdering 957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez public: 977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; 987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** Compute the permutation vector from a column-major sparse matrix */ 1007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template <typename MatrixType> 1017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void operator()(const MatrixType& /*mat*/, PermutationType& perm) 1027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez perm.resize(0); 1047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 1077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/** \ingroup OrderingMethods_Module 1097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \class COLAMDOrdering 1107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Functor computing the \em column \em approximate \em minimum \em degree ordering 1127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * The matrix should be in column-major and \b compressed format (see SparseMatrix::makeCompressed()). 1137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename Index> 1157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezclass COLAMDOrdering 1167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 1177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez public: 1187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; 1197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Index, Dynamic, 1> IndexVector; 1207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** Compute the permutation vector \a perm form the sparse matrix \a mat 1227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \warning The input sparse matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()). 1237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template <typename MatrixType> 1257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void operator() (const MatrixType& mat, PermutationType& perm) 1267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(mat.isCompressed() && "COLAMDOrdering requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to COLAMDOrdering"); 1287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index m = mat.rows(); 1307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index n = mat.cols(); 1317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index nnz = mat.nonZeros(); 1327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Get the recommended value of Alen to be used by colamd 1337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index Alen = internal::colamd_recommended(nnz, m, n); 1347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Set the default parameters 1357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double knobs [COLAMD_KNOBS]; 1367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index stats [COLAMD_STATS]; 1377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::colamd_set_defaults(knobs); 1387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez IndexVector p(n+1), A(Alen); 1407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(Index i=0; i <= n; i++) p(i) = mat.outerIndexPtr()[i]; 1417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(Index i=0; i < nnz; i++) A(i) = mat.innerIndexPtr()[i]; 1427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Call Colamd routine to compute the ordering 1437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index info = internal::colamd(m, n, Alen, A.data(), p.data(), knobs, stats); 1447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez EIGEN_UNUSED_VARIABLE(info); 1457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert( info && "COLAMD failed " ); 1467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez perm.resize(n); 1487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (Index i = 0; i < n; i++) perm.indices()(p(i)) = i; 1497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 1517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} // end namespace Eigen 1537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#endif 155