文档建议反馈控制台
首页
学习
活动
专区
工具
TVP
最新优惠活动
文章/答案/技术大牛
发布
社区首页 >问答首页 >OpenMP并行化GSL型常微分方程的计算

问OpenMP并行化GSL型常微分方程的计算
EN

Stack Overflow用户
提问于 2019-06-09 16:05:20
回答 2查看 155关注 0票数 0

我正在尝试并行化我的代码，但我得到了错误。我需要计算一个柯西问题(已经完成了)，但是我需要使用OpenMP库对其进行并行化。

我试着用OpenMP写了一些代码，但是它不能工作。

我已经创建了一个结构来收集结果。

代码语言：javascript
复制
struct Dots {
    double par;
    double x;
    double y;
};

这是我的带参数的目标函数。

代码语言：javascript
复制
int ode_func (double x, const double y[], double f[], void *params)
{

    double mu = *(int *)params;
    f[0] = x + 2 * y[0] / (1 + mu * mu);
    return GSL_SUCCESS;
}

这是主函数。我目前没有找到一种方法来创建结构数组的数组，但这不是主要问题。

代码语言：javascript
复制
void calc_cauchy_problem(struct Dots ArrayOfDots[], double x_start, double x_end, double y_start,
        int count) {

    int dim = 1;
    double x = x_start;
    double y[1] = {y_start};
    int mu = 5;
    int param = 0;
    gsl_odeiv2_system sys = {ode_func, NULL, dim, &param};
    gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new (&sys,
                                                               gsl_odeiv2_step_rkf45, 1e-6, 1e-6, 0.0);
    int status = 0;
#pragma omp parallel for shared(ArrayOfDots) private(sys, param, d, status)
    for (int param = 1; param < mu; param++) {

        gsl_odeiv2_system sys = {ode_func, NULL, dim, &param};
        gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new (&sys,
                                                               gsl_odeiv2_step_rkf45, 1e-6, 1e-6, 0.0);
        for (int i = 1; i <= count; i++)
        {
            double xi = x_start + i * (x_end - x_start) / count;

            int status = gsl_odeiv2_driver_apply(d, &x, xi, y);

            if (status != GSL_SUCCESS)
            {
                printf ("error, return value=%d\n", status);
                break;
            }
           // ArrayOfDots[i].par = mu;
           // ArrayOfDots[i].x = xi;
           // ArrayOfDots[i].y = y[0];
        }
        gsl_odeiv2_driver_free (d);
    }

}

main

代码语言：javascript
复制
int main() {
    double x_start = 0;
    double x_end = 10;
    double y_start = 0;
    int count = 10;
    struct Dots ArrayOfDots[count];
    calc_cauchy_problem(ArrayOfDots, x_start, x_end, y_start, count);
    return 0;
}

它用这个gcc main.c -o main -fopenmp -lgsl -std=gnu11编译成功，但当我启动它时，我得到了错误

代码语言：javascript
复制
gsl: driver.c:354: ERROR: integration limits and/or step direction not consistent
Default GSL error handler invoked.

我认为这是#pragma omp parallel for shared(ArrayOfDots) private(sys, param, d, status)的主要问题，但我不知道如何以其他方式重写它。感谢您的回复。

更新：

使用Kaveh帮助我的代码开始部分工作。这意味着我的for循环的一半开始工作了。UPD UPD:在另一次调查之后，我得到了以下代码:它是编译并运行的，但我得到的是Process finished with exit code 4和printf("Elapsed time = %f\n", omp_get_wtime() - start_time);不打印任何东西。

代码语言：javascript
复制
struct Dots {
    double par;
    double x;
    double y;
};

int ode_func (double x, const double y[], double f[], void *params)
{

    double mu = *(int *)params;
    f[0] = (x + 2 * y[0]) / (1 + mu * mu);
    return GSL_SUCCESS;
}
void calc_cauchy_problem(double x_start, double x_end, double y_start,
                         int count, int param1, int param2) {
    int dim = 1;
    double x = x_start;
    double y[1] = {y_start};
    int param = param1;
    int j = 0;
    int status = 0;
    char filename[10];

#pragma omp parallel for private(param, status, x, y)
    for (param = param1; param <= param2; param++) {
        struct Dots ArrayOfDots[count];
        gsl_odeiv2_system sys = {ode_func, NULL, dim, &param};
        gsl_odeiv2_driver * d =
                gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rkf45, 1e-6, 1e-6, 0.0);
        for (int i = 1; i <= count; i++) {
            double xi = x_start + i * (x_end - x_start) / count;

            int status = gsl_odeiv2_driver_apply(d, &x, xi, y);
            if (status != GSL_SUCCESS)
            {
                printf ("error, return value=%d\n", status);
                break;
            }

            ArrayOfDots[i].par = param;
            ArrayOfDots[i].x = xi;
            ArrayOfDots[i].y = y[0];
        }
        gsl_odeiv2_driver_free (d);

    }
}
int main() {
    double start_time = omp_get_wtime();
    double x_start = 0;
    double x_end = 10;
    double y_start = 0;
    const int count = 500;
    int param1 = 1;
    int param2 = 10;
    calc_cauchy_problem(x_start, x_end, y_start, count, param1, param2);
    printf("Elapsed time = %f\n", omp_get_wtime() - start_time);
    return 0;
}
openmp
gsl
c
EN

回答 2

Stack Overflow用户

回答已采纳

发布于 2019-06-09 20:26:52

将x添加到私有循环vars：private(sys, param, d, status, x)。如果您仍然遇到问题，请联系我。

代码语言：javascript
复制
void calc_cauchy_problem(double x_start, double x_end, double y_start,
                         int count, int param1, int param2) {

  int dim = 1;
  double x = x_start;
  double y[1] = {y_start};
  int param = param1;
  int j = 0;
  int status = 0;
  char filename[10];

#pragma omp parallel for private(param, status, x, y)
  for (param = param1; param <= param2; param++) {
    struct Dots ArrayOfDots[count];
    gsl_odeiv2_system sys = {ode_func, NULL, dim, &param};
    gsl_odeiv2_driver * d =
      gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rkf45, 1e-6, 1e-6, 0.0);
    for (int i = 1; i <= count; i++) {
      double xi = x_start + i * (x_end - x_start) / count;

      int status = gsl_odeiv2_driver_apply(d, &x, xi, y);
      if (status != GSL_SUCCESS)
        {
          printf ("error, return value=%d\n", status);
          break;
        }
      ArrayOfDots[i].par = param;
      ArrayOfDots[i].x = xi;
      ArrayOfDots[i].y = y[0];
    }
    //write_data_to_file(param, count, ArrayOfDots);                                                                                        
    for (int i = 0; i < count; ++i) {
      printf ("%d: %f, %f, %f\n", omp_get_thread_num(),
              ArrayOfDots[i].par, ArrayOfDots[i].x, ArrayOfDots[i].y);
    }
    gsl_odeiv2_driver_free (d);
  }
}
票数 1
EN

Stack Overflow用户

发布于 2019-10-10 17:02:58

看起来这个版本运行得很好。我认为问题出在这个结构体Dots ArrayOfDots[count];上，当我试图将值推入这个结构体时。

代码语言：javascript
复制
      ArrayOfDots[i].par = param;
      ArrayOfDots[i].x = xi;
      ArrayOfDots[i].y = y[0];

下面是完整的代码。

代码语言：javascript
复制
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>
// GSL lib includes
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>


int ode_func (double x, const double y[], double f[], void *params)
{

    double mu = *(int *)params;
    f[0] = (x + 2 * y[0]) / (1 + mu * mu);
    return GSL_SUCCESS;
}

void calc_cauchy_problem(double x_start, double x_end, double y_start,
                         int count, int param1, int param2) {

#pragma omp parallel for
    for(int param = param1; param < param2; param++) {
        gsl_odeiv2_system sys = {ode_func, NULL, 1, &param};

        gsl_odeiv2_driver * d =
                gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk8pd,
                                               1e-6, 1e-6, 0.0);
        int i;
        double x = x_start, x1 = x_end;
        double y[1] = { y_start };

        for (i = 1; i <= count; i++)
        {
            double xi = i * x1 / count;
            int status = gsl_odeiv2_driver_apply (d, &x, xi, y);

            if (status != GSL_SUCCESS)
            {
                printf ("error, return value=%d\n", status);
                break;
            }

//            printf ("%d %d %.5e %.5e\n", omp_get_thread_num(), param, x, y[0]);
        }

        gsl_odeiv2_driver_free (d);
        }
    }

int main() {
    double start_time = omp_get_wtime();
    double x_start = 0;
    double x_end = 10;
    double y_start = 0;
    const int count = 100000;
    int param1 = 1;
    int param2 = 20;
    calc_cauchy_problem(x_start, x_end, y_start, count, param1, param2);
    printf("Elapsed time = %f\n", omp_get_wtime() - start_time);
    return 0;
}

真的要感谢卡维·瓦赫迪普尔。

票数 0
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：

https://stackoverflow.com/questions/56512937

复制
相关文章

相似问题

领券

腾讯云开发者

扫码关注腾讯云开发者

扫码关注腾讯云开发者

领取腾讯云代金券

热门产品

热门推荐

更多推荐

Copyright © 2013 - 2024 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有 

深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059 深公网安备号 44030502008569

腾讯云计算（北京）有限责任公司 京ICP证150476号 |  京ICP备11018762号 | 京公网安备号11010802020287

问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档