Initial commit with benchmarking code for square matrix product

1 files changed, 73 insertions, 0 deletions

diff --git a/matrixproduct.f90 b/matrixproduct.f90
new file mode 100644
index 0000000..feb4cf1
--- /dev/null
+++ b/matrixproduct.f90
@@ -0,0 +1,73 @@
+program matrixproduct
+
+    !> Performs a matrix-matrix multiplication using various methods
+    !> and compares the performance of each. The following are used:
+    !>  1. Basic triple-loop (iterative algorithm)
+    !>  2. Fortran native matmul routine
+    !>  3. LAPACK/BLAS library call
+    
+    use, intrinsic :: iso_fortran_env
+    implicit none
+    external :: dgemm   !> double-precision general matrix-matrix multiplication
+    
+    real(real64), allocatable, dimension(:,:) :: A, B, C
+    real(real64) :: start, end, loop_time, matmul_time, blas_time
+    integer(int32) :: n
+
+    do n = 1000,100000,1000
+        call prep_mats(A,B,C,n)
+
+        call cpu_time(start)
+        !call triple_loop_mul(A,B,C,n)
+        call cpu_time(end)
+        loop_time = end-start
+        C = 0
+
+        call cpu_time(start)
+        C = matmul(A,B)
+        call cpu_time(end)
+        matmul_time = end-start
+        C = 0
+
+        call cpu_time(start)
+        call dgemm('N', 'N', n, n, n, 1.0_real64, A, n, B, n, 0.0_real64, C, n)
+        call cpu_time(end)
+        blas_time = end-start
+
+        write(*,'((i5),3(e16.8))') n, loop_time, matmul_time, blas_time
+    enddo
+
+contains
+
+    subroutine prep_mats(A,B,C,n)
+        real(real64), dimension(:,:), allocatable, intent(out) :: A,B,C
+        integer(int32), intent(in) :: n
+
+        allocate(A(n,n))
+        allocate(B(n,n))
+        allocate(C(n,n))
+
+        call random_number(A)
+        call random_number(B)
+        C = 0
+
+    end subroutine prep_mats
+    
+    subroutine triple_loop_mul(A,B,C,n)
+        real(real64), dimension(:,:), intent(in) :: A,B
+        real(real64), dimension(:,:), intent(out) :: C
+        integer(int32), intent(in) :: n
+        integer(int32) :: i, j, k
+
+        row: do i = 1,n
+            col: do j = 1,n
+                sum: do k = 1,n
+                    C(i,j) = C(i,j) + A(i,k)*B(k,j)
+                end do sum
+            end do col
+        end do row
+
+    end subroutine triple_loop_mul
+
+
+end program matrixproduct