GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-04-22 01:29:35 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add a draft (not clean ) version of the COLAMD ordering implementation

Browse Source
This commit is contained in:
Desire NUENTSA 2012-07-18 16:59:00 +02:00
parent 773804691a
commit b0cba2d988
5 changed files with 2604 additions and 56 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2515
Eigen/src/OrderingMethods/Eigen_Colamd.h
View File

File diff suppressed because it is too large Load Diff

93
Eigen/src/OrderingMethods/Ordering.h
View File

@ -27,6 +27,7 @@
#define EIGEN_ORDERING_H #define EIGEN_ORDERING_H
#include "Amd.h" #include "Amd.h"
#include "Eigen_Colamd.h"
namespace Eigen { namespace Eigen {
namespace internal { namespace internal {
@ -112,54 +113,50 @@ class NaturalOrdering
* Get the column approximate minimum degree ordering * Get the column approximate minimum degree ordering
* The matrix should be in column-major format * The matrix should be in column-major format
*/ */
// template<typename Scalar, typename Index> template<typename Index>
// class COLAMDOrdering: public OrderingBase< ColamdOrdering<Scalar, Index> > class COLAMDOrdering;
// { #include "Eigen_Colamd.h"
// public:
// typedef OrderingBase< ColamdOrdering<Scalar, Index> > Base; template<typename Index>
// typedef SparseMatrix<Scalar,ColMajor,Index> MatrixType; class COLAMDOrdering
// {
// public: public:
// COLAMDOrdering():Base() {} typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
// typedef Matrix<Index, Dynamic, 1> IndexVector;
// COLAMDOrdering(const MatrixType& matrix):Base() /** Compute the permutation vector form a sparse matrix */
// {
// compute(matrix);
// }
// COLAMDOrdering(const MatrixType& mat, PermutationType& perm_c):Base() template <typename MatrixType>
// { void operator() (const MatrixType& mat, PermutationType& perm)
// compute(matrix); {
// perm_c = this.get_perm(); int m = mat.rows();
// } int n = mat.cols();
// void compute(const MatrixType& mat) int nnz = mat.nonZeros();
// { // Get the recommended value of Alen to be used by colamd
// // Test if the matrix is column major... int Alen = eigen_colamd_recommended(nnz, m, n);
// // Set the default parameters
// int m = mat.rows(); double knobs [EIGEN_COLAMD_KNOBS];
// int n = mat.cols(); int stats [EIGEN_COLAMD_STATS];
// int nnz = mat.nonZeros(); eigen_colamd_set_defaults(knobs);
// // Get the recommended value of Alen to be used by colamd
// int Alen = colamd_recommended(nnz, m, n); int info;
// // Set the default parameters IndexVector p(n+1), A(Alen);
// double knobs[COLAMD_KNOBS]; for(int i=0; i <= n; i++) p(i) = mat.outerIndexPtr()[i];
// colamd_set_defaults(knobs); for(int i=0; i < nnz; i++) A(i) = mat.innerIndexPtr()[i];
// // Call Colamd routine to compute the ordering
// int info; info = eigen_colamd(m, n, Alen, A.data(), p.data(), knobs, stats);
// VectorXi p(n), A(nnz); eigen_assert( info && "COLAMD failed " );
// for(int i=0; i < n; i++) p(i) = mat.outerIndexPtr()(i);
// for(int i=0; i < nnz; i++) A(i) = mat.innerIndexPtr()(i); perm.resize(n);
// // Call Colamd routine to compute the ordering for (int i = 0; i < n; i++) perm.indices()(p(i)) = i;
// info = colamd(m, n, Alen, A,p , knobs, stats)
// eigen_assert( (info != FALSE)&& "COLAMD failed " ); }
//
// m_P.resize(n); private:
// for (int i = 0; i < n; i++) m_P(p(i)) = i;
// m_isInitialized = true;
// } };
// protected:
// using Base::m_isInitialized;
// using Base m_P;
// };
} // end namespace Eigen } // end namespace Eigen
#endif #endif

28
Eigen/src/SparseLU/SparseLU.h
View File

@ -99,8 +99,29 @@ class SparseLU
{ {
m_diagpivotthresh = thresh; m_diagpivotthresh = thresh;
} }
/** Return the number of nonzero elements in the L factor */
int nnzL()
{
if (m_factorizationIsOk)
return m_nnzL;
else
{
std::cerr<<"Numerical factorization should be done before\n";
return 0;
}
}
/** Return the number of nonzero elements in the U factor */
int nnzU()
{
if (m_factorizationIsOk)
return m_nnzU;
else
{
std::cerr<<"Numerical factorization should be done before\n";
return 0;
}
}
/** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A.
* *
* \sa compute() * \sa compute()
@ -325,7 +346,8 @@ void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat)
ord(mat,m_perm_c); ord(mat,m_perm_c);
//FIXME Check the right semantic behind m_perm_c //FIXME Check the right semantic behind m_perm_c
// that is, column j of mat goes to column m_perm_c(j) of mat * m_perm_c; // that is, column j of mat goes to column m_perm_c(j) of mat * m_perm_c;
// Apply the permutation to the column of the input matrix // Apply the permutation to the column of the input matrix
m_mat = mat * m_perm_c.inverse(); //FIXME It should be less expensive here to permute only the structural pattern of the matrix m_mat = mat * m_perm_c.inverse(); //FIXME It should be less expensive here to permute only the structural pattern of the matrix

2
Eigen/src/SuperLUSupport/SuperLUSupport.h
View File

@ -628,7 +628,7 @@ void SuperLU<MatrixType>::factorize(const MatrixType& a)
this->initFactorization(a); this->initFactorization(a);
//DEBUG //DEBUG
m_sluOptions.ColPerm = NATURAL; // m_sluOptions.ColPerm = COLAMD;
m_sluOptions.Equil = NO; m_sluOptions.Equil = NO;
int info = 0; int info = 0;
RealScalar recip_pivot_growth, rcond; RealScalar recip_pivot_growth, rcond;

22
bench/spbench/test_sparseLU.cpp
View File

@ -6,6 +6,7 @@
#include <iomanip> #include <iomanip>
#include <unsupported/Eigen/SparseExtra> #include <unsupported/Eigen/SparseExtra>
#include <Eigen/SparseLU> #include <Eigen/SparseLU>
#include <bench/BenchTimer.h>
using namespace std; using namespace std;
using namespace Eigen; using namespace Eigen;
@ -17,10 +18,12 @@ int main(int argc, char **args)
typedef Matrix<double, Dynamic, Dynamic> DenseMatrix; typedef Matrix<double, Dynamic, Dynamic> DenseMatrix;
typedef Matrix<double, Dynamic, 1> DenseRhs; typedef Matrix<double, Dynamic, 1> DenseRhs;
VectorXd b, x, tmp; VectorXd b, x, tmp;
SparseLU<SparseMatrix<double, ColMajor>, AMDOrdering<int> > solver; // SparseLU<SparseMatrix<double, ColMajor>, AMDOrdering<int> > solver;
SparseLU<SparseMatrix<double, ColMajor>, COLAMDOrdering<int> > solver;
ifstream matrix_file; ifstream matrix_file;
string line; string line;
int n; int n;
BenchTimer timer;
// Set parameters // Set parameters
/* Fill the matrix with sparse matrix stored in Matrix-Market coordinate column-oriented format */ /* Fill the matrix with sparse matrix stored in Matrix-Market coordinate column-oriented format */
@ -53,13 +56,26 @@ int main(int argc, char **args)
/* Compute the factorization */ /* Compute the factorization */
// solver.isSymmetric(true); // solver.isSymmetric(true);
solver.compute(A); timer.start();
// solver.compute(A);
solver.analyzePattern(A);
timer.stop();
cout << "Time to analyze " << timer.value() << std::endl;
timer.reset();
timer.start();
solver.factorize(A);
timer.stop();
cout << "Factorize Time " << timer.value() << std::endl;
timer.reset();
timer.start();
solver._solve(b, x); solver._solve(b, x);
timer.stop();
cout << "solve time " << timer.value() << std::endl;
/* Check the accuracy */ /* Check the accuracy */
VectorXd tmp2 = b - A*x; VectorXd tmp2 = b - A*x;
double tempNorm = tmp2.norm()/b.norm(); double tempNorm = tmp2.norm()/b.norm();
cout << "Relative norm of the computed solution : " << tempNorm <<"\n"; cout << "Relative norm of the computed solution : " << tempNorm <<"\n";
cout << "Number of nonzeros in the factor : " << solver.nnzL() + solver.nnzU() << std::endl;
return 0; return 0;
} }