Gaussian elimination¶
We consider a linear system of size \(n\):
\[A_{ij} x_j = b_i,\]
or explicitly
\[\begin{split}\newcommand{\x}[2]{a_{#1,#2}}
\newcommand{\y}[1]{x_{#1}}
\newcommand{\z}[1]{b_{#1}}
\begin{bmatrix}
\x{ 0}{ 0} & \x{ 0}{ 1} & \cdots & \x{ 0}{i-1} & \x{ 0}{i } & \x{ 0}{i+1} & \cdots & \x{ 0}{n-2} & \x{ 0}{n-1} \\
\x{ 1}{ 0} & \x{ 1}{ 1} & \cdots & \x{ 1}{i-1} & \x{ 1}{i } & \x{ 1}{i+1} & \cdots & \x{ 1}{n-2} & \x{ 1}{n-1} \\
\vdots & \vdots & \ddots & \vdots & \vdots & \vdots & & \vdots & \vdots \\
\x{i-1}{ 0} & \x{i-1}{ 1} & \cdots & \x{i-1}{i-1} & \x{i-1}{i } & \x{i-1}{i+1} & \cdots & \x{i-1}{n-2} & \x{i-1}{n-1} \\
\x{i }{ 0} & \x{i }{ 1} & \cdots & \x{i }{i-1} & \x{i }{i } & \x{i }{i+1} & \cdots & \x{i }{n-2} & \x{i }{n-1} \\
\x{i+1}{ 0} & \x{i+1}{ 1} & \cdots & \x{i+1}{i-1} & \x{i+1}{i } & \x{i+1}{i+1} & \cdots & \x{i+1}{n-2} & \x{i+1}{n-1} \\
\vdots & \vdots & & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
\x{n-2}{ 0} & \x{n-2}{ 1} & \cdots & \x{n-2}{i-1} & \x{n-2}{i } & \x{n-2}{i+1} & \cdots & \x{n-2}{n-2} & \x{n-2}{n-1} \\
\x{n-1}{ 0} & \x{n-1}{ 1} & \cdots & \x{n-1}{i-1} & \x{n-1}{i } & \x{n-1}{i+1} & \cdots & \x{n-1}{n-2} & \x{n-1}{n-1} \\
\end{bmatrix}
\begin{bmatrix}
\y{ 0} \\
\y{ 1} \\
\vdots \\
\y{i-1} \\
\y{i } \\
\y{i+1} \\
\vdots \\
\y{n-2} \\
\y{n-1} \\
\end{bmatrix}
=
\begin{bmatrix}
\z{ 0} \\
\z{ 1} \\
\vdots \\
\z{i-1} \\
\z{i } \\
\z{i+1} \\
\vdots \\
\z{n-2} \\
\z{n-1} \\
\end{bmatrix}.\end{split}\]
For notational simplicity, we consider the following augmented matrix:
\[\begin{split}\begin{bmatrix}
\x{ 0}{ 0} & \x{ 0}{ 1} & \cdots & \x{ 0}{i-1} & \x{ 0}{i } & \x{ 0}{i+1} & \cdots & \x{ 0}{n-2} & \x{ 0}{n-1} & \z{ 0} \\
\x{ 1}{ 0} & \x{ 1}{ 1} & \cdots & \x{ 1}{i-1} & \x{ 1}{i } & \x{ 1}{i+1} & \cdots & \x{ 1}{n-2} & \x{ 1}{n-1} & \z{ 1} \\
\vdots & \vdots & \ddots & \vdots & \vdots & \vdots & & \vdots & \vdots & \vdots \\
\x{i-1}{ 0} & \x{i-1}{ 1} & \cdots & \x{i-1}{i-1} & \x{i-1}{i } & \x{i-1}{i+1} & \cdots & \x{i-1}{n-2} & \x{i-1}{n-1} & \z{i-1} \\
\x{i }{ 0} & \x{i }{ 1} & \cdots & \x{i }{i-1} & \x{i }{i } & \x{i }{i+1} & \cdots & \x{i }{n-2} & \x{i }{n-1} & \z{i } \\
\x{i+1}{ 0} & \x{i+1}{ 1} & \cdots & \x{i+1}{i-1} & \x{i+1}{i } & \x{i+1}{i+1} & \cdots & \x{i+1}{n-2} & \x{i+1}{n-1} & \z{i+1} \\
\vdots & \vdots & & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\
\x{n-2}{ 0} & \x{n-2}{ 1} & \cdots & \x{n-2}{i-1} & \x{n-2}{i } & \x{n-2}{i+1} & \cdots & \x{n-2}{n-2} & \x{n-2}{n-1} & \z{n-2} \\
\x{n-1}{ 0} & \x{n-1}{ 1} & \cdots & \x{n-1}{i-1} & \x{n-1}{i } & \x{n-1}{i+1} & \cdots & \x{n-1}{n-2} & \x{n-1}{n-1} & \z{n-1} \\
\end{bmatrix}.\end{split}\]
Forward sweep¶
After forward-sweeping up to the \(i - 1\)-th row, we have
\[\begin{split}\begin{bmatrix}
\vdots & \vdots & \ddots & \vdots & \vdots & \vdots & & \vdots & \vdots & \vdots \\
0 & 0 & \cdots & 1 & a_{i-1,i } & a_{i-1,i+1} & \cdots & a_{i-1,n-2} & a_{i-1,n-1} & b_{i-1} \\
0 & 0 & \cdots & 0 & a_{i ,i } & a_{i ,i+1} & \cdots & a_{i ,n-2} & a_{i ,n-1} & b_{i } \\
0 & 0 & \cdots & 0 & a_{i+1,i } & a_{i+1,i+1} & \cdots & a_{i+1,n-2} & a_{i+1,n-1} & b_{i+1} \\
\vdots & \vdots & & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\
\end{bmatrix}.\end{split}\]
The \(i\)-th row is normalized (divided by \(\x{i}{i}\)) to obtain
\[\begin{split}\begin{bmatrix}
\vdots & \vdots & \ddots & \vdots & \vdots & \vdots & & \vdots & \vdots & \vdots \\
0 & 0 & \cdots & 1 & a_{i-1,i } & a_{i-1,i+1} & \cdots & a_{i-1,n-2} & a_{i-1,n-1} & b_{i-1} \\
0 & 0 & \cdots & 0 & 1 & a_{i ,i+1} := a_{i ,i+1} / a_{i ,i } & \cdots & a_{i ,n-2} := a_{i ,n-2} / a_{i ,i } & a_{i ,n-1} := a_{i ,n-1} / a_{i ,i } & b_{i } := b_{i } / a_{i ,i } \\
0 & 0 & \cdots & 0 & a_{i+1,i } & a_{i+1,i+1} & \cdots & a_{i+1,n-2} & a_{i+1,n-1} & b_{i+1} \\
\vdots & \vdots & & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\
\end{bmatrix}.\end{split}\]
src/gaussian_elimination.c¶
13const double f = 1. / a[i * nitems + i];
14for (size_t j = i; j < nitems; j++) {
15 a[i * nitems + j] *= f;
16}
17x[i] *= f;
The new \(i\)-th row is utilized to eliminate the \(i\)-th column:
\[\begin{split}\begin{bmatrix}
\vdots & \vdots & \ddots & \vdots & \vdots & \vdots & & \vdots & \vdots & \vdots \\
0 & 0 & \cdots & 1 & a_{i-1,i } & a_{i-1,i+1} & \cdots & a_{i-1,n-2} & a_{i-1,n-1} & b_{i-1} \\
0 & 0 & \cdots & 0 & 1 & a_{i ,i+1} & \cdots & a_{i ,n-2} & a_{i ,n-1} & b_{i } \\
0 & 0 & \cdots & 0 & 0 & a_{i+1,i+1} := a_{i+1,i+1} - a_{i+1,i }a_{i ,i+1} & \cdots & a_{i+1,n-2} := a_{i+1,n-2} - a_{i+1,i }a_{i ,n-2} & a_{i+1,n-1} := a_{i+1,n-1} - a_{i+1,i }a_{i ,n-1} & b_{i+1} := b_{i+1} - a_{i+1,i }b_{i } \\
\vdots & \vdots & & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\
\end{bmatrix}.\end{split}\]
src/gaussian_elimination.c¶
19for (size_t ii = i + 1; ii < nitems; ii++) {
20 const double f = a[ii * nitems + i];
21 for (size_t j = i + 1; j < nitems; j++) {
22 a[ii * nitems + j] -= f * a[i * nitems + j];
23 }
24 x[ii] -= f * x[i];
25}
Backward substitution¶
After the forward sweep, we have
\[\begin{split}\begin{bmatrix}
1 & \x{ 0}{ 1} & \cdots & \x{ 0}{i-1} & \x{ 0}{i } & \x{ 0}{i+1} & \cdots & \x{ 0}{n-2} & \x{ 0}{n-1} & \z{ 0} \\
0 & 1 & \cdots & \x{ 1}{i-1} & \x{ 1}{i } & \x{ 1}{i+1} & \cdots & \x{ 1}{n-2} & \x{ 1}{n-1} & \z{ 1} \\
\vdots & \vdots & \ddots & \vdots & \vdots & \vdots & & \vdots & \vdots & \vdots \\
0 & 0 & \cdots & 1 & \x{i-1}{i } & \x{i-1}{i+1} & \cdots & \x{i-1}{n-2} & \x{i-1}{n-1} & \z{i-1} \\
0 & 0 & \cdots & 0 & 1 & \x{i }{i+1} & \cdots & \x{i }{n-2} & \x{i }{n-1} & \z{i } \\
0 & 0 & \cdots & 0 & 0 & 1 & \cdots & \x{i+1}{n-2} & \x{i+1}{n-1} & \z{i+1} \\
\vdots & \vdots & & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\
0 & 0 & \cdots & 0 & 0 & 0 & \cdots & 1 & \x{n-2}{n-1} & \z{n-2} \\
0 & 0 & \cdots & 0 & 0 & 0 & \cdots & 0 & 1 & \y{n-1} \\
\end{bmatrix}.\end{split}\]
Starting from the bottom (\(n - 1\)-th) row, the upper-triangular elements are eliminated.
After backward-substituting up to the \(i + 1\)-th row, we have
\[\begin{split}\begin{bmatrix}
\vdots & \vdots & \ddots & \vdots & \vdots & \vdots & & \vdots & \vdots & \vdots \\
0 & 0 & \cdots & 1 & \x{i-1}{i } & \x{i-1}{i+1} & \cdots & \x{i-1}{n-2} & \x{i-1}{n-1} & \z{i-1} \\
0 & 0 & \cdots & 0 & 1 & \x{i }{i+1} & \cdots & \x{i }{n-2} & \x{i }{n-1} & \z{i } \\
0 & 0 & \cdots & 0 & 0 & 1 & \cdots & \x{i+1}{n-2} & \x{i+1}{n-1} & \y{i+1} \\
\vdots & \vdots & & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\
\end{bmatrix}.\end{split}\]
src/gaussian_elimination.c¶
28for (size_t i = nitems - 1; ; i--) {
29 for (size_t j = i + 1; j < nitems; j++) {
30 x[i] -= a[i * nitems + j] * x[j];
31 }
32 if (0 == i) {
33 break;
34 }
35}