subroutine utility_zgemm_new(a, b, c, transa_opt, transb_opt)
!=============================================================!
! !
! Return matrix product of complex matrices a and b: !
! !
! C = Op(A) Op(B) !
! !
! transa = 'N' ==> Op(A) = A !
! transa = 'T' ==> Op(A) = transpose(A) !
! transa = 'C' ==> Op(A) = congj(transpose(A)) !
! !
! similarly for B !
! !
! Due to the use of assumed shape arrays, this routine is a !
! safer and more general replacement for the above routine !
! utility_zgemm. Consider removing utility_zgemm and using !
! utility_zgemm_new throughout. !
! !
!=============================================================!
use w90_constants, only: cmplx_0, cmplx_1
implicit none
complex(kind=dp), intent(in) :: a(:, :)
complex(kind=dp), intent(in) :: b(:, :)
complex(kind=dp), intent(out) :: c(:, :)
character(len=1), intent(in), optional :: transa_opt
character(len=1), intent(in), optional :: transb_opt
integer :: m, n, k
character(len=1) :: transa, transb
transa = 'N'
transb = 'N'
if (present(transa_opt)) transa = transa_opt
if (present(transb_opt)) transb = transb_opt
! m ... number of rows in Op(A) and C
! n ... number of columns in Op(B) and C
! k ... number of columns in Op(A) resp. rows in Op(B)
m = size(c, 1)
n = size(c, 2)
if (transa /= 'N') then
k = size(a, 1)
else
k = size(a, 2)
end if
call zgemm(transa, transb, m, n, k, cmplx_1, a, size(a, 1), b, size(b, 1), cmplx_0, c, m)
end subroutine utility_zgemm_new