Eigen  3.3.3
Ordering.h
00001  
00002 // This file is part of Eigen, a lightweight C++ template library
00003 // for linear algebra.
00004 //
00005 // Copyright (C) 2012  Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
00006 //
00007 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_ORDERING_H
00012 #define EIGEN_ORDERING_H
00013 
00014 namespace Eigen {
00015   
00016 #include "Eigen_Colamd.h"
00017 
00018 namespace internal {
00019     
00026 template<typename MatrixType> 
00027 void ordering_helper_at_plus_a(const MatrixType& A, MatrixType& symmat)
00028 {
00029   MatrixType C;
00030   C = A.transpose(); // NOTE: Could be  costly
00031   for (int i = 0; i < C.rows(); i++) 
00032   {
00033       for (typename MatrixType::InnerIterator it(C, i); it; ++it)
00034         it.valueRef() = 0.0;
00035   }
00036   symmat = C + A;
00037 }
00038     
00039 }
00040 
00041 #ifndef EIGEN_MPL2_ONLY
00042 
00051 template <typename StorageIndex>
00052 class AMDOrdering
00053 {
00054   public:
00055     typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType;
00056     
00060     template <typename MatrixType>
00061     void operator()(const MatrixType& mat, PermutationType& perm)
00062     {
00063       // Compute the symmetric pattern
00064       SparseMatrix<typename MatrixType::Scalar, ColMajor, StorageIndex> symm;
00065       internal::ordering_helper_at_plus_a(mat,symm); 
00066     
00067       // Call the AMD routine 
00068       //m_mat.prune(keep_diag());
00069       internal::minimum_degree_ordering(symm, perm);
00070     }
00071     
00073     template <typename SrcType, unsigned int SrcUpLo> 
00074     void operator()(const SparseSelfAdjointView<SrcType, SrcUpLo>& mat, PermutationType& perm)
00075     { 
00076       SparseMatrix<typename SrcType::Scalar, ColMajor, StorageIndex> C; C = mat;
00077       
00078       // Call the AMD routine 
00079       // m_mat.prune(keep_diag()); //Remove the diagonal elements 
00080       internal::minimum_degree_ordering(C, perm);
00081     }
00082 };
00083 
00084 #endif // EIGEN_MPL2_ONLY
00085 
00094 template <typename StorageIndex>
00095 class NaturalOrdering
00096 {
00097   public:
00098     typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType;
00099     
00101     template <typename MatrixType>
00102     void operator()(const MatrixType& /*mat*/, PermutationType& perm)
00103     {
00104       perm.resize(0); 
00105     }
00106     
00107 };
00108 
00117 template<typename StorageIndex>
00118 class COLAMDOrdering
00119 {
00120   public:
00121     typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType; 
00122     typedef Matrix<StorageIndex, Dynamic, 1> IndexVector;
00123     
00127     template <typename MatrixType>
00128     void operator() (const MatrixType& mat, PermutationType& perm)
00129     {
00130       eigen_assert(mat.isCompressed() && "COLAMDOrdering requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to COLAMDOrdering");
00131       
00132       StorageIndex m = StorageIndex(mat.rows());
00133       StorageIndex n = StorageIndex(mat.cols());
00134       StorageIndex nnz = StorageIndex(mat.nonZeros());
00135       // Get the recommended value of Alen to be used by colamd
00136       StorageIndex Alen = internal::colamd_recommended(nnz, m, n); 
00137       // Set the default parameters
00138       double knobs [COLAMD_KNOBS]; 
00139       StorageIndex stats [COLAMD_STATS];
00140       internal::colamd_set_defaults(knobs);
00141       
00142       IndexVector p(n+1), A(Alen); 
00143       for(StorageIndex i=0; i <= n; i++)   p(i) = mat.outerIndexPtr()[i];
00144       for(StorageIndex i=0; i < nnz; i++)  A(i) = mat.innerIndexPtr()[i];
00145       // Call Colamd routine to compute the ordering 
00146       StorageIndex info = internal::colamd(m, n, Alen, A.data(), p.data(), knobs, stats); 
00147       EIGEN_UNUSED_VARIABLE(info);
00148       eigen_assert( info && "COLAMD failed " );
00149       
00150       perm.resize(n);
00151       for (StorageIndex i = 0; i < n; i++) perm.indices()(p(i)) = i;
00152     }
00153 };
00154 
00155 } // end namespace Eigen
00156 
00157 #endif
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends