对称块矩阵乘法答案

【问题标题】：Symmetric Block Matrix Multiplication对称块矩阵乘法
【发布时间】：2023-03-14 14:53:01
【问题描述】：

我正在尝试将两个块对称矩阵相乘 (MATRIX_SIZExMATRIX_SIZE)。我想执行块矩阵乘法（将矩阵划分为多个 BLOCK_SIZExBLOCK_SIZE 矩阵并乘以相应的块）。我写了一些代码，但想改进它并存储主对角线上方的块，但我没有任何想法。可以的话请大家帮忙看看？

#define IND(A, x, y) A[y*MATRIX_SIZE+x]
void block_mult2(double*& A, double*& B, double*& C){
int i, j, k, i0, j0, k0;
for (i = 0; i < MATRIX_SIZE; i += BLOCK_SIZE)
for (j = 0; j < MATRIX_SIZE; j += BLOCK_SIZE)
for (k = 0; k < MATRIX_SIZE; k += BLOCK_SIZE)
    for (i0 = i; i0 < min(BLOCK_SIZE+i, MATRIX_SIZE); i0++)
        for (j0 = j; j0 < min(BLOCK_SIZE+j, MATRIX_SIZE); j0++)
            for (k0 = k; k0 < min(BLOCK_SIZE+k, MATRIX_SIZE); k0++)
                IND(C, i0, j0) += IND(A, i0, k0) * IND(B, k0, j0);
}

【问题讨论】：

标签： c++ algorithm matrix matrix-multiplication

【解决方案1】：

您可以使用现有的线性代数包吗？如果您正在处理原始类型，例如double BLAS 可能是最好的方法，但可能会有陡峭的学习曲线。对于高度优化但对用户非常友好的库，Eigen 是我在 c++ 中执行此类任务时最喜欢的选项之一。

我强烈建议使用现有的线性代数包（甚至不一定是我提到的那些）。由于实际的实现是从包中处理的，这将使您更容易充实您的想法。更不用说这样的软件包已经存在多年（在 BLAS 的情况下是几十年），并且应该非常非常非常擅长此类任务。除非你真的知道你在做什么（有一个非常非常具体的任务，你可以编码的特定优化）我怀疑你可以轻松地优化这些库（如果有的话）。即便如此，也需要考虑成本效益分析：与现有的好软件包相比，我自己做这件事要花多少时间？

虽然我强烈建议您不要自己做，但如果您绝对必须自己做，一个不清楚的问题是所有块的大小都相同吗？矩阵还有什么形式存储，列或行主要？假设块大小相同，并且您具有行主要形式，您可以做的草图是迭代块并将块-块乘法归为通用矩阵乘法函数。我正在删除double*& 并仅传递指针double*。 operator[] 应该注意引用正确的位置，但请检查我在 [] 中的算术是否正确：

编辑：如果 A 和 B 只存储上三角块，我更正了代码

//Assuming all blocks are the same size
//Assuming matrix stored in row major form

#define NUMBER_OF_BLOCKS = MATRIX_SIZE/BLOCK_SIZE

void block_mult2(double* A, double* B, double* C){
  for(size_t i=0; i<NUMBER_OF_BLOCKS; i++)
    for(size_t j=0; j<NUMBER_OF_BLOCKS; j++)
      for(size_t k=0; k<NUMBER_OF_BLOCKS; k++)
        mult2(A[min(i,j)*BLOCK_SIZE*NUMBER_OF_BLOCKS + max(i,j)*BLOCK_SIZE],
              B[min(j,k)*BLOCK_SIZE*NUMBER_OF_BLOCKS + max(j,k)*BLOCK_SIZE],
              C[i*BLOCK_SIZE*NUMBER_OF_BLOCKS + k*BLOCK_SIZE]);
  return;
}

void mult2(double* A, double* B, double* C){
  for(size_t i=0; i<BLOCK_SIZE; i++)
    for(size_t j=0; j<BLOCK_SIZE; j++)
      for(size_t k=0; k<BLOCK_SIZE; k++)
        C[i*BLOCK_SIZE+k] = A[min(i,j)*BLOCK_SIZE+max(i,j)]*B[min(j,k)*BLOCK_SIZE+max(j,k)];
  return;
}

我怎么强调都不为过，我建议你放弃所有这些，花点时间学习线性代数包。您将摆脱很多技术问题（例如刚刚提出的：我是否正确地进行了指针运算？）并且您可以使用该包来完成更多任务。我认为这将有利于您的整体工作。

【讨论】：

【解决方案2】：

for(int jj=0;jj<N;jj+= s){
    for(int kk=0;kk<N;kk+= s){
            for(int i=0;i<N;i++){
                    for(int j = jj; j<((jj+s)>N?N:(jj+s)); j++){
                            temp = 0;
                            for(int k = kk; k<((kk+s)>N?N:(kk+s)); k++){
                                    temp += a[i][k]*b[k][j];
                            }
                            c[i][j] += temp;
                    }
            }
     }
 }

我很抱歉这个虚拟代码，但你可以考虑 N 是你的 BLOCK_SIZE

【讨论】：