C++ Implementation#
Here’s an implementation of Gaussian elimination in C++, with partial pivoting.
This uses a two-dimensional multidimensional array class, Array defined
in the header array.H. This allows us to do things like:
Array A{{{1, 1, 1},
{-1, 2, 0},
{2, 0, 1}}};
To initialize a 2-d array. We can then get the number of rows and
columns via A.nrows() and A.ncols(). We can also index the
array using A(i,j) to get the i-th row and j-th column.
Some details of this class can be found here: https://zingale.github.io/phy504/contiguous_array.html (Although the implementation here is slightly different).
Here’s the class header array.H:
#ifndef ARRAY_H
#define ARRAY_H
#include <vector>
#include <iostream>
#include <cassert>
// a contiguous 2-d array
// here the data is stored in row-major order in a 1-d memory space
// managed as a vector. We overload () to allow us to index this as
// a(irow, icol)
struct Array {
std::size_t _rows;
std::size_t _cols;
std::vector<double> _data;
Array (std::size_t rows, std::size_t cols, double val=0.0)
: _rows{rows},
_cols{cols},
_data(rows * cols, val)
{
// we do the asserts here after the initialization of _data
// in the initialization list, but if the size is zero, we'll
// abort here anyway.
assert (rows > 0);
assert (cols > 0);
}
// constructor that allows us to initialize array
// via an initializer list of rows
Array (std::vector<std::vector<double>> v)
: _rows{v.size()},
_cols{v[0].size()},
_data(_rows * _cols, 0.0)
{
int idx = 0;
for (std::size_t i = 0; i < _rows; ++i) {
assert (v[i].size() == _cols);
for (std::size_t j = 0; j < _cols; ++j) {
_data[idx] = v[i][j];
++idx;
}
}
}
// note the "const" after the argument list here -- this means
// that this can be called on a const Array
inline std::size_t ncols() const { return _cols;}
inline std::size_t nrows() const { return _rows;}
inline double& operator()(int row, int col) {
assert (row >= 0 && row < static_cast<int>(_rows));
assert (col >= 0 && col < static_cast<int>(_cols));
return _data[row*_cols + col];
}
inline const double& operator()(int row, int col) const {
assert (row >= 0 && row < static_cast<int>(_rows));
assert (col >= 0 && col < static_cast<int>(_cols));
return _data[row*_cols + col];
}
};
// the << operator is not part of the of the class, so it is not a member
std::ostream& operator<< (std::ostream& os, const Array& a) {
for (std::size_t row = 0; row < a.nrows(); ++row) {
for (std::size_t col = 0; col < a.ncols(); ++col) {
os << a(row, col) << " ";
}
os << std::endl;
}
return os;
}
#endif
A function to do a matrix-vector multiply is contained in matmul.H:
#ifndef MATMUL_H
#define MATMUL_H
#include <vector>
#include <array.H>
inline
std::vector<double> matmul(const Array& A, const std::vector<double>& x) {
assert (A.ncols() == static_cast<int>(x.size()));
std::vector<double> b(A.nrows(), 0.0);
for (int m = 0; m < A.nrows(); ++m) {
for (int n = 0; n < A.ncols(); ++n) {
b[m] += A(m, n) * x[n];
}
}
return b;
}
#endif
The Gaussian elimination function is in the header gauss.H:
#ifndef GAUSS_H
#define GAUSS_H
#include <array.H>
#include <vector>
#include <limits>
#include <iostream>
#include <iomanip>
inline
void print_Ab(const Array& A, const std::vector<double>& b) {
// printout the matrix A and vector b in a pretty fashion. We don't use the
// numpy print here, because we want to make them side by side
const int N = b.size();
for (int irow = 0; irow < N; ++irow) {
std::cout.precision(4);
std::cout << std::fixed;
if (irow == 0) {
std::cout << "/ ";
} else if (irow == N-1) {
std::cout << "\\ ";
} else {
std::cout << "| ";
}
for (int jcol = 0; jcol < N; ++jcol) {
std::cout << std::setw(8) << A(irow, jcol);
}
std::cout << " | ";
std::cout << std::setw(8) << b[irow];
if (irow == 0) {
std::cout << " \\ ";
} else if (irow == N-1) {
std::cout << " / ";
} else {
std::cout << " | ";
}
std::cout << std::endl;
}
std::cout << std::endl;
}
inline
std::vector<double> gauss_elim(Array& A, std::vector<double>& b,
const bool quiet=false) {
// perform gaussian elimination with pivoting, solving A x = b.
//
// A is an NxN matrix, x and b are an N-element vectors. Note: A and b are
// changed upon exit to be in upper triangular (row echelon) form
const int N = b.size();
// A is square, with each dimension of length N
assert(static_cast<int>(A.nrows()) == N && static_cast<int>(A.ncols()) == N);
// allocate the solution array
std::vector<double> x(N, 0.0);
// find the scale factors for each row -- this is used when pivoting
std::vector<double> scales(N, 0.0);
for (int irow = 0; irow < N; ++irow) {
double maxval = std::numeric_limits<double>::lowest();
for (int jcol = 0; jcol < N; ++jcol) {
maxval = std::max(maxval, std::abs(A(irow, jcol)));
}
scales[irow] = maxval;
}
// keep track of the number of times we swapped rows
int num_row_swap{};
if (!quiet) {
print_Ab(A, b);
}
// main loop over rows
for (int krow = 0; krow < N; ++krow) {
// find the pivot row based on the size of column k -- only consider the
// rows beyond the current row
int row_max{krow};
double col_max = std::numeric_limits<double>::lowest();
for (int kk = krow; kk < N; ++kk) {
double scaled_col_val = std::abs(A(kk, krow)) / scales[kk];
if (scaled_col_val > col_max) {
col_max = scaled_col_val;
row_max = kk;
}
}
// swap the row with the largest scaled element in the current column
// with the current row (pivot) -- do this with b too!
if (row_max != krow) {
// swap row krow with row row_max
for (int jcol = 0; jcol < N; ++jcol) {
double tmp_val = A(krow, jcol);
A(krow, jcol) = A(row_max, jcol);
A(row_max, jcol) = tmp_val;
}
// now swap the same elements in b
double tmp_b = b[krow];
b[krow] = b[row_max];
b[row_max] = tmp_b;
// finally swap the scales
double tmp_scale = scales[krow];
scales[krow] = scales[row_max];
scales[row_max] = tmp_scale;
if (!quiet) {
std::cout << "pivoted" << std::endl;
}
++num_row_swap;
}
// do the forward-elimination for all rows below the current
for (int irow = krow+1; irow < N; ++irow) {
double coeff = A(irow, krow) / A(krow, krow);
for (int jcol = krow+1; jcol < N; ++jcol) {
A(irow, jcol) += -A(krow, jcol) * coeff;
}
A(irow, krow) = 0.0;
b[irow] += -coeff * b[krow];
}
if (!quiet) {
print_Ab(A, b);
}
}
// back-substitution
// last solution is easy
x[N-1] = b[N-1] / A(N-1, N-1);
for (int irow = N-2; irow >= 0; --irow) {
double bsum = b[irow];
for (int jcol = irow+1; jcol < N; ++jcol) {
bsum += -A(irow, jcol) * x[jcol];
}
x[irow] = bsum / A(irow, irow);
}
return x;
}
#endif
Finally, a test driver gauss_test.cpp is here:
#include <iostream>
#include <array.H>
#include <gauss.H>
int main() {
Array A{{{1, 1, 1},
{-1, 2, 0},
{2, 0, 1}}};
std::vector<double> b{6, 3, 5};
auto x = gauss_elim(A, b);
for (auto e : x) {
std::cout << e << std::endl;
}
std::cout << std::endl;
std::cout << std::endl;
// next example
Array A2{{{0, 0, 0, 4},
{0, 0, 3, 0},
{5, 6, 7, 8},
{0, 4, 3, 2}}};
std::vector<double> b2{5, 4, 9, 1};
auto x2 = gauss_elim(A2, b2);
for (auto e : x2) {
std::cout << e << std::endl;
}
std::cout << std::endl;
std::cout << std::endl;
}
To run this, we would do:
g++ -I . --std=c++17 -o gauss_test gauss_test.cpp
./gauss_test
Try it#
Let’s test this out on the system we used when we introduced pivoting:
\[\begin{alignat*}{4}
\epsilon x &+ y &+ z &= 6 \\
-x &+ 2y & &= 3 \\
2x &+ &+ z &= 5
\end{alignat*}\]