Wrapper routine for using the Intel MKL xGEMM matrix multiplication routines

mkl_matrix_multiply.f90

    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