This is the mail archive of the gcc-patches@gcc.gnu.org mailing list for the GCC project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

Re: [patch, fortran, RFC] First steps towards inlining matmul


I have done some timings

(1) with the test given below, before the patch I get (last column in Gflops)

[Book15] f90/bug% gfc -Ofast timing/matmul_tst_sys.f90 -framework Accelerate
[Book15] f90/bug% time a.out
 Time, MATMUL:    373.708008       373.69497100000001        4.2815668504139435     
 Time, MATMUL:    172.086609       23.117919000000001        69.210381782201068     
545.374u 0.537s 6:36.93 137.5%	0+0k 0+0io 0pf+0w
[Book15] f90/bug% gfc -Ofast -fexternal-blas timing/matmul_tst_sys.f90  -framework Accelerate
[Book15] f90/bug% time a.out
 Time, MATMUL:    176.327881       23.855111999999998        67.071577781735002     
 Time, MATMUL:    182.086746       24.453551000000001        65.430170039516952     
358.059u 0.471s 0:48.43 740.2%	0+0k 0+0io 0pf+0w

after the patch

[Book15] f90/bug% time a.out
 Time, MATMUL:    392.415436       392.41728799999999        4.0772923337669056     
 Time, MATMUL:    171.690399       22.905118000000002        69.853383859450091     
563.671u 0.551s 6:55.44 135.8%	0+0k 0+0io 0pf+0w
[Book15] f90/bug% gfc -Ofast -fexternal-blas timing/matmul_tst_sys.f90 -framework Accelerate
[Book15] f90/bug% time a.out
 Time, MATMUL:    392.850342       392.84190799999999        4.0728852177349673     
 Time, MATMUL:    174.534821       23.784797000000001        67.269861500184334     
566.824u 0.674s 6:56.74 136.1%	0+0k 0+0io 0pf+0w

which means that -fexternal-blas should disable the inlining.

program t2 
implicit none 
REAL time_begin, time_end 
integer, parameter :: n=2000; 
integer(8) :: ts, te, rate8, cmax8
real(8) :: elapsed
REAL(8) :: a(n,n), b(n,n), c(n,n) 
integer, parameter :: m = 100 
integer :: i 
call RANDOM_NUMBER(a) 
call RANDOM_NUMBER(b) 
call cpu_time(time_begin) 
call SYSTEM_CLOCK (ts, rate8, cmax8)
do i = 1,m 
    a(1,1) = a(1,1) + 0.1 
    c = MATMUL(a,b) 
enddo 
call SYSTEM_CLOCK (te, rate8, cmax8)
call cpu_time(time_end) 
elapsed = real(te-ts, kind=8)/real(rate8, kind=8)
PRINT *, 'Time, MATMUL: ',time_end-time_begin, elapsed , 2*m*real(n, kind=8)**3/(10**9*elapsed)
call cpu_time(time_begin) 
call SYSTEM_CLOCK (ts, rate8, cmax8)
do i = 1,m 
    a(1,1) = a(1,1) + 0.1 
    call dgemm('n','n',n, n, n, dble(1.0), a, n, b, n, dble(0.0), c, n) 
enddo 
call SYSTEM_CLOCK (te, rate8, cmax8)
call cpu_time(time_end) 
elapsed = real(te-ts, kind=8)/real(rate8, kind=8)
PRINT *, 'Time, MATMUL: ',time_end-time_begin, elapsed , 2*m*real(n, kind=8)**3/(10**9*elapsed)
end program 
! http://groups.google.com/group/comp.lang.fortran/browse_thread/thread/1cba8e6ce5080197

(2) The polyhedron test compiled with -fprotect-parens -Ofast -funroll-loops -ftree-loop-linear -fomit-frame-pointer -fwhole-program -flto runs in ~6.5s. After applying the patch below, it runs in ~66s with/without the patch.

Dominique

--- induct.f90	2005-10-11 22:53:32.000000000 +0200
+++ induct_vmc.f90	2015-04-05 19:06:30.000000000 +0200
@@ -1644,18 +1644,17 @@ contains
           coil_tmp_vector(1) = -sin(theta)
           coil_tmp_vector(2) = cos(theta)
           coil_tmp_vector(3) = 0.0_longreal
-          coil_current_vec(1) = dot_product(rotate_coil(1,:),coil_tmp_vector(:))
-          coil_current_vec(2) = dot_product(rotate_coil(2,:),coil_tmp_vector(:))
-          coil_current_vec(3) = dot_product(rotate_coil(3,:),coil_tmp_vector(:))
+          coil_current_vec = matmul(rotate_coil,coil_tmp_vector)
 !
           do j = 1, 9
               c_vector(3) = 0.5 * h_coil * z1gauss(j)
 !
 !       rotate coil vector into the global coordinate system and translate it
 !
-              rot_c_vector(1) = dot_product(rotate_coil(1,:),c_vector(:)) + dx
-              rot_c_vector(2) = dot_product(rotate_coil(2,:),c_vector(:)) + dy
-              rot_c_vector(3) = dot_product(rotate_coil(3,:),c_vector(:)) + dz
+              rot_c_vector = matmul(rotate_coil,c_vector)
+              rot_c_vector(1) = rot_c_vector(1) + dx
+              rot_c_vector(2) = rot_c_vector(2) + dy
+              rot_c_vector(3) = rot_c_vector(3) + dz
 !
               do k = 1, 9
                   q_vector(1) = 0.5_longreal * a * (x2gauss(k) + 1.0_longreal)
@@ -1664,9 +1663,7 @@ contains
 !
 !       rotate quad vector into the global coordinate system
 !
-                  rot_q_vector(1) = dot_product(rotate_quad(1,:),q_vector(:))
-                  rot_q_vector(2) = dot_product(rotate_quad(2,:),q_vector(:))
-                  rot_q_vector(3) = dot_product(rotate_quad(3,:),q_vector(:))
+                  rot_q_vector = matmul(rotate_quad,q_vector)
 !
 !       compute and add in quadrature term
 !
@@ -1756,18 +1753,17 @@ contains
           coil_tmp_vector(1) = -sin(theta)
           coil_tmp_vector(2) = cos(theta)
           coil_tmp_vector(3) = 0.0_longreal
-          coil_current_vec(1) = dot_product(rotate_coil(1,:),coil_tmp_vector(:))
-          coil_current_vec(2) = dot_product(rotate_coil(2,:),coil_tmp_vector(:))
-          coil_current_vec(3) = dot_product(rotate_coil(3,:),coil_tmp_vector(:))
+          coil_current_vec = matmul(rotate_coil,coil_tmp_vector)
 !
           do j = 1, 9
               c_vector(3) = 0.5 * h_coil * z1gauss(j)
 !
 !       rotate coil vector into the global coordinate system and translate it
 !
-              rot_c_vector(1) = dot_product(rotate_coil(1,:),c_vector(:)) + dx
-              rot_c_vector(2) = dot_product(rotate_coil(2,:),c_vector(:)) + dy
-              rot_c_vector(3) = dot_product(rotate_coil(3,:),c_vector(:)) + dz
+              rot_c_vector = matmul(rotate_coil,c_vector)
+              rot_c_vector(1) = rot_c_vector(1) + dx
+              rot_c_vector(2) = rot_c_vector(2) + dy
+              rot_c_vector(3) = rot_c_vector(3) + dz
 !
               do k = 1, 9
                   q_vector(1) = 0.5_longreal * a * (x2gauss(k) + 1.0_longreal)
@@ -1776,9 +1772,7 @@ contains
 !
 !       rotate quad vector into the global coordinate system
 !
-                  rot_q_vector(1) = dot_product(rotate_quad(1,:),q_vector(:))
-                  rot_q_vector(2) = dot_product(rotate_quad(2,:),q_vector(:))
-                  rot_q_vector(3) = dot_product(rotate_quad(3,:),q_vector(:))
+                  rot_q_vector = matmul(rotate_quad,q_vector)
 !
 !       compute and add in quadrature term
 !
@@ -2061,18 +2055,17 @@ contains
 !
 !       compute current vector for the coil in the global coordinate system
 !
-          coil_current_vec(1) = dot_product(rotate_coil(1,:),coil_tmp_vector(:))
-          coil_current_vec(2) = dot_product(rotate_coil(2,:),coil_tmp_vector(:))
-          coil_current_vec(3) = dot_product(rotate_coil(3,:),coil_tmp_vector(:))
+          coil_current_vec = matmul(rotate_coil,coil_tmp_vector)
 !
           do j = 1, 9
               c_vector(3) = 0.5 * h_coil * z1gauss(j)
 !
 !       rotate coil vector into the global coordinate system and translate it
 !
-              rot_c_vector(1) = dot_product(rotate_coil(1,:),c_vector(:)) + dx
-              rot_c_vector(2) = dot_product(rotate_coil(2,:),c_vector(:)) + dy
-              rot_c_vector(3) = dot_product(rotate_coil(3,:),c_vector(:)) + dz
+              rot_c_vector = matmul(rotate_coil,c_vector)
+              rot_c_vector(1) = rot_c_vector(1) + dx
+              rot_c_vector(2) = rot_c_vector(2) + dy
+              rot_c_vector(3) = rot_c_vector(3) + dz
 !
               do k = 1, 9
                   q_vector(1) = 0.5_longreal * a * (x2gauss(k) + 1.0_longreal)
@@ -2081,9 +2074,7 @@ contains
 !
 !       rotate quad vector into the global coordinate system
 !
-                  rot_q_vector(1) = dot_product(rotate_quad(1,:),q_vector(:))
-                  rot_q_vector(2) = dot_product(rotate_quad(2,:),q_vector(:))
-                  rot_q_vector(3) = dot_product(rotate_quad(3,:),q_vector(:))
+                  rot_q_vector = matmul(rotate_quad,q_vector)
 !
 !       compute and add in quadrature term
 !
@@ -2204,18 +2195,17 @@ contains
 !
 !       compute current vector for the coil in the global coordinate system
 !
-          coil_current_vec(1) = dot_product(rotate_coil(1,:),coil_tmp_vector(:))
-          coil_current_vec(2) = dot_product(rotate_coil(2,:),coil_tmp_vector(:))
-          coil_current_vec(3) = dot_product(rotate_coil(3,:),coil_tmp_vector(:))
+          coil_current_vec = matmul(rotate_coil,coil_tmp_vector)
 !
           do j = 1, 9
               c_vector(3) = 0.5 * h_coil * z1gauss(j)
 !
 !       rotate coil vector into the global coordinate system and translate it
 !
-              rot_c_vector(1) = dot_product(rotate_coil(1,:),c_vector(:)) + dx
-              rot_c_vector(2) = dot_product(rotate_coil(2,:),c_vector(:)) + dy
-              rot_c_vector(3) = dot_product(rotate_coil(3,:),c_vector(:)) + dz
+              rot_c_vector = matmul(rotate_coil,c_vector)
+              rot_c_vector(1) = rot_c_vector(1) + dx
+              rot_c_vector(2) = rot_c_vector(2) + dy
+              rot_c_vector(3) = rot_c_vector(3) + dz
 !
               do k = 1, 9
                   q_vector(1) = 0.5_longreal * a * (x2gauss(k) + 1.0_longreal)
@@ -2224,9 +2214,7 @@ contains
 !
 !       rotate quad vector into the global coordinate system
 !
-                  rot_q_vector(1) = dot_product(rotate_quad(1,:),q_vector(:))
-                  rot_q_vector(2) = dot_product(rotate_quad(2,:),q_vector(:))
-                  rot_q_vector(3) = dot_product(rotate_quad(3,:),q_vector(:))
+                  rot_q_vector = matmul(rotate_quad,q_vector)
 !
 !       compute and add in quadrature term
 !

> Le 5 avr. 2015 à 16:33, Thomas Koenig <tkoenig@netcologne.de> a écrit :
> 
> Hi Dominique,
> 
>> IMO the inlining of MATMUL should be restricted to small matrices (less than 4x4, 9x9
>> or 16x16 depending of your field!-)
> 
> The problem with the library function we have is that it is quite
> general; it can deal with all the complexity of assumed-shape array
> arguments.  Inlining allows the compiler to take advantage of
> contiguous memory, compile-time knowledge of array sizes, etc.
> to allow for vectorization or even complete unrolling.
> 
> Of course, if you use c=matmul(a,b) where all three are assumed-shape
> arrays, your advantage over the library function will be minimal at
> best.
> 
> It will be interesting to see at what size calling BLAS directly
> will be better.
> 
> 	Thomas


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]