function d_matmul_(A,B) result(C)
real(8),dimension(:,:),intent(in) :: A ![N,K]
real(8),dimension(:,:),intent(in) :: B ![K,M]
real(8),dimension(size(A,1),size(B,2)) :: C ![N,M]
integer :: N,K,M
!
! C = alfa*A*B + beta*C
!
N = size(A,1)
K = size(A,2) !==size(B,1)
M = size(B,2)
if(any(shape(B)/=[K,M]))stop "d_matmul error: B has illegal shape"
if(any(shape(C)/=[N,M]))stop "d_matmul error: C has illegal shape"
!
call DGEMM('N', 'N', N, M, K, 1d0, A, N, B, K, 0d0, C, M)
!
return
end function d_matmul_