Metody Numeryczne - Rozkład LU

Macierz trójkątna ( ang. triangular matrix ) jest specjalnym rodzajem macierzy kwadratowej. Występuje ona w dwóch odmianach:

Macierz trójkątna dolna (ang. lower triangular matrix, L ) posiada wszystkie wyrazy leżące ponad główną przekątna równe zero:

Macierz trójkątna górna ( ang. upper triangular matrix, U ) posiada wszystkie wyrazy leżące poniżej przekątnej głównej równe zero:

Rozkład LU ( ang. LU decomposition ) polega na rozłożeniu macierzy kwadratowej A na iloczyn dwóch macierzy trójkątnych, dolnej L o wyrazach przekątnej głównej równych 1 oraz górnej U:

Zgodnie z własnościami wyznaczników, wyznacznik iloczynu macierzy jest równy iloczynowi ich wyznaczników:

Z kolei wyznacznik macierzy trójkątnej jest równy iloczynowi wyrazów na przekątnej głównej, co wynika bezpośrednio z rozwinięcia Laplace'a. Zatem:

Wyznacznik macierzy A jest równy wyznacznikowi macierzy U. Z kolei wyznacznik macierzy U, która jest macierzą trójkątną, jest równy iloczynowi wyrazów na głównej przekątnej macierzy U:

I ostatecznie:

Wynika z tego, iż jeśli uda nam się znaleźć macierz U dla macierzy A, to wyznacznik macierzy A policzymy w prosty sposób przez wymnożenie przez siebie wyrazów głównej przekątnej macierzy U. Zajmijmy się zatem sposobem wyznaczenia macierzy U.

Na początek: podrozdziału strony

Ponieważ iloczyn macierzy L i U musi dać w wyniku macierz A, to każdy wyraz macierzy A jest sumą iloczynów wyrazów z odpowiedniego wiersza macierzy L i odpowiedniej kolumny macierzy U:

Na przykład wyraz a₁₁ jest sumą iloczynów kolejnych wyrazów z wiersza nr 1 macierzy L przez kolejne wyrazy z kolumny nr 1 macierzy U:

gdzie n jest stopniem macierzy

Ponieważ obie macierze L i U są macierzami trójkątnymi, to część ich elementów ma wartość zero. Np. w macierzy L w wierszu nr 1 tylko pierwszy wyraz jest niezerowy i równy 1, pozostałe wyrazy tego wiersza są równe 0, zatem:

Pozwala to wyznaczyć pierwszy wyraz macierzy U, który jest taki sam jak wyraz macierzy A.

Rozważmy teraz dwa możliwe przypadki dla wyrazu a_ij:

Indeksy: i ≤ j

Wyraz a_ij leży na głównej przekątnej macierzy A lub ponad nią. Ostatnim niezerowym elementem w macierzy L w wierszu i-tym jest l_ii. Wzór iloczynowy możemy więc zredukować:

Co więcej, element l_ii jest elementem przekątnej głównej macierzy L i ma wartość 1, zatem:

Wzór ten pozwala wyznaczyć element u_ij, jeśli znamy poprzednie elementy macierzy L i U:

Indeksy: i > j

Wyraz a_ij leży pod główną przekątną macierzy A. Ostatnim potencjalnie niezerowym elementem w kolumnie j-tej w macierzy U jest element u_jj. Wszystkie pozostałe elementy macierzy U w tej kolumnie mają wartość zero. Wzór iloczynowy redukuje się do:

Wyprowadźmy poza sumę jej ostatni wyraz:

Wzór ten pozwala wyliczyć wyraz l_ij macierzy L, jeśli znane są pozostałe wyrazy macierzy L i U:

Otrzymaliśmy wzory na elementy obu macierzy L i U rozkładu LU macierzy A. Wzory te wykorzystuje znany algorytm Doolitle'a, który wygląda następująco:

Tworzymy dwie macierze L i U o stopniu n macierzy A. Wypełniamy je zerami. W macierzy L główną przekątną wypełniamy elementami o wartości 1.

Dla j od 1 do n wykonujemy dwie poniższe operacje:

Dla i od 1 do j wyznaczamy u_ij wg wzoru:

Dla i od j = 1 do n wyznaczamy l_ij wg wzoru:

Uwaga: obliczenia powinny być wykonywane na liczbach zmiennoprzecinkowych.

Przyjrzyj się drugiemu ze wzorów w algorytmie Doolitle'a. Występuje tutaj dzielenie przez element u_jj. Jeśli u_jj = 0, dzielenie jest niewykonalne i nie można wyznaczyć rozkładu LU macierzy A. Problemem tym zajmiemy się dalej w rozdziale.

Metoda Doolitle'a ma złożoność obliczeniową klasy O(n³), co sprawia iż nadaje się ona do obliczeń wyznaczników nawet dużych macierzy.

Poniższy program wyznacza rozkład LU macierzy wejściowej A, a następnie oblicza jej wyznacznik mnożąc przez siebie elementy głównej przekątnej macierzy U. Macierz A przekazywana jest przez standardowe wejście w następującej postaci:

pierwsza liczba określa stopień macierzy, następne n × n liczb to zawartość kolejnych wierszy macierzy A. Program wyświetla odczytaną macierz oraz wartość jej wyznacznika.

Przykładowe dane ( przekopiuj je do schowka i wklej w programie ):

Macierz:

Dane wejściowe: 5 5 3 7 4 2 9 2 2 1 1 3 6 2 8 9 9 4 -2 -1 -3 0 5 3 -6 -11

Przykładowy program w języku C++

// Macierze
// (C)2019 mgr Jerzy Wałaszek
// Metody numeryczne
//---------------------------

#include <iostream>
#include <iomanip>

using namespace std;

// Definicja klasy, może być w pliku nagłówkowym
//----------------------------------------------

#ifndef _matrix_class
    #define _matrix_class

    template <class T> class matrix
    {
        private:
            T * A;                   // Wyrazy macierzy
        public:
            int n;                   // Stopień macierzy
            matrix();                // Konstruktor
            matrix(int nn);          // Konstruktor
            ~matrix();               // Destruktor
            T  getv(int r, int c);
            void setv(int r, int c, T v);
            void set1();             // ustawia elementy przekątnej na 1
            T det();                 // Oblicza wyznacznik
    };

    // Konstruktor - tworzy macierz ze strumienia cin
    //-----------------------------------------------
    template <class T> matrix<T>::matrix()
    {
        cin >> n;       // Stopień macierzy
        int nn = n * n;
        A = new T [nn]; // Rezerwujemy pamięć na macierz n x n

        for(int i = 0; i < nn; i++) cin >> A[i]; // Wczytujemy wyrazy macierzy
    }

    // Konstruktor - tworzy macierz stopnia n i wypełnia ją zerami
    //------------------------------------------------------------
    template <class T> matrix<T>::matrix(int nn)
    {
        n = nn;       // Stopień macierzy
        nn = n * n;
        A = new T [nn]; // Rezerwujemy pamięć na macierz n x n

        for(int i = 0; i < nn; i++) A[i] = (T)0;
    }

    // Destruktor - usuwa macierz i zwraca zajętą przez nią pamięć
    //------------------------------------------------------------
    template <class T> matrix<T>::~matrix()
    {
        delete [] A;
    }

    // Zwraca wartość elementu
    // r - numer wiersza
    // c - numer kolumny
    //-------------------------
    template <class T> T matrix<T>::getv(int r, int c)
    {
        return A[r * n = c];
    }

    // Ustawia wartość elementu
    // r - numer wiersza
    // c - numer kolumny
    // v - wartość elementu
    //-------------------------
    template <class T> void matrix<T>::setv(int r, int c, T v)
    {
        A[r * n + c] = v;
    }

    // Ustawia przekatną na 1
    //-----------------------
    template <class T> void matrix<T>::set1()
    {
        int i;
        for(i = 0; i < n; i++) setv(i,i,1);
    }

    // Oblicza wyznacznik macierzy
    template <class T> T matrix<T>::det()
    {
        matrix<T> L(n);
        matrix<T> U(n);
        L.set1();
        int i,j,k;
        T sum;
        for(j = 0; j < n; j++)
        {
            for(i = 0; i <= j; i++)
            {
                sum = (T) 0;
                for(k = 0; k < i; k++) sum += L.getv(i,k) * U.getv(k,j);
                U.setv(i,j,getv(i,j) - sum);
            }
            for(i = j + 1; i < n; i++)
            {
                sum = (T) 0;
                for(k = 0; k < j; k++) sum += L.getv(i,k) * U.getv(k,j);
                if(U.getv(j,j)) L.setv(i,j,(getv(i,j) - sum) / U.getv(j,j));
            }
        }
        sum = (T) 1;
        for(i = 0; i < n; i++) sum *= U.getv(i,i);
        return sum;
    }
#endif // _matrix_class

// Program główny
//---------------

int main()
{

    // Wczytujemy dane dla macierzy

    setlocale(LC_ALL,"");

    cout << setprecision(2) << fixed;

    cout << "Wpisz dane dla macierzy:" << endl << endl;

    matrix<double> A;

    cout << endl
         << "Wyznacznik macierzy metodą rozkładu LU Doolitle'a" << endl
         << "-------------------------------------------------" << endl << endl;
    for(int i = 0; i < A.n; i++)
    {
        for(int j = 0; j < A.n; j++) cout << setw(8) << A.getv(i,j);
        cout << endl;
    }

    cout << endl << "det A = " << A.det();

    cout << endl << endl;

    return 0;
}

Wynik

Wpisz dane dla macierzy:

5
5 3 7 4 2
9 2 2 1 1
3 6 2 8 9
9 4 -2 -1 -3
0 5 3 -6 -11

Wyznacznik macierzy metodą rozkładu LU Doolitle'a
-------------------------------------------------

    5.00    3.00    7.00    4.00    2.00
    9.00    2.00    2.00    1.00    1.00
    3.00    6.00    2.00    8.00    9.00
    9.00    4.00   -2.00   -1.00   -3.00
    0.00    5.00    3.00   -6.00 -11.00

det A = -9368.00

Jeśli przyjrzysz się dokładnie wzorom Doolitle'a, to powinieneś zauważyć, że wyrazy macierzy A są wykorzystywane tylko jeden raz w każdym obiegu wyznaczania wyrazów macierzy L i U. Pozwala to umieścić macierze L i U w macierzy A, przez co algorytm działa w miejscu, tzn. nie zużywa dodatkowej pamięci. Jest to możliwe, ponieważ przekątna główna macierzy dolnej L nie jest używana przez algorytm. Zatem obie macierze L i U wypełnią dokładnie macierz A.

Poniższy program wykorzystuje tę koncepcję do wyliczenia wyznacznika. W przypadku próby dzielenia przez zero ustawiany jest znacznik error na true w macierzy A. Pozwala to programowi sprawdzić, czy obliczony wynik jest poprawny.

Przykładowe dane ( przekopiuj je do schowka i wklej w programie ):

Macierz:

Dane wejściowe: 5 5 3 7 4 2 9 2 2 1 1 3 6 2 8 9 9 4 -2 -1 -3 0 5 3 -6 -11

Przykładowy program w języku C++

// Macierze
// (C)2019 mgr Jerzy Wałaszek
// Metody numeryczne
//---------------------------

#include <iostream>
#include <iomanip>

using namespace std;

// Definicja klasy, może być w pliku nagłówkowym
//----------------------------------------------

#ifndef _matrix_class
    #define _matrix_class

    template <class T> class matrix
    {
        private:
            T * A;                   // Wyrazy macierzy
            T eps;
        public:
            int n;                   // Stopień macierzy
            bool error;              // Błąd obliczania wyznacznika
            matrix();                // Konstruktor
            ~matrix();               // Destruktor
            T  getv(int r, int c);
            void setv(int r, int c, T v);
            T det();                 // Oblicza wyznacznik
    };

    // Konstruktor - tworzy macierz ze strumienia cin
    //-----------------------------------------------
    template <class T> matrix<T>::matrix()
    {
        error = false;
        eps = (T) 0.000000001;
        cin >> n;       // Stopień macierzy
        int nn = n * n;
        A = new T [nn]; // Rezerwujemy pamięć na macierz n x n

        for(int i = 0; i < nn; i++) cin >> A[i]; // Wczytujemy wyrazy macierzy
    }

    // Destruktor - usuwa macierz i zwraca zajętą przez nią pamięć
    //------------------------------------------------------------
    template <class T> matrix<T>::~matrix()
    {
        delete [] A;
    }

    // Zwraca wartość elementu
    // r - numer wiersza
    // c - numer kolumny
    //-------------------------
    template <class T> T matrix<T>::getv(int r, int c)
    {
        return A[r * n + c];
    }

    // Ustawia wartość elementu
    // r - numer wiersza
    // c - numer kolumny
    // v - wartość elementu
    //-------------------------
    template <class T> void matrix<T>::setv(int r, int c, T v)
    {
        A[r * n + c] = v;
    }

    // Oblicza wyznacznik macierzy
    template <class T> T matrix<T>::det()
    {
        int i,j,k;
        T sum;
        for(j = 0; j < n; j++)
        {
            if((getv(j,j) > - eps) && (getv(j,j) < eps))
            {
                error = true; // Błąd w obliczaniu wyznacznika
                return 0;
            }
            for(i = 0; i <= j; i++)
            {
                sum = (T) 0;
                for(k = 0; k < i; k++) sum += getv(i,k) * getv(k,j);
                setv(i,j,getv(i,j) - sum);
            }
            for(i = j + 1; i < n; i++)
            {
                sum = (T) 0;
                for(k = 0; k < j; k++) sum += getv(i,k) * getv(k,j);
                setv(i,j,(getv(i,j) - sum) / getv(j,j));
            }
        }
        sum = (T) 1;
        for(i = 0; i < n; i++) sum *= getv(i,i);
        return sum;
    }
#endif // _matrix_class

// Program główny
//---------------

int main()
{

    // Wczytujemy dane dla macierzy

    setlocale(LC_ALL,"");

    cout << setprecision(2) << fixed;

    cout << "Wpisz dane dla macierzy:" << endl << endl;

    matrix<double> A;

    cout << endl
         << "Wyznacznik macierzy metodą rozkładu LU Doolitle'a w miejscu" << endl
         << "-----------------------------------------------------------" << endl << endl;
    for(int i = 0; i < A.n; i++)
    {
        for(int j = 0; j < A.n; j++) cout << setw(8) << A.getv(i,j);
        cout << endl;
    }

    cout << endl << "det A = " << A.det();

    if(A.error) cout << " - BŁĄD W OBLICZENIACH !!!";

    cout << endl << endl;

    return 0;
}

Wynik

Wpisz dane dla macierzy:

5
5 3 7 4 2
9 2 2 1 1
3 6 2 8 9
9 4 -2 -1 -3
0 5 3 -6 -11

Wyznacznik macierzy metodą rozkładu LU Doolitle'a w miejscu
-----------------------------------------------------------

    5.00    3.00    7.00    4.00    2.00
    9.00    2.00    2.00    1.00    1.00
    3.00    6.00    2.00    8.00    9.00
    9.00    4.00   -2.00   -1.00   -3.00
    0.00    5.00    3.00   -6.00 -11.00

det A = -9368.00

Istnieje prostsze rozwiązanie problemu wyznaczania rozkładu LU macierzy. Polega ono na przechodzeniu przez kolejne elementy głównej przekątnej macierzy A od pierwszego do przedostatniego i dla każdego elementu A[k,k] wykonywaniu kolejno dwóch operacji:

Normalizacji kolumny – elementy w kolumnie pod elementem a_kk są dzielone przez ten element.
Modyfikacji podmacierzy, która zawiera wszystkie elementy macierzy A leżące poniżej i na prawo elementu przekątnej z pominięciem wiersza i kolumny zawierającego element przekątnej. Modyfikacja polega na odjęciu od każdego elementu a_ij tej podmacierzy iloczynu elementów a_ik × a_kj.

Dzięki takiemu podejściu algorytm rozkładu LU znacznie się upraszcza.

Poniższy program wykorzystuje ten algorytm do wyznaczenia rozkładu LU macierzy. Macierz jest wprowadzana tak samo, jak w programie poprzednim. Po dokonaniu rozkładu LU program oblicza wyznacznik macierzy i wyświetla wynik. Jeśli w trakcie pracy program napotka zerowy element na głównej przekątnej, to zakończy działanie z informacją o błędzie.

Przykładowe dane ( przekopiuj je do schowka i wklej w programie ):

Macierz:

Dane wejściowe: 5 5 3 7 4 2 9 2 2 1 1 3 6 2 8 9 9 4 -2 -1 -3 0 5 3 -6 -11

Przykładowy program w języku C++

// Macierze
// (C)2019 mgr Jerzy Wałaszek
// Metody numeryczne
//---------------------------

#include <iostream>
#include <iomanip>

using namespace std;

// Definicja klasy, może być w pliku nagłówkowym
//----------------------------------------------

#ifndef _matrix_class
    #define _matrix_class

    template <class T> class matrix
    {
        private:
            T * A;                   // Wyrazy macierzy
            T eps;
        public:
            int n;                   // Stopień macierzy
            bool error;              // Błąd obliczania wyznacznika
            matrix();                // Konstruktor
            ~matrix();               // Destruktor
            T  getv(int r, int c);
            void setv(int r, int c, T v);
            T det();                 // Oblicza wyznacznik
    };

    // Konstruktor - tworzy macierz ze strumienia cin
    //-----------------------------------------------
    template <class T> matrix<T>::matrix()
    {
        error = false;
        eps = (T) 0.000000001;
        cin >> n;       // Stopień macierzy
        int nn = n * n;
        A = new T [nn]; // Rezerwujemy pamięć na macierz n x n

        for(int i = 0; i < nn; i++) cin >> A[i]; // Wczytujemy wyrazy macierzy
    }

    // Destruktor - usuwa macierz i zwraca zajętą przez nią pamięć
    //------------------------------------------------------------
    template <class T> matrix<T>::~matrix()
    {
        delete [] A;
    }

    // Zwraca wartość elementu
    // r - numer wiersza
    // c - numer kolumny
    //-------------------------
    template <class T> T matrix<T>::getv(int r, int c)
    {
        return A[r * n + c];
    }

    // Ustawia wartość elementu
    // r - numer wiersza
    // c - numer kolumny
    // v - wartość elementu
    //-------------------------
    template <class T> void matrix<T>::setv(int r, int c, T v)
    {
        A[r * n + c] = v;
    }

    // Oblicza wyznacznik macierzy
    template <class T> T matrix<T>::det()
    {
        int i,j,k;
        T sum, akk;

        for(k = 0; k < n - 1; k++)
        {
            // Dzielnik
            akk = getv(k,k);
            if((akk > -eps) && (akk < eps))
            {
                error = true;
                return 0;
            }

            // Normalizujemy kolumnę
            for(i = k + 1; i < n; i++) setv(i,k,getv(i,k) / akk);

            // Modyfikujemy podmacierz
            for(i = k + 1; i < n; i++)
                for(j = k + 1; j < n; j++) setv(i,j,getv(i,j) - getv(i,k) *  getv(k,j));
        }

        // Obliczamy wyznacznik;
        sum = getv(0,0);
        for(k = 1; k < n; k++) sum *= getv(k,k);
        return sum;
    }
#endif // _matrix_class

// Program główny
//---------------

int main()
{

    // Wczytujemy dane dla macierzy

    setlocale(LC_ALL,"");

    cout << setprecision(2) << fixed;

    cout << "Wpisz dane dla macierzy:" << endl << endl;

    matrix<double> A;

    cout << endl
         << "Wyznacznik macierzy metodą normalizacji kolumn" << endl
         << "----------------------------------------------" << endl << endl;
    for(int i = 0; i < A.n; i++)
    {
        for(int j = 0; j < A.n; j++) cout << setw(8) << A.getv(i,j);
        cout << endl;
    }

    cout << endl << "det A = " << A.det();

    if(A.error) cout << " - BŁĄD W OBLICZENIACH !!!";

    cout << endl << endl;

    return 0;
}

Wynik

Wpisz dane dla macierzy:

5
5 3 7 4 2
9 2 2 1 1
3 6 2 8 9
9 4 -2 -1 -3
0 5 3 -6 -11

Wyznacznik macierzy metodą normalizacji kolumn
----------------------------------------------

    5.00    3.00    7.00    4.00    2.00
    9.00    2.00    2.00    1.00    1.00
    3.00    6.00    2.00    8.00    9.00
    9.00    4.00   -2.00   -1.00   -3.00
    0.00    5.00    3.00   -6.00 -11.00

det A = -9368.00

Na początek: podrozdziału strony

Podane powyżej metody rozkładu LU zawodzą w przypadku napotkania elementu zerowego przez który należy wykonać dzielenie. W takich przypadkach programy przerywały obliczenia. Tymczasem rozwiązanie istnieje i jest dosyć proste. Wystarczy w takim przypadku tak poprzestawiać wiersze rozkładanej macierzy, aby na jej przekątnej znalazły się wyrazy o największych modułach. Wprowadźmy pojęcie macierzy permutacji P. Jest to macierz jednostkowa o stopniu macierzy A, w której zamieniono miejscami wiersze lub kolumny. Jeśli pomnożymy macierz A przez tak zmienioną macierz permutacji P, to wynikowa macierz będzie miała przestawione wiersze lub kolumny tak samo, jak zostały one przestawione w macierzy P. Obowiązuje zasada: mnożenie P × A przestawia wiersze, A × P przestawia kolumny. Zobaczmy na przykładzie, jak to działa:

Mamy macierz A czwartego stopnia:

Tworzymy dla niej macierz permutacji P:

W takiej postaci permutacja nie zamienia miejscami wierszy/kolumn macierzy A:

Zamieńmy miejscami dwa środkowe wiersze ( lub dwie środkowe kolumny, co daje identyczny wynik ) w macierzy permutacji P:

Teraz mnożenia ( sprawdź to ) dadzą wyniki:

Pierwsze mnożenie dało w wyniku macierz A z zamienionymi dwoma środkowymi wierszami, a drugie mnożenie dało w wyniku macierz A z zamienionymi dwiema środkowymi kolumnami. Zatem macierz permutacji P ma za zadanie odpowiednio pozamieniać ze sobą miejscami wiersze lub kolumny. Taka zamiana nazywa się piwotem ( ang. pivoting ). Termin ten pochodzi z gry w kosza, gdzie oznacza obrót zawodnika z piłką.

Po co to wszystko? Otóż permutacje pozwolą nam tak poukładać wiersze ( lub kolumny ) macierzy, aby na jej przekątnej głównej znalazły się wyrazy o największym module. Jeśli macierz nie jest zdegenerowana, to jest to zawsze możliwe. Permutacje wykonujemy w trakcie wyznaczania rozkładu LU. Musimy jednak zapamiętywać ich liczbę, ponieważ każda zamiana wierszy/kolumn zmienia znak wyznacznika na przeciwny.

Praktycznie nie tworzy się macierzy P. Permutacje wierszy zapamiętujemy w wektorze wierszowym W, poprzez który odwołujemy się do wierszy macierzy A przy wyznaczaniu rozkładu LU. Podobnie postępowaliśmy w trakcie wyliczania wartości wyznacznika za pomocą rozwinięcia Laplace'a. Każdorazowo w pętli rozkładu LU wyznaczamy największy element w pierwszej kolumnie bieżącej podmacierzy i zamieniamy jego wiersz z pierwszym wierszem poprzez wektor kolumnowy. Algorytm jest następujący:

Algorytm rozkładu LU z piwotowaniem

Dane wejściowe:

ε	–	dokładność przyrównania do zera
n	–	stopień macierzy
A	–	macierz podlegająca rozkładowi LU

Dane wyjściowe

A	–	macierz zawierająca połączone macierze L i U
W	–	wektor wierszowy
md	–	mnożnik do korekcji wyznacznika

Zmienne pomocnicze

i,j,k	–	indeksy
maxw	–	numer wiersza z elementem maksymalnym
maxe	–	moduł elementu maksymalnego
akk	–	element macierzy A leżący na jej przekątnej po piwotowaniu

Lista kroków

K01:	md ← 1	mnożnik korekcyjny ustawiamy na 1
K02:	Dla i = 0,1,...,n – 1: wykonuj: W [ i ] ← i	ustawiamy wektor wierszowy
K03:	Dla k = 0,1,...,n – 1: wykonuj kroki K04...K13	rozpoczynamy rozkład LU
K04:	maxw ← k	wiersz z elementem maksymalnym
K05:	maxe ← \| A [ W [ k ], k ] \|	moduł elementu maksymalnego
K06:	Dla i = k + 1, k + 2,..., n – 1: wykonuj krok K07	szukamy elementu większego
K07:	Jeśli \| A [W [ i ], k ] \| > maxe, to maxw ← i maxe ← \| A [ W [ i ], k ] \|
K08:	Jeśli maxe ≤ ε, to zakończ z wynikiem false	macierz jest zdegenerowana
K09:	Jeśli maxw ≠ W [ k ], to md ← - md W [ k ] ↔ W [ maxw ]	znaleziono element maksymalny
K10:	akk ← A [ W [ k ], k ]	rozkład LU
K11:	Dla i = k + 1, k + 2,...,n – 1: wykonuj: A [ W [ i ], k ] ← A [ W [ i ], k ] / akk	normalizacja kolumny
K12:	Dla i = k + 1,k + 2,...,n – 1: wykonuj krok K13	modyfikujemy podmacierz
K13:	Dla j = k + 1,k + 2,...,n – 1: wykonuj: A [ W [ i ], j ] ← A [ W [ i ], j ] - A [ W [ i ], k ] × A [ W [ k ], j ]
K14:	Zakończ z wynikiem true

Przy obliczaniu wyznacznika do wierszy macierzy A należy odwoływać się poprzez wektor wierszowy W.

Podany algorytm posiada jeszcze jedną zaletę: redukuje błędy zaokrągleń, ponieważ gwarantuje dzielenie przez największy element.

Poniższy program wykorzystuje ten algorytm do wyliczania wartości wyznacznika. Jest to modyfikacja poprzedniego programu.

Macierz:

Dane wejściowe: 5 0 3 7 4 2 9 2 2 1 1 3 6 2 8 9 9 4 -2 -1 -3 0 5 3 -6 -11

Przykładowy program w języku C++

// Macierze
// (C)2019 mgr Jerzy Wałaszek
// Metody numeryczne
//---------------------------

#include <iostream>
#include <iomanip>

using namespace std;

// Definicja klasy, może być w pliku nagłówkowym
//----------------------------------------------

#ifndef _matrix_class
    #define _matrix_class

    template <class T> class matrix
    {
        private:
            T * A;                   // Wyrazy macierzy
            T eps;
        public:
            int n;                   // Stopień macierzy
            bool error;              // Błąd obliczania wyznacznika
            matrix();                // Konstruktor
            ~matrix();               // Destruktor
            T  getv(int r, int c);
            void setv(int r, int c, T v);
            T  absT(T x);
            T det();                 // Oblicza wyznacznik
    };

    // Konstruktor - tworzy macierz ze strumienia cin
    //-----------------------------------------------
    template <class T> matrix<T>::matrix()
    {
        error = false;
        eps = (T) 0.000000001;
        cin >> n;       // Stopień macierzy
        int nn = n * n;
        A = new T [nn]; // Rezerwujemy pamięć na macierz n x n

        for(int i = 0; i < nn; i++) cin >> A[i]; // Wczytujemy wyrazy macierzy
    }

    // Destruktor - usuwa macierz i zwraca zajętą przez nią pamięć
    //------------------------------------------------------------
    template <class T> matrix<T>::~matrix()
    {
        delete [] A;
    }

    // Zwraca wartość elementu
    // r - numer wiersza
    // c - numer kolumny
    //-------------------------
    template <class T> T matrix<T>::getv(int r, int c)
    {
        return A[r * n + c];
    }

    // Ustawia wartość elementu
    // r - numer wiersza
    // c - numer kolumny
    // v - wartość elementu
    //-------------------------
    template <class T> void matrix<T>::setv(int r, int c, T v)
    {
        A[r * n + c] = v;
    }

    template <class T> T matrix<T>::absT(T x)
    {
        if(x < 0) x = - x;
        return x;
    }

    // Oblicza wyznacznik macierzy
    template <class T> T matrix<T>::det()
    {
        int i, j, k, * W, md, maxw;
        T sum, akk, maxe;

        // Mnożnik korekcyjny dla wyznacznika
        md = 1;

        // Tworzymy wektor wierszowy
        W = new int [n];

        // Inicjujemy wektor wierszowy numerami wierszy
        for(i = 0; i < n; i++) W[i] = i;

        // Rozkład LU
        for(k = 0; k < n; k++)
        {
            maxw = k; // numer wiersza z elementem max
            maxe = absT(getv(W[k],k));

            // Szukamy elementu max
            for(i = k + 1; i < n; i++)
                if(absT(getv(W[i],k)) > maxe)
                {
                    maxw = i;
                    maxe = absT(getv(W[i],k));
                }
            // Jeśli max jest zero, to macierz jest zdegenerowana
            if(maxe <= eps)
            {
                error = false;
                delete [] W;	
                return 0;
            }

            // Jeśli max nie leży na przekątnej, zamieniamy wiersze
            if(maxw != k)
            {
                md = -md; // Modyfikujemy mnożnik
                swap(W[k], W[maxw]);
            }

            // Dzielnik
            akk = getv(W[k],k);

            // Normalizujemy kolumnę
            for(i = k + 1; i < n; i++) setv(W[i],k,getv(W[i],k) / akk);

            // Modyfikujemy podmacierz
            for(i = k + 1; i < n; i++)
                for(j = k + 1; j < n; j++) setv(W[i],j,getv(W[i],j) - getv(W[i],k) *  getv(W[k],j));
        }

        // Obliczamy wyznacznik;
        sum = md * getv(W[0],0);
        for(k = 1; k < n; k++) sum *= getv(W[k],k);
        delete [] W;
        return sum;
    }
#endif // _matrix_class

// Program główny
//---------------

int main()
{

    // Wczytujemy dane dla macierzy

    setlocale(LC_ALL,"");

    cout << setprecision(2) << fixed;

    cout << "Wpisz dane dla macierzy:" << endl << endl;

    matrix<double> A;

    cout << endl
         << "Wyznacznik metodą LU z częściowym wyborem elementu głównego" << endl
         << "-----------------------------------------------------------" << endl << endl;
    for(int i = 0; i < A.n; i++)
    {
        for(int j = 0; j < A.n; j++) cout << setw(8) << A.getv(i,j);
        cout << endl;
    }

    cout << endl << "det A = " << A.det();

    if(A.error) cout << " - BŁĄD W OBLICZENIACH !!!";

    cout << endl << endl;

    return 0;
}

Wynik

Wpisz dane dla macierzy:

5
0 3 7 4 2
9 2 2 1 1
3 6 2 8 9
9 4 -2 -1 -3
0 5 3 -6 -11

Wyznacznik metodą LU z częściowym wyborem elementu głównego
-----------------------------------------------------------

    0.00    3.00    7.00    4.00    2.00
    9.00    2.00    2.00    1.00    1.00
    3.00    6.00    2.00    8.00    9.00
    9.00    4.00   -2.00   -1.00   -3.00
    0.00    5.00    3.00   -6.00 -11.00

det A = -8598.00

Na początek: podrozdziału strony

Równania

Rozkład LU

Definicje

Rozkład LU Doolitle'a

Rozkład LU z częściowym wyborem elementu głównego

Algorytm rozkładu LU z piwotowaniem

Dane wejściowe:

Dane wyjściowe

Zmienne pomocnicze

Lista kroków