#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <cstdlib>

int main()
{
	//using namespace boost::numeric::ublas;
	//using namespace std;

    int size = 4;
  
//    matrix<float>  A(4, 4),  B(identity_matrix<float>(4));                     // OK with 1, 2, 3, 4.  Crash with 5, 6, ... (Floating exception)
    //matrix<float>  A(size, size),  B(identity_matrix<float>(size));          // I will get a compile time error here ***, even if size is a const int
    boost::numeric::ublas::matrix<float>  A(size, size);
	boost::numeric::ublas::vector<float> B(boost::numeric::ublas::scalar_vector<float>(size, 1));
   
    for (int i = 0; i < size; ++i)
        for (int j = 0; j < size; ++j)
            A(i,j) = std::rand() % 10;
       
    boost::numeric::ublas::permutation_matrix<> pm(size);
   
    int result = boost::numeric::ublas::lu_factorize<boost::numeric::ublas::matrix<float>, boost::numeric::ublas::permutation_matrix<> >(A, pm);
   
    if ( 0 == result ) {
        boost::numeric::ublas::lu_substitute(A, pm, B);           // *** Here is where the error is generated
    }
    else
        std::cout << "Error at row " << result << std::endl;
}