Каков наилучший способ реализовать массив 3d векторов?

Я решил использовать собственный библиотека в моем проекте.
Но из документации не ясно, как наиболее эффективно указать массив 3d векторов.

Как я предполагаю, первый способ

Eigen::Matrix<Eigen::Vector3d, Eigen::Dynamic, 1> array_of_v3d(size);

Но в таком случае, как я должен получить другой массив, элементы которого равны скалярным произведениям элементов array_of_v3d и какой-то другой случай Vector3d?
Другими словами, могу ли я переписать следующий цикл, используя Eigenфункции:

Eigen::Vector3d v3d(0.5, 0.5, 0.5);
Eigen::VectorXd other_array(size);
for (size_t i = 0; i < size; ++i)
other_array(i) = array_of_v3d(i).dot(v3d);

Второй способ — использовать матрицу, размер которой (3 x size) или же (size x 3),
Например, я могу объявить это так:

Eigen::Matrix<double, 3, Eigen::Dynamic> matrix;

Но я не понял из документации, как установить количество столбцов.
Кажется, что следующее работает, но я должен повторить количество строк 3 дважды:

Eigen::Matrix<double, 3, Eigen::Dynamic> matrix(3, size);

Тогда вышеуказанный цикл эквивалентен

other_array = v3d.transpose() * array_of_v3d;

Как показали мои эксперименты, это немного быстрее, чем

Eigen::Matrix<double, Eigen::Dynamic, 3> matrix(size, 3);
other_array = array_of_v3d * v3d;

К тому же:

Во всяком случае, мое использование Eigen кажется не очень оптимальным, так как одна и та же программа на простом C почти в 1,5 раза быстрее (на самом деле это не так, это зависит от size):

for (size_t i = 0; i < 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;
}

Я что-то пропустил? Есть ли какой-то другой способ в объеме Eigen библиотека, чтобы решить мою проблему?

Обновить:

Здесь я представляю результаты моих испытаний. Есть 5 случаев:

1. C-стиль для цикла, который частично развернут
2. Eigen::Matrix (rows x cols = 3 x size). В этом случае значения 3d векторов хранятся вместе в памяти, потому что по умолчанию Eigen хранит данные в столбце-мажоре. Или я мог бы установить Eigen::RowMajor а все остальное как в следующем случае.
3. Eigen::Matrix (rows x cols = size x 3).
4. Здесь каждый компонент трехмерного вектора хранится на отдельном VectorXd, Таким образом, есть три объекта VectorXd, которые объединены в Vector3d,
5. std::vector контейнер используется для хранения Vector3d объекты.

Это моя тестовая программа

#include <iostream>
#include <iomanip>
#include <vector>
#include <ctime>

#include <Eigen/Dense>

double for_loop(size_t rep, size_t size)
{
std::vector<double>  a(size), b(size), c(size);
double x = 1, y = 2, z = - 3;
std::vector<double>  s(size);
for(size_t i = 0; i < size; ++i) {
a[i] = i;
b[i] = i;
c[i] = i;
s[i] = 0;
}
double dtime = clock();
for(size_t j = 0; j < rep; j++)
for(size_t i = 0; i < 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 < size; ++i)
res += std::abs(s[i]);
assert(res == 0.);

return dtime;
}

double eigen_3_size(size_t rep, size_t size)
{
Eigen::Matrix<double, 3, Eigen::Dynamic> A(3, size);
Eigen::Matrix<double, 1, Eigen::Dynamic> S(size);
Eigen::Vector3d X(1, 2, -3);
for(size_t i = 0; i < 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 < 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<double, Eigen::Dynamic, 3> A(size, 3);
Eigen::Matrix<double, Eigen::Dynamic, 1> S(size);
Eigen::Vector3d X(1, 2, -3);
for(size_t i = 0; i < 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 < 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<Eigen::VectorXd, 3, 1> 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 < 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 < 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<                 Eigen::Vector3d,
Eigen::aligned_allocator<Eigen::Vector3d> > A(size);
std::vector<double> S(size);
Eigen::Vector3d X(1, 2, -3);
for(size_t i = 0; i < 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 < rep; j++)
for(size_t i = 0; i < size; i++)
S[i] += X.dot(A[i]);
dtime = (clock() - dtime) / CLOCKS_PER_SEC;

double res = 0;
for(size_t i = 0; i < size; i++)
res += std::abs(S[i]);
assert(res == 0.);

return dtime;
}

int main()
{
std::cout << "    size    |  for loop  |   Matrix   |   Matrix   | Vector3 of  | STL vector of \n"<< "            |            |  3 x size  |  size x 3  | Vector_size |  TinyVectors  \n"<< 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 < 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 << setw(8)  << size
<< setw(13) << t1 << setw(13) << t2 << setw(13) << t3
<< setw(14) << t4 << setw(15) << t5 << endl;
}
}

Программа была составлена gcc 4.6 с вариантами -march=native -O2 -msse2 -mfpmath=sse, На моем Athlon 64 X2 4600+ Я получил хороший стол:

  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

Таблица показывает, что хорошие представления массива 3d векторов Matrix (компоненты 3d векторов должны храниться вместе) и std::vector фиксированного размера Vector3d объекты. Это подтверждает ответ Якоба. Для больших векторов Eigen действительно показывает хорошие результаты.