title: OpenCL Matrix Multiplication on Radeon HD 7970 and Intel HD 4000 GPUs vs Eigen on CPU
body: <p>I am having some problems understanding the outcome of my comparison.</p>
  
  <p>Laptop with i7/intel hd4000 and server with 8Xeon 5400/7970 HDRadeon.</p>
  
  <p>I am multiplying:</p>
  
  <pre><code>int M =1024*2, N = 1024*6, P = 1024*2;
      //      N               P      
      //|-----------|   |-----------|
      //|           |   |           |
      //|M          | * |N          |
      //|           |   |           |
      //|-----------|   |-----------|
  </code></pre>
  
  <p>This is the kernel:</p>
  
  <pre><code>/*
   * Copyright 1993-2010 NVIDIA Corporation.  All rights reserved.
   *
   * Please refer to the NVIDIA end user license agreement (EULA) associated
   * with this source code for terms and conditions that govern your use of
   * this software. Any use, reproduction, disclosure, or distribution of
   * this software and related documentation outside the terms of the EULA
   * is strictly prohibited.
   *
   */
  
  /* Matrix multiplication: C = A * B.
   * Device code.
   */
  #ifndef BLOCK_SIZE
  #define BLOCK_SIZE 16
  #endif
  
  #define AS(i, j) As[j + i * BLOCK_SIZE]
  #define BS(i, j) Bs[j + i * BLOCK_SIZE]
  
  ///////////////////////////////////////////////////////////////////////////////
  //! Matrix multiplication on the device: C = A * B
  //! uiWA is A's width and uiWB is B's width
  ////////////////////////////////////////////////////////////////////////////////
  __kernel void
  m_m_mul( __global float* A, __global float* B, __global float* C, 
         /*__local float* As, __local float* Bs,*/ int uiWA, int uiWB, int trueLocalSize1)
  {
  
     __local float As[BLOCK_SIZE*BLOCK_SIZE];
     __local float Bs[BLOCK_SIZE*BLOCK_SIZE];
  
      // Block index
      int bx = get_group_id(0);
      int by = get_group_id(1);
  
      // Thread index
      int tx = get_local_id(0);
      int ty = get_local_id(1);
  
      // Index of the first sub-matrix of A processed by the block
      int aBegin = uiWA * BLOCK_SIZE * by;
  
      // Index of the last sub-matrix of A processed by the block
      int aEnd   = aBegin + uiWA - 1;
  
      // Step size used to iterate through the sub-matrices of A
      int aStep  = BLOCK_SIZE;
  
      // Index of the first sub-matrix of B processed by the block
      int bBegin = BLOCK_SIZE * bx;
  
      // Step size used to iterate through the sub-matrices of B
      int bStep  = BLOCK_SIZE * uiWB;
  
      // Csub is used to store the element of the block sub-matrix
      // that is computed by the thread
      float Csub = 0.0f;
  
      // Loop over all the sub-matrices of A and B
      // required to compute the block sub-matrix
      for (int a = aBegin, b = bBegin;
               a &lt;= aEnd;
               a += aStep, b += bStep) {
  
          // Load the matrices from device memory
          // to shared memory; each thread loads
          // one element of each matrix
          AS(ty, tx) = A[a + uiWA * ty + tx];
          BS(ty, tx) = B[b + uiWB * ty + tx];
  
          // Synchronize to make sure the matrices are loaded
          barrier(CLK_LOCAL_MEM_FENCE);
  
          // Multiply the two matrices together;
          // each thread computes one element
          // of the block sub-matrix        
          #pragma unroll
          for (int k = 0; k &lt; BLOCK_SIZE; ++k)
              Csub += AS(ty, k) * BS(k, tx);
  
          // Synchronize to make sure that the preceding
          // computation is done before loading two new
          // sub-matrices of A and B in the next iteration
          barrier(CLK_LOCAL_MEM_FENCE);
      }
  
      if (get_global_id(1) &lt; trueLocalSize1)
      // Write the block sub-matrix to device memory;
      // each thread writes one element
      C[get_global_id(1) * get_global_size(0) + get_global_id(0)] = Csub;
  
  }
  </code></pre>
  
  <p>I compare with <code>Eigen::Matrix&lt;float,-1,-1,Eigen::RowMajor&gt; m4 = m1 * m2</code>;</p>
  
  <p>On the server:</p>
  
  <pre><code>Creating matrices on GPU....... Done [0ms]
  Creating matrices on CPU....... Done [0ms]
  Filling GPU with random numbers....... Done [19ms]
  M3 = M1 * M2... on GPU (Loading Kernels)... Done [240ms]
  M3 = M1 * M2... on GPU (3 times)... Done [211ms]
  Loading M1, M2 on GPU... Done [93ms]
  M4 = M1 * M2 on CPU... Done [7775ms] Error:3.78049e-008
  Press any key to continue . . .
  Matlab: Elapsed time is 3.010626 seconds.
  </code></pre>
  
  <p>On the Laptop:</p>
  
  <pre><code>Creating matrices on GPU....... Done [22ms]
  Creating matrices on CPU....... Done [0ms]
  Filling GPU with random numbers....... Done [35ms]
  M3 = M1 * M2... on GPU (Loading Kernels)... Done [2975ms]
  M3 = M1 * M2... on GPU (3 times)... Done [6891ms]
  Loading M1, M2 on GPU... Done [80ms]
  M4 = M1 * M2 on CPU... Done [5966ms] Error:3.78049e-008
  Press any key to continue . . .
  Matlab: Elapsed time is 2.310626 seconds.
  </code></pre>
  
  <p>My questions are now.
  1) Why are the laptop faster with eigen then the 8 cores of xeon. Could it be that eigen only use one core on both systems and i7 has a higher clockspeed?2.0 vs 2.4?</p>
  
  <p>2) Almost 3 times speed up using Intel HD4000 vs Eigen on labtop, but Matlab takes 2.3sec for the same multiplication. Thats the same as the kernel on HD4000. (Can i do something to get Eigen to run at same speed as Matlab?)</p>