从参数列表生成C++函数

您试图解决的关键问题是，在给定问题域的情况下，计算中会被多次调用的"leaf"函数通常也不会起作用。希望冗余工作(即将一个值与编译时已知为零的数组元素相乘)可以作为编译时步骤的一部分进行折叠。

C++有处理这一问题的语言设施，即模板元编程。C++模板非常强大(即图灵完备)，允许基于编译时间常数进行递归计算。

下面是一个如何使用模板和模板专业化实现示例的示例(您也可以在这里找到我创建的一个可运行的示例http://ideone.com/BDtBt7)。代码背后的基本思想是生成一个带有静态函数的类型，该函数返回值的输入向量和编译时常数数组的点积。静态函数递归地调用自身的实例，在元素的输入/常量数组中移动时传递较低的索引值。它还以正在计算的编译时常数数组p中的值是否为零为模板。如果是，我们可以跳过计算，然后转到递归中的下一个值。最后，有一个基本情况，一旦我们到达数组中的第一个元素，就会停止递归。

#include <array>
#include <iostream>
#include <vector>
constexpr std::array<double, 5> p = { 1.0, 0.0, 3.0, 5.0, 0.0 };
template<size_t index, bool isZero>
struct DotProductCalculator
{
static double Calculate(const std::vector<double>& xArg)
{
return (xArg[index] * p[index]) 
+ DotProductCalculator<index - 1, p[index - 1] == 0.0>::Calculate(xArg);
}
};
template<>
struct DotProductCalculator<0, true>
{
static double Calculate(const std::vector<double>& xArg)
{
return 0.0;
}
};
template<>
struct DotProductCalculator<0, false>
{
static double Calculate(const std::vector<double>& xArg)
{
return xArg[0] * p[0];
}
};
template<size_t index>
struct DotProductCalculator<index, true>
{
static double Calculate(const std::vector<double>& xArg)
{
return 0.0 + DotProductCalculator<index - 1, p[index - 1] == 0.0>::Calculate(xArg);
}
};
template<typename ArrayType>
double f_p_driver(const std::vector<double>& xArg, const ArrayType& pAsArgument) 
{
return DotProductCalculator<std::tuple_size<ArrayType>::value - 1, 
p[std::tuple_size<ArrayType>::value -1] == 0.0>::Calculate(xArg); 
}
int main() 
{
std::vector<double> x = { 1.0, 2.0, 3.0, 4.0, 5.0 };
double result = f_p_driver(x, p);
std::cout << "Result: " << result;
return 0;
}

你在评论中说，p实际上是矩阵的一行或一列，并且矩阵是稀疏的。我不熟悉你正在解决的具体物理问题，但通常情况下，稀疏矩阵有某种固定的对角线"带状"结构，例如：

| a1 b1  0  0  0  0  0 d1 |
| c1 a2 b2  0  0  0  0  0 |
| 0  c2 a3 b3  0  0  0  0 |
| 0  0  c3 a4 b4  0  0  0 |
| 0  0  0  c4 a5 b5  0  0 |
| 0  0  0  0  c5 a6 b6  0 |
| 0  0  0  0  0  c6 a7 b7 |
| e1 0  0  0  0  0  c7 a8 |

存储此类矩阵的最有效方法往往是将对角线存储为数组/向量，因此：

A = [a1, a2, a3, a4, a5, a6, a7, a8]
B = [b1, b2, b3, b4, b5, b6, b7]
C = [c1, c2, c3, c4, c5, c6, c7]
D = [d1]
E = [e1]

将行矢量X = [x1, x2, x3, x4, x5, x6, x7, x8]乘以上述矩阵，就变成：

Y = X . M
Y[0] =               X[0] * A[0] + X[1] * C[0] + X[7] * E[0]
Y[1] = X[0] * B[0] + X[1] * A[1] + X[2] * C[1]

等等。

或者更一般地说：

Y[i] = X[i-7] * D[i] + X[i-1] * B[i] + X[i] * A[i] + X[i+1] * C[i] + X[i+7] * E[i]

其中超出范围的数组访问(< 0或>= 8应被视为求值为0。为了避免在任何地方都测试越界，实际上可以将每个对角线和向量本身存储在超大数组中，这些数组的前导元素和尾部元素都填充了零。

请注意，由于所有阵列访问都是线性的，因此这也将具有很高的缓存效率。

在给定约束的情况下，我将创建一个自定义函数对象，该对象存储矩阵p，并在其函数调用运算符中计算运算。我将实现该函数的两个版本：一个版本在构造时对矩阵进行预处理，以"知道"非零元素在哪里，另一个版本只执行所述操作，接受许多计算只产生0。引用的10%非零元素的数量听起来可能过于密集，不利于利用稀疏性来获得回报。

忽略p是一个矩阵，并将其用作向量，没有预处理的版本将是这样的：

class dotProduct {
std::vector<double> p;
public:
dotProduct(std::vector<double> const& p): p(p) {}
double operator()(std::vector<double> const& x) const {
return std::inner_product(p.begin(), p.end(), x.begin());
}
};
// ...
... RungeKutta(dotProduct(p), initial conditions IC, etc.);

当使用C++11时，可以使用lambda函数：

... RungeKutta([=](std::vector<double> const& x) {
return std::inner_product(p.begin(), p.end(), x.begin());
}, intitial conditions IC, etc.);

对于预处理版本，您将存储一个std::vector<std::pair<double, std::size_t>>，指示哪些索引实际上需要相乘：

class sparseDotProduct {
typedef std::vector<std::pair<double, std::size_t>> Vector;
Vector p;
public:
sparsedotProduct(std::vector<double> const& op) {
for (std::size_t i(0), s(op.size()); i != s; ++i) {
if (op[i]) {
p.push_back(std::make_pair(op[i], i));
}
}
}
double operator()(std::vector<double> const& x) {
double result(0);
for (Vector::const_iterator it(p.begin()), end(p.end()); it != end; ++it) {
result += it->first * x[it->second];
}
return result;
}
};

这个函数对象的使用是一样的，尽管如果p不变，保留这个对象可能是合理的。

我个人认为，如果有10%的非零值，非稀疏版本实际上会优于稀疏版本。然而，有了这两个版本，衡量不同方法的性能应该相对简单。我不认为自定义创建的代码会更好，尽管它可以改进计算。如果是这样的话，使用元编程技术来创建代码可能会，但我怀疑这太实用了。