Matrix inverse from PLU decomposition

Question

I am trying to compute the matrix inverse from plu decomposition, but the results I get are wrong. Using the debugger, I couldn't find a single step I could blame, so here is the code for the inversion function :

template<class T>
Matrix<T> Matrix<T>::inv() const {
    std::tuple<Matrix<T>, Matrix<T>, Matrix<T>> PLU(this->decomp_PLU());
    Matrix<T> P(std::get<0>(PLU));
    Matrix<T> L(std::get<1>(PLU));
    Matrix<T> U(std::get<2>(PLU));

    if (!this->lines() == this->cols()) { throw std::logic_error("Matrix must be square"); }
    unsigned N = this->lines();

    Matrix<T> Y(N, N);
    for (unsigned i = 0; i < N; i++) {
        Matrix<T> B(Matrix<T>::gen_col(P.transp().data[i])); // B is the ith column vector of P
        Y[i] = Matrix<T>::solve_climb(L, B).transp().data[0]; // Compute LY=P for each column vector

        Matrix<T> conf(L.dot(Matrix<T>::gen_col(Y[i])));
        if (!conf.allclose(B, 1e-9, 1e-9)) {
            throw std::logic_error("WHYY");
        }
    }
    Y = Y.transp();

    Matrix<T> X(N, N);
    for (unsigned i = 0; i < N; i++) {
        Matrix<T> B(Matrix<T>::gen_col(Y.transp().data[i])); // B is the ith column vector of Y
        X[i] = Matrix<T>::solve_descent(U, B).transp().data[0]; // Compute UX=Y for each column vector

        Matrix<T> conf(U.dot(Matrix<T>::gen_col(X[i])));
        if (!conf.allclose(B, 1e-9, 1e-9)) {
            throw std::logic_error("WHYY");
        }
    }
    X = X.transp();

    return X;
}

The function Matrix<T>::gen_col generates a column vector from its input, and solve_climb / solve_descent each compute the column vector that solves AX=B where A is a triangular matrix (inferior or superior depending on the case)

FYI, The code throws the logic errors 'WHYY' when trying to check the calculations are right for each vector.

Any clue of where the code could be wrong ?

Thank you

EDIT : The full code is here (matrix_def.h) and here (matrix.h)

Answer 1

As L is triangular inferior, you should use solve_descent, and not solve_climb. The same thing goes for U (triangular superior), you need to use solve_climb.