Использование LU-решателя boost, но избегать вычисления обратной матрицы?

Я решаю систему линейных уравнений в C ++, используя разложение LU, как это предусмотрено 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);
}

Похоже, что lu_substitute вычисляет обратную величину входного matix, что очень дорого (как обсуждалось выше) Вот).

Есть ли способ избежать этого, используя функцию повышения?