qp_plugins_scemama/fnmf/non_hermit.irp.f

subroutine lapack_g_non_sym_real(n, H, LDH, S, LDS, beta, &
      n_real, vl, LDVL, vr, LDVR)
  implicit none
  integer, intent(in) :: n, LDH, LDS, LDVL, LDVR
  double precision, intent(in) :: H(LDH,n), S(LDS,n)
  double precision, intent(out) :: VL(LDVL,n), VR(LDVR,n), beta(n)
  integer, intent(out) :: n_real

  integer :: lwork, info, i,j 
  double precision, allocatable :: work(:), Htmp(:,:), Stmp(:,:)
  double precision, allocatable :: alphar(:), alphai(:), vltmp(:,:), vrtmp(:,:)
  integer, allocatable :: iorder(:), list_good(:)


  lwork = -1
  allocate(work(1), alphar(n), alphai(n), vltmp(n,n), vrtmp(n,n), &
  Htmp(n,n), Stmp(n,n), list_good(n))

  Htmp(1:n,1:n) = H(1:n,1:n)
  Stmp(1:n,1:n) = S(1:n,1:n)
  call dggev('V', 'V', n, Htmp, size(Htmp,1), Stmp, size(Stmp,1), alphar, alphai, beta, &
             vltmp, size(vltmp,1), vrtmp, size(vrtmp,1), work, lwork, info)

  lwork = int(work(1))
  deallocate(work)
  allocate(work(lwork))
  call dggev('V', 'V', n, Htmp, size(Htmp,1), Stmp, size(Stmp,1), alphar, alphai, beta, &
             vltmp, size(vltmp,1), vrtmp, size(vrtmp,1), work, lwork, info)

  deallocate(work)
  if (info /= 0) then
    stop 'DGGEV Diagonalization failed'
  endif

  allocate(iorder(n))
  n_real = 0
  do i=1,n
    if (dabs(alphai(i)) < 1.d-10) then 
       n_real += 1
       list_good(n_real) = i
    else
       alphar(i) = dble(huge(1.0))
    endif
  enddo

  do i=1,n
    if (dabs(beta(i)/(alphar(i)+1.d-12)) < 1.d-10) then
      beta(i) = 0.d0
    else
      beta(i) = alphar(i)/beta(i)
    endif
  enddo

  do i=1,n_real
    iorder(i) = i
    do j=1,n
      vrtmp(j,i) = vrtmp(j,list_good(i))
      vltmp(j,i) = vltmp(j,list_good(i))
    end do
  end do

  call dsort(beta, iorder, n_real)

  do i=1,n_real
    do j=1,n
      vr(j,i) = vrtmp(j,iorder(i))
      vl(j,i) = vltmp(j,iorder(i))
    end do
  end do

  deallocate(vrtmp, vltmp, iorder, Htmp, Stmp)

end