的OpenMP#pragma omp并行速度较慢

Question

我正在努力改进我的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 ("ttx(t)tv(t)tE(t)tdEn");
    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]nn", 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] Euleron");
            fprintf(stderr, "t>> [2] EuleroCromer - RungeKutta (1 ordine)n");
            fprintf(stderr, "t>> [3] PuntoDiMezzon");
            fprintf(stderr, "t>> [4] Verletn");
            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, "nnStato: [%s%sDONE%s]nn", 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]nn", 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 differenzialen");
    fprintf(stderr, "==============================================nn");
}
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 1n");
        fprintf(stderr, "t>> [2] Esercizio 2n");
        fprintf(stderr, "t>> [3] Esercizio 3n");
        fprintf(stderr, "t>> [4] Esercizio 4n");
        fprintf(stderr, "t>> [5] Esercizio 5nn");
        fprintf(stderr, "t>> ");
        scanf("%d", &pars->choice);
    } while (( pars->choice <= 0 ));
    do {
        fprintf(stderr, "> Selezionare l'algoritmo voluto    :n");
        fprintf(stderr, "t>> [1] Euleron");
        fprintf(stderr, "t>> [2] EuleroCromer - RungeKutta (1 ordine)n");
        fprintf(stderr, "t>> [3] PuntoDiMezzon");
        fprintf(stderr, "t>> [4] Verletn");
        fprintf(stderr, "t>> [5] RungeKutta (2 ordine)nn");
        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, "dttE(t*)tE(0)tdE/E(0)t# passitittn");
    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, "%+.14lft%+.14lft%+.14lft%+.14lft%+.14dt%+.14dt%+.14lfn", 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, "%lft%lft", 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, "%lft%lfn", 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("%.14gt%.14gt%.14gt%.14gt%.14gn", 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("%.14gt%.14gt%.14gt%.14gt%.14gn", 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("%.14gt%.14gt%.14gt%.14gt%.14gn", 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("%.14gt%.14gt%.14gt%.14gt%.14gn", 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("%.14gt%.14gt%.14gt%.14gt%.14gn", 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);
}

------------------编辑------------------

尊敬的Gilles：，我向你解释我的程序是做什么的。这个程序旨在用不同的算法（algEulero，algEuler，algPuntoDiMezzo，algVerlet，algRungeKutta2）数值求解微分方程。它运行良好。物理方程为d^2x/dt^2=f（x，v，gamma，epsilon）。这个f（）正是你可以在我的代码中找到的f（）。我的简单C程序的"大"问题是他真的很慢：当我选择5'练习（pars.choice==5）时，程序会准确地执行：

1） 用pars.dx增量计算（getBasin（）中的两个）[-1，1]x[-1，1]的面积；2） 在x和y上的每一个增量都将启动算法，该算法以for的两个起始条件（x，y）数值求解微分方程。当算法达到渐近t*（pars.t_basin）时，他将在basin.dat中写入x（t*）和v（t*）的输出。3） （x，y）将改变并且再次转到点（1）。

现在，您可以使用以下参数进行测试：0.83、4、0.01、62、1、1、5、5。

我测试了你的代码，但并不比我的顺序代码快（例如，超过11分钟）。为了改进它：

1） basin.dat的输出顺序无关紧要，因为我将根据3'和4'列的值（在Gnuplot上有图像）绘制空间（x，y）给我的点着色。2） 为什么还要在函数getBasin（）之前创建线程，而不仅仅是为两个for（）创建线程？对于每个线程，不应该重复getBasin（）吗？

很抱歉我的英语和并行编程不好，我正在努力提高在线阅读教程。

Answer 1

好吧，代码中明显的主要问题是并行版本（非常）错误。您在parallel区域之外定义所有变量，但不声明其中任何一个private（甚至不声明循环索引）。

此外，由于您将getBasin()函数的所有参数都作为指针传递，因此之后将它们声明为私有参数会变得更加棘手。

然而，对代码的快速分析表明，尽管这些参数是作为指针传递的，但在例程退出时，您实际上并不关心它们的值。此外，您尝试并行化的循环似乎没有数据依赖性（尽管我没有对此进行全面的检查）。我发现的唯一明显的依赖关系是关于您打开的输出文件，在保持顺序的同时并行更新会很困难。。。

因此，为了修复代码的并行化，我做了以下操作：

1. 通过值而不是通过引用将参数传递给getBasin()。这将生成一些额外的数据副本，但由于这些数据需要声明为private，因此无论如何都需要一个副本
2. 将对getBasin()的调用封装在parallel区域内。事实上，在我看来，这样做比处理功能内部的数据私有化和IO问题更简单
3. 每个线程打开一个输出文件的私有副本，名为"basinXX.dat"，其中"XX"是当前线程的id，用0填充。这些文件将包含与当前线程对应的全局输出的共享。希望这件衣服适合你。否则，处理打印的顺序可能会有点棘手
4. 在函数内部使用单个孤立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 ("ttx(t)tv(t)tE(t)tdEn");
    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]nn", 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] Euleron");
            fprintf(stderr, "t>> [2] EuleroCromer - RungeKutta (1 ordine)n");
            fprintf(stderr, "t>> [3] PuntoDiMezzon");
            fprintf(stderr, "t>> [4] Verletn");
            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, "nnStato: [%s%sDONE%s]nn", 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]nn", 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 differenzialen");
    fprintf(stderr, "==============================================nn");
}
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 1n");
        fprintf(stderr, "t>> [2] Esercizio 2n");
        fprintf(stderr, "t>> [3] Esercizio 3n");
        fprintf(stderr, "t>> [4] Esercizio 4n");
        fprintf(stderr, "t>> [5] Esercizio 5nn");
        fprintf(stderr, "t>> ");
        scanf("%d", &pars->choice);
    } while (( pars->choice <= 0 ));
    do {
        fprintf(stderr, "> Selezionare l'algoritmo voluto    :n");
        fprintf(stderr, "t>> [1] Euleron");
        fprintf(stderr, "t>> [2] EuleroCromer - RungeKutta (1 ordine)n");
        fprintf(stderr, "t>> [3] PuntoDiMezzon");
        fprintf(stderr, "t>> [4] Verletn");
        fprintf(stderr, "t>> [5] RungeKutta (2 ordine)nn");
        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, "dttE(t*)tE(0)tdE/E(0)t# passitittn");
    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, "%+.14lft%+.14lft%+.14lft%+.14lft%+.14dt%+.14dt%+.14lfn", 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, "%lft%lft", 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, "%lft%lfn", 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("%.14gt%.14gt%.14gt%.14gt%.14gn", 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("%.14gt%.14gt%.14gt%.14gt%.14gn", 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("%.14gt%.14gt%.14gt%.14gt%.14gn", 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("%.14gt%.14gt%.14gt%.14gt%.14gn", 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("%.14gt%.14gt%.14gt%.14gt%.14gn", 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);
}

---

编辑

有了你给出的新输入参数，我可以测试代码，结果是：

1. 我在给您的代码中有一个小错误（用于打开输出文件，使用"fname"而不是fname）
2. 您的代码大部分（全部）时间都在printf()中

在修复了这两个之后，并像这样编译：

gcc -O3 -fopenmp -mtune=native -march=native -w bassin.c

我在2.12s的双核笔记本电脑上运行了它

代码已更新。请你自己试试。