title: What is the best way to implement an array of 3d vectors?
body: <p>I decided to use <a href="http://eigen.tuxfamily.org" rel="noreferrer">Eigen</a> library in my project.
  But it is not clear from documentation how the most efficiently one should specify an array of 3d vectors.</p>
  
  <p>As I suggest, the first way is  </p>
  
  <pre><code>Eigen::Matrix&lt;Eigen::Vector3d, Eigen::Dynamic, 1&gt; array_of_v3d(size);  
  </code></pre>
  
  <p>But in that case how should I get another array which elements are equal to scalar products of elements of <code>array_of_v3d</code> and some other instance of <code>Vector3d</code>?
  In other words, can I rewrite the following loop using <code>Eigen</code>'s functions:</p>
  
  <pre><code>Eigen::Vector3d v3d(0.5, 0.5, 0.5);  
  Eigen::VectorXd other_array(size);  
  for (size_t i = 0; i &lt; size; ++i)
      other_array(i) = array_of_v3d(i).dot(v3d);
  </code></pre>
  
  <p>The second way is to use a matrix which size is <code>(3 x size)</code> or <code>(size x 3)</code>.
  For example, I can declare it like this:</p>
  
  <pre><code>Eigen::Matrix&lt;double, 3, Eigen::Dynamic&gt; matrix;  
  </code></pre>
  
  <p>But I didn't get from the documentation how to set the number of columns.
  The following seems to work but I have to repeat the number of rows <code>3</code> twice:</p>
  
  <pre><code>Eigen::Matrix&lt;double, 3, Eigen::Dynamic&gt; matrix(3, size);  
  </code></pre>
  
  <p>Then the above loop is equivalent to</p>
  
  <pre><code>other_array = v3d.transpose() * array_of_v3d;
  </code></pre>
  
  <p>As my experiments shown that this is a little bit faster than </p>
  
  <pre><code>Eigen::Matrix&lt;double, Eigen::Dynamic, 3&gt; matrix(size, 3);  
  other_array = array_of_v3d * v3d;
  </code></pre>
  
  <p>In addition:</p>
  
  <p>Anyway, my use of <code>Eigen</code> seems to be not so optimal since the same program in plain <code>C</code> is almost 1.5 times faster (it is not true indeed, it depends on <code>size</code>): </p>
  
  <pre><code>for (size_t i = 0; i &lt; size; i+=4) {
      s[i]   += a[i]   * x + b[i]   * y  + c[i]   * z;
      s[i+1] += a[i+1] * x + b[i+1] * y  + c[i+1] * z;
      s[i+2] += a[i+2] * x + b[i+2] * y  + c[i+2] * z;
      s[i+3] += a[i+3] * x + b[i+3] * y  + c[i+3] * z;
  }
  </code></pre>
  
  <p>Have I missed something? Is there some other way in the scope of <code>Eigen</code> library to solve my problem?</p>
  
  <p>Update:</p>
  
  <p>Here I present results of my testings. There are 5 cases:</p>
  
  <ol>
  <li><code>C</code>-style for loop which is partially unrolled</li>
  <li><code>Eigen::Matrix</code> (<code>rows x cols = 3 x size</code>). In this case values of 3d vectors are stored together in memory because by default <code>Eigen</code> stores data in column-major. Or I could set <code>Eigen::RowMajor</code> and everything else like in next case.  </li>
  <li><code>Eigen::Matrix</code> (<code>rows x cols = size x 3</code>). </li>
  <li>Here each component of 3d vector is stored on a separate <code>VectorXd</code>. So there are three VectorXd objects which are put together on <code>Vector3d</code>.</li>
  <li>An <code>std::vector</code> container is used to store <code>Vector3d</code> objects.</li>
  </ol>
  
  <p>This is my test program</p>
  
  <pre><code>#include &lt;iostream&gt;
  #include &lt;iomanip&gt;
  #include &lt;vector&gt;
  #include &lt;ctime&gt;
  
  #include &lt;Eigen/Dense&gt;
  
  double for_loop(size_t rep, size_t size)
  {
      std::vector&lt;double&gt;  a(size), b(size), c(size);
      double x = 1, y = 2, z = - 3;
      std::vector&lt;double&gt;  s(size);
      for(size_t i = 0; i &lt; size; ++i) {
          a[i] = i;
          b[i] = i;
          c[i] = i;
          s[i] = 0;
      }
      double dtime = clock();
      for(size_t j = 0; j &lt; rep; j++) 
          for(size_t i = 0; i &lt; size; i += 8) {
              s[i] += a[i]   * x + b[i]   * y  + c[i]   * z;
              s[i] += a[i+1] * x + b[i+1] * y  + c[i+1] * z;
              s[i] += a[i+2] * x + b[i+2] * y  + c[i+2] * z;
              s[i] += a[i+3] * x + b[i+3] * y  + c[i+3] * z;
              s[i] += a[i+4] * x + b[i+4] * y  + c[i+4] * z;
              s[i] += a[i+5] * x + b[i+5] * y  + c[i+5] * z;
              s[i] += a[i+6] * x + b[i+6] * y  + c[i+6] * z;
              s[i] += a[i+7] * x + b[i+7] * y  + c[i+7] * z;
          }
      dtime = (clock() - dtime) / CLOCKS_PER_SEC;
  
      double res = 0;
      for(size_t i = 0; i &lt; size; ++i)
          res += std::abs(s[i]);
      assert(res == 0.);
  
      return dtime;
  }
  
  double eigen_3_size(size_t rep, size_t size)
  {
      Eigen::Matrix&lt;double, 3, Eigen::Dynamic&gt; A(3, size);
      Eigen::Matrix&lt;double, 1, Eigen::Dynamic&gt; S(size);
      Eigen::Vector3d X(1, 2, -3);
      for(size_t i = 0; i &lt; size; ++i) {
          A(0, i) = i;
          A(1, i) = i;
          A(2, i) = i;
          S(i)    = 0;
      }
      double dtime = clock();
      for (size_t j = 0; j &lt; rep; j++) 
          S.noalias() += X.transpose() * A;
      dtime = (clock() - dtime) / CLOCKS_PER_SEC;
  
      double res = S.array().abs().sum();
      assert(res == 0.);
  
      return dtime;
  }
  
  double eigen_size_3(size_t rep, size_t size)
  {
      Eigen::Matrix&lt;double, Eigen::Dynamic, 3&gt; A(size, 3);
      Eigen::Matrix&lt;double, Eigen::Dynamic, 1&gt; S(size);
      Eigen::Vector3d X(1, 2, -3);
      for(size_t i = 0; i &lt; size; ++i) {
          A(i, 0) = i;
          A(i, 1) = i;
          A(i, 2) = i;
          S(i)    = 0;
      }
      double dtime = clock();
      for (size_t j = 0; j &lt; rep; j++) 
          S.noalias() += A * X;
      dtime = (clock() - dtime) / CLOCKS_PER_SEC;
  
      double res = S.array().abs().sum();
      assert(res == 0.);
  
      return dtime;
  }
  
  double eigen_vector3_vector(size_t rep, size_t size)
  {
      Eigen::Matrix&lt;Eigen::VectorXd, 3, 1&gt; A;
      A(0).resize(size);
      A(1).resize(size);
      A(2).resize(size);
      Eigen::VectorXd S(size);
      Eigen::Vector3d X(1, 2, -3);
      for(size_t i = 0; i &lt; size; ++i) {
          A(0)(i) = i;
          A(1)(i) = i;
          A(2)(i) = i;
          S(i)    = 0;
      }
      double dtime = clock();
      for (size_t j = 0; j &lt; rep; j++) 
          S.noalias() += A(0) * X(0) + A(1) * X(1) + A(2) * X(2);
      dtime = (clock() - dtime) / CLOCKS_PER_SEC;
  
      double res = S.array().abs().sum();
      assert(res == 0.);
  
      return dtime;
  }
  
  double eigen_stlvector_vector3(size_t rep, size_t size)
  {
      std::vector&lt;                 Eigen::Vector3d,
          Eigen::aligned_allocator&lt;Eigen::Vector3d&gt; &gt; A(size);
      std::vector&lt;double&gt; S(size);
      Eigen::Vector3d X(1, 2, -3);
      for(size_t i = 0; i &lt; size; ++i) {
          A[i](0) = i;
          A[i](1) = i;
          A[i](2) = i;
          S[i]    = 0;
      }
      double dtime = clock();
      for (size_t j = 0; j &lt; rep; j++) 
          for(size_t i = 0; i &lt; size; i++) 
              S[i] += X.dot(A[i]);
      dtime = (clock() - dtime) / CLOCKS_PER_SEC;
  
      double res = 0;
      for(size_t i = 0; i &lt; size; i++) 
          res += std::abs(S[i]);
      assert(res == 0.);
  
      return dtime;
  }
  
  int main()
  {
      std::cout &lt;&lt; "    size    |  for loop  |   Matrix   |   Matrix   | Vector3 of  | STL vector of \n" 
                &lt;&lt; "            |            |  3 x size  |  size x 3  | Vector_size |  TinyVectors  \n" 
                &lt;&lt; std::endl;
  
      size_t n       = 10;
      size_t sizes[] = {16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192};
      int rep_all    = 1024 * 1024 * 1024;
  
      for (int i = 0; i &lt; n; ++i) {
          size_t size = sizes[i];
          size_t rep = rep_all / size;
  
          double t1 = for_loop                (rep, size);
          double t2 = eigen_3_size            (rep, size);
          double t3 = eigen_size_3            (rep, size);
          double t4 = eigen_vector3_vector    (rep, size);
          double t5 = eigen_stlvector_vector3 (rep, size);
  
          using namespace std;
          cout &lt;&lt; setw(8)  &lt;&lt; size 
               &lt;&lt; setw(13) &lt;&lt; t1 &lt;&lt; setw(13) &lt;&lt; t2 &lt;&lt; setw(13) &lt;&lt; t3
               &lt;&lt; setw(14) &lt;&lt; t4 &lt;&lt; setw(15) &lt;&lt; t5 &lt;&lt; endl;
      }
  }
  </code></pre>
  
  <p>The program was compiled by <code>gcc 4.6</code> with options <code>-march=native -O2 -msse2 -mfpmath=sse</code>. On my <code>Athlon 64 X2 4600+</code> I got a nice table:</p>
  
  <pre><code>  size    |  for loop  |   Matrix   |   Matrix   | Vector3 of  | STL vector of 
            |            |  3 x size  |  size x 3  | Vector_size |  TinyVectors  
  
      16         2.23          3.1         3.29          1.95           3.34
      32         2.12         2.72         3.51          2.25           2.95
      64         2.15         2.52         3.27          2.03           2.74
     128         2.22         2.43         3.14          1.92           2.66
     256         2.19         2.38         3.34          2.15           2.61
     512         2.17         2.36         3.54          2.28           2.59
    1024         2.16         2.35         3.52          2.28           2.58
    2048         2.16         2.36         3.43          2.42           2.59
    4096        11.57         5.35        20.29         13.88           5.23
    8192        11.55         5.31        16.17         13.79           5.24
  </code></pre>
  
  <p>The table shows that good representations of an array of 3d vectors are <code>Matrix</code> (components of 3d vectors should stored together) and <code>std::vector</code> of fixed size <code>Vector3d</code> objects. This confirms Jakob's answer. For big vectors <code>Eigen</code> indeed shows  good results.</p>