使用数组时的抽象与性能

Question

这是一个问题，我担心在处理数组时在更好的性能和更清晰的代码（更好的抽象）之间进行选择。我试图将其提炼成一个玩具示例。

C++ 特别擅长在不影响性能的情况下允许抽象。问题是这在类似于下面的示例中是否可行。

考虑一个使用连续行优先存储的普通任意大小的矩阵类：

#include <cmath>
#include <cassert>

class Matrix {
    int nrow, ncol;
    double *data;
public:
    Matrix(int nrow, int ncol) : nrow(nrow), ncol(ncol), data(new double[nrow*ncol]) { }
    ~Matrix() { delete [] data; }

    int rows() const { return nrow; }
    int cols() const { return ncol; }

    double & operator [] (int i) { return data[i]; }

    double & operator () (int i, int j) { return data[i*ncol + j]; }
};

它有一个二维索引operator ()，使其易于使用。它还有operator [] 用于连续访问，但更好抽象的矩阵可能没有这个。

让我们实现一个函数，它采用一个 n×2 矩阵，本质上是一个二维向量列表，并就地对每个向量进行归一化。

明确的方法：

inline double veclen(double x, double y) {
    return std::sqrt(x*x + y*y);
}

void normalize(Matrix &mat) {
    assert(mat.cols() == 2); // some kind of check for correct input
    for (int i=0; i < mat.rows(); ++i) {
        double norm = veclen(mat(i,0), mat(i,1));
        mat(i,0) /= norm;
        mat(i,1) /= norm;
    }
}

快速但不太清晰的方法：

void normalize2(Matrix &mat) {
    assert(mat.cols() == 2);
    for (int i=0; i < mat.rows(); ++i) {
        double norm = veclen(mat[2*i], mat[2*i+1]);
        mat[2*i] /= norm;
        mat[2*i+1] /= norm;
    }
}

第二个版本 (normalize2) 有可能更快，因为它的编写方式很明显，循环的第二次迭代不会访问在第一次迭代中计算的数据。因此，它可以更好地利用 SIMD 指令。 Looking at godbolt, this seems to be what happens（除非我误读了程序集）。

在第一个版本 (normalize) 中，编译器无法知道输入矩阵不是 n×1，这会导致数组访问重叠。

问题：是否有可能以某种方式告诉编译器输入矩阵在 normalize() 中实际上是 n×2 以使其优化到与 @987654330 中相同的水平@?

---

解决 cmets：

• John Zwinck：我去做了基准测试。 normalize2() 的速度要快得多（2.4 秒对 1.3 秒），但只有在我删除了assert 宏或者我定义了NDEBUG 时。这是-DNDEBUG 的一个相当违反直觉的效果，不是吗？它会降低性能而不是提高性能。

• Max：证据既是我链接到的神螺栓输出，也是上述基准。我也对不能内联这两个函数的情况感兴趣（例如，因为它们位于单独的翻译单元中）。

• Jarod42 和 bolov：这就是我一直在寻找的答案。由第一点中提到的基准确认。尽管如此，如果有人实现自己的assert（这正是我在我的应用程序中所做的），了解这一点很重要。

Answer 1

我相信模板可以让您同时获得性能和可读性。

通过模板化矩阵的大小（就像流行的数学库一样），您可以在编译时让编译器知道很多信息。

我修改了你的小课堂：

template<int R, int C>
class Matrix {
    double data[R * C] = {0.0};
public:
    Matrix() = default;

    int rows() const { return R; }
    int cols() const { return C; }
    int size() const { return R*C; }

    double & operator [] (int i) { return data[i]; }

    double & operator () (int row, int col) { return data[row*C + col]; }
};

inline double veclen(double x, double y) {
    return std::sqrt(x*x + y*y);
}

template<int R>
void normalize(Matrix<R, 2> &mat) {
    for (int i = 0; i < R; ++i) {
        double norm = veclen(mat(i, 0), mat(i, 1));
        mat(i, 0) /= norm;
        mat(i, 1) /= norm;
    }
}

template<int R>
void normalize2(Matrix<R, 2> &mat) {
    for (int i = 0; i < R; ++i) {
        double norm = veclen(mat[2 * i], mat[2 * i + 1]);
        mat[2 * i] /= norm;
        mat[2 * i + 1] /= norm;
    }
}

我也更喜欢将数据作为普通成员（=不带指针），因此您可以在矩阵构造期间选择内存所在的位置（堆栈或堆）。

额外的好处是您现在可以确定在编译时规范化函数只接受 n×2 矩阵。

我没有在编译器资源管理器上测试我的代码，因为老实说我无法破译 asm。所以，是的，我声称我的版本更好但不确定；）

最后一句话：不要滚动你自己的矩阵，使用一个库，比如 glm 或 eigen。

最后一句话²：如果您不知道该选择什么，请选择可读性。

Answer 2

我可以接受的答案基本上是由 cmets 中的 @bolov 和 @Jared42 给出的。由于他们没有发布，我会自己发布。

要让编译器知道矩阵的大小是 n×2，只要在函数的开头添加代码，当矩阵大小不正确时，其余代码就无法访问就足够了。

例如添加

if (mat.cols() != 2)
    throw std::runtime_error("Input array is not of expected shape.");

normalize() 的开头允许它与normalize2() 一样快地运行（在我使用 clang 5.0 的基准测试中是 1.3 而不是 2.4 秒）。

我们也可以添加assert(mat.cols() == 2)，但这会导致反直觉的效果，即在编译期间定义-DNDEBUG 会使函数变得相当慢（因为它删除了断言）。