eig_cmplx Subroutine

private subroutine eig_cmplx(matrix, e_values, e_vectors, t_left_ev)

Arguments

Type IntentOptional Attributes Name
complex(kind=dp), intent(in) :: matrix(:,:)
real(kind=dp), intent(out) :: e_values(size(matrix,1))
complex(kind=dp), intent(out), optional :: e_vectors(size(matrix,1),size(matrix,1))
logical, intent(in), optional :: t_left_ev

Contents

Source Code


Source Code

    subroutine eig_cmplx(matrix, e_values, e_vectors, t_left_ev)
        ! for very restricted matrices do a diag routine here!
        complex(dp), intent(in) :: matrix(:,:)
        real(dp), intent(out) :: e_values(size(matrix,1))
        complex(dp), intent(out), optional :: e_vectors(size(matrix,1),size(matrix,1))
        logical, intent(in), optional :: t_left_ev
        character(*), parameter :: this_routine = 'eig_cmplx'

        ! get the specifics for the eigenvectors still..
        ! i think i need a bigger work, and maybe also a flag for how many
        ! eigenvectors i want.. maybe also the number of eigenvalues..
        integer :: n, info
        complex(dp) :: work(4*size(matrix,1)), tmp_matrix(size(matrix,1),size(matrix,2))
        complex(dp) :: left_ev(size(matrix,1),size(matrix,1))
        real(dp), allocatable :: rwork(:)
        real(dp) :: right_ev(size(matrix,1),size(matrix,1))
        integer :: sort_ind(size(matrix,1))
        character(len=1) :: left, right


        ! and convention is: we only want the right eigenvectors!!
        ! and always assume real-only eigenvalues
        if (present(e_vectors)) then

            if (present(t_left_ev)) then
                if (t_left_ev) then
                    left = 'V'
                    right = 'N'
                else
                    left = 'N'
                    right = 'V'
                end if
            else
                left = 'N'
                right = 'V'
            end if

            n = size(matrix,1)

            tmp_matrix = matrix

            left = 'V'
            right = 'V'


            allocate(rwork(max(1,3*n-2)))
            call zheev(&
                 left, &
                 right, &
                 n, &
                 tmp_matrix, &
                 n, &
                 e_values, &
                 work, &
                 4*n, &
                 rwork, &
                 info)
            if (info /= 0) call stop_all(this_routine, 'Failed in BLAS call.')
            deallocate(rwork)

            sort_ind = [(n, n = 1, size(matrix,1))]

            call sort(e_values, sort_ind)

            if (present(t_left_ev)) then
                if (t_left_ev) then
                    e_vectors = left_ev(:,sort_ind)
                else
                    e_vectors = right_ev(:,sort_ind)
                end if
            else
                e_vectors = right_ev(:,sort_ind)
            end if

        else
            e_values = calc_eigenvalues(matrix)
        end if

    end subroutine eig_cmplx