张量积算法优化答案

【问题标题】：Tensor Product Algorithm Optimization张量积算法优化
【发布时间】：2011-09-20 19:00:15
【问题描述】：

double data[12] = {1, z, z^2, z^3, 1, y, y^2, y^3, 1, x, x^2, x^3};
double result[64] = {1, z, z^2, z^3, y, zy, (z^2)y, (z^3)y, y^2, z(y^2), (z^2)(y^2), (z^3)(y^2), y^3, z(y^3), (z^2)(y^3), (z^3)(y^3), x, zx, (z^2)x, (z^3)x, yx, zyx, (z^2)yx, (z^3)yx, (y^2)x, z(y^2)x, (z^2)(y^2)x, (z^3)(y^2)x, (y^3)x, z(y^3)x, (z^2)(y^3)x, (z^3)(y^3)x, x^2, z(x^2), (z^2)(x^2), (z^3)(x^2), y(x^2), zy(x^2), (z^2)y(x^2), (z^3)y(x^2), (y^2)(x^2), z(y^2)(x^2), (z^2)(y^2)(x^2), (z^3)(y^2)(x^2), (y^3)(x^2), z(y^3)(x^2), (z^2)(y^3)(x^2), (z^3)(y^3)(x^2), x^3, z(x^3), (z^2)(x^3), (z^3)(x^3), y(x^3), zy(x^3), (z^2)y(x^3), (z^3)y(x^3), (y^2)(x^3), z(y^2)(x^3), (z^2)(y^2)(x^3), (z^3)(y^2)(x^3), (y^3)(x^3), z(y^3)(x^3), (z^2)(y^3)(x^3), (z^3)(y^3)(x^3)};

在给定数据的情况下，产生结果的最快（最少执行）是什么？假设数据的大小是可变的，但始终是 4 倍（例如，4、8、12 等）。
没有提升。我试图保持我的依赖小。 STL 算法没问题。
提示：结果数组大小应始终为 4^（倍数大小）（例如，4、16、64 等）。
奖励：如果您可以计算 result 刚刚给定 x, y, z

其他示例：

double data[4] = {1, z, z^2, z^3};
double result[4] = {1, z, z^2, z^3};

double data[8] = {1, z, z^2, z^3, 1, y, y^2, y^3};
double result[16] = { ... };

我在运行此基准测试后选择了接受的答案代码：https://gist.github.com/1232406。基本上，前两个代码被运行，执行时间最短的一个获胜。

【问题讨论】：

你有算法太慢，还是需要算法？
@Als 不是作业；它的工作。
@Seth 我已经有一个解决方案，但是速度很慢，并且引入了外部依赖项（特别是 Eigen eigen.tuxfamily.org）来执行 Kronecker 产品。
如果您使用的是 Eigen，而且速度很慢，那么无论是您做错了还是做得更快，都意味着英特尔数学库的魔力。
我很可能使用 Eigen 做错了；我会同意的。基本上，我无法让 Eigen 在 Map<VectorXd> 上执行 Kronecker 乘积而不创建 Map 的副本，从而执行内存分配然后复制，从而减慢执行速度。

标签： c++ algorithm optimization

【解决方案1】：

void Tensor(std::vector<double>& result, double x, double y, double z) {
    result.resize(64); //almost noop if already right size
    double tz = z*z;
    double ty = y*y;
    double tx = x*x;
    std::array<double, 12> data = {0, 0, tz, tz*z, 1, y, ty, ty*y, 1, x, tx, tx*x};
    register std::vector<double>::iterator iter = result.begin();
    register int yi;
    register double xy;
    for(register int xi=0; xi<4; ++xi) {
        for(yi=0; yi<4; ++yi) {
            xy = data[4+yi]*data[8+xi];
            *iter = xy; //a smart compiler can do these four in parallell
            *(++iter) = z*xy;
            *(++iter) = data[2]*xy;
            *(++iter) = data[3]*xy;
            ++iter; //workaround for speed!
        }
    }        
}

这里的某个地方可能至少存在一个错误，但它应该很快，没有依赖关系（std::vector/std::array 之外），只需要 x,y,z。我避免了递归，所以它只适用于 3 in/64 out。不过，这个概念可以应用于任意数量的参数。你只需要实例化自己。

【讨论】：

在循环前使用register int xi, yi; register double xy;，在循环体中使用xy = data[4+yi]*data[8+xi]; *iter = xy; *(++iter) = z*xy; *(++iter) = data[2]*xy; *(++iter) = data[3]*xy; ++iter;。 iter++ 通常实现为 iterator temp(*this); ++(*this); return temp;。
我通常会忽略 register 关键字，因为我的编译器会这样做。我添加了它以防某些编译器仍然关心它。我知道++iter 与iter++，但由于它处于循环中，因此优化在这里不起作用。另外，由于vector::iterator 基本上是一个指针，所以它并不是真正的减速。
没有完全理解问题，也没有正确阅读答案，但为努力+1！
我已经运行了代码，它正确地完成了我上面展示的示例。所以我对答案投了赞成票。我没有接受它的原因是它还不适用于通用案例。我现在正在努力，当我扩展你的方法时，我会接受答案。
哪个通用案例？更多变数？不同类型的变量？ 0-N 的幂？

【解决方案2】：

~~一个好的编译器会自动向量化这个~~我猜我的编译器都不是好的：

void tensor(const double *restrict data,
            int dimensions,
            double *restrict result) {
  result[0] = 1.0;
  for (int i = 0; i < dimensions; i++) {
    for (int j = (1 << (i * 2)) - 1; j > -1; j--) {
      double alpha = result[j];
      {
        double *restrict dst = &result[j * 4];
        const double *restrict src = &data[(dimensions - 1 - i) * 4];
        for (int k = 0; k < 4; k++) dst[k] = alpha * src[k];
      }
    }
  }
}

【讨论】：

这有点用。我正在使用 VS 2010 并且 restrict 关键字无效；它使用__restrict。此外，由于某种原因，在使用编译器优化后运行代码时，它不起作用；虽然如果我在没有优化的情况下运行它，它确实可以工作。
去掉 VS 中的__restrict 关键字修复了编译器优化错误。
@Ryan：显然，对这些答案中的几个进行计时。此方法多次写入数组的某些元素，但代码太小，它可能在缓存中执行得更好，因此它可能比其他方法更快。
@Mooing Duck 这个会多次写入数组的一些元素 这个不会假设每个四边形的第一个元素是 1。这很容易在这种假设下避免无关的写入。
@以前叫 d 的家伙：实际上，是的。你说的对。我现在明白了。但是，考虑到每个四边形的第一个元素是 1，它可以更快。ideone.com/wRCaw

【解决方案3】：

你应该使用动态算法。也就是说，您可以使用以前的结果。例如，您保留 y^2 结果并在计算 (y^2)z 时使用它，而不是再次计算它。

【讨论】：

【解决方案4】：

#include <vector> 
#include <cstddef>
#include <cmath>

void Tensor(std::vector<double>& result, const std::vector<double>& variables, size_t index)    
{
    double p1 = variables[index];
    double p2 = p1*p1;
    double p3 = p1*p2;
    if (index == variables.size() - 1) {
        result.push_back(1);
        result.push_back(p1);
        result.push_back(p2);
        result.push_back(p3);
    } else {
        Tensor(result, variables, index+1);
        ptrdiff_t size = result.size();
        for(int j=0; j<size; ++j)
            result.push_back(result[j]*p1);
        for(int j=0; j<size; ++j)
            result.push_back(result[j]*p2);
        for(int j=0; j<size; ++j)
            result.push_back(result[j]*p3);
    }
}

std::vector<double> Tensor(const std::vector<double>& params) {
    std::vector<double> result;
    double rsize = (1<<(2*params.size());
    result.reserve(rsize);
    Tensor(result, params);
    return result;
}

int main() {
    std::vector<double> params;
    params.push_back(3.1415926535);
    params.push_back(2.7182818284);
    params.push_back(42);
    params.push_back(65536);
    std::vector<double> result = Tensor(params);
}

我验证了这个可以编译和运行 (http://ideone.com/IU1eQ)。它运行速度很快，没有依赖关系（在 std::vector 之外）。它还接受任意数量的参数。由于调用递归形式很尴尬，我做了一个包装器。它对每个参数进行一次函数调用，并对动态内存进行一次调用（在包装器中）。

【讨论】：

【解决方案5】：

您应该寻找Pascal's pyramid 以获得快速解决方案。 Useful link 1、useful link 2、useful link 3 和 useful link 4。

还有一件事：我认为它是有限元求解器的基础。通常编写自己的 BLAS 求解器不是一个好主意。不要重新发明轮子！我认为你应该使用像 intel MKL 或 Cuda base BLAS 这样的 BLAS 求解器。

【讨论】：