🌟 WARNING: If you want to run the test, add a header (Third party dependence: https://github.com/catchorg/Catch2) in the same directory;

Matrix.hpp

#ifndef MATRIX_H
#define MATRIX_H

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <fstream>
#include <initializer_list>
#include <iostream>
#include <random>
#include <sstream>
#include <type_traits>
#include <utility>
#include <vector>

class Matrix {
 public:
  using num = double;
  static Matrix zeros(size_t rows, size_t cols);
  static Matrix eye(size_t size);
  static Matrix randu(size_t rows, size_t cols);
  static Matrix randn(size_t rows, size_t cols);

  Matrix();
  Matrix(const Matrix& other);
  Matrix& operator=(Matrix other);
  Matrix(size_t rows, size_t cols, num init_value);
  Matrix(std::initializer_list<std::initializer_list<num>> mat);
  ~Matrix();

  Matrix operator+(num k) const;
  Matrix operator-(num k) const;
  Matrix operator*(num k) const;
  Matrix operator+(const Matrix& other) const;
  Matrix operator-(const Matrix& other) const;
  Matrix operator%(const Matrix& other) const;
  Matrix operator/(const Matrix& other) const;

  Matrix operator*(const Matrix& other) const;
  Matrix operator^(size_t k) const;

  bool operator==(const Matrix& other) const;
  bool operator!=(const Matrix& other) const;

  bool write(std::string filename) const;
  bool load(std::string filename);

  size_t get_rows() const { return rows_; }
  size_t get_cols() const { return cols_; }

 private:
  size_t rows_, cols_;
  num** mat_;

  void init_(size_t rows, size_t cols, num default_value);
  void despose_();
  void debug_print_() const;

  std::string serialize_() const;
  void deserialize_(std::ifstream& in);
  static Matrix identity_(size_t n);
  Matrix get_colvec_(size_t col) const;
  void simplify_row_elems_(size_t row, size_t pivot_col);
  void gauss_jordan_elimination_();
};

Matrix AugmentMatrix(const Matrix& a, const Matrix& b);
Matrix Inverse(const Matrix& a);
Matrix GaussianElimination(const Matrix& a, const Matrix& b);

#endif

void Matrix::init_(size_t rows, size_t cols, num default_value = num()) {
  mat_ = new num*[rows];
  for (int i = 0; i < rows; ++i) mat_[i] = new num[cols];

  rows_ = rows;
  cols_ = cols;
  for (int i = 0; i < rows_; i++) {
    for (int j = 0; j < cols_; j++) {
      mat_[i][j] = default_value;
    }
  }
}

void Matrix::despose_() {
  for (int i = 0; i < rows_; i++) delete[] mat_[i];
  delete[] mat_;
}

void Matrix::debug_print_() const {
  for (int i = 0; i < rows_; i++) {
    for (int j = 0; j < cols_; j++) {
      std::cout << mat_[i][j] << " ";
    }
    std::cout << std::endl;
  }
}

Matrix::Matrix() { init_(1, 1); }

Matrix::Matrix(const Matrix& other) {
  cols_ = other.cols_;
  rows_ = other.rows_;

  mat_ = new num*[rows_];
  for (int i = 0; i < rows_; ++i) mat_[i] = new num[cols_];

  for (int i = 0; i < rows_; ++i) {
    for (int j = 0; j < cols_; ++j) {
      mat_[i][j] = other.mat_[i][j];
    }
  }
}

Matrix& Matrix::operator=(Matrix other /* const Matrix& other*/) {
  cols_ = other.cols_;
  rows_ = other.rows_;
  std::swap(mat_, other.mat_);
  return *this;
}

Matrix::Matrix(size_t rows, size_t cols, num init_value = num()) {
  init_(rows, cols, init_value);
}

Matrix::Matrix(std::initializer_list<std::initializer_list<num>> mat) {
  if (mat.begin() == mat.end()) {
    init_(1, 1);
  } else {
    init_(mat.size(), mat.begin()->size());
    int i = 0;
    for (auto row : mat) {
      int j = 0;
      for (auto e : row) {
        mat_[i][j] = e;
        j++;
      }
      i++;
    }
  }
}

Matrix::~Matrix() { despose_(); }

Matrix Matrix::operator+(num k) const {
  Matrix m = *this;
  for (size_t i = 0; i < rows_; i++) {
    for (size_t j = 0; j < cols_; j++) {
      m.mat_[i][j] += k;
    }
  }
  return m;
}

Matrix Matrix::operator-(num k) const {
  Matrix m = *this;
  for (size_t i = 0; i < rows_; i++) {
    for (size_t j = 0; j < cols_; j++) {
      m.mat_[i][j] -= k;
    }
  }
  return m;
}

Matrix Matrix::operator*(num k) const {
  Matrix m = *this;
  for (size_t i = 0; i < rows_; i++) {
    for (size_t j = 0; j < cols_; j++) {
      m.mat_[i][j] *= k;
    }
  }
  return m;
}

Matrix Matrix::operator+(const Matrix& other) const {
  // NOTE: Why it will cause error: (SIG)abnormal error?
  // assert(*this == other);
  Matrix m = other;
  for (size_t i = 0; i < rows_; i++) {
    for (size_t j = 0; j < cols_; j++) {
      m.mat_[i][j] += mat_[i][j];
    }
  }
  return m;
}

Matrix Matrix::operator-(const Matrix& other) const {
  Matrix m = *this;
  for (size_t i = 0; i < rows_; i++) {
    for (size_t j = 0; j < cols_; j++) {
      m.mat_[i][j] -= other.mat_[i][j];
    }
  }
  return m;
}

Matrix Matrix::operator%(const Matrix& other) const {
  Matrix m = *this;
  for (size_t i = 0; i < rows_; i++) {
    for (size_t j = 0; j < cols_; j++) {
      m.mat_[i][j] *= other.mat_[i][j];
    }
  }
  return m;
}

Matrix Matrix::operator/(const Matrix& other) const {
  Matrix m = *this;
  for (size_t i = 0; i < rows_; i++) {
    for (size_t j = 0; j < cols_; j++) {
      m.mat_[i][j] /= other.mat_[i][j];
    }
  }
  return m;
}

Matrix Matrix::get_colvec_(size_t col) const {
  Matrix m(rows_, 1);
  for (size_t i = 0; i < rows_; i++) {
    m.mat_[i][0] = mat_[i][col];
  }
  return m;
}

Matrix Matrix::operator*(const Matrix& other) const {
  Matrix mat(rows_, other.cols_);
  int cnt = cols_;
  double tmp = 0;
  for (size_t i = 0; i < rows_; i++) {
    for (size_t j = 0; j < other.cols_; j++) {
      for (size_t n = 0; n < cnt; n++) {
        tmp += mat_[i][n] * other.mat_[n][j];
      }
      mat.mat_[i][j] = tmp;
      tmp = 0;
    }
  }
  return mat;
}

Matrix Matrix::operator^(size_t k) const {
  auto res = Matrix::identity_(rows_);

  for (Matrix t = *this; k;) {
    if (k & 1)  // k is odd
      res = res * t;
    t = t * t;
    k = k >> 1;  // k / 2
  }
  return res;
}

template <typename T>
static bool equal(T lhs, T rhs) {
  if (std::is_same<T, double>()) return std::fabs(lhs - rhs) <= 1e-8;
  return lhs == rhs;
}

bool Matrix::operator==(const Matrix& other) const {
  if (rows_ != other.rows_ || cols_ != other.cols_) return false;
  for (size_t i = 0; i < rows_; i++) {
    for (size_t j = 0; j < cols_; j++) {
      if (!equal<num>(mat_[i][j], other.mat_[i][j])) return false;
    }
  }
  return true;
}

bool Matrix::operator!=(const Matrix& other) const { return !(*this == other); }

Matrix Matrix::zeros(size_t rows, size_t cols) { return Matrix(rows, cols, 0); }

Matrix Matrix::eye(size_t size) {
  int n = std::sqrt(size);
  auto other = Matrix(n, n);
  for (int i = 0; i < n; i++) {
    other.mat_[i][i] = 1;
  }
  return other;
}

Matrix Matrix::randu(size_t rows, size_t cols) {
  std::random_device rd;
  std::mt19937 gen(rd());
  std::uniform_real_distribution<> dis(1.0, 2.0);

  auto other = Matrix(rows, cols);
  for (int i = 0; i < rows; i++) {
    for (int j = 0; j < cols; j++) {
      other.mat_[i][j] = dis(gen);
    }
  }
  return other;
}

Matrix Matrix::randn(size_t rows, size_t cols) {
  std::random_device rd;
  std::mt19937 gen(rd());
  std::normal_distribution<> dis{5, 2};

  auto other = Matrix(rows, cols);
  for (int i = 0; i < rows; i++) {
    for (int j = 0; j < cols; j++) {
      other.mat_[i][j] = dis(gen);
    }
  }
  return other;
}

Matrix Matrix::identity_(size_t n) {
  Matrix i(n, n, 0);
  for (size_t k = 0; k < n; k++) {
    i.mat_[k][k] = 1;
  }

  return i;
}

std::string Matrix::serialize_() const {
  std::string result;
  for (size_t i = 0; i < rows_; i++) {
    for (size_t j = 0; j < cols_; j++) {
      result += std::to_string(mat_[i][j]);
      result += ' ';
    }
    result += '\n';
  }
  return result;
}

bool Matrix::write(std::string filename) const {
  std::ofstream ofstrm(filename, std::ios::out);
  ofstrm << serialize_();
  return true;
}

void Matrix::deserialize_(std::ifstream& in) {
  std::vector<std::vector<num>> m;
  int i = -1;
  for (std::string row; std::getline(in, row);) {
    m.push_back({});
    i++;
    std::stringstream ss(row);
    for (std::string s; ss >> s;) {
      m[i].push_back(std::stod(s));
    }
  }

  despose_();
  init_(m.size(), m[0].size());
  for (size_t i = 0; i < rows_; i++) {
    for (size_t j = 0; j < cols_; j++) {
      mat_[i][j] = m[i][j];
    }
  }
}

bool Matrix::load(std::string filename) {
  std::ifstream ifstrm(filename, std::ios::in);
  deserialize_(ifstrm);
  return true;
}

Matrix AugmentMatrix(const Matrix& a, const Matrix& b) {
  auto rows = a.get_rows();
  auto cols = a.get_cols() + b.get_cols();
  Matrix m(rows, cols);
  for (size_t i = 0; i < rows; i++) {
    for (size_t j = 0; j < a.get_cols(); j++) {
      m.mat_[i][j] = a.mat_[i][j];
    }
  }

  for (size_t i = 0; i < rows; i++) {
    for (size_t j = 0; j < b.get_cols(); j++) {
      m.mat_[i][j + a.get_cols()] = b.mat_[i][j];
    }
  }
  return m;
}

void Matrix::simplify_row_elems_(size_t row, size_t pivot_col) {
  auto factor = mat_[row][pivot_col];
  for (size_t i = pivot_col; i < cols_; i++) {
    mat_[row][i] /= factor;
  }
}

static size_t find_pivot(Matrix::num* row_elems, size_t rows) {
  for (size_t i = 0; i < rows; i++) {
    if (!equal<Matrix::num>(row_elems[i], 0)) {
      return i;
    }
  }
  return SIZE_MAX;
}

void Matrix::gauss_jordan_elimination_() {
  num max_elem = num();
  for (size_t max_i = 0, h = 0, k = 0; h < rows_ && k < cols_;) {
    for (size_t i = h; i < rows_; i++) {
      if (max_elem < std::fabs(mat_[i][k])) {
        max_elem = std::fabs(mat_[i][k]);
        max_i = i;
      }
    }
    if (equal<num>(mat_[max_i][k], 0)) {
      k++;
    } else {
      std::swap(mat_[h], mat_[max_i]);
      for (size_t i = h + 1; i < rows_; i++) {
        num factor = mat_[i][k] / mat_[h][k];
        mat_[i][k] = 0;
        for (size_t j = k + 1; j < cols_; j++) {
          mat_[i][j] = mat_[i][j] - mat_[h][j] * factor;
        }
      }
      h++;
      k++;
    }
  }

  std::vector<std::pair<size_t, size_t>> order;

  for (size_t pivot, i = 0; i < rows_; i++) {
    pivot = find_pivot(mat_[i], rows_);
    if (pivot != SIZE_MAX) {
      order.push_back({i, pivot});
      simplify_row_elems_(i, pivot);
    }
  }

  for (size_t i = 0; i < order.size() - 1; i++) {
    for (size_t j = 0; j < order.size() - i - 1; j++) {
      if (order[j].second > order[j + 1].second) {
        std::swap(mat_[order[j].first], mat_[order[j + 1].first]);
        std::swap(order[j].second, order[j + 1].second);
      }
    }
  }

  for (size_t pivot = SIZE_MAX, i = rows_ - 1; i >= 1; i--, pivot = SIZE_MAX) {
    for (size_t j = 0; j < cols_; j++) {
      if (equal<num>(mat_[i][j], 1)) {
        pivot = j;
        break;
      };
    }

    // NOTE: underflow error!
    // solution(now): use implicit conversion: size_t(`unsigned`) to
    // long/int(`signed type`)
    for (long j = i - 1; j >= 0; j--) {
      if (pivot != SIZE_MAX && !equal<num>(mat_[j][pivot], 0)) {
        auto factor = mat_[j][pivot];
        for (size_t k = pivot; k < cols_; k++) {
          mat_[j][k] -= mat_[i][k] * factor;
        }
      }
    }
  }
}

Matrix Inverse(const Matrix& a) {
  Matrix aI = AugmentMatrix(a, Matrix::identity_(a.rows_));

  aI.gauss_jordan_elimination_();
  Matrix m(a.rows_, a.cols_);

  for (size_t i = 0; i < aI.rows_; i++) {
    for (size_t j = 0, offset = a.cols_; j + offset < aI.cols_; j++) {
      m.mat_[i][j] = aI.mat_[i][j + offset];
    }
  }
  return m;
}

Matrix GaussianElimination(const Matrix& a, const Matrix& b) {
  auto ab = AugmentMatrix(a, b);
  ab.gauss_jordan_elimination_();

  return ab.get_colvec_(ab.cols_ - 1);
}

test.cpp

#define CATCH_CONFIG_MAIN  // This tells Catch to provide a main() - only do
                           // this in one cpp file
#include "catch.hpp"

// Test private members:
// - Use compiler `define(D)` option:
//   - like: `$ g++ -Dprivate=public -std=c++14 test.cc`
// - Or just write `#define private public` in your test file.
#define private public  // WARNING: Don't move me
#include "Matrix.hpp"

TEST_CASE("2) 矩阵类的构造函数", "[Matrix]") {
  SECTION("a) 默认构造函数，生成一个 1,1 列的矩阵") {
    Matrix mat;
    REQUIRE(mat.get_cols() == 1);
    REQUIRE(mat.get_rows() == 1);
    REQUIRE(mat.mat_[0][0] == 0.0);
  }

  SECTION("b) 生成行列分别为 rows, cols 的矩阵") {
    Matrix mat(5, 5);
    REQUIRE(mat.get_cols() == 5);
    REQUIRE(mat.get_rows() == 5);
    REQUIRE(mat.mat_[0][0] == 0.0);
  }

  SECTION("c) 生成行列分别为 rows, cols,初始值为 init_value 的矩阵") {
    const double INIT_VALUE = 6.66;
    Matrix mat(5, 5, INIT_VALUE);
    REQUIRE(mat.get_cols() == 5);
    REQUIRE(mat.get_rows() == 5);
    REQUIRE(mat.mat_[0][0] == INIT_VALUE);
  }

  SECTION("d) 利用初始化列表生成矩阵") {
    SECTION("空列表") {
      Matrix mat = {};

      REQUIRE(mat.get_cols() == 1);
      REQUIRE(mat.get_rows() == 1);
      REQUIRE(mat.mat_[0][0] == 0.0);
    }

    SECTION("{{1,2,3,4}, {5,6,7,8}}") {
      Matrix mat = {
          {1, 2, 3, 4},
          {5, 6, 7, 8},
      };
      REQUIRE(mat.get_cols() == 4);
      REQUIRE(mat.get_rows() == 2);
      REQUIRE(mat.mat_[0][0] == 1);
      REQUIRE(mat.mat_[1][0] == 5);
    }
  }
}

TEST_CASE("4) 拷贝构造函数与赋值运算符函数(实现深拷贝)", "[Matrix]") {
  const size_t R = 5, C = 5;
  Matrix m(R, C);

  SECTION("拷贝构造函数") {
    Matrix cc(m);
    REQUIRE(cc.get_cols() == R);
    REQUIRE(cc.get_rows() == C);
    REQUIRE(cc.mat_ != m.mat_);
  }

  SECTION("赋值运算符函数") {
    Matrix as = m;
    REQUIRE(as.get_cols() == R);
    REQUIRE(as.get_rows() == C);
    REQUIRE(as.mat_ != m.mat_);
  }
}

TEST_CASE("5) 实现基本运算符重载函数", "[Matrix]") {
  SECTION("A + k 矩阵对应元素加上标量") {
    Matrix a = {
        {1, 1, 1},
        {1, 1, 1},
    };
    Matrix::num k = 1;
    Matrix actual = a + k;
    Matrix expected = {
        {2, 2, 2},
        {2, 2, 2},
    };
    REQUIRE(actual == expected);
  }

  SECTION("A - k 矩阵对应元素减去标量") {
    Matrix a = {
        {1, 1, 1},
        {1, 1, 1},
    };
    Matrix::num k = 1;
    Matrix actual = a - k;
    Matrix expected = {
        {0, 0, 0},
        {0, 0, 0},
    };
    REQUIRE(actual == expected);
  }

  SECTION("A * k 矩阵对应元素乘以标量") {
    Matrix a = {
        {1, 1, 1},
        {1, 1, 1},
    };
    Matrix::num k = 2;
    Matrix actual = a * k;
    Matrix expected = {
        {2, 2, 2},
        {2, 2, 2},
    };
    REQUIRE(actual == expected);
  }

  SECTION("A + B 矩阵对应元素加法") {
    Matrix a = {
        {1, 1, 1},
        {1, 1, 1},
    };
    Matrix b = {
        {1, 1, 1},
        {0, 0, 0},
    };
    Matrix actual = a + b;
    Matrix expected = {
        {2, 2, 2},
        {1, 1, 1},
    };
    REQUIRE(actual == expected);
  }

  SECTION("A - B 矩阵对应元素减法") {
    Matrix a = {
        {1, 1, 1},
        {1, 1, 1},
    };
    Matrix b = {
        {1, 1, 1},
        {0, 0, 0},
    };
    Matrix actual = a - b;
    Matrix expected = {
        {0, 0, 0},
        {1, 1, 1},
    };
    REQUIRE(actual == expected);
  }

  SECTION("A % B 矩阵对应元素乘法") {
    Matrix a = {
        {1, 1, 1},
        {1, 1, 1},
    };
    Matrix b = {
        {2, 2, 2},
        {0, 0, 0},
    };
    Matrix actual = a % b;
    Matrix expected = {
        {2, 2, 2},
        {0, 0, 0},
    };
    REQUIRE(actual == expected);
  }

  SECTION("A / B 矩阵对应元素除法") {
    Matrix a = {
        {4, 4, 4},
        {2, 2, 2},
    };
    Matrix b = {
        {2, 2, 2},
        {2, 2, 2},
    };
    Matrix actual = a / b;
    Matrix expected = {
        {2, 2, 2},
        {1, 1, 1},
    };
    REQUIRE(actual == expected);
  }

  SECTION("A * B 矩阵乘法") {
    Matrix a = {
        {1, 2, 3},
        {4, 5, 6},
    };

    Matrix b = {
        {100.000, 10.000},
        {10.000, 100.000},
        {1.000, 1000.000},
    };
    Matrix actual = a * b;
    Matrix expected = {
        {123, 3210},
        {456, 6540},
    };
    REQUIRE(actual == expected);
  }

  SECTION("A ^ k 矩阵k次幂") {
    Matrix a = {
        {1, 2, 3},
        {4, 5, 6},
        {7, 8, 9},
    };

    Matrix actual = a ^ 3;
    Matrix expected = {
        {468, 576, 684},
        {1062, 1305, 1548},
        {1656, 2034, 2412},
    };

    REQUIRE(actual == expected);
  }
}

TEST_CASE("6) 静态函数生成基本矩阵", "[Matrix]") {
  SECTION("a) 零矩阵 Matrix::zeros(rows, cols);") {
    auto zerom = Matrix::zeros(3, 3);
    for (int i = 0; i < 3; i++) {
      for (int j = 0; j < 3; j++) {
        REQUIRE(zerom.mat_[i][j] == 0);
      }
    }
  }

  SECTION("b) 单位矩阵 Matrix::eye(size);") {
    auto eyem = Matrix::eye(9);
    for (int i = 0; i < 3; i++) {
      REQUIRE(eyem.mat_[i][i] == 1);
    }
  }

  SECTION("c) 均匀分布随机矩阵 1 Matrix::randu(rows, cols);") {
    auto randum = Matrix::randu(3, 3);
    // randum.debug_print_();
    WARN("我暂时不知道如何测试这个 QwQ");
  }

  SECTION("d) 高斯分布随机矩阵 2 Matrix::randn(rows, cols);") {
    auto randum = Matrix::randn(3, 3);
    // randum.debug_print_();
    WARN("我暂时不知道如何测试这个 QwQ");
  }
}

TEST_CASE("7) 应用于方程组的矩阵函数", "[Matrix]") {
  SECTION("a) 增广矩阵") {
    Matrix a = {
        {1, 2, 3},
        {4, 5, 6},
    };
    Matrix b = {
        {4, 5},
        {7, 8},
    };

    Matrix expected = {
        {1, 2, 3, 4, 5},
        {4, 5, 6, 7, 8},
    };

    auto actual = AugmentMatrix(a, b);
    REQUIRE(actual == expected);
  }

  SECTION("b) 逆矩阵") {
    Matrix a = {
        {1, 2},
        {3, 4},
    };

    auto actual = Inverse(a);
    Matrix expected = {
        {-2, 1},
        {1.5, -0.5},
    };
    REQUIRE(actual == expected);
  }

  SECTION("c) 实现高斯消元法求解线性方程组") {
    Matrix a = {
        {2, 1, -1},
        {-3, -1, 2},
        {-2, 1, 2},
    };
    Matrix b = {
        {8},
        {-11},
        {-3},
    };

    auto actual = GaussianElimination(a, b);
    Matrix expected = {
        {2},
        {3},
        {-1},
    };

    REQUIRE(actual == expected);
  }
}

TEST_CASE("9) 实现成员函数实现将矩阵读写函数", "[Matrix]") {
  Matrix actual = {
      {1, 2, 3},
      {4, 5, 6},
      {7, 8, 9},
  };
  actual.write("tmp.txt");
  Matrix expected;
  expected.load("tmp.txt");
  REQUIRE(actual == expected);
}

TEST_CASE("额外的", "[Matrix]") {
  SECTION("操纵列向量") {
    Matrix m = {
        {1, 0, 0},
        {2, 0, 0},
        {3, 0, 0},
    };

    auto actual = m.get_colvec_(0);
    Matrix expected = {
        {1},
        {2},
        {3},
    };
    REQUIRE(actual == expected);
  }

  SECTION("单位矩阵 I") {
    auto actual = Matrix::identity_(3);
    Matrix expected = {
        {1, 0, 0},
        {0, 1, 0},
        {0, 0, 1},
    };
    REQUIRE(actual == expected);
  }

  SECTION("矩阵判等 (==,!=)") {
    Matrix actual = {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}};
    Matrix rows_not_match = {{1, 2, 3}, {4, 5, 6}};
    Matrix elem_not_match = {{1, 2, 3}, {41, 5, 6}, {7, 8, 9}};
    Matrix expected = {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}};
    REQUIRE(actual == expected);
    REQUIRE(actual != rows_not_match);
    REQUIRE(actual != elem_not_match);
  }

  SECTION("化简指定行元素") {
    Matrix m = {{2, 4, 6}};
    SECTION("指定列") {
      Matrix actual(m);
      actual.simplify_row_elems_(0, 0);
      Matrix expected = {{1, 2, 3}};

      REQUIRE(actual == expected);
    }
  }

  SECTION("化最简阶梯型 (高斯－若尔当消元法) ") {
    Matrix actual = {
        {2, 1, -1, 8},
        {-3, -1, 2, -11},
        {-2, 1, 2, -3},
    };

    actual.gauss_jordan_elimination_();
    Matrix expected = {
        {1, 0, 0, 2},
        {0, 1, 0, 3},
        {0, 0, 1, -1},
    };

    REQUIRE(actual == expected);
  }
}