我决定使用 Eigen我项目中的库。 但是文档中并不清楚如何最有效地指定 3d vector 数组。

按照我的建议，第一种方式是

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 部分展开的 for 循环样式
2. Eigen::Matrix(行 x 列 = 3 x 大小)。在这种情况下，3d vector 的值一起存储在内存中，因为默认情况下 Eigen 将数据存储在列优先中。或者我可以设置 Eigen::RowMajor 和下一个案例中的其他所有内容。
3. Eigen::Matrix(行 x 列 = 大小 x 3)。
4. 此处 3d vector 的每个分量都存储在单独的 VectorXd 中。所以在 Vector3d 上有三个 VectorXd 对象。
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+ 上，我得到了一张漂亮的 table :

  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 vector 数组的良好表示是 Matrix(3d vector 的组件应存储在一起)和固定大小的 std::vector Vector3d 对象。这证实了雅各布的回答。对于大 vector Eigen 确实显示了很好的结果。

最佳答案

如果您只想保存一个 Vector3d 对象数组，通常使用 std::vector 就可以了，尽管您需要注意 alignment问题。

您描述的第二种使用 size x 3 矩阵的方法也非常可行，而且通常速度更快。除了性能问题之外，我不确定您是否在问这个问题。

至于性能，我假设您确实在编译器中使用了完全优化，因为 Eigen 在关闭优化时表现不佳。在任何情况下，您都可以通过使用 aligned types 获得一些性能。可以使用优化的 SIMD 指令进行处理。如果您使用例如，Eigen 应该自动执行此操作Vector4d 代替。

关于c++ - 实现 3d vector 数组的最佳方法是什么？，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/12998410/