OpenMP矩阵乘法实现

前面已经介绍了MPI和Pthreads矩阵乘法的实现，下面将继续介绍OpenMP矩阵乘法的实现方式。其实OpenMP矩阵乘法的实现与前面的Pthreads的实现方式有共同之处，都是基于线程的并行矩阵乘法的实现，因此如果Pthreads那章的代码看明白的话，本章就会变得非常简单，代码实现也和上一章差不了多少。

串行化思路

首先通过I/O操作将两个矩阵读入数组中去，再对第二个矩阵进行转置操作，之后将计算结果存入c.txt文件中，整个过程中使用clock函数计算时间。

并行化思路

假设矩阵乘法为A * B = C。让每个线程计算乘积向量 x 的不同部分,特别地,p 个线程中的每一个线程计算x中的 1000/P 个连续的元素。这个算法通过执行以下代码来启动thread_count个线程：

# pragma omp parallel num_threads(thread_count)
omp_mult();

其中omp_mult函数是用来计算每个线程块的矩阵乘法结果。
在主线程中将矩阵 A 按行分配给每个线程块，并将矩阵B定义为全局变量，使每个线程能计算相应部分的矩阵相乘结果，然后子线程将计算结果返回给主线程，当所有线程都完成了乘法操作，从线程将终止，主线程等到所有从线程结束，最后在主线程里将矩阵返回结果存储到文件中。整个过程中使用gettimeofday计算程序运行总时间。
这里提一点，gettimeofday这个函数的用法可以上网查一下，它计算的时间在这里要比clock函数更加精确一些。

转置操作

考虑到Cache命中率的问题，本实验中对于第二个矩阵进行转置操作，提高Cache命中率来进行并行的优化。转置时只考虑对角线以上元素(不包括对角线元素)，比如对于a[i][j]，矩阵转置时就将a[i][j]与a[j][i]交换，由于对角线元素交换后不变，因此就不考虑对角线元素了。

OpenMP矩阵乘法实现

#include<stdio.h>
#include<stdlib.h>
#include<omp.h>

int thread_count;
int size, local_size;
int *a, *b, *c;
FILE *fp;

int* transpose_matrix(int *m, int size);

void Init();

void omp_mult();

int main(int argc, char* argv[]){
    int i, j;

    float time_use = 0;
    struct timeval start;
    struct timeval end;

    size = 1000;

    gettimeofday(&start, NULL);

    thread_count = strtol(argv[1], NULL, 10);

    local_size = size/thread_count;

    Init();

# pragma omp parallel num_threads(thread_count)
    omp_mult();

    fp=fopen("c.txt","w");//打开文件
    for(i=0;i<1000;i++) {//写数据
        for(j=0;j<1000;j++)
            fprintf(fp,"%d ",c[i*size+j]); 
        fputc('\n',fp); 
    }
    fclose(fp);//关闭文件

    gettimeofday(&end, NULL);
    time_use = (end.tv_sec-start.tv_sec)*1000000+(end.tv_usec-start.tv_usec);

    printf("time_use is %f\n", time_use/1000000);

    free(a);
    free(b);
    free(c);

    return 0;
}

int* transpose_matrix(int *m, int size){
    int i, j;
    for(i=0; i<size; i++){
        for(j=i+1; j<size; j++){
            int temp = m[i*size+j];
            m[i*size+j] = m[j*size+i];
            m[j*size+i] = temp;
        }
    }
    return m;
}

void Init(){

    int i, j;

    a = (int*)malloc(sizeof(int)*size*size);
    b = (int*)malloc(sizeof(int)*size*size);
    c = (int*)malloc(sizeof(int)*size*size);

    //从文件中读入矩阵
    fp=fopen("a.txt","r");//打开文件
    for(i=0;i<1000;i++) //读数据
        for(j=0;j<1000;j++)
            fscanf(fp,"%d",&a[i*size+j]);
    fclose(fp);//关闭文件

    fp=fopen("b.txt","r");

    for(i=0;i<1000;i++)
        for(j=0;j<1000;j++)
            fscanf(fp,"%d",&b[i*size+j]);
    fclose(fp);

    b = transpose_matrix(b, size);
}

void omp_mult(){
    int my_rank = omp_get_thread_num();
    //int thread_count = omp_get_num_threads();
    int i, j, k, temp;

    int my_first_row = my_rank*local_size;
    int my_last_row = (my_rank+1)*local_size - 1;

    for(i = my_first_row; i <= my_last_row; i++){
        for(j = 0; j<size; j++){
            temp = 0;
            for(k = 0; k<size; k++)
                temp += a[i*size+k] * b[j*size+k];
            c[i*size+j] = temp;
        }

    }
}

结果分析

① 执行时间分析：
开始时由1个线程增长为2个线程的过程中，执行时间接近于减半，较符合并行计算的情况，但之后随着线程数目的增多，并行计算的时间再也没有减半，基本稳定在2秒左右。
② 加速比分析：
在1个线程增加到2个线程的过程中，加速比几乎增加了两倍，符合并行计算的情况，但之后随着线程数的增加，加速比基本稳定在1.6~1.7左右，再也没有加倍的情况出现。
③ 效率分析：
程序执行效率随线程数增加在不断下降，只是下降的趋势在不断减小。
④ 原因分析：
OpenMP并行程序的测试平台为Intel Core i5 CPU，为双核CPU，即在一个处理器上集成两个运算核心，提高了运算效率，因此会比串行的执行时间要短。由于一个进程只能在一个核上执行，因此只能有两个进程并行执行，又因为多进程运行在两个CPU上，会有进程切换等操作，所以才会出现进程数增加而执行时间增加的情况。