张量积算法优化

Question

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等）。
• 奖金：如果你可以计算结果刚刚给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。基本上，前两个代码是运行的，执行时间最短的代码获胜。

Answer 1

一个好的编译器会autovectorize这 我想没有我的编译器是好的：

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]; 
     } 
    } 
    } 
} 

Answer 2

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

Answer 3

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

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

Answer 4

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进/ 64出。这个概念可以应用于任何数量的参数。你只需要实例化自己。

Answer 5

#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之外）。它也需要任何数量的参数。由于调用递归表单很尴尬，我做了一个包装器。它为每个参数调用一个函数，并且调用一次动态内存（在包装器中）。