числовой — генерация функции C ++ из списка аргументов

С учетом этих ограничений я бы создал объект пользовательской функции, который хранит матрицу 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 вместо этого может использоваться лямбда-функция:

... 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% ненулевых значений. Тем не менее, с этими двумя версиями должно быть относительно просто измерить производительность различных подходов. Я бы не ожидал, что созданный на заказ код будет значительно лучше, хотя это может улучшить вычисления. Если так, то это может работать, чтобы использовать методы мета-программирования для создания кода, но я сомневаюсь, что это было бы слишком практично.

В комментариях вы говорите, что 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. Чтобы избежать повсеместного тестирования на выход за пределы, вы можете хранить каждую диагональ и сам вектор в массивах больших размеров, в которых начальные и конечные элементы заполнены нулями.

Обратите внимание, что это также будет очень эффективно кешировать, поскольку все обращения к массивам являются линейными