OpenMP求π和实现矩阵乘法

OpenMP积分法求解π

• 使用sum数组保存各线程的局部和
• 将n次计算平均分配给NUM_THREADS个线程，各线程并行计算其分配到的计算任务

#include <omp.h>
#include <stdio.h>
#define NUM_THREADS 4
static long num_steps = 100000;
double step = 0;
int main()
{
    int i =0;
    double pi = 0;
    double sum[NUM_THREADS];
    double start_time,end_time;
    step = 1.0/(double)num_steps;
    omp_set_num_threads(NUM_THREADS);
    start_time = omp_get_wtime();
# pragma omp parallel

    {
        int id;
        double x;
        id = omp_get_thread_num();
        for (i =id,sum[id]=0;i<num_steps;i=i+NUM_THREADS)
        {
            x = (i+0.5)*step;
            sum[id] +=4.0/(1.0+x*x);
        }
    }
    for (i=0,pi=0;i<NUM_THREADS;i++) pi+=sum[i]*step;
    end_time = omp_get_wtime();
    printf("PI = %f\n",pi);
    printf("Running time = %lf\n",end_time-start_time);
    return 0;
}

OpenMP实现矩阵乘法

矩阵相乘是线性代数中最常见的问题之一，它在数值计算中有广泛的应用，在计算机的世界里，矩阵相乘扮演着一个不可或缺的角色。设A和B是2个
矩阵，它们的乘积AB同样是一个 矩阵。A和B乘积矩阵C中元素

实验参考~~

https://blog.csdn.net/qq_34149581/article/details/93882400

矩阵生成函数（利用C++11生成随机数的方法）

void init_matrix() {
    default_random_engine engine;
    uniform_real_distribution<double> u(0.0, 1.0);
    for (int i = 0; i < M; ++i) {
        for (int j = 0; j < N; ++j) {
            matrix_a[i][j] = u(engine);
        }
    }
    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < M; ++j) {
            matrix_b[i][j] = u(engine);
        }
    }
}

串行算法描述

若以上面的定义来计算A和B的乘积矩阵C，则每计算C的一个元素 ，需要作n次乘法运算和n－1次加法运算。因此，算出矩阵C的 个元素所需的计算时间为 。从程序运行的效率上来分析，如果n比较大的时候，会耗费大量的空间和时间。

并行算法

下面使用了OpenMP，修改后的矩阵相乘程序使最外面的迭代在多个线程间静态分割.
#pragma omp paralllel shared(a,b,c,dim) num_threads(2)
       #pragma omp for schedule(static) private(j,k)
for(i=0;i<dim;i++)
{  for(j=0;j<dim;j++)
    {  c(i,j)=0;
       for(k=0;k<dim;k++)
       {
         c(i ,j)+=a(i ,k)*b(k,j);
       } 
    } 
}

Static调度类的一般形式为schedule（static[,chunk-size]）。它将迭代空间分割成大小为chunk-size的相等的块，并将这些块以循环的方式分配给线程。如果没有指定chunk-size，迭代空间就分成与线程数目一样多的块，每个块分配给每个线程。

OpenMP实现矩阵乘法的相关说明

矩阵乘法代码实现

#include <iostream>
#include "cstdlib"
#include "random"
#include "ctime"
#include "chrono"
#include "omp.h"

using namespace std;
using namespace std::chrono;

#define M 1000
#define N 1000
#define ThreadNumber 4

double matrix_a[M][N], matrix_b[N][M], result[M][M];

void init_matrix();

void serial_multiplication();

void parallel_multiplication();

int main() {
    init_matrix();
    auto start = system_clock::now();
    serial_multiplication();
    auto end = system_clock::now();
    auto duration = duration_cast<microseconds>(end - start);
    // 关闭了输入输出，使得测试时间更准确
    //for (int i = 0; i < M; i++)
    //{
    //    for (int j = 0; j < M; j++)
    //        cout << result[i][j] << " ";
    //    cout << endl;
    //}
    cout << "serial multiplication takes "
        << double(duration.count()) * microseconds::period::num / microseconds::period::den << " seconds" << endl;
    start = system_clock::now();
    parallel_multiplication();
    end = system_clock::now();
    duration = duration_cast<microseconds>(end - start);
    cout << "parallel multiplication takes "
        << double(duration.count()) * microseconds::period::num / microseconds::period::den << " seconds" << endl;
    return 0;
}

//generate the same matrix everytime
void init_matrix() {
    default_random_engine engine;
    uniform_real_distribution<double> u(0.0, 1.0);
    for (int i = 0; i < M; ++i) {
        for (int j = 0; j < N; ++j) {
            matrix_a[i][j] = u(engine);
        }
    }
    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < M; ++j) {
            matrix_b[i][j] = u(engine);
        }
    }
}


// 普通矩阵乘法
void serial_multiplication() {
    for (int i = 0; i < M; ++i) {
        for (int j = 0; j < M; ++j) {
            double temp = 0;
            for (int k = 0; k < N; ++k) {
                temp += matrix_a[i][k] * matrix_b[k][j];
            }
            result[i][j] = temp;
        }
    }
}

// OpenMP计算矩阵乘法
void parallel_multiplication() {
#pragma omp paralllel default(private) shared(matrix_a,matrix_b, result,N) num_threads(ThreadNumber)
    #pragma omp for schedule(static)
    for (int i = 0; i < M; ++i) {
        for (int j = 0; j < M; ++j) {
            double temp = 0;
            for (int k = 0; k < N; ++k) {
                temp += matrix_a[i][k] * matrix_b[k][j];
            }
            result[i][j] = temp;
        }
    }
}

串行和并行实现矩阵乘法的实验结果

实验结果

硬件平台：AMD Athlon(tm) 64 X2 Dual Core Processor 4400+ 2.21GHz ,
2. 00GB的内存
编译器：Visual Studio2019
操作系统：Microsoft Windows Server 2003 Enterprise Edition Service Pack 1
测试数据集合：由随机数函数产生0－100的矩阵元素。

测试结果可视化：