1 // Small bench routine for Eigen available in Eigen
2 // (C) Desire NUENTSA WAKAM, INRIA
3
4 #include <iostream>
5 #include <fstream>
6 #include <iomanip>
7 #include <Eigen/Jacobi>
8 #include <Eigen/Householder>
9 #include <Eigen/IterativeLinearSolvers>
10 #include <Eigen/LU>
11 #include <unsupported/Eigen/SparseExtra>
12 //#include <Eigen/SparseLU>
13 #include <Eigen/SuperLUSupport>
14 // #include <unsupported/Eigen/src/IterativeSolvers/Scaling.h>
15 #include <bench/BenchTimer.h>
16 #include <unsupported/Eigen/IterativeSolvers>
17 using namespace std;
18 using namespace Eigen;
19
main(int argc,char ** args)20 int main(int argc, char **args)
21 {
22 SparseMatrix<double, ColMajor> A;
23 typedef SparseMatrix<double, ColMajor>::Index Index;
24 typedef Matrix<double, Dynamic, Dynamic> DenseMatrix;
25 typedef Matrix<double, Dynamic, 1> DenseRhs;
26 VectorXd b, x, tmp;
27 BenchTimer timer,totaltime;
28 //SparseLU<SparseMatrix<double, ColMajor> > solver;
29 // SuperLU<SparseMatrix<double, ColMajor> > solver;
30 ConjugateGradient<SparseMatrix<double, ColMajor>, Lower,IncompleteCholesky<double,Lower> > solver;
31 ifstream matrix_file;
32 string line;
33 int n;
34 // Set parameters
35 // solver.iparm(IPARM_THREAD_NBR) = 4;
36 /* Fill the matrix with sparse matrix stored in Matrix-Market coordinate column-oriented format */
37 if (argc < 2) assert(false && "please, give the matrix market file ");
38
39 timer.start();
40 totaltime.start();
41 loadMarket(A, args[1]);
42 cout << "End charging matrix " << endl;
43 bool iscomplex=false, isvector=false;
44 int sym;
45 getMarketHeader(args[1], sym, iscomplex, isvector);
46 if (iscomplex) { cout<< " Not for complex matrices \n"; return -1; }
47 if (isvector) { cout << "The provided file is not a matrix file\n"; return -1;}
48 if (sym != 0) { // symmetric matrices, only the lower part is stored
49 SparseMatrix<double, ColMajor> temp;
50 temp = A;
51 A = temp.selfadjointView<Lower>();
52 }
53 timer.stop();
54
55 n = A.cols();
56 // ====== TESTS FOR SPARSE TUTORIAL ======
57 // cout<< "OuterSize " << A.outerSize() << " inner " << A.innerSize() << endl;
58 // SparseMatrix<double, RowMajor> mat1(A);
59 // SparseMatrix<double, RowMajor> mat2;
60 // cout << " norm of A " << mat1.norm() << endl; ;
61 // PermutationMatrix<Dynamic, Dynamic, int> perm(n);
62 // perm.resize(n,1);
63 // perm.indices().setLinSpaced(n, 0, n-1);
64 // mat2 = perm * mat1;
65 // mat.subrows();
66 // mat2.resize(n,n);
67 // mat2.reserve(10);
68 // mat2.setConstant();
69 // std::cout<< "NORM " << mat1.squaredNorm()<< endl;
70
71 cout<< "Time to load the matrix " << timer.value() <<endl;
72 /* Fill the right hand side */
73
74 // solver.set_restart(374);
75 if (argc > 2)
76 loadMarketVector(b, args[2]);
77 else
78 {
79 b.resize(n);
80 tmp.resize(n);
81 // tmp.setRandom();
82 for (int i = 0; i < n; i++) tmp(i) = i;
83 b = A * tmp ;
84 }
85 // Scaling<SparseMatrix<double> > scal;
86 // scal.computeRef(A);
87 // b = scal.LeftScaling().cwiseProduct(b);
88
89 /* Compute the factorization */
90 cout<< "Starting the factorization "<< endl;
91 timer.reset();
92 timer.start();
93 cout<< "Size of Input Matrix "<< b.size()<<"\n\n";
94 cout<< "Rows and columns "<< A.rows() <<" " <<A.cols() <<"\n";
95 solver.compute(A);
96 // solver.analyzePattern(A);
97 // solver.factorize(A);
98 if (solver.info() != Success) {
99 std::cout<< "The solver failed \n";
100 return -1;
101 }
102 timer.stop();
103 float time_comp = timer.value();
104 cout <<" Compute Time " << time_comp<< endl;
105
106 timer.reset();
107 timer.start();
108 x = solver.solve(b);
109 // x = scal.RightScaling().cwiseProduct(x);
110 timer.stop();
111 float time_solve = timer.value();
112 cout<< " Time to solve " << time_solve << endl;
113
114 /* Check the accuracy */
115 VectorXd tmp2 = b - A*x;
116 double tempNorm = tmp2.norm()/b.norm();
117 cout << "Relative norm of the computed solution : " << tempNorm <<"\n";
118 // cout << "Iterations : " << solver.iterations() << "\n";
119
120 totaltime.stop();
121 cout << "Total time " << totaltime.value() << "\n";
122 // std::cout<<x.transpose()<<"\n";
123
124 return 0;
125 }