Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Wrapper routine for using the Intel MKL xGEMM matrix multiplication routines
subroutine matrix_multiply(first_matrix, first_transposed, second_matrix, second_transposed, result_matrix)
!This subroutine multiplies two matrices using the Intel MKL xGEMM multiplication routines
!The routines can be a little tricky, so I've made it a little easier to use them consistently
!first_matrix and second_matrix are the input matrices
!first_transposed and second_transposed are integers indicating (with a 1 or 0) whether or not
!the respective matrix is transposed
!The resulting matrix is stored in result_matrix, but you'll have to allocate that outside the routine
!The sGEMM below can be changed to dGEMM for double precision without any other changes to the
!subroutine call. You'll then have to change the 'real' variables to 'double precision'.
real, intent(in) :: first_matrix(:,:), second_matrix(:,:)
real, intent(inout) :: result_matrix(:,:)
integer, intent(in) :: first_transposed, second_transposed
integer :: a1,a2,b1,b2
a1 = size(first_matrix,1)
a2 = size(first_matrix,2)
b1 = size(second_matrix,1)
b2 = size(second_matrix,2)
if (first_transposed .AND. second_transposed) then
call sGEMM('T','T',a2,b1,b2,1.0,first_matrix,a1,second_matrix,b1,0.0,result_matrix,a2)
else if (.NOT. first_transposed .AND. second_transposed) then
call sGEMM('N','T',a1,b1,b2,1.0,first_matrix,a1,second_matrix,b1,0.0,result_matrix,a1)
else if (first_transposed .AND. .NOT. second_transposed) then
call sGEMM('T','N',a2,b2,b1,1.0,first_matrix,a1,second_matrix,b1,0.0,result_matrix,a2)
else if (.NOT. first_transposed .AND. .NOT. second_transposed) then
call sGEMM('N','N',a1,b2,b1,1.0,first_matrix,a1,second_matrix,b1,0.0,result_matrix,a1)
end if
end subroutine matrix_multiply
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment