(C++) boost::numeric::ublas::matrix example 7: determinant for matrix of any size

Technical facts

Application type(s)

Operating system(s) or programming environment(s)

Libraries used:

• STL: GNU ISO C++ Library, version 4.7.2

Qt project file: CppUblasMatrixExample7.pro

 ``` TEMPLATE = app CONFIG += console CONFIG -= qt QMAKE_CXXFLAGS += -std=c++11 SOURCES += main.cpp ```

main.cpp

 ``` #ifdef _WIN32 #undef __STRICT_ANSI__ #endif #include #include #include #include #include ///Calculate the determinant //Adapted from the code Maik Beckmann posted at //  http://boost.2283326.n4.nabble.com/How-to-compute-determinant-td2710896.html double CalcDeterminant(boost::numeric::ublas::matrix m) {   assert(m.size1() == m.size2() && "Can only calculate the determinant of square matrices");   boost::numeric::ublas::permutation_matrix pivots(m.size1() );   const int is_singular = boost::numeric::ublas::lu_factorize(m, pivots);   if (is_singular) return 0.0;   double d = 1.0;   const std::size_t sz = pivots.size();   for (std::size_t i=0; i != sz; ++i)   {     if (pivots(i) != i)     {       d *= -1.0;     }     d *= m(i,i);   }   return d; } double CalcDeterminantSmall(const boost::numeric::ublas::matrix& m) {   assert(m.size1() == m.size2() && "Can only calculate the determinant of square matrices");   switch(m.size1())   {     case 0: return 1.0;     case 1:     {       return m(0,0);     }     case 2:     {       const double a = m(0,0);       const double b = m(0,1);       const double c = m(1,0);       const double d = m(1,1);       const double determinant = (a * d) - (b * c);       return determinant;     }     case 3:     {       assert(m.size1() == 3 && m.size2() == 3 && "Only for 3x3 matrices");       const double a = m(0,0);       const double b = m(0,1);       const double c = m(0,2);       const double d = m(1,0);       const double e = m(1,1);       const double f = m(1,2);       const double g = m(2,0);       const double h = m(2,1);       const double k = m(2,2);       const double determinant         = (a * ((e*k) - (f*h)))         - (b * ((k*d) - (f*g)))         + (c * ((d*h) - (e*g)));       return determinant;     }     default:       assert("Cannot handle matrix bigger than 3x3");       throw std::logic_error("Cannot handle matrix bigger than 3x3");   } } ///Chop returns a std::vector of sub-matrices //[ A at [0]   B at [1] ] //[ C at [2]   D at [4] ] const std::vector > Chop(   const boost::numeric::ublas::matrix& m) {   using boost::numeric::ublas::range;   using boost::numeric::ublas::matrix;   using boost::numeric::ublas::matrix_range;   std::vector > v;   v.reserve(4);   const int midy = m.size1() / 2;   const int midx = m.size2() / 2;   const matrix_range > top_left(    m,range(0   ,midy     ),range(0   ,midx     ));   const matrix_range > bottom_left( m,range(midy,m.size1()),range(0   ,midx     ));   const matrix_range > top_right(   m,range(0   ,midy     ),range(midx,m.size2()));   const matrix_range > bottom_right(m,range(midy,m.size1()),range(midx,m.size2()));   v.push_back(matrix(top_left));   v.push_back(matrix(top_right));   v.push_back(matrix(bottom_left));   v.push_back(matrix(bottom_right));   return v; } const boost::numeric::ublas::matrix CreateRandomMatrix(const std::size_t n_rows, const std::size_t n_cols) {   boost::numeric::ublas::matrix m(n_rows,n_cols);   for (std::size_t row=0; row!=n_rows; ++row)   {     for (std::size_t col=0; col!=n_cols; ++col)     {       m(row,col) = static_cast(std::rand()) / static_cast(RAND_MAX);     }   }   return m; } int main() {   using boost::numeric::ublas::matrix;   //Compare CalcDeterminant for smaller matrices   for (std::size_t sz = 0; sz!=4; ++sz)   {     const matrix m = CreateRandomMatrix(sz,sz);     const double determinant = CalcDeterminant(m);     const double determinant_small = CalcDeterminantSmall(m);     const double epsilon = 0.0001; //Rounding error     assert(std::abs(determinant - determinant_small) < epsilon);   }   //Calculate larger and larger matrices their determinants   //Notice that the computation takes longer and longer   for (std::size_t sz = 0; sz!=10; ++sz)   {     const matrix m = CreateRandomMatrix(sz,sz);     const double determinant = CalcDeterminant(m);     std::cout << determinant << '\n';   } } /* Screen output: 1 0.95223 -0.326259 -0.0324396 -0.0537625 0.0014962 0.023256 0.0172371 -0.0193395 0.00152402 */ ```