how to do 3d sum using openmp

Viewed 49

I am a freshman in openmp. I have some trouble in a 3d sum, and I don't know how to improve my code. Here's the code I want to improve in openmp. My aim is to speed up the calculation of this 3d sum. What should I add in my code according to the rules of openmp? I add #pragma omp parallel for reduction(+:integral) in my code. But an error happens which says the initialization of 'for' is not correct. This is the information of this error:enter image description here I am a chinese, so the language of my IDE is chinese. I use Visual Studio 2019.

 #include<omp.h>
#include<stdio.h>
#include<math.h>
int main()
{
    double a = 0.3291;
    double d_title = 2.414;
    double b = 3.8037;
    double c = 4086;
    double nu_start = 0;
    double mu_start = 0;
    double z_start = 0;
    double step_nu = 2 * 3.1415926 / 100;
    double step_mu = 3.1415926 / 100;
    double step_z = 0;
    double nu = 0;
    double mu = 0;
    double z = 0;
    double integral=0;
    double d_uv = 0;
    int i = 0;
    int j = 0;
    int k = 0;
#pragma omp parallel for default(none) shared(a, d_title, b, c, nu_start, mu_start, z_start, step_nu, step_mu) private( j,k,mu, nu, step_z, z, d_uv) reduction(+:integral)
    for (i = 0; i < 100; i++)
    {
        mu = mu_start + (i + 1) * step_mu;
        for (j = 0; j < 100; j++)
        {
            
            nu = nu_start + (j + 1) * step_nu;
            for (k = 0; k < 500; k++)
            {
                d_uv = (sin(mu) * sin(mu) * cos(nu) * cos(nu) + sin(mu) * sin(mu) * (a * sin(nu) - d_title * cos(nu)) * (a * sin(nu) - d_title * cos(nu)) + b * b * cos(mu) * cos(mu)) / (c * c);
                step_z = 20 / (d_uv * 500);
                z = z_start + (k + 1) * step_z;
                
                integral = integral + sin(mu) * (1 - 3 * sin(mu) * sin(mu) * cos(nu) * cos(nu)) * exp(-d_uv * z) * log(1 + z * z) * step_z * step_mu * step_nu;
            }

        }

    }

    double out = 0;
    out = integral / (c * c);
    return 0;
}
1 Answers

Solutions (UPDATE: It is an answer to the original question:)

  1. To do the least typing you just have to add the following line before for(int i=..)

#pragma omp parallel for private( mu, nu, step_z, z, d_uv) reduction(+:integral)

Here you define which variables have to be private to avoid data race. Note that variables are shared by default, so variable integral also shared, but all threads update its value, which is a data race. To avoid it, you have 2 possibilities: use atomic operation, or a much better option is to use use reduction (add reduction(+:integral) clause).

  1. As you mentioned that you are beginner in OpenMP it is recommended to use default(none) clause in the #pragma omp parallel for directive, so you have to explicitly define sharing attributes. If you forget a variable you will get an error, so you have to consider all variables involved in your parallel region and can think about possible data races:

#pragma omp parallel for default(none) shared(a, d_title, b, c, nu_start, mu_start, z_start, step_nu, step_mu) private( mu, nu, step_z, z, d_uv) reduction(+:integral)

  1. Generally, it is recommended to define your variables in their minimum required scope, so variables defined inside the for loop to parallelize will be private. In this case you just have to add #pragma omp parallel for reduction(+:integral) before your outermost for loop, so your code will be:
    #pragma omp parallel for reduction(+:integral)
    for (int i = 0; i < 100; i++)
    {
        double mu = mu_start + (i + 1) * step_mu;
        for (int j = 0; j < 100; j++)
        {
            //int id = omp_get_thread_num();
            double nu = nu_start + (j + 1) * step_nu;
            for (int k = 0; k < 500; k++)
            {
                double d_uv = (sin(mu) * sin(mu) * cos(nu) * cos(nu) + sin(mu) * sin(mu) * (a * sin(nu) - d_title * cos(nu)) * (a * sin(nu) - d_title * cos(nu)) + b * b * cos(mu) * cos(mu)) / (c * c);
                double step_z = 20 / (d_uv * 500);
                double z = z_start + (k + 1) * step_z;
                //int id = omp_get_thread_num();
                integral = integral + sin(mu) * (1 - 3 * sin(mu) * sin(mu) * cos(nu) * cos(nu)) * exp(-d_uv * z) * log(1 + z * z) * step_z * step_mu * step_nu;
            }

        }

    } 

Runtimes: 44 ms (1 thread) and 11 ms (4 threads) on my computer (g++ -O3 -mavx2 -fopenmp).

Related