在fortran下进行openmp并行计算编程

最近写水动力的程序，体系太大，必须用并行才能算的动，无奈只好找了并行编程的资料学习了。我想我没有必要在博客里开一个什么并行编程的教程之类，因为网上到处都是，我就随手记点重要的笔记吧。这里主要是openmp的~

1 临界与归约
   在涉及到openmp的并行时，最需要注意的就是被并行的区域中的公共变量，对于需要reduce的变量，尤其要注意，比如这段代码：

program main
implicit none
include 'omp_lib.h'
integer N,M,i
real(kind=8) t
N=20000
t=0.0
!$OMP PARALLEL DO
do i=1,N
t=t+float(i);
M=OMP_get_num_threads()
enddo
write(*, "('t =  ', F20.5, ' running on ', I3, ' threads.')") t,M
pause
stop
end

 串行代码可以很容易的得到正确结果：

t = 200010000.00000 running on   1 threads.
不幸的是，如果是并行的话，可能每次都得到一个不同的结果：
t = 54821260.00000 running on   8 threads.
t = 54430262.00000 running on   8 threads.
....
原因很简单，假设do被并行了两个线程，A1，A2，则每个线程都可以t，在其中一个线程访问t的时候，另一个线程修改了t，导致t的某些值“丢了”。解决方法有两种，第一种就是“临界”，就是锁定t：

!$OMP PARALLEL DO
do i = 1, N
!$OMP CRITICAL
t = t+float(i)
            !$OMP END CRITICAL
           M = OMP_get_num_threads()
enddo

这样每个时刻只有一个线程能访问这个变量。显然，这种方法会遇到“短木板瓶颈”，更高效的方法是使用“归约”：

!$OMP PARALLEL DO REDUCTION(+:t)
do i = 1, N
t = t+float(i)
           M = OMP_get_num_threads()
enddo

 此时程序会自动在内部实现储存部分和之类的操作。这个方法比临界要高效的多，这是我这里运行的结果：临界0.005s， 归约0.003s。对于大任务，速度会更快。

2 条件并行
有时，对于小的循环，多线程的消耗超过了并行的节省时间，显然这是就不值得并行了。比如

do i = 1, N
t = t+(sin(float(i))+2.0)**0.3+abs(cos(log(float(i))))**0.7
M = OMP_get_num_threads()
enddo

发现：

N              20000         5000
tserial       0.027s        0.003
tparallel    0.013s       0.004

推断在N>5000时应该并行更有效，可以加上条件编译：

!$OMP PARALLEL DO REDUCTION(+:t) if(N > 5000)

3 负载平衡
不同线程间的工作量“不平等”是个很麻烦的问题，他会大大降低程序并行效率，比如这个程序：

N = 5000
!$OMP PARALLEL DO PRIVATE(j)
do i = 1, N
do j = i, N
a(j, i) = fun(i, j)
enddo
enddo

其中fun是个费时的函数，串行与8核CPU并行的时间比较：
serial：3m28.007s；paralle：49.940s  加速比 4.1 太低了

      这个显然与CPU个数无关。分析上面的循环发现，i=1时内层需要N个循环，而i=2500时候内部仅仅N/2个循环，极其不平衡，因此可以显式指定其调动模式，改进负载平衡。NAMD中有个LDB模块就是干这个的。SCHEDULE一般格式：
SCHEDULE(type, chunk)

可以比较一下：

!$OMP PARALLEL DO SCHEDULE(static,1)         34.955s
!$OMP PARALLEL DO SCHEDULE(dynamic,1)   29.773s
     !$OMP PARALLEL DO SCHEDULE(guided,1)      53.116s
     !$OMP PARALLEL DO SCHEDULE(static,500)      48.822s
     !$OMP PARALLEL DO SCHEDULE(dynamic,500)  50.485s
!$OMP PARALLEL DO SCHEDULE(guided,500)     51.611s

      需要注意的是，实际中很难一下看出那种调度方式最好。通常需要实际试验，这还与你调用的CPU数目有关。SCHEDULE中，增大chunk可以提高缓存命中率，但是以降低负载平衡为代价的