lapack zheev example program

module matgen
    implicit none
    contains

    ! Creates pseudo-random positive-definite hermetian matrix
    function create_random_hermitian_matrix(N) result(A)
    complex(8), allocatable, dimension(:,:) :: A, temp
    real(8)                                 :: rv, iv
    integer                                 :: i, j, N

    if (.not. allocated(A))    allocate(A(N,N))
    if (.not. allocated(temp)) allocate(temp(N,N))

    ! Create general hermitian temp
    do j = 1, N
        do i = 1, N
            if (i > j) then
                call random_number(rv)
                call random_number(iv)
                temp(i,j) = cmplx(rv, iv, 8)
                temp(j,i) = conjg(temp(i,j))
            else if (i == j) then
                call random_number(rv)
                temp(i,j) = rv
            end if
        end do
    end do

    ! Multiply temp by conjugate transpose of temp to get positive definite hermitian A
    call zgemm('N', 'C', N, N, N, cmplx(1.0, 0.0, 8), temp, N, temp, N, cmplx(0.0, 0.0, 8), A, N)

    deallocate(temp)

    end function create_random_hermitian_matrix

end module matgen


program main
    use matgen
    implicit none

    integer                                 :: N, i, k, lda, istat
    integer                                 :: lwork, lrwork, liwork
    character(len=20)                       :: arg
    complex(8), dimension(:,:), allocatable :: A, Aref
    complex(8), dimension(:,:), allocatable :: B, Bref
    complex(8), dimension(:), allocatable   :: work
    real(8), dimension(:), allocatable      :: w, rwork
    integer, dimension(:), allocatable      :: iwork

    ! Parse command line arguments
    i = command_argument_count()

    if (i >= 1) then
        ! If N is provided, generate random hermetian matrices for A and B
        print*, "Using randomly-generated matrices..."
        call get_command_argument(1, arg)
        read(arg, *) N
        lda = N
    else
        print*, "Usage error: please provide at least one integer argument (matrix dimension)."
        stop
    endif

    print*, "Running with N = ", N

    ! size of temp. work arrays from lapack
    lwork = 2*N + N*N
    lrwork = 1 + 5*N + 2*N*N
    liwork = 3 + 5*N

    ! deallocate first (just in case)
    if (allocated(Aref))  deallocate(Aref)
    if (allocated(Bref))  deallocate(Bref)
    if (allocated(A))     deallocate(A)
    if (allocated(B))     deallocate(B)
    if (allocated(w))     deallocate(w)
    if (allocated(iwork)) deallocate(iwork)
    if (allocated(rwork)) deallocate(rwork)
    if (allocated(work))  deallocate(work)
    if (allocated(work))  deallocate(work)
    if (allocated(rwork)) deallocate(rwork)
    if (allocated(iwork)) deallocate(iwork)

    ! dynamic memory allocation
    if (.not. allocated(Aref))  allocate(Aref(N,N))
    if (.not. allocated(Bref))  allocate(Bref(N,N))
    if (.not. allocated(A))     allocate(A(N,N))
    if (.not. allocated(B))     allocate(B(N,N))
    if (.not. allocated(w))     allocate(w(N))
    if (.not. allocated(iwork)) allocate(iwork(liwork))
    if (.not. allocated(rwork)) allocate(rwork(lrwork))
    if (.not. allocated(work))  allocate(work(lwork))
    if (.not. allocated(work))  allocate(work(lwork))
    if (.not. allocated(rwork)) allocate(rwork(lrwork))
    if (.not. allocated(iwork)) allocate(iwork(liwork))


    ! Create random positive-definite hermetian matrices on host
    Aref = create_random_hermitian_matrix(N)

    ! Bref is the identity matrix
    Bref = 0
    forall(k = 1:N) Bref(k,k) = cmplx(1.0, 0.0, 8)    ! Set the diagonal

    A(:,:) = Aref(:,:)
    B(:,:) = Bref(:,:)

    !call zhegvd(1, 'V', 'U', N, A, lda, B, lda, w, work, lwork, rwork, lrwork, iwork, liwork, istat)
    !call zheevd('V', 'U', N, A, lda, w, work, lwork, rwork, lrwork, iwork, liwork, istat)
    call zheev( 'V', 'U', N, A, lda, w, work, lwork, rwork, istat)
    if (istat /= 0) write(*,*) 'zhexxx failed. istat = ', istat

    print*, 'First eigenvalue = ', w(1)
    print*, 'Last eigenvalue = ', w(N)


end program main