向量和矩阵的表达式模板

不是我折腾矩阵和向量, 是矩阵和向量折腾我.
早期为速度, 使用 cblas 比较多. cblas 非常不好用, double, float 函数各一 还好说, 关键是大部分函数的参数太多了, 如

void cblas_dgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA,
                 const enum CBLAS_TRANSPOSE TransB, const int M, const int N,
                 const int K, const double alpha, const double *A,
                 const int lda, const double *B, const int ldb,
                 const double beta, double *C, const int ldc);

基本上, 每次使用, 我都测试一下 M, N, K 几个值差得比较大的情形, 看看有没 有段错误, 这个方法比认真检查参数还要管用.
后来, 用 C++ 比较多, 仿照 STL, 希望尽量泛化, 不要 double, float 各写一 遍. 还希望进一步简化, 对所有矩阵, dense 的也好, csr, csc 也好, 我 们应该只需要写 gemm(A, B).
还要将计算和存储分开, 比如 gemm(A, B) 只表示一个表达式, 并不做任何实 际运算, 需要时, 可以将结果并行地存到另一个区域. 这样, 比如写 C=gemm(A, B) * alpha + beta * C 就比上面的 cblas 简明易懂多了. 还可以 用 C.parallel_assign(gemm(A, B) * alpha + beta * C) 表示并行地计算
$$C = \alpha A B + \beta C.$$
要想实现上面的表达式, 我们需要一个 proxy 抽象, vector_proxy, matrix_proxy. 矩阵, 向量的所有公有操作, 在 proxy 里实现. 具体的向量, 矩阵继承此 proxy 即可.
对所有 vector_proxy, 定义

template <class Derived, class Value, class Reference>
struct vector_proxy {
...
size_t size() const;
value_type get(size_t idx) const;
reference_type get(size_t idx);
...
};

其中, Derived 为具体的向量, 需要定义 get, size, 并继承 vector_proxy<Derived>.
再为 vector_proxy 定义 加, 减 等基本运算, 结果也为一个 proxy, 此 proxy 的 get(size_t idx) 函数实现计算.
从 proxy 到 proxy 的赋值如下.

template <class Derived, class Value, class Reference>
struct vector_proxy {
...
 
template <class R>
void assign(const vector_proxy<R> &r) {
  size_t n = size();
  assert(r.size() == n);
  for (size_t i = 0; i < n; ++i) {
     get(i) = r.get(i);
  }
}
 
template <class R>
void parallel_assign(const vector_proxy<R> &r) {
  size_t n = size();
  assert(r.size() == n);
#pragma omp parallel for 
  for (size_t i = 0; i < n; ++i) {
    get(i) = r.get(i);
  }
}
 
...
 
}

matrix_proxy 和 vector_proxy 完全类似.
具体到 gemm_proxy 上, 如下实现 get 方法即可.

value_type get(size_t i, size_t j) {
  return vec_inner_prod(m_A.row(i), m_B.col(j));
}

矩阵的 row, col 返回的也是 vector_proxy. 在 vec_inner_prod 里, 根据向量是 dense 还是 sparse, 采用不同的计算方法即可.
gemm_proxy 里的 row(i), col(j) 可以采用矩阵和向量的乘积 gemv(m_B, m_A.row(i)), gemv(m_A, m_B.col(j)). 在 gemv_proxy 里, 还是用 vec_inner_prod 来自动区分 dense/sparse 的向量.
gemm(A, B) 在 A, B 都是稀疏矩阵时, 情况比较复杂, 使 gemv 返回一个稀 疏向量的实现效率比较低, 我直接将结果当成 dense 了.