# Gaussian elimination algorithm in C++

Given:

j-g-h-i=0
a+b-c-j=0
c+i-d-e=0
e+g-f=0


And known:

a=10
b=7
d=3
e=2
f=3
j=14


I want to solve this (or similar equations) with the Gaussian elimination algorithm, and wrote the following C++ code:

#include <string>
#include <vector>
#include <sstream>
#include <iostream>
using namespace std;
const int my_magic1 = 999;
const int my_magic2 = -999;

void printMatrix1(vector<vector<int>> &matrix1)
{
printf("Matrix1: (%d x %d)\n", matrix1.size(), matrix1[0].size());
for (int i = 0; i < matrix1.size(); i++)
{
for (int j = 0; j < matrix1[i].size(); j++)
{
cout << matrix1[i][j] << " ";
}
cout << endl;
}
cout << endl;
}

void printMatrix2(vector<int> &matrix2)
{
printf("Matrix2: (%d)\n", matrix2.size());
for (int i = 0; i < matrix2.size(); i++)
{
cout << "  " << (char)('a' + i) << " ";
}
cout << endl;
for (int i = 0; i < matrix2.size(); i++)
{
printf("% 3d ", matrix2[i]);
}
cout << endl;
cout << endl;
}

void solveAndPrint(vector<vector<int>> &matrix1, vector<int> &matrix2)
{
vector<vector<int>> result(matrix1);
vector<int> result2(matrix2);
for (int i = 0; i < result.size(); i++)
{
for (int j = 0; j < result[i].size(); j++)
{
if (result[i][j] != 0)
{
if (result2[j] != my_magic1)
{
result[i][j] = result[i][j] * result2[j];
}
else if (result[i][j] < 0)
{
result[i][j] = my_magic2;
}
else
{
result[i][j] = my_magic1;
}
}
}
}

printMatrix1(result);
printMatrix2(result2);

bool changed = true;
while (changed)
{
changed = false;
for (int i = 0; i < result.size(); i++)
{
int magic_n = 0;
int magic_i = 0;
bool isneg = false;
for (int j = 0; j < result[i].size(); j++)
{
if (result[i][j] == my_magic1 || result[i][j] == my_magic2)
{
magic_n++;
magic_i = j;
isneg = result[i][j] == my_magic2;
}
}
if (magic_n == 1)
{
int sum = 0;
for (int j = 0; j < result[i].size(); j++)
{
if (j != magic_i)
{
sum += result[i][j];
}
}
for (int k = 0; k < result.size(); k++)
{
if (result[k][magic_i] == my_magic1)
{
result[k][magic_i] = isneg ? sum : -sum;
}
else if (result[k][magic_i] == my_magic2)
{
result[k][magic_i] = isneg ? -sum : sum;
}
}
if (result2[magic_i] == my_magic1)
{
result2[magic_i] = isneg ? sum : -sum;
}
changed = true;
break;
}
}
}

printMatrix1(result);
printMatrix2(result2);
}

int main(int argc, char const *argv[])
{
vector<string> args({"0 0 0 0 0 0 -1 -1 -1 1",
"1 1 -1 0 0 0 0 0 0 -1",
"0 0 1 -1 -1 0 0 0 1 0",
"0 0 0 0 1 -1 1 0 0 0",
"0",
"10 7 - 3 2 3 - - - 14"});
vector<vector<int>> matrix1;
vector<int> matrix2;
for (int i = 1;; i++)
{
printf("Coefficients %d. row, or 0:\n", i);
string l;
// getline(cin, l);
l = args[i - 1];
cout << l << endl;

if (l == "0")
{
break;
}
stringstream ss;
ss << l;
string segment;
vector<int> v;
while (getline(ss, segment, ' '))
{
v.push_back(stoi(segment));
}
matrix1.push_back(v);
}
{
cout << "Known, or just - if unknown:" << endl;
string l;
// getline(cin, l);
l = args[args.size() - 1];
cout << l << endl;

stringstream ss;
ss << l;
string segment;
while (getline(ss, segment, ' '))
{
if (segment == "-")
{
matrix2.push_back(my_magic1);
}
else
{
matrix2.push_back(stoi(segment));
}
}
}
solveAndPrint(matrix1, matrix2);
return 0;
}


I would like to get a code review, and how I can avoid my_magic1 and my_magic2.

• vector<vector<T>> is an inefficient way to hold a matrix. If you still want to be able reorder whole rows by just swapping pointers, you can do that with pointers into one contiguous allocation so at least it's not so bad for cache access pattern. If you don't need that, just do 2D indexing (maybe with a wrapper class like boost::multi_array (example) to keep the [i][j] syntax instead of flat_vector[i*n + j]. Commented Nov 26, 2023 at 19:55
• @PeterCordes I wanted to solve this with C++(17) "standard" language features. So, the boost lib isn't an option here. Based on the reviews to date, I have do a refactoring, and the code can be found here. I'm still using vectors there, because sorting or swapping isn't needed (yet). Commented Nov 27, 2023 at 14:53
• You have that backwards. If sorting and swapping of whole rows is not needed, you're not getting any benefit from vector<vector<>> to make up for the downsides. (Except I guess convenience in your I/O routine which IIRC didn't check that all rows are the same length, or find out the total dimensions before starting I/O.) Use anything else, like perhaps one flat std::vector<int> which still can grow easily during I/O without caring about rows vs. cols unless you want to do 2D indexing. Commented Nov 27, 2023 at 18:50

## 2 Answers

Adding to pacmaninbw's answer:

# Naming things

matrix2 is not a matrix, it's a vector. I would just call it a vector to avoid any confusion. Of course, since you used using namespace std, you can't name a variable vector anymore. Either don't use using namespace std, or rename matrix2 to vec, or instead of naming variables after its type, name it after what they represents. For example, matrix1 could be renamed to coefficients, and matrix2 to constant_terms (as that is how they are named in the system of linear equations).

Another thing you can do is create new types or type aliases that represent matrices and vectors. For example:

using Matrix = std::vector<std::vector<int>>;
using Vector = std::vector<int>;


That way you can write:

void solveAndPrint(const Matrix& system, const Vector& constant_terms) {
…
}


# Group related data

The coefficients and the constant terms together form the system of linear equations. You could group them into a struct:

struct System {
Matrix coefficients;
Vector constant_terms;
};


That way you only have to pass one variable to solveAndPrint():

void solveAndPrint(const System& system) {
…
}


Alternatively, consider that in mathematics, a system of linear equations is often written as an augmented matrix. You could do the same in your code.

# Printing

Instead of having to call a function to print a matrix, wouldn't it be nice if you could just write something like this?

std::cout << matrix1;


This is actually easy to do. You just need to create an overload for the <<-operator. You can reuse most of your printMatrix1() for this:

std::ostream& operator<<(std::ostream& out, const Matrix& matrix) {
out << "Matrix: (" << matrix.size() << " x " << matrix[0].size() << ")\n";
for (auto& row: matrix) {
for (auto& value: row) {
out << value << ' ';
}
out << '\n';
}
return out << '\n';
}


Note that this function is no longer hardcoded to print to std::cout, so now you can also print to std::cerr, to a std::stringstream or to a file. Of course, you can also do the same for printMatrix2().

# Return values from functions

pacmaninbw already mentioned splitting up solveAndPrint(). Make a solve() function that returns the result. Of course, you can only return one thing from a function, but as shown above, you can always group things into a struct. If you reuse System, then you can write:

System solve(const System& input) {
System result = input;

for (auto& row: result.coefficients) {
…
}
…

return result;
}


Also consider creating a function that reads the system of linear equations from the command line or standard input. Then your main() can be simplified to:

int main(int argc, char* argv[]) {
System input = readInput(…);
System result = solve(input);
std::cout << result.coefficients;
std::cout << result.constant_terms;
}


# Getting rid of the magic constants

You are correct in wanting to avoid the magic constants my_magic1 and my_magic2. What if the input contained those numbers? The correct way to solve this is to create a new type that can store either a number or the state "unknown". If it's unknown, you also want to track if the sign was negative or positive during the algorithm. I see two ways to solve this:

1. Use std::variant together with an enum class that indicates positive or negative unknown:
enum class Unknown {
POSITIVE,
NEGATIVE,
};

using Value = std::variant<int, Unknown>;
using Matrix = std::vector<std::vector<Value>>;

Variants are a bit cumbersome to work with in C++ though. The alternative is:
2. Use a struct that combines both a number and a flag to indicate whether this is an unknown or not. If it is an unknown, you could keep the original value, or store 1 or -1 in the number to indicate positive/negative:
struct Value {
int value;
bool unknown;
};

using Matrix = std::vector<std::vector<Value>>;

This way the code in solve() could look like:
for (int i = 0; i < result.size(); i++) {
for (int j = 0; j < result[i].size(); j++) {
if (result[i][j].value != 0) {
if (!result2[j].unknown) {
result[i][j].value *= result2[j].value;
} else {
// keep original value, but mark it as unknown
result[i][j].unknown = true;
}
}
}
}
…


# Other minor issues

• Use const references where appropriate.

# Don't Mix C Language I/O with C++ I/O

The code is using both printf() statements combined with std::cout, this is generally a bad idea. C++ does allow for formatted output.

# Prefer "\n" Over std::endl

For performance reasons, it is better to use std::cout << "\n"; to output new lines rather than std::cout << std::endl;. std::endl flushes the output in addition to adding a newline character. This means that std::endl includes a system call every time you user it. This is especially true when you are outputting information in a loop.

# Avoid using namespace std;

If you are coding professionally you probably should get out of the habit of using the using namespace std; statement. The code will more clearly define where cout and other identifiers are coming from (std::cin, std::cout). As you start using namespaces in your code it is better to identify where each function comes from because there may be function name collisions from different namespaces. The identifier cout you may override within your own classes, and you may override the operator << in your own classes as well. This stack overflow question discusses this in more detail.

# Use Iterators

The code has many loops of the form

    for (int index = 0; index < matrix1.size(); index++)
{
Do something
}


The code would be more performant and easier to modify if it was using iterators rather than indexes.

void printMatrix1(vector<vector<int>>& matrix1)
{
printf("Matrix1: (%d x %d)\n", matrix1.size(), matrix1[0].size());
for (auto submatrix:matrix1)
{
for (auto element:submatrix)
{
cout << element << " ";
}
cout << "\n";
}
cout << "\n";
}


Just FYI, std::vector.size() returns std::size_t which is an unsigned integer rather than an integer.

# Complexity

The function main() and the function solveAndPrint() are too complex (do too much). As programs grow in size the use of main() should be limited to calling functions that parse the command line, calling functions that set up for processing, calling functions that execute the desired function of the program, and calling functions to clean up after the main portion of the program. In the function solveAndPrint() the while loop after the calls to printMatrix1() and printMatrix2() should be a separate function.

There is also a programming principle called the Single Responsibility Principle that applies here. The Single Responsibility Principle states:

that every module, class, or function should have responsibility over a single part of the functionality provided by the software, and that responsibility should be entirely encapsulated by that module, class or function.

• Those for loops should take vectors by reference (for (auto &submatrix:matrix1) or even better for (const auto &submatrix:matrix1)), otherwise you'll be copying every inner vector Commented Nov 26, 2023 at 22:55