2015-11-15 34 views
0

我试图改进我的C源代码以并行执行。我有一个四核CPU,所以我认为4是很多线程(一个用于CPU)来运行我的程序,而不是像顺序代码那样进行优化。OpenMP #pragma omp parallel for slow

但它不工作。我的代码没有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程序的“大”问题是他非常慢:当我选择5'练习(pars.choice == 5)时,程序将完成:

1)计算在getBasin())中,一个面积为[-1,1] x [-1,1]的pars.dx增量; 2)在x和y上的每个增量将启动算法,该算法用for的两个起始条件(x,y)数值地求解微分方程。当算法达到asynthotic 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)盆地输出的顺序是微不足道的,因为我将绘制依赖于3'和4'列值的图像(x,y)在Gnuplot上)。 2)为什么在函数getBasin()之前创建线程,而不是仅为两个for()?创建线程。这不应该重复,每个线程,getBasin()?

对不起,我的英语不好,并行编程,我试图改进它在线阅读教程。

回答

4

那么,你的代码的主要问题是并行版本是(非常)错误的。您可以在parallel区域之外定义所有变量,但不要声明它们中的任何一个private(甚至不包括循环索引)。

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

但是,对代码的快速分析表明,尽管这些参数是作为指针传递的,但实际上并不在意它们在退出例程时的值。此外,它似乎没有数据依赖于你试图平行的循环(虽然我没有做一个全面的检查)。我发现的唯一明确的依赖关系是关于您打开的输出文件,这将在并行更新中保持顺序顺序并不容易...

因此,为了解决你的代码的并行化,这是我做的:

  1. 通过值而不是通过引用传递的参数getBasin()。这会生成一些额外的数据副本,但由于这些数据需要声明为private,因此无论如何都需要副本。
  2. parallel区域内呼叫getBasin()。事实上,在我看来,这样做比较简单,而不是在函数本身内部处理数据私有化和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 ("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); 
} 

编辑

有了新的输入参数你给的,我是能够测试的代码,事实证明:

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

固定这两个和编译这样后:

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

2.12s

跑在我的双核笔记本电脑的代码更新。请自行尝试。

+0

太棒了!谢谢 ! 'printf()'非常糟糕。代码是完美的,但最后一个问题是,bassin00,bassin01,bassin02,bassin03产生的是完全相同的数据。如何告诉omp在不同索引的不同线程中分割for()'? –

+0

'bassin00.dat'比我的输出中的其他人略有些,但是,其他三个完全相同。也就是说,代码确实与线程间分布的'j'循环迭代并行。简单地说,我没有在代码中看到对'j'的任何依赖。所以要么这里有问题,要么代码本身是迭代的,不能(很容易)并行化。无论如何,我的笔记本电脑上的顺序版本现在是6.07s。 – Gilles

+0

在我的桌面上,顺序也是〜20秒,没关系。但我不知道为什么第一个for()(j一)不是从'0'到'j_max',而是从'0'到几分之一'j_max',每个线程都这样做。尝试绘制bassin01-4的输出,第一列的函数为2,颜色由3列确定:您将看到三个输出是相同的。问题是脚本正在执行顺序代码中的对应工作,但在并行版本中并没有这样做。 –