最近想写一个自己用的库来计算带状存储矩阵。主要需求就是1.带状矩阵A * 向量x；2.向量x * 向量y；3.迭代法求解线性方程组Ax=y。这过程中发现了一些疑惑特来求教。（PS：如果有带状存储矩阵可以用的库那其实更好了）
我首先先尝试写向量点积的方法来积累一些经验，看了一些用C来举例的文章，觉得矩阵乘法可能应该按,直接循环->交换循环顺序->SIMD向量化->循环展开->blocking,的顺序来优化。但是我发现手写SIMD256根本不如或几乎和编译器优化的一样。而SIMD向量化之后再循环展开更慢了，一重循环再展开倒是比直接循环快，几乎接近nalgebra了。
这把我搞郁闷了。难道rust写矩阵乘法，居然是直接循环更快么，更离谱的是，貌似nalgebra也是直接写循环？
这是因为rust的循环优化太强了么？
测试代码如下：
use std::arch::x86_64::{
__m128, _mm_setr_ps, _mm_mul_ps,
_mm_hadd_ps, _mm_setzero_ps, __m256, _mm256_setzero_ps, _mm256_setr_ps, _mm256_mul_ps,
_mm256_hadd_ps, _mm256_load_ps, _mm_load_ps, _mm256_fmadd_ps, _mm256_add_ps, _mm_fmadd_ps
};
use rand::{thread_rng, Rng};
use nalgebra as na;
fn main() {
let mut x = vec![0.0f32; 10000000];
let mut y = vec![0.0f32; 10000000];
thread_rng().fill(&mut x[..]);
thread_rng().fill(&mut y[..]);
// time dot_product
let start = std::time::Instant::now();
let sum = dot_product(&x, &y);
let elapsed = start.elapsed();
println!("dot_product: {} in {:?}", sum, elapsed);
// time unrolling_dot_product
let start = std::time::Instant::now();
let sum = unrolling_dot_product(&x, &y);
let elapsed = start.elapsed();
println!("unrolling_dot_product: {} in {:?}", sum, elapsed);
// time simd_dot_product
let start = std::time::Instant::now();
let simd_sum = simd_dot_product(&x, &y);
let simd_elapsed = start.elapsed();
println!("simd_dot_product: {} in {:?}", simd_sum, simd_elapsed);
// time simd256_dot_product
let start = std::time::Instant::now();
let simd_sum = simd256_dot_product(&x, &y);
let simd_elapsed = start.elapsed();
println!("simd256_dot_product: {} in {:?}", simd_sum, simd_elapsed);
// time simd256_unrolling_dot_product
let start = std::time::Instant::now();
let simd_sum = simd256_unrolling_dot_product(&x, &y);
let simd_elapsed = start.elapsed();
println!("simd256_unrolling_dot_product: {} in {:?}", simd_sum, simd_elapsed);
// time na_dot_product
let vx = na::DVector::from_vec(x);
let vy = na::DVector::from_vec(y);
let start = std::time::Instant::now();
let na_sum = vx.dot(&vy);
let na_elapsed = start.elapsed();
println!("na_dot_product: {} in {:?}", na_sum, na_elapsed);
}
// 简单点积
fn dot_product(x: &Vec<f32>, y: &Vec<f32>) -> f32 {
assert!(x.len() == y.len());
let mut sum = 0.0f32;
for i in 0..x.len() {
sum += x[i] * y[i];
}
sum
}
// 循环展开点积
fn unrolling_dot_product(x: &Vec<f32>, y: &Vec<f32>) -> f32 {
assert!(x.len() == y.len());
let mut sum = 0.0f32;
let mut i = 0;
let mut acc0 = 0.0f32;
let mut acc1 = 0.0f32;
let mut acc2 = 0.0f32;
let mut acc3 = 0.0f32;
let mut acc4 = 0.0f32;
let mut acc5 = 0.0f32;
let mut acc6 = 0.0f32;
let mut acc7 = 0.0f32;
while x.len() - i >= 8 {
acc0 += unsafe { x.get_unchecked(i) * y.get_unchecked(i) };
acc1 += unsafe { x.get_unchecked(i + 1) * y.get_unchecked(i + 1) };
acc2 += unsafe { x.get_unchecked(i + 2) * y.get_unchecked(i + 2) };
acc3 += unsafe { x.get_unchecked(i + 3) * y.get_unchecked(i + 3) };
acc4 += unsafe { x.get_unchecked(i + 4) * y.get_unchecked(i + 4) };
acc5 += unsafe { x.get_unchecked(i + 5) * y.get_unchecked(i + 5) };
acc6 += unsafe { x.get_unchecked(i + 6) * y.get_unchecked(i + 6) };
acc7 += unsafe { x.get_unchecked(i + 7) * y.get_unchecked(i + 7) };
i += 8;
}
sum += acc0 + acc4;
sum += acc1 + acc5;
sum += acc2 + acc6;
sum += acc3 + acc7;
sum
}
// SIMD128点积
fn simd_dot_product(x: &Vec<f32>, y: &Vec<f32>) -> f32 {
assert!(x.len() == y.len());
let mut sum = 0.0f32;
unsafe {
let mut sum4 = _mm_setzero_ps();
for i in 0..x.len() / 4 {
let x4 = _mm_load_ps(&x[i * 4]);
let y4 = _mm_load_ps(&y[i * 4]);
sum4 = _mm_fmadd_ps(x4, y4, sum4);
}
let sum4: [f32; 4] = std::mem::transmute(sum4);
sum += sum4[0] + sum4[1] + sum4[2] + sum4[3];
}
for i in (x.len() / 4) * 4..x.len() {
sum += x[i] * y[i];
}
sum
}
// SIMD256点积
fn simd256_dot_product(x: &Vec<f32>, y: &Vec<f32>) -> f32 {
assert!(x.len() == y.len());
let mut sum = 0.0f32;
unsafe {
let mut sum8 = _mm256_setzero_ps();
for i in 0..x.len() / 8 {
let x8 = _mm256_load_ps(x.get_unchecked(i * 8));
let y8 = _mm256_load_ps(y.get_unchecked(i * 8));
sum8 = _mm256_fmadd_ps(x8, y8, sum8);
}
let sum8: [f32; 8] = std::mem::transmute(sum8);
sum += sum8[0] + sum8[1] + sum8[2] + sum8[3] + sum8[4] + sum8[5] + sum8[6] + sum8[7];
}
for i in (x.len() / 8) * 8..x.len() {
sum += x[i] * y[i];
}
sum
}
// SIMD256循环展开点积
fn simd256_unrolling_dot_product(x: &Vec<f32>, y: &Vec<f32>) -> f32 {
assert!(x.len() == y.len());
let mut sum = 0.0f32;
unsafe {
let mut sum8_0 = _mm256_setzero_ps();
let mut sum8_1 = _mm256_setzero_ps();
let mut sum8_2 = _mm256_setzero_ps();
let mut sum8_3 = _mm256_setzero_ps();
let mut sum8_4 = _mm256_setzero_ps();
let mut sum8_5 = _mm256_setzero_ps();
let mut sum8_6 = _mm256_setzero_ps();
let mut sum8_7 = _mm256_setzero_ps();
let mut i = 0;
while x.len() - i >= 64 {
let x8_0 = _mm256_load_ps(&x[i]);
let y8_0 = _mm256_load_ps(&y[i]);
sum8_0 = _mm256_fmadd_ps(x8_0, y8_0, sum8_0);
let x8_1 = _mm256_load_ps(&x[i + 8]);
let y8_1 = _mm256_load_ps(&y[i + 8]);
sum8_1 = _mm256_fmadd_ps(x8_1, y8_1, sum8_1);
let x8_2 = _mm256_load_ps(&x[i + 16]);
let y8_2 = _mm256_load_ps(&y[i + 16]);
sum8_2 = _mm256_fmadd_ps(x8_2, y8_2, sum8_2);
let x8_3 = _mm256_load_ps(&x[i + 24]);
let y8_3 = _mm256_load_ps(&y[i + 24]);
sum8_3 = _mm256_fmadd_ps(x8_3, y8_3, sum8_3);
let x8_4 = _mm256_load_ps(&x[i + 32]);
let y8_4 = _mm256_load_ps(&y[i + 32]);
sum8_4 = _mm256_fmadd_ps(x8_4, y8_4, sum8_4);
let x8_5 = _mm256_load_ps(&x[i + 40]);
let y8_5 = _mm256_load_ps(&y[i + 40]);
sum8_5 = _mm256_fmadd_ps(x8_5, y8_5, sum8_5);
let x8_6 = _mm256_load_ps(&x[i + 48]);
let y8_6 = _mm256_load_ps(&y[i + 48]);
sum8_6 = _mm256_fmadd_ps(x8_6, y8_6, sum8_6);
let x8_7 = _mm256_load_ps(&x[i + 56]);
let y8_7 = _mm256_load_ps(&y[i + 56]);
sum8_7 = _mm256_fmadd_ps(x8_7, y8_7, sum8_7);
i += 64;
}
sum8_0 = _mm256_add_ps(sum8_0, sum8_4);
sum8_1 = _mm256_add_ps(sum8_1, sum8_5);
sum8_2 = _mm256_add_ps(sum8_2, sum8_6);
sum8_3 = _mm256_add_ps(sum8_3, sum8_7);
sum8_0 = _mm256_add_ps(sum8_0, sum8_2);
sum8_1 = _mm256_add_ps(sum8_1, sum8_3);
sum8_0 = _mm256_add_ps(sum8_0, sum8_1);
let sum8_0: [f32; 8] = std::mem::transmute(sum8_0);
sum += sum8_0[0] + sum8_0[1] + sum8_0[2] + sum8_0[3] + sum8_0[4] + sum8_0[5] + sum8_0[6] + sum8_0[7];
}
for i in (x.len() / 64) * 64..x.len() {
sum += x[i] * y[i];
}
sum
}
1
共 12 条评论, 1 页
评论区
写评论应该是并行度,你可以调出来一个unroll的最优值.然后测试性能时候,应该多跑几次取平均值.
--
👇
eweca-d: 大佬，感觉是不是直接循环+unroll有什么奇怪的优化我没发现吧。
--
👇
wangbyby: 还有一件事, 你需要查看下CPU型号以及对应的指令集扩展. 比如我的笔记本是i5-7300HQ, 指令集扩展是Intel® SSE4.1, Intel® SSE4.2, Intel® AVX2,但没查到具体寄存器个数. 运行结果是 unrolling_dot_product: 24575894 in 42.7105ms simd256_dot_product: 24575892 in 126.9367ms
然后查看汇编
即一个用的是avx指令集,一个是SSE指令集. 我不太了解intel给这两个指令集分配的寄存器情况, 所以我猜测性能差距来自于寄存器和计算单元数量的差异.
另外，还有unalign的，unrolling还有movups有24处，simd256有movups有4处。
--
👇
eweca-d: 我的是台式，amd tr 2950x，扩展基本全有，包括SSE4.1, SSE4.2, SSE 4A，AVX，AVX2，FMA3。寄存器有几个也没看到。我也不太懂汇编，按你这个方法查了下：unrolling那个，vmovaps没有，movaps有26处；simd256那个，vmovaps有6处，movaps有42-6=36处。
大佬，感觉是不是直接循环+unroll有什么奇怪的优化我没发现吧。
--
👇
wangbyby: 还有一件事, 你需要查看下CPU型号以及对应的指令集扩展. 比如我的笔记本是i5-7300HQ, 指令集扩展是Intel® SSE4.1, Intel® SSE4.2, Intel® AVX2,但没查到具体寄存器个数. 运行结果是 unrolling_dot_product: 24575894 in 42.7105ms simd256_dot_product: 24575892 in 126.9367ms
然后查看汇编
即一个用的是avx指令集,一个是SSE指令集. 我不太了解intel给这两个指令集分配的寄存器情况, 所以我猜测性能差距来自于寄存器和计算单元数量的差异.
我的是台式，amd tr 2950x，扩展基本全有，包括SSE4.1, SSE4.2, SSE 4A，AVX，AVX2，FMA3。寄存器有几个也没看到。我也不太懂汇编，按你这个方法查了下：unrolling那个，vmovaps没有，movaps有26处；simd256那个，vmovaps有6处，movaps有42-6=36处。
--
👇
wangbyby: 还有一件事, 你需要查看下CPU型号以及对应的指令集扩展. 比如我的笔记本是i5-7300HQ, 指令集扩展是Intel® SSE4.1, Intel® SSE4.2, Intel® AVX2,但没查到具体寄存器个数. 运行结果是 unrolling_dot_product: 24575894 in 42.7105ms simd256_dot_product: 24575892 in 126.9367ms
然后查看汇编
即一个用的是avx指令集,一个是SSE指令集. 我不太了解intel给这两个指令集分配的寄存器情况, 所以我猜测性能差距来自于寄存器和计算单元数量的差异.
还有一件事, 你需要查看下CPU型号以及对应的指令集扩展. 比如我的笔记本是i5-7300HQ, 指令集扩展是Intel® SSE4.1, Intel® SSE4.2, Intel® AVX2,但没查到具体寄存器个数. 运行结果是 unrolling_dot_product: 24575894 in 42.7105ms simd256_dot_product: 24575892 in 126.9367ms
然后查看汇编
即一个用的是avx指令集,一个是SSE指令集. 我不太了解intel给这两个指令集分配的寄存器情况, 所以我猜测性能差距来自于寄存器和计算单元数量的差异.
--
👇
eweca-d: 感谢大佬，其实我看了IR之后隐隐约约感觉到确实不太对啊。貌似unroll8的时候只是调用了simd的add和mul？那按道理直接调用simd258至少一样快吧，甚至更快（我合并了mul和add到fma指令集的fmadd）。
--
👇
wangbyby: 看起来是编译器自动做了循环向量化(可以通过查看x86汇编验证) 不过为啥unroll8和simd258有差别我还得看下
感谢大佬的回复，嗯嗯，已经改成&[T]。不过那个sum8貌似影响不大，事实上我后来尝试过某个StackOverflow的技巧：
感觉跟最后累加区别也不大，写的时候也是估计全程没转换回来，等到最后才累加这么一次。整个x和y的向量的长度很大很大，所以最后影响就不太大的样子。
--
👇
Neutron3529: 首先给一个建议：
所有函数签名，比如
(x: &Vec<f32>, y: &Vec<f32>)，里面的&Vec<T>都可以改成&[T]另外，simd256_dot_product里面的
可能有点慢，你或许应该
let mut sum8: [f32; 8] = [0f32;8];，然后用sum8存累加结果首先给一个建议：
所有函数签名，比如
(x: &Vec<f32>, y: &Vec<f32>)，里面的&Vec<T>都可以改成&[T]另外，simd256_dot_product里面的
可能有点慢，你或许应该
let mut sum8: [f32; 8] = [0f32;8];，然后用sum8存累加结果感谢大佬，其实我看了IR之后隐隐约约感觉到确实不太对啊。貌似unroll8的时候只是调用了simd的add和mul？那按道理直接调用simd258至少一样快吧，甚至更快（我合并了mul和add到fma指令集的fmadd）。
--
👇
wangbyby: 看起来是编译器自动做了循环向量化(可以通过查看x86汇编验证) 不过为啥unroll8和simd258有差别我还得看下
看起来是编译器自动做了循环向量化(可以通过查看x86汇编验证) 不过为啥unroll8和simd258有差别我还得看下
👇
eweca-d: 忘记分行了，重来下：
dot_product: 2471604.8 in 7.6143ms
unrolling_dot_product: 2499994.3 in 4.4126ms
simd_dot_product: 2498251.8 in 12.8602ms
simd256_dot_product: 2499994.3 in 7.727ms
simd256_unrolling_dot_product: 2500656.8 in 7.7704ms
na_dot_product: 2499994.3 in 4.3216ms
--
👇
wangbyby: 1. 可否把时间贴出来 2. 编译选项是否为release 3. 查看编译出来的汇编/LLVMIR了吗? 建议为每个函数加上noinline
忘记分行了，重来下：
dot_product: 2471604.8 in 7.6143ms
unrolling_dot_product: 2499994.3 in 4.4126ms
simd_dot_product: 2498251.8 in 12.8602ms
simd256_dot_product: 2499994.3 in 7.727ms
simd256_unrolling_dot_product: 2500656.8 in 7.7704ms
na_dot_product: 2499994.3 in 4.3216ms
--
👇
wangbyby: 1. 可否把时间贴出来 2. 编译选项是否为release 3. 查看编译出来的汇编/LLVMIR了吗? 建议为每个函数加上noinline
大佬，你要的运行时间，release下大概率如下：
dot_product: 2471604.8 in 7.6143ms unrolling_dot_product: 2499994.3 in 4.4126ms simd_dot_product: 2498251.8 in 12.8602ms simd256_dot_product: 2499994.3 in 7.727ms simd256_unrolling_dot_product: 2500656.8 in 7.7704ms na_dot_product: 2499994.3 in 4.3216ms
加上noinline是指
#[inline(never)]这个么，加了感觉没明显差别啊。IR我着实不懂。稍微贴一点我觉得可能重要的把。
循环展开（展开8个）：
直接循环：
手动SIMD256：
--
👇
wangbyby: 1. 可否把时间贴出来 2. 编译选项是否为release 3. 查看编译出来的汇编/LLVMIR了吗? 建议为每个函数加上noinline