稀疏矩阵与向量相乘的并行算法介绍

稀疏矩阵与向量相乘在CG法中很重要，其计算时间所占的比重与稀疏矩阵平均每行中的非零元个数有关，比如在图13.1所示的算例中，稀疏矩阵平均每行中的非零元个数大约为81个，所以所占比重很大，其并行算法的设计是决定整个算法并行效率的关键因素之一。

在进行稀疏矩阵A与稠密向量x的乘法，以得到结果y时，在每个进程上都只具有局部x向量，同时也只需要存储局部y向量。设将x、y都分成为p块，进程k上分配第k块，则在进程k上计算yk时，将需要用到其他进程上的x分量，这些分量的标号只与稀疏矩阵的结构有关。

为确定y＝Ax的通信结构，首先在每个进程k上分析yk依赖于x中的哪些分量，这些分量又处于哪些进程上。这些数据可以采用类似于系数矩阵的存储方法存储在idxr中，利用iptrr（j+1）记录依赖于第j个进程上的首分量在idxr中的位置，并利用lengr（j+1）记录所依赖的进程j上的分量个数。

之后，利用mpi_alltoall就可以得知每个进程k上需要发送到其他进程上的分量个数，同时，以此为基础，计算出进程k需要发送给其他所有进程的分量在idxs中的首地址，这里idxs记录需要发送到其他进程上的分量的标号。(www.xing528.com)

最后，利用mpi_alltoallv就可以将idxs中的对应分量标号分散出去，并在每个进程上收集得到idxs。这样，就确定了y＝Ax的通信结构。

在进行正式的y＝Ax操作时，需要将每个进程上需要用到的x分量收集到该进程上，这可以这样来完成。首先，利用idxs将需要分散的分量收集到ws中，之后通过mpi_alltoallv操作进行分散，每个进程上收集得到wr，并将wr中的每个元素存储到一个临时数组w中，使得w[idxr（j）]＝wr（j）。这样，就可以直接进行yi＝Aw的计算了。

由于与刚度矩阵的装配一样，稀疏矩阵在与向量进行并行乘法时，所用到的各个处理器之间的依赖关系，以及一个处理器上结果向量中某分量的计算需要用到的其他处理器上的分量标号等信息都将一直不变，所以可以在第一次进行刚度矩阵装配时，先收集该操作中有关通信时所需要的整型信息，再在具体矩阵向量乘法时利用这些信息来进行相关浮点数据的通信即可。这可以大大节约整型信息的计算次数、通信次数与通信时间，从而也将大大提高并行计算的效率。