Ole*_*yan 9 c++ linear-algebra eigen
我决定在我的项目中使用Eigen库.但是从文档中不清楚如何最有效地指定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种情况:
C-style for loop,部分展开Eigen::Matrix(rows x cols = 3 x size).在这种情况下,3d向量的值一起存储在内存中,因为默认情况下Eigen将数据存储在column-major中.或者我可以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
该表显示了3d矢量数组的良好表示(3d矢量的Matrix组件应该存储在一起)和std::vector固定大小的Vector3d对象.这证实了雅各布的答案.对于大载体Eigen确实显示出良好的结果.
如果你只想保存一个Vector3d对象数组，通常使用std::vector就可以了，尽管你需要注意对齐问题。
您描述的第二种方法使用size x 3矩阵也非常可行，而且通常是更快的方法。除了性能问题之外，我不确定您是否在问这个问题。
至于性能，我假设您确实在编译器中使用了完全优化，因为当优化关闭时，Eigen 的性能不佳。在任何情况下，您都可以通过使用可以使用优化的 SIMD 指令进行处理的对齐类型来获得一些性能。如果您使用egVector4d代替，Eigen 应该自动执行此操作。
| 归档时间： |
|
| 查看次数： |
4173 次 |
| 最近记录： |