title: multi-thread matrix multiplication
body: <p>I'm writing a code that does N X N matrix multiplication using thread level parallelism.<br>
  <br>
  To get C = A X B, <br>
  first I transposed matrix B, divided matrices into blocks.<br></p>
  
  <p>A thread takes a block from A and B and multiply them,<br>
  and then adds the result to the corresponding block in C.<br>
  All matrices are allocated in heap memory using malloc().<br></p>
  
  <p>But the problem is that for the same input,<br>
  the answer is sometimes incorrect and sometimes correct.<br></p>
  
  <p>I'm not quite sure why this happens, <br>but I guess the code need to be improved in terms of thread safety.<br></p>
  
  <p>I post some part of my code.<br>
  <em>blocks</em>  is the number of blocks in row and column, i.e. N / block size.<br>
  So the total number of threads is blocks^3.</p>
  
  <pre><code>while (thread_loaded &lt; total_thread)
   {
    if (thread_count &lt; MAX_THREAD)
    {
     p[thread_loaded].idx = thread_idx;
     p[thread_loaded].matinfo = &amp;mi;   
     threads[thread_loaded] = CreateThread(NULL, 0, matmul_tlp, &amp;p[thread_loaded], 0, &amp;threadid[thread_loaded]);
     if (threads[thread_loaded] != NULL)
     {     
      thread_loaded++;
      thread_count++;
      thread_idx += blocks;
      if (thread_idx &gt;= total_thread)
       thread_idx = (thread_idx % total_thread) + 1;
     }          
    }
   }
  </code></pre>
  
  <p>for the thread function,</p>
  
  <pre><code>int i, j, k;  
  param* p = (param*)arg; 
  
  int blocks = BLOCKS;
  int block_size = BLOCK_SIZE;
  int Ar = p-&gt;idx / (blocks * blocks); 
  int Ac = p-&gt;idx % blocks;
  int Br = (p-&gt;idx / blocks) % blocks; 
  int Bc = p-&gt;idx % blocks;
  int Cr = p-&gt;idx / (blocks * blocks);
  int Cc = (p-&gt;idx % (blocks * blocks)) / blocks;
  
  double** A = p-&gt;matinfo-&gt;A;
  double** B = p-&gt;matinfo-&gt;B;
  double** C = p-&gt;matinfo-&gt;C;
  
  
  DWORD res = WaitForSingleObject(mutex[Cr * blocks + Cc], INFINITE);
  if (res != WAIT_OBJECT_0)
   perror("cannot acquire mutex.");
  
  
  for (i = 0; i &lt; block_size; i++)
  {
   for (j = 0; j &lt; block_size; j++)
   {
    for (k = 0; k &lt; block_size; k++)
    {
     C[Cr * block_size + i][Cc * block_size + j] +=
      A[Ar * block_size + i][Ac * block_size + k] *
      B[Br * block_size + j][Bc * block_size + k];
    }
   }
  }
  
  
  ReleaseMutex(mutex[Cr * blocks + Cc]);
  
  
  thread_count--;
  return NULL;
  </code></pre>
  
  <p>It would be appreciated if anyone can find a hole. :)</p>