1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// g++ -I.. sparse_lu.cpp -O3 -g0 -I /usr/include/superlu/ -lsuperlu -lgfortran -DSIZE=1000 -DDENSITY=.05 && ./a.out 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_SUPERLU_SUPPORT 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_UMFPACK_SUPPORT 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include <Eigen/Sparse> 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define NOGMM 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define NOMTL 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef SIZE 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define SIZE 10 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef DENSITY 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define DENSITY 0.01 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef REPEAT 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define REPEAT 1 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include "BenchSparseUtil.h" 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef MINDENSITY 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define MINDENSITY 0.0004 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef NBTRIES 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define NBTRIES 10 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define BENCH(X) \ 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath timer.reset(); \ 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int _j=0; _j<NBTRIES; ++_j) { \ 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath timer.start(); \ 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int _k=0; _k<REPEAT; ++_k) { \ 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath X \ 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } timer.stop(); } 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtypedef Matrix<Scalar,Dynamic,1> VectorX; 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include <Eigen/LU> 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<int Backend> 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid doEigen(const char* name, const EigenSparseMatrix& sm1, const VectorX& b, VectorX& x, int flags = 0) 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cout << name << "..." << std::flush; 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath BenchTimer timer; timer.start(); 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseLU<EigenSparseMatrix,Backend> lu(sm1, flags); 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath timer.stop(); 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (lu.succeeded()) 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cout << ":\t" << timer.value() << endl; 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cout << ":\t FAILED" << endl; 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return; 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool ok; 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath timer.reset(); timer.start(); 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ok = lu.solve(b,&x); 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath timer.stop(); 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (ok) 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cout << " solve:\t" << timer.value() << endl; 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cout << " solve:\t" << " FAILED" << endl; 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath //std::cout << x.transpose() << "\n"; 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathint main(int argc, char *argv[]) 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int rows = SIZE; 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int cols = SIZE; 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath float density = DENSITY; 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath BenchTimer timer; 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorX b = VectorX::Random(cols); 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorX x = VectorX::Random(cols); 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool densedone = false; 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath //for (float density = DENSITY; density>=MINDENSITY; density*=0.5) 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// float density = 0.5; 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenSparseMatrix sm1(rows, cols); 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fillMatrix(density, rows, cols, sm1); 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // dense matrices 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath #ifdef DENSEMATRIX 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (!densedone) 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath densedone = true; 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cout << "Eigen Dense\t" << density*100 << "%\n"; 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseMatrix m1(rows,cols); 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eiToDense(sm1, m1); 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath BenchTimer timer; 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath timer.start(); 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath FullPivLU<DenseMatrix> lu(m1); 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath timer.stop(); 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cout << "Eigen/dense:\t" << timer.value() << endl; 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath timer.reset(); 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath timer.start(); 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lu.solve(b,&x); 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath timer.stop(); 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cout << " solve:\t" << timer.value() << endl; 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// std::cout << b.transpose() << "\n"; 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// std::cout << x.transpose() << "\n"; 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath #endif 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath #ifdef EIGEN_UMFPACK_SUPPORT 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x.setZero(); 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath doEigen<Eigen::UmfPack>("Eigen/UmfPack (auto)", sm1, b, x, 0); 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath #endif 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath #ifdef EIGEN_SUPERLU_SUPPORT 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x.setZero(); 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath doEigen<Eigen::SuperLU>("Eigen/SuperLU (nat)", sm1, b, x, Eigen::NaturalOrdering); 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// doEigen<Eigen::SuperLU>("Eigen/SuperLU (MD AT+A)", sm1, b, x, Eigen::MinimumDegree_AT_PLUS_A); 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// doEigen<Eigen::SuperLU>("Eigen/SuperLU (MD ATA)", sm1, b, x, Eigen::MinimumDegree_ATA); 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath doEigen<Eigen::SuperLU>("Eigen/SuperLU (COLAMD)", sm1, b, x, Eigen::ColApproxMinimumDegree); 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath #endif 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 133