2413445047

Committed 27 Mar 2026 12:30PM UTC coverage: 6.984% (+7.0%) from 0.0%

Build Type

push

gitlab-ci

Committed by David Dickinson

Commit Message

Add missing side and trans arguments

Coverage Stats

743 of 10639 relevant lines covered (6.98%)

16175.3 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

36.21

/matrix_multiply.fpp

!> Provides a consistent wrapper to different matrix multiply implementations
module matrix_multiply
  implicit none

  private

  public :: matmul_wrapper, multiply_methods, recommend_multiplication_method
  public :: matrix_multiply_method_type, matmul_lapack, matmul_intrinsic, matmul_custom

  !> A simple wrapper type around an integer to represent different
  !> multiply methods.
  type matrix_multiply_method_type
     private
     integer :: flag = 0
     character(len=16) :: name = "UNSET NAME"
   contains
     procedure, public :: get_flag => matrix_multiply_method_get_flag
     procedure, public :: get_name => matrix_multiply_method_get_name
  end type matrix_multiply_method_type

  !> Method calling blas' gemm directly
  type(matrix_multiply_method_type), parameter :: matmul_lapack = &
       matrix_multiply_method_type(flag = 1, name = 'Lapack')

  !> Method simply calling intrinsic matmul
  type(matrix_multiply_method_type), parameter :: matmul_intrinsic = &
       matrix_multiply_method_type(flag = 2, name = 'Intrinsic')

  !> Method using explicit loops with OpenMP parallelisation. Could do with a better name.
  type(matrix_multiply_method_type), parameter :: matmul_custom = &
       matrix_multiply_method_type(flag = 3, name = 'Custom')

  !> Array of all available methods
  type(matrix_multiply_method_type), dimension(*), parameter :: multiply_methods = [ &
#ifdef LAPACK
       matmul_lapack, &
#endif
       matmul_intrinsic, &
       matmul_custom &
       ]

  !> Provides a wrapper to different matrix multiplication methods
  interface matmul_wrapper
     module procedure matmul_wrapper_complex_2d
  end interface matmul_wrapper

  !> Used to store the default method to use if not passed to matmul_wrapper.
  !> This way we can override the default (e.g. after timing performance)
  !> whilst still having a reasonable default
  type(matrix_multiply_method_type) :: matmul_default_method = matmul_intrinsic

contains

  !> Helper routine to get access to the matrix_multiply_method's flag
  !> in a read-only way. Primarily used for testing.
  elemental integer function matrix_multiply_method_get_flag(self) result(flag)
    implicit none
    class(matrix_multiply_method_type), intent(in) :: self
    flag = self%flag
  end function matrix_multiply_method_get_flag

  !> Helper routine to get access to the matrix_multiply_method's name
  !> in a read-only way. Primarily used for testing.
  elemental character(len=16) function matrix_multiply_method_get_name(self) result(name)
    implicit none
    class(matrix_multiply_method_type), intent(in) :: self
    name = self%name
  end function matrix_multiply_method_get_name

  !> Times each available method for the given problem size and returns
  !> the flag of the fastest one.
  function recommend_multiplication_method(trial_size, repeats, display_times) result(method)
    use job_manage, only: timer_local
    use file_utils, only: error_unit
    use ran, only: ranf
    use optionals, only: get_option_with_default
    implicit none
    integer, intent(in) :: trial_size
    integer, intent(in), optional :: repeats
    logical, intent(in), optional :: display_times
    type(matrix_multiply_method_type) :: method
    real :: start_time
    integer :: method_index, j, k
    complex, dimension(:, :), allocatable :: a, b, c
    integer :: error, number_of_methods, number_of_repeats, repeat
    real, dimension(:), allocatable :: times
    logical :: should_display_times

    ! Initialise to invalid method flag
    method = matrix_multiply_method_type(flag = -1)

    ! Try to allocate matrix of trial_size, if it doesn't work report error
    ! and return.
    allocate(a(trial_size, trial_size), stat = error)
    allocate(b(trial_size, trial_size), stat = error)
    allocate(c(trial_size, trial_size), stat = error)

    if (error > 0) then
       write(error_unit(), '("Unable to allocated trial array with size ",I0)') trial_size
       return
    end if

    ! Initialise arrays
    do j = 1, trial_size
       do k = 1, trial_size
          a(j,k) = cmplx(ranf(), ranf())
          b(j,k) = a(j,k)
          c(j,k) = 0.0
       end do
    end do

    number_of_methods = size(multiply_methods)

    allocate(times(number_of_methods))
    times = 0.0

    number_of_repeats = get_option_with_default(repeats, 1)
    should_display_times = get_option_with_default(display_times, .false.)

    ! Time each method
    do method_index = 1, number_of_methods
       start_time = timer_local()
       do repeat = 1, number_of_repeats
          c = matmul_wrapper(a, b, multiply_methods(method_index))
       end do
       times(method_index) = timer_local() - start_time
    end do

    ! Report the results if requested
    if (should_display_times) then
       write(*,'("Matrix multiplication benchmark : ",I0,&
            &" repeats of square matrices with size ",I0)') number_of_repeats, trial_size
       write(*,'("Method",T18,"Total time (s)")')
       do method_index = 1, number_of_methods
          method = multiply_methods(method_index)
          write(*,'(A16,T18,F12.8)') trim(method%get_name()), times(method_index)
       end do
    end if

    ! Recommend the method with the smallest time
    method = multiply_methods( minloc(times, dim = 1))

  end function recommend_multiplication_method

  !> Wrapper to complex matrix multiplication with two 2D matrices
  function matmul_wrapper_complex_2d(a, b, method) result(output)
    use mp, only: mp_abort
#ifdef LAPACK
    use lapack_wrapper, only: gemm
#endif
    implicit none
    complex, dimension(:, :), intent(in) :: a, b
    type(matrix_multiply_method_type), optional, intent(in) :: method
    type(matrix_multiply_method_type) :: internal_method
    complex, dimension(:, :), allocatable :: output
    complex, parameter :: one = 1.0, zero = 0.0
    complex :: tmp
    integer :: d1, d2, d3
    integer :: i, j, k
    d1 = size(a, dim = 1)
    d2 = size(b, dim = 2)
    d3 = size(a, dim = 2) ! Should match size(b, dim = 1)
#ifdef GK_DEBUG
    if (d3 /= size(b, dim = 1)) error stop
#endif
    allocate(output(d1, d2))

    internal_method = matmul_default_method
    if (present(method)) internal_method = method

    select case (internal_method%flag)
    case(matmul_intrinsic%flag)
       output = matmul(a, b)
#ifdef LAPACK
    case(matmul_lapack%flag)
       call gemm('n', 'n', d1, d2, d3, one, a, d1, b, d3, zero, output, d1)
#endif
    case(matmul_custom%flag)
       !$OMP PARALLEL DO DEFAULT(none)&
       !$OMP PRIVATE(i, j, k, tmp) &
       !$OMP SHARED(output, d1, d2, d3, a, b) &
       !$OMP SCHEDULE(static)
       do j = 1, d2
          output(:, j) = zero
          do k = 1, d3
             tmp = b(k, j)
             do i = 1, d1
                output(i, j) = output(i, j) + tmp * a(i, k)
             end do
          end do
       end do
       !$OMP END PARALLEL DO
    case default
       call mp_abort("Unknown multiply method passed to matrix_multiply::matmul_wrapper", .true.)
    end select

  end function matmul_wrapper_complex_2d

end module matrix_multiply

1	!> Provides a consistent wrapper to different matrix multiply implementations
2	module matrix_multiply
3	implicit none
4
5	private
6
7	public :: matmul_wrapper, multiply_methods, recommend_multiplication_method
8	public :: matrix_multiply_method_type, matmul_lapack, matmul_intrinsic, matmul_custom
9
10	!> A simple wrapper type around an integer to represent different
11	!> multiply methods.
12	type matrix_multiply_method_type
13	private
14	integer :: flag = 0
15	character(len=16) :: name = "UNSET NAME"
16	contains
17	procedure, public :: get_flag => matrix_multiply_method_get_flag
18	procedure, public :: get_name => matrix_multiply_method_get_name
19	end type matrix_multiply_method_type
20
21	!> Method calling blas' gemm directly
22	type(matrix_multiply_method_type), parameter :: matmul_lapack = &
23	matrix_multiply_method_type(flag = 1, name = 'Lapack')
24
25	!> Method simply calling intrinsic matmul
26	type(matrix_multiply_method_type), parameter :: matmul_intrinsic = &
27	matrix_multiply_method_type(flag = 2, name = 'Intrinsic')
28
29	!> Method using explicit loops with OpenMP parallelisation. Could do with a better name.
30	type(matrix_multiply_method_type), parameter :: matmul_custom = &
31	matrix_multiply_method_type(flag = 3, name = 'Custom')
32
33	!> Array of all available methods
34	type(matrix_multiply_method_type), dimension(*), parameter :: multiply_methods = [ &
35	#ifdef LAPACK
36	matmul_lapack, &
37	#endif
38	matmul_intrinsic, &
39	matmul_custom &
40	]
41
42	!> Provides a wrapper to different matrix multiplication methods
43	interface matmul_wrapper
44	module procedure matmul_wrapper_complex_2d
45	end interface matmul_wrapper
46
47	!> Used to store the default method to use if not passed to matmul_wrapper.
48	!> This way we can override the default (e.g. after timing performance)
49	!> whilst still having a reasonable default
50	type(matrix_multiply_method_type) :: matmul_default_method = matmul_intrinsic
51
52	contains
53
54	!> Helper routine to get access to the matrix_multiply_method's flag
55	!> in a read-only way. Primarily used for testing.
56	elemental integer function matrix_multiply_method_get_flag(self) result(flag)	96✔
57	implicit none
58	class(matrix_multiply_method_type), intent(in) :: self
59	flag = self%flag	96✔
60	end function matrix_multiply_method_get_flag	96✔
61
62	!> Helper routine to get access to the matrix_multiply_method's name
63	!> in a read-only way. Primarily used for testing.
64	elemental character(len=16) function matrix_multiply_method_get_name(self) result(name)	×
65	implicit none
66	class(matrix_multiply_method_type), intent(in) :: self
67	name = self%name	×
68	end function matrix_multiply_method_get_name	×
69
70	!> Times each available method for the given problem size and returns
71	!> the flag of the fastest one.
72	function recommend_multiplication_method(trial_size, repeats, display_times) result(method)	×
73	use job_manage, only: timer_local
74	use file_utils, only: error_unit
75	use ran, only: ranf
76	use optionals, only: get_option_with_default
77	implicit none
78	integer, intent(in) :: trial_size
79	integer, intent(in), optional :: repeats
80	logical, intent(in), optional :: display_times
81	type(matrix_multiply_method_type) :: method
82	real :: start_time
83	integer :: method_index, j, k
84	complex, dimension(:, :), allocatable :: a, b, c	×
85	integer :: error, number_of_methods, number_of_repeats, repeat
86	real, dimension(:), allocatable :: times	×
87	logical :: should_display_times
88
89	! Initialise to invalid method flag
90	method = matrix_multiply_method_type(flag = -1)	×
91
92	! Try to allocate matrix of trial_size, if it doesn't work report error
93	! and return.
94	allocate(a(trial_size, trial_size), stat = error)	×
95	allocate(b(trial_size, trial_size), stat = error)	×
96	allocate(c(trial_size, trial_size), stat = error)	×
97
98	if (error > 0) then	×
99	write(error_unit(), '("Unable to allocated trial array with size ",I0)') trial_size	×
100	return	×
101	end if
102
103	! Initialise arrays
104	do j = 1, trial_size	×
105	do k = 1, trial_size	×
106	a(j,k) = cmplx(ranf(), ranf())	×
107	b(j,k) = a(j,k)	×
108	c(j,k) = 0.0	×
109	end do
110	end do
111
112	number_of_methods = size(multiply_methods)	×
113
114	allocate(times(number_of_methods))	×
115	times = 0.0	×
116
117	number_of_repeats = get_option_with_default(repeats, 1)	×
118	should_display_times = get_option_with_default(display_times, .false.)	×
119
120	! Time each method
121	do method_index = 1, number_of_methods	×
122	start_time = timer_local()	×
123	do repeat = 1, number_of_repeats	×
124	c = matmul_wrapper(a, b, multiply_methods(method_index))	×
125	end do
126	times(method_index) = timer_local() - start_time	×
127	end do
128
129	! Report the results if requested
130	if (should_display_times) then	×
131	write(*,'("Matrix multiplication benchmark : ",I0,&
132	&" repeats of square matrices with size ",I0)') number_of_repeats, trial_size	×
133	write(*,'("Method",T18,"Total time (s)")')	×
134	do method_index = 1, number_of_methods	×
135	method = multiply_methods(method_index)	×
136	write(*,'(A16,T18,F12.8)') trim(method%get_name()), times(method_index)	×
137	end do
138	end if
139
140	! Recommend the method with the smallest time
141	method = multiply_methods( minloc(times, dim = 1))	×
142
143	end function recommend_multiplication_method	×
144
145	!> Wrapper to complex matrix multiplication with two 2D matrices
146	function matmul_wrapper_complex_2d(a, b, method) result(output)	12✔
147	use mp, only: mp_abort
148	#ifdef LAPACK
149	use lapack_wrapper, only: gemm
150	#endif
151	implicit none
152	complex, dimension(:, :), intent(in) :: a, b
153	type(matrix_multiply_method_type), optional, intent(in) :: method
154	type(matrix_multiply_method_type) :: internal_method
155	complex, dimension(:, :), allocatable :: output
156	complex, parameter :: one = 1.0, zero = 0.0
157	complex :: tmp
158	integer :: d1, d2, d3
159	integer :: i, j, k
160	d1 = size(a, dim = 1)	12✔
161	d2 = size(b, dim = 2)	12✔
162	d3 = size(a, dim = 2) ! Should match size(b, dim = 1)	12✔
163	#ifdef GK_DEBUG
164	if (d3 /= size(b, dim = 1)) error stop
165	#endif
166	allocate(output(d1, d2))	12✔
167
168	internal_method = matmul_default_method	12✔
169	if (present(method)) internal_method = method	12✔
170
171	select case (internal_method%flag)	16✔
172	case(matmul_intrinsic%flag)
173	output = matmul(a, b)	4✔
174	#ifdef LAPACK
175	case(matmul_lapack%flag)
176	call gemm('n', 'n', d1, d2, d3, one, a, d1, b, d3, zero, output, d1)	4✔
177	#endif
178	case(matmul_custom%flag)
179	!$OMP PARALLEL DO DEFAULT(none)&
180	!$OMP PRIVATE(i, j, k, tmp) &
181	!$OMP SHARED(output, d1, d2, d3, a, b) &
182	!$OMP SCHEDULE(static)
183	do j = 1, d2	164✔
184	output(:, j) = zero	19,186✔
185	do k = 1, d3	19,190✔
186	tmp = b(k, j)	19,030✔
187	do i = 1, d1	2,594,644✔
188	output(i, j) = output(i, j) + tmp * a(i, k)	2,594,488✔
189	end do
190	end do
191	end do
192	!$OMP END PARALLEL DO
193	case default
194	call mp_abort("Unknown multiply method passed to matrix_multiply::matmul_wrapper", .true.)	12✔
195	end select
196
197	end function matmul_wrapper_complex_2d	12✔
198
199	end module matrix_multiply	×

gyrokinetics / utils / 2413445047

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous