![]() |
Eigen
3.3.3
|
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