program main

    character :: uplo
    integer :: n
    double complex :: alphah, alphal
    double complex, allocatable :: ah(:,:), al(:,:)
    integer :: lda
    double complex, allocatable :: xh(:), xl(:)
    integer :: incx
    double complex :: betah, betal
    double complex, allocatable :: yh(:), yl(:)
    integer :: incy

    integer :: i, j

    uplo = 'u'
    n = 3
    incx = 1
    incy = 1

    lda = max(1, n)
    allocate(ah(lda,n))
    allocate(al(lda,n))
    allocate(xh(1+(n-1)*abs(incx)))
    allocate(xl(1+(n-1)*abs(incx)))
    allocate(yh(1+(n-1)*abs(incy)))
    allocate(yl(1+(n-1)*abs(incy)))

    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(xh,1)
       call random_complex(xh(i))
    end do
    xl = 0d0

    do i = 1, size(xh,1)
       call random_complex(yh(i))
    end do
    yl = 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_vector("x", xh, xl)
    call print_scalar("beta", betah, betal)
    call print_vector("y", yh, yl)

    call gzzhemv_host_multigpu(uplo, n, alphah, alphal, ah, al, lda, xh, xl, incx, &
        & betah, betal, yh, yl, incy)

    print *, "OUTPUT"
    call print_vector("y", yh, yl)

    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_vector(vname, xh, xl)

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

        do i = 1, size(xh, 1)
            print '(a,"(",i0,") = (",e23.15,", ",e23.15,")")', &
                & vname//'h', i, dble(xh(i)), dimag(xh(i))
            print '(a,"(",i0,") = (",e23.15,", ",e23.15,")")', &
                & vname//'l', i, dble(xl(i)), dimag(xl(i))
        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