OpenMP #pragma omp parallel for 速度较慢
OpenMP #pragma omp parallel for is slower
我正在尝试改进我的 C 源代码以并行执行。我有一个四核 CPU，所以我认为 4 个线程（一个用于 CPU）对 运行 我的程序来说真的比像顺序代码那样优化更快。
但它不起作用。我的没有 OpenMP 的代码需要 11 分钟才能执行，而并行代码需要超过 11 分钟。我报告了我的所有来源，但并行代码仅在 getBasin() 函数中。
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <omp.h>
#define BLD "\x1B[1m"
#define RED "\x1B[31m"
#define GRN "\x1B[32m"
#define RST "\x1B[0m"
struct point{
double x;
double v;
double t;
double E0;
double E;
double dE;
} typedef point;
struct parametri {
double gamma;
double epsilon;
double dt;
double t_star;
double t_basin;
double t_max;
double x0;
double v0;
int choice;
int alg[1];
double dx;
} typedef parametri;
// Prototipi delle funzioni
void Setup();
void getParams(parametri *pars);
void getError(point *xv, parametri *pars, int *i, int k, int *n_passi, double *xn1, double *vn1);
void getBasin(point *xv, parametri *pars, int *h, int *n_passi, double *xn1, double *vn1);
double f(double x, double v, parametri *pars);
void algEulero(point *xv, parametri *pars, double *xn1, double *vn1);
void algEuleroCromer(point *xv, parametri *pars, double *xn1, double *vn1);
void algPuntoDiMezzo(point *xv, parametri *pars, double *xn1, double *vn1);
void algVerlet(point *xv, parametri *pars, double *xn1, double *vn1);
void algRungeKutta2(point *xv, parametri *pars, double *xn1, double *vn1);
int main(void) {
// Inizializzo il display per l'utente. Maggiori informazioni vedere funzione Setup();
Setup();
// Dichiaro e recupero i parametri richiesti per il corretto funzionamento del programma;
parametri pars;
getParams(&pars);
// Dichiaro e ricavo i parametri essenziali, annesse variabili di supporto;
int i, n_passi = pars.t_max/pars.dt;
double dt0, xn1, vn1, t_star = 5.;
point xv;
// Imposto i parametri iniziali;
xv.x = pars.x0;
xv.v = pars.v0;
xv.E0 = 0.5*(xv.v)*(xv.v) - (xv.x)*(xv.x)/8. + (xv.x)*(xv.x)*(xv.x)*(xv.x)/4.;
xv.E = xv.E0;
xv.dE = 0;
xv.t = 0;
pars.t_star = 5;
pars.t_basin = 60;
dt0 = pars.dt;
// Formato dell'output;
printf ("t\tx(t)\tv(t)\tE(t)\tdE\n");
if ((pars.choice == 1) || (pars.choice == 3) || (pars.choice == 4)) { // L'utente ha deciso di affrontare il primo/terzo/quarto esercizio;
// Avverto l'utente che il processo sta per iniziare. Può risultare inutile su tempi brevi, ma efficace su tempi molto lunghi;
fprintf(stderr, "\nAvvio integrazione numerica... ");
if (pars.alg[0] == 1) { // L'utente ha selezionato l'algoritmo di Eulero;
for (i=0; i<n_passi; i++){
algEulero(&xv, &pars, &xn1, &vn1);
}
} else if (pars.alg[0] == 2) { // L'utente ha selezionato l'algoritmo di EuleroCromer;
for (i=0; i<n_passi; i++){
algEuleroCromer(&xv, &pars, &xn1, &vn1);
}
} else if (pars.alg[0] == 3) { // L'utente ha selezionato l'algoritmo di PuntoDiMezzo;
for (i=0; i<n_passi; i++){
algPuntoDiMezzo(&xv, &pars, &xn1, &vn1);
}
} else if (pars.alg[0] == 4) { // L'utente ha selezionato l'algoritmo di Verlet;
for (i=0; i<n_passi; i++){
algVerlet(&xv, &pars, &xn1, &vn1);
}
} else if (pars.alg[0] == 5) { // L'utente ha selezionato l'algoritmo di RungeKutta di ordine 2;
for (i=0; i<n_passi; i++) {
algRungeKutta2(&xv, &pars, &xn1, &vn1);
}
}
// Algoritmo terminato;
fprintf(stderr, "[%s%sDONE%s]\n\n", BLD, GRN, RST);
} else if (pars.choice == 2) { // L'utente ha deciso di affrontare il secondo esercizio;
// Seleziono il secondo algoritmo da confrontare
do {
fprintf(stderr, "> Selezionare il secondo algoritmo:\n");
fprintf(stderr, "\t>> [1] Eulero\n");
fprintf(stderr, "\t>> [2] EuleroCromer - RungeKutta (1 ordine)\n");
fprintf(stderr, "\t>> [3] PuntoDiMezzo\n");
fprintf(stderr, "\t>> [4] Verlet\n");
fprintf(stderr, "\t>> [5] RungeKutta (2 ordine)\n");
fprintf(stderr, "\t>> ");
scanf("%d", &pars.alg[1]);
} while (( pars.alg[1] <= 0 ));
// Avverto l'utente che il processo sta per iniziare. Può risultare inutile su tempi brevi, ma efficace su tempi molto lunghi;
fprintf(stderr, "\nAvvio calcolo errori... ");
// Eseguo lo studio degli errori d'integrazione mediante il primo e secondo algoritmo scelto, rispettivamente nei file error_1.dat, error_2.dat;
getError(&xv, &pars, &i, 0, &n_passi, &xn1, &vn1);
// Resetto le variabili e richiamo l'algoritmo per calcolare gli errori;
pars.dt = dt0;
n_passi = pars.t_max/pars.dt;
xv.x = pars.x0;
xv.v = pars.v0;
xv.E = xv.E0;
xv.dE = 0;
xv.t = 0;
getError(&xv, &pars, &i, 1, &n_passi, &xn1, &vn1);
// Processo terminato;
fprintf(stderr, "\n\nStato: [%s%sDONE%s]\n\n", BLD, GRN, RST);
} else if (pars.choice == 5) { // L'utente ha deciso di affrontare il quinto esercizio;
// Avverto l'utente che il processo sta per iniziare. Può risultare inutile su tempi brevi, ma efficace su tempi molto lunghi;
fprintf(stderr, "\nAvvio calcolo griglia... ");
getBasin(&xv, &pars, &i, &n_passi, &xn1, &vn1);
// Processo terminato;
fprintf(stderr, "[%s%sDONE%s]\n\n", BLD, GRN, RST);
} else { // L'utente non ha selezionato un esercizio valido;
fprintf(stderr, "\n[%s%sFAILED%s] Esercizio non disponibile.\n", BLD, RED, RST);
exit(EXIT_FAILURE);
}
return 0;
}
void Setup() {
fprintf(stderr, "\nAnalisi numerica di un'equazione differenziale\n");
fprintf(stderr, "==============================================\n\n");
}
void getParams(parametri *pars) {
do {
fprintf(stderr, "> Inserire un valore per gamma : ");
scanf("%lf", &pars->gamma);
} while (pars->gamma < 0);
do {
fprintf(stderr, "> Inserire un valore per epsilon : ");
scanf("%lf", &pars->epsilon);
} while (pars->epsilon < 0);
do {
fprintf(stderr, "> Inserire un valore per dt : ");
scanf("%lf", &pars->dt);
} while (pars->dt < 0);
do {
fprintf(stderr, "> Inserire un valore per tmax t.c. :\n");
fprintf(stderr, " >> (tmax > 5) per I eserc.\n");
fprintf(stderr, " >> (tmax > 60) per V eserc. : ");
scanf("%lf", &pars->t_max);
} while (pars->t_max < 0);
do {
fprintf(stderr, "> Inserire un valore per x(0) : ");
scanf("%lf", &pars->x0);
} while (pars->x0 < -1);
do {
fprintf(stderr, "> Inserire un valore per v(0) : ");
scanf("%lf", &pars->v0);
} while (pars->v0 < -1);
do {
fprintf(stderr, "> Selezionare l'esercizio richiesto :\n");
fprintf(stderr, "\t>> [1] Esercizio 1\n");
fprintf(stderr, "\t>> [2] Esercizio 2\n");
fprintf(stderr, "\t>> [3] Esercizio 3\n");
fprintf(stderr, "\t>> [4] Esercizio 4\n");
fprintf(stderr, "\t>> [5] Esercizio 5\n\n");
fprintf(stderr, "\t>> ");
scanf("%d", &pars->choice);
} while (( pars->choice <= 0 ));
do {
fprintf(stderr, "> Selezionare l'algoritmo voluto :\n");
fprintf(stderr, "\t>> [1] Eulero\n");
fprintf(stderr, "\t>> [2] EuleroCromer - RungeKutta (1 ordine)\n");
fprintf(stderr, "\t>> [3] PuntoDiMezzo\n");
fprintf(stderr, "\t>> [4] Verlet\n");
fprintf(stderr, "\t>> [5] RungeKutta (2 ordine)\n\n");
fprintf(stderr, "\t>> ");
scanf("%d", &pars->alg[0]);
} while (( pars->alg[0] <= 0 ));
}
void getError(point *xv, parametri *pars, int *i, int k, int *n_passi, double *xn1, double *vn1) {
void (*algF)(point *, parametri *, double *, double *);
int j, n = 0;
FILE *fp;
// Questo if controlla se il codice sta eseguendo lo studio degli errori per il primo o secondo algoritmo (pars->alg[k]);
if (k == 0) {
fp = fopen("error_1.dat", "w+");
} else {
fp = fopen("error_2.dat", "w+");
}
// Assegno il puntatore corretto a seconda dell'algoritmo scelto;
if (pars->alg[k] == 1) {
algF = algEulero;
} else if (pars->alg[k] == 2) {
algF = algEuleroCromer;
} else if (pars->alg[k] == 3) {
algF = algPuntoDiMezzo;
} else if (pars->alg[k] == 4) {
algF = algVerlet;
} else if (pars->alg[k] == 5) {
algF = algRungeKutta2;
} else {
fprintf(stderr, "\n[%s%sFAILED%s] E' stato selezionato un algoritmo non valido.\n", BLD, RED, RST);
exit(EXIT_FAILURE);
}
// Informo l'utente che il processo sta iniziando;
fprintf(stderr, "\n>> Avvio %d algoritmo... ", k+1);
// Formattazione dell'output del file contenente gli errori;
fprintf(fp, "dt\tE(t*)\tE(0)\tdE/E(0)\t# passi\ti\tt\n");
for (j=0; j<8; j++) {
for ((*i)=0; (*i)<(*n_passi); (*i)++){
(*algF)(xv, pars, xn1, vn1);
if (((pars->t_star - pars->dt/2. <= xv->t) && (xv->t >= pars->t_star + pars->dt/2.)) && (n == 0)) {
fprintf(fp, "%+.14lf\t%+.14lf\t%+.14lf\t%+.14lf\t%+.14d\t%+.14d\t%+.14lf\n", pars->dt, xv->E, xv->E0, (xv->E - xv->E0)/xv->E0, (*n_passi), (*i), xv->t);
n = 1;
}
}
// Resetto le variabili per rilanciare l'algoritmo con dt/2^j
n = 0;
xv->t = 0;
xv->x = pars->x0;
xv->v = pars->v0;
pars->dt = pars->dt/2.;
(*n_passi) = pars->t_max/pars->dt;
(*xn1) = 0;
(*vn1) = 0;
xv->E = xv->E0;
}
fclose(fp);
fprintf(stderr, "[%s%sDONE%s]", BLD, GRN, RST);
}
void getBasin(point *xv, parametri *pars, int *h, int *n_passi, double *xn1, double *vn1) {
// Dichiaro e setto i parametri che delimitano la griglia;
point xv1;
pars->dx = 0.01;
xv1.x = -1;
xv1.v = -1;
// Dichiaro le variabili necessarie per il bacino d'attrazione;
int i, j, i_max = 2./pars->dx, j_max = 2./pars->dx;
void (*algF)(point *, parametri *, double *, double *);
FILE *fp = fopen("basin.dat", "w+");
// Assegno il puntatore corretto a seconda dell'algoritmo scelto;
if (pars->alg[0] == 1) {
algF = algEulero;
} else if (pars->alg[0] == 2) {
algF = algEuleroCromer;
} else if (pars->alg[0] == 3) {
algF = algPuntoDiMezzo;
} else if (pars->alg[0] == 4) {
algF = algVerlet;
} else if (pars->alg[0] == 5) {
algF = algRungeKutta2;
} else {
fprintf(stderr, "\n[%s%sFAILED%s] E' stato selezionato un algoritmo non valido.\n", BLD, RED, RST);
exit(EXIT_FAILURE);
}
// Formattazione output file basin.dat;
fprintf(fp, "x(0)\tx'(0)\tx(t*)\tv(t*)\n");
omp_set_num_threads(4);
#pragma omp parallel for
// Eseguo il for della griglia sull'asse x';
for (j=0; j<=j_max; j++) {
// Eseguo il for della griglia sull'asse x;
for (i=0; i<=i_max; i++) {
fprintf(fp, "%lf\t%lf\t", xv1.x, xv1.v);
xv->x = xv1.x;
xv->v = xv1.v;
// Eseguo l'integrazione numerica
for ((*h)=0; (*h)<(*n_passi); (*h)++) {
(*algF)(xv, pars, xn1, vn1);
// Entro in t = t*, stampo v(t*) ed esco dal ciclo;
if ((pars->t_basin - pars->dt/2. <= xv->t) && (xv->t >= pars->t_basin + pars->dt/2.)) {
fprintf(fp, "%lf\t%lf\n", xv->x, xv->v);
break;
}
}
xv1.x += pars->dx;
xv->t = 0;
(*xn1) = 0;
(*vn1) = 0;
}
// Resetto la x e incremento la x';
xv1.x = -1;
xv1.v += pars->dx;
}
}
double f(double x, double v, parametri *pars) {
return 0.25*x - x*x*x + (pars->gamma - pars->epsilon*(x*x))*v;
}
void algEulero(point *xv, parametri *pars, double *xn1, double *vn1){
printf("%.14g\t%.14g\t%.14g\t%.14g\t%.14g\n", xv->t, xv->x, xv->v, xv->E, xv->dE);
*xn1 = xv->x + (xv->v)*(pars->dt);
*vn1 = xv->v + f(xv->x, xv->v, pars)*(pars->dt);
xv->E = 0.5*(xv->v)*(xv->v) - (xv->x)*(xv->x)/8. + (xv->x)*(xv->x)*(xv->x)*(xv->x)/4.;
xv->dE = fabs(xv->E-xv->E0)/xv->E0;
xv->t += (pars->dt);
xv->x = *xn1;
xv->v = *vn1;
}
void algEuleroCromer(point *xv, parametri *pars, double *xn1, double *vn1){
printf("%.14g\t%.14g\t%.14g\t%.14g\t%.14g\n", xv->t, xv->x, xv->v, xv->E, xv->dE);
xv->v = xv->v + f(xv->x, xv->v, pars)*(pars->dt);
xv->x = xv->x + (xv->v)*(pars->dt);
xv->E = 0.5*(xv->v)*(xv->v) - (xv->x)*(xv->x)/8. + (xv->x)*(xv->x)*(xv->x)*(xv->x)/4.;
xv->dE = fabs(xv->E-xv->E0)/xv->E0;
xv->t += (pars->dt);
}
void algPuntoDiMezzo(point *xv, parametri *pars, double *xn1, double *vn1) {
printf("%.14g\t%.14g\t%.14g\t%.14g\t%.14g\n", xv->t, xv->x, xv->v, xv->E, xv->dE);
*vn1 = xv->v + f(xv->x, xv->v, pars)*(pars->dt);
xv->x = xv->x + (0.5*(xv->v + (*vn1)))*(pars->dt);
xv->E = 0.5*(xv->v)*(xv->v) - (xv->x)*(xv->x)/8. + (xv->x)*(xv->x)*(xv->x)*(xv->x)/4.;
xv->dE = fabs(xv->E-xv->E0)/xv->E0;
xv->t += (pars->dt);
xv->v = *vn1;
}
void algVerlet(point *xv, parametri *pars, double *xn1, double *vn1) {
printf("%.14g\t%.14g\t%.14g\t%.14g\t%.14g\n", xv->t, xv->x, xv->v, xv->E, xv->dE);
*xn1 = xv->x + xv->v*pars->dt + 0.5*(f(xv->x, xv->v, pars))*pars->dt*pars->dt;
*vn1 = xv->v + 0.5*(f(xv->x, xv->v, pars) + f((*xn1), xv->v, pars))*pars->dt;
xv->x = *xn1;
xv->v = *vn1;
xv->E = 0.5*(xv->v)*(xv->v) - (xv->x)*(xv->x)/8. + (xv->x)*(xv->x)*(xv->x)*(xv->x)/4.;
xv->dE = fabs(xv->E-xv->E0)/xv->E0;
xv->t += (pars->dt);
}
void algRungeKutta2(point *xv, parametri *pars, double *xn1, double *vn1) {
printf("%.14g\t%.14g\t%.14g\t%.14g\t%.14g\n", xv->t, xv->x, xv->v, xv->E, xv->dE);
*xn1 = xv->x + xv->v*pars->dt + 0.5*f(xv->x, xv->v, pars)*pars->dt*pars->dt;
*vn1 = xv->v + f(xv->x + 0.5*xv->v*pars->dt, xv->v + 0.5*f(xv->x, xv->v, pars)*pars->dt, pars)*pars->dt;
xv->x = *xn1;
xv->v = *vn1;
xv->E = 0.5*(xv->v)*(xv->v) - (xv->x)*(xv->x)/8. + (xv->x)*(xv->x)*(xv->x)*(xv->x)/4.;
xv->dE = fabs(xv->E-xv->E0)/xv->E0;
xv->t += (pars->dt);
}
----------------编辑----------------
亲爱的吉尔斯，
我向您解释我的程序是做什么的。该程序旨在用不同的算法（algEulero、algEuleroCromer、algPuntoDiMezzo、algVerlet、algRungeKutta2）对微分方程进行数值求解。它工作正常。物理方程是 d^2x/dt^2 = f(x, v, gamma, epsilon)。这个 f() 正是您可以在我的代码中找到的 f() 。我的简单 C 程序的 "big" 问题是他真的很慢：当我选择 5' 练习 (pars.choice == 5) 时，该程序将完全执行：
1) 计算（getBasin()中的两个for）[-1, 1]x[-1, 1]的面积，增量为pars.dx；
2) 在 x 和 y 上的每个增量都将启动算法，该算法以 for 的两个起始条件 (x, y) 对微分方程进行数值求解。当算法达到一个异步t*(pars.t_basin)时，他会把x(t*)和v(t*)的输出写到basin.dat中。
3) (x, y) 将改变并再次转到 (1) 点。
现在，您可以使用以下参数进行测试：0.83、4、0.01、62、1、1、5、5。
我测试了你的代码，但并不比我的顺序代码快（例如超过 11 分钟）。为了改进它：
1) basin.dat 的输出顺序无关紧要，因为我将根据 3' 和 4' 列值绘制 space (x, y) 着色我的点（与Gnuplot 上的图像）。
2) 为什么还要在函数 getBasin() 之前创建线程，而不仅仅是为两个 for() 创建线程？对于每个线程，这不应该重复 getBasin() 吗？
抱歉我的英语和并行编程不好，我正在努力改进它在线阅读教程。
好吧，您的代码中明显的主要问题是并行版本（非常）错误。您在 parallel 区域之外定义了所有变量，但不声明其中任何一个 private（甚至循环索引也不声明）。
此外，由于您将 getBasin() 函数的所有参数都作为指针传递，之后将它们声明为私有变得更加棘手。
然而，对代码的快速分析表明，尽管这些参数作为指针传递，但您实际上并不关心它们在例程退出时的值。此外，看起来您尝试并行化的循环没有数据依赖性（尽管我没有对此进行全面检查）。我发现的唯一明确的依赖关系是关于你打开的输出文件，在保持顺序的同时并行更新会很棘手...
所以为了修复你的代码的并行化，这就是我所做的：
- 通过值而不是通过引用将参数传递给
getBasin()。这将生成一些额外的数据副本，但由于这些数据需要声明 private，因此无论如何都需要一份副本。
- 在
parallel 区域内包含对 getBasin() 的调用。事实上，在我看来，这样做比在函数本身内部处理数据私有化和 IO 问题更简单。
- 为每个线程打开一个输出文件的私有副本，名为 "basinXX.dat"，"XX" 是当前线程的 ID，左填充 0。这些文件将包含共享当前线程对应的全局输出。希望这适合您。否则，处理打印顺序可能会有点棘手。
- 在函数内使用一个孤立的
omp for 指令来并行化循环。
有了这个，代码应该（希望）扩展得更好。但是，由于您没有指明用于计算的输入参数，因此我无法对其进行测试，无论是正确性还是性能。所以它可能会惨败...
无论如何，这是我想出的：
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <omp.h>
#define BLD "\x1B[1m"
#define RED "\x1B[31m"
#define GRN "\x1B[32m"
#define RST "\x1B[0m"
struct point{
double x;
double v;
double t;
double E0;
double E;
double dE;
} typedef point;
struct parametri {
double gamma;
double epsilon;
double dt;
double t_star;
double t_basin;
double t_max;
double x0;
double v0;
int choice;
int alg[1];
double dx;
} typedef parametri;
// Prototipi delle funzioni
void Setup();
void getParams(parametri *pars);
void getError(point *xv, parametri *pars, int *i, int k, int *n_passi, double *xn1, double *vn1);
void getBasin(point xv, parametri pars, int h, int n_passi, double xn1, double vn1);
double f(double x, double v, parametri *pars);
void algEulero(point *xv, parametri *pars, double *xn1, double *vn1);
void algEuleroCromer(point *xv, parametri *pars, double *xn1, double *vn1);
void algPuntoDiMezzo(point *xv, parametri *pars, double *xn1, double *vn1);
void algVerlet(point *xv, parametri *pars, double *xn1, double *vn1);
void algRungeKutta2(point *xv, parametri *pars, double *xn1, double *vn1);
int main(void) {
// Inizializzo il display per l'utente. Maggiori informazioni vedere funzione Setup();
Setup();
// Dichiaro e recupero i parametri richiesti per il corretto funzionamento del programma;
parametri pars;
getParams(&pars);
// Dichiaro e ricavo i parametri essenziali, annesse variabili di supporto;
int i, n_passi = pars.t_max/pars.dt;
double dt0, xn1, vn1, t_star = 5.;
point xv;
// Imposto i parametri iniziali;
xv.x = pars.x0;
xv.v = pars.v0;
xv.E0 = 0.5*(xv.v)*(xv.v) - (xv.x)*(xv.x)/8. + (xv.x)*(xv.x)*(xv.x)*(xv.x)/4.;
xv.E = xv.E0;
xv.dE = 0;
xv.t = 0;
pars.t_star = 5;
pars.t_basin = 60;
dt0 = pars.dt;
// Formato dell'output;
printf ("t\tx(t)\tv(t)\tE(t)\tdE\n");
if ((pars.choice == 1) || (pars.choice == 3) || (pars.choice == 4)) { // L'utente ha deciso di affrontare il primo/terzo/quarto esercizio;
// Avverto l'utente che il processo sta per iniziare. Può risultare inutile su tempi brevi, ma efficace su tempi molto lunghi;
fprintf(stderr, "\nAvvio integrazione numerica... ");
if (pars.alg[0] == 1) { // L'utente ha selezionato l'algoritmo di Eulero;
for (i=0; i<n_passi; i++){
algEulero(&xv, &pars, &xn1, &vn1);
}
} else if (pars.alg[0] == 2) { // L'utente ha selezionato l'algoritmo di EuleroCromer;
for (i=0; i<n_passi; i++){
algEuleroCromer(&xv, &pars, &xn1, &vn1);
}
} else if (pars.alg[0] == 3) { // L'utente ha selezionato l'algoritmo di PuntoDiMezzo;
for (i=0; i<n_passi; i++){
algPuntoDiMezzo(&xv, &pars, &xn1, &vn1);
}
} else if (pars.alg[0] == 4) { // L'utente ha selezionato l'algoritmo di Verlet;
for (i=0; i<n_passi; i++){
algVerlet(&xv, &pars, &xn1, &vn1);
}
} else if (pars.alg[0] == 5) { // L'utente ha selezionato l'algoritmo di RungeKutta di ordine 2;
for (i=0; i<n_passi; i++) {
algRungeKutta2(&xv, &pars, &xn1, &vn1);
}
}
// Algoritmo terminato;
fprintf(stderr, "[%s%sDONE%s]\n\n", BLD, GRN, RST);
} else if (pars.choice == 2) { // L'utente ha deciso di affrontare il secondo esercizio;
// Seleziono il secondo algoritmo da confrontare
do {
fprintf(stderr, "> Selezionare il secondo algoritmo:\n");
fprintf(stderr, "\t>> [1] Eulero\n");
fprintf(stderr, "\t>> [2] EuleroCromer - RungeKutta (1 ordine)\n");
fprintf(stderr, "\t>> [3] PuntoDiMezzo\n");
fprintf(stderr, "\t>> [4] Verlet\n");
fprintf(stderr, "\t>> [5] RungeKutta (2 ordine)\n");
fprintf(stderr, "\t>> ");
scanf("%d", &pars.alg[1]);
} while (( pars.alg[1] <= 0 ));
// Avverto l'utente che il processo sta per iniziare. Può risultare inutile su tempi brevi, ma efficace su tempi molto lunghi;
fprintf(stderr, "\nAvvio calcolo errori... ");
// Eseguo lo studio degli errori d'integrazione mediante il primo e secondo algoritmo scelto, rispettivamente nei file error_1.dat, error_2.dat;
getError(&xv, &pars, &i, 0, &n_passi, &xn1, &vn1);
// Resetto le variabili e richiamo l'algoritmo per calcolare gli errori;
pars.dt = dt0;
n_passi = pars.t_max/pars.dt;
xv.x = pars.x0;
xv.v = pars.v0;
xv.E = xv.E0;
xv.dE = 0;
xv.t = 0;
getError(&xv, &pars, &i, 1, &n_passi, &xn1, &vn1);
// Processo terminato;
fprintf(stderr, "\n\nStato: [%s%sDONE%s]\n\n", BLD, GRN, RST);
} else if (pars.choice == 5) { // L'utente ha deciso di affrontare il quinto esercizio;
// Avverto l'utente che il processo sta per iniziare. Può risultare inutile su tempi brevi, ma efficace su tempi molto lunghi;
fprintf(stderr, "\nAvvio calcolo griglia... ");
#pragma omp parallel num_threads( 4 )
getBasin(xv, pars, i, n_passi, xn1, vn1);
// Processo terminato;
fprintf(stderr, "[%s%sDONE%s]\n\n", BLD, GRN, RST);
} else { // L'utente non ha selezionato un esercizio valido;
fprintf(stderr, "\n[%s%sFAILED%s] Esercizio non disponibile.\n", BLD, RED, RST);
exit(EXIT_FAILURE);
}
return 0;
}
void Setup() {
fprintf(stderr, "\nAnalisi numerica di un'equazione differenziale\n");
fprintf(stderr, "==============================================\n\n");
}
void getParams(parametri *pars) {
do {
fprintf(stderr, "> Inserire un valore per gamma : ");
scanf("%lf", &pars->gamma);
} while (pars->gamma < 0);
do {
fprintf(stderr, "> Inserire un valore per epsilon : ");
scanf("%lf", &pars->epsilon);
} while (pars->epsilon < 0);
do {
fprintf(stderr, "> Inserire un valore per dt : ");
scanf("%lf", &pars->dt);
} while (pars->dt < 0);
do {
fprintf(stderr, "> Inserire un valore per tmax t.c. :\n");
fprintf(stderr, " >> (tmax > 5) per I eserc.\n");
fprintf(stderr, " >> (tmax > 60) per V eserc. : ");
scanf("%lf", &pars->t_max);
} while (pars->t_max < 0);
do {
fprintf(stderr, "> Inserire un valore per x(0) : ");
scanf("%lf", &pars->x0);
} while (pars->x0 < -1);
do {
fprintf(stderr, "> Inserire un valore per v(0) : ");
scanf("%lf", &pars->v0);
} while (pars->v0 < -1);
do {
fprintf(stderr, "> Selezionare l'esercizio richiesto :\n");
fprintf(stderr, "\t>> [1] Esercizio 1\n");
fprintf(stderr, "\t>> [2] Esercizio 2\n");
fprintf(stderr, "\t>> [3] Esercizio 3\n");
fprintf(stderr, "\t>> [4] Esercizio 4\n");
fprintf(stderr, "\t>> [5] Esercizio 5\n\n");
fprintf(stderr, "\t>> ");
scanf("%d", &pars->choice);
} while (( pars->choice <= 0 ));
do {
fprintf(stderr, "> Selezionare l'algoritmo voluto :\n");
fprintf(stderr, "\t>> [1] Eulero\n");
fprintf(stderr, "\t>> [2] EuleroCromer - RungeKutta (1 ordine)\n");
fprintf(stderr, "\t>> [3] PuntoDiMezzo\n");
fprintf(stderr, "\t>> [4] Verlet\n");
fprintf(stderr, "\t>> [5] RungeKutta (2 ordine)\n\n");
fprintf(stderr, "\t>> ");
scanf("%d", &pars->alg[0]);
} while (( pars->alg[0] <= 0 ));
}
void getError(point *xv, parametri *pars, int *i, int k, int *n_passi, double *xn1, double *vn1) {
void (*algF)(point *, parametri *, double *, double *);
int j, n = 0;
FILE *fp;
// Questo if controlla se il codice sta eseguendo lo studio degli errori per il primo o secondo algoritmo (pars->alg[k]);
if (k == 0) {
fp = fopen("error_1.dat", "w+");
} else {
fp = fopen("error_2.dat", "w+");
}
// Assegno il puntatore corretto a seconda dell'algoritmo scelto;
if (pars->alg[k] == 1) {
algF = algEulero;
} else if (pars->alg[k] == 2) {
algF = algEuleroCromer;
} else if (pars->alg[k] == 3) {
algF = algPuntoDiMezzo;
} else if (pars->alg[k] == 4) {
algF = algVerlet;
} else if (pars->alg[k] == 5) {
algF = algRungeKutta2;
} else {
fprintf(stderr, "\n[%s%sFAILED%s] E' stato selezionato un algoritmo non valido.\n", BLD, RED, RST);
exit(EXIT_FAILURE);
}
// Informo l'utente che il processo sta iniziando;
fprintf(stderr, "\n>> Avvio %d algoritmo... ", k+1);
// Formattazione dell'output del file contenente gli errori;
fprintf(fp, "dt\tE(t*)\tE(0)\tdE/E(0)\t# passi\ti\tt\n");
for (j=0; j<8; j++) {
for ((*i)=0; (*i)<(*n_passi); (*i)++){
(*algF)(xv, pars, xn1, vn1);
if (((pars->t_star - pars->dt/2. <= xv->t) && (xv->t >= pars->t_star + pars->dt/2.)) && (n == 0)) {
fprintf(fp, "%+.14lf\t%+.14lf\t%+.14lf\t%+.14lf\t%+.14d\t%+.14d\t%+.14lf\n", pars->dt, xv->E, xv->E0, (xv->E - xv->E0)/xv->E0, (*n_passi), (*i), xv->t);
n = 1;
}
}
// Resetto le variabili per rilanciare l'algoritmo con dt/2^j
n = 0;
xv->t = 0;
xv->x = pars->x0;
xv->v = pars->v0;
pars->dt = pars->dt/2.;
(*n_passi) = pars->t_max/pars->dt;
(*xn1) = 0;
(*vn1) = 0;
xv->E = xv->E0;
}
fclose(fp);
fprintf(stderr, "[%s%sDONE%s]", BLD, GRN, RST);
}
void getBasin(point xv, parametri pars, int h, int n_passi, double xn1, double vn1) {
// Dichiaro e setto i parametri che delimitano la griglia;
point xv1;
pars.dx = 0.01;
xv1.x = -1;
xv1.v = -1;
// Dichiaro le variabili necessarie per il bacino d'attrazione;
int i, j, i_max = 2./pars.dx, j_max = 2./pars.dx;
void (*algF)(point *, parametri *, double *, double *);
char fname[13];
sprintf( fname, "bassin%02d.dat", omp_get_thread_num() );
FILE *fp = fopen( fname, "w+");
// Assegno il puntatore corretto a seconda dell'algoritmo scelto;
if (pars.alg[0] == 1) {
algF = algEulero;
} else if (pars.alg[0] == 2) {
algF = algEuleroCromer;
} else if (pars.alg[0] == 3) {
algF = algPuntoDiMezzo;
} else if (pars.alg[0] == 4) {
algF = algVerlet;
} else if (pars.alg[0] == 5) {
algF = algRungeKutta2;
} else {
fprintf(stderr, "\n[%s%sFAILED%s] E' stato selezionato un algoritmo non valido.\n", BLD, RED, RST);
exit(EXIT_FAILURE);
}
// Formattazione output file basin.dat;
fprintf(fp, "x(0)\tx'(0)\tx(t*)\tv(t*)\n");
#pragma omp for
// Eseguo il for della griglia sull'asse x';
for (j=0; j<=j_max; j++) {
// Eseguo il for della griglia sull'asse x;
for (i=0; i<=i_max; i++) {
fprintf(fp, "%lf\t%lf\t", xv1.x, xv1.v);
xv.x = xv1.x;
xv.v = xv1.v;
// Eseguo l'integrazione numerica
for (h=0; h<n_passi; h++) {
(*algF)(&xv, &pars, &xn1, &vn1);
// Entro in t = t*, stampo v(t*) ed esco dal ciclo;
if ((pars.t_basin - pars.dt/2. <= xv.t) && (xv.t >= pars.t_basin + pars.dt/2.)) {
fprintf(fp, "%lf\t%lf\n", xv.x, xv.v);
break;
}
}
xv1.x += pars.dx;
xv.t = 0;
xn1 = 0;
vn1 = 0;
}
// Resetto la x e incremento la x';
xv1.x = -1;
xv1.v += pars.dx;
}
}
double f(double x, double v, parametri *pars) {
return 0.25*x - x*x*x + (pars->gamma - pars->epsilon*(x*x))*v;
}
void algEulero(point *xv, parametri *pars, double *xn1, double *vn1){
//printf("%.14g\t%.14g\t%.14g\t%.14g\t%.14g\n", xv->t, xv->x, xv->v, xv->E, xv->dE);
*xn1 = xv->x + (xv->v)*(pars->dt);
*vn1 = xv->v + f(xv->x, xv->v, pars)*(pars->dt);
xv->E = 0.5*(xv->v)*(xv->v) - (xv->x)*(xv->x)/8. + (xv->x)*(xv->x)*(xv->x)*(xv->x)/4.;
xv->dE = fabs(xv->E-xv->E0)/xv->E0;
xv->t += (pars->dt);
xv->x = *xn1;
xv->v = *vn1;
}
void algEuleroCromer(point *xv, parametri *pars, double *xn1, double *vn1){
//printf("%.14g\t%.14g\t%.14g\t%.14g\t%.14g\n", xv->t, xv->x, xv->v, xv->E, xv->dE);
xv->v = xv->v + f(xv->x, xv->v, pars)*(pars->dt);
xv->x = xv->x + (xv->v)*(pars->dt);
xv->E = 0.5*(xv->v)*(xv->v) - (xv->x)*(xv->x)/8. + (xv->x)*(xv->x)*(xv->x)*(xv->x)/4.;
xv->dE = fabs(xv->E-xv->E0)/xv->E0;
xv->t += (pars->dt);
}
void algPuntoDiMezzo(point *xv, parametri *pars, double *xn1, double *vn1) {
//printf("%.14g\t%.14g\t%.14g\t%.14g\t%.14g\n", xv->t, xv->x, xv->v, xv->E, xv->dE);
*vn1 = xv->v + f(xv->x, xv->v, pars)*(pars->dt);
xv->x = xv->x + (0.5*(xv->v + (*vn1)))*(pars->dt);
xv->E = 0.5*(xv->v)*(xv->v) - (xv->x)*(xv->x)/8. + (xv->x)*(xv->x)*(xv->x)*(xv->x)/4.;
xv->dE = fabs(xv->E-xv->E0)/xv->E0;
xv->t += (pars->dt);
xv->v = *vn1;
}
void algVerlet(point *xv, parametri *pars, double *xn1, double *vn1) {
//printf("%.14g\t%.14g\t%.14g\t%.14g\t%.14g\n", xv->t, xv->x, xv->v, xv->E, xv->dE);
*xn1 = xv->x + xv->v*pars->dt + 0.5*(f(xv->x, xv->v, pars))*pars->dt*pars->dt;
*vn1 = xv->v + 0.5*(f(xv->x, xv->v, pars) + f((*xn1), xv->v, pars))*pars->dt;
xv->x = *xn1;
xv->v = *vn1;
xv->E = 0.5*(xv->v)*(xv->v) - (xv->x)*(xv->x)/8. + (xv->x)*(xv->x)*(xv->x)*(xv->x)/4.;
xv->dE = fabs(xv->E-xv->E0)/xv->E0;
xv->t += (pars->dt);
}
void algRungeKutta2(point *xv, parametri *pars, double *xn1, double *vn1) {
//printf("%.14g\t%.14g\t%.14g\t%.14g\t%.14g\n", xv->t, xv->x, xv->v, xv->E, xv->dE);
*xn1 = xv->x + xv->v*pars->dt + 0.5*f(xv->x, xv->v, pars)*pars->dt*pars->dt;
*vn1 = xv->v + f(xv->x + 0.5*xv->v*pars->dt, xv->v + 0.5*f(xv->x, xv->v, pars)*pars->dt, pars)*pars->dt;
xv->x = *xn1;
xv->v = *vn1;
xv->E = 0.5*(xv->v)*(xv->v) - (xv->x)*(xv->x)/8. + (xv->x)*(xv->x)*(xv->x)*(xv->x)/4.;
xv->dE = fabs(xv->E-xv->E0)/xv->E0;
xv->t += (pars->dt);
}
编辑
使用您提供的新输入参数，我能够测试代码，结果是：
- 我给你的代码中有一个小错误（用于打开输出文件，使用
"fname" 而不是 fname）
- 您的代码大部分（全部）时间都花在
printf()。
修复这两个后，编译如下：
gcc -O3 -fopenmp -mtune=native -march=native -w bassin.c
我 运行 它在我的双核笔记本电脑上 2.12s
代码已更新。请自行尝试。
我正在尝试改进我的 C 源代码以并行执行。我有一个四核 CPU，所以我认为 4 个线程（一个用于 CPU）对 运行 我的程序来说真的比像顺序代码那样优化更快。
但它不起作用。我的没有 OpenMP 的代码需要 11 分钟才能执行，而并行代码需要超过 11 分钟。我报告了我的所有来源，但并行代码仅在 getBasin() 函数中。
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <omp.h>
#define BLD "\x1B[1m"
#define RED "\x1B[31m"
#define GRN "\x1B[32m"
#define RST "\x1B[0m"
struct point{
double x;
double v;
double t;
double E0;
double E;
double dE;
} typedef point;
struct parametri {
double gamma;
double epsilon;
double dt;
double t_star;
double t_basin;
double t_max;
double x0;
double v0;
int choice;
int alg[1];
double dx;
} typedef parametri;
// Prototipi delle funzioni
void Setup();
void getParams(parametri *pars);
void getError(point *xv, parametri *pars, int *i, int k, int *n_passi, double *xn1, double *vn1);
void getBasin(point *xv, parametri *pars, int *h, int *n_passi, double *xn1, double *vn1);
double f(double x, double v, parametri *pars);
void algEulero(point *xv, parametri *pars, double *xn1, double *vn1);
void algEuleroCromer(point *xv, parametri *pars, double *xn1, double *vn1);
void algPuntoDiMezzo(point *xv, parametri *pars, double *xn1, double *vn1);
void algVerlet(point *xv, parametri *pars, double *xn1, double *vn1);
void algRungeKutta2(point *xv, parametri *pars, double *xn1, double *vn1);
int main(void) {
// Inizializzo il display per l'utente. Maggiori informazioni vedere funzione Setup();
Setup();
// Dichiaro e recupero i parametri richiesti per il corretto funzionamento del programma;
parametri pars;
getParams(&pars);
// Dichiaro e ricavo i parametri essenziali, annesse variabili di supporto;
int i, n_passi = pars.t_max/pars.dt;
double dt0, xn1, vn1, t_star = 5.;
point xv;
// Imposto i parametri iniziali;
xv.x = pars.x0;
xv.v = pars.v0;
xv.E0 = 0.5*(xv.v)*(xv.v) - (xv.x)*(xv.x)/8. + (xv.x)*(xv.x)*(xv.x)*(xv.x)/4.;
xv.E = xv.E0;
xv.dE = 0;
xv.t = 0;
pars.t_star = 5;
pars.t_basin = 60;
dt0 = pars.dt;
// Formato dell'output;
printf ("t\tx(t)\tv(t)\tE(t)\tdE\n");
if ((pars.choice == 1) || (pars.choice == 3) || (pars.choice == 4)) { // L'utente ha deciso di affrontare il primo/terzo/quarto esercizio;
// Avverto l'utente che il processo sta per iniziare. Può risultare inutile su tempi brevi, ma efficace su tempi molto lunghi;
fprintf(stderr, "\nAvvio integrazione numerica... ");
if (pars.alg[0] == 1) { // L'utente ha selezionato l'algoritmo di Eulero;
for (i=0; i<n_passi; i++){
algEulero(&xv, &pars, &xn1, &vn1);
}
} else if (pars.alg[0] == 2) { // L'utente ha selezionato l'algoritmo di EuleroCromer;
for (i=0; i<n_passi; i++){
algEuleroCromer(&xv, &pars, &xn1, &vn1);
}
} else if (pars.alg[0] == 3) { // L'utente ha selezionato l'algoritmo di PuntoDiMezzo;
for (i=0; i<n_passi; i++){
algPuntoDiMezzo(&xv, &pars, &xn1, &vn1);
}
} else if (pars.alg[0] == 4) { // L'utente ha selezionato l'algoritmo di Verlet;
for (i=0; i<n_passi; i++){
algVerlet(&xv, &pars, &xn1, &vn1);
}
} else if (pars.alg[0] == 5) { // L'utente ha selezionato l'algoritmo di RungeKutta di ordine 2;
for (i=0; i<n_passi; i++) {
algRungeKutta2(&xv, &pars, &xn1, &vn1);
}
}
// Algoritmo terminato;
fprintf(stderr, "[%s%sDONE%s]\n\n", BLD, GRN, RST);
} else if (pars.choice == 2) { // L'utente ha deciso di affrontare il secondo esercizio;
// Seleziono il secondo algoritmo da confrontare
do {
fprintf(stderr, "> Selezionare il secondo algoritmo:\n");
fprintf(stderr, "\t>> [1] Eulero\n");
fprintf(stderr, "\t>> [2] EuleroCromer - RungeKutta (1 ordine)\n");
fprintf(stderr, "\t>> [3] PuntoDiMezzo\n");
fprintf(stderr, "\t>> [4] Verlet\n");
fprintf(stderr, "\t>> [5] RungeKutta (2 ordine)\n");
fprintf(stderr, "\t>> ");
scanf("%d", &pars.alg[1]);
} while (( pars.alg[1] <= 0 ));
// Avverto l'utente che il processo sta per iniziare. Può risultare inutile su tempi brevi, ma efficace su tempi molto lunghi;
fprintf(stderr, "\nAvvio calcolo errori... ");
// Eseguo lo studio degli errori d'integrazione mediante il primo e secondo algoritmo scelto, rispettivamente nei file error_1.dat, error_2.dat;
getError(&xv, &pars, &i, 0, &n_passi, &xn1, &vn1);
// Resetto le variabili e richiamo l'algoritmo per calcolare gli errori;
pars.dt = dt0;
n_passi = pars.t_max/pars.dt;
xv.x = pars.x0;
xv.v = pars.v0;
xv.E = xv.E0;
xv.dE = 0;
xv.t = 0;
getError(&xv, &pars, &i, 1, &n_passi, &xn1, &vn1);
// Processo terminato;
fprintf(stderr, "\n\nStato: [%s%sDONE%s]\n\n", BLD, GRN, RST);
} else if (pars.choice == 5) { // L'utente ha deciso di affrontare il quinto esercizio;
// Avverto l'utente che il processo sta per iniziare. Può risultare inutile su tempi brevi, ma efficace su tempi molto lunghi;
fprintf(stderr, "\nAvvio calcolo griglia... ");
getBasin(&xv, &pars, &i, &n_passi, &xn1, &vn1);
// Processo terminato;
fprintf(stderr, "[%s%sDONE%s]\n\n", BLD, GRN, RST);
} else { // L'utente non ha selezionato un esercizio valido;
fprintf(stderr, "\n[%s%sFAILED%s] Esercizio non disponibile.\n", BLD, RED, RST);
exit(EXIT_FAILURE);
}
return 0;
}
void Setup() {
fprintf(stderr, "\nAnalisi numerica di un'equazione differenziale\n");
fprintf(stderr, "==============================================\n\n");
}
void getParams(parametri *pars) {
do {
fprintf(stderr, "> Inserire un valore per gamma : ");
scanf("%lf", &pars->gamma);
} while (pars->gamma < 0);
do {
fprintf(stderr, "> Inserire un valore per epsilon : ");
scanf("%lf", &pars->epsilon);
} while (pars->epsilon < 0);
do {
fprintf(stderr, "> Inserire un valore per dt : ");
scanf("%lf", &pars->dt);
} while (pars->dt < 0);
do {
fprintf(stderr, "> Inserire un valore per tmax t.c. :\n");
fprintf(stderr, " >> (tmax > 5) per I eserc.\n");
fprintf(stderr, " >> (tmax > 60) per V eserc. : ");
scanf("%lf", &pars->t_max);
} while (pars->t_max < 0);
do {
fprintf(stderr, "> Inserire un valore per x(0) : ");
scanf("%lf", &pars->x0);
} while (pars->x0 < -1);
do {
fprintf(stderr, "> Inserire un valore per v(0) : ");
scanf("%lf", &pars->v0);
} while (pars->v0 < -1);
do {
fprintf(stderr, "> Selezionare l'esercizio richiesto :\n");
fprintf(stderr, "\t>> [1] Esercizio 1\n");
fprintf(stderr, "\t>> [2] Esercizio 2\n");
fprintf(stderr, "\t>> [3] Esercizio 3\n");
fprintf(stderr, "\t>> [4] Esercizio 4\n");
fprintf(stderr, "\t>> [5] Esercizio 5\n\n");
fprintf(stderr, "\t>> ");
scanf("%d", &pars->choice);
} while (( pars->choice <= 0 ));
do {
fprintf(stderr, "> Selezionare l'algoritmo voluto :\n");
fprintf(stderr, "\t>> [1] Eulero\n");
fprintf(stderr, "\t>> [2] EuleroCromer - RungeKutta (1 ordine)\n");
fprintf(stderr, "\t>> [3] PuntoDiMezzo\n");
fprintf(stderr, "\t>> [4] Verlet\n");
fprintf(stderr, "\t>> [5] RungeKutta (2 ordine)\n\n");
fprintf(stderr, "\t>> ");
scanf("%d", &pars->alg[0]);
} while (( pars->alg[0] <= 0 ));
}
void getError(point *xv, parametri *pars, int *i, int k, int *n_passi, double *xn1, double *vn1) {
void (*algF)(point *, parametri *, double *, double *);
int j, n = 0;
FILE *fp;
// Questo if controlla se il codice sta eseguendo lo studio degli errori per il primo o secondo algoritmo (pars->alg[k]);
if (k == 0) {
fp = fopen("error_1.dat", "w+");
} else {
fp = fopen("error_2.dat", "w+");
}
// Assegno il puntatore corretto a seconda dell'algoritmo scelto;
if (pars->alg[k] == 1) {
algF = algEulero;
} else if (pars->alg[k] == 2) {
algF = algEuleroCromer;
} else if (pars->alg[k] == 3) {
algF = algPuntoDiMezzo;
} else if (pars->alg[k] == 4) {
algF = algVerlet;
} else if (pars->alg[k] == 5) {
algF = algRungeKutta2;
} else {
fprintf(stderr, "\n[%s%sFAILED%s] E' stato selezionato un algoritmo non valido.\n", BLD, RED, RST);
exit(EXIT_FAILURE);
}
// Informo l'utente che il processo sta iniziando;
fprintf(stderr, "\n>> Avvio %d algoritmo... ", k+1);
// Formattazione dell'output del file contenente gli errori;
fprintf(fp, "dt\tE(t*)\tE(0)\tdE/E(0)\t# passi\ti\tt\n");
for (j=0; j<8; j++) {
for ((*i)=0; (*i)<(*n_passi); (*i)++){
(*algF)(xv, pars, xn1, vn1);
if (((pars->t_star - pars->dt/2. <= xv->t) && (xv->t >= pars->t_star + pars->dt/2.)) && (n == 0)) {
fprintf(fp, "%+.14lf\t%+.14lf\t%+.14lf\t%+.14lf\t%+.14d\t%+.14d\t%+.14lf\n", pars->dt, xv->E, xv->E0, (xv->E - xv->E0)/xv->E0, (*n_passi), (*i), xv->t);
n = 1;
}
}
// Resetto le variabili per rilanciare l'algoritmo con dt/2^j
n = 0;
xv->t = 0;
xv->x = pars->x0;
xv->v = pars->v0;
pars->dt = pars->dt/2.;
(*n_passi) = pars->t_max/pars->dt;
(*xn1) = 0;
(*vn1) = 0;
xv->E = xv->E0;
}
fclose(fp);
fprintf(stderr, "[%s%sDONE%s]", BLD, GRN, RST);
}
void getBasin(point *xv, parametri *pars, int *h, int *n_passi, double *xn1, double *vn1) {
// Dichiaro e setto i parametri che delimitano la griglia;
point xv1;
pars->dx = 0.01;
xv1.x = -1;
xv1.v = -1;
// Dichiaro le variabili necessarie per il bacino d'attrazione;
int i, j, i_max = 2./pars->dx, j_max = 2./pars->dx;
void (*algF)(point *, parametri *, double *, double *);
FILE *fp = fopen("basin.dat", "w+");
// Assegno il puntatore corretto a seconda dell'algoritmo scelto;
if (pars->alg[0] == 1) {
algF = algEulero;
} else if (pars->alg[0] == 2) {
algF = algEuleroCromer;
} else if (pars->alg[0] == 3) {
algF = algPuntoDiMezzo;
} else if (pars->alg[0] == 4) {
algF = algVerlet;
} else if (pars->alg[0] == 5) {
algF = algRungeKutta2;
} else {
fprintf(stderr, "\n[%s%sFAILED%s] E' stato selezionato un algoritmo non valido.\n", BLD, RED, RST);
exit(EXIT_FAILURE);
}
// Formattazione output file basin.dat;
fprintf(fp, "x(0)\tx'(0)\tx(t*)\tv(t*)\n");
omp_set_num_threads(4);
#pragma omp parallel for
// Eseguo il for della griglia sull'asse x';
for (j=0; j<=j_max; j++) {
// Eseguo il for della griglia sull'asse x;
for (i=0; i<=i_max; i++) {
fprintf(fp, "%lf\t%lf\t", xv1.x, xv1.v);
xv->x = xv1.x;
xv->v = xv1.v;
// Eseguo l'integrazione numerica
for ((*h)=0; (*h)<(*n_passi); (*h)++) {
(*algF)(xv, pars, xn1, vn1);
// Entro in t = t*, stampo v(t*) ed esco dal ciclo;
if ((pars->t_basin - pars->dt/2. <= xv->t) && (xv->t >= pars->t_basin + pars->dt/2.)) {
fprintf(fp, "%lf\t%lf\n", xv->x, xv->v);
break;
}
}
xv1.x += pars->dx;
xv->t = 0;
(*xn1) = 0;
(*vn1) = 0;
}
// Resetto la x e incremento la x';
xv1.x = -1;
xv1.v += pars->dx;
}
}
double f(double x, double v, parametri *pars) {
return 0.25*x - x*x*x + (pars->gamma - pars->epsilon*(x*x))*v;
}
void algEulero(point *xv, parametri *pars, double *xn1, double *vn1){
printf("%.14g\t%.14g\t%.14g\t%.14g\t%.14g\n", xv->t, xv->x, xv->v, xv->E, xv->dE);
*xn1 = xv->x + (xv->v)*(pars->dt);
*vn1 = xv->v + f(xv->x, xv->v, pars)*(pars->dt);
xv->E = 0.5*(xv->v)*(xv->v) - (xv->x)*(xv->x)/8. + (xv->x)*(xv->x)*(xv->x)*(xv->x)/4.;
xv->dE = fabs(xv->E-xv->E0)/xv->E0;
xv->t += (pars->dt);
xv->x = *xn1;
xv->v = *vn1;
}
void algEuleroCromer(point *xv, parametri *pars, double *xn1, double *vn1){
printf("%.14g\t%.14g\t%.14g\t%.14g\t%.14g\n", xv->t, xv->x, xv->v, xv->E, xv->dE);
xv->v = xv->v + f(xv->x, xv->v, pars)*(pars->dt);
xv->x = xv->x + (xv->v)*(pars->dt);
xv->E = 0.5*(xv->v)*(xv->v) - (xv->x)*(xv->x)/8. + (xv->x)*(xv->x)*(xv->x)*(xv->x)/4.;
xv->dE = fabs(xv->E-xv->E0)/xv->E0;
xv->t += (pars->dt);
}
void algPuntoDiMezzo(point *xv, parametri *pars, double *xn1, double *vn1) {
printf("%.14g\t%.14g\t%.14g\t%.14g\t%.14g\n", xv->t, xv->x, xv->v, xv->E, xv->dE);
*vn1 = xv->v + f(xv->x, xv->v, pars)*(pars->dt);
xv->x = xv->x + (0.5*(xv->v + (*vn1)))*(pars->dt);
xv->E = 0.5*(xv->v)*(xv->v) - (xv->x)*(xv->x)/8. + (xv->x)*(xv->x)*(xv->x)*(xv->x)/4.;
xv->dE = fabs(xv->E-xv->E0)/xv->E0;
xv->t += (pars->dt);
xv->v = *vn1;
}
void algVerlet(point *xv, parametri *pars, double *xn1, double *vn1) {
printf("%.14g\t%.14g\t%.14g\t%.14g\t%.14g\n", xv->t, xv->x, xv->v, xv->E, xv->dE);
*xn1 = xv->x + xv->v*pars->dt + 0.5*(f(xv->x, xv->v, pars))*pars->dt*pars->dt;
*vn1 = xv->v + 0.5*(f(xv->x, xv->v, pars) + f((*xn1), xv->v, pars))*pars->dt;
xv->x = *xn1;
xv->v = *vn1;
xv->E = 0.5*(xv->v)*(xv->v) - (xv->x)*(xv->x)/8. + (xv->x)*(xv->x)*(xv->x)*(xv->x)/4.;
xv->dE = fabs(xv->E-xv->E0)/xv->E0;
xv->t += (pars->dt);
}
void algRungeKutta2(point *xv, parametri *pars, double *xn1, double *vn1) {
printf("%.14g\t%.14g\t%.14g\t%.14g\t%.14g\n", xv->t, xv->x, xv->v, xv->E, xv->dE);
*xn1 = xv->x + xv->v*pars->dt + 0.5*f(xv->x, xv->v, pars)*pars->dt*pars->dt;
*vn1 = xv->v + f(xv->x + 0.5*xv->v*pars->dt, xv->v + 0.5*f(xv->x, xv->v, pars)*pars->dt, pars)*pars->dt;
xv->x = *xn1;
xv->v = *vn1;
xv->E = 0.5*(xv->v)*(xv->v) - (xv->x)*(xv->x)/8. + (xv->x)*(xv->x)*(xv->x)*(xv->x)/4.;
xv->dE = fabs(xv->E-xv->E0)/xv->E0;
xv->t += (pars->dt);
}
----------------编辑----------------
亲爱的吉尔斯， 我向您解释我的程序是做什么的。该程序旨在用不同的算法（algEulero、algEuleroCromer、algPuntoDiMezzo、algVerlet、algRungeKutta2）对微分方程进行数值求解。它工作正常。物理方程是 d^2x/dt^2 = f(x, v, gamma, epsilon)。这个 f() 正是您可以在我的代码中找到的 f() 。我的简单 C 程序的 "big" 问题是他真的很慢：当我选择 5' 练习 (pars.choice == 5) 时，该程序将完全执行：
1) 计算（getBasin()中的两个for）[-1, 1]x[-1, 1]的面积，增量为pars.dx； 2) 在 x 和 y 上的每个增量都将启动算法，该算法以 for 的两个起始条件 (x, y) 对微分方程进行数值求解。当算法达到一个异步t*(pars.t_basin)时，他会把x(t*)和v(t*)的输出写到basin.dat中。 3) (x, y) 将改变并再次转到 (1) 点。
现在，您可以使用以下参数进行测试：0.83、4、0.01、62、1、1、5、5。
我测试了你的代码，但并不比我的顺序代码快（例如超过 11 分钟）。为了改进它：
1) basin.dat 的输出顺序无关紧要，因为我将根据 3' 和 4' 列值绘制 space (x, y) 着色我的点（与Gnuplot 上的图像）。 2) 为什么还要在函数 getBasin() 之前创建线程，而不仅仅是为两个 for() 创建线程？对于每个线程，这不应该重复 getBasin() 吗？
抱歉我的英语和并行编程不好，我正在努力改进它在线阅读教程。
好吧，您的代码中明显的主要问题是并行版本（非常）错误。您在 parallel 区域之外定义了所有变量，但不声明其中任何一个 private（甚至循环索引也不声明）。
此外，由于您将 getBasin() 函数的所有参数都作为指针传递，之后将它们声明为私有变得更加棘手。
然而，对代码的快速分析表明，尽管这些参数作为指针传递，但您实际上并不关心它们在例程退出时的值。此外，看起来您尝试并行化的循环没有数据依赖性（尽管我没有对此进行全面检查）。我发现的唯一明确的依赖关系是关于你打开的输出文件，在保持顺序的同时并行更新会很棘手...
所以为了修复你的代码的并行化，这就是我所做的：
- 通过值而不是通过引用将参数传递给
getBasin()。这将生成一些额外的数据副本，但由于这些数据需要声明private，因此无论如何都需要一份副本。 - 在
parallel区域内包含对getBasin()的调用。事实上，在我看来，这样做比在函数本身内部处理数据私有化和 IO 问题更简单。 - 为每个线程打开一个输出文件的私有副本，名为 "basinXX.dat"，"XX" 是当前线程的 ID，左填充 0。这些文件将包含共享当前线程对应的全局输出。希望这适合您。否则，处理打印顺序可能会有点棘手。
- 在函数内使用一个孤立的
omp for指令来并行化循环。
有了这个，代码应该（希望）扩展得更好。但是，由于您没有指明用于计算的输入参数，因此我无法对其进行测试，无论是正确性还是性能。所以它可能会惨败...
无论如何，这是我想出的：
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <omp.h>
#define BLD "\x1B[1m"
#define RED "\x1B[31m"
#define GRN "\x1B[32m"
#define RST "\x1B[0m"
struct point{
double x;
double v;
double t;
double E0;
double E;
double dE;
} typedef point;
struct parametri {
double gamma;
double epsilon;
double dt;
double t_star;
double t_basin;
double t_max;
double x0;
double v0;
int choice;
int alg[1];
double dx;
} typedef parametri;
// Prototipi delle funzioni
void Setup();
void getParams(parametri *pars);
void getError(point *xv, parametri *pars, int *i, int k, int *n_passi, double *xn1, double *vn1);
void getBasin(point xv, parametri pars, int h, int n_passi, double xn1, double vn1);
double f(double x, double v, parametri *pars);
void algEulero(point *xv, parametri *pars, double *xn1, double *vn1);
void algEuleroCromer(point *xv, parametri *pars, double *xn1, double *vn1);
void algPuntoDiMezzo(point *xv, parametri *pars, double *xn1, double *vn1);
void algVerlet(point *xv, parametri *pars, double *xn1, double *vn1);
void algRungeKutta2(point *xv, parametri *pars, double *xn1, double *vn1);
int main(void) {
// Inizializzo il display per l'utente. Maggiori informazioni vedere funzione Setup();
Setup();
// Dichiaro e recupero i parametri richiesti per il corretto funzionamento del programma;
parametri pars;
getParams(&pars);
// Dichiaro e ricavo i parametri essenziali, annesse variabili di supporto;
int i, n_passi = pars.t_max/pars.dt;
double dt0, xn1, vn1, t_star = 5.;
point xv;
// Imposto i parametri iniziali;
xv.x = pars.x0;
xv.v = pars.v0;
xv.E0 = 0.5*(xv.v)*(xv.v) - (xv.x)*(xv.x)/8. + (xv.x)*(xv.x)*(xv.x)*(xv.x)/4.;
xv.E = xv.E0;
xv.dE = 0;
xv.t = 0;
pars.t_star = 5;
pars.t_basin = 60;
dt0 = pars.dt;
// Formato dell'output;
printf ("t\tx(t)\tv(t)\tE(t)\tdE\n");
if ((pars.choice == 1) || (pars.choice == 3) || (pars.choice == 4)) { // L'utente ha deciso di affrontare il primo/terzo/quarto esercizio;
// Avverto l'utente che il processo sta per iniziare. Può risultare inutile su tempi brevi, ma efficace su tempi molto lunghi;
fprintf(stderr, "\nAvvio integrazione numerica... ");
if (pars.alg[0] == 1) { // L'utente ha selezionato l'algoritmo di Eulero;
for (i=0; i<n_passi; i++){
algEulero(&xv, &pars, &xn1, &vn1);
}
} else if (pars.alg[0] == 2) { // L'utente ha selezionato l'algoritmo di EuleroCromer;
for (i=0; i<n_passi; i++){
algEuleroCromer(&xv, &pars, &xn1, &vn1);
}
} else if (pars.alg[0] == 3) { // L'utente ha selezionato l'algoritmo di PuntoDiMezzo;
for (i=0; i<n_passi; i++){
algPuntoDiMezzo(&xv, &pars, &xn1, &vn1);
}
} else if (pars.alg[0] == 4) { // L'utente ha selezionato l'algoritmo di Verlet;
for (i=0; i<n_passi; i++){
algVerlet(&xv, &pars, &xn1, &vn1);
}
} else if (pars.alg[0] == 5) { // L'utente ha selezionato l'algoritmo di RungeKutta di ordine 2;
for (i=0; i<n_passi; i++) {
algRungeKutta2(&xv, &pars, &xn1, &vn1);
}
}
// Algoritmo terminato;
fprintf(stderr, "[%s%sDONE%s]\n\n", BLD, GRN, RST);
} else if (pars.choice == 2) { // L'utente ha deciso di affrontare il secondo esercizio;
// Seleziono il secondo algoritmo da confrontare
do {
fprintf(stderr, "> Selezionare il secondo algoritmo:\n");
fprintf(stderr, "\t>> [1] Eulero\n");
fprintf(stderr, "\t>> [2] EuleroCromer - RungeKutta (1 ordine)\n");
fprintf(stderr, "\t>> [3] PuntoDiMezzo\n");
fprintf(stderr, "\t>> [4] Verlet\n");
fprintf(stderr, "\t>> [5] RungeKutta (2 ordine)\n");
fprintf(stderr, "\t>> ");
scanf("%d", &pars.alg[1]);
} while (( pars.alg[1] <= 0 ));
// Avverto l'utente che il processo sta per iniziare. Può risultare inutile su tempi brevi, ma efficace su tempi molto lunghi;
fprintf(stderr, "\nAvvio calcolo errori... ");
// Eseguo lo studio degli errori d'integrazione mediante il primo e secondo algoritmo scelto, rispettivamente nei file error_1.dat, error_2.dat;
getError(&xv, &pars, &i, 0, &n_passi, &xn1, &vn1);
// Resetto le variabili e richiamo l'algoritmo per calcolare gli errori;
pars.dt = dt0;
n_passi = pars.t_max/pars.dt;
xv.x = pars.x0;
xv.v = pars.v0;
xv.E = xv.E0;
xv.dE = 0;
xv.t = 0;
getError(&xv, &pars, &i, 1, &n_passi, &xn1, &vn1);
// Processo terminato;
fprintf(stderr, "\n\nStato: [%s%sDONE%s]\n\n", BLD, GRN, RST);
} else if (pars.choice == 5) { // L'utente ha deciso di affrontare il quinto esercizio;
// Avverto l'utente che il processo sta per iniziare. Può risultare inutile su tempi brevi, ma efficace su tempi molto lunghi;
fprintf(stderr, "\nAvvio calcolo griglia... ");
#pragma omp parallel num_threads( 4 )
getBasin(xv, pars, i, n_passi, xn1, vn1);
// Processo terminato;
fprintf(stderr, "[%s%sDONE%s]\n\n", BLD, GRN, RST);
} else { // L'utente non ha selezionato un esercizio valido;
fprintf(stderr, "\n[%s%sFAILED%s] Esercizio non disponibile.\n", BLD, RED, RST);
exit(EXIT_FAILURE);
}
return 0;
}
void Setup() {
fprintf(stderr, "\nAnalisi numerica di un'equazione differenziale\n");
fprintf(stderr, "==============================================\n\n");
}
void getParams(parametri *pars) {
do {
fprintf(stderr, "> Inserire un valore per gamma : ");
scanf("%lf", &pars->gamma);
} while (pars->gamma < 0);
do {
fprintf(stderr, "> Inserire un valore per epsilon : ");
scanf("%lf", &pars->epsilon);
} while (pars->epsilon < 0);
do {
fprintf(stderr, "> Inserire un valore per dt : ");
scanf("%lf", &pars->dt);
} while (pars->dt < 0);
do {
fprintf(stderr, "> Inserire un valore per tmax t.c. :\n");
fprintf(stderr, " >> (tmax > 5) per I eserc.\n");
fprintf(stderr, " >> (tmax > 60) per V eserc. : ");
scanf("%lf", &pars->t_max);
} while (pars->t_max < 0);
do {
fprintf(stderr, "> Inserire un valore per x(0) : ");
scanf("%lf", &pars->x0);
} while (pars->x0 < -1);
do {
fprintf(stderr, "> Inserire un valore per v(0) : ");
scanf("%lf", &pars->v0);
} while (pars->v0 < -1);
do {
fprintf(stderr, "> Selezionare l'esercizio richiesto :\n");
fprintf(stderr, "\t>> [1] Esercizio 1\n");
fprintf(stderr, "\t>> [2] Esercizio 2\n");
fprintf(stderr, "\t>> [3] Esercizio 3\n");
fprintf(stderr, "\t>> [4] Esercizio 4\n");
fprintf(stderr, "\t>> [5] Esercizio 5\n\n");
fprintf(stderr, "\t>> ");
scanf("%d", &pars->choice);
} while (( pars->choice <= 0 ));
do {
fprintf(stderr, "> Selezionare l'algoritmo voluto :\n");
fprintf(stderr, "\t>> [1] Eulero\n");
fprintf(stderr, "\t>> [2] EuleroCromer - RungeKutta (1 ordine)\n");
fprintf(stderr, "\t>> [3] PuntoDiMezzo\n");
fprintf(stderr, "\t>> [4] Verlet\n");
fprintf(stderr, "\t>> [5] RungeKutta (2 ordine)\n\n");
fprintf(stderr, "\t>> ");
scanf("%d", &pars->alg[0]);
} while (( pars->alg[0] <= 0 ));
}
void getError(point *xv, parametri *pars, int *i, int k, int *n_passi, double *xn1, double *vn1) {
void (*algF)(point *, parametri *, double *, double *);
int j, n = 0;
FILE *fp;
// Questo if controlla se il codice sta eseguendo lo studio degli errori per il primo o secondo algoritmo (pars->alg[k]);
if (k == 0) {
fp = fopen("error_1.dat", "w+");
} else {
fp = fopen("error_2.dat", "w+");
}
// Assegno il puntatore corretto a seconda dell'algoritmo scelto;
if (pars->alg[k] == 1) {
algF = algEulero;
} else if (pars->alg[k] == 2) {
algF = algEuleroCromer;
} else if (pars->alg[k] == 3) {
algF = algPuntoDiMezzo;
} else if (pars->alg[k] == 4) {
algF = algVerlet;
} else if (pars->alg[k] == 5) {
algF = algRungeKutta2;
} else {
fprintf(stderr, "\n[%s%sFAILED%s] E' stato selezionato un algoritmo non valido.\n", BLD, RED, RST);
exit(EXIT_FAILURE);
}
// Informo l'utente che il processo sta iniziando;
fprintf(stderr, "\n>> Avvio %d algoritmo... ", k+1);
// Formattazione dell'output del file contenente gli errori;
fprintf(fp, "dt\tE(t*)\tE(0)\tdE/E(0)\t# passi\ti\tt\n");
for (j=0; j<8; j++) {
for ((*i)=0; (*i)<(*n_passi); (*i)++){
(*algF)(xv, pars, xn1, vn1);
if (((pars->t_star - pars->dt/2. <= xv->t) && (xv->t >= pars->t_star + pars->dt/2.)) && (n == 0)) {
fprintf(fp, "%+.14lf\t%+.14lf\t%+.14lf\t%+.14lf\t%+.14d\t%+.14d\t%+.14lf\n", pars->dt, xv->E, xv->E0, (xv->E - xv->E0)/xv->E0, (*n_passi), (*i), xv->t);
n = 1;
}
}
// Resetto le variabili per rilanciare l'algoritmo con dt/2^j
n = 0;
xv->t = 0;
xv->x = pars->x0;
xv->v = pars->v0;
pars->dt = pars->dt/2.;
(*n_passi) = pars->t_max/pars->dt;
(*xn1) = 0;
(*vn1) = 0;
xv->E = xv->E0;
}
fclose(fp);
fprintf(stderr, "[%s%sDONE%s]", BLD, GRN, RST);
}
void getBasin(point xv, parametri pars, int h, int n_passi, double xn1, double vn1) {
// Dichiaro e setto i parametri che delimitano la griglia;
point xv1;
pars.dx = 0.01;
xv1.x = -1;
xv1.v = -1;
// Dichiaro le variabili necessarie per il bacino d'attrazione;
int i, j, i_max = 2./pars.dx, j_max = 2./pars.dx;
void (*algF)(point *, parametri *, double *, double *);
char fname[13];
sprintf( fname, "bassin%02d.dat", omp_get_thread_num() );
FILE *fp = fopen( fname, "w+");
// Assegno il puntatore corretto a seconda dell'algoritmo scelto;
if (pars.alg[0] == 1) {
algF = algEulero;
} else if (pars.alg[0] == 2) {
algF = algEuleroCromer;
} else if (pars.alg[0] == 3) {
algF = algPuntoDiMezzo;
} else if (pars.alg[0] == 4) {
algF = algVerlet;
} else if (pars.alg[0] == 5) {
algF = algRungeKutta2;
} else {
fprintf(stderr, "\n[%s%sFAILED%s] E' stato selezionato un algoritmo non valido.\n", BLD, RED, RST);
exit(EXIT_FAILURE);
}
// Formattazione output file basin.dat;
fprintf(fp, "x(0)\tx'(0)\tx(t*)\tv(t*)\n");
#pragma omp for
// Eseguo il for della griglia sull'asse x';
for (j=0; j<=j_max; j++) {
// Eseguo il for della griglia sull'asse x;
for (i=0; i<=i_max; i++) {
fprintf(fp, "%lf\t%lf\t", xv1.x, xv1.v);
xv.x = xv1.x;
xv.v = xv1.v;
// Eseguo l'integrazione numerica
for (h=0; h<n_passi; h++) {
(*algF)(&xv, &pars, &xn1, &vn1);
// Entro in t = t*, stampo v(t*) ed esco dal ciclo;
if ((pars.t_basin - pars.dt/2. <= xv.t) && (xv.t >= pars.t_basin + pars.dt/2.)) {
fprintf(fp, "%lf\t%lf\n", xv.x, xv.v);
break;
}
}
xv1.x += pars.dx;
xv.t = 0;
xn1 = 0;
vn1 = 0;
}
// Resetto la x e incremento la x';
xv1.x = -1;
xv1.v += pars.dx;
}
}
double f(double x, double v, parametri *pars) {
return 0.25*x - x*x*x + (pars->gamma - pars->epsilon*(x*x))*v;
}
void algEulero(point *xv, parametri *pars, double *xn1, double *vn1){
//printf("%.14g\t%.14g\t%.14g\t%.14g\t%.14g\n", xv->t, xv->x, xv->v, xv->E, xv->dE);
*xn1 = xv->x + (xv->v)*(pars->dt);
*vn1 = xv->v + f(xv->x, xv->v, pars)*(pars->dt);
xv->E = 0.5*(xv->v)*(xv->v) - (xv->x)*(xv->x)/8. + (xv->x)*(xv->x)*(xv->x)*(xv->x)/4.;
xv->dE = fabs(xv->E-xv->E0)/xv->E0;
xv->t += (pars->dt);
xv->x = *xn1;
xv->v = *vn1;
}
void algEuleroCromer(point *xv, parametri *pars, double *xn1, double *vn1){
//printf("%.14g\t%.14g\t%.14g\t%.14g\t%.14g\n", xv->t, xv->x, xv->v, xv->E, xv->dE);
xv->v = xv->v + f(xv->x, xv->v, pars)*(pars->dt);
xv->x = xv->x + (xv->v)*(pars->dt);
xv->E = 0.5*(xv->v)*(xv->v) - (xv->x)*(xv->x)/8. + (xv->x)*(xv->x)*(xv->x)*(xv->x)/4.;
xv->dE = fabs(xv->E-xv->E0)/xv->E0;
xv->t += (pars->dt);
}
void algPuntoDiMezzo(point *xv, parametri *pars, double *xn1, double *vn1) {
//printf("%.14g\t%.14g\t%.14g\t%.14g\t%.14g\n", xv->t, xv->x, xv->v, xv->E, xv->dE);
*vn1 = xv->v + f(xv->x, xv->v, pars)*(pars->dt);
xv->x = xv->x + (0.5*(xv->v + (*vn1)))*(pars->dt);
xv->E = 0.5*(xv->v)*(xv->v) - (xv->x)*(xv->x)/8. + (xv->x)*(xv->x)*(xv->x)*(xv->x)/4.;
xv->dE = fabs(xv->E-xv->E0)/xv->E0;
xv->t += (pars->dt);
xv->v = *vn1;
}
void algVerlet(point *xv, parametri *pars, double *xn1, double *vn1) {
//printf("%.14g\t%.14g\t%.14g\t%.14g\t%.14g\n", xv->t, xv->x, xv->v, xv->E, xv->dE);
*xn1 = xv->x + xv->v*pars->dt + 0.5*(f(xv->x, xv->v, pars))*pars->dt*pars->dt;
*vn1 = xv->v + 0.5*(f(xv->x, xv->v, pars) + f((*xn1), xv->v, pars))*pars->dt;
xv->x = *xn1;
xv->v = *vn1;
xv->E = 0.5*(xv->v)*(xv->v) - (xv->x)*(xv->x)/8. + (xv->x)*(xv->x)*(xv->x)*(xv->x)/4.;
xv->dE = fabs(xv->E-xv->E0)/xv->E0;
xv->t += (pars->dt);
}
void algRungeKutta2(point *xv, parametri *pars, double *xn1, double *vn1) {
//printf("%.14g\t%.14g\t%.14g\t%.14g\t%.14g\n", xv->t, xv->x, xv->v, xv->E, xv->dE);
*xn1 = xv->x + xv->v*pars->dt + 0.5*f(xv->x, xv->v, pars)*pars->dt*pars->dt;
*vn1 = xv->v + f(xv->x + 0.5*xv->v*pars->dt, xv->v + 0.5*f(xv->x, xv->v, pars)*pars->dt, pars)*pars->dt;
xv->x = *xn1;
xv->v = *vn1;
xv->E = 0.5*(xv->v)*(xv->v) - (xv->x)*(xv->x)/8. + (xv->x)*(xv->x)*(xv->x)*(xv->x)/4.;
xv->dE = fabs(xv->E-xv->E0)/xv->E0;
xv->t += (pars->dt);
}
编辑
使用您提供的新输入参数，我能够测试代码，结果是：
- 我给你的代码中有一个小错误（用于打开输出文件，使用
"fname"而不是fname） - 您的代码大部分（全部）时间都花在
printf()。
修复这两个后，编译如下：
gcc -O3 -fopenmp -mtune=native -march=native -w bassin.c
我 运行 它在我的双核笔记本电脑上 2.12s
代码已更新。请自行尝试。