我面临的问题是,在矩阵严格为实且向量复杂的情况下,尽快计算矩阵向量乘积。更准确地说,实矩阵是(反)对称的,并且具有非零密集块的稀疏结构。
到目前为止,我的策略是将向量分为实部和虚部,并为每个密集块计算实部和虚部的矩阵向量乘积。由于矩阵是(反)对称的,因此我打算同时计算一个块及其转置的乘积,因此我可以重复使用矩阵所在的缓存行。因此,对于每个块,我将计算4个矩阵块的实矢量积及其每个实部和虚部的转置。
我为单个块计算这4个乘积的代码最终看起来像这样:
#define no_alias __restrict__
template <typename VecType, typename MatType>
void trans_mul( const VecType* const no_alias re_in,
VecType* const no_alias re_out,
const VecType* const no_alias im_in,
VecType* const no_alias im_out,
const VecType* const no_alias re_tin,
VecType* const no_alias re_tout,
const VecType* const no_alias im_tin,
VecType* const no_alias im_tout,
const MatType* no_alias mat, // the matrix block
const int rows,
const int cols)
{
for(int j = 0; j < cols; ++j) {
for(int i = 0; i < rows; ++i) {
const auto m = *mat++; // this is mat[i, j]
re_tout[j] += m * re_tin[i]; // transposed
im_tout[j] += m * im_tin[i]; // transposed
re_out[i] -= m * re_in[j];
im_out[i] -= m * im_in[j];
}
}
}
典型的矩阵大小约为10 ^ 2。我使用带有-Ofast -march=native
的GCC 9.2.1编译我的代码。从assembly output,我可以看到编译器正在自动向量化并使用SIMD指令。
我正在与使用Fortran编写的类似代码竞争,该代码仍以大约25%的速度运行。当然,我的代码是非常幼稚的,但是,我仍然想不出比这更快的速度,因为积极的优化似乎非常有效。我还尝试使用四个cblas_dgemv
调用,但是这比我的幼稚方法要慢得多。还有什么我可以做的,还是有任何适合我情况的有用的BLAS例程?
对于相当大的矩阵(例如> = 1k),可以使用寄存器分块来提高性能(与算术运算相比,减少了内存加载/存储的数量)。对于小型矩阵,很难做得比原始代码更好。
这是带有寄存器阻止的结果代码:
#define no_alias __restrict__
#define mini(a, b) (((a)<(b)) ? (a) : (b))
template <typename VecType, typename MatType>
void trans_mul_v2( const VecType* const no_alias re_in,
VecType* const no_alias re_out,
const VecType* const no_alias im_in,
VecType* const no_alias im_out,
const VecType* const no_alias re_tin,
VecType* const no_alias re_tout,
const VecType* const no_alias im_tin,
VecType* const no_alias im_tout,
const MatType* no_alias mat, // the matrix block
const int rows,
const int cols)
{
// Block size (tuned for Clang/GCC on Intel Skylake processors)
const int si = 16;
const int sj = 8;
for(int bj = 0; bj < cols; bj+=sj) {
for(int bi = 0; bi < rows; bi+=si) {
if(bi+si <= rows && bj+sj <= cols)
{
// The underlying loops are expected to be unrolled by the compiler
for(int j = bj; j < bj+sj; ++j) {
for(int i = bi; i < bi+si; ++i) {
const auto m = mat[j*rows+i]; // Assume a column major ordering
re_tout[j] += m * re_tin[i];
im_tout[j] += m * im_tin[i];
re_out[i] -= m * re_in[j];
im_out[i] -= m * im_in[j];
}
}
}
else
{
// General case (borders)
for(int j = bj; j < mini(bj+sj,cols); ++j) {
for(int i = bi; i < mini(bi+si,rows); ++i) {
const auto m = mat[j*rows+i];
re_tout[j] += m * re_tin[i];
im_tout[j] += m * im_tin[i];
re_out[i] -= m * re_in[j];
im_out[i] -= m * im_in[j];
}
}
}
}
}
}
请注意,si
和sj
的值对执行时间有很大的影响。最佳值取决于编译器和基础架构。您可能应该针对目标计算机进行调整(如果需要性能可移植性,请将它们保持较小,尽管性能可能并非最佳)。
这里是结果(在GCC 9中使用双精度类型:]
With a row=100 and cols=100:
trans_mul_v1: 2.438 us
trans_mul_v2: 2.842 us
With a row=1k and cols=1k:
trans_mul_v1: 452 us
trans_mul_v2: 296 us
With a row=10k and cols=10k:
trans_mul_v1: 71.2 ms
trans_mul_v2: 35.9 ms