title: My matrix multiplication doesn't have the same precision as BLAS
body: <p>I wrote a matrix multiplication routine in C++ and used it in an eigensolver. Then I realized that my implementation has precision problems even though I already use double-precision float-points. It seems the result of my implementation doesn't have enough precision. I tested it with BLAS and found its result doesn't match with BLAS's. The simple code below illustrates the problem and has the following output.</p>
  
  <pre><code>0: 0
  1: 0
  2: -1.42109e-14
  test: test.cpp:59: int main(): Assertion `C1[i] == C2[i]' failed.
  </code></pre>
  
  <p>I can simply call BLAS routines for matrix multiplication, but I want to know why this happens. I don't know if this is an architecture-related problem. I use Intel Xeon E5-4620.</p>
  
  <p>BTW, I tried to define sum with "long double" and still got different results.</p>
  
  <p>Thanks,
  Da</p>
  
  <pre><code>#include &lt;stdlib.h&gt;
  #include &lt;stdio.h&gt;
  #include &lt;stddef.h&gt;
  #include &lt;cblas.h&gt;
  #include &lt;assert.h&gt;
  
  #include &lt;vector&gt;
  
  void init_vec(std::vector&lt;double&gt; &amp;vec)
  {
      for (size_t i = 0; i &lt; vec.size(); i++)
          vec[i] = ((double) random()) / ((double) random());
  }
  
  void mv(const std::vector&lt;double&gt; &amp;A, const std::vector&lt;double&gt; &amp;B,
          std::vector&lt;double&gt; &amp;C)
  {
      size_t num_rows = A.size() / B.size();
      size_t num_cols = B.size();
      for (size_t i = 0; i &lt; num_rows; i++) {
          double sum = 0;
          for (size_t j = 0; j &lt; num_cols; j++)
              sum += A[j * num_rows + i] * B[j];
          C[i] = sum;
      }
  }
  
  int main()
  {
      std::vector&lt;double&gt; A(100 * 10);
      std::vector&lt;double&gt; B(10);
      std::vector&lt;double&gt; C1(100);
      std::vector&lt;double&gt; C2(100);
      init_vec(A);
      init_vec(B);
  
      mv(A, B, C1);
      cblas_dgemv(CblasColMajor, CblasNoTrans, 100, 10, 1, A.data(), 100, B.data(), 1, 0, C2.data(), 1);
      for (size_t i = 0; i &lt; C1.size(); i++) {
          printf("%ld: %g\n", i, C1[i] - C2[i]);
          assert(C1[i] == C2[i]);
      }
  }
  </code></pre>