diff options
author | Desire NUENTSA <desire.nuentsa_wakam@inria.fr> | 2012-06-14 18:45:04 +0200 |
---|---|---|
committer | Desire NUENTSA <desire.nuentsa_wakam@inria.fr> | 2012-06-14 18:45:04 +0200 |
commit | 0c9b08e46e7507d9f13200f0702bc57ed6aae52c (patch) | |
tree | bc319fe32b4bdfa8dd601082e25c508606f53853 /bench/spbench | |
parent | f8a0745cb0426eb3095dbea24288a64eddab04f0 (diff) |
build complete... almost
Diffstat (limited to 'bench/spbench')
-rw-r--r-- | bench/spbench/test_sparseLU.cpp | 64 |
1 files changed, 64 insertions, 0 deletions
diff --git a/bench/spbench/test_sparseLU.cpp b/bench/spbench/test_sparseLU.cpp new file mode 100644 index 000000000..0bbbb0627 --- /dev/null +++ b/bench/spbench/test_sparseLU.cpp @@ -0,0 +1,64 @@ +// 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> + +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; + SparseLU<SparseMatrix<double, ColMajor>, AMDOrdering<double, int> > solver; + ifstream matrix_file; + string line; + int n; + + // 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<double, 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.compute(A); + + solver._solve(b, x); + /* Check the accuracy */ + VectorXd tmp2 = b - A*x; + double tempNorm = tmp2.norm()/b.norm(); + cout << "Relative norm of the computed solution : " << tempNorm <<"\n"; + + return 0; +}
\ No newline at end of file |