Openmp 多线程并行计算。为什么结果每次run出来不一样？

[janeyxm]

各位，我遇到下面这样一个问题，困扰很久没有解决。不知哪位可以慷慨赐教？
我用了两个线程把计算任务分为相互独立的两个部分。分别完成之后把结果合并成为完整的结果。但两个部分内部的语句，如 do while...end do中的语句不是相互独立的，即有上下关联。
结果是，如果我计算简单的问题，最终得到正确结果，但如果计算比较复杂的例子，得到的计算结果不完整。看上去像是有些计算循环没有完成。
请问这是do while引起的吗？还是可能其他的原因？ 可是我并不知道循环次数，只能用do while。而且，我只是各用一thread完成一半的工作，至于每一半里面我不需要并行。

program test
use omp_lib
......

call omp_set_num_threads (2)
......

!$OMP PARALLEL
if (omp_get_thread_num()==0) then
do while(条件)
这里调用了一些subroutine.完成计算并储存结果到矩阵X1
end do
else
do while(条件)
调用同样的subroutine，但用了不同的变量。完成计算并储存结果到矩阵X2。
end do
end do
!$OMP BARRIER
!$OMP END PARALLEL
合并X1和X2，输出结果。

[youxia000]

do while 不适合 用omp

能不能改成for的？ 并行的要求就是 每次计算先后不相关

[janeyxm]

谢谢楼上～～！！

可好像是Fortran里面没有for循环。不过有forall。用法和for很像。forall(i<N,条件)。我的问题在于不知道循环次数的上限，而且很有可能我的循环次数N会远远超过integer的最大上限。所以我用了do while。

[intel_www]

lz的例子是按功能划分的并行，不是按数据集划分的并行。这里建议使用OpenMP中的SECTIONS字句，最好不要轻易使用OpenMP runtime library中的函数，因为这样的函数调用容易引起可移植性和缩放性方面的问题。

这里结果不正确可能是由于两部分循环操作的数据集有冲突。对于存在数据竞争的地方应该加以保护。

[intel_www]

建议你用Intel Thread Checker检查一下吧。像这种变量较多，数据共享情况不易看清楚地程序，Thread Checker可以有效的帮助检查出有没有冲突。至于怎样获取和使用Thread Checker可以参见顶置的Intel工具索引贴。

[janeyxm]

谢谢楼上～～
你是说 可能我在两个部分都给相同的变量赋值，所以造成冲突吗？ 可是我在每个部分都使用了不同的变量，如果用到相同的变量，也只是读取不赋值。
以下是部分代码。其中矩阵C，A，B和常数n,k只是读取不赋值。
另外，我想问，call omp_set_num_threads (2) ,然后code里面分配了thread 0做什么什么，thread 1做什么什么，那每个thread里面做的其实就是 按照serial program来做的对吧？ 无所谓有几个cpu？

!$OMP PARALLEL
if (omp_get_thread_num()==0) then

call eqsets(C,Y1,n,k,A,B,slctM1(1,:))

call AGGJE(C,k,Y1,L1,JS1)
m1=0
if(L1.EQ.1) then
Prod1=matmul(A,Y1)
p1=0
do j=1,n+k
if((Prod1(j).GE.B(j)).OR.((B(j)-Prod1(j)).LT.0.001)) then
p1=p1+1
endif
end do
if(p1.EQ.n+k) then
m1=m1+1
X1(:,m1)=Y1
endif
end if

s1=1
do i=0,k-3
s1=s1*slctM1(1,n+k-i)
end do

do while (s1==0)
call fslct(slctM1,n+k)
call eqsets(C,Y1,n,k,A,B,slctM1(2,:))
call AGGJE(C,k,Y1,L1,JS1)

if(L1.EQ.1) then
Prod1=matmul(A,Y1)
L1=0
call check(Prod1,B,n,k,m1,Y1,X1,L1)
if(L1==1) then
m1=m1+1
X1(:,m1)=Y1(:)
end if
end if
slctM1(1,:)=slctM1(2,:)

s1=1
do i=0,k-3
s1=s1*slctM1(1,n+k-i)
end do
end do

else
call eqsets(C,Y2,n,k,A,B,slctM2(1,:))

call AGGJE(C,k,Y2,L2,JS2)
m2=0

if(L2.EQ.1) then
Prod2=matmul(A,Y2)
p2=0
do j=1,n+k
if((Prod2(j).GE.B(j)).OR.((B(j)-Prod2(j)).LT.0.001)) then
p2=p2+1
endif
end do

if(p2.EQ.n+k) then
m2=m2+1
X2(:,m2)=Y2
endif
end if

s2=1
do i=0,k-2
s2=s2*slctM2(1,n+k-i)
end do

do while (s2==0)
call fslct(slctM2,n+k)
call eqsets(C,Y2,n,k,A,B,slctM2(2,:))

call AGGJE(C,k,Y2,L2,JS2)
if(L2.EQ.1) then
Prod2=matmul(A,Y2)
L2=0
call check(Prod2,B,n,k,m2,Y2,X2,L2)


if(L2==1) then
m2=m2+1
X2(:,m2)=Y2(:)
end if
end if

slctM2(1,:)=slctM2(2,:)

s2=1
do i=0,k-2
s2=s2*slctM2(1,n+k-i)
end do
end do
end if
!$OMP BARRIER
!$OMP End parallel

[janeyxm]

上面贴的格式不对，不方便阅读。这个好一些。

!$OMP PARALLEL
if (omp_get_thread_num()==0) then

call eqsets(C,Y1,n,k,A,B,slctM1(1,:))

call AGGJE(C,k,Y1,L1,JS1)
m1=0
if(L1.EQ.1) then
Prod1=matmul(A,Y1)
p1=0
do j=1,n+k
if((Prod1(j).GE.B(j)).OR.((B(j)-Prod1(j)).LT.0.001)) then
p1=p1+1
endif
end do
if(p1.EQ.n+k) then
m1=m1+1
X1(:,m1)=Y1
endif
end if

s1=1
do i=0,k-3
s1=s1*slctM1(1,n+k-i)
end do

do while (s1==0)
call fslct(slctM1,n+k)
call eqsets(C,Y1,n,k,A,B,slctM1(2,:))
call AGGJE(C,k,Y1,L1,JS1)

if(L1.EQ.1) then
Prod1=matmul(A,Y1)
L1=0
call check(Prod1,B,n,k,m1,Y1,X1,L1)
if(L1==1) then
m1=m1+1
X1(:,m1)=Y1(:)
end if
end if
slctM1(1,:)=slctM1(2,:)

s1=1
do i=0,k-3
s1=s1*slctM1(1,n+k-i)
end do
end do

else
call eqsets(C,Y2,n,k,A,B,slctM2(1,:))

call AGGJE(C,k,Y2,L2,JS2)
m2=0

if(L2.EQ.1) then
Prod2=matmul(A,Y2)
p2=0
do j=1,n+k
if((Prod2(j).GE.B(j)).OR.((B(j)-Prod2(j)).LT.0.001)) then
p2=p2+1
endif
end do

if(p2.EQ.n+k) then
m2=m2+1
X2(:,m2)=Y2
endif
end if

s2=1
do i=0,k-2
s2=s2*slctM2(1,n+k-i)
end do

do while (s2==0)
call fslct(slctM2,n+k)
call eqsets(C,Y2,n,k,A,B,slctM2(2,:))

call AGGJE(C,k,Y2,L2,JS2)
if(L2.EQ.1) then
Prod2=matmul(A,Y2)
L2=0
call check(Prod2,B,n,k,m2,Y2,X2,L2)


if(L2==1) then
m2=m2+1
X2(:,m2)=Y2(:)
end if
end if

slctM2(1,:)=slctM2(2,:)

s2=1
do i=0,k-2
s2=s2*slctM2(1,n+k-i)
end do
end do
end if
!$OMP BARRIER
!$OMP End parallel

[intel_www]

[Quote=引用 7 楼 janeyxm 的回复:]
谢谢 intel_www，谢谢你推荐的工具！！终于找到问题了！
原来是我定义的矩阵C，在subroutine里面不只是读取而是赋值所以冲突了～～！太感谢了！！
[/Quote]

呵呵，不用谢。对了，忘了这几个问题：“call omp_set_num_threads (2) ,然后code里面分配了thread 0做什么什么，thread 1做什么什么，那每个thread里面做的其实就是 按照serial program来做的对吧？ 无所谓有几个cpu？”

"call omp_set_num_threads(2)"表示让OpenMP的runtime library创建两个工作线程来执行并行区域内的代码，在你的例子中就是!$OMP PARALLEL和 !$OMP End parallel 之间的代码段。也就是说两个工作线程执行的代码段是同一段代码。然而，你在这段代码中用了一个if语句来根据omp_get_thread_num()的返回值来做判断。这个调用返回的是当前工作线程在OpenMP所创建的线程组中的ID(注意：不要和操作系统的线程ID相混淆)。由于不同的工作线程有不同的ID，他们会转向不同的if-else分支，一个执行if，一个执行else。于是，这两个分支会被系统调用到不同的CPU Core上并行执行。

[janeyxm]

谢谢 intel_www，谢谢你推荐的工具！！终于找到问题了！
原来是我定义的矩阵C，在subroutine里面不只是读取而是赋值所以冲突了～～！太感谢了！！

[risca]

学习了

[intel_iclifort]

又一个很有意思的老帖，顶上来