矩阵设计：从基础到应用

# 矩阵设计：从基础到应用 ## 1. 引言 在编程中，孤立的语言特性往往既无趣又无用。本文将展示如何结合多种语言特性和编程技巧，设计并实现一个通用的 N 维矩阵。虽然完美的矩阵类可能并不存在，但我们可以实现一个简单且实用的 N 维密集矩阵。 ### 1.1 矩阵的基本使用 `Matrix<T,N>` 是一个 N 维矩阵，其中 `T` 是元素类型。以下是一些使用示例： ```cpp Matrix<double, 0> m0 {1}; // 零维矩阵：标量 Matrix<double, 1> m1 {1, 2, 3, 4}; // 一维矩阵：向量（4 个元素） Matrix<double, 2> m2 { {00, 01, 02, 03}, // 第 0 行 {10, 11, 12, 13}, // 第 1 行 {20, 21, 22, 23} // 第 2 行 }; Matrix<double, 3> m3(4, 7, 9); // 三维矩阵（4*7*9 个元素），全部初始化为 0 Matrix<complex<double>, 17> m17; // 17 维矩阵（目前无元素） ``` 元素类型可以是任何可存储的类型，但矩阵运算的数学性质与整数或浮点运算有所不同，例如矩阵乘法不满足交换律。在初始化矩阵时，我们使用 `()` 指定维度大小，使用 `{}` 指定元素值。同时，维度和元素数量必须匹配，否则会导致错误。 矩阵的一些基本属性包括： - `order()`：返回矩阵的维度数。 - `extent(n)`：返回第 `n` 维的元素数量。 - `size()`：返回矩阵的总元素数量。 我们可以通过多种下标方式访问矩阵元素，例如： ```cpp Matrix<double, 2> m { {00, 01, 02, 03}, {10, 11, 12, 13}, {20, 21, 22, 23} }; double d1 = m(1, 2); // d1 == 12 double d2 = m[1][2]; // d2 == 12 Matrix<double, 1> m1 = m[1]; // 第 1 行：{10, 11, 12, 13} double d3 = m1[2]; // d3 == 12 ``` 为了方便调试，我们可以定义一个输出函数： ```cpp template<typename M> Enable_if<Matrix_type<M>(), ostream&> operator<<(ostream& os, const M& m) { os << '{'; for (size_t i = 0; i != rows(m); ++i) { os << m[i]; if (i + 1 != rows(m)) os << ','; } return os << '}'; } ``` ### 1.2 矩阵的需求 在实现矩阵之前，我们需要考虑以下需求： - **维度灵活性**：支持 0 到多个维度，无需为每个维度编写专门代码。 - **元素类型通用性**：元素类型可以是任何可存储的类型。 - **数学运算适用性**：数学运算适用于任何可被视为数字的类型，包括矩阵。 - **下标方式**：支持 Fortran 风格（如 `m(1, 2, 3)`）和 C 风格（如 `m[7]`）的下标访问。 - **下标性能和范围检查**：下标访问应快速且可进行范围检查。 - **移动语义**：支持移动赋值和移动构造，以提高矩阵结果的传递效率，减少临时对象的开销。 - **数学运算**：支持一些基本的矩阵运算，如 `+` 和 `*=`。 - **子矩阵引用**：支持读取、写入和传递子矩阵的引用。 - **资源管理**：避免资源泄漏，提供基本的资源保证。 - **融合操作**：支持融合关键操作，如 `m*v + v2` 作为单个函数调用。 为了满足这些需求，我们将结合以下语言特性和编程技巧： - 类 - 数字和类型参数化 - 移动构造和赋值 - RAII（资源获取即初始化） - 可变参数模板 - 初始化列表 - 运算符重载 - 函数对象 - 简单的模板元编程 - 实现继承 ## 2. 矩阵模板 以下是 `Matrix` 类的声明，包含了一些重要的操作： ```cpp template<typename T, size_t N> class Matrix { public: static constexpr size_t order = N; using value_type = T; using iterator = typename std::vector<T>::iterator; using const_iterator = typename std::vector<T>::const_iterator; Matrix() = default; Matrix(Matrix&&) = default; // 移动构造 Matrix& operator=(Matrix&&) = default; Matrix(Matrix const&) = default; // 拷贝构造 Matrix& operator=(Matrix const&) = default; ~Matrix() = default; template<typename U> Matrix(const Matrix_ref<U, N>&); // 从 Matrix_ref 构造 template<typename U> Matrix& operator=(const Matrix_ref<U, N>&); // 从 Matrix_ref 赋值 template<typename... Exts> // 指定维度大小 explicit Matrix(Exts... exts); Matrix(Matrix_initializer<T, N>); // 从初始化列表构造 Matrix& operator=(Matrix_initializer<T, N>); // 从初始化列表赋值 template<typename U> Matrix(initializer_list<U>) = delete; // 禁止使用 {} 除元素初始化外的初始化 template<typename U> Matrix& operator=(initializer_list<U>) = delete; static constexpr size_t order() { return N; } // 维度数 size_t extent(size_t n) const { return desc.extents[n]; } // 第 n 维的元素数量 size_t size() const { return elems.size(); } // 总元素数量 const Matrix_slice<N>& descriptor() const { return desc; } // 切片描述符 T* data() { return elems.data(); } // 元素数据指针 const T* data() const { return elems.data(); } // ... private: Matrix_slice<N> desc; // 切片描述符 vector<T> elems; // 元素 }; ``` 使用 `vector<T>` 存储元素可以避免内存管理和异常安全问题。`Matrix_slice` 用于定义矩阵的形状和下标映射，`Matrix_ref` 用于表示子矩阵的引用，`Matrix_initializer` 是一个嵌套的初始化列表。 ### 2.1 构造和赋值 `Matrix` 的默认拷贝和移动操作具有正确的语义，即对 `desc` 和元素进行逐成员拷贝或移动。以下是构造函数的实现： ```cpp // 接受维度大小的构造函数 template<typename T, size_t N> template<typename... Exts> Matrix<T, N>::Matrix(Exts... exts) : desc{exts...}, // 拷贝维度大小 elems(desc.size) // 分配并初始化元素 { } // 接受初始化列表的构造函数 template<typename T, size_t N> Matrix<T, N>::Matrix(Matrix_initializer<T, N> init) { Matrix_impl::derive_extents(init, desc.extents); // 从初始化列表推导维度大小 elems.reserve(desc.size); // 预留空间 Matrix_impl::insert_flat(init, elems); // 从初始化列表插入元素 assert(elems.size() == desc.size); } ``` 为了确保 `{}` 仅用于元素列表的初始化，我们删除了简单的初始化列表构造函数。同时，我们还支持从 `Matrix_ref` 构造和赋值： ```cpp template<typename T, size_t N> template<typename U> Matrix<T, N>::Matrix(const Matrix_ref<U, N>& x) : desc{x.desc}, elems{x.begin(), x.end()} // 拷贝描述符和元素 { static_assert(Convertible<U, T>(), "Matrix constructor: incompatible element types"); } template<typename T, size_t N> template<typename U> Matrix<T, N>& Matrix<T, N>::operator=(const Matrix_ref<U, N>& x) { static_assert(Convertible<U, T>(), "Matrix =: incompatible element types"); desc = x.desc; elems.assign(x.begin(), x.end()); return *this; } ``` ### 2.2 下标和切片 矩阵可以通过下标（元素或行）、行和列、切片（行或列的部分）进行访问。以下是一些访问方式： | 访问方式 | 描述 | | ---- | ---- | | `m.row(i)` | 矩阵 `m` 的第 `i` 行，返回 `Matrix_ref<T, N - 1>` | | `m.column(i)` | 矩阵 `m` 的第 `i` 列，返回 `Matrix_ref<T, N - 1>` | | `m[i]` | C 风格下标访问，等同于 `m.row(i)` | | `m(i, j)` | Fortran 风格元素访问，返回 `T&`，下标数量必须为 `N` | | `m(slice(i, n), slice(j))` | 切片访问子矩阵，返回 `Matrix_ref<T, N>` | 这些访问操作都是成员函数，以下是部分实现： ```cpp template<typename T, size_t N> Matrix_ref<T, N - 1> Matrix<T, N>::operator[](size_t n) { return row(n); } template<typename T, size_t N> template<typename... Args> Enable_if<Matrix_impl::Requesting_element<Args...>(), T&> Matrix<T, N>::operator()(Args... args) { assert(Matrix_impl::check_bounds(desc, args...)); return *(data() + desc(args...)); } ``` 切片操作使用 `slice` 结构体，它描述了从整数下标到元素位置的映射： ```cpp struct slice { slice() : start(-1), length(-1), stride(1) { } explicit slice(size_t s) : start(s), length(-1), stride(1) { } slice(size_t s, size_t l, size_t n = 1) : start(s), length(l), stride(n) { } size_t operator()(size_t i) const { return start + i * stride; } static slice all; size_t start; // 起始索引 size_t length; // ```