diff options
author | Desire NUENTSA <desire.nuentsa_wakam@inria.fr> | 2012-09-07 13:18:16 +0200 |
---|---|---|
committer | Desire NUENTSA <desire.nuentsa_wakam@inria.fr> | 2012-09-07 13:18:16 +0200 |
commit | fdd0f0c5fc807f3b90e0442cfe4f3fa1f624b50a (patch) | |
tree | 63a7f166f41a90c9048b1d743fc6fb5b4b47a54a /bench/spbench | |
parent | 288e6aab14cc12e604bd1a12f0cba20d88edf54f (diff) | |
parent | 063705b5be5a41e324773887d3d5ae065321a719 (diff) |
merge Sparse LU branch
Diffstat (limited to 'bench/spbench')
-rw-r--r-- | bench/spbench/CMakeLists.txt | 12 | ||||
-rw-r--r-- | bench/spbench/sp_solver.cpp | 124 | ||||
-rw-r--r-- | bench/spbench/test_sparseLU.cpp | 93 |
3 files changed, 229 insertions, 0 deletions
diff --git a/bench/spbench/CMakeLists.txt b/bench/spbench/CMakeLists.txt index 079912266..5451843b9 100644 --- a/bench/spbench/CMakeLists.txt +++ b/bench/spbench/CMakeLists.txt @@ -63,3 +63,15 @@ endif(RT_LIBRARY) add_executable(spbenchsolver spbenchsolver.cpp) target_link_libraries (spbenchsolver ${SPARSE_LIBS}) +add_executable(spsolver sp_solver.cpp) +target_link_libraries (spsolver ${SPARSE_LIBS}) + +if(METIS_FOUND) + include_directories(${METIS_INCLUDES}) + set (SPARSE_LIBS ${SPARSE_LIBS} ${METIS_LIBRARIES}) + add_definitions("-DEIGEN_METIS_SUPPORT") +endif(METIS_FOUND) + +add_executable(test_sparseLU test_sparseLU.cpp) +target_link_libraries (test_sparseLU ${SPARSE_LIBS}) + diff --git a/bench/spbench/sp_solver.cpp b/bench/spbench/sp_solver.cpp new file mode 100644 index 000000000..e18f2d1c3 --- /dev/null +++ b/bench/spbench/sp_solver.cpp @@ -0,0 +1,124 @@ +// Small bench routine for Eigen available in Eigen +// (C) Desire NUENTSA WAKAM, INRIA + +#include <iostream> +#include <fstream> +#include <iomanip> +#include <Eigen/Jacobi> +#include <Eigen/Householder> +#include <Eigen/IterativeLinearSolvers> +#include <Eigen/LU> +#include <unsupported/Eigen/SparseExtra> +//#include <Eigen/SparseLU> +#include <Eigen/SuperLUSupport> +// #include <unsupported/Eigen/src/IterativeSolvers/Scaling.h> +#include <bench/BenchTimer.h> + +using namespace std; +using namespace Eigen; + +int main(int argc, char **args) +{ + SparseMatrix<double, ColMajor> A; + typedef SparseMatrix<double, ColMajor>::Index Index; + typedef Matrix<double, Dynamic, Dynamic> DenseMatrix; + typedef Matrix<double, Dynamic, 1> DenseRhs; + VectorXd b, x, tmp; + BenchTimer timer,totaltime; + //SparseLU<SparseMatrix<double, ColMajor> > solver; + SuperLU<SparseMatrix<double, ColMajor> > solver; + ifstream matrix_file; + string line; + int n; + // Set parameters +// solver.iparm(IPARM_THREAD_NBR) = 4; + /* Fill the matrix with sparse matrix stored in Matrix-Market coordinate column-oriented format */ + if (argc < 2) assert(false && "please, give the matrix market file "); + + timer.start(); + totaltime.start(); + loadMarket(A, args[1]); + cout << "End charging matrix " << endl; + bool iscomplex=false, isvector=false; + int sym; + getMarketHeader(args[1], sym, iscomplex, isvector); + if (iscomplex) { cout<< " Not for complex matrices \n"; return -1; } + if (isvector) { cout << "The provided file is not a matrix file\n"; return -1;} + if (sym != 0) { // symmetric matrices, only the lower part is stored + SparseMatrix<double, ColMajor> temp; + temp = A; + A = temp.selfadjointView<Lower>(); + } + timer.stop(); + + n = A.cols(); + // ====== TESTS FOR SPARSE TUTORIAL ====== +// cout<< "OuterSize " << A.outerSize() << " inner " << A.innerSize() << endl; +// SparseMatrix<double, RowMajor> mat1(A); +// SparseMatrix<double, RowMajor> mat2; +// cout << " norm of A " << mat1.norm() << endl; ; +// PermutationMatrix<Dynamic, Dynamic, int> perm(n); +// perm.resize(n,1); +// perm.indices().setLinSpaced(n, 0, n-1); +// mat2 = perm * mat1; +// mat.subrows(); +// mat2.resize(n,n); +// mat2.reserve(10); +// mat2.setConstant(); +// std::cout<< "NORM " << mat1.squaredNorm()<< endl; + + cout<< "Time to load the matrix " << timer.value() <<endl; + /* Fill the right hand side */ + +// solver.set_restart(374); + if (argc > 2) + loadMarketVector(b, args[2]); + else + { + b.resize(n); + tmp.resize(n); +// tmp.setRandom(); + for (int i = 0; i < n; i++) tmp(i) = i; + b = A * tmp ; + } +// Scaling<SparseMatrix<double> > scal; +// scal.computeRef(A); +// b = scal.LeftScaling().cwiseProduct(b); + + /* Compute the factorization */ + cout<< "Starting the factorization "<< endl; + timer.reset(); + timer.start(); + cout<< "Size of Input Matrix "<< b.size()<<"\n\n"; + cout<< "Rows and columns "<< A.rows() <<" " <<A.cols() <<"\n"; + solver.compute(A); +// solver.analyzePattern(A); +// solver.factorize(A); + if (solver.info() != Success) { + std::cout<< "The solver failed \n"; + return -1; + } + timer.stop(); + float time_comp = timer.value(); + cout <<" Compute Time " << time_comp<< endl; + + timer.reset(); + timer.start(); + x = solver.solve(b); +// x = scal.RightScaling().cwiseProduct(x); + timer.stop(); + float time_solve = timer.value(); + cout<< " Time to solve " << time_solve << endl; + + /* Check the accuracy */ + VectorXd tmp2 = b - A*x; + double tempNorm = tmp2.norm()/b.norm(); + cout << "Relative norm of the computed solution : " << tempNorm <<"\n"; +// cout << "Iterations : " << solver.iterations() << "\n"; + + totaltime.stop(); + cout << "Total time " << totaltime.value() << "\n"; +// std::cout<<x.transpose()<<"\n"; + + return 0; +}
\ No newline at end of file diff --git a/bench/spbench/test_sparseLU.cpp b/bench/spbench/test_sparseLU.cpp new file mode 100644 index 000000000..c6511a9bc --- /dev/null +++ b/bench/spbench/test_sparseLU.cpp @@ -0,0 +1,93 @@ +// Small bench routine for Eigen available in Eigen +// (C) Desire NUENTSA WAKAM, INRIA + +#include <iostream> +#include <fstream> +#include <iomanip> +#include <unsupported/Eigen/SparseExtra> +#include <Eigen/SparseLU> +#include <bench/BenchTimer.h> +#ifdef EIGEN_METIS_SUPPORT +#include <Eigen/MetisSupport> +#endif + +using namespace std; +using namespace Eigen; + +int main(int argc, char **args) +{ +// typedef complex<double> scalar; + typedef double scalar; + SparseMatrix<scalar, ColMajor> A; + typedef SparseMatrix<scalar, ColMajor>::Index Index; + typedef Matrix<scalar, Dynamic, Dynamic> DenseMatrix; + typedef Matrix<scalar, Dynamic, 1> DenseRhs; + Matrix<scalar, Dynamic, 1> b, x, tmp; +// SparseLU<SparseMatrix<scalar, ColMajor>, AMDOrdering<int> > solver; +// #ifdef EIGEN_METIS_SUPPORT +// SparseLU<SparseMatrix<scalar, ColMajor>, MetisOrdering<int> > solver; +// std::cout<< "ORDERING : METIS\n"; +// #else + SparseLU<SparseMatrix<scalar, ColMajor>, COLAMDOrdering<int> > solver; + std::cout<< "ORDERING : COLAMD\n"; +// #endif + + ifstream matrix_file; + string line; + int n; + BenchTimer timer; + + // Set parameters + /* Fill the matrix with sparse matrix stored in Matrix-Market coordinate column-oriented format */ + if (argc < 2) assert(false && "please, give the matrix market file "); + loadMarket(A, args[1]); + cout << "End charging matrix " << endl; + bool iscomplex=false, isvector=false; + int sym; + getMarketHeader(args[1], sym, iscomplex, isvector); +// if (iscomplex) { cout<< " Not for complex matrices \n"; return -1; } + if (isvector) { cout << "The provided file is not a matrix file\n"; return -1;} + if (sym != 0) { // symmetric matrices, only the lower part is stored + SparseMatrix<scalar, ColMajor> temp; + temp = A; + A = temp.selfadjointView<Lower>(); + } + n = A.cols(); + /* Fill the right hand side */ + + if (argc > 2) + loadMarketVector(b, args[2]); + else + { + b.resize(n); + tmp.resize(n); +// tmp.setRandom(); + for (int i = 0; i < n; i++) tmp(i) = i; + b = A * tmp ; + } + + /* Compute the factorization */ +// solver.isSymmetric(true); + 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(); + x = solver.solve(b); + timer.stop(); + cout << "solve time " << timer.value() << std::endl; + /* Check the accuracy */ + Matrix<scalar, Dynamic, 1> tmp2 = b - A*x; + scalar tempNorm = tmp2.norm()/b.norm(); + cout << "Relative norm of the computed solution : " << tempNorm <<"\n"; + cout << "Number of nonzeros in the factor : " << solver.nnzL() + solver.nnzU() << std::endl; + + return 0; +}
\ No newline at end of file |