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 矩阵乘法实现

矩阵乘法相对复杂,我们实现三种形式:

  1. 矩阵与矩阵相乘
  2. 矩阵与标量相乘
  3. 标量与矩阵相乘
// 矩阵乘法
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. 性能优化与最佳实践

在实现高性能矩阵运算时,有几个关键点需要考虑:

  1. 循环顺序优化 :按行优先顺序访问内存通常更快
  2. SIMD指令 :使用现代CPU的向量指令
  3. 表达式模板 :延迟计算以避免临时对象
  4. 并行计算 :利用多线程加速

这里展示一个简单的循环展开优化示例:

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;
}
Logo

Agent 垂直技术社区,欢迎活跃、内容共建。

更多推荐