对称块矩阵乘法

Question

我正在尝试将两个块对称矩阵相乘 (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);
}

Answer 1

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

Answer 2

你能使用现有的线性代数包吗？如果您正在处理原始类型，例如 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;
}

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

对称块矩阵乘法

Symmetric Block Matrix Multiplication

c++

algorithm

matrix

matrix-multiplication