Using boost's LU-solver but avoid computing a matrix inverse?

Question

I am solving system of linear equations in C++ using the LU decomposition as provided by Boost.

static void boostLU(const boost::numeric::ublas::matrix<double> &A, const boost::numeric::ublas::matrix<double> &y)
{
    boost::numeric::ublas::matrix<double> Afactorized = A;  
    boost::numeric::ublas::matrix<double> Ainv = boost::numeric::ublas::identity_matrix<float>(A.size1());
    boost::numeric::ublas::vector<double> x_boost(Afactorized.size1(), 1);
    boost::numeric::ublas::permutation_matrix<size_t> pm(Afactorized.size1());
    boost::numeric::ublas::matrix<double> result = boost::numeric::ublas::identity_matrix<float>(A.size1());

    int singular = boost::numeric::ublas::lu_factorize(Afactorized,pm);
    if (singular)
    {
        throw std::runtime_error("[LinearSolver<LU>::solve()] A is singular.");
    } 
    result = y;
    boost::numeric::ublas::lu_substitute(Afactorized, pm, result);
}

It seems that lu_substitute calculates the inverse of the input matix, which is computationall expensive (as discussed here).

Is there any way to avoid that using boost functionality?

Answer 1

lu_substitute does not compute the inverse.

Look at the source code (http://www.boost.org/doc/libs/1_55_0/boost/numeric/ublas/lu.hpp): lu_substitute calls inplace_solve, and inplace_solve (defined here: http://www.boost.org/doc/libs/1_55_0/boost/numeric/ublas/triangular.hpp) does inplace forward/backwards substitution. So everything is as efficient as it gets.