【C Fortran】使用SIMD加速4x4矩阵Jacobi迭代计算

watermelon

上次自然辩证法上看知乎，忽然又想起来了SIMD，当初大约两年前，小弟我还是大三上半学期，初次见到站长写的SIMD加速计算，是使用泰勒展开法计算的cos函数（精度为16阶），当初仔细研读了代码并写了注释进行学习。
具体的代码以及SIMD的相关介绍可以见该贴：传送门->https://www.0xaa55.com/forum.php ... mp;highlight=0xAA55

这次是自己进行写一个简单的SIMD加速计算，是使用的4x4的Jacobi迭代进行练手，之前想过一些矩阵的直接解法使用SIMD计算的，但是碍于小弟水平有限，感觉不是那么好用SIMD进行编写，所以我就找了一个软柿子捏。
并且注意这次的这个SIMD计算的程序只能进行4x4的矩阵运算，因为__m128可以存储4个float(32bit)，若矩阵的维度多于4，感觉需要分块进行SIMD的计算；若矩阵维度少于4，则直接使用ANSI C进行计算即可。

直接上代码：
C语言代码

1. /*
   
2. * Use the SIMD to accelerate the calculation speed of the 4X4 matrix in Jacobi Iteration
   
3. * Author       : Zhou Xingguang
   
4. * Organization : Xi'an Jiaotong University - NuTHeL
   
5. * Date         : 2020/11/13
   
6. */
   
7. 
   
8. #include <stdio.h>
   
9. #include <stdlib.h>
   
10. #include <string.h>
    
11. #include <xmmintrin.h>
    
12. #include <smmintrin.h>
    
13. #include <Windows.h>
    
14. 
    
15. #define STOP_STEP 50
    
16. #define EPSILON 1.0E-3
    
17. 
    
18. 
    
19. // 4x4 matrix
    
20. void JacobiIter(float (*matrix)[4], float *b, float *x)
    
21. {
    
22.         int i;
    
23.         float *x_front = NULL;
    
24.         __m128 _matrix_0 = _mm_loadu_ps(matrix[0]);
    
25.         __m128 _matrix_1 = _mm_loadu_ps(matrix[1]);
    
26.         __m128 _matrix_2 = _mm_loadu_ps(matrix[2]);
    
27.         __m128 _matrix_3 = _mm_loadu_ps(matrix[3]);
    
28.         __m128 _x = _mm_set1_ps(0);
    
29.         __m128 _b = _mm_loadu_ps(b);
    
30.         __m128 _x_front = _mm_set1_ps(0);
    
31.         __m128 _diagonal = _mm_set_ps(matrix[3][3],
    
32.                 matrix[2][2],
    
33.                 matrix[1][1],
    
34.                 matrix[0][0]);
    
35. 
    
36.         // DEBUG USE
    
37.         __m128 temp;
    
38.         __m128 _sum0, _sum1, _sum2, _sum3;
    
39.         __m128 _sum_final;
    
40. 
    
41.         x_front = (float *)malloc(sizeof(float)*4);
    
42.         memset(x_front, 0, sizeof(float)* 4);
    
43. 
    
44.         _x = _mm_loadu_ps(x);
    
45.         _x_front = _mm_loadu_ps(x_front);
    
46. 
    
47.         for (i = 0; i < STOP_STEP; i++)
    
48.         {
    
49.                 // only the SSE4.1 can use the dot multiply between two _m128 vectors
    
50.                 temp = _mm_mul_ps(_matrix_0, _x_front);
    
51.                 _sum0 = _mm_dp_ps(temp, _mm_set1_ps(1.0), 0xFF);
    
52. 
    
53.                 temp = _mm_mul_ps(_matrix_1, _x_front);
    
54.                 _sum1 = _mm_dp_ps(temp, _mm_set1_ps(1.0), 0xFF);
    
55. 
    
56.                 temp = _mm_mul_ps(_matrix_2, _x_front);
    
57.                 _sum2 = _mm_dp_ps(temp, _mm_set1_ps(1.0), 0xFF);
    
58. 
    
59.                 temp = _mm_mul_ps(_matrix_3, _x_front);
    
60.                 _sum3 = _mm_dp_ps(temp, _mm_set1_ps(1.0), 0xFF);
    
61. 
    
62.                 // make the multiple result.
    
63.                 __m128 m000 = _mm_move_ss(_sum1, _sum0);
    
64.                 __m128 m111 = _mm_movelh_ps(m000, _sum2);
    
65.                 __m128 m222 = _mm_shuffle_ps(m111, m111, _MM_SHUFFLE(0, 1, 2, 3));
    
66.                 _sum_final = _mm_move_ss(m222, _sum3);
    
67.                 _sum_final = _mm_shuffle_ps(_sum_final, _sum_final, _MM_SHUFFLE(0, 1, 2, 3));
    
68. 
    
69.                 temp = _mm_add_ps(_mm_sub_ps(_b, _sum_final), _mm_mul_ps(_diagonal, _x_front));
    
70.                 _x = _mm_div_ps(temp, _diagonal);
    
71. 
    
72.                 // begin to judge the 2-Norm of the subtract vector between the solution vector and its front vector
    
73.                 __m128 square = _mm_mul_ps(_mm_sub_ps(_x, _x_front), _mm_sub_ps(_x, _x_front));
    
74.                 temp = _mm_dp_ps(square, _mm_set1_ps(1.0), 0xFF);
    
75.                 __m128 sqrt = _mm_sqrt_ps(temp);
    
76.                 // if the 2-Norm is less than 1E-3, then we get the solution, return the solution.
    
77.                 if (sqrt.m128_f32[0] < EPSILON)
    
78.                 {
    
79.                         printf("[*] 2-Norm is less than 1E-3...\n");
    
80.                         goto end;
    
81.                 }
    
82. 
    
83.                 // record the front value
    
84.                 _x_front = _x;
    
85. 
    
86.                 // print the result of each iteration
    
87.                 printf("Iteration : %d\n", i + 1);
    
88.                 printf("x = %f, %f, %f, %f\n", _x.m128_f32[0],
    
89.                         _x.m128_f32[1],
    
90.                         _x.m128_f32[2],
    
91.                         _x.m128_f32[3]);
    
92.         }
    
93. end:
    
94.         _mm_storeu_ps(x, _x);
    
95.         _mm_sfence();
    
96.         _ReadWriteBarrier();
    
97.         free(x_front);
    
98.         x_front = NULL;
    
99.         return;
    
100. }
     
101. 
     
102. 
     
103. int main(void)
     
104. {
     
105.         int i, j;
     
106.         float matrix[4][4] = { { 2.52, 0.95, 1.25, -0.85 },
     
107.         { 0.39, 1.69, -0.45, 0.49 },
     
108.         { 0.55, -1.25, 1.96, -0.98 },
     
109.         { 0.23, -1.15, -0.45, 2.31 } };
     
110. 
     
111.         float b[4] = { 1.38, -0.34, 0.67, 1.52 };
     
112.         float x[4] = { 0 };   // make the initial value to zero
     
113. 
     
114.         // time consumption test
     
115.         SetThreadAffinityMask(GetCurrentThread(), 1);
     
116.         LARGE_INTEGER begin_JacobiIter;
     
117.         LARGE_INTEGER end_JacobiIter;
     
118.         LARGE_INTEGER freq;
     
119. 
     
120.         QueryPerformanceFrequency(&freq);
     
121.         QueryPerformanceCounter(&begin_JacobiIter);
     
122.        
     
123.         for (i = 0; i < 100; i++)
     
124.         {
     
125.                 for (j = 0; j < 1048576; j++)
     
126.                 {
     
127.                         JacobiIter(matrix, b, x);
     
128.                 }
     
129.         }
     
130.        
     
131.         QueryPerformanceCounter(&end_JacobiIter);
     
132.        
     
133.         // print the result
     
134.         printf("The final result: x = %f, %f, %f, %f\n", x[0], x[1], x[2], x[3]);
     
135. 
     
136.         printf("Time consumption:%f\n",
     
137.                         (float)(end_JacobiIter.QuadPart - begin_JacobiIter.QuadPart) / (float)freq.QuadPart);
     
138. 
     
139.         return 0;
     
140. }

复制代码

运行结果：
The final result: x = 0.883346, -0.513747, -0.084159, 0.298170
Time consumption:22.806086
请按任意键继续. . .

Fortran代码

1. ! Jacobi Iteration
   
2.     module JACOBI
   
3.     implicit none
   
4.         real(8), parameter :: epsilon = 1E-3
   
5.         integer, parameter :: stop_step = 1e3
   
6.     contains
   
7.         ! Iteration solver
   
8.         subroutine IterSolver(matrix, b, x)
   
9.         implicit none
   
10.             real(8), intent(inout) :: matrix(:, :)
    
11.             real(8), intent(inout) :: b(:)
    
12.             real(8), intent(inout) :: x(:)
    
13.             real(8), allocatable   :: x_front(:)    ! try to record the value of the front calculation.
    
14.             integer                :: N
    
15.             integer                :: i, j
    
16.             N = size(matrix, 1)
    
17.             allocate(x_front(N))
    
18.             ! We should to try to calculate the initial value of x
    
19.             ! and, we try to initial the x = 0 here.
    
20.             x = 0
    
21.             x_front = x
    
22.             ! begin to iteration
    
23.             do i=1, stop_step
    
24.                 !write(*, *) "step", i
    
25.                 !write(*, *) x
    
26.                 do j=1, N
    
27.                     x(j) = (b(j)-sum(matrix(j, :)*x_front)+matrix(j, j)*x_front(j)) / matrix(j, j)
    
28.                 end do
    
29.                 ! Stop condition
    
30.                 ! Try to judge the 2-Norm of the subtract vector
    
31.                 ! between the solution vector and its front vector.
    
32.                 if(sqrt(sum((x - x_front) * (x - x_front))) .LT. epsilon) then
    
33.                     !write(*, *) "iteration stop:", x
    
34.                     exit
    
35.                 end if
    
36.                 ! give the value to x_front
    
37.                 x_front = x
    
38.             end do
    
39.             deallocate(x_front)
    
40.             return
    
41.         end subroutine
    
42.     end module
    
43.     !
    
44.     ! MAIN
    
45.     !
    
46.     program main
    
47.     use JACOBI
    
48.     implicit none
    
49.         real(8) :: matrix(4, 4)
    
50.         real(8) :: b(4) = (/ 1.38, -0.34, 0.67, 1.52 /)
    
51.         real(8) :: x(4)
    
52.         real(8) :: time_start, time_end
    
53.         integer :: i, j
    
54.         matrix(1, :) = (/  2.52, 0.95, 1.25, -0.85 /)
    
55.         matrix(2, :) = (/ 0.39, 1.69, -0.45, 0.49 /)
    
56.         matrix(3, :) = (/ 0.55, -1.25, 1.96, -0.98 /)
    
57.         matrix(4, :) = (/ 0.23, -1.15, -0.45, 2.31 /)
    
58.         
    
59.         call cpu_time(time_start)
    
60.         do i=1, 100
    
61.             do j=1, 1048576
    
62.                 call IterSolver(matrix, b, x)
    
63.             end do
    
64.         end do
    
65.         call cpu_time(time_end)
    
66.         write(*, *) "The final result:", x
    
67.         write(*, *) "Time consumption:", time_end - time_start
    
68.     end program
    

复制代码

The final result:  0.883345862741093      -0.513746862467949
-8.415921190582545E-002  0.298170337224724
Time consumption:   36.2500000000000
请按任意键继续. . .

这两个程序中逻辑结构是一模一样的，并且两个程序均为单线程程序，可以看到该C代码与Fortran代码计算的结果相同，不同的只有运行时间上的差别，其中C代码中使用了高精度计时函数QPC，Fortran代码中也使用了高精度计时函数cpu_time，
最终我们可以看到在循环进行100x1048576次Jacobi迭代后，我们可以看到SIMD的代码要比Fortran代码快了将近13秒多。

【注】
1.本程序C代码是使用msvc2013的速度最优优化，同时Fortran代码使用了优化能力较强的Intel Visual Fortran2013进行了速度最优优化，具体的项目配置可以见下面打包的项目，直接用vs2013打开项目文件即可。
2.程序进行测速时，将程序中的打印语句均注释掉了，防止因为打印语句而影响了最终函数计算的测速；若要观察每一步的迭代结果，只需要将循环次数改为1，同时，将C代码中第79，87-91行的注释取消，Fortran代码中第24，25,33行注释取消即可观察到每一步运行的结果。
3.本SIMD程序使用了_mm_dp_ps的向量内积函数，需要支持SSE4.1的指令集，不过目前稍微“现代”一些的CPU都支持SSE4.1了吧（起码教研室的电脑是都支持的）！
希望各位老师同学多多批评指正！

0xAA55

再次阅读你的代码的时候我严重怀疑你的Fortran的real实数类型其实是C的double双精度浮点数类型。你这不科学啊，要不要考虑写一个对双精度浮点数做加速运算的方法来测试测试？

说起来，AVX指令集就是做这个的，4个double同时算。你可以试试AVX。

0xAA55

watermelon 发表于 2020-11-16 18:43
哦哦好的，小弟我是这样想的，最后我有一个_mm_storeu_ps(x, _x)，然后我感觉就需要有一个_mm_sfence()来 ...

我当时那份代码加这两条语句是为了保证ABI层面的兼容性，也就是把浮点数从SSE寄存器里抠出来，然后用FPU栈顶来返回。

之所以这么做，是因为我当时考虑在运行时通过判断CPUID来决定是否使用SIMD，以及使用哪个版本的SSE。

后来一方面不爽函数指针对优化效果的副作用，另一方面每次返回浮点数都要走RAM来把SSE存储的浮点数倒进FPU。这些低效率的操作抵消了我的数学库的优化作用。于是我更新了我的mathutil库，默认做法就是只进行编译时的指令集选择。

而且事实上作为数学加速库我并不考虑支持2008年以前的CPU的加速效果了。所以要么FPU，要么SSE4.2。

你不需要考虑用FPU传递浮点数运行结果，自然不需要_mm_sfence();和_ReadWriteBarrier();这两条语句了。

0xAA55

_mm_sfence();和_ReadWriteBarrier();是不是不需要了？我感觉是不需要的

watermelon

0xAA55 发表于 2020-11-16 16:03
_mm_sfence();和_ReadWriteBarrier();是不是不需要了？我感觉是不需要的

哦哦好的，小弟我是这样想的，最后我有一个_mm_storeu_ps(x, _x)，然后我感觉就需要有一个_mm_sfence()来保证内容写入以后再返回，然后我实际上对_ReadWriteBarrier()不是特别了解，然后我看原先学习您的那个帖子，是_mm_sfence()和_ReadWriteBarrier()写在一起的，所以我就写在一起了。
对于_ReadWriteBarrier()这个函数，只找到了一个繁体字版本的MSDN：https://docs.microsoft.com/zh-tw ... rrier?view=msvc-160

watermelon

0xAA55 发表于 2020-11-17 05:45
我当时那份代码加这两条语句是为了保证ABI层面的兼容性，也就是把浮点数从SSE寄存器里抠出来，然后用FPU ...

可以可以！
原来如此！
一开始我也考虑要不要使用cpuid来判断用户的指令集是否支持SSE4.2，后来因为时间问题就没有深究；
学习到了！

watermelon

0xAA55 发表于 2020-11-17 18:08
再次阅读你的代码的时候我严重怀疑你的Fortran的real实数类型其实是C的double双精度浮点数类型。你这不科学 ...

草，我忘了，的确，real(8)是一个double类型的，real对应float，好的好的，我了解一下AVX，感谢A5大佬指导！

0xAA55

watermelon 发表于 2020-11-17 19:17
草，我忘了，的确，real(8)是一个double类型的，real对应float，好的好的，我了解一下AVX，感谢A5大佬指 ...

以及SSE可以做到2个double同时算。虽然听起来不是很“并行”。