Я решил использовать в своем проекте библиотеку Eigen. Но из документации не ясно, как наиболее эффективно указывать массив трехмерных векторов.
Как я предлагаю, первый способ
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 случаев:
C-стиль для частично развернутого циклаEigen::Matrix(rows x cols = 3 x size). В этом случае значения трехмерных векторов хранятся в памяти вместе, потому что по умолчаниюEigenхранит данные в основном столбце. Или я мог бы установитьEigen::RowMajorи все остальное, как в следующем случае.Eigen::Matrix(rows x cols = size x 3).- Здесь каждый компонент трехмерного вектора хранится в отдельном
VectorXd. Таким образом, есть три объекта VectorXd, которые объединены вVector3d. - Контейнер
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
Таблица показывает, что хорошими представлениями массива трехмерных векторов являются Matrix (компоненты трехмерных векторов должны храниться вместе) и std::vector объектов фиксированного размера Vector3d. Это подтверждает ответ Якоба. Для больших векторов Eigen действительно показывает хорошие результаты.
