sparse_trisolver.cpp
Go to the documentation of this file.
00001 
00002 //g++ -O3 -g0 -DNDEBUG  sparse_product.cpp -I.. -I/home/gael/Coding/LinearAlgebra/mtl4/ -DDENSITY=0.005 -DSIZE=10000 && ./a.out
00003 //g++ -O3 -g0 -DNDEBUG  sparse_product.cpp -I.. -I/home/gael/Coding/LinearAlgebra/mtl4/ -DDENSITY=0.05 -DSIZE=2000 && ./a.out
00004 // -DNOGMM -DNOMTL
00005 // -I /home/gael/Coding/LinearAlgebra/CSparse/Include/ /home/gael/Coding/LinearAlgebra/CSparse/Lib/libcsparse.a
00006 
00007 #ifndef SIZE
00008 #define SIZE 10000
00009 #endif
00010 
00011 #ifndef DENSITY
00012 #define DENSITY 0.01
00013 #endif
00014 
00015 #ifndef REPEAT
00016 #define REPEAT 1
00017 #endif
00018 
00019 #include "BenchSparseUtil.h"
00020 
00021 #ifndef MINDENSITY
00022 #define MINDENSITY 0.0004
00023 #endif
00024 
00025 #ifndef NBTRIES
00026 #define NBTRIES 10
00027 #endif
00028 
00029 #define BENCH(X) \
00030   timer.reset(); \
00031   for (int _j=0; _j<NBTRIES; ++_j) { \
00032     timer.start(); \
00033     for (int _k=0; _k<REPEAT; ++_k) { \
00034         X  \
00035   } timer.stop(); }
00036 
00037 typedef SparseMatrix<Scalar,UpperTriangular> EigenSparseTriMatrix;
00038 typedef SparseMatrix<Scalar,RowMajorBit|UpperTriangular> EigenSparseTriMatrixRow;
00039 
00040 void fillMatrix(float density, int rows, int cols,  EigenSparseTriMatrix& dst)
00041 {
00042   dst.startFill(rows*cols*density);
00043   for(int j = 0; j < cols; j++)
00044   {
00045     for(int i = 0; i < j; i++)
00046     {
00047       Scalar v = (internal::random<float>(0,1) < density) ? internal::random<Scalar>() : 0;
00048       if (v!=0)
00049         dst.fill(i,j) = v;
00050     }
00051     dst.fill(j,j) = internal::random<Scalar>();
00052   }
00053   dst.endFill();
00054 }
00055 
00056 int main(int argc, char *argv[])
00057 {
00058   int rows = SIZE;
00059   int cols = SIZE;
00060   float density = DENSITY;
00061   BenchTimer timer;
00062   #if 1
00063   EigenSparseTriMatrix sm1(rows,cols);
00064   typedef Matrix<Scalar,Dynamic,1> DenseVector;
00065   DenseVector b = DenseVector::Random(cols);
00066   DenseVector x = DenseVector::Random(cols);
00067 
00068   bool densedone = false;
00069 
00070   for (float density = DENSITY; density>=MINDENSITY; density*=0.5)
00071   {
00072     EigenSparseTriMatrix sm1(rows, cols);
00073     fillMatrix(density, rows, cols, sm1);
00074 
00075     // dense matrices
00076     #ifdef DENSEMATRIX
00077     if (!densedone)
00078     {
00079       densedone = true;
00080       std::cout << "Eigen Dense\t" << density*100 << "%\n";
00081       DenseMatrix m1(rows,cols);
00082       Matrix<Scalar,Dynamic,Dynamic,Dynamic,Dynamic,RowMajorBit> m2(rows,cols);
00083       eiToDense(sm1, m1);
00084       m2 = m1;
00085 
00086       BENCH(x = m1.marked<UpperTriangular>().solveTriangular(b);)
00087       std::cout << "   colmajor^-1 * b:\t" << timer.value() << endl;
00088 //       std::cerr << x.transpose() << "\n";
00089 
00090       BENCH(x = m2.marked<UpperTriangular>().solveTriangular(b);)
00091       std::cout << "   rowmajor^-1 * b:\t" << timer.value() << endl;
00092 //       std::cerr << x.transpose() << "\n";
00093     }
00094     #endif
00095 
00096     // eigen sparse matrices
00097     {
00098       std::cout << "Eigen sparse\t" << density*100 << "%\n";
00099       EigenSparseTriMatrixRow sm2 = sm1;
00100 
00101       BENCH(x = sm1.solveTriangular(b);)
00102       std::cout << "   colmajor^-1 * b:\t" << timer.value() << endl;
00103 //       std::cerr << x.transpose() << "\n";
00104 
00105       BENCH(x = sm2.solveTriangular(b);)
00106       std::cout << "   rowmajor^-1 * b:\t" << timer.value() << endl;
00107 //       std::cerr << x.transpose() << "\n";
00108 
00109 //       x = b;
00110 //       BENCH(sm1.inverseProductInPlace(x);)
00111 //       std::cout << "   colmajor^-1 * b:\t" << timer.value() << " (inplace)" << endl;
00112 //       std::cerr << x.transpose() << "\n";
00113 //
00114 //       x = b;
00115 //       BENCH(sm2.inverseProductInPlace(x);)
00116 //       std::cout << "   rowmajor^-1 * b:\t" << timer.value() << " (inplace)" << endl;
00117 //       std::cerr << x.transpose() << "\n";
00118     }
00119 
00120 
00121 
00122     // CSparse
00123     #ifdef CSPARSE
00124     {
00125       std::cout << "CSparse \t" << density*100 << "%\n";
00126       cs *m1;
00127       eiToCSparse(sm1, m1);
00128 
00129       BENCH(x = b; if (!cs_lsolve (m1, x.data())){std::cerr << "cs_lsolve failed\n"; break;}; )
00130       std::cout << "   colmajor^-1 * b:\t" << timer.value() << endl;
00131     }
00132     #endif
00133 
00134     // GMM++
00135     #ifndef NOGMM
00136     {
00137       std::cout << "GMM++ sparse\t" << density*100 << "%\n";
00138       GmmSparse m1(rows,cols);
00139       gmm::csr_matrix<Scalar> m2;
00140       eiToGmm(sm1, m1);
00141       gmm::copy(m1,m2);
00142       std::vector<Scalar> gmmX(cols), gmmB(cols);
00143       Map<Matrix<Scalar,Dynamic,1> >(&gmmX[0], cols) = x;
00144       Map<Matrix<Scalar,Dynamic,1> >(&gmmB[0], cols) = b;
00145 
00146       gmmX = gmmB;
00147       BENCH(gmm::upper_tri_solve(m1, gmmX, false);)
00148       std::cout << "   colmajor^-1 * b:\t" << timer.value() << endl;
00149 //       std::cerr << Map<Matrix<Scalar,Dynamic,1> >(&gmmX[0], cols).transpose() << "\n";
00150 
00151       gmmX = gmmB;
00152       BENCH(gmm::upper_tri_solve(m2, gmmX, false);)
00153       timer.stop();
00154       std::cout << "   rowmajor^-1 * b:\t" << timer.value() << endl;
00155 //       std::cerr << Map<Matrix<Scalar,Dynamic,1> >(&gmmX[0], cols).transpose() << "\n";
00156     }
00157     #endif
00158 
00159     // MTL4
00160     #ifndef NOMTL
00161     {
00162       std::cout << "MTL4\t" << density*100 << "%\n";
00163       MtlSparse m1(rows,cols);
00164       MtlSparseRowMajor m2(rows,cols);
00165       eiToMtl(sm1, m1);
00166       m2 = m1;
00167       mtl::dense_vector<Scalar> x(rows, 1.0);
00168       mtl::dense_vector<Scalar> b(rows, 1.0);
00169 
00170       BENCH(x = mtl::upper_trisolve(m1,b);)
00171       std::cout << "   colmajor^-1 * b:\t" << timer.value() << endl;
00172 //       std::cerr << x << "\n";
00173 
00174       BENCH(x = mtl::upper_trisolve(m2,b);)
00175       std::cout << "   rowmajor^-1 * b:\t" << timer.value() << endl;
00176 //       std::cerr << x << "\n";
00177     }
00178     #endif
00179 
00180 
00181     std::cout << "\n\n";
00182   }
00183   #endif
00184 
00185   #if 0
00186     // bench small matrices (in-place versus return bye value)
00187     {
00188       timer.reset();
00189       for (int _j=0; _j<10; ++_j) {
00190         Matrix4f m = Matrix4f::Random();
00191         Vector4f b = Vector4f::Random();
00192         Vector4f x = Vector4f::Random();
00193         timer.start();
00194         for (int _k=0; _k<1000000; ++_k) {
00195           b = m.inverseProduct(b);
00196         }
00197         timer.stop();
00198       }
00199       std::cout << "4x4 :\t" << timer.value() << endl;
00200     }
00201 
00202     {
00203       timer.reset();
00204       for (int _j=0; _j<10; ++_j) {
00205         Matrix4f m = Matrix4f::Random();
00206         Vector4f b = Vector4f::Random();
00207         Vector4f x = Vector4f::Random();
00208         timer.start();
00209         for (int _k=0; _k<1000000; ++_k) {
00210           m.inverseProductInPlace(x);
00211         }
00212         timer.stop();
00213       }
00214       std::cout << "4x4 IP :\t" << timer.value() << endl;
00215     }
00216   #endif
00217 
00218   return 0;
00219 }
00220 


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:33:32