Matrix class in C++ with addition, multiplication, concatenation and reduced echelon form operations.

matrix.cpp

#include <string>
#include <iostream>
#include <iomanip>
#include <sstream>
#include <vector>
#include <cstdint>
#include <assert.h>

template <typename T>
T typed_abs(T value) {
	if (value < 0) {
		return value * (-1);
	}
	return value;
}

template <typename T>
class Matrix {
	public:
		Matrix(uint32_t m, uint32_t n) : rows(m), columns(n) {
			matrix.resize(m, std::vector<T>(n, 0));
		}

		Matrix(std::initializer_list<std::vector<T>> m) : matrix(m) {
			rows = matrix.size();
			if (rows > 0) {
				columns = matrix[0].size();
				for (const auto& line : matrix) {
					assert(line.size() == columns);
				}
			}
		}

		class MatrixRow {
			public:
				MatrixRow(std::vector<T>& row) : row(row) {}
				T& operator[](uint32_t j) {
					return row.at(j);
				}

				void swap(MatrixRow&& other) {
					assert(row.size() == other.row.size());

					auto it = row.begin();
					auto it2 = other.row.begin();
					while (it != row.end()) {
						T tmp = *it;
						*it = *it2;
						*it2 = tmp;
						it++;
						it2++;
					}
				}

				void divideRow(T n) {
					auto it = row.begin();
					while (it != row.end()) {
						*it = (*it / n);
						it++;
					}
				}

				void addMultipleRow(T multiple, MatrixRow&& other) {
					assert(row.size() == other.row.size());

					auto it = row.begin();
					auto it2 = other.row.begin();
					while (it != row.end()) {
						*it = *it + (*it2 * multiple);
						it++;
						it2++;
					}
				}

			private:
				std::vector<T>& row;
		};

		MatrixRow operator[](uint32_t i) {
			return MatrixRow(matrix.at(i));
		}

		Matrix<T> operator+(const Matrix<T>& other) {
			assert(columns == other.columns && rows == other.rows);

			Matrix<T> ret(other);
			int i = 0;
			for (const auto& line : matrix) {
				int j = 0;
				for (const auto& value : line) {
					ret[i][j] += value;
					j++;
				}
				i++;
			}

			return ret;
		}

		Matrix<T> operator*(Matrix<T>& other) {
			assert(columns == other.rows);

			Matrix<T> ret(rows, other.columns);
			for (uint32_t i = 0; i < rows; i++) {
				for (uint32_t j = 0; j < other.columns; j++) {
					int sum = 0;
					for (uint32_t k = 0; k < columns; k++) {
						T a = matrix[i][k];
						T b = other[k][j];
						sum += a * b;
					}
					ret[i][j] = sum;
				}
			}
			return ret;
		}

		Matrix<T> reducedEchelon() {
			Matrix<T> ret(*this);
			uint32_t lead = 0;
			for (uint32_t r = 0; r < rows; r++) {
				if (columns <= lead) {
					break;
				}
				uint32_t i = r;
				while (ret[i][lead] == 0) {
					if (++i >= rows) {
						i = r;
						if (++lead >= columns) {
							return ret;
						}
					}
				}
				ret[i].swap(ret[r]);

				T tmp = ret[r][lead];
				
				for (uint32_t k = 0; k < columns; k++) {
					ret[r][k] /= tmp;
				}

				for (uint32_t i  = 0; i < rows; i++) {
					if (i != r) {
						ret[i].addMultipleRow(-ret[i][lead], ret[r]);
					}
				}
			}
			return ret;
		}

		Matrix<T> concat(const Matrix<T>& other) {
			assert(rows == other.rows);

			Matrix<T> ret(rows, columns + other.columns);
			for (uint32_t i = 0; i < rows; i++) {
				for (uint32_t j = 0; j < columns; j++) {
					ret[i][j] = matrix[i][j];
				}
			}
			for (uint32_t i = 0; i < rows; i++) {
				for (uint32_t j = 0; j < other.columns; j++) {
					ret[i][columns + j] = other.matrix[i][j];
				}
			}
			return ret;
		}

		void print(bool fraction = false) {
			auto doubleToFraction = [](double x, double error = 0.000001) {
				std::stringstream ss;

				int64_t n = x;
				x -= n;
				if (x < error) {
					ss << n;
					return ss.str();
				} else if (1 - error < x) {
					ss << n + 1;
					return ss.str();
				}

				int64_t lower_n = 0, lower_d = 1, upper_n = 1, upper_d = 1, middle_n, middle_d;

				while (1) {
					middle_n = lower_n + upper_n;
					middle_d = lower_d + upper_d;

					if (middle_d * (x + error) < middle_n) {
						upper_n = middle_n;
						upper_d = middle_d;
					} else if (middle_n < (x - error) * middle_d) {
						lower_n = middle_n;
						lower_d = middle_d;
					} else {
						ss << n * middle_d + middle_n << '/' << middle_d;
						return ss.str();
					}
				}
			};

			uint64_t max = 3;
			std::stringstream ss;
			for (const auto& line : matrix) {
				for (const auto& value : line) {
					ss.str(std::string());
					if (fraction) {
							ss << doubleToFraction(value);
					} else {
						ss << value;
					}
					ss.seekg(0, std::ios::end);
					max = std::max<uint64_t>(max, ss.tellg());
				}
			}

			for (const auto& line : matrix) {
				std::cout << std::left << std::setfill(' ');
				for (const auto& value : line) {
					std::cout << std::setw(max + 2);
					if (fraction) {
						std::cout << doubleToFraction(value);
					} else {
						std::cout << value;
					}
					std::cout << ' ';
				}
				std::cout << std::endl;
			}
		}

		uint32_t getRows() {
			return rows;
		}

		uint32_t getColumns() {
			return columns;
		}

	private:
		std::vector<std::vector<T>> matrix;
		uint32_t rows;
		uint32_t columns;
};

int main() {
	Matrix<double> A {
		{1, 1, 1},
		{0, 1, 1},
		{1, 0, 1}
	};
	
	Matrix<double> B {
		{1, 1, 0},
		{1, 0, 1},
		{0, 1, 1}
	};

	Matrix<double> v1 {
		{1},
		{0},
		{1}
	};

	Matrix<double> v2_A {
		{1},
		{1},
		{1}
	};

	Matrix<double> v3_B {
		{1},
		{0},
		{0}
	};

	// v2 na base canonica
	auto v2_C = A * v2_A;
	// v3 na base canonica
	auto v3_C = B * v3_B;

	auto s_B = B.concat(A).concat(v1).concat(v2_C).reducedEchelon();
	auto s_A = A.concat(B).concat(v1).concat(v3_C).reducedEchelon();

	std::cout << "a, d, e)" << std::endl;
	s_B.print(true);
	std::cout << std::endl;
	std::cout << "b, c, g)" << std::endl;
	s_A.print(true);
	std::cout << std::endl;

	std::cout << "f)" << std::endl;
	v2_C.print(true);
	std::cout << std::endl;
	std::cout << "h)" << std::endl;
	v3_C.print(true);
	std::cout << std::endl;

	return 0;
}

You have crucial error in operator* (line 106) there should be
T sum = 0;
otherwise it leads to rounding and ultimately wrong solutions (error increases dramatically)