[ublas] Solving non-triangular Ax=b

Subject: [ublas] Solving non-triangular Ax=b
From: Kenneth Porter (shiva_at_[hidden])
Date: 2009-04-30 16:54:52

In figuring out how to do this I came across this list message:

<http://lists.boost.org/MailArchives/ublas/2005/01/0062.php>

In digging through the ublas code I figured out that permutation_matrix is 
a matrix of indices, not values, so I was getting warnings about loss of 
precision between doubles and unsigneds in some internal std::swap 
operations. The correct code should look like this:

namespace ublas = boost::numeric::ublas;

template<class T>
void nonTriangularSolve(ublas::matrix<T>& A,
                        ublas::vector<T>& x,
                        const ublas::vector<T>& b)
{
   ublas::permutation_matrix<ublas::matrix<T>::size_type> P(b.size());
   ublas::lu_factorize(A,P);
   x = b;
   ublas::lu_substitute(A,P,x);
}

Note the template parameter for the permutation_matrix should be size_type, 
not double.

The matrix A is logically const but lu_factorize requires it to be mutable, 
presumably for performance.