blas_matmul Function

    function blas_matmul(A, B) result(C)
        ! a basic wrapper to the most fundamental matrix mult with blas
        HElement_t(dp), intent(in) :: A(:,:), B(:,:)
        HElement_t(dp) :: C(size(A,1),size(A,2))

        integer :: n

        n = size(A,1)
#ifdef CMPLX_
        call zgemm('N','N', n, n, n, cmplx(1.0_dp,0.0_dp,kind=dp), A, n, B, n, &
            cmplx(1.0_dp,0.0_dp,kind=dp), C, n)
#else
        call dgemm('N', 'N', n, n, n, 1.0_dp, A, n, B, n, 0.0_dp, C, n)
#endif
    end function blas_matmul