我试图使用OpenMP编写一个可重入算法，但成效有限。在github上有一个C版本的FORTRAN LSODA例程。它用于求解一阶常微分方程。它经历了几个版本，但最近的一个版本是Simon Frost，可以在这里找到：
https://github.com/sdwfrost/liblsoda
该库附带了一个简单的测试示例，我使用OpenMP更新为线程。关于并行实现，我有几个悬而未决的问题。我最“成功”的尝试是在下面（解决的经过时间没有减少，解决方案有时会有所不同）：
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include "common.h"
#include "lsoda.h"
#include "omp.h"
#define BILLION 1000000000L
#define NUM_THREADS 2 // Number of threads for OpenMP
struct timespec start, end;
long long int timeValues[2];
long long unsigned int deltaTimeArray;
double msTime = 0.;
int fex(double t, double *y, double *ydot, void *data)
{
ydot[0] = 1.0E4 * y[1] * y[2] - .04E0 * y[0];
ydot[2] = 3.0E7 * y[1] * y[1];
ydot[1] = -1.0 * (ydot[0] + ydot[2]);
return(0);
}
int test(void)
{
# ifdef _OPENMP
printf("Compiled by an OpenMP-compliant implementation.\n");
# endif
// Begin timing the algorithm
clock_gettime(CLOCK_MONOTONIC, &start);
timeValues[0] = (BILLION*start.tv_sec) + start.tv_nsec;
omp_set_dynamic(0);
omp_set_num_threads(NUM_THREADS);
int nThreads = 0;
#pragma omp parallel
{
double atol[3], rtol[3], t, tout, y[3];
int neq = 3;
int iout;
y[0] = 1.0E0;
y[1] = 0.0E0;
y[2] = 0.0E0;
t = 0.0E0;
tout = 0.4E0;
struct lsoda_opt_t opt = {0};
opt.ixpr = 0;
opt.rtol = rtol;
opt.atol = atol;
opt.itask = 1;
rtol[0] = rtol[2] = 1.0E-4;
rtol[1] = 1.0E-4;
atol[0] = 1.0E-6;
atol[1] = 1.0E-10;
atol[2] = 1.0E-6;
struct lsoda_context_t ctx = {
.function = fex,
.neq = neq,
.data = NULL,
.state = 1,
};
lsoda_prepare(&ctx, &opt);
#pragma omp master
nThreads = omp_get_num_threads();
#pragma omp for
for (iout = 1; iout <= 12; iout++)
{
lsoda(&ctx, y, &t, tout);
printf(" at t= %12.4e y= %14.6e %14.6e %14.6e\n", t, y[0], y[1], y[2]);
if (ctx.state <= 0)
{
printf("error istate = %d\n", ctx.state);
exit(0);
}
/*
if (iout == 1) tout = 4.0E-1 * 10.0E0;
if (iout == 2) tout = 4.0E0 * 10.0E0;
if (iout == 3) tout = 4.0E1 * 10.0E0;
if (iout == 4) tout = 4.0E2 * 10.0E0;
if (iout == 5) tout = 4.0E3 * 10.0E0;
if (iout == 6) tout = 4.0E4 * 10.0E0;
if (iout == 7) tout = 4.0E5 * 10.0E0;
if (iout == 8) tout = 4.0E6 * 10.0E0;
if (iout == 9) tout = 4.0E7 * 10.0E0;
if (iout == 10) tout = 4.0E8 * 10.0E0;
if (iout == 11) tout = 4.0E9 * 10.0E0;
if (iout == 12) tout = 4.0E10 * 10.0E0;
*/
tout = tout * 10.0E0;
}
lsoda_free(&ctx);
}
if (nThreads == NUM_THREADS)
{
printf("The expected number of threads, %d, were used.\n", NUM_THREADS);
}
else
{
printf("Expected %d OpenMP threads, but %d were used.\n", NUM_THREADS, nThreads);
}
// End timing the algorithm
clock_gettime(CLOCK_MONOTONIC, &end);
timeValues[1] = (BILLION*end.tv_sec) + end.tv_nsec;
deltaTimeArray = timeValues[1] - timeValues[0];
timeValues[0] = timeValues[1];
msTime = deltaTimeArray * pow(10, -6);
const char *elapsed = "elapsed";
printf("%13s\n%13.06lf\n", elapsed, msTime);
return(0);
}
int main(void)
{
for(int i = 0; i < 1; i++)
{
test();
}
return(0);
}
/*
The correct answer (up to certain precision):
at t= 4.0000e-01 y= 9.851712e-01 3.386380e-05 1.479493e-02
at t= 4.0000e+00 y= 9.055333e-01 2.240655e-05 9.444430e-02
at t= 4.0000e+01 y= 7.158403e-01 9.186334e-06 2.841505e-01
at t= 4.0000e+02 y= 4.505250e-01 3.222964e-06 5.494717e-01
at t= 4.0000e+03 y= 1.831976e-01 8.941773e-07 8.168015e-01
at t= 4.0000e+04 y= 3.898729e-02 1.621940e-07 9.610125e-01
at t= 4.0000e+05 y= 4.936362e-03 1.984221e-08 9.950636e-01
at t= 4.0000e+06 y= 5.161833e-04 2.065787e-09 9.994838e-01
at t= 4.0000e+07 y= 5.179804e-05 2.072027e-10 9.999482e-01
at t= 4.0000e+08 y= 5.283675e-06 2.113481e-11 9.999947e-01
at t= 4.0000e+09 y= 4.658667e-07 1.863468e-12 9.999995e-01
at t= 4.0000e+10 y= 1.431100e-08 5.724404e-14 1.000000e+00
*/
似乎有几个“omp.h”的变种，我在这里可以找到我使用的版本：
https://sites.uclouvain.be/SystInfo/usr/include/omp.h.html
我不确定这个算法是如何被线程化的。第一个突出的问题是调用lsoda（）例程的for循环是否应该是线程化的，因为tout具有循环依赖性。如果它是一个增量（就像在我的实际实现中那样），它可以简单地使用逗号运算符在for循环中添加为第二个增量语句（编辑：显然这不能用OpenMP完成 - for循环必须是规范的形式，不能有两个增量表达式）。这是注释掉的部分代码根据iout的值更新tout的原因。即使使用那些笨重的代码块，有时也会导致错误的解决方案。
提供的测试示例是使用嵌入在for循环中的test（）函数编写的，并在main中调用（在一次迭代中）。我不确定这背后的逻辑，但在我的情况下，我正在做类似的事情，在main中的for循环中多次调用test（）例程。对循环进行线程化的问题再次是循环依赖性，其中每个连续的解决方案（y []值）依赖于先前的解决方案。我知道这个算法能够被线程化，但我似乎无法成功实现它。任何指针将不胜感激。
似乎有几个“omp.h”的变种，我使用的版本可以在这里找到：[...]
不可以。您无权选择一些随机的OMP标头。您应该使用您正在使用的OMP实现提供的那个，否则您将面临未定义的行为风险。
我不确定这个算法是如何被线程化的。第一个突出的问题是调用lsoda（）例程的for循环是否应该是线程化的，因为tout具有循环依赖性。
你很关心依赖性。 OpenMP没有任何魔力 - 程序员可以处理依赖关系。但是，在这种情况下，您可以选择打破涉及tout的依赖关系，其中包括：
tout值的数组并传递该数组的元素：
double touts[12] = { 0.4 };
for (int i = 1; i < 12; i++) touts[i] = 10 * touts[i - 1];
// ...
lsoda(&ctx, y, &t, touts[i - 1]);
tout：
lsoda(&ctx, y, &t, tout * pow(10, i - 1));
// ... and avoid modifying tout later in the loop ...
但请注意其他依赖项。特别是，如果lsoda()既读取又修改了其他参数所指向的数据，那么这可能会引入更难的 - 也许是不可能的 - 来处理。
你似乎在说这确实是这样的：
对循环进行线程化的问题再次是循环依赖性，其中每个连续的解决方案（y []值）依赖于先前的解决方案。
如果lsoda()计算并存储y元素的新值，并且这些值依赖于原始值的难以预测的方式，则那就是游戏结束。你不会打破这种依赖性，如果不打破它，你最好的结果是你得到正确的输出而没有任何加速。并不是说你可以放心地看到最好的情况。
我知道这个算法能够被线程化，但我似乎无法成功实现它。
并非所有算法都可并行化。您可以使用OpenMP和lsoda()同时解决几个单独的问题，这是合理的，可以想象lsoda()可以内部并行化（我没有评估这种可能性），但我认为没有理由认为您的特定测试用例可以有效并行化。