C++分而治之的方阵乘法问题

Posted

技术标签:

【中文标题】C++分而治之的方阵乘法问题【英文标题】:C++ Divide and conquer square matrix multiplication problems 【发布时间】:2019-08-02 07:02:20 【问题描述】:

我正在尝试做分而治之的矩阵乘法,这样我就可以并行化它,但是我得到了一半的随机垃圾数和一半的 0,例如在 2x2 矩阵“[[15909360,0][15909360,0]]”上。这是我迄今为止基于the algorithm on Wikipedia 所拥有的,但我真的不知道从这里去哪里。我还没有将指针用于分区或线程。顺便说一句,这是家庭作业。

void partition(const std::vector<std::vector<IntElement> >& m, std::vector<std::vector<IntElement> >& m11, std::vector<std::vector<IntElement> >& m12,
                 std::vector<std::vector<IntElement> >& m21, std::vector<std::vector<IntElement> >& m22, int n)
for(int i=0;i<n/2;i++)
 for(int j=0;j<n/2;j++)
  m11[i][j] = m[i][j]; // top left
  m12[i][j] = m[i][j + n / 2]; // top right
  m21[i][j] = m[i + n / 2][j]; // bottom left
  m22[i][j] = m[i + n / 2][j + n / 2]; // bottom right
 
;
void add(std::vector<std::vector<IntElement> >& C, std::vector<std::vector<IntElement> >& T, int n)
if(n==1)
  C[0][0] += C[0][0] + T[0][0];

else
std::vector<std::vector<IntElement> > c11(n/2, std::vector<IntElement>(n/2)), c12(n/2, std::vector<IntElement>(n/2)),
  c21(n/2, std::vector<IntElement>(n/2)), c22(n/2, std::vector<IntElement>(n/2));
std::vector<std::vector<IntElement> > t11(n/2, std::vector<IntElement>(n/2)), t12(n/2, std::vector<IntElement>(n/2)),
  t21(n/2, std::vector<IntElement>(n/2)), t22(n/2, std::vector<IntElement>(n/2));
partition(C, c11, c12, c21, c22, n);
partition(T, t11, t12, t21, t22, n);

add(c11, t11, n/2);
add(c12, t12, n/2);
add(c21, t21, n/2);
add(c22, t22, n/2);
  
;
void multiply(std::vector<std::vector<IntElement> >& C, const std::vector<std::vector<IntElement> >& A,
          const std::vector<std::vector<IntElement> >& B, int n)
if(n==1)
    C[0][0] += A[0][0] * B[0][0];
else
  std::vector<std::vector<IntElement> > T(n, std::vector<IntElement>(n));
  std::vector<std::vector<IntElement> > a11(n/2, std::vector<IntElement>(n/2)), a12(n/2, std::vector<IntElement>(n/2)),
    a21(n/2, std::vector<IntElement>(n/2)), a22(n/2, std::vector<IntElement>(n/2));
  std::vector<std::vector<IntElement> > b11(n/2, std::vector<IntElement>(n/2)), b12(n/2, std::vector<IntElement>(n/2)),
    b21(n/2, std::vector<IntElement>(n/2)), b22(n/2, std::vector<IntElement>(n/2));
  std::vector<std::vector<IntElement> > c11(n/2, std::vector<IntElement>(n/2)), c12(n/2, std::vector<IntElement>(n/2)),
    c21(n/2, std::vector<IntElement>(n/2)), c22(n/2, std::vector<IntElement>(n/2));
  std::vector<std::vector<IntElement> > t11(n/2, std::vector<IntElement>(n/2)), t12(n/2, std::vector<IntElement>(n/2)),
    t21(n/2, std::vector<IntElement>(n/2)), t22(n/2, std::vector<IntElement>(n/2));

  partition(A, a11, a12, a21, a22, n);
  partition(B, b11, b12, b21, b22, n);
  partition(C, c11, c12, c21, c22, n);
  partition(T, t11, t12, t21, t22, n);

  multiply(c11, a11, b11, n/2);
  multiply(c12, a11, b12, n/2);
  multiply(c21, a21, b11, n/2);
  multiply(c22, a21, b12, n/2);

  multiply(t11, a12, b21, n/2);
  multiply(t12, a12, b22, n/2);
  multiply(t21, a22, b21, n/2);
  multiply(t22, a22, b22, n/2);

  add(C, T, n);

return;
;
SquareMatrix& SquareMatrix::operator*=(const SquareMatrix& m)
  std::vector<std::vector<IntElement> > C(n, std::vector<IntElement>(n));
  multiply(C, elements, m.elements, n);
  elements = C;
  return *this;


SquareMatrix operator*(const SquareMatrix& a, const SquareMatrix& b)
  SquareMatrix c = a;
  c *= b;
  return c;

编辑:我改变了 C[0][0] += C[0][0] + T[0][0];在 add() 到 C[0][0] += T[0][0];我还做了一个 unpartition 函数,它基本上做相反的事情,并在相乘和相加后将分区放回 C 和 T:

void unpartition(std::vector<std::vector<IntElement> >& m,std::vector<std::vector<IntElement> >& m11, std::vector<std::vector<IntElement> >& m12,
           std::vector<std::vector<IntElement> >& m21, std::vector<std::vector<IntElement> >& m22, int n)
for(int i=0;i<n/2;i++)
for(int j=0;j<n/2;j++)
  m[i][j] = m11[i][j]; // top left
  m[i][j + n / 2] = m12[i][j]; // top right
  m[i + n / 2][j] = m21[i][j]; // bottom left
  m[i + n / 2][j + n / 2] = m22[i][j]; // bottom right
 

在我修复了 IntElement 类的默认构造函数后,我的向量被正确初始化。

【问题讨论】:

你试过调试你的代码吗?也许在纸上手工完成步骤并与程序进行比较? 仅供参考:您可以通过在文件顶部使用 using elem_mat_t = std::vector&lt;std::vector&lt;IntElement&gt; &gt;; 之类的矢量类型别名来使您的代码(对我们和您自己)更具可读性。 欢迎来到 Stack Overflow。请在minimal complete examples 上查看我们的页面;产生错误的最小、最简单的乘法示例是什么? 但乍一看,您似乎从未初始化那些本地向量s T,a11,b11,c11,d11。它们不是零初始化的。 我建议您不要将向量的向量用作矩阵,因为它们会导致数据组织效率非常低。最好的方法可能是编写一个将数据存储在一维向量中的包装类。 【参考方案1】:

您的partition 函数会复制源向量中的数据。在multiplyadd 中,您最后的一系列调用(multiply(c11, a11, b11, n/2)add(c11, t11, n/2))将修改这些子矩阵对象。这不会修改原始的CT 矩阵。 T 将保持填充零,C 将偶尔从 1x1 矩阵中获得更新。您需要将结果“取消分区”回适当的矩阵中。

add 中的单元素情况是错误的。 C[0][0] += C[0][0] + T[0][0]; 应该是 C[0][0] += T[0][0];C[0][0] = C[0][0] + T[0][0];(但不能同时是两者)。

【讨论】:

【参考方案2】:

变量T,a11,b11,c11,d11 中的内部向量永远不会被初始化。它们包含在这 4 个 parititons 调用中使用的垃圾。

编辑:上述陈述是错误的,正如 Alan Birtles 指出的那样,它们是由IntElement 的默认构造函数初始化的。我仍然相信它可能只是 int 的 typedef,在这种情况下我的论点成立。

【讨论】:

Vector 元素将由IntElement 的默认构造函数初始化,不管是什么类型,大概是做一些有意义的事情,比如初始化为零 @AlanBirtles 是的,你是对的,我的错。如果IntElement 只是int 的typedef,那么它不会。我相信可能是这种情况,所以我会在这里留下答案,但请随意投反对票,因为它不正确。 我尝试初始化,我得到了一个全为零的矩阵,所以如果需要它,还有其他东西。 IntElement 是一个包含 int 的类。 确实不需要确保你的类默认构造函数分配一个默认值 即使一个整数向量会被初始化为 0,一个向量的所有元素都会被初始化

以上是关于C++分而治之的方阵乘法问题的主要内容,如果未能解决你的问题,请参考以下文章

矩阵乘法 - 分而治之 vs Strassen,分而治之更快?

矩阵乘法 - 分而治与斯特拉森,分而治之更快?

分治法:多项式乘法时间复杂度

算法导论分治思想—大数乘法矩阵相乘残缺棋盘

C++探索之旅第一部分第七课:函数效应,分而治之

C++ 模板参数推导/替换失败: