subroutine tran_transfer(tot, tott, h_00, h_01, e_scan_cmp, nxx, stdout, error, comm)
!================================================!
!
! iterative construction of the transfer matrix
! as Lopez-Sancho^2&Rubio, J.Phys.F:Met.Phys., v.14, 1205 (1984)
! and ibid. v.15, 851 (1985)
!
!================================================
use w90_constants, only: dp, cmplx_0, cmplx_1, eps7
use w90_error, only: w90_error_type, set_error_alloc, set_error_fatal, set_error_unconv, &
set_error_dealloc
implicit none
type(w90_error_type), allocatable, intent(out) :: error
type(w90_comm_type), intent(in) :: comm
integer, intent(in) :: nxx
integer, intent(in) :: stdout
complex(kind=dp), intent(in) :: e_scan_cmp
complex(kind=dp), intent(out) :: tot(nxx, nxx)
complex(kind=dp), intent(out) :: tott(nxx, nxx)
real(kind=dp), intent(in) :: h_00(nxx, nxx)
real(kind=dp), intent(in) :: h_01(nxx, nxx)
integer :: ierr, info
integer :: i, j, n, nxx2
integer, allocatable :: ipiv(:)
real(kind=dp) :: conver, conver2
complex(kind=dp), allocatable, dimension(:, :) :: tsum, tsumt
complex(kind=dp), allocatable, dimension(:, :) :: t11, t12
complex(kind=dp), allocatable, dimension(:, :) :: s1, s2
complex(kind=dp), allocatable, dimension(:, :, :) :: tau, taut
allocate (ipiv(nxx), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating ipiv in tran_transfer', comm)
return
end if
allocate (tsum(nxx, nxx), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating tsum in tran_transfer', comm)
return
end if
allocate (tsumt(nxx, nxx), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating tsumt in tran_transfer', comm)
return
end if
allocate (t11(nxx, nxx), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating t11 in tran_transfer', comm)
return
end if
allocate (t12(nxx, nxx), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating t12 in tran_transfer', comm)
return
end if
allocate (s1(nxx, nxx), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating s1 in tran_transfer', comm)
return
end if
allocate (s2(nxx, nxx), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating s2 in tran_transfer', comm)
return
end if
allocate (tau(nxx, nxx, 2), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating tau in tran_transfer', comm)
return
end if
allocate (taut(nxx, nxx, 2), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating taut in tran_transfer', comm)
return
end if
nxx2 = nxx*nxx
tot = cmplx_0
tott = cmplx_0
! construction of the transfer matrix
! t12 = e - h_00
t12(:, :) = cmplx(-h_00(:, :), kind=dp)
do i = 1, nxx
t12(i, i) = e_scan_cmp + t12(i, i)
end do
! compute (e - h_00)^-1 and store it in t11
t11 = cmplx_0
do i = 1, nxx
t11(i, i) = cmplx_1
end do
! inverse of t12 -> t11
call ZGESV(nxx, nxx, t12, nxx, ipiv, t11, nxx, info)
if (info .ne. 0) then
write (stdout, *) 'ERROR: IN ZGESV IN tran_transfer, INFO=', info
call set_error_fatal(error, 'tran_transfer: problem in ZGESV 1', comm)
return
end if
! compute intermediate t-matrices (defined as tau(nxx,nxx,niter)
! and taut(...)):
tau = cmplx_0
taut = cmplx_0
! t_0:
t12(:, :) = cmplx(h_01(:, :), kind=dp)
! tau = ( e - H_00 )^-1 * H_01^+
call ZGEMM('N', 'C', nxx, nxx, nxx, cmplx_1, t11, nxx, t12, nxx, cmplx_0, tau(1, 1, 1), nxx)
! taut = ( e - H_00 )^-1 * H_01
call ZGEMM('N', 'N', nxx, nxx, nxx, cmplx_1, t11, nxx, t12, nxx, cmplx_0, taut(1, 1, 1), nxx)
! initialize T:
tot(:, :) = tau(:, :, 1)
tsum(:, :) = taut(:, :, 1)
! initialize T^bar:
tott(:, :) = taut(:, :, 1)
tsumt(:, :) = tau(:, :, 1)
! main loop:
do n = 1, nterx
t11 = cmplx_0
t12 = cmplx_0
call ZGEMM('N', 'N', nxx, nxx, nxx, cmplx_1, tau(1, 1, 1), nxx, taut(1, 1, 1), nxx, cmplx_0, &
t11, nxx)
call ZGEMM('N', 'N', nxx, nxx, nxx, cmplx_1, taut(1, 1, 1), nxx, tau(1, 1, 1), nxx, cmplx_0, &
t12, nxx)
s1(:, :) = -t11(:, :) - t12(:, :)
do i = 1, nxx
s1(i, i) = cmplx_1 + s1(i, i)
end do
s2 = cmplx_0
do i = 1, nxx
s2(i, i) = cmplx_1
end do
call ZGESV(nxx, nxx, s1, nxx, ipiv, s2, nxx, info)
if (info .ne. 0) then
write (stdout, *) 'ERROR: IN ZGESV IN tran_transfer, INFO=', info
call set_error_fatal(error, 'tran_transfer: problem in ZGESV 2', comm)
return
end if
t11 = cmplx_0
t12 = cmplx_0
call ZGEMM('N', 'N', nxx, nxx, nxx, cmplx_1, tau(1, 1, 1), nxx, tau(1, 1, 1), nxx, cmplx_0, &
t11, nxx)
call ZGEMM('N', 'N', nxx, nxx, nxx, cmplx_1, taut(1, 1, 1), nxx, taut(1, 1, 1), nxx, &
cmplx_0, t12, nxx)
call ZGEMM('N', 'N', nxx, nxx, nxx, cmplx_1, s2, nxx, t11, nxx, cmplx_0, tau(1, 1, 2), nxx)
call ZGEMM('N', 'N', nxx, nxx, nxx, cmplx_1, s2, nxx, t12, nxx, cmplx_0, taut(1, 1, 2), nxx)
! put the transfer matrices together
t11 = cmplx_0
s1 = cmplx_0
call ZGEMM('N', 'N', nxx, nxx, nxx, cmplx_1, tsum, nxx, tau(1, 1, 2), nxx, cmplx_0, t11, nxx)
call ZGEMM('N', 'N', nxx, nxx, nxx, cmplx_1, tsum, nxx, taut(1, 1, 2), nxx, cmplx_0, s1, nxx)
call ZCOPY(nxx2, t11, 1, s2, 1)
call ZAXPY(nxx2, cmplx_1, tot, 1, s2, 1)
tot(:, :) = s2(:, :)
tsum(:, :) = s1(:, :)
t11 = cmplx_0
s1 = cmplx_0
call ZGEMM('N', 'N', nxx, nxx, nxx, cmplx_1, tsumt, nxx, taut(1, 1, 2), nxx, cmplx_0, t11, &
nxx)
call ZGEMM('N', 'N', nxx, nxx, nxx, cmplx_1, tsumt, nxx, tau(1, 1, 2), nxx, cmplx_0, s1, nxx)
call ZCOPY(nxx2, t11, 1, s2, 1)
call ZAXPY(nxx2, cmplx_1, tott, 1, s2, 1)
tott(:, :) = s2(:, :)
tsumt(:, :) = s1(:, :)
tau(:, :, 1) = tau(:, :, 2)
taut(:, :, 1) = taut(:, :, 2)
! convergence check on the t-matrices
conver = 0.0_dp
conver2 = 0.0_dp
do j = 1, nxx
do i = 1, nxx
conver = conver + sqrt(real(tau(i, j, 2), dp)**2 + aimag(tau(i, j, 2))**2)
conver2 = conver2 + sqrt(real(taut(i, j, 2), dp)**2 + aimag(taut(i, j, 2))**2)
end do
end do
if (conver .lt. eps7 .and. conver2 .lt. eps7) return
end do
if (conver .gt. eps7 .or. conver2 .gt. eps7) then
call set_error_unconv(error, 'Error in converging transfer matrix in tran_transfer', comm)
return
end if
deallocate (taut, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating taut in tran_transfer', comm)
return
end if
deallocate (tau, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating tau in tran_transfer', comm)
return
end if
deallocate (s2, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating s2 in tran_transfer', comm)
return
end if
deallocate (s1, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating s1 in tran_transfer', comm)
return
end if
deallocate (t12, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating t12 in tran_transfer', comm)
return
end if
deallocate (t11, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating t11 in tran_transfer', comm)
return
end if
deallocate (tsumt, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating tsumt in tran_transfer', comm)
return
end if
deallocate (tsum, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating tsum in tran_transfer', comm)
return
end if
deallocate (ipiv, stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in deallocating ipiv in tran_transfer', comm)
return
end if
return
end subroutine tran_transfer
