方阵行列式并行化计算（OpenMP，MPI），并计算加速比

以下内容为本人并行计算课程的期末作业，有不足的地方，请多多指教！

1 实验目的

本实验的目的主要有以下三点：

1、 实现方阵行列式的计算。

2、 实现方阵行列式的并行计算，分别基于 OpenMP和 MPI。

3、 比较以上三种算法的运行时间,计算加速比。

2 实验设计

2.1 生成方阵

为方便，本实验的方阵不采取手动输入的方式，而是使用随机数来生成矩阵元素。

我定义了一个全局方阵变量——int p[100][100]。在创建方阵时，方阵的阶数N（N<100)由外部输入。然后用两层“for循环”来给方阵 p左上角 N×N个位置赋值。具体实现如下：

/*
 * 定义矩阵阶数N
 */

int N;

/*
 * 定义一个全局矩阵
 */

int p[100][100];

/*
 * 用随机数生成矩阵
 */

void create(){
	int i,j;
	for(i=0;i<N;i++)
	{
		for(j=0;j<N;j++)
	
	     {
           int a=rand()%15;//产生随机数，并赋给数组
	       p[i][j]=a;
		 }
	}
}

2.2 打印矩阵

将生成的矩阵输出，以便验算其计算行列式的正确性。具体实现如下：

/*
 * 输出矩阵
 */

void print()
{
	int i,j;
	for(i=0;i<N;i++)
    {
		for(j=0;j<N;j++)
			printf("%d ",p[i][j]);
		printf("
");
	}
}

2.3 计算矩阵行列式

计算矩阵行列式的方法有很多，本实验选择的方法是：行列式按行展开法。行列式等于它任一行的各元素与其对应的代数余子式乘积之和。代数余子式：A(ij)=（-1)^(i+j)M(ij).  (ij)为下标。某个元素的余子式等于原行列式划去该元素所在的行和列。本实验采取按第一行展开的方法。即：将高阶的行列式按第一行展开，一直重复展开行为，直到阶数为 1。上述过程可用递归完成。

2.3.1 递归实现代码

根据上面的理论，我们容易得出如下的实现方法：

/*
 * 计算行列式的函数
 */

long long mydet(int  p [100][100],int n){
	if(n==1)  //n=1返回矩阵的唯一数，停止递归
		return p[0][0];
	else{
		long long sum=0;
		for(int i=0;i<n;i++)
		{
			    int pp[100][100];//用于存放少一维的矩阵，为方便直接定义为100×100.
			for(int j=1,j1=0;j<n;j++)//去掉第一行
			{
				for(int k=0,k1=0;k<n;k++)
				{
					if(k==i)
						;//去掉对应的列
					else
					{

				     pp[j1][k1]=p[j][k];//pp为余子式	
					 k1++;
					}
				}
			    j1++;
			}
			if(i%2)
                sum+=(-1)*p[0][i]*mydet(pp,n-1);
			else
				sum+=p[0][i]*mydet(pp,n-1);
		}
		return sum;
	}
}

2.4  实现串行OpenMPMPI计算

我这里的并行主要是放在第一次的按行展开那，具体实现看代码吧。

2.4.1  串行代码

/*************************************************************************
    > File Name: matrix_det.c
    > Author: surecheun
    > Mail: surecheun@163.com
    > Created Time: 2017年12月06日 星期三 17时28分00秒
 ************************************************************************/
#include<stdlib.h>
#include<stdio.h>
#include<math.h>
#include<time.h>
#include<time.h>
/*
 * 定义矩阵阶数N
 */

int N;

/*
 * 定义一个全局矩阵
 */

int p[100][100];

/*
 * 用随机数生成矩阵
 */

void create(){
	int i,j;
	for(i=0;i<N;i++)
	{
		for(j=0;j<N;j++)
	
	     {
           int a=rand()%15;//产生随机数，并赋给数组
	       p[i][j]=a;
		 }
	}
}
      
/*
 * 输出矩阵
 */

void print()
{
	int i,j;
	for(i=0;i<N;i++)
    {
		for(j=0;j<N;j++)
			printf("%d ",p[i][j]);
		printf("
");
	}
}

/*
 * 计算行列式的函数
 */

long long mydet(int  p [100][100],int n){
	if(n==1)  //n=1返回矩阵的唯一数，停止递归
		return p[0][0];
	else{
		long long sum=0;
		for(int i=0;i<n;i++)
		{
			    int pp[100][100];//用于存放少一维的矩阵，为方便直接定义为100×100.
			for(int j=1,j1=0;j<n;j++)//去掉第一行
			{
				for(int k=0,k1=0;k<n;k++)
				{
					if(k==i)
						;//去掉对应的列
					else
					{

				     pp[j1][k1]=p[j][k];//pp为余子式	
					 k1++;
					}
				}
			    j1++;
			}
			if(i%2)
                sum+=(-1)*p[0][i]*mydet(pp,n-1);
			else
				sum+=p[0][i]*mydet(pp,n-1);
		}
		return sum;
	}
}

int main(){
	printf("N= ");
	scanf("%d",&N);
    while(N){       //如果输入N就可以继续算下去，这个设计主要为了方便获取时间数据来计算平均用时
		create();
		print();
		clock_t start_t=clock();  //开始计时
		printf("the sum of 串行 is %lld .
",mydet(p,N));
	         clock_t  end_t=clock();  //结束计时
		double   runing_t =(double)(end_t-start_t)/CLOCKS_PER_SEC;
		printf("the runing time of 串行 is %f s.",runing_t);  //输出时间
		printf("
");
		printf("N= ");
		scanf("%d",&N);
	}
 
	return 0;
}

2.4.2 OpenMP代码

/*************************************************************************
    > File Name: matrix_det_omp.c
    > Author: surecheun
    > Mail: surecheun@163.com
    > Created Time: 2017年12月07日 星期四 17时23分51秒
 ************************************************************************/

#include<stdlib.h>
#include<stdio.h>
#include<math.h>
#include<vector>
#include<time.h>
#include<omp.h>

/*
 * 定义线程数
 */
#define n_threads 2

 /*
 *定义矩阵的阶数为全局变量
 */

int N;

/*
 * 定义一个全局矩阵
 */

int p[100][100];

/*
 * 用随机数生成矩阵
 */

void create(){
	int i,j;
	for(i=0;i<N;i++)
	{
		for(j=0;j<N;j++)
	
	     {
           int a=rand()%15;//产生随机数，并赋给数组
	       p[i][j]=a;
	     }
	}
}
      
/*
 * 输出矩阵
 */

void print()
{
	int i,j;
	for(i=0;i<N;i++)
    {
		for(j=0;j<N;j++)
			printf("%d ",p[i][j]);
		printf("
");
	}
}

/*
 * 计算行列式的函数
 */

long long mydet(int  p [100][100],int n){
	if(n==1)  //n=1返回矩阵的唯一数，停止递归
		return p[0][0];
	else{
		long long sum=0;
		for(int i=0;i<n;i++)
		{
		    int pp[100][100];
			for(int j=1,j1=0;j<n;j++)//去掉第一行
			{
				for(int k=0,k1=0;k<n;k++)
				{
					if(k==i)//去掉和改数相同的列
						;
					else
					{

				     pp[j1][k1]=p[j][k];	//pp为代数余子式
					 k1++;
					}
				}
			    j1++;
			}
			if(i%2)
                                          sum+=(-1)*p[0][i]*mydet(pp,n-1);
			else
				sum+=p[0][i]*mydet(pp,n-1);
		}
		return sum;
	}
}

 
int main(){
	printf("N= ");
	scanf("%d",&N);
    while(N){       //如果输入的N>0，则继续计算
		create();  //创建矩阵
		print();   //打印创建的矩阵
		  double start1,finish1;
                    start1=omp_get_wtime();  //开始计算时间
                    long long sum=0;
             omp_set_num_threads(n_threads);//设置线程数
           #pragma omp parallel for reduction(+:sum)//并行化
       for(int i=0;i<N;i++)
		{
		    int pp[100][100];
			for(int j=1,j1=0;j<N;j++)//去掉第一行
			{
				for(int k=0,k1=0;k<N;k++)
				{
					if(k==i)//去掉和i相同的列
						;
					else
					{

				     pp[j1][k1]=p[j][k];	//pp为余子式
					 k1++;
					}
				}
			    j1++;
			}
			if(i%2)
				sum+=(-1)*p[0][i]*mydet(pp,N-1);
			else
				sum+=p[0][i]*mydet(pp,N-1);
		}
        printf("the sum of omp is %lld .
",sum);//输出结果
	    finish1=omp_get_wtime();   //结束计算时间
		double   runing_t =finish1-start1;
		printf("the runing time of opm is %f s.",runing_t);//输出时间
		printf("
");
		printf("N= ");
		scanf("%d",&N);
	}
	return 0;
	}

2.4.3 MPI实现代码

/*************************************************************************
    > File Name: matrix_det_mpi.c
    > Author: surecheun
    > Mail: surecheun@163.com
    > Created Time: 2017年12月07日 星期四 16时24分03秒
 ************************************************************************/

#include<stdlib.h>
#include<stdio.h>
#include<math.h>
#include<time.h>
#include<mpi.h>

/*
 *定义矩阵的阶数为全局变量
 */

int N;

/*
 * 定义一个全局矩阵
 */

int p[100][100];

/*
 * 用随机数生成矩阵
 */

void create(){
	int i,j;
	for(i=0;i<N;i++)
	{
		for(j=0;j<N;j++)
	
	     {
           int a=rand()%15;//产生随机数，并赋给数组
	       p[i][j]=a;
	     }
	}
}
      
/*
 * 输出矩阵
 */

void print()
{
	int i,j;
	for(i=0;i<N;i++)
    {
		for(j=0;j<N;j++)
			printf("%d ",p[i][j]);
		printf("
");
	}
}

/*
 * 计算行列式的函数
 */

long long mydet(int  p [100][100],int n){
	if(n==1)  //n=1返回矩阵的唯一数，停止递归
		return p[0][0];
	else{
		long long sum=0;
		for(int i=0;i<n;i++)
		{
		    int pp[100][100];
			for(int j=1,j1=0;j<n;j++)//去掉第一行
			{
				for(int k=0,k1=0;k<n;k++)
				{
					if(k==i)
						;
					else
					{

				     pp[j1][k1]=p[j][k];	
					 k1++;
					}
				}
			    j1++;
			}
			if(i%2)
                sum+=(-1)*p[0][i]*mydet(pp,n-1);
			else
				sum+=p[0][i]*mydet(pp,n-1);
		}
		return sum;
	}
}



int main(int argc,char *argv[]){  
	  scanf("%d",&N);
      int num_procs,my_rank;    
      double start = 0.0, stop = 0.0;  //记录时间的变量
      long long  per_procs = 0.0;  //记录每个进程算的和
      long long  result = 0.0;  //矩阵行列式结果
      MPI_Init(&argc, &argv);  
      MPI_Comm_size(MPI_COMM_WORLD, &num_procs);  
      MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);  
      if (my_rank == 0)
	  {//0号线程创建矩阵
        create(); 
        print(); //打印创建的矩阵
       }  
     start = MPI_Wtime();   //开始计算时间
 
      MPI_Bcast(&N, 1, MPI_INT, 0, MPI_COMM_WORLD);    //将矩阵大小广播给所有进程  
      for (int i = 0; i <N; i++)
	  {  
      MPI_Bcast(p[i],N, MPI_INT, 0, MPI_COMM_WORLD);  
      }           //将矩阵广播给所有进程  
   
      for (int i = my_rank; i <N; i += num_procs){ //每个线程处理不同的行和列 
            long long sum_i=0;
		    int pp[100][100];
			for(int j=1,j1=0;j<N;j++)//去掉第一行
			{
		    	for(int k=0,k1=0;k<N;k++)
				{
					if(k==i)
						;
					else
					{
				     pp[j1][k1]=p[j][k];	
					 k1++;
					}
				}
			    j1++;
			}
       if(i%2)
                sum_i=(-1)*p[0][i]*mydet(pp,N-1);
      else
		sum_i=p[0][i]*mydet(pp,N-1);
        per_procs += sum_i;  //记录每个进程的和
    }  
    MPI_Reduce(&per_procs, &result, 1, MPI_LONG_LONG_INT, MPI_SUM, 0, MPI_COMM_WORLD);//在0号进程求总和  
    if (my_rank == 0){  
        printf("the sum of mpi is %lld .
",result) ; 
        stop = MPI_Wtime();  //结束计算时间
        printf("the time of mpi is %f s
", stop - start);  
        fflush(stdout);  
    }  
    MPI_Finalize();
    return 0;  
}  

4 实验结果

4.1 正确性

4.1.1串行

结果分析，以 n=4为例，输出的矩阵为：

13	1	12	10
8	10	1	12
9	1	2	7
5	4	8	1

输出结果为：3875

和matlab计算结果一致！

4.1.2 OpenMP

结果分析，以 n=4为例，输出的矩阵为：

5	0	8	1
1	5	11	3
2	5	1	1
0	0	14	12

输出结果为：-2710

和matlab计算结果一致！

4.1.3 MPI

结果分析，以n=4为例，输出矩阵为：

9	1	2	7
5	4	8	1
0	6	7	1
11	8	12	9

输出结果为：-202

和matlab计算结果一致！

4.2 加速比

通过多次求平均，得到三种计算实现方法的计算时间（保留 3 位有效数字）如下：

N(阶数）	串行	OpenMP	MPI
9	0.0239s	0.0117s	0.0117s
10	0.195s	0.105s	0.100s

柱状图如下：