14、MIMD计算中的矩阵向量乘法与FFT算法解析

week9

于 2025-08-04 15:16:47 发布

阅读量64

点赞数

CC 4.0 BY-SA版权

分类专栏：并行计算导论：从基础到应用文章标签： MIMD计算矩阵向量乘法 FFT算法

本文链接：https://blog.youkuaiyun.com/week9/article/details/150085108

并行计算导论：从基础到应用专栏收录该内容

19 篇文章 ¥499.90

订阅专栏¥69.90

会员秒杀 ¥9.9 重磅福利

超级会员免费看

MIMD计算中的矩阵向量乘法与FFT算法解析

1. MIMD计算基础与并行执行时间

在MIMD（多指令多数据）计算中，我们需要思考计算结果的存储位置，即最终的和 s 应该存储在哪一个处理器上。为了简化问题，假设所有处理器都恰好持有 n/p 个元素。那么，在 p 个处理器上执行 saxpy 操作的并行执行时间可以通过特定公式计算，以实现理想的加速比。

这里涉及到几个重要的时间参数：
- Treduce(p)(l) ：表示在 p 个处理器上对一个元素进行归约操作所需的时间。
- Tstartup ：通信启动时间。
- Tword ：传输单个数据项所需的时间。

一般来说，为了获得良好的加速比， n 必须远大于 logp ，因为 Tstartup > Tflop （ Tflop 为执行一次浮点运算的平均时间）。

2. 矩阵向量乘法的不同实现方式

2.1 MPI实现矩阵向量乘法

在分布式内存环境中，对于矩阵向量乘法 y = Ax ，我们首先要确定矩阵 A 、向量 x 和 y 的分布方式。假设 x 是一个 N 维向量， y 是一个 M 维向量，那么矩阵 A 是一个 M x N 的矩阵。这两个向量将按块分布，行块大小 br = ⌈M/p⌉ ，列块大小 bc = ⌈N/p⌉ 。

矩阵 A 按行块分布，每个处理器持有 br 行连续的 A 矩阵元素，这种数据布局被称为一维（1D）块行分布或条带挖掘。这样， A 的行与 y 的元素对齐，即对于所有的 k ， y 的第 k 个元素和 A 的第 k 行位于同一个处理器上。

具体的实现步骤如下：
1. 初始化MPI环境 ：

#include <stdio.h>
#include "mpi.h"
int N=374, M=53, one=1; 
double dzero=0.0, done=1.0;
main (int argc, char* argv[]) {
    int myid, p; 
    int m, mb, n, nb, i, i0, j, j0;
    double *A, *x_global, *x_local, *y;
    MPI_Init(&argc, &argv); 
    MPI_Comm_rank(MPI_COMM_WORLD, &myid); 
    MPI_Comm_size(MPI_COMM_WORLD, &p);

确定块大小和实际本地大小 ：

    mb = (M - 1)/p + 1; 
    nb = (N - 1)/p + 1; 
    m = mb; n = nb;
    if (myid == p-1) m = M - myid*mb;
    if (myid == p-1) n = N - myid*nb;

分配内存空间 ：

    A = (double*) malloc(mb*N*sizeof(double));
    y = (double*) malloc(mb*sizeof(double));
    x_local = (double*) malloc(nb*sizeof(double));
    x_global = (double*) malloc(N*sizeof(double));

初始化矩阵 A 和向量 x ：

    for(j=0;j<N;j++)
        for(i = 0; i < m; i++){
            A[i+j*mb] = 0.0;
            if (j == myid*mb+i) A[i+j*mb] = 1.0;
        }
    for(j=0;j<n;j++) x_local[j] = (double)(j+myid*nb);

并行矩阵向量乘法 ：

    MPI_Allgather(x_local, n, MPI_DOUBLE, x_global, nb, MPI_DOUBLE, MPI_COMM_WORLD);
    dgemv_("N", &m, &N, &done, A, &mb, x_global, &one, &dzero, y, &one);

输出结果并结束MPI环境 ：

    for(i=0;i<m;i++)
        printf ("y[%.3d] = %.10.2f\n", myid*mb+i,y[i]);
    MPI_Finalize(); 
}

需要注意的是，如果矩阵 A 按列块分布，矩阵向量乘法的计算和通信过程会有很大的变化，此时局部计算会先于通信阶段，通信阶段变为向量归约操作。

2.2 PBLAS实现矩阵向量乘法

使用PBLAS（并行基本线性代数子程序）来实现矩阵向量乘法时，矩阵 A 、向量 x 和 y 将分布在二维进程网格上。矩阵 A 采用二维块循环分布， x 作为一个 1 x M 的矩阵存储在第一行进程中， y 作为一个 M x 1 的矩阵存储在第一列进程中。

具体步骤如下：
1. 确定块大小并分配内存空间 ：

int M=15, N=20, ZERO=0, ONE=1;
int br = (M - 1)/pr + 1; 
int bc = (N - 1)/pc + 1; 
double *x = (double*) malloc(bc*sizeof(double));
double *y = (double*) malloc(br*sizeof(double));
double *A = (double*) malloc(br*bc*sizeof(double));

确定本地矩阵大小和基索引 ：

int i0 = myrow*br; j0 = mycol*bc;
int m = br; n = bc;
if (myrow == pr-1) m = M - i0;
if (mycol == pc-1) n = N - j0;

定义数组描述符 ：
在调用PBLAS或ScaLAPACK例程之前，需要定义数组描述符来提供数据的全局视图。数组描述符包含以下参数：
- 分布式矩阵的行数 M 。
- 分布式矩阵的列数 N 。
- 行块大小 br 。
- 列块大小 bc 。
- 矩阵第一行分布所在的进程行。
- 矩阵第一列分布所在的进程列。
- BLAGS上下文。
- 存储本地块的本地数组的领先维度。

descinit(descA, &M, &N, &br, &bc, &ZERO, &ZERO, &ctxt, &br, &info);
descinit(descx, &ONE, &N, &ONE, &bc, &ZERO, &ZERO, &ctxt, &ONE, &info);
descinit(descy, &M, &ONE, &br, &ONE, &ZERO, &ZERO, &ctxt, &br, &info);

执行矩阵向量乘法 ：

double alpha = 1.0; double beta = 0.0;
pdgemv("N", &M, &N, &alpha, A, &ONE, &ONE, descA, x, &ONE, &ONE, descx, &ONE, &beta, y, &ONE, &ONE, descy, &ONE);

3. ScaLAPACK的应用

ScaLAPACK例程是LAPACK例程的并行化版本，其LAPACK块大小等于二维块循环矩阵分布的块大小 b = bT = bc 。以块LU分解为例，其主要步骤如下：
1. 实际面板的LU分解 ：此任务仅需在一个进程列中进行通信，用于确定主元以及交换主元行与其他行。
2. 实际块行的计算 ：需要在主元行上进行一次集体通信，广播 b x b 的因子 LU ，同时结合广播 U21 以减少通信启动次数。
3. 剩余矩阵的秩 b 更新 ：这是计算密集型部分，在计算开始之前，需要在进程列上广播主元块行 U12 。

分布式LU分解的执行时间可以通过特定公式估算，但该公式是在一些简化假设下推导出来的。在实际应用中，需要注意以下几点：
- Tflop 并非常数，而是取决于块大小 b 。
- 对于小的块大小 b ，预测的浮点性能可能过于乐观，因为此时 Tword 低于平均水平。
- 当块大小 b 较大时，负载平衡很难实现，可能导致严重的负载不均衡。

通过在Beowulf集群上的实验，我们发现选择进程列数 pc 略大于进程行数 pr 可能是一个不错的选择。不同问题规模和处理器数量下，ScaLAPACK例程 pdgesv 的执行时间和加速比如下表所示：

问题规模 `n`	处理器数量 `p`	时间 `t(s)`	加速比 `S(p)`
500	1	0.959	1
500	2	0.686	1.4
500	4	0.788	1.2
500	8	0.684	1.4
500	16	1.12	0.86
500	32	1.12	0.86
1000	1	8.42	1.0
1000	2	4.92	1.7
1000	4	3.16	2.7
1000	8	2.31	3.7
1000	16	2.45	3.4
1000	32	2.27	3.7
2000	1	121	1
2000	2	47.3	2.7
2000	4	17.7	6.9
2000	8	10.8	11
2000	16	7.43	16
2000	32	6.53	19
5000	1	2220	1
5000	2	1262	1.8
5000	4	500	4.4
5000	8	303	7.3
5000	16	141	15
5000	32	48	46

从表中可以看出，对于小问题规模，Beowulf集群的加速比较低；而对于大问题规模（ n > 2000 ），会出现超线性加速比，这是由于随着处理器数量的增加，主存的累积大小增大，减少了对主存的访问流量。

4. MPI二维FFT示例

为了进一步说明并行计算的应用，我们来看一个MPI实现的二维FFT（快速傅里叶变换）示例。为了简化，我们考虑 n = 2^m 的二进制基数情况。二维FFT的核心思想是利用行和列变换的独立性，先对行进行独立变换，再对列进行独立变换。

具体步骤如下：
1. 按行进行一维FFT变换 ：

void FFT2D(float *a, float *w, float sign, int ny, int n) {
    int i, j, off;
    float *pa;
    for(i=0; i<ny; i++) {
        off = 2*i*n;
        pa = a + off;
        cfft2(n, pa, w, sign);
    }

进行矩阵转置 ：

    Xpose(a, n);

按列进行一维FFT变换 ：

    for(i=0; i<ny; i++) {
        off = 2*i*n;
        pa = a + off;
        cfft2(n, pa, w, sign);
    }

再次进行矩阵转置以恢复原始顺序 ：

    Xpose(a, n);
}

其中，矩阵转置函数 Xpose 的实现如下：

void Xpose(float *a, int n) {
    float t0, t1;
    static float *buf_io;
    int i, ij, is, j, step, n2, nn, size, rk, other;
    static int init = -1;
    MPI_Status stat;
    MPI_Comm_size(MPI_COMM_WORLD, &size);
    MPI_Comm_rank(MPI_COMM_WORLD, &rk);
    nn = n/size;
    n2 = 2*nn;
    if(init != n) {
        buf_io = (float *)malloc(nn*n2*sizeof(float));
        init = n;
    }
    // 本地转置第一个块
    for(j = 0; j < nn; j++) {
        for(i = 0; i < j; i++) {
            t0 = a[rk*n2 + i*2*n + j*2];
            t1 = a[rk*n2 + i*2*n + j*2 + 1];
            a[rk*n2 + i*2*n + j*2] = a[rk*n2 + j*2*n + 2*i];
            a[rk*n2 + i*2*n + j*2 + 1] = a[rk*n2 + j*2*n + 2*i + 1];
            a[rk*n2 + j*2*n + 2*i] = t0;
            a[rk*n2 + j*2*n + 2*i + 1] = t1;
        }
    }
    // 进行size - 1次通信步骤
    for (step = 1; step < size; step++) {
        other = rk ^ step; 
        ij = 0;
        for(i = 0; i < nn; i++) {
            is = other*n2 + i*2*n;
            for(j = 0; j < n2; j++) {
                buf_io[ij++] = a[is + j];
            }
        }
        MPI_Sendrecv_replace(buf_io, 2*nn*nn, MPI_FLOAT, other, rk, other, other, MPI_COMM_WORLD, &stat);
        for(i = 0; i < nn; i++) {
            for(j = 0; j < nn; j++) {
                a[other*n2 + j*2*n + i*2] = buf_io[i*n2 + j*2];
                a[other*n2 + j*2*n + i*2 + 1] = buf_io[i*n2 + j*2 + 1];
            }
        }
    }
}

需要注意的是，该算法仅适用于处理器数量 p = 2^q （ q < m - log2(n) ）的情况。对于更一般的情况，需要使用索引数字置换或 MPI_Alltoall 命令，但 MPI_Alltoall 命令的效率可能会受到实现的影响。

综上所述，通过MPI和PBLAS等工具，我们可以在分布式内存环境中实现高效的矩阵向量乘法和FFT变换。在实际应用中，需要根据具体的问题规模和硬件环境选择合适的算法和参数，以达到最佳的性能。同时，要注意负载平衡和通信开销等因素对计算效率的影响。

MIMD计算中的矩阵向量乘法与FFT算法解析

5. 不同实现方式的性能分析与对比

前面介绍了MPI和PBLAS两种实现矩阵向量乘法的方式，接下来我们对它们的性能进行分析与对比。

5.1 MPI实现的性能特点

通信开销 ：在MPI实现中， MPI_Allgather 函数用于收集所有处理器上的向量 x ，这会带来一定的通信开销。特别是当处理器数量较多或者向量 x 的规模较大时，通信时间可能会成为性能瓶颈。
负载平衡 ：假设所有处理器持有相同数量的元素，但在实际情况中，最后一个处理器可能会处理不同数量的元素，这可能会导致一定的负载不均衡。不过，通过合理调整块大小，可以在一定程度上缓解这个问题。

5.2 PBLAS实现的性能特点

全局视图 ：PBLAS实现需要定义数组描述符来提供数据的全局视图，这增加了一定的编程复杂度，但也使得在分布式环境中更方便地进行矩阵向量乘法。
通信与计算的协同 ：PBLAS通过二维块循环分布矩阵，使得通信和计算可以更好地协同进行。例如，在广播和归约操作中，可以更高效地利用网络带宽。

5.3 性能对比表格

为了更直观地对比两种实现方式的性能，我们可以列出以下表格：

实现方式	编程复杂度	通信开销	负载平衡	适用场景
MPI	较低	可能较高，取决于数据规模和处理器数量	可能存在一定不均衡	小规模问题或对编程复杂度要求较低的场景
PBLAS	较高	相对较低，通过二维块循环分布优化	较好，通过合理分布数据	大规模问题或对性能要求较高的场景

6. FFT算法的复杂度分析与优化思路

6.1 复杂度分析

二维FFT算法的复杂度主要取决于一维FFT的复杂度和矩阵转置的复杂度。
- 一维FFT复杂度 ：一维FFT算法 cfft2 的复杂度通常为$O(n log n)$，其中 n 是向量的长度。
- 矩阵转置复杂度 ：矩阵转置函数 Xpose 的复杂度与处理器数量 p 和矩阵规模 n 有关。在理想情况下，当 p = 2^q 时，转置操作的复杂度可以得到优化。

6.2 优化思路

减少通信开销 ：在矩阵转置过程中，通过合理安排通信顺序和数据布局，可以减少通信开销。例如，使用更高效的通信算法或者减少不必要的通信。
并行化优化 ：可以进一步挖掘FFT算法中的并行性，例如对多个行或列同时进行一维FFT变换。
数据局部性优化 ：提高数据的局部性可以减少内存访问时间。例如，在进行一维FFT变换时，尽量让数据在缓存中保持连续访问。

7. 实际应用中的注意事项

7.1 硬件环境的影响

处理器性能 ：不同的处理器具有不同的浮点运算性能，这会影响到 Tflop 的值。例如，在Beowulf集群和Hewlett - Packard Superdome上，单个处理器的浮点性能不同，从而导致整体性能的差异。
网络带宽 ：网络带宽决定了 Tword 的值，较高的网络带宽可以减少通信时间。在分布式环境中，网络带宽的差异会对性能产生显著影响。

7.2 参数选择的重要性

块大小 ：在矩阵向量乘法和LU分解中，块大小 b 的选择非常重要。较小的块大小可能会导致频繁的通信启动，而较大的块大小可能会导致负载不均衡。需要根据具体的问题规模和硬件环境进行调整。
处理器数量和布局 ：在分布式计算中，处理器数量和布局（如进程行和列的数量）会影响通信和计算的效率。例如，在ScaLAPACK的LU分解中，选择合适的进程列数 pc 和进程行数 pr 可以提高性能。

8. 总结与展望

8.1 总结

本文详细介绍了MIMD计算中矩阵向量乘法和FFT算法的不同实现方式，包括MPI和PBLAS实现矩阵向量乘法，以及MPI实现二维FFT算法。通过具体的代码示例和步骤说明，展示了如何在分布式内存环境中进行这些计算。同时，对不同实现方式的性能进行了分析和对比，指出了各自的优缺点和适用场景。此外，还对FFT算法的复杂度进行了分析，并提出了一些优化思路。

8.2 展望

新算法的探索 ：随着计算机技术的不断发展，可能会出现更高效的矩阵向量乘法和FFT算法。例如，基于量子计算的算法可能会带来性能的巨大提升。
硬件与软件的协同优化 ：未来可以进一步探索硬件和软件的协同优化，例如设计专门的硬件加速器来加速矩阵向量乘法和FFT算法。
跨平台应用 ：随着云计算和边缘计算的发展，需要将这些算法应用到不同的平台上。因此，如何在不同的硬件和软件环境中实现高效的计算将是一个重要的研究方向。

9. 代码示例的使用说明

为了方便读者使用本文中的代码示例，以下是一些使用说明：

9.1 MPI矩阵向量乘法代码

#include <stdio.h>
#include "mpi.h"
int N=374, M=53, one=1; 
double dzero=0.0, done=1.0;
main (int argc, char* argv[]) {
    int myid, p; 
    int m, mb, n, nb, i, i0, j, j0;
    double *A, *x_global, *x_local, *y;
    MPI_Init(&argc, &argv); 
    MPI_Comm_rank(MPI_COMM_WORLD, &myid); 
    MPI_Comm_size(MPI_COMM_WORLD, &p); 
    mb = (M - 1)/p + 1; 
    nb = (N - 1)/p + 1; 
    m = mb; n = nb;
    if (myid == p-1) m = M - myid*mb;
    if (myid == p-1) n = N - myid*nb;
    A = (double*) malloc(mb*N*sizeof(double));
    y = (double*) malloc(mb*sizeof(double));
    x_local = (double*) malloc(nb*sizeof(double));
    x_global = (double*) malloc(N*sizeof(double));
    for(j=0;j<N;j++)
        for(i = 0; i < m; i++){
            A[i+j*mb] = 0.0;
            if (j == myid*mb+i) A[i+j*mb] = 1.0;
        }
    for(j=0;j<n;j++) x_local[j] = (double)(j+myid*nb);
    MPI_Allgather(x_local, n, MPI_DOUBLE, x_global, nb, MPI_DOUBLE, MPI_COMM_WORLD);
    dgemv_("N", &m, &N, &done, A, &mb, x_global, &one, &dzero, y, &one);
    for(i=0;i<m;i++)
        printf ("y[%.3d] = %.10.2f\n", myid*mb+i,y[i]);
    MPI_Finalize(); 
}

使用步骤：
1. 确保MPI环境已经正确安装和配置。
2. 将上述代码保存为一个 .c 文件，例如 mpi_matvec.c 。
3. 使用MPI编译器进行编译，例如 mpicc mpi_matvec.c -o mpi_matvec 。
4. 运行编译后的程序，例如 mpirun -np 4 ./mpi_matvec ，其中 4 是处理器数量，可以根据需要进行调整。

9.2 MPI二维FFT代码

void Xpose(float *a, int n) {
    float t0, t1;
    static float *buf_io;
    int i, ij, is, j, step, n2, nn, size, rk, other;
    static int init = -1;
    MPI_Status stat;
    MPI_Comm_size(MPI_COMM_WORLD, &size);
    MPI_Comm_rank(MPI_COMM_WORLD, &rk);
    nn = n/size;
    n2 = 2*nn;
    if(init != n) {
        buf_io = (float *)malloc(nn*n2*sizeof(float));
        init = n;
    }
    for(j = 0; j < nn; j++) {
        for(i = 0; i < j; i++) {
            t0 = a[rk*n2 + i*2*n + j*2];
            t1 = a[rk*n2 + i*2*n + j*2 + 1];
            a[rk*n2 + i*2*n + j*2] = a[rk*n2 + j*2*n + 2*i];
            a[rk*n2 + i*2*n + j*2 + 1] = a[rk*n2 + j*2*n + 2*i + 1];
            a[rk*n2 + j*2*n + 2*i] = t0;
            a[rk*n2 + j*2*n + 2*i + 1] = t1;
        }
    }
    for (step = 1; step < size; step++) {
        other = rk ^ step; 
        ij = 0;
        for(i = 0; i < nn; i++) {
            is = other*n2 + i*2*n;
            for(j = 0; j < n2; j++) {
                buf_io[ij++] = a[is + j];
            }
        }
        MPI_Sendrecv_replace(buf_io, 2*nn*nn, MPI_FLOAT, other, rk, other, other, MPI_COMM_WORLD, &stat);
        for(i = 0; i < nn; i++) {
            for(j = 0; j < nn; j++) {
                a[other*n2 + j*2*n + i*2] = buf_io[i*n2 + j*2];
                a[other*n2 + j*2*n + i*2 + 1] = buf_io[i*n2 + j*2 + 1];
            }
        }
    }
}

void FFT2D(float *a, float *w, float sign, int ny, int n) {
    int i, j, off;
    float *pa;
    for(i=0; i<ny; i++) {
        off = 2*i*n;
        pa = a + off;
        cfft2(n, pa, w, sign);
    }
    Xpose(a, n);
    for(i=0; i<ny; i++) {
        off = 2*i*n;
        pa = a + off;
        cfft2(n, pa, w, sign);
    }
    Xpose(a, n);
}

使用步骤：
1. 确保MPI环境已经正确安装和配置。
2. 将上述代码保存为一个 .c 文件，例如 mpi_fft2d.c 。
3. 编写一个主函数来调用 FFT2D 函数，并进行必要的初始化。
4. 使用MPI编译器进行编译，例如 mpicc mpi_fft2d.c -o mpi_fft2d 。
5. 运行编译后的程序，例如 mpirun -np 4 ./mpi_fft2d ，其中 4 是处理器数量，可以根据需要进行调整。

通过以上的介绍和说明，希望读者能够更好地理解和应用矩阵向量乘法和FFT算法在分布式内存环境中的实现。在实际应用中，需要根据具体情况进行调整和优化，以达到最佳的性能。