! { dg-do compile }
! { dg-options "-std=gnu" }
!
! Minimal implementation of the BLAS ?gemm subroutines for testing
! purposes. Very inefficient, but very straightforward.
subroutine sgemm (transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
  character(len=1) :: transa, transb
  integer :: m, n, k, lda, ldb, ldc
  real :: alpha, beta, a(lda,*), b(ldb,*), c(ldc,*)

  integer :: i, j
  logical :: nota, notb

  nota  = (transa == "N" .or. transa == "N")
  notb  = (transb == "N" .or. transb == "N")

  if (alpha /= 1 .or. beta /= 0) call abort
  if (all(transa /= (/"N", "n", "T", "t"/)) .or. &
      all(transa /= (/"N", "n", "T", "t"/))) call abort
  if (any((/m, n, k/) < 0)) call abort
  if ((nota .and. lda < max(1,m)) .or. ((.not. nota) .and. lda < max(1,k))) &
    call abort
  if ((notb .and. ldb < max(1,k)) .or. ((.not. notb) .and. ldb < max(1,n))) &
    call abort
  if (ldc < max(1,m)) call abort
  if (any((/m, n/) == 0)) return

  do i = 1, m
    do j = 1, n
      if (notb) then
        if (nota) then
          c(i,j) = sum(a(i,1:k)*b(1:k,j))
        else
          c(i,j) = sum(a(1:k,i)*b(1:k,j))
        end if
      else
        if (nota) then
          c(i,j) = sum(a(i,1:k)*b(j,1:k))
        else
          c(i,j) = sum(a(1:k,i)*b(j,1:k))
        end if
      end if
    end do
  end do
end subroutine


subroutine dgemm (transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
  character(len=1) :: transa, transb
  integer :: m, n, k, lda, ldb, ldc
  double precision :: alpha, beta, a(lda,*), b(ldb,*), c(ldc,*)

  integer :: i, j
  logical :: nota, notb

  nota  = (transa == "N" .or. transa == "N")
  notb  = (transb == "N" .or. transb == "N")

  if (alpha /= 1 .or. beta /= 0) call abort
  if (all(transa /= (/"N", "n", "T", "t"/)) .or. &
      all(transa /= (/"N", "n", "T", "t"/))) call abort
  if (any((/m, n, k/) < 0)) call abort
  if ((nota .and. lda < max(1,m)) .or. ((.not. nota) .and. lda < max(1,k))) &
    call abort
  if ((notb .and. ldb < max(1,k)) .or. ((.not. notb) .and. ldb < max(1,n))) &
    call abort
  if (ldc < max(1,m)) call abort
  if (any((/m, n/) == 0)) return

  do i = 1, m
    do j = 1, n
      if (notb) then
        if (nota) then
          c(i,j) = sum(a(i,1:k)*b(1:k,j))
        else
          c(i,j) = sum(a(1:k,i)*b(1:k,j))
        end if
      else
        if (nota) then
          c(i,j) = sum(a(i,1:k)*b(j,1:k))
        else
          c(i,j) = sum(a(1:k,i)*b(j,1:k))
        end if
      end if
    end do
  end do
end subroutine


subroutine cgemm (transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
  character(len=1) :: transa, transb
  integer :: m, n, k, lda, ldb, ldc
  complex :: alpha, beta, a(lda,*), b(ldb,*), c(ldc,*)

  integer :: i, j
  logical :: nota, notb

  nota  = (transa == "N" .or. transa == "N")
  notb  = (transb == "N" .or. transb == "N")

  if (alpha /= 1 .or. beta /= 0) call abort
  if (all(transa /= (/"N", "n", "T", "t"/)) .or. &
      all(transa /= (/"N", "n", "T", "t"/))) call abort
  if (any((/m, n, k/) < 0)) call abort
  if ((nota .and. lda < max(1,m)) .or. ((.not. nota) .and. lda < max(1,k))) &
    call abort
  if ((notb .and. ldb < max(1,k)) .or. ((.not. notb) .and. ldb < max(1,n))) &
    call abort
  if (ldc < max(1,m)) call abort
  if (any((/m, n/) == 0)) return

  do i = 1, m
    do j = 1, n
      if (notb) then
        if (nota) then
          c(i,j) = sum(a(i,1:k)*b(1:k,j))
        else
          c(i,j) = sum(a(1:k,i)*b(1:k,j))
        end if
      else
        if (nota) then
          c(i,j) = sum(a(i,1:k)*b(j,1:k))
        else
          c(i,j) = sum(a(1:k,i)*b(j,1:k))
        end if
      end if
    end do
  end do
end subroutine


subroutine zgemm (transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
  character(len=1) :: transa, transb
  integer :: m, n, k, lda, ldb, ldc
  double complex :: alpha, beta, a(lda,*), b(ldb,*), c(ldc,*)

  integer :: i, j
  logical :: nota, notb

  nota  = (transa == "N" .or. transa == "N")
  notb  = (transb == "N" .or. transb == "N")

  if (alpha /= 1 .or. beta /= 0) call abort
  if (all(transa /= (/"N", "n", "T", "t"/)) .or. &
      all(transa /= (/"N", "n", "T", "t"/))) call abort
  if (any((/m, n, k/) < 0)) call abort
  if ((nota .and. lda < max(1,m)) .or. ((.not. nota) .and. lda < max(1,k))) &
    call abort
  if ((notb .and. ldb < max(1,k)) .or. ((.not. notb) .and. ldb < max(1,n))) &
    call abort
  if (ldc < max(1,m)) call abort
  if (any((/m, n/) == 0)) return

  do i = 1, m
    do j = 1, n
      if (notb) then
        if (nota) then
          c(i,j) = sum(a(i,1:k)*b(1:k,j))
        else
          c(i,j) = sum(a(1:k,i)*b(1:k,j))
        end if
      else
        if (nota) then
          c(i,j) = sum(a(i,1:k)*b(j,1:k))
        else
          c(i,j) = sum(a(1:k,i)*b(j,1:k))
        end if
      end if
    end do
  end do
end subroutine