// no automatic singular check :
#define BOOST_UBLAS_TYPE_CHECK 0
#include <boost/numeric/ublas/io.hpp> 
#include <boost/numeric/ublas/lu.hpp> 

int main () {
  using namespace boost::numeric::ublas; 
  using std::cout;
  using std::endl;
  typedef double mtype; 
  matrix<mtype> A(3,3), B(identity_matrix<mtype>(3)); 
  for (size_t i=0; i<A.size1(); ++i) 
    for (size_t j=0; j<A.size2(); ++j) {
      //A(i,j)=  A(i,j)=3*(i+4)/(7-i-j);  // inv. ok
      A(i,j) = 3*i*(j+3) + 5*j + 2; // inv. fails, singular
    }     
  cout << "  A = "<< A << endl; 
  cout << "  B = "<< B << endl; 
  permutation_matrix<> pm(3);     
  int result = lu_factorize<matrix<mtype>, permutation_matrix<> >(A, pm); 
  cout << " LU = " << A << endl; 
  cout << " pm = " << pm << endl; 
  if ( result==0 ) { 
    lu_substitute(A, pm, B);      
    cout << "inv = " << B << endl; 
  } else { 
    cout << "singular, at row " << result << endl; 
  } 
  return EXIT_SUCCESS; 
}