LU факторизация в программе (C++)

Question

Мне нужна функция или класс, который выполняет факторизацию LU (разложение, в чем разница?) В моей программе на C++. Я использую в качестве компилятора последнюю версию Dev-C (на Windows). Я подумал, как установить броненосец и лапак, но это кажется действительно сложным и с небольшим количеством проблем http://icl.cs.utk.edu/lapack-for-windows/lapack/ (как сказано в нижней части страницы), Поэтому я хотел бы иметь библиотеку (на C++), которая работает хорошо, и ее не так сложно установить. Я нашел кое-что об Эйгене в примере, как оно? Есть еще предложения?

Спасибо

PS матрица плотная на диагонали и около нее, в остальном разреженная, а часть под углом (NE) пустая.

Answer 1

Я написал немного кода на C для вас. Может быть, это может помочь. Вот процедура разложения плотной матрицы на L, U, где L*U=A, L - нижняя треугольная, U - верхняя треугольная, L[i,i]=U[i,i] (диагональные элементы равны), Такое разложение также известно как LU(sq).

#include <math.h>
// A, L, U each allocates at least N*N doubles
// A contains elements of given matrix, written row by row
void decompose(unsigned N, double *A, double *L, double *U) {
#define _(M,i,j) M[N*i + j]
#define _A(i,j) _(A,i,j)
#define _L(i,j) _(L,i,j)
#define _U(i,j) _(U,i,j)
    _L(0,0) = sqrt(_A(0,0));
    _U(0,0) = sqrt(_A(0,0));
    for ( int i = 0; i < N; ++i ) {
        _L(i,0) = _A(i,0) / _A(0,0);
        _U(0,i) = _A(0,i) / _A(0,0);
        _L(0,i) = _U(i,0) = 0.0;

        double s = 0.0;
        for ( int k = 0; k < i; ++k ) {
            s += _L(i,k) * _U(k,i);
        }
        _L(i,i) = _U(i,i) = sqrt(_A(i,i) - s);

        for ( int j = 1; j < i; ++j ) {
            double s = 0.0;
            for ( int k = 0; k < j; ++k ) {
                s += _L(i,k) * _U(k,j);
            }
            _L(i,j) = (_A(i,j) - s) / _L(i,i);
            _L(j,i) = 0.0;

            double s = 0.0;
            for ( int k = 0; k < i; ++k ) {
                s += _L(i,k) * _U(k,j);
            }
            _U(j,i) = (_A(i,j) - s) / _U(i,i);
            _U(i,j) = 0.0;
        }
    }
}

К сожалению, сейчас у меня нет времени проверять код на наличие ошибок. У меня есть некоторые формулы, в которых я уверен (использовался несколько лет назад):

Мой код основан на этих формулах. Конечно, если вы хотите использовать какой-то конкретный подходящий для вашего разреженного формата данных, вам нужно изменить процедуру, но не формулы.