Матрица обратная от разложения PLU - PullRequest
Новая
       92

Матрица обратная от разложения PLU

2 голосов
/ 11 апреля 2019

Я пытаюсь вычислить матрицу, обратную от разложения по плюсу, но полученные результаты неверны. Используя отладчик, я не смог найти ни одного шага, который мог бы винить, поэтому вот код для функции инверсии: 

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;
}

Функция Matrix<T>::gen_col генерирует вектор столбца из своих входных данных, и solve_climb / solve_descent каждый вычисляет вектор столбца, который решает AX = B, где A - треугольная матрица (младшая или старшая в зависимости от случая)

FYI, код выдает логические ошибки «WHYY» при попытке проверить правильность вычислений для каждого вектора.

Любая подсказка, где код может быть неправильным?

Спасибо

РЕДАКТИРОВАТЬ: полный код здесь (matrix_def.h) и здесь (matrix.h)

1 Ответ

2 голосов
/ 11 апреля 2019

Поскольку L имеет треугольную форму, вы должны использовать solve_descent, а не solve_climb. То же самое касается U (треугольный супериор), вам нужно использовать solve_climb.

...