矩阵乘法的并行化改造

在嵌套循环中，如果外层循环迭代次数较少时，如果将来CPU核数增加到一定程度时，创建的线程数将可能小于CPU核数。另外如果内层循环存在负载平衡的情况下，很难调度外层循环使之达到负载平衡。

       下面以矩阵乘法作为例子来讲述如何将嵌套循环并行化，以满足上述扩展性和负载平衡需求。

       一个串行的矩阵乘法的函数代码如下：

/**  矩阵串行乘法函数

    @param    int *a -  指向要相乘的第个矩阵的指针

    @param    int row_a -  矩阵a的行数

    @param    int col_a -  矩阵a的列数

    @param    int *b -  指向要相乘的第个矩阵的指针 

    @param    int row_b -  矩阵b的行数

    @param    int col_b -  矩阵b的列数

    @param    int *c -  计算结果的矩阵的指针

    @param    int c_size -  矩阵c的空间大小（总元素个数）

    @return   void -  无

*/

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++;

上面这种优化，省去了除法，效率高，但是只能在串行代码中使用，因为它存在循环迭代间的依赖关系，无法将其正确地并行化。