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

//#include <boost/numeric/ublas/vector_proxy.hpp>

namespace ag { namespace linalg {
   namespace ublas = boost::numeric::ublas;

   void solve( ublas::matrix< float >& A,
               ublas::vector< float >& b )
   {
      //ublas::permutation_matrix<> pivot( 3 );
      
      ublas::matrix< float > triangular( A );
      const int singular = 0;//ublas::lu_factorize( triangular/*, pivot*/ );
      if( singular == 0 )
      {
         ublas::lu_substitute< ublas::matrix< float >, ublas::vector< float > >( triangular, /*pivot, */b );
      }
      else
      {
         std::cerr << "singular: " << singular << std::endl;
         // do something with A
      }
   }
} }