对称块矩阵乘法

Question

我试图乘以两个块对称矩阵（矩阵大小矩阵大小）。我想执行块矩阵乘法（将一个矩阵分成多个块大小矩阵，并将相应的块相乘）。我已经写了一些代码，但想改进它，并存储主对角线以上的块，但我没有任何想法。如果可能的话，你们能帮忙吗？

#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是你的块大小。

Answer 2

你能使用现有的线性代数软件包吗？如果你要对付灵长类动物，比如doubleBLAS可能是最理想的方法，但可能会有一个陡峭的学习曲线。对于高度优化但非常用户友好的库，Eigen是我最喜欢的c语言任务之一。

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

虽然我强烈建议你不要自己动手，但如果你绝对必须自己动手，有一个问题是，所有的积木大小都一样吗？矩阵是以什么形式存储的，主要是列还是行？假设块大小相同，并且有行主形式，可以做的一个草图是迭代块，并将块乘法降级为通用矩阵乘法函数。我要把加倍*

编辑：如果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;
}

我非常强调我建议你放弃所有这些，花点时间学习线性代数软件包。你将摆脱很多技术问题（例如：我做的指针算术正确吗？）你可以用这个软件包完成更多的任务。我认为这将有利于你的整体工作。