program main

    character :: side, uplo
    integer :: m, n
    double complex :: alphah, alphal
    double complex, allocatable :: ah(:,:), al(:,:)
    integer :: lda
    double complex, allocatable :: bh(:,:), bl(:,:)
    integer :: ldb
    double complex :: betah, betal
    double complex, allocatable :: ch(:,:), cl(:,:)
    integer :: ldc

    integer :: i, j

    side = 'l'
    uplo = 'u'
    m = 3
    n = 2

    lda = max(1, m)
    ldb = max(1, m)
    ldc = max(1, m)
    allocate(ah(lda,m))
    allocate(al(lda,m))
    allocate(bh(ldb,n))
    allocate(bl(ldb,n))
    allocate(ch(ldc,n))
    allocate(cl(ldc,n))

    do i = 1, size(ah, 1)
        do j = i, size(ah, 2)
            call random_complex(ah(i,j))
        end do
        ah(i,i) = dble(ah(i,i))
    end do
    al = 0d0

    do i = 1, size(bh,1)
       do j = 1, size(bh, 2)
           call random_complex(bh(i,j))
       end do
    end do
    bl = 0d0

    do i = 1, size(ch,1)
       do j = 1, size(ch, 2)
           call random_complex(ch(i,j))
       end do
    end do
    cl = 0d0

    call random_complex(alphah)
    alphal = 0d0

    call random_complex(betah)
    betal = 0d0

    print *, "INPUT"
    call print_scalar("alpha", alphah, alphal)
    call print_matrix("a", ah, al)
    call print_matrix("b", bh, bl)
    call print_scalar("beta", betah, betal)
    call print_matrix("c", ch, cl)

    call zzhemm(side, uplo, m, n, alphah, alphal, ah, al, lda, &
        & bh, bl, ldb, betah, betal, ch, cl, ldc)

    print *, "OUTPUT"
    call print_matrix("c", ch, cl)

    contains

    subroutine random_complex(x)

        double complex, intent(out) :: x
        double precision :: re, im

        call random_number(re)
        call random_number(im)
        x = dcmplx(re, im)

    end subroutine

    subroutine print_matrix(vname, ah, al)

        character(*), intent(in) :: vname
        double complex, intent(in) :: ah(:,:), al(:,:)
        integer :: i, j

        do i = 1, size(ah, 1)
            do j = 1, size(ah, 2)
                print '(a,"(",i0,",",i0,") = (",e23.15,", ",e23.15,")")', &
                    & vname//'h', i, j, dble(ah(i,j)), dimag(ah(i,j))
                print '(a,"(",i0,",",i0,") = (",e23.15,", ",e23.15,")")', &
                    & vname//'l', i, j, dble(al(i,j)), dimag(al(i,j))
            end do
        end do

    end subroutine

    subroutine print_scalar(vname, xh, xl)

        character(*), intent(in) :: vname
        double complex, intent(in) :: xh, xl

        print '(a," = (",e23.15,", ",e23.15,")")', &
            & vname//'h', dble(xh), dimag(xh)
        print '(a," = (",e23.15,", ",e23.15,")")', &
            & vname//'l', dble(xl), dimag(xl)

    end subroutine

end program