Home | History | Annotate | Download | only in spbench 
      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 
     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 }