并行计算实验报告(高性能计算与网格技术)Word文档下载推荐.doc
《并行计算实验报告(高性能计算与网格技术)Word文档下载推荐.doc》由会员分享,可在线阅读,更多相关《并行计算实验报告(高性能计算与网格技术)Word文档下载推荐.doc(15页珍藏版)》请在冰豆网上搜索。
1、独立完成实验内容;
2、了解并行算法的设计基础;
3、熟悉OpenMP和MPI的编程环境以及运行环境;
4、理解不同线程数,进程数对于加速比的影响。
三、实验内容
3.1、矩阵LU分解算法的设计:
参考文档sy6.doc所使用的并行算法:
在LU分解的过程中,主要的计算是利用主行i对其余各行j,(j>
i)作初等行变换,各行计算之间没有数据相关关系,因此可以对矩阵A按行划分来实现并行计算。
考虑到在计算过程中处理器之间的负载均衡,对A采用行交叉划分:
设处理器个数为p,矩阵A的阶数为n,,对矩阵A行交叉划分后,编号为i(i=0,1,…,p-1)的处理器存有A的第i,i+p,…,i+(m-1)p行。
然后依次以第0,1,…,n-1行作为主行,将其广播给所有处理器,各处理器利用主行对其部分行向量做行变换,这实际上是各处理器轮流选出主行并广播。
若以编号为my_rank的处理器的第i行元素作为主行,并将它广播给所有处理器,则编号大于等于my_rank的处理器利用主行元素对其第i+1,…,m-1行数据做行变换,其它处理器利用主行元素对其第i,…,m-1行数据做行变换。
根据上述算法原理用代码表示如下(关键代码):
for(k=0;
k<
N;
k++)
{
for(i=0;
i<
THREADS_NUM;
i++){
thread_data_arrray[i].thread_id=i;
thread_data_arrray[i].K_number=k;
thread_data_arrray[i].chushu=a[k][k];
//创建线程
rc=pthread_create(&
pid[i],NULL,work,(void*)&
thread_data_arrray[i]);
…
}
i++){
//等待线程同步
rc=pthread_join(pid[i],&
ret);
…
}
void*work(void*arg)
{
structthread_data*my_data;
my_data=(structthread_data*)arg;
intmyid=my_data->
thread_id;
//线程ID
intmyk=my_data->
K_number;
//外层循环计数K
floatmychushu=my_data->
chushu;
//对角线的值
ints,e;
inti,j;
s=(N-myk-1)*myid/THREADS_NUM;
//确定起始循环的行数的相对位置
e=(N-myk-1)*(myid+1)/THREADS_NUM;
//确定终止循环的行数的相对位置
for(i=s+myk+1;
e+myk+1;
i++) //由于矩阵规模在缩小,找到偏移位置
a[i][myk]=a[i][myk]/mychushu;
for(j=myk+1;
j<
N;
j++)
a[i][j]=a[i][j]-a[i][myk]*a[myk][j];
//printMatrix(a);
returnNULL;
}
第一部分为入口函数,其创建指定的线程数,并根据不同的线程id按行划分矩阵,将矩阵的不同部分作为参数传递给线程,在多处理器电脑上,不同的线程并行执行,实现并行计算LU分解。
在LU分解的过程中,主要的计算是利用主行i对其余各行j,(j)i)做初等行变换,由于各行计算之间没有数据相关关系,因此可以对矩阵按行划分来实现并行算法。
考虑到计算过程中处理器负载的均衡,对矩阵采用行交叉划分;
假设处理器个数为p,矩阵的阶数为n,则每个处理器处理的行数为。
由于在OpenMP和MPI中并行算法的实现不太一样,所以接下来的两小节中我将分别针对两个编程环境设计LU分解的并行实现。
3.2、OpenMP编程
因为OpenMP是基于线程的编程模型,所以设计了一个基于多线程的OpenMP的LU分解算法,关键代码如下:
for(k=0;
omp_set_num_threads(THREADS_NUM);
#pragmaompparallelprivate(tid)
tid=omp_get_thread_num();
//当前线程ID
intmyid=tid;
printf("
helloworldfromOMPthread%d\n"
tid);
intmyk=k;
floatmychushu=A[k][k];
ints,e;
inti,j;
s=(N-myk-1)*myid/THREADS_NUM;
//确定起始循环的行数的相对位置
e=(N-myk-1)*(myid+1)/THREADS_NUM;
for(i=s+myk+1;
i++) //由于矩阵规模在缩小,找到偏移位置
{
A[i][myk]=A[i][myk]/mychushu;
for(j=myk+1;
A[i][j]=A[i][j]-A[i][myk]*A[myk][j];
//对行进行初等行变换
其主要思想为:
外层设置一个列循环,在每次循环中开设THREAD_NUMS个线程,每个线程处理的矩阵A的行为上述的m,一次循环过后则完成对应列的变换,这样在N此循环过后便可完成矩阵A的LU分解。
即L为A[k][j]中k>
j的元素,其对角线上元素为1.0,其它为0,U为A[k][j]中k<
=j的元素,其余为0。
这里如果我们使用的是一般的多线程编程,则在开启THREAD_NUMS个线程后,在下次循环开始之前,需要手动配置等待线程同步,不然可能出现错误。
但由于OpenMP使用Fork-Join并行执行模型,其会在线程队执行完以后才转到主线程执行,所以不需要等待线程同步。
详细的代码请参看附带源程序。
3.3、MPI编程
设处理器个数为p,矩阵A的阶数为n,,对矩阵A行交叉划分后,编号为i(i=0,1,…,p-1)的处理器存有A的第i,i+p,…,i+(m-1)p行。
若以编号为my_rank的处理器的第i行元素作为主行,并将它广播给所有处理器,则编号大于等于my_rank的处理器利用主行元素对其第i+1,…,m-1行数据做行变换,其它处理器利用主行元素对其第i,…,m-1行数据做行变换,计算完成后,编号为0的处理器收集各处理器中的计算结果,并从经过初等行变换的矩阵A中分离出下三角矩阵L和上三角矩阵U。
关键代码如下:
/*0号进程采用行交叉划分将矩阵A划分为大小m*M的p块子矩阵,依次发送给1至p-1号进程*/
if(my_rank==0)
{
for(i=0;
i<
m;
i++)
for(j=0;
j<
M;
j++)
a(i,j)=A((i*p),j);
if((i%p)!
=0)
{
i1=i%p;
i2=i/p+1;
MPI_Send(&
A(i,0),M,MPI_FLOAT,i1,i2,MPI_COMM_WORLD);
}
}
else
MPI_Recv(&
a(i,0),M,MPI_FLOAT,0,i+1,MPI_COMM_WORLD,&
status);
for(i=0;
for(j=0;
p;
/*j号进程负责广播主行元素*/
if(my_rank==j)
v=i*p+j;
for(k=v;
f[k]=a(i,k);
MPI_Bcast(f,M,MPI_FLOAT,my_rank,MPI_COMM_WORLD);
else
MPI_Bcast(f,M,MPI_FLOAT,j,MPI_COMM_WORLD);
/*编号小于my_rank的进程(包括my_rank本身)利用主行对其第i+1,…,m-1行数据做行变换*/
if(my_rank<
=j){
for(k=i+1;
{
a(k,v)=a(k,v)/f[v];
for(w=v+1;
w<
w++)
a(k,w)=a(k,w)-f[w]*a(k,v);
}
/*编号大于my_rank的进程利用主行对其第i,…,m-1行数据做行变换*/
if(my_rank>
j){
for(k=i;
/*0号进程从其余各进程中接收子矩阵a,得到经过变换的矩阵A*/
A(i*p,j)=a(i,j);
if(my_rank!
MPI_Send(&
a(i,0),M,MPI_FLOAT,0,i,MPI_COMM_WORLD);
for(i=1;
MPI_Recv(&
a(j,0),M,MPI_FLOAT,i,j,MPI_COMM_WORLD,&
for(k=0;
A((j*p+i),k)=a(j,k);