我决定在我的项目中使用Eigen库。但是从文档中不清楚如何最有效地指定一个 3d 向量数组。
正如我所建议的，第一种方法是
Eigen::Matrix<Eigen::Vector3d, Eigen::Dynamic, 1> array_of_v3d(size);
但是在那种情况下，我应该如何获得另一个数组，其中元素等于 的元素array_of_v3d和其他实例的标量积Vector3d？换句话说，我可以使用Eigen's 函数重写以下循环吗：
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- 部分展开的 for 循环Eigen::Matrix(rows x cols = 3 x size)。在这种情况下，3d 向量的值一起存储在内存中，因为默认情况下Eigen以列为主存储数据。或者我可以Eigen::RowMajor在下一个案例中设置和其他所有内容。Eigen::Matrix(rows x cols = size x 3)。- 这里 3d 向量的每个分量都存储在一个单独的
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.6with options编译-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确实显示出很好的结果。