C++ Divide and conquer square matrix multiplication problems

Question

I'm trying to do divide and conquer matrix multiplication so i can parallelize it, but I'm getting half random garbage numbers and half 0's in the result, e.g. on a 2x2 matrix "[[15909360,0][15909360,0]]". This is what I have so far based on the algorithm on Wikipedia but I don't really know where to go from here. I'm not using pointers yet for the partitioning or threads. This is homework btw.

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;
}

EDIT: I changed C[0][0] += C[0][0] + T[0][0]; in add() to C[0][0] += T[0][0]; Also I made a unpartition function that basically does the reverse and puts the partitions back into C and T after multiplying and adding:

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
 }
}

My vectors get initialized correctly after I fixed the default constructor for my IntElement class.

Did you try debugging your code? Perhaps do the steps by hand on paper and compare with the program? — Quimby, Mar 11 '19 at 21:43
FYI: you can make your code much more readable (to us, and to yourself) by aliasing your vector type with something like `using elem_mat_t = std::vector >;` at the top of your file. — , Mar 11 '19 at 21:45
Welcome to Stack Overflow. Please take a look at our page on [minimal complete examples](https://stackoverflow.com/help/mcve); what is the smallest, simplest multiplication example that produces the error? — Beta, Mar 11 '19 at 21:47
But just at first glance it looks like you never initialize those local vectors `T,a11,b11,c11,d11`. They are NOT zero-initialized. — Quimby, Mar 11 '19 at 21:47
I would suggest you to not use vectors of vectors as matrices because they result in very inefficient data organization. The best way probably is to write a wrapper class that stores the data in a 1D vector. — eesiraed, Mar 12 '19 at 05:08

score 0 · Answer 1 · answered Mar 11 '19 at 22:15

Your partition function makes copies of the data in the source vectors. In multiply and add, you final series of calls (multiply(c11, a11, b11, n/2), add(c11, t11, n/2)) will modify those submatrix objects. This does not modify the original C, or T matrices. T will stay filled with zeroes, and C will get occasional updates from 1x1 matrices. You need to "unpartition" the results back into the appropriate matrices.

And the single element case in add is wrong. C[0][0] += C[0][0] + T[0][0]; should either be C[0][0] += T[0][0]; or C[0][0] = C[0][0] + T[0][0]; (but not both).

Quimby · Answer 2 · 2019-03-11T22:14:56.493

-1

The inner vectors in variables T,a11,b11,c11,d11 are never initialized. They contain the garbage which is used in those 4 parititons calls.

EDIT: The above statement is false, as Alan Birtles pointed out they are initialized by default constructor of IntElement. I still believe that it might be just a typedef for int in which case my argument stands.

edited Mar 11 '19 at 22:14

answered Mar 11 '19 at 21:50

Quimby

17,735
4
35
55

Vector elements will be initialised by the default constructor of `IntElement` whatever type that is, presumably that does something sensible like initialise to zero – Alan Birtles Mar 11 '19 at 22:00
@AlanBirtles Yes, you are right, my mistake. If `IntElement` is a just typedef for `int` then it won't. I believe that might be the case so I will leave the answer here, but feel free to downvote as it's not correct. – Quimby Mar 11 '19 at 22:13
i tried initializing and i got a matrix with all zeroes so if it is needed there's still something else too. IntElement is a class that holds an int. – Sasha Mar 11 '19 at 22:28
it indeed isn't needed just make sure your classes default constructor assigns a default value – Sasha Mar 11 '19 at 23:10
Even a vector of ints will be initialised to 0, all elements of a vector will always be initialised – Alan Birtles Mar 12 '19 at 06:52

C++ Divide and conquer square matrix multiplication problems

2 Answers2