Boost logo

Ublas :

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:


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());
   x = b;

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.