C++17 Matrix 类模板进阶:5个关键成员函数与运算符重载实战
·
C++17 Matrix 类模板进阶:5个关键成员函数与运算符重载实战
在科学计算和工程应用中,矩阵运算是最基础也是最重要的数学工具之一。本文将带领你从零开始构建一个功能完备的矩阵类模板,不仅支持基本的构造和拷贝控制,还能实现算术运算、流输出等高级功能。这个实现将充分利用C++17的特性,为中级C++开发者提供一个可复用的高质量数学库。
1. 矩阵类模板基础架构
我们先从最基础的矩阵存储开始。一个通用的矩阵类需要能够存储任意数据类型的元素,并管理内存的生命周期。以下是我们的基础实现:
template <typename T>
class Matrix {
private:
size_t rows_;
size_t cols_;
std::unique_ptr<T[]> data_; // 使用智能指针管理内存
public:
// 基础构造函数
Matrix(size_t rows, size_t cols)
: rows_(rows), cols_(cols), data_(new T[rows * cols]) {}
// 拷贝构造函数
Matrix(const Matrix& other)
: rows_(other.rows_), cols_(other.cols_), data_(new T[rows_ * cols_]) {
std::copy(other.data_.get(), other.data_.get() + rows_ * cols_, data_.get());
}
// 移动构造函数
Matrix(Matrix&& other) noexcept
: rows_(other.rows_), cols_(other.cols_), data_(std::move(other.data_)) {
other.rows_ = 0;
other.cols_ = 0;
}
// 析构函数 - 由unique_ptr自动处理
~Matrix() = default;
// 获取行列数
size_t rows() const { return rows_; }
size_t cols() const { return cols_; }
// 元素访问
T& operator()(size_t row, size_t col) {
return data_[row * cols_ + col];
}
const T& operator()(size_t row, size_t col) const {
return data_[row * cols_ + col];
}
};
这个基础版本已经提供了:
- 安全的动态内存管理(使用
unique_ptr) - 完整的拷贝控制(构造、拷贝、移动)
- 方便的元素访问接口
2. 核心运算符重载实现
让矩阵支持基本的算术运算是提升实用性的关键。我们将实现加法、减法和乘法运算符。
2.1 加法和减法运算符
// 矩阵加法
template <typename T>
Matrix<T> operator+(const Matrix<T>& lhs, const Matrix<T>& rhs) {
if (lhs.rows() != rhs.rows() || lhs.cols() != rhs.cols()) {
throw std::invalid_argument("Matrix dimensions mismatch for addition");
}
Matrix<T> result(lhs.rows(), lhs.cols());
for (size_t i = 0; i < lhs.rows(); ++i) {
for (size_t j = 0; j < lhs.cols(); ++j) {
result(i, j) = lhs(i, j) + rhs(i, j);
}
}
return result;
}
// 矩阵减法
template <typename T>
Matrix<T> operator-(const Matrix<T>& lhs, const Matrix<T>& rhs) {
if (lhs.rows() != rhs.rows() || lhs.cols() != rhs.cols()) {
throw std::invalid_argument("Matrix dimensions mismatch for subtraction");
}
Matrix<T> result(lhs.rows(), lhs.cols());
for (size_t i = 0; i < lhs.rows(); ++i) {
for (size_t j = 0; j < lhs.cols(); ++j) {
result(i, j) = lhs(i, j) - rhs(i, j);
}
}
return result;
}
2.2 矩阵乘法实现
矩阵乘法相对复杂,我们实现三种形式:
- 矩阵与矩阵相乘
- 矩阵与标量相乘
- 标量与矩阵相乘
// 矩阵乘法
template <typename T>
Matrix<T> operator*(const Matrix<T>& lhs, const Matrix<T>& rhs) {
if (lhs.cols() != rhs.rows()) {
throw std::invalid_argument("Matrix dimensions mismatch for multiplication");
}
Matrix<T> result(lhs.rows(), rhs.cols());
for (size_t i = 0; i < lhs.rows(); ++i) {
for (size_t j = 0; j < rhs.cols(); ++j) {
T sum{};
for (size_t k = 0; k < lhs.cols(); ++k) {
sum += lhs(i, k) * rhs(k, j);
}
result(i, j) = sum;
}
}
return result;
}
// 标量乘法(矩阵 * 标量)
template <typename T>
Matrix<T> operator*(const Matrix<T>& mat, T scalar) {
Matrix<T> result(mat.rows(), mat.cols());
for (size_t i = 0; i < mat.rows(); ++i) {
for (size_t j = 0; j < mat.cols(); ++j) {
result(i, j) = mat(i, j) * scalar;
}
}
return result;
}
// 标量乘法(标量 * 矩阵)
template <typename T>
Matrix<T> operator*(T scalar, const Matrix<T>& mat) {
return mat * scalar; // 重用上面的实现
}
3. 流输出与输入运算符
为了方便调试和使用,我们重载 << 和 >> 运算符:
// 输出运算符
template <typename T>
std::ostream& operator<<(std::ostream& os, const Matrix<T>& mat) {
for (size_t i = 0; i < mat.rows(); ++i) {
for (size_t j = 0; j < mat.cols(); ++j) {
os << mat(i, j);
if (j != mat.cols() - 1) os << ' ';
}
if (i != mat.rows() - 1) os << '\n';
}
return os;
}
// 输入运算符
template <typename T>
std::istream& operator>>(std::istream& is, Matrix<T>& mat) {
for (size_t i = 0; i < mat.rows(); ++i) {
for (size_t j = 0; j < mat.cols(); ++j) {
is >> mat(i, j);
}
}
return is;
}
4. 五个核心成员函数实现
除了基本运算,一个实用的矩阵类还需要提供一些核心数学运算功能。
4.1 转置函数
template <typename T>
Matrix<T> Matrix<T>::transpose() const {
Matrix<T> result(cols_, rows_);
for (size_t i = 0; i < rows_; ++i) {
for (size_t j = 0; j < cols_; ++j) {
result(j, i) = (*this)(i, j);
}
}
return result;
}
4.2 行列式计算
行列式计算使用递归的拉普拉斯展开法:
template <typename T>
T Matrix<T>::determinant() const {
if (rows_ != cols_) {
throw std::logic_error("Determinant is only defined for square matrices");
}
if (rows_ == 1) return (*this)(0, 0);
if (rows_ == 2) {
return (*this)(0, 0) * (*this)(1, 1) - (*this)(0, 1) * (*this)(1, 0);
}
T det{};
for (size_t col = 0; col < cols_; ++col) {
Matrix<T> minor(rows_ - 1, cols_ - 1);
for (size_t i = 1; i < rows_; ++i) {
for (size_t j = 0, k = 0; j < cols_; ++j) {
if (j == col) continue;
minor(i - 1, k++) = (*this)(i, j);
}
}
T sign = (col % 2 == 0) ? 1 : -1;
det += sign * (*this)(0, col) * minor.determinant();
}
return det;
}
4.3 逆矩阵计算
使用伴随矩阵法计算逆矩阵:
template <typename T>
Matrix<T> Matrix<T>::inverse() const {
T det = determinant();
if (det == T{}) {
throw std::runtime_error("Matrix is singular, cannot compute inverse");
}
Matrix<T> adj(rows_, cols_);
for (size_t i = 0; i < rows_; ++i) {
for (size_t j = 0; j < cols_; ++j) {
Matrix<T> minor(rows_ - 1, cols_ - 1);
for (size_t m = 0, p = 0; m < rows_; ++m) {
if (m == i) continue;
for (size_t n = 0, q = 0; n < cols_; ++n) {
if (n == j) continue;
minor(p, q++) = (*this)(m, n);
}
++p;
}
T sign = ((i + j) % 2 == 0) ? 1 : -1;
adj(j, i) = sign * minor.determinant(); // 注意这里是(j,i)不是(i,j)
}
}
return adj * (1 / det);
}
4.4 迹(trace)计算
template <typename T>
T Matrix<T>::trace() const {
if (rows_ != cols_) {
throw std::logic_error("Trace is only defined for square matrices");
}
T tr{};
for (size_t i = 0; i < rows_; ++i) {
tr += (*this)(i, i);
}
return tr;
}
4.5 单位矩阵生成
template <typename T>
Matrix<T> Matrix<T>::identity(size_t n) {
Matrix<T> mat(n, n);
for (size_t i = 0; i < n; ++i) {
for (size_t j = 0; j < n; ++j) {
mat(i, j) = (i == j) ? T{1} : T{};
}
}
return mat;
}
5. 模板特化与优化
对于某些特定类型,我们可以进行优化。例如,对于bool类型的矩阵,可以使用更紧凑的存储方式。
5.1 bool矩阵特化
template <>
class Matrix<bool> {
private:
size_t rows_;
size_t cols_;
std::vector<std::bitset<64>> data_; // 每64个bool用1个bitset存储
public:
Matrix(size_t rows, size_t cols)
: rows_(rows), cols_(cols), data_((rows * cols + 63) / 64) {}
class BoolProxy {
private:
std::bitset<64>& chunk;
size_t pos;
public:
BoolProxy(std::bitset<64>& c, size_t p) : chunk(c), pos(p) {}
operator bool() const { return chunk[pos]; }
BoolProxy& operator=(bool val) { chunk[pos] = val; return *this; }
};
BoolProxy operator()(size_t row, size_t col) {
size_t index = row * cols_ + col;
return BoolProxy(data_[index / 64], index % 64);
}
bool operator()(size_t row, size_t col) const {
size_t index = row * cols_ + col;
return data_[index / 64][index % 64];
}
// 其他成员函数...
};
5.2 针对小型矩阵的优化
对于小型矩阵(如3x3或4x4),可以使用栈存储而非堆存储:
template <typename T, size_t Rows = 3, size_t Cols = 3>
class SmallMatrix {
private:
std::array<T, Rows * Cols> data_;
public:
SmallMatrix() = default;
T& operator()(size_t row, size_t col) {
return data_[row * Cols + col];
}
const T& operator()(size_t row, size_t col) const {
return data_[row * Cols + col];
}
// 其他成员函数...
};
6. 实际应用示例
让我们看一个使用这个矩阵类解决实际问题的例子:解线性方程组。
template <typename T>
std::vector<T> solveLinearSystem(const Matrix<T>& A, const std::vector<T>& b) {
if (A.rows() != A.cols() || A.rows() != b.size()) {
throw std::invalid_argument("Invalid dimensions for linear system");
}
Matrix<T> invA = A.inverse();
std::vector<T> x(b.size());
for (size_t i = 0; i < A.rows(); ++i) {
x[i] = {};
for (size_t j = 0; j < A.cols(); ++j) {
x[i] += invA(i, j) * b[j];
}
}
return x;
}
这个实现虽然简单,但对于小型方程组已经足够。对于大型系统,你可能需要实现更高效的算法如LU分解。
7. 性能优化与最佳实践
在实现高性能矩阵运算时,有几个关键点需要考虑:
- 循环顺序优化 :按行优先顺序访问内存通常更快
- SIMD指令 :使用现代CPU的向量指令
- 表达式模板 :延迟计算以避免临时对象
- 并行计算 :利用多线程加速
这里展示一个简单的循环展开优化示例:
template <typename T>
Matrix<T> operator*(const Matrix<T>& lhs, const Matrix<T>& rhs) {
if (lhs.cols() != rhs.rows()) {
throw std::invalid_argument("Matrix dimensions mismatch");
}
Matrix<T> result(lhs.rows(), rhs.cols());
constexpr size_t blockSize = 4;
for (size_t i = 0; i < lhs.rows(); ++i) {
for (size_t j = 0; j < rhs.cols(); j += blockSize) {
for (size_t k = 0; k < lhs.cols(); ++k) {
for (size_t b = 0; b < blockSize && (j + b) < rhs.cols(); ++b) {
result(i, j + b) += lhs(i, k) * rhs(k, j + b);
}
}
}
}
return result;
}
更多推荐
所有评论(0)