OpenMP程序设计的两个小技巧
http://blog.csdn.net/drzhouweiming/archive/2008/05/23/2472454.aspx
1、动态设置并行循环的线程数量
在实际情况中，程序可能运行在不同的机器环境里，有些机器是双核，有些机器是4核甚至更多核。并且未来硬件存在升级的可能，CPU核数会变得越来越多。如何根据机器硬件的不同来自动设置合适的线程数量就显得很重要了，否则硬件升级后程序就得进行修改，那将是一件很麻烦的事情。
比如刚开始在双核系统中开发的软件，线程数量缺省都设成2，那么当机器升级到4核或8核以后，线程数量就不能满足要求了，除非修改程序。
线程数量的设置除了要满足机器硬件升级的可扩展性外，还需要考虑程序的可扩展性，当程序运算量增加或减少后，设置的线程数量仍然能够满足要求。显然这也不能通过设置静态的线程数量来解决。
在具体计算需要使用多少线程时，主要需要考虑以下两点：
1） 当循环次数比较少时，如果分成过多数量的线程来执行，可能会使得总运行时间高于较少线程或一个线程执行的情况。并且会增加能耗。
2） 如果设置的线程数量远大于CPU核数的话，那么存在着大量的任务切换和调度等开销，也会降低整体效率。
那么如何根据循环的次数和CPU核数来动态地设置线程的数量呢？下面以一个例子来说明动态设置线程数量的算法，假设一个需要动态设置线程数的需求为：
1、 以多个线程运行时的每个线程运行的循环次数不低于4次
2、 总的运行线程数最大不超过2倍CPU核数
下面代码便是一个实现上述需求的动态设置线程数量的例子
const int MIN_ITERATOR_NUM = 4;
int ncore = omp_get_num_procs(); //获取执行核的数量
int max_tn = n / MIN_ITERATOR_NUM;
int tn = max_tn > 2*ncore ? 2*ncore : max_tn; //tn表示要设置的线程数量
#pragma omp parallel for if( tn > 1) num_threads(tn)
for ( i = 0; i < n; i++ )
{
printf("Thread Id = %ld/n", omp_get_thread_num());
//Do some work here
}
在上面代码中，根据每个线程运行的循环次数不低于4次，先计算出最大可能的线程数max_tn，然后计算需要的线程数量tn，tn的值等于max_tn和2倍CPU核数中的较小值。
然后在parallel for构造中使用if子句来判断tn是否大于1，大于1时使用单个线程，否则使用tn个线程，，这样就使得设置的线程数量满足了需求中的条件。
比如在一个双核CPU上，n=64，最终会以2倍CPU核数（4个）线程运行，而不会以max_tn = 64/4＝16个线程运行。
在实际情况中，当然不能每个循环都象上面一样写几行代码来计算一遍，可以将其写成一个独立的功能函数如下：
const int g_ncore = omp_get_num_procs(); //获取执行核的数量
int dtn(int n, int min_n)
{
int max_tn = n / min_n;
int tn = max_tn > g_ncore ? g_ncore : max_tn; //tn表示要设置的线程数量
if ( tn < 1 )
{
tn = 1;
}
return tn;
}
这样每次并行化循环时就可以直接使用函数dtn()来获取合适的线程数量，前面的代码可以简写成如下形式：
#pragma omp parallel for num_threads(dtn(n, MIN_ITERATOR_NUM))
for ( i = 0; i < n; i++ )
{
printf("Thread Id = %ld/n", omp_get_thread_num());
//Do some work here
}
当然具体设置多少线程要视情况而定的，一般情况下线程数量刚好等于CPU核数可以取得比较好的性能，因为线程数等于CPU核数时，每个核执行一个任务，没有任务切换开销。
2、嵌套循环的并行化
在嵌套循环中，如果外层循环迭代次数较少时，如果将来CPU核数增加到一定程度时，创建的线程数将可能小于CPU核数。另外如果内层循环存在负载平衡的情况下，很难调度外层循环使之达到负载平衡。
下面以矩阵乘法作为例子来讲述如何将嵌套循环并行化，以满足上述扩展性和负载平衡需求。
一个串行的矩阵乘法的函数代码如下：
void Matrix_Multiply(int *a, int row_a, int col_a,
int *b, int row_b,int col_b,
int *c, int c_size)
{
if ( col_a != row_b || c_size < row_a * col_b )
{
return;
}
int i, j, k;
//#pragma omp for private(i, j, k)
for ( i = 0; i < row_a; i++ )
{
int row_i = i * col_a;
int row_c = i * col_b;
for ( j = 0; j < col_b; j++ )
{
c[row_c + j] = 0;
for ( k = 0; k < row_b; k++ )
{
c[row_c + j] += a[row_i + k] * b[k * col_b + j];
}
}
}
}
如果在外层循环前加上OpenMP的for语句时，它就变成了一个并行的矩阵乘法函数，但是这样简单地将其并行化显然无法满足前面所述的扩展性需求。
其实可以采用一个简单的方法将最外层循环和第2层循环合并成一个循环，下面便是采用合并循环后的并行实现。
void Parallel_Matrix_Multiply(int *a, int row_a, int col_a,
int *b, int row_b,int col_b,
int *c, int c_size )
{
if ( col_a != row_b )
{
return;
}
int i, j, k;
int index;
int border = row_a * col_b;
i = 0;
j = 0;
#pragma omp parallel private(i,j,k) num_threads(dtn(border, 1))
for ( index = 0; index < border; index++ )
{
i = index / col_b;
j = index % col_b;
int row_i = i * col_a;
int row_c = i * col_b;
c[row_c+j] = 0;
for ( k = 0; k < row_b; k++ )
{
c[row_c + j] += a[row_i+k] * b[k*col_b+j];
}
}
}
从上面代码可以看出，合并后的循环边界border = row_a * col_b;即等于原来两个循环边界之积，然后在循环中计算出原来的外层循环和第2层循环的迭代变量i和j，采用除法和取余来求出i和j的值。
需要注意的是，上面求i和j的值必须要保证循环迭代的独立性，即不能有循环迭代间的依赖关系。不能将求i和j值的过程优化成如下的形式：
if ( j == col_b )
{
j = 0;
i++;
}
// …… 此处代表实际的矩阵乘法代码
j++;
上面这种优化，省去了除法，效率高，但是只能在串行代码中使用，因为它存在循环迭代间的依赖关系，无法将其正确地并行化。
==========================================
OpenMP中的任务调度
http://blog.csdn.net/drzhouweiming/archive/2007/10/26/1844762.aspx
OpenMP中，任务调度主要用于并行的for循环中，当循环中每次迭代的计算量不相等时，如果简单地给各个线程分配相同次数的迭代的话，会造成各个线程计算负载不均衡，这会使得有些线程先执行完，有些后执行完，造成某些CPU核空闲，影响程序性能。例如以下代码：
int i, j;
int a[100][100] = {0};
for ( i =0; i < 100; i++)
{
for( j = i; j < 100; j++ )
{
a[i][j] = i*j;
}
}
如果将最外层循环并行化的话，比如使用4个线程，如果给每个线程平均分配25次循环迭代计算的话，显然i＝0和i＝99的计算量相差了100倍，那么各个线程间可能出现较大的负载不平衡情况。为了解决这些问题，OpenMP中提供了几种对for循环并行化的任务调度方案。
在OpenMP中，对for循环并行化的任务调度使用schedule子句来实现，下面介绍schedule字句的用法。
1.1.1 Schedule子句用法
schedule子句的使用格式为：
schedule(type[,size])
schedule有两个参数：type和size，size参数是可选的。
1.type参数
表示调度类型，有四种调度类型如下：
· dynamic
· guided
· runtime
· static
这四种调度类型实际上只有static、dynamic、guided三种调度方式，runtime实际上是根据环境变量来选择前三种中的某中类型。
run-sched-var
2. size参数 (可选)
size参数表示循环迭代次数，size参数必须是整数。static、dynamic、guided三种调度方式都可以使用size参数，也可以不使用size参数。当type参数类型为runtime时，size参数是非法的（不需要使用，如果使用的话编译器会报错）。
1.1.2静态调度(static)
当parallel for编译指导语句没有带schedule子句时，大部分系统中默认采用static调度方式，这种调度方式非常简单。假设有n次循环迭代，t个线程，那么给每个线程静态分配大约n/t次迭代计算。这里为什么说大约分配n/t次呢？因为n/t不一定是整数，因此实际分配的迭代次数可能存在差1的情况，如果指定了size参数的话，那么可能相差一个size。
静态调度时可以不使用size参数，也可以使用size参数。
3.不使用size参数
不使用size参数时，分配给每个线程的是n/t次连续的迭代，不使用size参数的用法如下：
schedule(static)
例如以下代码:
#pragma omp parallel for schedule(static)
for(i = 0; i < 10; i++ )
{
printf("i=%d, thread_id=%d/n", i, omp_get_thread_num());
}
上面代码执行时打印的结果如下：
i=0, thread_id=0
i=1, thread_id=0
i=2, thread_id=0
i=3, thread_id=0
i=4, thread_id=0
i=5, thread_id=1
i=6, thread_id=1
i=7, thread_id=1
i=8, thread_id=1
i=9, thread_id=1
可以看出线程0得到了0～4次连续迭代，线程1得到5～9次连续迭代。注意由于多线程执行时序的随机性，每次执行时打印的结果顺序可能存在差别，后面的例子也一样。
4.使用size参数
使用size参数时，分配给每个线程的size次连续的迭代计算，用法如下：
schedule(static, size)
例如以下代码:
#pragma omp parallel for schedule(static, 2)
for(i = 0; i < 10; i++ )
{
printf("i=%d, thread_id=%d/n", i, omp_get_thread_num());
}
执行时会打印以下结果：
i=0, thread_id=0
i=1, thread_id=0
i=4, thread_id=0
i=5, thread_id=0
i=8, thread_id=0
i=9, thread_id=0
i=2, thread_id=1
i=3, thread_id=1
i=6, thread_id=1
i=7, thread_id=1
从打印结果可以看出，0、1次迭代分配给线程0，2、3次迭代分配给线程1，4、5次迭代分配给线程0，6、7次迭代分配给线程1，…。每个线程依次分配到2次连续的迭代计算。
1.1.3动态调度(dynamic)
动态调度是动态地将迭代分配到各个线程，动态调度可以使用size参数也可以不使用size参数，不使用size参数时是将迭代逐个地分配到各个线程，使用size参数时，每次分配给线程的迭代次数为指定的size次。
下面为使用动态调度不带size参数的例子：
#pragma omp parallel for schedule(dynamic)
for(i = 0; i < 10; i++ )
{
printf("i=%d, thread_id=%d/n", i, omp_get_thread_num());
}
打印结果如下：
i=0, thread_id=0
i=1, thread_id=1
i=2, thread_id=0
i=3, thread_id=1
i=5, thread_id=1
i=6, thread_id=1
i=7, thread_id=1
i=8, thread_id=1
i=4, thread_id=0
i=9, thread_id=1
下面为动态调度使用size参数的例子：
#pragma omp parallel for schedule(dynamic, 2)
for(i = 0; i < 10; i++ )
{
printf("i=%d, thread_id=%d/n", i, omp_get_thread_num());
}
打印结果如下：
i=0, thread_id=0
i=1, thread_id=0
i=4, thread_id=0
i=2, thread_id=1
i=5, thread_id=0
i=3, thread_id=1
i=6, thread_id=0
i=8, thread_id=1
i=7, thread_id=0
i=9, thread_id=1
从打印结果可以看出第0、1，4、5，6、7次迭代被分配给了线程0，第2、3，8、9次迭代则分配给了线程1，每次分配的迭代次数为2。
1.1.4guided调度（guided）
guided调度是一种采用指导性的启发式自调度方法。开始时每个线程会分配到较大的迭代块，之后分配到的迭代块会逐渐递减。迭代块的大小会按指数级下降到指定的size大小，如果没有指定size参数，那么迭代块大小最小会降到1。
例如以下代码：
#pragma omp parallel for schedule(guided,2)
for(i = 0; i < 10; i++ )
{
printf("i=%d, thread_id=%d/n", i, omp_get_thread_num());
}
打印结果如下：
i=0, thread_id=0
i=1, thread_id=0
i=2, thread_id=0
i=3, thread_id=0
i=4, thread_id=0
i=8, thread_id=0
i=9, thread_id=0
i=5, thread_id=1
i=6, thread_id=1
i=7, thread_id=1
第0、1、2、3、4次迭代被分配给线程0，第5、6、7次迭代被分配给线程1，第8、9次迭代被分配给线程0，分配的迭代次数呈递减趋势，最后一次递减到2次。
1.1.5runtime调度（rumtime）
runtime调度并不是和前面三种调度方式似的真实调度方式，它是在运行时根据环境变量OMP_SCHEDULE来确定调度类型，最终使用的调度类型仍然是上述三种调度方式中的某种。
例如在unix系统中，可以使用setenv命令来设置OMP_SCHEDULE环境变量：
setenv OMP_SCHEDULE “dynamic, 2”
上述命令设置调度类型为动态调度，动态调度的迭代次数为2。
在windows环境中，可以在”系统属性|高级|环境变量”对话框中进行设置环境变量。