计算科学 - 矩阵的平行组装 - 吾爱随笔录

矩阵的平行组装

计算科学矩阵并行计算布拉斯

2021-11-30 02:55:42

我有一个我想快速组装的矩阵，它是块状的：

A = (\begin{matrix} A_{11} & A_{12} & A_{13} \\ A_{21} & A_{22} & A_{23} \\ A_{31} & A_{32} & A_{33} \end{matrix})

$A = \pmatrix{ A_{11} & A_{12} & A_{13} \\ A_{21} & A_{22} & A_{23} \\ A_{31} & A_{32} & A_{33}}$ 和：

A_{i j} = \sum_{k = 1}^{N} α_{i j k} v_{k} v_{k}^{T}

$A_{ij} = \sum_{k=1}^N \alpha_{ijk} v_k v_k^T$ 而且我目前的实施速度很慢，所以我正在寻找有关如何加快速度的建议。对于给定的

k

$k$ ，矩阵

v_{k} v_{k}^{T}

$v_k v_k^T$ 出现在每个区块中

A_{i j}

$A_{ij}$ 所以它可以被预先计算，然后复制到每个块中。我使用 OpenMP 来分配计算

v_{k} v_{k}^{T}

$v_k v_k^T$ 跨多个线程，这非常快。代码的瓶颈是当我需要更新矩阵时

A

$A$ . 我给每个线程自己的矩阵 A 进行更新，然后在最后将它们全部相加以获得最终的 A。我的伪代码如下：

# pragma omp parallel for
loop from k = 1 to N
  v_k = buildv(k);   // v_k and G_k are stored in thread-specific storage
  G_k = v_k v_k^T;   // using DSYR - this step is fast
  loop i = 1 to 3    // this part is really SLOW
    loop j = 1 to 3
      alpha_ijk = ...; // this is fast
      A[thread_id](i,j) += alpha_ijk * G_k; // I wrote a matrix DAXPY routine for this

sum A[thread_id] over all threads to build final A

当我写代码的时候，我以为计算 $v_k$ 和 $G_k$ 将是缓慢的部分，因为它们涉及计算（BLAS 调用）。但是我很惊讶缓慢的部分只是将已经计算的结果复制到矩阵 $A$ （即使每个线程都有自己的 $A$ ）。所以我想我需要专注于以某种方式并行复制操作。有没有人有什么建议？

注意 $N$ 相当大（大约 500 万），这就是为什么我专注于并行计算 $v_k v_k^T$ .

编辑

根据要求，我编写了一个最小的工作示例（可在此处获得）。我编译了下面的代码：

gcc -g -Wall -O2 -fopenmp -o matrix matrix.c -lm -lgsl -lcblas -latlas

请注意，运行此程序的先决条件是 GSL 库和 ATLAS BLAS - 尽管您可以将 ATLAS 替换为任何 cblas 库（例如 -lgslcblas）。

我添加了一个#if 标志来启用/禁用矩阵复制操作，这是缓慢的部分。当矩阵被复制到更大的线程特定矩阵中时，我机器上的总运行时间是81 seconds。禁用此代码时，总运行时间为4.1 秒。我的机器上有 24 个内核。

为简单起见，我使用随机数生成器来指定向量 $v_k$ .

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>
#include <sys/time.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_rng.h>

static void
random_vector(const double a, const double b, gsl_vector * v, gsl_rng * r)
{
  size_t i;

  for (i = 0; i < v->size; ++i)
    {
      double vi = (b - a) * gsl_rng_uniform(r) + a;
      gsl_vector_set(v, i, vi);
    }
}

static int
matrix_daxpy(const double alpha, const gsl_matrix * X, gsl_matrix * Y)
{
  const size_t M = Y->size1;
  const size_t N = Y->size2;

  if (X->size1 != M || X->size2 != N)
    {
      GSL_ERROR ("matrices must have same dimensions", GSL_EBADLEN);
    }
  else
    {
      if (alpha != 0.0)
        {
          const size_t tda_X = X->tda;
          const size_t tda_Y = Y->tda;
          size_t i, j;

          for (i = 0; i < M; ++i)
            {
              for (j = 0; j < N; ++j)
                {
                  Y->data[tda_Y * i + j] += alpha * X->data[tda_X * i + j];
                }
            }
        }

      return GSL_SUCCESS;
    }
}

int
main(int argc, char *argv[])
{
  const size_t max_threads = (size_t) omp_get_max_threads();
  const size_t N = 500000;  /* my real N is around 5 million */
  const size_t m = 300;
  const size_t nblocks = 3; /* number of block matrices in A */
  const double a = -10.0;
  const double b = 10.0;
  gsl_vector **v = malloc(max_threads * sizeof(gsl_vector *));
  gsl_matrix **GTG = malloc(max_threads * sizeof(gsl_matrix *));
  gsl_rng **r = malloc(max_threads * sizeof(gsl_rng *));
  gsl_matrix **A = malloc(max_threads * sizeof(gsl_matrix *));
  size_t k;
  struct timeval tv0, tv1;

  for (k = 0; k < max_threads; ++k)
    {
      v[k] = gsl_vector_alloc(m);
      GTG[k] = gsl_matrix_calloc(m, m);
      A[k] = gsl_matrix_calloc(nblocks * m, nblocks * m);
      r[k] = gsl_rng_alloc(gsl_rng_default);
    }

  gettimeofday(&tv0, NULL);

#pragma omp parallel for
  for (k = 0; k < N; ++k)
    {
      int thread_id = omp_get_thread_num();
      size_t i, j;

      /* construct v_k */
      random_vector(a, b, v[thread_id], r[thread_id]);

      /* GTG += v_k v_k^T */
      gsl_blas_dsyr(CblasLower, 1.0, v[thread_id], GTG[thread_id]);

      /* copy lower triangle to upper triangle */
      gsl_matrix_transpose_tricpy('L', 0, GTG[thread_id], GTG[thread_id]);

#if 1
      /* SLOW PART: now copy alpha_{ijk} v_k v_k^T into thread-specific final A */
      for (i = 0; i < nblocks; ++i)
        {
          for (j = 0; j <= i; ++j)
            {
              double alpha_ijk = 1.0; /* for simplicity */
              gsl_matrix_view Av = gsl_matrix_submatrix(A[thread_id], i * m, j * m, m, m); /* A_{ij} m-by-m block */

              /* copy alpha_ijk v_k v_k^T into appropriate spot */
              matrix_daxpy(alpha_ijk, GTG[thread_id], &Av.matrix);
            }
        }
#endif
    }

  gettimeofday(&tv1, NULL);
  fprintf(stderr, "elapsed time = %.4f [sec]\n",
          (tv1.tv_sec + tv1.tv_usec * 1.0e-6) - (tv0.tv_sec + tv0.tv_usec * 1.0e-6));

  return 0;
}

2个回答

第一个建议（不看代码）是检查您的内存访问模式是否合理，或者最好是对代码的关键部分进行优化。

需要牢记数据是如何存储在 RAM 中的，并在操作过程中最大限度地减少缓存未命中。在你的情况下，你应该检查

无论你的大矩阵 $A$ , 它的子块 $A_{ij}$ 按列或按行存储。无论它们以何种方式存储，都应该反映在它们被访问以进行读/写操作的方式上。
尽量避免DSYR和DAXPY来电有INCX和INCY不同1。如果存储间距增量不同于 1（非连续内存访问模式），则重新组织数据可能会有所帮助。这不太可能是问题，DSYR因为您提到它很快。
您提到您已经DAXPY为矩阵编写了自己的代码。我不确定您是否需要自己的一个，因为您可以调用 regular DAXPY，只是N现在不是向量的长度，而是矩阵的大小。您的自定义实现DAXPY很可能会减慢您的计算速度。
更改操作顺序也可能有所帮助。现在，你已经完成了循环 $k$ 在上面。我会尝试先循环块，然后有一个 OpenMP 并行循环 $k$ 里面。块的数量非常少，并且预计算 $G_k$ 和 $v_k$ （由 9 个块共享）可以交易数据局部性。
您可以尝试手动展开“块循环” $i,j$ . 很难说你是如何存储各个块的以及使用了哪些编译器选项，但是从二维索引改变 $i,j = 1,\ldots,3$ 一维索引 $b=1,\ldots,9$ 可以是有益的。它肯定会毫不费力地删除一层折叠循环。

你在更新 $m^2\times\mathit{nblocks}^2$ 每次迭代的内存位置，每个 8 字节，总共 $N$ 次，所以你得到了内存带宽 $N\times m^2\times\mathit{nblocks}^2\times 8\mathrm{B}\,/\,81\mathrm{s} = 40\mathrm{GB}/\mathrm{s}$ . （当您禁用慢速部分时，它的公式与 $\mathit{nblocks}=1$ , $t=4.1\mathrm{s}$ , 给 $82\mathrm{GB}/\mathrm{s}$ .) 这不算更新到 $G$ ，所以可能会偏离一个小因素，但这并不重要。在我看来，您的代码可能是内存有限的，并且您已经在最大可能内存带宽的一个小整数因子内（无论您的机器上是什么）。

在这种情况下，一种可能的解决方案是进一步减少矩阵块，直到它们适合缓存，并特别注意仅操作已经在缓存中的值。gemm这与典型例程对矩阵块所做的非常相似。所以让块更小，这样 $8\mathrm{B}\times m^2$ 适合您的每线程缓存：如果它是 256KB，那就是 $\sqrt{256\times 2^{10}/8}\sim 180$ 给予或接受（直接尝试手动设置此值而不是尝试自己计算它通常更容易）。对于每个这样的块，你仍然会累积 $\alpha_{ijk}v_k v_k^\top$ 但只使用每个的一部分 $v_k$ ，取决于块，所以它ger是非对称 rank-1 更新。注意：我没有在代码中测试过这个想法。

PS我认为你也可以通过只在最后进行三角转置复制来节省两倍，每个块一次，而不是每次迭代。

其它你可能感兴趣的问题

上一篇是否存在关于计算科学实验如何声称有效性的“框架”？计算科学的科学方法？下一篇牛顿法的计算复杂度