Algorytmy i Struktury Danych

Rozważmy macierz kwadratową \(A\) o wymiarze \({n\times n}\). Załóżmy, iż jest ona iloczynem dwóch macierzy trójkątnych \(L\) i \(U\):

\[
A_{n\times n} =
\begin{bmatrix}
a_{1,1} & a_{1,2} & a_{1,3} & \dots & a_{1,n-1} & a_{1,n} \\
a_{2,1} & a_{2,2} & a_{2,3} & \dots & a_{2,n-1} & a_{2,n} \\
a_{3,1} & a_{3,2} & a_{3,3} & \dots & a_{3,n-1} & a_{3,n} \\
\dots & \dots & \dots & \dots & \dots & \dots \\
a_{n-1,1} & a_{n-1,2} & a_{n-1,3} & \dots & a_{n-1,n-1} & a_{n-1,n} \\
a_{n,1} & a_{n,2} & a_{n,3} & \dots & a_{n,n-1} & a_{n,n}
\end{bmatrix}
= L_{n\times n} \times U_{n\times n} =
\\ \\
=
\begin{bmatrix}
1 & 0 & 0 & \dots & 0 & 0 \\
l_{2,1} & 1 & 0 & \dots & 0 & 0 \\
l_{3,1} & l_{3,2} & 1 & \dots & 0 & 0 \\
\dots & \dots & \dots & \dots & \dots & \dots \\
l_{n-1,1} & l_{n-1,2} & l_{n-1,3} & \dots & 1 & 0 \\
l_{n,1} & l_{n,2} & l_{n,3} & \dots & l_{n,n-1} & 1
\end{bmatrix}
\times
\begin{bmatrix}
u_{1,1} & u_{1,2} & u_{1,3} & \dots & u_{1,n-1} & u_{1,n} \\
0 & u_{2,2} & u_{2,3} & \dots & u_{2,n-1} & u_{2,n} \\
0 & 0 & u_{3,3} & \dots & u_{3,n-1} & u_{3,n} \\
\dots & \dots & \dots & \dots & \dots & \dots \\
0 & 0 & 0 & \dots & u_{n-1,n-1} & u_{n-1,n} \\
0 & 0 & 0 & \dots & 0 & u_{n,n}
\end{bmatrix}
\]

Wykorzystując wzory na mnożenie macierzy otrzymujemy dla każdego elementu \(a_{i,j}\) macierzy \(A\) następujące równanie:

\[
A_{n\times n} = L_{n\times n} \times U_{n\times n}
\\ \\
a_{i,j} = \sum_{k=1}^{n} l_{i,k} \cdot u_{k,j}
\]

Ponieważ macierze \(L\) i \(U\) są macierzami trójkątnymi, to część ich elementów jest równa zero. Otrzymujemy stąd dwa przypadki w zależności od indeksów \(i\) oraz \(j\):

\(i \le j\), czyli element macierzy \(A\): \(a_{i,j}\) leży na głównej przekątnej lub ponad nią. W takim przypadku ostatnim, niezerowym elementem w wierszu \({i{-}\text{tym}}\) macierzy \(L\) jest element \(l_{i,i}\). Skoro tak, to podany wzór możemy zredukować do postaci:

\[
a_{i,j} = \sum\limits_{k=1}^{i} l_{i,k} \cdot u_{k,j}
\]

Dodatkowo z postaci macierzy \(L\) wiemy, iż element \(l_{i,i}\), leżący na głównej przekątnej, jest zawsze równy 1. Stąd:

\[
a_{i,j} = \left( \sum\limits_{k=1}^{i-1} l_{i,k} \cdot u_{k,j} \right) + u_{i,j}
\]

Wzór ten pozwala nam wyliczyć element \(u_{i,j}\) o ile są znane poprzednie elementy macierzy \(L\) i \(U\):

\[
u_{i,j} = a_{i,j} - \sum\limits_{k=1}^{i-1} l_{i,k} \cdot u_{k,j}
\]

\(i > j\), czyli element macierzy \(A\): \(a_{i,j}\) leży pod główną przekątną. W tym przypadku ostatnim, niezerowym elementem w \({j{-}\text{tej}}\) kolumnie macierzy \(U\) jest \(u_{j,j}\). Zatem wzór wyjściowy redukuje się do postaci:

\[
a_{i,j} = \sum\limits_{k=1}^{j} l_{i,k} \cdot u_{k,j}
\\ \\
a_{i,j} = \left( \sum\limits_{k=1}^{j-1} l_{i,k} \cdot u_{k,j} \right) + l_{i,j} \cdot u_{j,j}
\]

Otrzymany wzór pozwala wyliczyć element \(l_{i,j}\) o ile znane są poprzednie elementy macierzy \(L\) i \(U\):

\[
l_{i,j} = \frac{1}{u_{j,j}} \cdot \left( a_{i,j} - \sum\limits_{k=1}^{i-1} l_{i,k} \cdot u_{k,j} \right)
\]

W ten sposób doszliśmy do algorytmu Doolitle'a wyznaczania rozkładu macierzy \(A\) na iloczyn dwóch macierzy \(L\) i \(U\) (ang. LU matrix decomposition). Zasada działania tego algorytmu jest następująca:

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

Dla \(j = 1, 2, \dots, n\): wykonujemy kolejno dwie poniższe operacje:

Dla \(i = 1, 2, \dots, j\): wykorzystujemy wzór z punktu 1) do wyznaczenia elementów \(u_{i,j}\).

Dla \(i = j+1, j+2, \dots, n\): wykorzystujemy wzór z punktu 2) do wyznaczenia elementów \(l_{i,j}\).

Przykład:

W celu zobrazowania pracy algorytmu Doolitle'a wyznaczmy wzory na rozkład LU macierzy \(A_{3\times3}\):

\[
A_{3\times 3} =
\begin{bmatrix}
a_{1,1} & a_{1,2} & a_{1,3} \\
a_{2,1} & a_{2,2} & a_{2,3} \\
a_{3,1} & a_{3,2} & a_{3,3}
\end{bmatrix}
= L_{3\times 3} \times U_{3\times 3} =
\\ \\
=
\begin{bmatrix}
1 & 0 & 0 \\
l_{2,1} & 1 & 0 \\
l_{3,1} & l_{3,2} & 1
\end{bmatrix}
\times
\begin{bmatrix}
u_{1,1} & u_{1,2} & u_{1,3} \\
0 & u_{2,2} & u_{2,3} \\
0 & 0 & u_{3,3}
\end{bmatrix}
\]

\(j = 1\)

\(i = 1, i = j\), zatem wyznaczamy \(u_{1,1}\):

\[
u_{1,1} = a_{1,1} - \sum\limits_{k=1}^{0} l_{1,k} \cdot u_{k,1} = a_{1,1}
\]

\(i = 2, i > j\), wyznaczamy \(l_{2,1}\):

\[
l_{2,1} = \frac{1}{u_{1,1}} \cdot \left( a_{2,1} - \sum\limits_{k=1}^{0} l_{2,k} \cdot u_{k,1} \right) = \frac{1}{u_{1,1}} \cdot a_{2,1}
\]

\(i = 3, i > j\), wyznaczamy \(l_{3,1}\):

\[
l_{3,1} = \frac{1}{u_{1,1}} \cdot \left( a_{3,1} - \sum\limits_{k=1}^{0} l_{3,k} \cdot u_{k,1} \right) = \frac{1}{u_{1,1}} \cdot a_{3,1}
\]

\(j = 2\)

\(i = 1, i < j\), wyznaczamy \(u_{1,2}\):

\[
u_{1,2} = a_{1,2} - \sum\limits_{k=1}^{0} l_{1,k} \cdot u_{k,2} = a_{1,2}
\]

\(i = 2, i = j\), wyznaczamy \(u_{2,2}\):

\[
u_{2,2} = a_{2,2} - \sum\limits_{k=1}^{1} l_{2,k} \cdot u_{k,2} = a_{2,2} - l_{2,1} \cdot u_{1,2}
\]

\(i = 3, i > j\), wyznaczamy \(l_{3,2}\):

\[
l_{3,2} = \frac{1}{u_{2,2}} \cdot \left( a_{3,2} - \sum\limits_{k=1}^{1} l_{3,k} \cdot u_{k,2} \right) = \frac{1}{u_{2,2}} \cdot \left( a_{3,2} - l_{3,1} \cdot u_{1,2} \right)
\]

\(j = 3\)

\(i = 1, i < j\), wyznaczamy \(u_{1,3}\):

\[
u_{1,3} = a_{1,3} - \sum\limits_{k=1}^{0} l_{1,k} \cdot u_{k,3} = a_{1,3}
\]

\(i = 2, i < j\), wyznaczamy \(u_{2, 3}\):

\[
u_{2,3} = a_{2,3} - \sum\limits_{k=1}^{1} l_{2,k} \cdot u_{k,3} = a_{2,3} - l_{2,1} \cdot u_{1,3}
\]

\(i = 3, i = j\), wyznaczamy \(u_{3,3}\):

\[
u_{3,3} = a_{3,3} - \sum\limits_{k=1}^{2} l_{3,k} \cdot u_{k,3} = a_{3,3} - l_{3,1} \cdot u_{1,3} - l_{3,2} \cdot u_{2,3}
\]

Dla przykładowych danych liczbowych otrzymujemy ze wzorów:

\[
\begin{bmatrix}
1 & 4 & 7 \\
2 & 5 & 8 \\
3 & 6 & 9
\end{bmatrix}
=
\begin{bmatrix}
1 & 0 & 0 \\
2 & 1 & 0 \\
3 & 2 & 1
\end{bmatrix}
\times
\begin{bmatrix}
1 & 4 & 7 \\
0 & -3 & -6 \\
0 & 0 & 0
\end{bmatrix}
\]

W powyższym przykładzie zwróć uwagę na następujące fakty:

Potrzebne we wzorach elementy macierzy \(L\) i \(U\) są zawsze wyliczone wcześniej.
We wszystkich wzorach elementy macierzy \(A\) są używane tylko jeden raz.

Prowadzi to do wniosku, iż macierze \(L\) i \(U\) mogą zostać umieszczone w miejscu macierzy \(A\), czyli operacja będzie wykonana in situ (łac. w miejscu) i zaoszczędzimy pamięć. W macierzy \(L\) istotna jest tylko część leżąca pod główną przekątną. Elementy przekątnej są zawsze równe 1, a elementy leżące nad przekątną są równe zero i nie musimy ich wcale zapamiętywać. Z drugiej strony w macierzy \(U\) istotna staje się tylko główna przekątna oraz elementy leżące nad nią – pozostałe elementy są równe zero i również nie musimy ich zapamiętywać. Zatem algorytm Doolitle'a może przekształcać macierz \(A\) na połączone macierze \(L\) i \(U\):

\[
\begin{bmatrix}
a_{1,1} & a_{1,2} & a_{1,3} & \dots & a_{1,n-1} & a_{1,n} \\
a_{2,1} & a_{2,2} & a_{2,3} & \dots & a_{2,n-1} & a_{2,n} \\
a_{3,1} & a_{3,2} & a_{3,3} & \dots & a_{3,n-1} & a_{3,n} \\
\dots & \dots & \dots & \dots & \dots & \dots \\
a_{n-1,1} & a_{n-1,2} & a_{n-1,3} & \dots & a_{n-1,n-1} & a_{n-1,n} \\
a_{n,1} & a_{n,2} & a_{n,3} & \dots & a_{n,n-1} & a_{n,n}
\end{bmatrix}
\rightarrow
\\ \\
\rightarrow
\begin{bmatrix}
u_{1,1} & u_{1,2} & u_{1,3} & \dots & u_{1,n-1} & u_{1,n} \\
l_{2,1} & u_{2,2} & u_{2,3} & \dots & u_{2,n-1} & u_{2,n} \\
l_{3,1} & l_{3,2} & u_{3,3} & \dots & u_{3,n-1} & u_{3,n} \\
\dots & \dots & \dots & \dots & \dots & \dots \\
l_{n-1,1} & l_{n-1,2} & l_{n-1,3} & \dots & u_{n-1,n-1} & u_{n-1,n} \\
l_{n,1} & l_{n,2} & l_{n,3} & \dots & l_{n,n-1} & u_{n,n}
\end{bmatrix}
\]

Algorytm Doolitle'a rozkładu LU w miejscu

Wejście:

\(n\) : stopień macierzy; \({n\in\mathbb{N}}\).
\(A\) : macierz stopnia \(n\); \({elementy \; A\in\mathbb{R}}\).

Wyjście:

\(r = \text{true}\), macierz \(A\) zawiera połączone macierze \(L\) i \(U\).
\(r = \text{false}\), błąd dzielenia przez zero.

Elementy pomocnicze:

\(i, j, k\) : indeksy; i, j, k ∈ N.
\(s\) : suma iloczynów; s ∈ R.
\(\varepsilon\) : otoczenie zera; \({\varepsilon = 10^{-12}, \varepsilon\in\mathbb{R}}\).

Lista kroków:

K01: r ← false ; zakładamy porażkę
K02: Dla j = 1, 2, …, n:
     wykonuj kroki K03…K11
K03:   Jeśli |A[j,j]| < ε, ; rozkład niemożliwy, 
       to zakończ ; gdyż dzielnik równy zero
K04:   Dla i = 1, 2, …, j: ; najpierw wyznaczamy w kolumnie wyrazy u
       wykonuj kroki K05…K07
K05:     s ← 0
K06:     Dla k = 1, 2, …, i-1: ; wyznaczamy sumę iloczynów wyrazów z L i U
         wykonuj s ← s+A[i,k]×A[k,j]
K07:     A[i,j] ← A[i,j]-s ; obliczamy U[i,j]
K08:   Dla i = j+1, j+2, …, n:
       wykonuj kroki K09…K11
K09:     s ← 0
K10:     Dla k = 1, 2, …, j-1: ; wyznaczamy sumę iloczynów wyrazów z L i U
         wykonuj s ← s+A[i,k]×A[k,j]
K11:     A[i,j] ← (A[i,j]-s):A[j,j] ; obliczamy L[i,j]
K12: r ← true ; zaznaczamy sukces
K13: Zakończ

Przykładowe programy

Uwaga:

Zanim uruchomisz program, przeczytaj wstęp do tego artykułu, w którym wyjaśniamy funkcje tych programów oraz sposób korzystania z nich.

Program odczytuje ze standardowego wejścia macierz \(A\) (najpierw liczba n określająca stopień macierzy, a następnie n² liczb dla kolejnych elementów), wyznacza jej rozkład LU za pomocą algorytmu Doolitle'a i wyświetla wynikową macierz \(A\).

Dane przykładowe (przekopiuj do schowka i wklej do okna konsoli):

3
1 4 7
2 5 8
3 6 9

Pascal

// Rozkład LU algorytmem Doolitle'a
// Data: 21.03.2010
// (C)2020 mgr Jerzy Wałaszek
//-----------------------------

program ludoolitle;

type
  NReal = array of double;  // typ tablicy wierszy
  MReal = array of NReal;   // typ tablicy wskaźników

const

  eps = 1e-12;

// Funkcja realizuje algorytm Doolitle'a rozkładu LU
//--------------------------------------------------
function Doolitle(n : integer; A : MReal) : boolean;

var
  i, j, k : integer;
  s     : double;

begin

  for j := 0 to n-1 do
  begin
    if abs(A[j][j]) < eps then Exit(false);
    for i := 0 to j do
    begin
      s := 0;
      for k := 0 to i-1 do
        s := s+A[i][k]*A[k][j];
      A[i][j] := A[i][j]-s;
    end;
    for i := j+1 to n-1 do
    begin
      s := 0;
      for k := 0 to j-1 do
        s := s+A[i][k]*A[k][j];
      A[i][j] := (A[i][j]-s)/A[j][j];
    end;
  end;
  Doolitle := true;
end;

// Program główny

var
  A     : MReal;
  n, i, j : integer;

begin

  // odczytujemy stopień macierzy A
  read(n);

  // tworzymy macierz A o odpowiednich rozmiarach
  SetLength(A, n);
  for i := 0 to n-1 do SetLength(A[i], n);

  // odczytujemy dane dla macierzy A
  for i := 0 to n-1 do
    for j := 0 to n-1 do read(A[i][j]);

  // wyznaczamy rozkład LU macierzy A
  if Doolitle(n, A) then
    for i := 0 to n-1 do
    begin
      for j := 0 to n-1 do
        write(A[i][j]:8:3, ' ');
      writeln;
    end
  else
    writeln('DZIELNIK ZERO');

  // usuwamy macierz z pamięci
  for i := 0 to n-1 do
    SetLength(A[i], 0);
  SetLength(A, 0);

end.

C++

// Rozkład LU algorytmem Doolitle'a
// Data: 21.03.2010
// (C)2020 mgr Jerzy Wałaszek
//-----------------------------

#include <iostream>
#include <iomanip>
#include <cmath>

using namespace std;

const double eps = 1e-12;

// Funkcja realizuje algorytm
// Doolitle'a rozkładu LU
//---------------------------
bool Doolitle(int n, 
              double ** A)
{
  int i, j, k;
  double s;
 
  for(j = 0; j < n; j++)
  {
    if(fabs(A[j][j]) < eps)
      return false;
    for(i = 0; i <= j; i++)
    {
      s = 0;
      for(k = 0; k < i; k++)
        s +=  A[i][k]*A[k][j];
      A[i][j] -= s;
    }
    for(i = j+1; i < n; i++)
    {
      s = 0;
      for(k = 0; k < j; k++)
        s += A[i][k]*A[k][j];
      A[i][j] = (A[i][j]-s)/A[j][j];
    }
  }
  return true;
}

// Program główny

int main()
{
  double **A;
  int n, i, j;
 
  cout << setprecision(3) << fixed;
 
  // odczytujemy stopień macierzy A
  cin >> n;

  // tworzymy macierz A o odpowiednich rozmiarach
  A = new double * [n];
  for(i = 0; i < n; i++)
    A[i] = new double[n];

  // odczytujemy dane dla macierzy A
  for(i = 0; i < n; i++)
    for(j = 0; j < n; j++)
      cin >> A[i][j];

  // wyznaczamy rozkład LU macierzy A
  if(Doolitle(n, A))
    for(i = 0; i < n; i++)
    {
      for(j = 0; j < n; j++)
        cout << setw(8) << A[i][j] << " ";
      cout << endl;
    }
  else
    cout << "DZIELNIK ZERO\n";

  // usuwamy macierz z pamięci
  for(i = 0; i < n; i++)
    delete [] A[i];
  delete [] A;
  return 0;
}

Basic

' Rozkład LU algorytmem Doolitle'a
' Data: 21.03.2010
' (C)2020 mgr Jerzy Wałaszek
'-----------------------------

Const eps As Double = 1e-12

' Funkcja realizuje algorytm
' Doolitle'a rozkładu LU
'---------------------------
Function doolitle(n As Integer, _
                  A As Double Ptr Ptr) _
                  As Integer
  Dim As Integer i, j, k
  Dim As Double s
 
  For j = 0 To n-1
    If Abs(A[j][j]) < eps Then
      Return 0
    For i = 0 To j
      s = 0
      For k = 0 To i-1
        s += A[i][k]*A[k][j]
      Next
      A[i][j] -= s
    Next
    For i = j+1 To n-1
      s = 0
      For k = 0 To j-1
        s += A[i][k]*A[k][j]
      Next
      A[i][j] = (A[i][j]-s)/A[j][j]
    Next
  Next
  doolitle = 1
End Function

Dim As Integer i, j, n
Dim A As Double Ptr Ptr

Open Cons For Input As #1

' odczytujemy stopień macierzy A
Input #1, n

' tworzymy macierz A o odpowiednich rozmiarach
A = New Double Ptr [n]
For i = 0 To n-1
  A[i] = New Double [n]
Next

' odczytujemy dane dla macierzy A
For i = 0 To n-1
  For j = 0 To n-1
    Input #1, A[i][j]
  Next
Next

Close #1

' wyznaczamy rozkład LU macierzy A
If doolitle(n, A) = 1 Then
  For i = 0 To n-1
    For j = 0 To n-1
      Print Using "####.### ";A[i][j];
    Next
    Print
  Next
Else
  Print "DZIELNIK ZERO"
End If

' usuwamy macierz z pamięci
For i = 0 To n-1
  Delete [] A[i]
Next
Delete [] A

End

Python (dodatek)

# Rozkład LU algorytmem Doolitle'a
# Data: 23.04.2024
# (C)2024 mgr Jerzy Wałaszek
#-----------------------------

EPS = 1e-12


# Funkcja realizuje algorytm
# Doolitle'a rozkładu LU
#---------------------------
def doolitle(a):
    n = len(a)
    for j in range(n):
        if abs(a[j][j]) < EPS:
            return False
        for i in range(j+1):
            s = 0
            for k in range(i):
                s += a[i][k]*a[k][j]
            a[i][j] -= s
        for i in range(j+1, n):
            s = 0
            for k in range(j):
                s += a[i][k]*a[k][j]
            a[i][j] = (a[i][j]-s)/a[j][j]
    return True


# Program główny

# odczytujemy stopień macierzy A
n = int(input())
# odczytujemy macierz A
a = []
for i in range(n):
    arr = input().split()
    arr = [float(arr[j]) for j in range(n)]
    a.append(arr)
# wyznaczamy rozkład LU macierzy A
if doolitle(a):
    for i in range(n):
        for j in range(n):
            print("%8.3f " % a[i][j], end="")
        print()
else:
    print("DZIELNIK ZERO")
# usuwamy macierz z pamięci
a = None

Wynik:

3
1 4 7
2 5 8
3 6 9
   1.000    4.000    7.000 
   2.000   -3.000   -6.000 
   3.000    2.000    0.000

do podrozdziału do strony

W literaturze informatycznej można znaleźć prostsze rozwiązanie 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[k, k]}\) 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[i, j]}\) tej podmacierzy iloczynu elementów \({A[i, k]\cdot A[k, j]}\).

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

Algorytm rozkładu LU w miejscu

Wejście:

\(n\) : stopień macierzy; \({n\in\mathbb{N}}\).
\(A\) : macierz stopnia \(n\); \({elementy \; A\in\mathbb{R}}\).

Wyjście:

\(r = \text{true}\), macierz \(A\) zawiera połączone macierze\(L\) i \(U\).
\(r = \text{false}\), błąd dzielenia przez zero.

Elementy pomocnicze:

\(i, j, k\) : indeksy; \({i, j, k\in\mathbb{N}}\).
\(\varepsilon\) : otoczenie zera; \({\varepsilon = 10^{-12}, \varepsilon\in\mathbb{R}}\).

Lista kroków:

K01: r ← false ; zakładamy porażkę
K02: Dla k = 1, 2, …, n-1:
     wykonuj kroki K03…K06
K03:   Jeśli |A[k, k]| < ε, ; element zerowy nie pozwoli
       to zakończ          ; wykonać dzielenia
K04:   Dla i = k+1, k+2, …, n: ; normalizujemy kolumnę
       wykonuj: A[i, j] ← A[i, k]:A[k, k]
K05:   Dla i = k+1, k+2, …, n:
       wykonuj krok K06
K06:   Dla j = k+1, k+2, …, n: ; modyfikujemy podmacierz
       wykonuj: A[i, j] ← A[i, j]-A[i, k]×A[k, j]
K07: r ← true ; sukces
K08: Zakończ

Przykładowe programy

Uwaga:

Zanim uruchomisz program, przeczytaj wstęp do tego artykułu, w którym wyjaśniamy funkcje tych programów oraz sposób korzystania z nich.

Pascal

// Rozkład LU
// Data: 30.04.2010
// (C)2020 mgr Jerzy Wałaszek
//-----------------------------

program ludecomp;

type
  NReal = array of double; // typ tablicy wierszy
  MReal = array of NReal;  // typ tablicy wskaźników

const
  eps = 1e-12;

// Funkcja realizuje algorytm rozkładu LU
//---------------------------------------
function lu(n : integer;
            A : MReal)
            : boolean;
var
  i, j, k : integer;
begin
  for k := 0 to n-2 do
  begin
    if abs(A[k][k]) < eps then
      Exit(false);
    for i := k+1 to n-1 do
      A[i][k] := A[i][k]/A[k][k];
    for i := k+1 to n-1 do
      for j := k+1 to n-1 do
        A[i][j] := A[i][j]-A[i][k]*A[k][j];
  end;
  lu := true;
end;

// Program główny

var
  A     : MReal;
  n, i, j : integer;

begin

  // odczytujemy stopień macierzy A
  read(n);

  // tworzymy macierz A o odpowiednich rozmiarach
  SetLength(A, n);
  for i := 0 to n-1 do
    SetLength(A[i], n);

  // odczytujemy dane dla macierzy A
  for i := 0 to n-1 do
    for j := 0 to n-1 do
      read(A[i][j]);

  // wyznaczamy rozkład LU macierzy A
  if lu(n, A) then
    for i := 0 to n-1 do
    begin
      for j := 0 to n-1 do
        write(A[i][j]:8:3, ' ');
      writeln;
    end
  else
    writeln('DZIELNIK ZERO');

  // usuwamy macierz z pamięci
  for i := 0 to n-1 do
    SetLength(A[i], 0);
  SetLength(A, 0);

end.

C++

// Rozkład LU
// Data: 30.04.2010
// (C)2020 mgr Jerzy Wałaszek
//-----------------------------

#include <iostream>
#include <iomanip>
#include <cmath>

using namespace std;

const double eps = 1e-12;

// Funkcja realizuje algorytm rozkładu LU
//---------------------------------------
bool lu(int n, double ** A)
{
  int i, j, k;

  for(k = 0; k < n-1; k++)
  {
    if(fabs(A[k][k]) < eps)
      return false;
    for(i = k+1; i < n; i++)
      A[i][k] /= A[k][k];
    for(i = k+1; i < n; i++)
      for(j = k+1; j < n; j++)
        A[i][j] -= A[i][k]*A[k][j];
  }
  return true;
}

// Program główny

int main()
{
  double **A;
  int n, i, j;
 
  cout << setprecision(3) << fixed;
 
  // odczytujemy stopień macierzy A
  cin >> n;

  // tworzymy macierz A
  // o odpowiednich rozmiarach
  A = new double * [n];
  for(i = 0; i < n; i++)
    A[i] = new double[n];

  // odczytujemy dane dla macierzy A
  for(i = 0; i < n; i++)
    for(j = 0; j < n; j++)
      cin >> A[i][j];

  // wyznaczamy rozkład LU macierzy A
  if(lu(n, A))
    for(i = 0; i < n; i++)
    {
      for(j = 0; j < n; j++)
        cout << setw(8) << A[i][j] << " ";
      cout << endl;
    }
  else
    cout << "DZIELNIK ZERO\n";

  // usuwamy macierz z pamięci
  for(i = 0; i < n; i++)
    delete [] A[i];
  delete [] A;

  return 0;
}

Basic

' Rozkład LU
' Data: 30.04.2010
' (C)2020 mgr Jerzy Wałaszek
'-----------------------------

Const eps As Double = 1e-12

' Funkcja realizuje algorytm rozkładu LU
'---------------------------------------
Function lu(n As Integer, _
            A As Double Ptr Ptr) _
            As Integer

  Dim As Integer i, j, k
 
  For k = 0 To n-2
    If Abs(A[k][k]) < eps Then _
      Return 0
    For i = k+1 To n-1
      A[i][k] /= A[k][k]
    Next
    For i = k+1 To n-1
      For j = k+1 To n-1
        A[i][j] -= A[i][k]*A[k][j]
      Next
    Next
  Next

  lu = 1

End Function

Dim A As Double Ptr Ptr
Dim As Integer i, j, n

Open Cons For Input As #1

' odczytujemy stopień macierzy A
Input #1, n

' tworzymy macierz A
' o odpowiednich rozmiarach
A = New Double Ptr[n]
For i = 0 To n-1
  A[i] = New Double[n]
Next

' odczytujemy dane dla macierzy A
For i = 0 To n-1
  For j = 0 To n-1
    Input #1, A[i][j]
  Next
Next

Close #1

' wyznaczamy rozkład LU macierzy A
If lu(n, A) = 1 Then
  For i = 0 To n-1
    For j = 0 To n-1
      Print Using "####.### ";A[i][j];
    Next
    Print
  Next
Else
    Print "DZIELNIK ZERO"
End If

' Usuwamy macierz A

For i = 0 To n-1
  Delete [] A[i]
Next
Delete [] A

End

Python (dodatek)

# Rozkład LU
# Data: 23.04.2024
# (C)2024 mgr Jerzy Wałaszek
#---------------------------

EPS = 1e-12


# Funkcja realizuje algorytm rozkładu LU
#---------------------------------------
def lu(a):
    n = len(a)
    for k in range(n-1):
        if abs(a[k][k]) < EPS:
            return False
        for i in range(k+1, n):
            a[i][k] /= a[k][k]
        for i in range(k+1, n):
            for j in range(k+1, n):
                a[i][j] -= a[i][k]*a[k][j]
    return True


# Program główny
#---------------

# odczytujemy stopień macierzy A
n = int(input())
# odczytujemy macierz A
a = []
for i in range(n):
    arr = input().split()
    arr = [float(arr[j]) for j in range(n)]
    a.append(arr)
# wyznaczamy rozkład LU macierzy A
if lu(a):
    for i in range(n):
        for j in range(n):
            print("%8.3f" % a[i][j], end=" ")
        print()
else:
    print("DZIELNIK ZERO")
# usuwamy macierz z pamięci
a = None

do podrozdziału do strony

Rozkład LU

Problem

Rozwiązanie nr 1

\(j = 1\)

\(j = 2\)

\(j = 3\)

Algorytm Doolitle'a rozkładu LU w miejscu

Wejście:

Wyjście:

Elementy pomocnicze:

Lista kroków:

Przykładowe programy

Rozwiązanie nr 2

Algorytm rozkładu LU w miejscu

Wejście:

Wyjście:

Elementy pomocnicze:

Lista kroków:

Przykładowe programy