"""
Fortran module KrylovOps: August 2017 (dj)
Contains the subroutines for propagating globally or locally the
wave-function with the Krylov algorithm.
**Authors**
* D. Jaschke
* M. L. Wall
**Details**
Contains algorithms for Lanczos approach (Hamiltonians) and
Arnoldi (Liouville operator for Lindblad master equation).
The following subroutines / functions are defined for the
applicable data type
+--------------------+-------------+---------+
| procedure | include.f90 | mpi.f90 |
+====================+=============+=========+
| krylov_method | X | |
+--------------------+-------------+---------+
"""
[docs]def krylov_method_mpo_complex():
"""
fortran-subroutine - August 2017 (dj, updated)
Propagate the wave function Psi with the Krylov/Lanczos algorithm.
**Arguments**
Psi : TYPE(mpsc), inout
This is the wave function to be propagated.
deltat : complex, in
The time step for the evolution.
Ham : TYPE(mpo), inout
Evolve the wave function under this Hamiltonian.
converged : LOGICAL, out
On exit, it indicates if the algorithm is converged within then
given tolerances.
cerr : REAL, inout
The cumulated error. The error done during this subroutine is
added to the incoming value.
Cp : TYPE(ConvParam), inout
Stores all the convergence parameters necessary for the algorithm.
init_random : CHARACTER, inout
Determines the initialization in the fit of H | psi >.
renorm : LOGICAL, in
Flag if state vector should be renormalized to 1 (true).
'N' : do not normalize (default); 'M' : normalize for MPS
1 / sqrt(norm);
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
Currently uses complete reorthogonalization for stability.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_method_mpo_complex(Psi, deltat, Ham, converged, &
cerr, Cp, init_random, renorm, pbc, errst)
type(mpsc), intent(inout) :: Psi
complex(KIND=rKind), intent(in) :: deltat
type(mpo), intent(inout) :: Ham
logical, intent(out) :: converged
real(KIND=rKind), intent(inout) :: cerr
type(ConvParam), intent(in) :: Cp
character, intent(in) :: init_random
character, intent(in) :: renorm
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: jj, kk, pp
! for looping / unused index
integer :: ii
! diagonal and 1st offdiagonal of tridiagonal upper hessenberg
real(KIND=rKind), allocatable :: beta(:), alpha(:)
! number of Krylov vectors and Krylov basis
integer :: nkrylov
type(mpsc), pointer :: Krylovuv(:)
! MPS to be modified in iteration process
type(mpsc) :: Rs
! the overlap between two states
complex(KIND=rKind) :: overlap
! error of a single call; added to cumulative error
real(KIND=rKind) :: err
! left right overlap
type(tensorlistc), pointer :: LR(:)
! The exponential in the Krylov basis, including the previous step
complex(KIND=rKind), allocatable :: prop(:, :)
complex(KIND=rKind), allocatable :: propsave(:)
! estimate of convergence based on propagator in Krylov space
real(KIND=rKind) :: estimate
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_method_mpo_complex:'//&
! 'Krylov does not have PBC.', 99, 'KrylovOps_include.f90:134', &
! errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylovuv(nkrylov), LR(nkrylov), alpha(nkrylov), &
beta(0:nkrylov), propsave(nkrylov))
propsave(1) = 10.0_rkind
call copy(Rs, Psi)
beta(0) = sqrt(norm(Rs))
! First iteration
! ---------------
jj = 1
call copy(Krylovuv(jj), Rs, scalar=zone / beta(jj - 1))
call destroy(Rs)
call fit_hpsi(Rs, Ham, Krylovuv(jj), Cp, err, converged, init_random, &
'V', 'N', errst=errst)
!if(prop_error('krylov_method_mpo: fit_hpsi (1) failed', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method: failed on jj =', &
jj, 'Hpsi with error', err
alpha(jj) = real(dot(Rs, Krylovuv(jj)), KIND=rKind)
! Setup LR overlaps from j to r - use them to get alpha
do kk = 1, jj
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc)
end do
! Orthogonalize
call orth(Rs, Krylovuv, 1, 1, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_method_mpo : orth (1) failed.', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method:, failed on jj =', &
jj, 'three-term recurrence with error', err
! Reorthogonalize
call orth(Rs, Krylovuv, 1, jj, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_method_mpo : orth (2) failed.', &
! errst=errst)) return
if(.not. converged) write(slog, *) 'krylov_method:, failed on jj =', &
jj, 'reorthogonalization', err
do kk = 1, jj
do pp = 1, size(LR(kk)%Li, 1)
call destroy(LR(kk)%Li(pp))
end do
deallocate(LR(kk)%Li)
end do
beta(jj) = sqrt(norm(Rs))
if(beta(1) < Cp%lanczos_tol) then
! Psi is an eigenstate of H to within the given tolerance
call scale(exp(deltat * alpha(1)), Psi%Aa(max(Psi%oc, 1)))
call destroy(Rs)
call destroy(Krylovuv(1))
deallocate(Krylovuv, alpha, beta, propsave)
return
end if
call destroy(Psi)
! Further iterations to build Krylov basis
! ----------------------------------------
do jj = 2, nkrylov
call copy(Krylovuv(jj), Rs, scalar=zone / beta(jj - 1))
call destroy(Rs)
call fit_hpsi(Rs, Ham, Krylovuv(jj), Cp, err, converged, &
init_random, 'V', 'N', errst=errst)
!if(prop_error('krylov_method_mpo: fit_hpsi (2) '//&
! 'failed', errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method: faild on jj =', &
jj, 'Hpsi with error', err
alpha(jj) = real(dot(Rs, Krylovuv(jj)), KIND=rKind)
do kk = (jj - 1), jj
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc)
end do
! Three-term recurrence
call orth(Rs, Krylovuv, (jj - 1), jj, LR, Cp, err, converged, &
errst=errst)
!if(prop_error('krylov_method_mpo : orth (3) failed.', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method:, failed '//&
'on jj =', jj, 'three-term recurrence with error', err
do kk = 1, (jj - 2)
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc)
end do
! Reorthogonalize
call orth(Rs, Krylovuv, 1, jj, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_method_mpo : orth (4) failed.', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method:, failed '//&
'on jj =', jj, 'reorthogonalization', err
do kk = 1, jj
do pp = 1, size(LR(kk)%Li, 1)
call destroy(LR(kk)%Li(pp))
end do
deallocate(LR(kk)%Li)
end do
beta(jj) = sqrt(norm(Rs))
allocate(prop(jj, jj))
call expmh(prop, deltat, alpha(1:jj), beta(1:jj - 1), jj, &
errst=errst)
!if(prop_error('krylov_method_MPOM_TYPE_complex: '//&
! 'expmh failed.', 'KrylovOps_include.f90:274', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * beta(jj - 1) * abs(prop(jj, 1)))
if(abs(estimate) <= Cp%lanczos_tol) exit
propsave(1:jj - 1) = prop(1:jj - 1, 1) - propsave(1:jj - 1)
estimate = real(dot_product(propsave(1:jj - 1) , &
propsave(1:jj - 1)), KIND=rKind)
if(estimate <= Cp%psi_tol) exit
propsave(1:jj) = prop(:, 1)
! Avoid increment to jj = nkrylov + 1
if(jj == nkrylov) exit
deallocate(prop)
end do
! To-Do : seems to be repetitive
call expmh(prop, deltat, alpha(1:jj), beta(1:jj - 1), jj, errst=errst)
!if(prop_error('krylov_method_MPOM_TYPE_complex: '//&
! 'expmh failed.', 'KrylovOps_include.f90:298', &
! errst=errst)) then
! return
!end if
call destroy(Rs)
! location of largest weight
kk = maxloc(abs(prop(1:jj, 1)), 1)
! Initialize variational ansatz to be the largest weight MPS
call copy(Psi, Krylovuv(kk))
call fitmps(Psi, prop(1:jj, 1), Krylovuv, Cp, err, converged)
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method: failed on '//&
'fitting with error', err
do ii = 1, jj
call destroy(Krylovuv(ii))
end do
deallocate(Krylovuv, LR, alpha, beta, prop, propsave)
end subroutine krylov_method_mpo_complex
"""
return
[docs]def krylov_method_mpoc_complex():
"""
fortran-subroutine - August 2017 (dj, updated)
Propagate the wave function Psi with the Krylov/Lanczos algorithm.
**Arguments**
Psi : TYPE(mpsc), inout
This is the wave function to be propagated.
deltat : complex, in
The time step for the evolution.
Ham : TYPE(mpoc), inout
Evolve the wave function under this Hamiltonian.
converged : LOGICAL, out
On exit, it indicates if the algorithm is converged within then
given tolerances.
cerr : REAL, inout
The cumulated error. The error done during this subroutine is
added to the incoming value.
Cp : TYPE(ConvParam), inout
Stores all the convergence parameters necessary for the algorithm.
init_random : CHARACTER, inout
Determines the initialization in the fit of H | psi >.
renorm : LOGICAL, in
Flag if state vector should be renormalized to 1 (true).
'N' : do not normalize (default); 'M' : normalize for MPS
1 / sqrt(norm);
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
Currently uses complete reorthogonalization for stability.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_method_mpoc_complex(Psi, deltat, Ham, converged, &
cerr, Cp, init_random, renorm, pbc, errst)
type(mpsc), intent(inout) :: Psi
complex(KIND=rKind), intent(in) :: deltat
type(mpoc), intent(inout) :: Ham
logical, intent(out) :: converged
real(KIND=rKind), intent(inout) :: cerr
type(ConvParam), intent(in) :: Cp
character, intent(in) :: init_random
character, intent(in) :: renorm
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: jj, kk, pp
! for looping / unused index
integer :: ii
! diagonal and 1st offdiagonal of tridiagonal upper hessenberg
real(KIND=rKind), allocatable :: beta(:), alpha(:)
! number of Krylov vectors and Krylov basis
integer :: nkrylov
type(mpsc), pointer :: Krylovuv(:)
! MPS to be modified in iteration process
type(mpsc) :: Rs
! the overlap between two states
complex(KIND=rKind) :: overlap
! error of a single call; added to cumulative error
real(KIND=rKind) :: err
! left right overlap
type(tensorlistc), pointer :: LR(:)
! The exponential in the Krylov basis, including the previous step
complex(KIND=rKind), allocatable :: prop(:, :)
complex(KIND=rKind), allocatable :: propsave(:)
! estimate of convergence based on propagator in Krylov space
real(KIND=rKind) :: estimate
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_method_mpoc_complex:'//&
! 'Krylov does not have PBC.', 99, 'KrylovOps_include.f90:134', &
! errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylovuv(nkrylov), LR(nkrylov), alpha(nkrylov), &
beta(0:nkrylov), propsave(nkrylov))
propsave(1) = 10.0_rkind
call copy(Rs, Psi)
beta(0) = sqrt(norm(Rs))
! First iteration
! ---------------
jj = 1
call copy(Krylovuv(jj), Rs, scalar=zone / beta(jj - 1))
call destroy(Rs)
call fit_hpsi(Rs, Ham, Krylovuv(jj), Cp, err, converged, init_random, &
'V', 'N', errst=errst)
!if(prop_error('krylov_method_mpoc: fit_hpsi (1) failed', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method: failed on jj =', &
jj, 'Hpsi with error', err
alpha(jj) = real(dot(Rs, Krylovuv(jj)), KIND=rKind)
! Setup LR overlaps from j to r - use them to get alpha
do kk = 1, jj
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc)
end do
! Orthogonalize
call orth(Rs, Krylovuv, 1, 1, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_method_mpoc : orth (1) failed.', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method:, failed on jj =', &
jj, 'three-term recurrence with error', err
! Reorthogonalize
call orth(Rs, Krylovuv, 1, jj, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_method_mpoc : orth (2) failed.', &
! errst=errst)) return
if(.not. converged) write(slog, *) 'krylov_method:, failed on jj =', &
jj, 'reorthogonalization', err
do kk = 1, jj
do pp = 1, size(LR(kk)%Li, 1)
call destroy(LR(kk)%Li(pp))
end do
deallocate(LR(kk)%Li)
end do
beta(jj) = sqrt(norm(Rs))
if(beta(1) < Cp%lanczos_tol) then
! Psi is an eigenstate of H to within the given tolerance
call scale(exp(deltat * alpha(1)), Psi%Aa(max(Psi%oc, 1)))
call destroy(Rs)
call destroy(Krylovuv(1))
deallocate(Krylovuv, alpha, beta, propsave)
return
end if
call destroy(Psi)
! Further iterations to build Krylov basis
! ----------------------------------------
do jj = 2, nkrylov
call copy(Krylovuv(jj), Rs, scalar=zone / beta(jj - 1))
call destroy(Rs)
call fit_hpsi(Rs, Ham, Krylovuv(jj), Cp, err, converged, &
init_random, 'V', 'N', errst=errst)
!if(prop_error('krylov_method_mpoc: fit_hpsi (2) '//&
! 'failed', errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method: faild on jj =', &
jj, 'Hpsi with error', err
alpha(jj) = real(dot(Rs, Krylovuv(jj)), KIND=rKind)
do kk = (jj - 1), jj
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc)
end do
! Three-term recurrence
call orth(Rs, Krylovuv, (jj - 1), jj, LR, Cp, err, converged, &
errst=errst)
!if(prop_error('krylov_method_mpoc : orth (3) failed.', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method:, failed '//&
'on jj =', jj, 'three-term recurrence with error', err
do kk = 1, (jj - 2)
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc)
end do
! Reorthogonalize
call orth(Rs, Krylovuv, 1, jj, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_method_mpoc : orth (4) failed.', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method:, failed '//&
'on jj =', jj, 'reorthogonalization', err
do kk = 1, jj
do pp = 1, size(LR(kk)%Li, 1)
call destroy(LR(kk)%Li(pp))
end do
deallocate(LR(kk)%Li)
end do
beta(jj) = sqrt(norm(Rs))
allocate(prop(jj, jj))
call expmh(prop, deltat, alpha(1:jj), beta(1:jj - 1), jj, &
errst=errst)
!if(prop_error('krylov_method_MPOM_TYPE_complex: '//&
! 'expmh failed.', 'KrylovOps_include.f90:274', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * beta(jj - 1) * abs(prop(jj, 1)))
if(abs(estimate) <= Cp%lanczos_tol) exit
propsave(1:jj - 1) = prop(1:jj - 1, 1) - propsave(1:jj - 1)
estimate = real(dot_product(propsave(1:jj - 1) , &
propsave(1:jj - 1)), KIND=rKind)
if(estimate <= Cp%psi_tol) exit
propsave(1:jj) = prop(:, 1)
! Avoid increment to jj = nkrylov + 1
if(jj == nkrylov) exit
deallocate(prop)
end do
! To-Do : seems to be repetitive
call expmh(prop, deltat, alpha(1:jj), beta(1:jj - 1), jj, errst=errst)
!if(prop_error('krylov_method_MPOM_TYPE_complex: '//&
! 'expmh failed.', 'KrylovOps_include.f90:298', &
! errst=errst)) then
! return
!end if
call destroy(Rs)
! location of largest weight
kk = maxloc(abs(prop(1:jj, 1)), 1)
! Initialize variational ansatz to be the largest weight MPS
call copy(Psi, Krylovuv(kk))
call fitmps(Psi, prop(1:jj, 1), Krylovuv, Cp, err, converged)
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method: failed on '//&
'fitting with error', err
do ii = 1, jj
call destroy(Krylovuv(ii))
end do
deallocate(Krylovuv, LR, alpha, beta, prop, propsave)
end subroutine krylov_method_mpoc_complex
"""
return
[docs]def krylov_method_qmpo_complex():
"""
fortran-subroutine - August 2017 (dj, updated)
Propagate the wave function Psi with the Krylov/Lanczos algorithm.
**Arguments**
Psi : TYPE(qmpsc), inout
This is the wave function to be propagated.
deltat : complex, in
The time step for the evolution.
Ham : TYPE(qmpo), inout
Evolve the wave function under this Hamiltonian.
converged : LOGICAL, out
On exit, it indicates if the algorithm is converged within then
given tolerances.
cerr : REAL, inout
The cumulated error. The error done during this subroutine is
added to the incoming value.
Cp : TYPE(ConvParam), inout
Stores all the convergence parameters necessary for the algorithm.
init_random : CHARACTER, inout
Determines the initialization in the fit of H | psi >.
renorm : LOGICAL, in
Flag if state vector should be renormalized to 1 (true).
'N' : do not normalize (default); 'M' : normalize for MPS
1 / sqrt(norm);
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
Currently uses complete reorthogonalization for stability.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_method_qmpo_complex(Psi, deltat, Ham, converged, &
cerr, Cp, init_random, renorm, pbc, errst)
type(qmpsc), intent(inout) :: Psi
complex(KIND=rKind), intent(in) :: deltat
type(qmpo), intent(inout) :: Ham
logical, intent(out) :: converged
real(KIND=rKind), intent(inout) :: cerr
type(ConvParam), intent(in) :: Cp
character, intent(in) :: init_random
character, intent(in) :: renorm
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: jj, kk, pp
! for looping / unused index
integer :: ii
! diagonal and 1st offdiagonal of tridiagonal upper hessenberg
real(KIND=rKind), allocatable :: beta(:), alpha(:)
! number of Krylov vectors and Krylov basis
integer :: nkrylov
type(qmpsc), pointer :: Krylovuv(:)
! MPS to be modified in iteration process
type(qmpsc) :: Rs
! the overlap between two states
complex(KIND=rKind) :: overlap
! error of a single call; added to cumulative error
real(KIND=rKind) :: err
! left right overlap
type(qtensorclist), pointer :: LR(:)
! The exponential in the Krylov basis, including the previous step
complex(KIND=rKind), allocatable :: prop(:, :)
complex(KIND=rKind), allocatable :: propsave(:)
! estimate of convergence based on propagator in Krylov space
real(KIND=rKind) :: estimate
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_method_qmpo_complex:'//&
! 'Krylov does not have PBC.', 99, 'KrylovOps_include.f90:134', &
! errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylovuv(nkrylov), LR(nkrylov), alpha(nkrylov), &
beta(0:nkrylov), propsave(nkrylov))
propsave(1) = 10.0_rkind
call copy(Rs, Psi)
beta(0) = sqrt(norm(Rs))
! First iteration
! ---------------
jj = 1
call copy(Krylovuv(jj), Rs, scalar=zone / beta(jj - 1))
call destroy(Rs)
call fit_hpsi(Rs, Ham, Krylovuv(jj), Cp, err, converged, init_random, &
'V', 'N', errst=errst)
!if(prop_error('krylov_method_qmpo: fit_hpsi (1) failed', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method: failed on jj =', &
jj, 'Hpsi with error', err
alpha(jj) = real(dot(Rs, Krylovuv(jj)), KIND=rKind)
! Setup LR overlaps from j to r - use them to get alpha
do kk = 1, jj
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc)
end do
! Orthogonalize
call orth(Rs, Krylovuv, 1, 1, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_method_qmpo : orth (1) failed.', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method:, failed on jj =', &
jj, 'three-term recurrence with error', err
! Reorthogonalize
call orth(Rs, Krylovuv, 1, jj, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_method_qmpo : orth (2) failed.', &
! errst=errst)) return
if(.not. converged) write(slog, *) 'krylov_method:, failed on jj =', &
jj, 'reorthogonalization', err
do kk = 1, jj
do pp = 1, size(LR(kk)%Li, 1)
call destroy(LR(kk)%Li(pp))
end do
deallocate(LR(kk)%Li)
end do
beta(jj) = sqrt(norm(Rs))
if(beta(1) < Cp%lanczos_tol) then
! Psi is an eigenstate of H to within the given tolerance
call scale(exp(deltat * alpha(1)), Psi%Aa(max(Psi%oc, 1)))
call destroy(Rs)
call destroy(Krylovuv(1))
deallocate(Krylovuv, alpha, beta, propsave)
return
end if
call destroy(Psi)
! Further iterations to build Krylov basis
! ----------------------------------------
do jj = 2, nkrylov
call copy(Krylovuv(jj), Rs, scalar=zone / beta(jj - 1))
call destroy(Rs)
call fit_hpsi(Rs, Ham, Krylovuv(jj), Cp, err, converged, &
init_random, 'V', 'N', errst=errst)
!if(prop_error('krylov_method_qmpo: fit_hpsi (2) '//&
! 'failed', errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method: faild on jj =', &
jj, 'Hpsi with error', err
alpha(jj) = real(dot(Rs, Krylovuv(jj)), KIND=rKind)
do kk = (jj - 1), jj
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc)
end do
! Three-term recurrence
call orth(Rs, Krylovuv, (jj - 1), jj, LR, Cp, err, converged, &
errst=errst)
!if(prop_error('krylov_method_qmpo : orth (3) failed.', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method:, failed '//&
'on jj =', jj, 'three-term recurrence with error', err
do kk = 1, (jj - 2)
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc)
end do
! Reorthogonalize
call orth(Rs, Krylovuv, 1, jj, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_method_qmpo : orth (4) failed.', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method:, failed '//&
'on jj =', jj, 'reorthogonalization', err
do kk = 1, jj
do pp = 1, size(LR(kk)%Li, 1)
call destroy(LR(kk)%Li(pp))
end do
deallocate(LR(kk)%Li)
end do
beta(jj) = sqrt(norm(Rs))
allocate(prop(jj, jj))
call expmh(prop, deltat, alpha(1:jj), beta(1:jj - 1), jj, &
errst=errst)
!if(prop_error('krylov_method_MPOM_TYPE_complex: '//&
! 'expmh failed.', 'KrylovOps_include.f90:274', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * beta(jj - 1) * abs(prop(jj, 1)))
if(abs(estimate) <= Cp%lanczos_tol) exit
propsave(1:jj - 1) = prop(1:jj - 1, 1) - propsave(1:jj - 1)
estimate = real(dot_product(propsave(1:jj - 1) , &
propsave(1:jj - 1)), KIND=rKind)
if(estimate <= Cp%psi_tol) exit
propsave(1:jj) = prop(:, 1)
! Avoid increment to jj = nkrylov + 1
if(jj == nkrylov) exit
deallocate(prop)
end do
! To-Do : seems to be repetitive
call expmh(prop, deltat, alpha(1:jj), beta(1:jj - 1), jj, errst=errst)
!if(prop_error('krylov_method_MPOM_TYPE_complex: '//&
! 'expmh failed.', 'KrylovOps_include.f90:298', &
! errst=errst)) then
! return
!end if
call destroy(Rs)
! location of largest weight
kk = maxloc(abs(prop(1:jj, 1)), 1)
! Initialize variational ansatz to be the largest weight MPS
call copy(Psi, Krylovuv(kk))
call fitmps(Psi, prop(1:jj, 1), Krylovuv, Cp, err, converged)
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method: failed on '//&
'fitting with error', err
do ii = 1, jj
call destroy(Krylovuv(ii))
end do
deallocate(Krylovuv, LR, alpha, beta, prop, propsave)
end subroutine krylov_method_qmpo_complex
"""
return
[docs]def krylov_method_qmpoc_complex():
"""
fortran-subroutine - August 2017 (dj, updated)
Propagate the wave function Psi with the Krylov/Lanczos algorithm.
**Arguments**
Psi : TYPE(qmpsc), inout
This is the wave function to be propagated.
deltat : complex, in
The time step for the evolution.
Ham : TYPE(qmpoc), inout
Evolve the wave function under this Hamiltonian.
converged : LOGICAL, out
On exit, it indicates if the algorithm is converged within then
given tolerances.
cerr : REAL, inout
The cumulated error. The error done during this subroutine is
added to the incoming value.
Cp : TYPE(ConvParam), inout
Stores all the convergence parameters necessary for the algorithm.
init_random : CHARACTER, inout
Determines the initialization in the fit of H | psi >.
renorm : LOGICAL, in
Flag if state vector should be renormalized to 1 (true).
'N' : do not normalize (default); 'M' : normalize for MPS
1 / sqrt(norm);
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
Currently uses complete reorthogonalization for stability.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_method_qmpoc_complex(Psi, deltat, Ham, converged, &
cerr, Cp, init_random, renorm, pbc, errst)
type(qmpsc), intent(inout) :: Psi
complex(KIND=rKind), intent(in) :: deltat
type(qmpoc), intent(inout) :: Ham
logical, intent(out) :: converged
real(KIND=rKind), intent(inout) :: cerr
type(ConvParam), intent(in) :: Cp
character, intent(in) :: init_random
character, intent(in) :: renorm
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: jj, kk, pp
! for looping / unused index
integer :: ii
! diagonal and 1st offdiagonal of tridiagonal upper hessenberg
real(KIND=rKind), allocatable :: beta(:), alpha(:)
! number of Krylov vectors and Krylov basis
integer :: nkrylov
type(qmpsc), pointer :: Krylovuv(:)
! MPS to be modified in iteration process
type(qmpsc) :: Rs
! the overlap between two states
complex(KIND=rKind) :: overlap
! error of a single call; added to cumulative error
real(KIND=rKind) :: err
! left right overlap
type(qtensorclist), pointer :: LR(:)
! The exponential in the Krylov basis, including the previous step
complex(KIND=rKind), allocatable :: prop(:, :)
complex(KIND=rKind), allocatable :: propsave(:)
! estimate of convergence based on propagator in Krylov space
real(KIND=rKind) :: estimate
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_method_qmpoc_complex:'//&
! 'Krylov does not have PBC.', 99, 'KrylovOps_include.f90:134', &
! errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylovuv(nkrylov), LR(nkrylov), alpha(nkrylov), &
beta(0:nkrylov), propsave(nkrylov))
propsave(1) = 10.0_rkind
call copy(Rs, Psi)
beta(0) = sqrt(norm(Rs))
! First iteration
! ---------------
jj = 1
call copy(Krylovuv(jj), Rs, scalar=zone / beta(jj - 1))
call destroy(Rs)
call fit_hpsi(Rs, Ham, Krylovuv(jj), Cp, err, converged, init_random, &
'V', 'N', errst=errst)
!if(prop_error('krylov_method_qmpoc: fit_hpsi (1) failed', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method: failed on jj =', &
jj, 'Hpsi with error', err
alpha(jj) = real(dot(Rs, Krylovuv(jj)), KIND=rKind)
! Setup LR overlaps from j to r - use them to get alpha
do kk = 1, jj
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc)
end do
! Orthogonalize
call orth(Rs, Krylovuv, 1, 1, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_method_qmpoc : orth (1) failed.', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method:, failed on jj =', &
jj, 'three-term recurrence with error', err
! Reorthogonalize
call orth(Rs, Krylovuv, 1, jj, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_method_qmpoc : orth (2) failed.', &
! errst=errst)) return
if(.not. converged) write(slog, *) 'krylov_method:, failed on jj =', &
jj, 'reorthogonalization', err
do kk = 1, jj
do pp = 1, size(LR(kk)%Li, 1)
call destroy(LR(kk)%Li(pp))
end do
deallocate(LR(kk)%Li)
end do
beta(jj) = sqrt(norm(Rs))
if(beta(1) < Cp%lanczos_tol) then
! Psi is an eigenstate of H to within the given tolerance
call scale(exp(deltat * alpha(1)), Psi%Aa(max(Psi%oc, 1)))
call destroy(Rs)
call destroy(Krylovuv(1))
deallocate(Krylovuv, alpha, beta, propsave)
return
end if
call destroy(Psi)
! Further iterations to build Krylov basis
! ----------------------------------------
do jj = 2, nkrylov
call copy(Krylovuv(jj), Rs, scalar=zone / beta(jj - 1))
call destroy(Rs)
call fit_hpsi(Rs, Ham, Krylovuv(jj), Cp, err, converged, &
init_random, 'V', 'N', errst=errst)
!if(prop_error('krylov_method_qmpoc: fit_hpsi (2) '//&
! 'failed', errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method: faild on jj =', &
jj, 'Hpsi with error', err
alpha(jj) = real(dot(Rs, Krylovuv(jj)), KIND=rKind)
do kk = (jj - 1), jj
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc)
end do
! Three-term recurrence
call orth(Rs, Krylovuv, (jj - 1), jj, LR, Cp, err, converged, &
errst=errst)
!if(prop_error('krylov_method_qmpoc : orth (3) failed.', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method:, failed '//&
'on jj =', jj, 'three-term recurrence with error', err
do kk = 1, (jj - 2)
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc)
end do
! Reorthogonalize
call orth(Rs, Krylovuv, 1, jj, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_method_qmpoc : orth (4) failed.', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method:, failed '//&
'on jj =', jj, 'reorthogonalization', err
do kk = 1, jj
do pp = 1, size(LR(kk)%Li, 1)
call destroy(LR(kk)%Li(pp))
end do
deallocate(LR(kk)%Li)
end do
beta(jj) = sqrt(norm(Rs))
allocate(prop(jj, jj))
call expmh(prop, deltat, alpha(1:jj), beta(1:jj - 1), jj, &
errst=errst)
!if(prop_error('krylov_method_MPOM_TYPE_complex: '//&
! 'expmh failed.', 'KrylovOps_include.f90:274', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * beta(jj - 1) * abs(prop(jj, 1)))
if(abs(estimate) <= Cp%lanczos_tol) exit
propsave(1:jj - 1) = prop(1:jj - 1, 1) - propsave(1:jj - 1)
estimate = real(dot_product(propsave(1:jj - 1) , &
propsave(1:jj - 1)), KIND=rKind)
if(estimate <= Cp%psi_tol) exit
propsave(1:jj) = prop(:, 1)
! Avoid increment to jj = nkrylov + 1
if(jj == nkrylov) exit
deallocate(prop)
end do
! To-Do : seems to be repetitive
call expmh(prop, deltat, alpha(1:jj), beta(1:jj - 1), jj, errst=errst)
!if(prop_error('krylov_method_MPOM_TYPE_complex: '//&
! 'expmh failed.', 'KrylovOps_include.f90:298', &
! errst=errst)) then
! return
!end if
call destroy(Rs)
! location of largest weight
kk = maxloc(abs(prop(1:jj, 1)), 1)
! Initialize variational ansatz to be the largest weight MPS
call copy(Psi, Krylovuv(kk))
call fitmps(Psi, prop(1:jj, 1), Krylovuv, Cp, err, converged)
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method: failed on '//&
'fitting with error', err
do ii = 1, jj
call destroy(Krylovuv(ii))
end do
deallocate(Krylovuv, LR, alpha, beta, prop, propsave)
end subroutine krylov_method_qmpoc_complex
"""
return
[docs]def krylov_method_mpo_real():
"""
fortran-subroutine - August 2017 (dj, updated)
Propagate the wave function Psi with the Krylov/Lanczos algorithm.
**Arguments**
Psi : TYPE(mps), inout
This is the wave function to be propagated.
deltat : real, in
The time step for the evolution.
Ham : TYPE(mpo), inout
Evolve the wave function under this Hamiltonian.
converged : LOGICAL, out
On exit, it indicates if the algorithm is converged within then
given tolerances.
cerr : REAL, inout
The cumulated error. The error done during this subroutine is
added to the incoming value.
Cp : TYPE(ConvParam), inout
Stores all the convergence parameters necessary for the algorithm.
init_random : CHARACTER, inout
Determines the initialization in the fit of H | psi >.
renorm : LOGICAL, in
Flag if state vector should be renormalized to 1 (true).
'N' : do not normalize (default); 'M' : normalize for MPS
1 / sqrt(norm);
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
Currently uses complete reorthogonalization for stability.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_method_mpo_real(Psi, deltat, Ham, converged, &
cerr, Cp, init_random, renorm, pbc, errst)
type(mps), intent(inout) :: Psi
real(KIND=rKind), intent(in) :: deltat
type(mpo), intent(inout) :: Ham
logical, intent(out) :: converged
real(KIND=rKind), intent(inout) :: cerr
type(ConvParam), intent(in) :: Cp
character, intent(in) :: init_random
character, intent(in) :: renorm
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: jj, kk, pp
! for looping / unused index
integer :: ii
! diagonal and 1st offdiagonal of tridiagonal upper hessenberg
real(KIND=rKind), allocatable :: beta(:), alpha(:)
! number of Krylov vectors and Krylov basis
integer :: nkrylov
type(mps), pointer :: Krylovuv(:)
! MPS to be modified in iteration process
type(mps) :: Rs
! the overlap between two states
real(KIND=rKind) :: overlap
! error of a single call; added to cumulative error
real(KIND=rKind) :: err
! left right overlap
type(tensorlist), pointer :: LR(:)
! The exponential in the Krylov basis, including the previous step
real(KIND=rKind), allocatable :: prop(:, :)
real(KIND=rKind), allocatable :: propsave(:)
! estimate of convergence based on propagator in Krylov space
real(KIND=rKind) :: estimate
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_method_mpo_real:'//&
! 'Krylov does not have PBC.', 99, 'KrylovOps_include.f90:134', &
! errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylovuv(nkrylov), LR(nkrylov), alpha(nkrylov), &
beta(0:nkrylov), propsave(nkrylov))
propsave(1) = 10.0_rkind
call copy(Rs, Psi)
beta(0) = sqrt(norm(Rs))
! First iteration
! ---------------
jj = 1
call copy(Krylovuv(jj), Rs, scalar=done / beta(jj - 1))
call destroy(Rs)
call fit_hpsi(Rs, Ham, Krylovuv(jj), Cp, err, converged, init_random, &
'V', 'N', errst=errst)
!if(prop_error('krylov_method_mpo: fit_hpsi (1) failed', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method: failed on jj =', &
jj, 'Hpsi with error', err
alpha(jj) = real(dot(Rs, Krylovuv(jj)), KIND=rKind)
! Setup LR overlaps from j to r - use them to get alpha
do kk = 1, jj
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc)
end do
! Orthogonalize
call orth(Rs, Krylovuv, 1, 1, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_method_mpo : orth (1) failed.', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method:, failed on jj =', &
jj, 'three-term recurrence with error', err
! Reorthogonalize
call orth(Rs, Krylovuv, 1, jj, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_method_mpo : orth (2) failed.', &
! errst=errst)) return
if(.not. converged) write(slog, *) 'krylov_method:, failed on jj =', &
jj, 'reorthogonalization', err
do kk = 1, jj
do pp = 1, size(LR(kk)%Li, 1)
call destroy(LR(kk)%Li(pp))
end do
deallocate(LR(kk)%Li)
end do
beta(jj) = sqrt(norm(Rs))
if(beta(1) < Cp%lanczos_tol) then
! Psi is an eigenstate of H to within the given tolerance
call scale(exp(deltat * alpha(1)), Psi%Aa(max(Psi%oc, 1)))
call destroy(Rs)
call destroy(Krylovuv(1))
deallocate(Krylovuv, alpha, beta, propsave)
return
end if
call destroy(Psi)
! Further iterations to build Krylov basis
! ----------------------------------------
do jj = 2, nkrylov
call copy(Krylovuv(jj), Rs, scalar=done / beta(jj - 1))
call destroy(Rs)
call fit_hpsi(Rs, Ham, Krylovuv(jj), Cp, err, converged, &
init_random, 'V', 'N', errst=errst)
!if(prop_error('krylov_method_mpo: fit_hpsi (2) '//&
! 'failed', errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method: faild on jj =', &
jj, 'Hpsi with error', err
alpha(jj) = real(dot(Rs, Krylovuv(jj)), KIND=rKind)
do kk = (jj - 1), jj
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc)
end do
! Three-term recurrence
call orth(Rs, Krylovuv, (jj - 1), jj, LR, Cp, err, converged, &
errst=errst)
!if(prop_error('krylov_method_mpo : orth (3) failed.', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method:, failed '//&
'on jj =', jj, 'three-term recurrence with error', err
do kk = 1, (jj - 2)
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc)
end do
! Reorthogonalize
call orth(Rs, Krylovuv, 1, jj, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_method_mpo : orth (4) failed.', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method:, failed '//&
'on jj =', jj, 'reorthogonalization', err
do kk = 1, jj
do pp = 1, size(LR(kk)%Li, 1)
call destroy(LR(kk)%Li(pp))
end do
deallocate(LR(kk)%Li)
end do
beta(jj) = sqrt(norm(Rs))
allocate(prop(jj, jj))
call expmh(prop, deltat, alpha(1:jj), beta(1:jj - 1), jj, &
errst=errst)
!if(prop_error('krylov_method_MPOM_TYPE_real: '//&
! 'expmh failed.', 'KrylovOps_include.f90:274', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * beta(jj - 1) * abs(prop(jj, 1)))
if(abs(estimate) <= Cp%lanczos_tol) exit
propsave(1:jj - 1) = prop(1:jj - 1, 1) - propsave(1:jj - 1)
estimate = real(dot_product(propsave(1:jj - 1) , &
propsave(1:jj - 1)), KIND=rKind)
if(estimate <= Cp%psi_tol) exit
propsave(1:jj) = prop(:, 1)
! Avoid increment to jj = nkrylov + 1
if(jj == nkrylov) exit
deallocate(prop)
end do
! To-Do : seems to be repetitive
call expmh(prop, deltat, alpha(1:jj), beta(1:jj - 1), jj, errst=errst)
!if(prop_error('krylov_method_MPOM_TYPE_real: '//&
! 'expmh failed.', 'KrylovOps_include.f90:298', &
! errst=errst)) then
! return
!end if
call destroy(Rs)
! location of largest weight
kk = maxloc(abs(prop(1:jj, 1)), 1)
! Initialize variational ansatz to be the largest weight MPS
call copy(Psi, Krylovuv(kk))
call fitmps(Psi, prop(1:jj, 1), Krylovuv, Cp, err, converged)
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method: failed on '//&
'fitting with error', err
do ii = 1, jj
call destroy(Krylovuv(ii))
end do
deallocate(Krylovuv, LR, alpha, beta, prop, propsave)
end subroutine krylov_method_mpo_real
"""
return
[docs]def krylov_method_mpoc_real():
"""
fortran-subroutine - August 2017 (dj, updated)
Propagate the wave function Psi with the Krylov/Lanczos algorithm.
**Arguments**
Psi : TYPE(mpsc), inout
This is the wave function to be propagated.
deltat : real, in
The time step for the evolution.
Ham : TYPE(mpoc), inout
Evolve the wave function under this Hamiltonian.
converged : LOGICAL, out
On exit, it indicates if the algorithm is converged within then
given tolerances.
cerr : REAL, inout
The cumulated error. The error done during this subroutine is
added to the incoming value.
Cp : TYPE(ConvParam), inout
Stores all the convergence parameters necessary for the algorithm.
init_random : CHARACTER, inout
Determines the initialization in the fit of H | psi >.
renorm : LOGICAL, in
Flag if state vector should be renormalized to 1 (true).
'N' : do not normalize (default); 'M' : normalize for MPS
1 / sqrt(norm);
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
Currently uses complete reorthogonalization for stability.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_method_mpoc_real(Psi, deltat, Ham, converged, &
cerr, Cp, init_random, renorm, pbc, errst)
type(mpsc), intent(inout) :: Psi
real(KIND=rKind), intent(in) :: deltat
type(mpoc), intent(inout) :: Ham
logical, intent(out) :: converged
real(KIND=rKind), intent(inout) :: cerr
type(ConvParam), intent(in) :: Cp
character, intent(in) :: init_random
character, intent(in) :: renorm
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: jj, kk, pp
! for looping / unused index
integer :: ii
! diagonal and 1st offdiagonal of tridiagonal upper hessenberg
real(KIND=rKind), allocatable :: beta(:), alpha(:)
! number of Krylov vectors and Krylov basis
integer :: nkrylov
type(mpsc), pointer :: Krylovuv(:)
! MPS to be modified in iteration process
type(mpsc) :: Rs
! the overlap between two states
real(KIND=rKind) :: overlap
! error of a single call; added to cumulative error
real(KIND=rKind) :: err
! left right overlap
type(tensorlistc), pointer :: LR(:)
! The exponential in the Krylov basis, including the previous step
real(KIND=rKind), allocatable :: prop(:, :)
real(KIND=rKind), allocatable :: propsave(:)
! estimate of convergence based on propagator in Krylov space
real(KIND=rKind) :: estimate
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_method_mpoc_real:'//&
! 'Krylov does not have PBC.', 99, 'KrylovOps_include.f90:134', &
! errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylovuv(nkrylov), LR(nkrylov), alpha(nkrylov), &
beta(0:nkrylov), propsave(nkrylov))
propsave(1) = 10.0_rkind
call copy(Rs, Psi)
beta(0) = sqrt(norm(Rs))
! First iteration
! ---------------
jj = 1
call copy(Krylovuv(jj), Rs, scalar=zone / beta(jj - 1))
call destroy(Rs)
call fit_hpsi(Rs, Ham, Krylovuv(jj), Cp, err, converged, init_random, &
'V', 'N', errst=errst)
!if(prop_error('krylov_method_mpoc: fit_hpsi (1) failed', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method: failed on jj =', &
jj, 'Hpsi with error', err
alpha(jj) = real(dot(Rs, Krylovuv(jj)), KIND=rKind)
! Setup LR overlaps from j to r - use them to get alpha
do kk = 1, jj
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc)
end do
! Orthogonalize
call orth(Rs, Krylovuv, 1, 1, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_method_mpoc : orth (1) failed.', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method:, failed on jj =', &
jj, 'three-term recurrence with error', err
! Reorthogonalize
call orth(Rs, Krylovuv, 1, jj, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_method_mpoc : orth (2) failed.', &
! errst=errst)) return
if(.not. converged) write(slog, *) 'krylov_method:, failed on jj =', &
jj, 'reorthogonalization', err
do kk = 1, jj
do pp = 1, size(LR(kk)%Li, 1)
call destroy(LR(kk)%Li(pp))
end do
deallocate(LR(kk)%Li)
end do
beta(jj) = sqrt(norm(Rs))
if(beta(1) < Cp%lanczos_tol) then
! Psi is an eigenstate of H to within the given tolerance
call scale(exp(deltat * alpha(1)), Psi%Aa(max(Psi%oc, 1)))
call destroy(Rs)
call destroy(Krylovuv(1))
deallocate(Krylovuv, alpha, beta, propsave)
return
end if
call destroy(Psi)
! Further iterations to build Krylov basis
! ----------------------------------------
do jj = 2, nkrylov
call copy(Krylovuv(jj), Rs, scalar=zone / beta(jj - 1))
call destroy(Rs)
call fit_hpsi(Rs, Ham, Krylovuv(jj), Cp, err, converged, &
init_random, 'V', 'N', errst=errst)
!if(prop_error('krylov_method_mpoc: fit_hpsi (2) '//&
! 'failed', errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method: faild on jj =', &
jj, 'Hpsi with error', err
alpha(jj) = real(dot(Rs, Krylovuv(jj)), KIND=rKind)
do kk = (jj - 1), jj
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc)
end do
! Three-term recurrence
call orth(Rs, Krylovuv, (jj - 1), jj, LR, Cp, err, converged, &
errst=errst)
!if(prop_error('krylov_method_mpoc : orth (3) failed.', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method:, failed '//&
'on jj =', jj, 'three-term recurrence with error', err
do kk = 1, (jj - 2)
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc)
end do
! Reorthogonalize
call orth(Rs, Krylovuv, 1, jj, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_method_mpoc : orth (4) failed.', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method:, failed '//&
'on jj =', jj, 'reorthogonalization', err
do kk = 1, jj
do pp = 1, size(LR(kk)%Li, 1)
call destroy(LR(kk)%Li(pp))
end do
deallocate(LR(kk)%Li)
end do
beta(jj) = sqrt(norm(Rs))
allocate(prop(jj, jj))
call expmh(prop, deltat, alpha(1:jj), beta(1:jj - 1), jj, &
errst=errst)
!if(prop_error('krylov_method_MPOM_TYPE_real: '//&
! 'expmh failed.', 'KrylovOps_include.f90:274', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * beta(jj - 1) * abs(prop(jj, 1)))
if(abs(estimate) <= Cp%lanczos_tol) exit
propsave(1:jj - 1) = prop(1:jj - 1, 1) - propsave(1:jj - 1)
estimate = real(dot_product(propsave(1:jj - 1) , &
propsave(1:jj - 1)), KIND=rKind)
if(estimate <= Cp%psi_tol) exit
propsave(1:jj) = prop(:, 1)
! Avoid increment to jj = nkrylov + 1
if(jj == nkrylov) exit
deallocate(prop)
end do
! To-Do : seems to be repetitive
call expmh(prop, deltat, alpha(1:jj), beta(1:jj - 1), jj, errst=errst)
!if(prop_error('krylov_method_MPOM_TYPE_real: '//&
! 'expmh failed.', 'KrylovOps_include.f90:298', &
! errst=errst)) then
! return
!end if
call destroy(Rs)
! location of largest weight
kk = maxloc(abs(prop(1:jj, 1)), 1)
! Initialize variational ansatz to be the largest weight MPS
call copy(Psi, Krylovuv(kk))
call fitmps(Psi, prop(1:jj, 1), Krylovuv, Cp, err, converged)
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method: failed on '//&
'fitting with error', err
do ii = 1, jj
call destroy(Krylovuv(ii))
end do
deallocate(Krylovuv, LR, alpha, beta, prop, propsave)
end subroutine krylov_method_mpoc_real
"""
return
[docs]def krylov_method_qmpo_real():
"""
fortran-subroutine - August 2017 (dj, updated)
Propagate the wave function Psi with the Krylov/Lanczos algorithm.
**Arguments**
Psi : TYPE(qmps), inout
This is the wave function to be propagated.
deltat : real, in
The time step for the evolution.
Ham : TYPE(qmpo), inout
Evolve the wave function under this Hamiltonian.
converged : LOGICAL, out
On exit, it indicates if the algorithm is converged within then
given tolerances.
cerr : REAL, inout
The cumulated error. The error done during this subroutine is
added to the incoming value.
Cp : TYPE(ConvParam), inout
Stores all the convergence parameters necessary for the algorithm.
init_random : CHARACTER, inout
Determines the initialization in the fit of H | psi >.
renorm : LOGICAL, in
Flag if state vector should be renormalized to 1 (true).
'N' : do not normalize (default); 'M' : normalize for MPS
1 / sqrt(norm);
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
Currently uses complete reorthogonalization for stability.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_method_qmpo_real(Psi, deltat, Ham, converged, &
cerr, Cp, init_random, renorm, pbc, errst)
type(qmps), intent(inout) :: Psi
real(KIND=rKind), intent(in) :: deltat
type(qmpo), intent(inout) :: Ham
logical, intent(out) :: converged
real(KIND=rKind), intent(inout) :: cerr
type(ConvParam), intent(in) :: Cp
character, intent(in) :: init_random
character, intent(in) :: renorm
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: jj, kk, pp
! for looping / unused index
integer :: ii
! diagonal and 1st offdiagonal of tridiagonal upper hessenberg
real(KIND=rKind), allocatable :: beta(:), alpha(:)
! number of Krylov vectors and Krylov basis
integer :: nkrylov
type(qmps), pointer :: Krylovuv(:)
! MPS to be modified in iteration process
type(qmps) :: Rs
! the overlap between two states
real(KIND=rKind) :: overlap
! error of a single call; added to cumulative error
real(KIND=rKind) :: err
! left right overlap
type(qtensorlist), pointer :: LR(:)
! The exponential in the Krylov basis, including the previous step
real(KIND=rKind), allocatable :: prop(:, :)
real(KIND=rKind), allocatable :: propsave(:)
! estimate of convergence based on propagator in Krylov space
real(KIND=rKind) :: estimate
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_method_qmpo_real:'//&
! 'Krylov does not have PBC.', 99, 'KrylovOps_include.f90:134', &
! errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylovuv(nkrylov), LR(nkrylov), alpha(nkrylov), &
beta(0:nkrylov), propsave(nkrylov))
propsave(1) = 10.0_rkind
call copy(Rs, Psi)
beta(0) = sqrt(norm(Rs))
! First iteration
! ---------------
jj = 1
call copy(Krylovuv(jj), Rs, scalar=done / beta(jj - 1))
call destroy(Rs)
call fit_hpsi(Rs, Ham, Krylovuv(jj), Cp, err, converged, init_random, &
'V', 'N', errst=errst)
!if(prop_error('krylov_method_qmpo: fit_hpsi (1) failed', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method: failed on jj =', &
jj, 'Hpsi with error', err
alpha(jj) = real(dot(Rs, Krylovuv(jj)), KIND=rKind)
! Setup LR overlaps from j to r - use them to get alpha
do kk = 1, jj
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc)
end do
! Orthogonalize
call orth(Rs, Krylovuv, 1, 1, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_method_qmpo : orth (1) failed.', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method:, failed on jj =', &
jj, 'three-term recurrence with error', err
! Reorthogonalize
call orth(Rs, Krylovuv, 1, jj, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_method_qmpo : orth (2) failed.', &
! errst=errst)) return
if(.not. converged) write(slog, *) 'krylov_method:, failed on jj =', &
jj, 'reorthogonalization', err
do kk = 1, jj
do pp = 1, size(LR(kk)%Li, 1)
call destroy(LR(kk)%Li(pp))
end do
deallocate(LR(kk)%Li)
end do
beta(jj) = sqrt(norm(Rs))
if(beta(1) < Cp%lanczos_tol) then
! Psi is an eigenstate of H to within the given tolerance
call scale(exp(deltat * alpha(1)), Psi%Aa(max(Psi%oc, 1)))
call destroy(Rs)
call destroy(Krylovuv(1))
deallocate(Krylovuv, alpha, beta, propsave)
return
end if
call destroy(Psi)
! Further iterations to build Krylov basis
! ----------------------------------------
do jj = 2, nkrylov
call copy(Krylovuv(jj), Rs, scalar=done / beta(jj - 1))
call destroy(Rs)
call fit_hpsi(Rs, Ham, Krylovuv(jj), Cp, err, converged, &
init_random, 'V', 'N', errst=errst)
!if(prop_error('krylov_method_qmpo: fit_hpsi (2) '//&
! 'failed', errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method: faild on jj =', &
jj, 'Hpsi with error', err
alpha(jj) = real(dot(Rs, Krylovuv(jj)), KIND=rKind)
do kk = (jj - 1), jj
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc)
end do
! Three-term recurrence
call orth(Rs, Krylovuv, (jj - 1), jj, LR, Cp, err, converged, &
errst=errst)
!if(prop_error('krylov_method_qmpo : orth (3) failed.', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method:, failed '//&
'on jj =', jj, 'three-term recurrence with error', err
do kk = 1, (jj - 2)
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc)
end do
! Reorthogonalize
call orth(Rs, Krylovuv, 1, jj, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_method_qmpo : orth (4) failed.', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method:, failed '//&
'on jj =', jj, 'reorthogonalization', err
do kk = 1, jj
do pp = 1, size(LR(kk)%Li, 1)
call destroy(LR(kk)%Li(pp))
end do
deallocate(LR(kk)%Li)
end do
beta(jj) = sqrt(norm(Rs))
allocate(prop(jj, jj))
call expmh(prop, deltat, alpha(1:jj), beta(1:jj - 1), jj, &
errst=errst)
!if(prop_error('krylov_method_MPOM_TYPE_real: '//&
! 'expmh failed.', 'KrylovOps_include.f90:274', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * beta(jj - 1) * abs(prop(jj, 1)))
if(abs(estimate) <= Cp%lanczos_tol) exit
propsave(1:jj - 1) = prop(1:jj - 1, 1) - propsave(1:jj - 1)
estimate = real(dot_product(propsave(1:jj - 1) , &
propsave(1:jj - 1)), KIND=rKind)
if(estimate <= Cp%psi_tol) exit
propsave(1:jj) = prop(:, 1)
! Avoid increment to jj = nkrylov + 1
if(jj == nkrylov) exit
deallocate(prop)
end do
! To-Do : seems to be repetitive
call expmh(prop, deltat, alpha(1:jj), beta(1:jj - 1), jj, errst=errst)
!if(prop_error('krylov_method_MPOM_TYPE_real: '//&
! 'expmh failed.', 'KrylovOps_include.f90:298', &
! errst=errst)) then
! return
!end if
call destroy(Rs)
! location of largest weight
kk = maxloc(abs(prop(1:jj, 1)), 1)
! Initialize variational ansatz to be the largest weight MPS
call copy(Psi, Krylovuv(kk))
call fitmps(Psi, prop(1:jj, 1), Krylovuv, Cp, err, converged)
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method: failed on '//&
'fitting with error', err
do ii = 1, jj
call destroy(Krylovuv(ii))
end do
deallocate(Krylovuv, LR, alpha, beta, prop, propsave)
end subroutine krylov_method_qmpo_real
"""
return
[docs]def krylov_method_qmpoc_real():
"""
fortran-subroutine - August 2017 (dj, updated)
Propagate the wave function Psi with the Krylov/Lanczos algorithm.
**Arguments**
Psi : TYPE(qmpsc), inout
This is the wave function to be propagated.
deltat : real, in
The time step for the evolution.
Ham : TYPE(qmpoc), inout
Evolve the wave function under this Hamiltonian.
converged : LOGICAL, out
On exit, it indicates if the algorithm is converged within then
given tolerances.
cerr : REAL, inout
The cumulated error. The error done during this subroutine is
added to the incoming value.
Cp : TYPE(ConvParam), inout
Stores all the convergence parameters necessary for the algorithm.
init_random : CHARACTER, inout
Determines the initialization in the fit of H | psi >.
renorm : LOGICAL, in
Flag if state vector should be renormalized to 1 (true).
'N' : do not normalize (default); 'M' : normalize for MPS
1 / sqrt(norm);
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
Currently uses complete reorthogonalization for stability.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_method_qmpoc_real(Psi, deltat, Ham, converged, &
cerr, Cp, init_random, renorm, pbc, errst)
type(qmpsc), intent(inout) :: Psi
real(KIND=rKind), intent(in) :: deltat
type(qmpoc), intent(inout) :: Ham
logical, intent(out) :: converged
real(KIND=rKind), intent(inout) :: cerr
type(ConvParam), intent(in) :: Cp
character, intent(in) :: init_random
character, intent(in) :: renorm
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: jj, kk, pp
! for looping / unused index
integer :: ii
! diagonal and 1st offdiagonal of tridiagonal upper hessenberg
real(KIND=rKind), allocatable :: beta(:), alpha(:)
! number of Krylov vectors and Krylov basis
integer :: nkrylov
type(qmpsc), pointer :: Krylovuv(:)
! MPS to be modified in iteration process
type(qmpsc) :: Rs
! the overlap between two states
real(KIND=rKind) :: overlap
! error of a single call; added to cumulative error
real(KIND=rKind) :: err
! left right overlap
type(qtensorclist), pointer :: LR(:)
! The exponential in the Krylov basis, including the previous step
real(KIND=rKind), allocatable :: prop(:, :)
real(KIND=rKind), allocatable :: propsave(:)
! estimate of convergence based on propagator in Krylov space
real(KIND=rKind) :: estimate
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_method_qmpoc_real:'//&
! 'Krylov does not have PBC.', 99, 'KrylovOps_include.f90:134', &
! errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylovuv(nkrylov), LR(nkrylov), alpha(nkrylov), &
beta(0:nkrylov), propsave(nkrylov))
propsave(1) = 10.0_rkind
call copy(Rs, Psi)
beta(0) = sqrt(norm(Rs))
! First iteration
! ---------------
jj = 1
call copy(Krylovuv(jj), Rs, scalar=zone / beta(jj - 1))
call destroy(Rs)
call fit_hpsi(Rs, Ham, Krylovuv(jj), Cp, err, converged, init_random, &
'V', 'N', errst=errst)
!if(prop_error('krylov_method_qmpoc: fit_hpsi (1) failed', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method: failed on jj =', &
jj, 'Hpsi with error', err
alpha(jj) = real(dot(Rs, Krylovuv(jj)), KIND=rKind)
! Setup LR overlaps from j to r - use them to get alpha
do kk = 1, jj
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc)
end do
! Orthogonalize
call orth(Rs, Krylovuv, 1, 1, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_method_qmpoc : orth (1) failed.', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method:, failed on jj =', &
jj, 'three-term recurrence with error', err
! Reorthogonalize
call orth(Rs, Krylovuv, 1, jj, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_method_qmpoc : orth (2) failed.', &
! errst=errst)) return
if(.not. converged) write(slog, *) 'krylov_method:, failed on jj =', &
jj, 'reorthogonalization', err
do kk = 1, jj
do pp = 1, size(LR(kk)%Li, 1)
call destroy(LR(kk)%Li(pp))
end do
deallocate(LR(kk)%Li)
end do
beta(jj) = sqrt(norm(Rs))
if(beta(1) < Cp%lanczos_tol) then
! Psi is an eigenstate of H to within the given tolerance
call scale(exp(deltat * alpha(1)), Psi%Aa(max(Psi%oc, 1)))
call destroy(Rs)
call destroy(Krylovuv(1))
deallocate(Krylovuv, alpha, beta, propsave)
return
end if
call destroy(Psi)
! Further iterations to build Krylov basis
! ----------------------------------------
do jj = 2, nkrylov
call copy(Krylovuv(jj), Rs, scalar=zone / beta(jj - 1))
call destroy(Rs)
call fit_hpsi(Rs, Ham, Krylovuv(jj), Cp, err, converged, &
init_random, 'V', 'N', errst=errst)
!if(prop_error('krylov_method_qmpoc: fit_hpsi (2) '//&
! 'failed', errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method: faild on jj =', &
jj, 'Hpsi with error', err
alpha(jj) = real(dot(Rs, Krylovuv(jj)), KIND=rKind)
do kk = (jj - 1), jj
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc)
end do
! Three-term recurrence
call orth(Rs, Krylovuv, (jj - 1), jj, LR, Cp, err, converged, &
errst=errst)
!if(prop_error('krylov_method_qmpoc : orth (3) failed.', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method:, failed '//&
'on jj =', jj, 'three-term recurrence with error', err
do kk = 1, (jj - 2)
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc)
end do
! Reorthogonalize
call orth(Rs, Krylovuv, 1, jj, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_method_qmpoc : orth (4) failed.', &
! errst=errst)) return
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method:, failed '//&
'on jj =', jj, 'reorthogonalization', err
do kk = 1, jj
do pp = 1, size(LR(kk)%Li, 1)
call destroy(LR(kk)%Li(pp))
end do
deallocate(LR(kk)%Li)
end do
beta(jj) = sqrt(norm(Rs))
allocate(prop(jj, jj))
call expmh(prop, deltat, alpha(1:jj), beta(1:jj - 1), jj, &
errst=errst)
!if(prop_error('krylov_method_MPOM_TYPE_real: '//&
! 'expmh failed.', 'KrylovOps_include.f90:274', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * beta(jj - 1) * abs(prop(jj, 1)))
if(abs(estimate) <= Cp%lanczos_tol) exit
propsave(1:jj - 1) = prop(1:jj - 1, 1) - propsave(1:jj - 1)
estimate = real(dot_product(propsave(1:jj - 1) , &
propsave(1:jj - 1)), KIND=rKind)
if(estimate <= Cp%psi_tol) exit
propsave(1:jj) = prop(:, 1)
! Avoid increment to jj = nkrylov + 1
if(jj == nkrylov) exit
deallocate(prop)
end do
! To-Do : seems to be repetitive
call expmh(prop, deltat, alpha(1:jj), beta(1:jj - 1), jj, errst=errst)
!if(prop_error('krylov_method_MPOM_TYPE_real: '//&
! 'expmh failed.', 'KrylovOps_include.f90:298', &
! errst=errst)) then
! return
!end if
call destroy(Rs)
! location of largest weight
kk = maxloc(abs(prop(1:jj, 1)), 1)
! Initialize variational ansatz to be the largest weight MPS
call copy(Psi, Krylovuv(kk))
call fitmps(Psi, prop(1:jj, 1), Krylovuv, Cp, err, converged)
cerr = cerr + err
if(.not. converged) write(slog, *) 'krylov_method: failed on '//&
'fitting with error', err
do ii = 1, jj
call destroy(Krylovuv(ii))
end do
deallocate(Krylovuv, LR, alpha, beta, prop, propsave)
end subroutine krylov_method_qmpoc_real
"""
return
[docs]def krylov_arnoldi_mpoc():
"""
fortran-subroutine - August 2017 (dj)
Propagate Psi, possibly representing an MPDO, with the Arnoldi
algorithm.
**Arguments**
Psi : TYPE(mpsc), inout
Propagate this MPS under the Hamiltonian passed as MPO and
its parameters.
deltat : COMPLEX, in
Time step for the evolution step (including -i)
H : TYPE(mpoc), inout
MPO representing the Hamiltonian.
maxbondd : INTEGER, inout
Maximal bond dimension :math:`\\chi`.
converged : LOGICAL, inout
Convergence flag of last subroutine called (FitMPS).
cerr : REAL, inout
The cumulated error. The error done during this subroutine is
added to the incoming value.
random_init : CHARACTER, in
Flag for initial state of Fit_Hpsi.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
Currently uses complete reorthogonalization for stability.
(defined in TimeEvolutionOps_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_arnoldi_mpoc(Psi, deltat, Ham, converged, &
cerr, Cp, init_random, pbc, errst)
type(mpsc), intent(inout) :: Psi
complex(KIND=rKind), intent(in) :: deltat
type(mpoc), intent(inout) :: Ham
logical, intent(out) :: converged
real(KIND=rKind), intent(inout) :: cerr
type(ConvParam), intent(in) :: Cp
character, intent(in) :: init_random
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: jj, kk, pp
! number of krylov vectors (copy of maxnlanczositerations)
integer :: nkrylov
! Upper Hessenberg matrix / upper Hessenberg weighted with dt
complex(KIND=rKind), dimension(:, :), allocatable :: hessen, hessendt
! intial norm
real(KIND=rKind) :: beta0
! vector for Krylov algorithm - modification carried out on r
type(mpsc) :: Rs
! the krylov vectors as MPS representation
type(mpsc), pointer :: krylovuv(:)
! Overlaps
type(tensorlistc), pointer :: LR(:)
! Estimate for error to check tolerance and error of called functions
real(KIND=rKind) :: estimate, err
! propagator and first column of last propagator
complex(KIND=rKind), ALLOCATABLE :: prop(:, :), propsave(:)
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_arnoldi_mpoc:'//&
! 'Krylov-Arnoldi does not have PBC.', 99, 'KrylovOps_include.f90:452', &
! errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylovuv(nkrylov), LR(nkrylov), &
hessen(nkrylov + 1, nkrylov + 1), propsave(nkrylov))
propsave(1) = 10.0_rKind
hessen = 0.0_rKind
call copy(Rs, Psi)
beta0 = sqrt(norm(Rs, errst=errst))
!if(prop_error('krylov_arnoldi_mpoc : norm (1) failed.', &
! errst=errst)) return
! First data point - initialize values
! ------------------------------------
jj = 1
call copy(Krylovuv(jj), Rs, scalar=zone / beta0)
call destroy(Rs)
call Fit_Hpsi(Rs, Ham, Krylovuv(jj), Cp, err, converged, &
init_random, 'V', 'N', errst=errst)
!if(prop_error('krylov_arnoldi_mpoc : fit_hpsi (1) failed.', &
! errst=errst)) return
cerr = cerr + err
hessen(jj + 1, jj) = sqrt(norm(Rs, errst=errst))
!if(prop_error('krylov_arnoldi_mpoc : norm (2) failed.', &
! errst=errst)) return
if(real(hessen(jj + 1, jj), KIND=rKind) < Cp%lanczos_tol) then
! Psi is an eigenstate of H to within the given tolerance
call scale(exp(deltat * hessen(1 , 1)), Psi%Aa(max(Psi%oc, 1)))
call destroy(Rs)
call destroy(Krylovuv(1))
deallocate(Krylovuv, hessen, propsave, LR)
return
end if
call destroy(Psi)
! Iteration over vectors 2, ...
! -----------------------------
do jj = 2, nkrylov
call copy(Krylovuv(jj), Rs)
call scale(1.0_rKind / real(hessen(jj, jj - 1), KIND=rKind), &
Krylovuv(jj))
call destroy(Rs)
call fit_hpsi(Rs, Ham, Krylovuv(jj), Cp, err, converged, &
init_random, 'V', 'N', errst=errst)
!if(prop_error('krylov_arnoldi_mpoc : fit_hpsi (2) '//&
! 'failed.', errst=errst)) return
cerr = cerr + err
if((.not. converged) .and. (verbose > 0)) write(slog, *) &
'krylov_arnoldi : failed on jj = ', jj, &
' Hpsi with error ', err
do kk = 1, (jj - 1)
hessen(kk, jj) = dot(Rs, Krylovuv(kk))
! exactdiag orthogonalizes directly after this step
end do
do kk = (jj - 1), jj
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc, errst=errst)
!if(prop_error('krylov_arnoldi_mpoc : setuplr (2) '//&
! 'failed.', errst=errst)) return
end do
call orth(Rs, Krylovuv, jj - 1, jj, LR, Cp, err, converged, &
errst=errst)
!if(prop_error('krylov_arnoldi_mpoc : orth (3) '//&
! 'failed.', errst=errst)) return
cerr = cerr + err
if((.not. converged) .and. (verbose > 0)) write(slog, *) &
'krylov_arnoldi : failed on jj = ', jj, &
' three-term recurrence with error', err
do kk = 1, jj - 2
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc, errst=errst)
!if(prop_error('krylov_arnoldi_mpoc : setuplr (3) '//&
! 'failed.', errst=errst)) return
end do
! Reorthogonalize
call orth(Rs, Krylovuv, 1, jj, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_arnoldi_mpoc : orth (4) '//&
! 'failed.', errst=errst)) return
cerr = cerr + err
if((.not. converged) .and. (verbose > 0)) write(slog, *) &
'krylov_arnoldi : failed on jj = ', jj, &
' reorthogonalization with error', err
do kk = 1, jj
do pp = 1, size(LR(kk)%Li)
call destroy(LR(kk)%Li(pp))
end do
deallocate(LR(kk)%Li)
end do
hessen(jj + 1, jj) = sqrt(norm(Rs, errst=errst))
!if(prop_error('krylov_arnoldi_mpoc : norm (3) failed.', &
! errst=errst)) return
! Exponentiate matrix
!call Exponentiate(prop(1:jj, 1:jj), deltat, hessen(1:jj, 1:jj))
allocate(prop(jj, jj), hessendt(jj, jj))
hessendt = hessen(1:jj, 1:jj)
call expm(prop, deltat, hessendt, jj, errst=errst)
!if(prop_error('krylov_arnoldi_MPOM_TYPE: '//&
! 'expm failed.', 'KrylovOps_include.f90:574', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * hessen(jj, jj - 1) * abs(prop(jj, 1)))
if(abs(estimate) < Cp%lanczos_tol) exit
propsave(1:jj - 1) = prop(1:jj - 1, 1) - propsave(1:jj - 1)
estimate = real(dot_product(propsave(1:jj - 1), &
propsave(1:jj - 1)), KIND=rKind)
!if(estimate < Cp%psi_tol) exit
if(abs(hessen(jj + 1, jj)) < Cp%Psi_tol) exit
propsave(1:jj) = prop(1:jj, 1)
! Avoid increment to jj = nkrylov + 1
if(jj == nkrylov) exit
deallocate(prop, hessendt)
end do
! Get the final propagated state
! ------------------------------
call destroy(Rs)
! location of largest and 2nd largest
kk = maxloc(abs(prop(1:jj, 1)), 1)
! Initialize variational ansatz to be the largest weight MPS
call copy(Psi, Krylovuv(kk))
prop(1:jj, 1) = prop(1:jj, 1) * beta0 !!! NEW
call fitmps(Psi, prop(1:jj, 1), Krylovuv, Cp, err, converged, &
errst=errst)
!if(prop_error('krylov_arnoldi_mpoc : fitmps failed.', &
! errst=errst)) return
cerr = cerr + err
if((.not. converged) .and. (verbose > 0)) write(slog, *) &
'krylov_arnoldi : failed on fitting with error ', err
do kk = 1, jj
call destroy(Krylovuv(kk))
end do
deallocate(Krylovuv, LR, hessen, prop, propsave, hessendt)
end subroutine krylov_arnoldi_mpoc
"""
return
[docs]def krylov_arnoldi_qmpoc():
"""
fortran-subroutine - August 2017 (dj)
Propagate Psi, possibly representing an MPDO, with the Arnoldi
algorithm.
**Arguments**
Psi : TYPE(qmpsc), inout
Propagate this MPS under the Hamiltonian passed as MPO and
its parameters.
deltat : COMPLEX, in
Time step for the evolution step (including -i)
H : TYPE(qmpoc), inout
MPO representing the Hamiltonian.
maxbondd : INTEGER, inout
Maximal bond dimension :math:`\\chi`.
converged : LOGICAL, inout
Convergence flag of last subroutine called (FitMPS).
cerr : REAL, inout
The cumulated error. The error done during this subroutine is
added to the incoming value.
random_init : CHARACTER, in
Flag for initial state of Fit_Hpsi.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
Currently uses complete reorthogonalization for stability.
(defined in TimeEvolutionOps_include.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_arnoldi_qmpoc(Psi, deltat, Ham, converged, &
cerr, Cp, init_random, pbc, errst)
type(qmpsc), intent(inout) :: Psi
complex(KIND=rKind), intent(in) :: deltat
type(qmpoc), intent(inout) :: Ham
logical, intent(out) :: converged
real(KIND=rKind), intent(inout) :: cerr
type(ConvParam), intent(in) :: Cp
character, intent(in) :: init_random
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: jj, kk, pp
! number of krylov vectors (copy of maxnlanczositerations)
integer :: nkrylov
! Upper Hessenberg matrix / upper Hessenberg weighted with dt
complex(KIND=rKind), dimension(:, :), allocatable :: hessen, hessendt
! intial norm
real(KIND=rKind) :: beta0
! vector for Krylov algorithm - modification carried out on r
type(qmpsc) :: Rs
! the krylov vectors as MPS representation
type(qmpsc), pointer :: krylovuv(:)
! Overlaps
type(qtensorclist), pointer :: LR(:)
! Estimate for error to check tolerance and error of called functions
real(KIND=rKind) :: estimate, err
! propagator and first column of last propagator
complex(KIND=rKind), ALLOCATABLE :: prop(:, :), propsave(:)
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_arnoldi_qmpoc:'//&
! 'Krylov-Arnoldi does not have PBC.', 99, 'KrylovOps_include.f90:452', &
! errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylovuv(nkrylov), LR(nkrylov), &
hessen(nkrylov + 1, nkrylov + 1), propsave(nkrylov))
propsave(1) = 10.0_rKind
hessen = 0.0_rKind
call copy(Rs, Psi)
beta0 = sqrt(norm(Rs, errst=errst))
!if(prop_error('krylov_arnoldi_qmpoc : norm (1) failed.', &
! errst=errst)) return
! First data point - initialize values
! ------------------------------------
jj = 1
call copy(Krylovuv(jj), Rs, scalar=zone / beta0)
call destroy(Rs)
call Fit_Hpsi(Rs, Ham, Krylovuv(jj), Cp, err, converged, &
init_random, 'V', 'N', errst=errst)
!if(prop_error('krylov_arnoldi_qmpoc : fit_hpsi (1) failed.', &
! errst=errst)) return
cerr = cerr + err
hessen(jj + 1, jj) = sqrt(norm(Rs, errst=errst))
!if(prop_error('krylov_arnoldi_qmpoc : norm (2) failed.', &
! errst=errst)) return
if(real(hessen(jj + 1, jj), KIND=rKind) < Cp%lanczos_tol) then
! Psi is an eigenstate of H to within the given tolerance
call scale(exp(deltat * hessen(1 , 1)), Psi%Aa(max(Psi%oc, 1)))
call destroy(Rs)
call destroy(Krylovuv(1))
deallocate(Krylovuv, hessen, propsave, LR)
return
end if
call destroy(Psi)
! Iteration over vectors 2, ...
! -----------------------------
do jj = 2, nkrylov
call copy(Krylovuv(jj), Rs)
call scale(1.0_rKind / real(hessen(jj, jj - 1), KIND=rKind), &
Krylovuv(jj))
call destroy(Rs)
call fit_hpsi(Rs, Ham, Krylovuv(jj), Cp, err, converged, &
init_random, 'V', 'N', errst=errst)
!if(prop_error('krylov_arnoldi_qmpoc : fit_hpsi (2) '//&
! 'failed.', errst=errst)) return
cerr = cerr + err
if((.not. converged) .and. (verbose > 0)) write(slog, *) &
'krylov_arnoldi : failed on jj = ', jj, &
' Hpsi with error ', err
do kk = 1, (jj - 1)
hessen(kk, jj) = dot(Rs, Krylovuv(kk))
! exactdiag orthogonalizes directly after this step
end do
do kk = (jj - 1), jj
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc, errst=errst)
!if(prop_error('krylov_arnoldi_qmpoc : setuplr (2) '//&
! 'failed.', errst=errst)) return
end do
call orth(Rs, Krylovuv, jj - 1, jj, LR, Cp, err, converged, &
errst=errst)
!if(prop_error('krylov_arnoldi_qmpoc : orth (3) '//&
! 'failed.', errst=errst)) return
cerr = cerr + err
if((.not. converged) .and. (verbose > 0)) write(slog, *) &
'krylov_arnoldi : failed on jj = ', jj, &
' three-term recurrence with error', err
do kk = 1, jj - 2
call setuplr(LR(kk), Rs, Krylovuv(kk), Rs%oc, errst=errst)
!if(prop_error('krylov_arnoldi_qmpoc : setuplr (3) '//&
! 'failed.', errst=errst)) return
end do
! Reorthogonalize
call orth(Rs, Krylovuv, 1, jj, LR, Cp, err, converged, errst=errst)
!if(prop_error('krylov_arnoldi_qmpoc : orth (4) '//&
! 'failed.', errst=errst)) return
cerr = cerr + err
if((.not. converged) .and. (verbose > 0)) write(slog, *) &
'krylov_arnoldi : failed on jj = ', jj, &
' reorthogonalization with error', err
do kk = 1, jj
do pp = 1, size(LR(kk)%Li)
call destroy(LR(kk)%Li(pp))
end do
deallocate(LR(kk)%Li)
end do
hessen(jj + 1, jj) = sqrt(norm(Rs, errst=errst))
!if(prop_error('krylov_arnoldi_qmpoc : norm (3) failed.', &
! errst=errst)) return
! Exponentiate matrix
!call Exponentiate(prop(1:jj, 1:jj), deltat, hessen(1:jj, 1:jj))
allocate(prop(jj, jj), hessendt(jj, jj))
hessendt = hessen(1:jj, 1:jj)
call expm(prop, deltat, hessendt, jj, errst=errst)
!if(prop_error('krylov_arnoldi_MPOM_TYPE: '//&
! 'expm failed.', 'KrylovOps_include.f90:574', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * hessen(jj, jj - 1) * abs(prop(jj, 1)))
if(abs(estimate) < Cp%lanczos_tol) exit
propsave(1:jj - 1) = prop(1:jj - 1, 1) - propsave(1:jj - 1)
estimate = real(dot_product(propsave(1:jj - 1), &
propsave(1:jj - 1)), KIND=rKind)
!if(estimate < Cp%psi_tol) exit
if(abs(hessen(jj + 1, jj)) < Cp%Psi_tol) exit
propsave(1:jj) = prop(1:jj, 1)
! Avoid increment to jj = nkrylov + 1
if(jj == nkrylov) exit
deallocate(prop, hessendt)
end do
! Get the final propagated state
! ------------------------------
call destroy(Rs)
! location of largest and 2nd largest
kk = maxloc(abs(prop(1:jj, 1)), 1)
! Initialize variational ansatz to be the largest weight MPS
call copy(Psi, Krylovuv(kk))
prop(1:jj, 1) = prop(1:jj, 1) * beta0 !!! NEW
call fitmps(Psi, prop(1:jj, 1), Krylovuv, Cp, err, converged, &
errst=errst)
!if(prop_error('krylov_arnoldi_qmpoc : fitmps failed.', &
! errst=errst)) return
cerr = cerr + err
if((.not. converged) .and. (verbose > 0)) write(slog, *) &
'krylov_arnoldi : failed on fitting with error ', err
do kk = 1, jj
call destroy(Krylovuv(kk))
end do
deallocate(Krylovuv, LR, hessen, prop, propsave, hessendt)
end subroutine krylov_arnoldi_qmpoc
"""
return
[docs]def krylov_full_local_tensor_complex():
"""
fortran-subroutine - September 2017 (dj, updated)
Propagate 1 or 2 site Hamiltonian with LR-overlaps and Krylov.
**Arguments**
Psi : TYPE(tensorc), inout
Tensor representing one or two sites of the MPS. Tensors
is evolved under the MPO plus left-right-overlap.
LL : , inout
Left overlap of MPS and MPO.
Wm : TYPE(sr_matrix_tensor), inout
The MPO matrix for one or two sites represented as sparse
matrix.
RR : TYPE(tensorlistc), inout
Right overlap of MPS and MPO.
deltat : complex, inout
Time step for the evolution. Complex for real time evolution
and real for imaginary time evolution.
leftmost : LOGICAL, inout
Flag if subroutine is acting on the leftmost site.
rightmost : LOGICAL, inout
Flag if subroutine is acting on the rightmost site.
Cp : TYPE(ConvParam), inout
Convergence parameters for the algorithm.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
This version uses the Krylov-Lanczos algorithm for
hermitian Hamiltonians building the tridiagonal matrix in the
Krylov subspace.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_full_local_tensor_complex(Psi, LL, Wn, RR, deltat, &
leftmost, rightmost, Cp, pbc, errst)
type(tensorc), intent(inout) :: Psi
type(tensorlistc), intent(inout) :: LL, RR
type(sr_matrix_tensor), intent(inout) :: Wn
complex(KIND=rKind), intent(in) :: deltat
logical, intent(in) :: leftmost, rightmost
type(ConvParam), intent(in) :: Cp
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk, pp
! diagonal and 1st offdiagonal of tridiagonal upper hessenberg
real(KIND=rKind), allocatable :: beta(:), alpha(:)
! number of Krylov vectors and Krylov basis
integer :: nkrylov
type(tensorc), pointer :: Krylov(:)
! MPS at tensor w/ orthogonality center to be modified in
! iteration process
type(tensorc) :: Rs
! The exponential in the Krylov basis, including the previous step
complex(KIND=rKind), allocatable :: prop(:, :)
REAL(KIND=rKind) :: tol, estimate
complex(KIND=rKind) :: overlap
INTEGER :: k, n,p,q, chunk, i,j,top
LOGICAL :: exitSwitch
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_full_local_MPO_TYPE_complex:'//&
! 'Krylov does not have PBC.', 99, 'KrylovOps_include.f90:751', &
! errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylov(nkrylov), alpha(nkrylov), beta(0:nkrylov))
call copy(Rs, Psi)
beta(0) = sqrt(norm(Rs))
! First iteration
! ----------------
jj = 1
call pointto(Krylov(jj), Rs)!, scalar=zone / beta(jj - 1))
call scale(zone / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo_dot_psi(Rs, LL, Wn, RR, Krylov(jj), leftmost, rightmost)
alpha(jj) = real(dot(Krylov(jj), Rs), KIND=rKind)
call gaxpy(Rs, -alpha(jj), Krylov(jj))
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs)
call gaxpy(Rs, -overlap, Krylov(ii))
end do
beta(jj) = sqrt(norm(Rs))
if(verbose .ge. 2) write(slog, *) 'beta', jj, beta(jj)
if(beta(1) < Cp%lanczos_tol) then
call scale(exp(deltat * alpha(1)), Psi)
call destroy(Rs)
call destroy(Krylov(1))
deallocate(Krylov, alpha, beta)
return
end if
estimate = beta(1)
! Further iterations to build Krylov basis
! ----------------------------------------
do jj = 2, nkrylov
if(verbose > 2) write(slog, *) 'jj', jj, nkrylov, estimate
call pointto(Krylov(jj), Rs)!, scalar=zone / beta(jj - 1))
call scale(zone / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo_dot_psi(Rs, LL, Wn, RR, Krylov(jj), leftmost, rightmost)
alpha(jj) = real(dot(Krylov(jj), Rs), KIND=rKind)
call gaxpy(Rs, -alpha(jj), Krylov(jj))
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs)
call gaxpy(Rs, -overlap, Krylov(ii))
end do
beta(jj) = sqrt(norm(Rs))
! Exponentiate tridiagonal matrix
allocate(prop(jj, jj))
call expmh(prop, deltat, alpha(1:jj), beta(1:jj - 1), jj, &
errst=errst)
!if(prop_error('krylov_full_local_tensor_complex: '//&
! 'expmh failed.', 'KrylovOps_include.f90:826', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * beta(jj - 1) * abs(prop(jj, 1)))
if(verbose.ge.2) write(slog, *) 'beta', jj, beta(jj), estimate, &
abs(prop(jj, 1))
deallocate(prop)
if(abs(estimate) <= Cp%lanczos_tol) exit
if(abs(beta(jj)) <= 10.0_rKind * numzero) exit
end do
top = min(jj, nkrylov)
allocate(prop(top, top))
call expmh(prop, deltat, alpha(1:top), beta(1:top-1), top, errst=errst)
!if(prop_error('krylov_full_local_tensor_complex: '//&
! 'expmh failed.', 'KrylovOps_include.f90:846', &
! errst=errst)) then
! return
!end if
call scale(0.0_rKind, Psi)
do ii = 1, top
call gaxpy(Psi, prop(ii, 1), Krylov(ii))
call destroy(Krylov(ii))
end do
call destroy(Rs)
deallocate(Krylov, alpha, beta, prop)
end subroutine krylov_full_local_tensor_complex
"""
return
[docs]def krylov_full_local_tensorc_complex():
"""
fortran-subroutine - September 2017 (dj, updated)
Propagate 1 or 2 site Hamiltonian with LR-overlaps and Krylov.
**Arguments**
Psi : TYPE(tensorc), inout
Tensor representing one or two sites of the MPS. Tensors
is evolved under the MPO plus left-right-overlap.
LL : , inout
Left overlap of MPS and MPO.
Wm : TYPE(sr_matrix_tensorc), inout
The MPO matrix for one or two sites represented as sparse
matrix.
RR : TYPE(tensorlistc), inout
Right overlap of MPS and MPO.
deltat : complex, inout
Time step for the evolution. Complex for real time evolution
and real for imaginary time evolution.
leftmost : LOGICAL, inout
Flag if subroutine is acting on the leftmost site.
rightmost : LOGICAL, inout
Flag if subroutine is acting on the rightmost site.
Cp : TYPE(ConvParam), inout
Convergence parameters for the algorithm.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
This version uses the Krylov-Lanczos algorithm for
hermitian Hamiltonians building the tridiagonal matrix in the
Krylov subspace.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_full_local_tensorc_complex(Psi, LL, Wn, RR, deltat, &
leftmost, rightmost, Cp, pbc, errst)
type(tensorc), intent(inout) :: Psi
type(tensorlistc), intent(inout) :: LL, RR
type(sr_matrix_tensorc), intent(inout) :: Wn
complex(KIND=rKind), intent(in) :: deltat
logical, intent(in) :: leftmost, rightmost
type(ConvParam), intent(in) :: Cp
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk, pp
! diagonal and 1st offdiagonal of tridiagonal upper hessenberg
real(KIND=rKind), allocatable :: beta(:), alpha(:)
! number of Krylov vectors and Krylov basis
integer :: nkrylov
type(tensorc), pointer :: Krylov(:)
! MPS at tensor w/ orthogonality center to be modified in
! iteration process
type(tensorc) :: Rs
! The exponential in the Krylov basis, including the previous step
complex(KIND=rKind), allocatable :: prop(:, :)
REAL(KIND=rKind) :: tol, estimate
complex(KIND=rKind) :: overlap
INTEGER :: k, n,p,q, chunk, i,j,top
LOGICAL :: exitSwitch
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_full_local_MPO_TYPE_complex:'//&
! 'Krylov does not have PBC.', 99, 'KrylovOps_include.f90:751', &
! errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylov(nkrylov), alpha(nkrylov), beta(0:nkrylov))
call copy(Rs, Psi)
beta(0) = sqrt(norm(Rs))
! First iteration
! ----------------
jj = 1
call pointto(Krylov(jj), Rs)!, scalar=zone / beta(jj - 1))
call scale(zone / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo_dot_psi(Rs, LL, Wn, RR, Krylov(jj), leftmost, rightmost)
alpha(jj) = real(dot(Krylov(jj), Rs), KIND=rKind)
call gaxpy(Rs, -alpha(jj), Krylov(jj))
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs)
call gaxpy(Rs, -overlap, Krylov(ii))
end do
beta(jj) = sqrt(norm(Rs))
if(verbose .ge. 2) write(slog, *) 'beta', jj, beta(jj)
if(beta(1) < Cp%lanczos_tol) then
call scale(exp(deltat * alpha(1)), Psi)
call destroy(Rs)
call destroy(Krylov(1))
deallocate(Krylov, alpha, beta)
return
end if
estimate = beta(1)
! Further iterations to build Krylov basis
! ----------------------------------------
do jj = 2, nkrylov
if(verbose > 2) write(slog, *) 'jj', jj, nkrylov, estimate
call pointto(Krylov(jj), Rs)!, scalar=zone / beta(jj - 1))
call scale(zone / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo_dot_psi(Rs, LL, Wn, RR, Krylov(jj), leftmost, rightmost)
alpha(jj) = real(dot(Krylov(jj), Rs), KIND=rKind)
call gaxpy(Rs, -alpha(jj), Krylov(jj))
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs)
call gaxpy(Rs, -overlap, Krylov(ii))
end do
beta(jj) = sqrt(norm(Rs))
! Exponentiate tridiagonal matrix
allocate(prop(jj, jj))
call expmh(prop, deltat, alpha(1:jj), beta(1:jj - 1), jj, &
errst=errst)
!if(prop_error('krylov_full_local_tensorc_complex: '//&
! 'expmh failed.', 'KrylovOps_include.f90:826', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * beta(jj - 1) * abs(prop(jj, 1)))
if(verbose.ge.2) write(slog, *) 'beta', jj, beta(jj), estimate, &
abs(prop(jj, 1))
deallocate(prop)
if(abs(estimate) <= Cp%lanczos_tol) exit
if(abs(beta(jj)) <= 10.0_rKind * numzero) exit
end do
top = min(jj, nkrylov)
allocate(prop(top, top))
call expmh(prop, deltat, alpha(1:top), beta(1:top-1), top, errst=errst)
!if(prop_error('krylov_full_local_tensorc_complex: '//&
! 'expmh failed.', 'KrylovOps_include.f90:846', &
! errst=errst)) then
! return
!end if
call scale(0.0_rKind, Psi)
do ii = 1, top
call gaxpy(Psi, prop(ii, 1), Krylov(ii))
call destroy(Krylov(ii))
end do
call destroy(Rs)
deallocate(Krylov, alpha, beta, prop)
end subroutine krylov_full_local_tensorc_complex
"""
return
[docs]def krylov_full_local_qtensor_complex():
"""
fortran-subroutine - September 2017 (dj, updated)
Propagate 1 or 2 site Hamiltonian with LR-overlaps and Krylov.
**Arguments**
Psi : TYPE(qtensorc), inout
Tensor representing one or two sites of the MPS. Tensors
is evolved under the MPO plus left-right-overlap.
LL : , inout
Left overlap of MPS and MPO.
Wm : TYPE(sr_matrix_qtensor), inout
The MPO matrix for one or two sites represented as sparse
matrix.
RR : TYPE(qtensorclist), inout
Right overlap of MPS and MPO.
deltat : complex, inout
Time step for the evolution. Complex for real time evolution
and real for imaginary time evolution.
leftmost : LOGICAL, inout
Flag if subroutine is acting on the leftmost site.
rightmost : LOGICAL, inout
Flag if subroutine is acting on the rightmost site.
Cp : TYPE(ConvParam), inout
Convergence parameters for the algorithm.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
This version uses the Krylov-Lanczos algorithm for
hermitian Hamiltonians building the tridiagonal matrix in the
Krylov subspace.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_full_local_qtensor_complex(Psi, LL, Wn, RR, deltat, &
leftmost, rightmost, Cp, pbc, errst)
type(qtensorc), intent(inout) :: Psi
type(qtensorclist), intent(inout) :: LL, RR
type(sr_matrix_qtensor), intent(inout) :: Wn
complex(KIND=rKind), intent(in) :: deltat
logical, intent(in) :: leftmost, rightmost
type(ConvParam), intent(in) :: Cp
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk, pp
! diagonal and 1st offdiagonal of tridiagonal upper hessenberg
real(KIND=rKind), allocatable :: beta(:), alpha(:)
! number of Krylov vectors and Krylov basis
integer :: nkrylov
type(qtensorc), pointer :: Krylov(:)
! MPS at tensor w/ orthogonality center to be modified in
! iteration process
type(qtensorc) :: Rs
! The exponential in the Krylov basis, including the previous step
complex(KIND=rKind), allocatable :: prop(:, :)
REAL(KIND=rKind) :: tol, estimate
complex(KIND=rKind) :: overlap
INTEGER :: k, n,p,q, chunk, i,j,top
LOGICAL :: exitSwitch
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_full_local_MPO_TYPE_complex:'//&
! 'Krylov does not have PBC.', 99, 'KrylovOps_include.f90:751', &
! errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylov(nkrylov), alpha(nkrylov), beta(0:nkrylov))
call copy(Rs, Psi)
beta(0) = sqrt(norm(Rs))
! First iteration
! ----------------
jj = 1
call pointto(Krylov(jj), Rs)!, scalar=zone / beta(jj - 1))
call scale(zone / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo_dot_psi(Rs, LL, Wn, RR, Krylov(jj), leftmost, rightmost)
alpha(jj) = real(dot(Krylov(jj), Rs), KIND=rKind)
call gaxpy(Rs, -alpha(jj), Krylov(jj))
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs)
call gaxpy(Rs, -overlap, Krylov(ii))
end do
beta(jj) = sqrt(norm(Rs))
if(verbose .ge. 2) write(slog, *) 'beta', jj, beta(jj)
if(beta(1) < Cp%lanczos_tol) then
call scale(exp(deltat * alpha(1)), Psi)
call destroy(Rs)
call destroy(Krylov(1))
deallocate(Krylov, alpha, beta)
return
end if
estimate = beta(1)
! Further iterations to build Krylov basis
! ----------------------------------------
do jj = 2, nkrylov
if(verbose > 2) write(slog, *) 'jj', jj, nkrylov, estimate
call pointto(Krylov(jj), Rs)!, scalar=zone / beta(jj - 1))
call scale(zone / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo_dot_psi(Rs, LL, Wn, RR, Krylov(jj), leftmost, rightmost)
alpha(jj) = real(dot(Krylov(jj), Rs), KIND=rKind)
call gaxpy(Rs, -alpha(jj), Krylov(jj))
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs)
call gaxpy(Rs, -overlap, Krylov(ii))
end do
beta(jj) = sqrt(norm(Rs))
! Exponentiate tridiagonal matrix
allocate(prop(jj, jj))
call expmh(prop, deltat, alpha(1:jj), beta(1:jj - 1), jj, &
errst=errst)
!if(prop_error('krylov_full_local_qtensor_complex: '//&
! 'expmh failed.', 'KrylovOps_include.f90:826', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * beta(jj - 1) * abs(prop(jj, 1)))
if(verbose.ge.2) write(slog, *) 'beta', jj, beta(jj), estimate, &
abs(prop(jj, 1))
deallocate(prop)
if(abs(estimate) <= Cp%lanczos_tol) exit
if(abs(beta(jj)) <= 10.0_rKind * numzero) exit
end do
top = min(jj, nkrylov)
allocate(prop(top, top))
call expmh(prop, deltat, alpha(1:top), beta(1:top-1), top, errst=errst)
!if(prop_error('krylov_full_local_qtensor_complex: '//&
! 'expmh failed.', 'KrylovOps_include.f90:846', &
! errst=errst)) then
! return
!end if
call scale(0.0_rKind, Psi)
do ii = 1, top
call gaxpy(Psi, prop(ii, 1), Krylov(ii))
call destroy(Krylov(ii))
end do
call destroy(Rs)
deallocate(Krylov, alpha, beta, prop)
end subroutine krylov_full_local_qtensor_complex
"""
return
[docs]def krylov_full_local_qtensorc_complex():
"""
fortran-subroutine - September 2017 (dj, updated)
Propagate 1 or 2 site Hamiltonian with LR-overlaps and Krylov.
**Arguments**
Psi : TYPE(qtensorc), inout
Tensor representing one or two sites of the MPS. Tensors
is evolved under the MPO plus left-right-overlap.
LL : , inout
Left overlap of MPS and MPO.
Wm : TYPE(sr_matrix_qtensorc), inout
The MPO matrix for one or two sites represented as sparse
matrix.
RR : TYPE(qtensorclist), inout
Right overlap of MPS and MPO.
deltat : complex, inout
Time step for the evolution. Complex for real time evolution
and real for imaginary time evolution.
leftmost : LOGICAL, inout
Flag if subroutine is acting on the leftmost site.
rightmost : LOGICAL, inout
Flag if subroutine is acting on the rightmost site.
Cp : TYPE(ConvParam), inout
Convergence parameters for the algorithm.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
This version uses the Krylov-Lanczos algorithm for
hermitian Hamiltonians building the tridiagonal matrix in the
Krylov subspace.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_full_local_qtensorc_complex(Psi, LL, Wn, RR, deltat, &
leftmost, rightmost, Cp, pbc, errst)
type(qtensorc), intent(inout) :: Psi
type(qtensorclist), intent(inout) :: LL, RR
type(sr_matrix_qtensorc), intent(inout) :: Wn
complex(KIND=rKind), intent(in) :: deltat
logical, intent(in) :: leftmost, rightmost
type(ConvParam), intent(in) :: Cp
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk, pp
! diagonal and 1st offdiagonal of tridiagonal upper hessenberg
real(KIND=rKind), allocatable :: beta(:), alpha(:)
! number of Krylov vectors and Krylov basis
integer :: nkrylov
type(qtensorc), pointer :: Krylov(:)
! MPS at tensor w/ orthogonality center to be modified in
! iteration process
type(qtensorc) :: Rs
! The exponential in the Krylov basis, including the previous step
complex(KIND=rKind), allocatable :: prop(:, :)
REAL(KIND=rKind) :: tol, estimate
complex(KIND=rKind) :: overlap
INTEGER :: k, n,p,q, chunk, i,j,top
LOGICAL :: exitSwitch
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_full_local_MPO_TYPE_complex:'//&
! 'Krylov does not have PBC.', 99, 'KrylovOps_include.f90:751', &
! errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylov(nkrylov), alpha(nkrylov), beta(0:nkrylov))
call copy(Rs, Psi)
beta(0) = sqrt(norm(Rs))
! First iteration
! ----------------
jj = 1
call pointto(Krylov(jj), Rs)!, scalar=zone / beta(jj - 1))
call scale(zone / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo_dot_psi(Rs, LL, Wn, RR, Krylov(jj), leftmost, rightmost)
alpha(jj) = real(dot(Krylov(jj), Rs), KIND=rKind)
call gaxpy(Rs, -alpha(jj), Krylov(jj))
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs)
call gaxpy(Rs, -overlap, Krylov(ii))
end do
beta(jj) = sqrt(norm(Rs))
if(verbose .ge. 2) write(slog, *) 'beta', jj, beta(jj)
if(beta(1) < Cp%lanczos_tol) then
call scale(exp(deltat * alpha(1)), Psi)
call destroy(Rs)
call destroy(Krylov(1))
deallocate(Krylov, alpha, beta)
return
end if
estimate = beta(1)
! Further iterations to build Krylov basis
! ----------------------------------------
do jj = 2, nkrylov
if(verbose > 2) write(slog, *) 'jj', jj, nkrylov, estimate
call pointto(Krylov(jj), Rs)!, scalar=zone / beta(jj - 1))
call scale(zone / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo_dot_psi(Rs, LL, Wn, RR, Krylov(jj), leftmost, rightmost)
alpha(jj) = real(dot(Krylov(jj), Rs), KIND=rKind)
call gaxpy(Rs, -alpha(jj), Krylov(jj))
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs)
call gaxpy(Rs, -overlap, Krylov(ii))
end do
beta(jj) = sqrt(norm(Rs))
! Exponentiate tridiagonal matrix
allocate(prop(jj, jj))
call expmh(prop, deltat, alpha(1:jj), beta(1:jj - 1), jj, &
errst=errst)
!if(prop_error('krylov_full_local_qtensorc_complex: '//&
! 'expmh failed.', 'KrylovOps_include.f90:826', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * beta(jj - 1) * abs(prop(jj, 1)))
if(verbose.ge.2) write(slog, *) 'beta', jj, beta(jj), estimate, &
abs(prop(jj, 1))
deallocate(prop)
if(abs(estimate) <= Cp%lanczos_tol) exit
if(abs(beta(jj)) <= 10.0_rKind * numzero) exit
end do
top = min(jj, nkrylov)
allocate(prop(top, top))
call expmh(prop, deltat, alpha(1:top), beta(1:top-1), top, errst=errst)
!if(prop_error('krylov_full_local_qtensorc_complex: '//&
! 'expmh failed.', 'KrylovOps_include.f90:846', &
! errst=errst)) then
! return
!end if
call scale(0.0_rKind, Psi)
do ii = 1, top
call gaxpy(Psi, prop(ii, 1), Krylov(ii))
call destroy(Krylov(ii))
end do
call destroy(Rs)
deallocate(Krylov, alpha, beta, prop)
end subroutine krylov_full_local_qtensorc_complex
"""
return
[docs]def krylov_full_local_tensor_real():
"""
fortran-subroutine - September 2017 (dj, updated)
Propagate 1 or 2 site Hamiltonian with LR-overlaps and Krylov.
**Arguments**
Psi : TYPE(tensor), inout
Tensor representing one or two sites of the MPS. Tensors
is evolved under the MPO plus left-right-overlap.
LL : , inout
Left overlap of MPS and MPO.
Wm : TYPE(sr_matrix_tensor), inout
The MPO matrix for one or two sites represented as sparse
matrix.
RR : TYPE(tensorlist), inout
Right overlap of MPS and MPO.
deltat : real, inout
Time step for the evolution. Complex for real time evolution
and real for imaginary time evolution.
leftmost : LOGICAL, inout
Flag if subroutine is acting on the leftmost site.
rightmost : LOGICAL, inout
Flag if subroutine is acting on the rightmost site.
Cp : TYPE(ConvParam), inout
Convergence parameters for the algorithm.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
This version uses the Krylov-Lanczos algorithm for
hermitian Hamiltonians building the tridiagonal matrix in the
Krylov subspace.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_full_local_tensor_real(Psi, LL, Wn, RR, deltat, &
leftmost, rightmost, Cp, pbc, errst)
type(tensor), intent(inout) :: Psi
type(tensorlist), intent(inout) :: LL, RR
type(sr_matrix_tensor), intent(inout) :: Wn
real(KIND=rKind), intent(in) :: deltat
logical, intent(in) :: leftmost, rightmost
type(ConvParam), intent(in) :: Cp
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk, pp
! diagonal and 1st offdiagonal of tridiagonal upper hessenberg
real(KIND=rKind), allocatable :: beta(:), alpha(:)
! number of Krylov vectors and Krylov basis
integer :: nkrylov
type(tensor), pointer :: Krylov(:)
! MPS at tensor w/ orthogonality center to be modified in
! iteration process
type(tensor) :: Rs
! The exponential in the Krylov basis, including the previous step
real(KIND=rKind), allocatable :: prop(:, :)
REAL(KIND=rKind) :: tol, estimate
real(KIND=rKind) :: overlap
INTEGER :: k, n,p,q, chunk, i,j,top
LOGICAL :: exitSwitch
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_full_local_MPO_TYPE_real:'//&
! 'Krylov does not have PBC.', 99, 'KrylovOps_include.f90:751', &
! errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylov(nkrylov), alpha(nkrylov), beta(0:nkrylov))
call copy(Rs, Psi)
beta(0) = sqrt(norm(Rs))
! First iteration
! ----------------
jj = 1
call pointto(Krylov(jj), Rs)!, scalar=done / beta(jj - 1))
call scale(done / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo_dot_psi(Rs, LL, Wn, RR, Krylov(jj), leftmost, rightmost)
alpha(jj) = real(dot(Krylov(jj), Rs), KIND=rKind)
call gaxpy(Rs, -alpha(jj), Krylov(jj))
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs)
call gaxpy(Rs, -overlap, Krylov(ii))
end do
beta(jj) = sqrt(norm(Rs))
if(verbose .ge. 2) write(slog, *) 'beta', jj, beta(jj)
if(beta(1) < Cp%lanczos_tol) then
call scale(exp(deltat * alpha(1)), Psi)
call destroy(Rs)
call destroy(Krylov(1))
deallocate(Krylov, alpha, beta)
return
end if
estimate = beta(1)
! Further iterations to build Krylov basis
! ----------------------------------------
do jj = 2, nkrylov
if(verbose > 2) write(slog, *) 'jj', jj, nkrylov, estimate
call pointto(Krylov(jj), Rs)!, scalar=done / beta(jj - 1))
call scale(done / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo_dot_psi(Rs, LL, Wn, RR, Krylov(jj), leftmost, rightmost)
alpha(jj) = real(dot(Krylov(jj), Rs), KIND=rKind)
call gaxpy(Rs, -alpha(jj), Krylov(jj))
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs)
call gaxpy(Rs, -overlap, Krylov(ii))
end do
beta(jj) = sqrt(norm(Rs))
! Exponentiate tridiagonal matrix
allocate(prop(jj, jj))
call expmh(prop, deltat, alpha(1:jj), beta(1:jj - 1), jj, &
errst=errst)
!if(prop_error('krylov_full_local_tensor_real: '//&
! 'expmh failed.', 'KrylovOps_include.f90:826', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * beta(jj - 1) * abs(prop(jj, 1)))
if(verbose.ge.2) write(slog, *) 'beta', jj, beta(jj), estimate, &
abs(prop(jj, 1))
deallocate(prop)
if(abs(estimate) <= Cp%lanczos_tol) exit
if(abs(beta(jj)) <= 10.0_rKind * numzero) exit
end do
top = min(jj, nkrylov)
allocate(prop(top, top))
call expmh(prop, deltat, alpha(1:top), beta(1:top-1), top, errst=errst)
!if(prop_error('krylov_full_local_tensor_real: '//&
! 'expmh failed.', 'KrylovOps_include.f90:846', &
! errst=errst)) then
! return
!end if
call scale(0.0_rKind, Psi)
do ii = 1, top
call gaxpy(Psi, prop(ii, 1), Krylov(ii))
call destroy(Krylov(ii))
end do
call destroy(Rs)
deallocate(Krylov, alpha, beta, prop)
end subroutine krylov_full_local_tensor_real
"""
return
[docs]def krylov_full_local_tensorc_real():
"""
fortran-subroutine - September 2017 (dj, updated)
Propagate 1 or 2 site Hamiltonian with LR-overlaps and Krylov.
**Arguments**
Psi : TYPE(tensorc), inout
Tensor representing one or two sites of the MPS. Tensors
is evolved under the MPO plus left-right-overlap.
LL : , inout
Left overlap of MPS and MPO.
Wm : TYPE(sr_matrix_tensorc), inout
The MPO matrix for one or two sites represented as sparse
matrix.
RR : TYPE(tensorlistc), inout
Right overlap of MPS and MPO.
deltat : real, inout
Time step for the evolution. Complex for real time evolution
and real for imaginary time evolution.
leftmost : LOGICAL, inout
Flag if subroutine is acting on the leftmost site.
rightmost : LOGICAL, inout
Flag if subroutine is acting on the rightmost site.
Cp : TYPE(ConvParam), inout
Convergence parameters for the algorithm.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
This version uses the Krylov-Lanczos algorithm for
hermitian Hamiltonians building the tridiagonal matrix in the
Krylov subspace.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_full_local_tensorc_real(Psi, LL, Wn, RR, deltat, &
leftmost, rightmost, Cp, pbc, errst)
type(tensorc), intent(inout) :: Psi
type(tensorlistc), intent(inout) :: LL, RR
type(sr_matrix_tensorc), intent(inout) :: Wn
real(KIND=rKind), intent(in) :: deltat
logical, intent(in) :: leftmost, rightmost
type(ConvParam), intent(in) :: Cp
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk, pp
! diagonal and 1st offdiagonal of tridiagonal upper hessenberg
real(KIND=rKind), allocatable :: beta(:), alpha(:)
! number of Krylov vectors and Krylov basis
integer :: nkrylov
type(tensorc), pointer :: Krylov(:)
! MPS at tensor w/ orthogonality center to be modified in
! iteration process
type(tensorc) :: Rs
! The exponential in the Krylov basis, including the previous step
real(KIND=rKind), allocatable :: prop(:, :)
REAL(KIND=rKind) :: tol, estimate
complex(KIND=rKind) :: overlap
INTEGER :: k, n,p,q, chunk, i,j,top
LOGICAL :: exitSwitch
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_full_local_MPO_TYPE_real:'//&
! 'Krylov does not have PBC.', 99, 'KrylovOps_include.f90:751', &
! errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylov(nkrylov), alpha(nkrylov), beta(0:nkrylov))
call copy(Rs, Psi)
beta(0) = sqrt(norm(Rs))
! First iteration
! ----------------
jj = 1
call pointto(Krylov(jj), Rs)!, scalar=zone / beta(jj - 1))
call scale(zone / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo_dot_psi(Rs, LL, Wn, RR, Krylov(jj), leftmost, rightmost)
alpha(jj) = real(dot(Krylov(jj), Rs), KIND=rKind)
call gaxpy(Rs, -alpha(jj), Krylov(jj))
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs)
call gaxpy(Rs, -overlap, Krylov(ii))
end do
beta(jj) = sqrt(norm(Rs))
if(verbose .ge. 2) write(slog, *) 'beta', jj, beta(jj)
if(beta(1) < Cp%lanczos_tol) then
call scale(exp(deltat * alpha(1)), Psi)
call destroy(Rs)
call destroy(Krylov(1))
deallocate(Krylov, alpha, beta)
return
end if
estimate = beta(1)
! Further iterations to build Krylov basis
! ----------------------------------------
do jj = 2, nkrylov
if(verbose > 2) write(slog, *) 'jj', jj, nkrylov, estimate
call pointto(Krylov(jj), Rs)!, scalar=zone / beta(jj - 1))
call scale(zone / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo_dot_psi(Rs, LL, Wn, RR, Krylov(jj), leftmost, rightmost)
alpha(jj) = real(dot(Krylov(jj), Rs), KIND=rKind)
call gaxpy(Rs, -alpha(jj), Krylov(jj))
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs)
call gaxpy(Rs, -overlap, Krylov(ii))
end do
beta(jj) = sqrt(norm(Rs))
! Exponentiate tridiagonal matrix
allocate(prop(jj, jj))
call expmh(prop, deltat, alpha(1:jj), beta(1:jj - 1), jj, &
errst=errst)
!if(prop_error('krylov_full_local_tensorc_real: '//&
! 'expmh failed.', 'KrylovOps_include.f90:826', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * beta(jj - 1) * abs(prop(jj, 1)))
if(verbose.ge.2) write(slog, *) 'beta', jj, beta(jj), estimate, &
abs(prop(jj, 1))
deallocate(prop)
if(abs(estimate) <= Cp%lanczos_tol) exit
if(abs(beta(jj)) <= 10.0_rKind * numzero) exit
end do
top = min(jj, nkrylov)
allocate(prop(top, top))
call expmh(prop, deltat, alpha(1:top), beta(1:top-1), top, errst=errst)
!if(prop_error('krylov_full_local_tensorc_real: '//&
! 'expmh failed.', 'KrylovOps_include.f90:846', &
! errst=errst)) then
! return
!end if
call scale(0.0_rKind, Psi)
do ii = 1, top
call gaxpy(Psi, prop(ii, 1), Krylov(ii))
call destroy(Krylov(ii))
end do
call destroy(Rs)
deallocate(Krylov, alpha, beta, prop)
end subroutine krylov_full_local_tensorc_real
"""
return
[docs]def krylov_full_local_qtensor_real():
"""
fortran-subroutine - September 2017 (dj, updated)
Propagate 1 or 2 site Hamiltonian with LR-overlaps and Krylov.
**Arguments**
Psi : TYPE(qtensor), inout
Tensor representing one or two sites of the MPS. Tensors
is evolved under the MPO plus left-right-overlap.
LL : , inout
Left overlap of MPS and MPO.
Wm : TYPE(sr_matrix_qtensor), inout
The MPO matrix for one or two sites represented as sparse
matrix.
RR : TYPE(qtensorlist), inout
Right overlap of MPS and MPO.
deltat : real, inout
Time step for the evolution. Complex for real time evolution
and real for imaginary time evolution.
leftmost : LOGICAL, inout
Flag if subroutine is acting on the leftmost site.
rightmost : LOGICAL, inout
Flag if subroutine is acting on the rightmost site.
Cp : TYPE(ConvParam), inout
Convergence parameters for the algorithm.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
This version uses the Krylov-Lanczos algorithm for
hermitian Hamiltonians building the tridiagonal matrix in the
Krylov subspace.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_full_local_qtensor_real(Psi, LL, Wn, RR, deltat, &
leftmost, rightmost, Cp, pbc, errst)
type(qtensor), intent(inout) :: Psi
type(qtensorlist), intent(inout) :: LL, RR
type(sr_matrix_qtensor), intent(inout) :: Wn
real(KIND=rKind), intent(in) :: deltat
logical, intent(in) :: leftmost, rightmost
type(ConvParam), intent(in) :: Cp
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk, pp
! diagonal and 1st offdiagonal of tridiagonal upper hessenberg
real(KIND=rKind), allocatable :: beta(:), alpha(:)
! number of Krylov vectors and Krylov basis
integer :: nkrylov
type(qtensor), pointer :: Krylov(:)
! MPS at tensor w/ orthogonality center to be modified in
! iteration process
type(qtensor) :: Rs
! The exponential in the Krylov basis, including the previous step
real(KIND=rKind), allocatable :: prop(:, :)
REAL(KIND=rKind) :: tol, estimate
real(KIND=rKind) :: overlap
INTEGER :: k, n,p,q, chunk, i,j,top
LOGICAL :: exitSwitch
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_full_local_MPO_TYPE_real:'//&
! 'Krylov does not have PBC.', 99, 'KrylovOps_include.f90:751', &
! errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylov(nkrylov), alpha(nkrylov), beta(0:nkrylov))
call copy(Rs, Psi)
beta(0) = sqrt(norm(Rs))
! First iteration
! ----------------
jj = 1
call pointto(Krylov(jj), Rs)!, scalar=done / beta(jj - 1))
call scale(done / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo_dot_psi(Rs, LL, Wn, RR, Krylov(jj), leftmost, rightmost)
alpha(jj) = real(dot(Krylov(jj), Rs), KIND=rKind)
call gaxpy(Rs, -alpha(jj), Krylov(jj))
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs)
call gaxpy(Rs, -overlap, Krylov(ii))
end do
beta(jj) = sqrt(norm(Rs))
if(verbose .ge. 2) write(slog, *) 'beta', jj, beta(jj)
if(beta(1) < Cp%lanczos_tol) then
call scale(exp(deltat * alpha(1)), Psi)
call destroy(Rs)
call destroy(Krylov(1))
deallocate(Krylov, alpha, beta)
return
end if
estimate = beta(1)
! Further iterations to build Krylov basis
! ----------------------------------------
do jj = 2, nkrylov
if(verbose > 2) write(slog, *) 'jj', jj, nkrylov, estimate
call pointto(Krylov(jj), Rs)!, scalar=done / beta(jj - 1))
call scale(done / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo_dot_psi(Rs, LL, Wn, RR, Krylov(jj), leftmost, rightmost)
alpha(jj) = real(dot(Krylov(jj), Rs), KIND=rKind)
call gaxpy(Rs, -alpha(jj), Krylov(jj))
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs)
call gaxpy(Rs, -overlap, Krylov(ii))
end do
beta(jj) = sqrt(norm(Rs))
! Exponentiate tridiagonal matrix
allocate(prop(jj, jj))
call expmh(prop, deltat, alpha(1:jj), beta(1:jj - 1), jj, &
errst=errst)
!if(prop_error('krylov_full_local_qtensor_real: '//&
! 'expmh failed.', 'KrylovOps_include.f90:826', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * beta(jj - 1) * abs(prop(jj, 1)))
if(verbose.ge.2) write(slog, *) 'beta', jj, beta(jj), estimate, &
abs(prop(jj, 1))
deallocate(prop)
if(abs(estimate) <= Cp%lanczos_tol) exit
if(abs(beta(jj)) <= 10.0_rKind * numzero) exit
end do
top = min(jj, nkrylov)
allocate(prop(top, top))
call expmh(prop, deltat, alpha(1:top), beta(1:top-1), top, errst=errst)
!if(prop_error('krylov_full_local_qtensor_real: '//&
! 'expmh failed.', 'KrylovOps_include.f90:846', &
! errst=errst)) then
! return
!end if
call scale(0.0_rKind, Psi)
do ii = 1, top
call gaxpy(Psi, prop(ii, 1), Krylov(ii))
call destroy(Krylov(ii))
end do
call destroy(Rs)
deallocate(Krylov, alpha, beta, prop)
end subroutine krylov_full_local_qtensor_real
"""
return
[docs]def krylov_full_local_qtensorc_real():
"""
fortran-subroutine - September 2017 (dj, updated)
Propagate 1 or 2 site Hamiltonian with LR-overlaps and Krylov.
**Arguments**
Psi : TYPE(qtensorc), inout
Tensor representing one or two sites of the MPS. Tensors
is evolved under the MPO plus left-right-overlap.
LL : , inout
Left overlap of MPS and MPO.
Wm : TYPE(sr_matrix_qtensorc), inout
The MPO matrix for one or two sites represented as sparse
matrix.
RR : TYPE(qtensorclist), inout
Right overlap of MPS and MPO.
deltat : real, inout
Time step for the evolution. Complex for real time evolution
and real for imaginary time evolution.
leftmost : LOGICAL, inout
Flag if subroutine is acting on the leftmost site.
rightmost : LOGICAL, inout
Flag if subroutine is acting on the rightmost site.
Cp : TYPE(ConvParam), inout
Convergence parameters for the algorithm.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
This version uses the Krylov-Lanczos algorithm for
hermitian Hamiltonians building the tridiagonal matrix in the
Krylov subspace.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_full_local_qtensorc_real(Psi, LL, Wn, RR, deltat, &
leftmost, rightmost, Cp, pbc, errst)
type(qtensorc), intent(inout) :: Psi
type(qtensorclist), intent(inout) :: LL, RR
type(sr_matrix_qtensorc), intent(inout) :: Wn
real(KIND=rKind), intent(in) :: deltat
logical, intent(in) :: leftmost, rightmost
type(ConvParam), intent(in) :: Cp
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk, pp
! diagonal and 1st offdiagonal of tridiagonal upper hessenberg
real(KIND=rKind), allocatable :: beta(:), alpha(:)
! number of Krylov vectors and Krylov basis
integer :: nkrylov
type(qtensorc), pointer :: Krylov(:)
! MPS at tensor w/ orthogonality center to be modified in
! iteration process
type(qtensorc) :: Rs
! The exponential in the Krylov basis, including the previous step
real(KIND=rKind), allocatable :: prop(:, :)
REAL(KIND=rKind) :: tol, estimate
complex(KIND=rKind) :: overlap
INTEGER :: k, n,p,q, chunk, i,j,top
LOGICAL :: exitSwitch
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_full_local_MPO_TYPE_real:'//&
! 'Krylov does not have PBC.', 99, 'KrylovOps_include.f90:751', &
! errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylov(nkrylov), alpha(nkrylov), beta(0:nkrylov))
call copy(Rs, Psi)
beta(0) = sqrt(norm(Rs))
! First iteration
! ----------------
jj = 1
call pointto(Krylov(jj), Rs)!, scalar=zone / beta(jj - 1))
call scale(zone / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo_dot_psi(Rs, LL, Wn, RR, Krylov(jj), leftmost, rightmost)
alpha(jj) = real(dot(Krylov(jj), Rs), KIND=rKind)
call gaxpy(Rs, -alpha(jj), Krylov(jj))
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs)
call gaxpy(Rs, -overlap, Krylov(ii))
end do
beta(jj) = sqrt(norm(Rs))
if(verbose .ge. 2) write(slog, *) 'beta', jj, beta(jj)
if(beta(1) < Cp%lanczos_tol) then
call scale(exp(deltat * alpha(1)), Psi)
call destroy(Rs)
call destroy(Krylov(1))
deallocate(Krylov, alpha, beta)
return
end if
estimate = beta(1)
! Further iterations to build Krylov basis
! ----------------------------------------
do jj = 2, nkrylov
if(verbose > 2) write(slog, *) 'jj', jj, nkrylov, estimate
call pointto(Krylov(jj), Rs)!, scalar=zone / beta(jj - 1))
call scale(zone / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo_dot_psi(Rs, LL, Wn, RR, Krylov(jj), leftmost, rightmost)
alpha(jj) = real(dot(Krylov(jj), Rs), KIND=rKind)
call gaxpy(Rs, -alpha(jj), Krylov(jj))
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs)
call gaxpy(Rs, -overlap, Krylov(ii))
end do
beta(jj) = sqrt(norm(Rs))
! Exponentiate tridiagonal matrix
allocate(prop(jj, jj))
call expmh(prop, deltat, alpha(1:jj), beta(1:jj - 1), jj, &
errst=errst)
!if(prop_error('krylov_full_local_qtensorc_real: '//&
! 'expmh failed.', 'KrylovOps_include.f90:826', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * beta(jj - 1) * abs(prop(jj, 1)))
if(verbose.ge.2) write(slog, *) 'beta', jj, beta(jj), estimate, &
abs(prop(jj, 1))
deallocate(prop)
if(abs(estimate) <= Cp%lanczos_tol) exit
if(abs(beta(jj)) <= 10.0_rKind * numzero) exit
end do
top = min(jj, nkrylov)
allocate(prop(top, top))
call expmh(prop, deltat, alpha(1:top), beta(1:top-1), top, errst=errst)
!if(prop_error('krylov_full_local_qtensorc_real: '//&
! 'expmh failed.', 'KrylovOps_include.f90:846', &
! errst=errst)) then
! return
!end if
call scale(0.0_rKind, Psi)
do ii = 1, top
call gaxpy(Psi, prop(ii, 1), Krylov(ii))
call destroy(Krylov(ii))
end do
call destroy(Rs)
deallocate(Krylov, alpha, beta, prop)
end subroutine krylov_full_local_qtensorc_real
"""
return
[docs]def krylov_arnoldi_full_local_tensor():
"""
fortran-subroutine - November 2017 (dj)
Propagate 1 or 2 site Hamiltonian with LR-overlaps and Krylov-Arnoldi.
**Arguments**
Psi : TYPE(tensorc), inout
Tensor representing one or two sites of the MPS. Tensors
is evolved under the MPO plus left-right-overlap.
LL : , inout
Left overlap of MPS and MPO.
Wm : TYPE(sr_matrix_tensor), inout
The MPO matrix for one or two sites represented as sparse
matrix.
RR : TYPE(tensorlistc), inout
Right overlap of MPS and MPO.
deltat : DT_TYPE, inout
Time step for the evolution. Complex for real time evolution
and real for imaginary time evolution.
leftmost : LOGICAL, inout
Flag if subroutine is acting on the leftmost site.
rightmost : LOGICAL, inout
Flag if subroutine is acting on the rightmost site.
Cp : TYPE(ConvParam), inout
Convergence parameters for the algorithm.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
This version uses the Krylov-Arnoldi algorithm for
general operators building the upper Hessenberg matrix in the
Krylov subspace.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_arnoldi_full_local_tensor(Psi, LL, Wn, RR, deltat, &
leftmost, rightmost, Cp, pbc, errst)
type(tensorc), intent(inout) :: Psi
type(tensorlistc), intent(inout) :: LL, RR
type(sr_matrix_tensor), intent(inout) :: Wn
complex(KIND=rKind), intent(in) :: deltat
logical, intent(in) :: leftmost, rightmost
type(ConvParam), intent(in) :: Cp
logical, intent(in) :: pbc
integer, intent(inout), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk
! number of Krylov vectors / Krylov basis
integer :: nkrylov
type(tensorc), pointer :: Krylov(:)
! Upper Hessenberg matrix / copy of it
complex(KIND=rKind), dimension(:, :), allocatable :: hessen, hessencopy
! initial norm
real(KIND=rKind) :: beta0
real(KIND=rKind) :: estimate
! vector for Krylov algorithm - modification carried out on r
type(tensorc) :: Rs, Tmp
! propagator
complex(KIND=rKind), allocatable :: prop(:, :)
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_arnoldi_full_local_MPO_TYPE:'//&
! 'Krylov does not have PBC.', 99, 'KrylovOps_include.f90:984', &
! errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
allocate(Krylov(nkrylov), hessen(nkrylov + 1, nkrylov + 1))
hessen = 0.0_rKind
call copy(Rs, Psi)
beta0 = sqrt(norm(Rs))
! First data point - initialize values
! ------------------------------------
jj = 1
call pointto(Krylov(jj), Rs)!, scalar=zone / beta0)
call scale(zone / beta0, Krylov(jj))
!call destroy(Rs)
call mpo_dot_psi(Rs, LL, Wn, RR, Krylov(jj), leftmost, rightmost)
hessen(jj + 1, jj) = sqrt(norm(Rs))
if(verbose >=2) write(slog, *) 'beta', jj, hessen(jj + 1, jj)
if(real(hessen(jj + 1, jj), KIND=rKind) < Cp%lanczos_tol) then
call scale(exp(deltat * hessen(1 , 1)), Psi)
call destroy(Rs)
call destroy(Krylov(1))
deallocate(Krylov, hessen)
return
end if
estimate = abs(hessen(2, 1))
! Iteration over vectors 2, ...
! -----------------------------
do jj = 2, nkrylov
if(verbose >= 2) write(slog, *) 'jj', jj, nkrylov, estimate
call pointto(Krylov(jj), Rs)!, scalar=zone / hessen(jj, jj -1))
call scale(zone / hessen(jj, jj -1), Krylov(jj))
!call destroy(Rs)
call mpo_dot_psi(Rs, LL, Wn, RR, Krylov(jj), leftmost, rightmost)
do kk = 1, (jj - 1)
hessen(kk, jj) = dot(Rs, Krylov(kk))
call gaxpy(Rs, - hessen(kk, jj), Krylov(kk))
end do
hessen(jj + 1, jj) = sqrt(norm(Rs))
! Exponentiate matrix
allocate(hessencopy(jj, jj), prop(jj, jj))
hessencopy = hessen(:jj, :jj)
call expm(prop, deltat, hessencopy, jj, errst=errst)
!if(prop_error('krylov_arnoldi_full_local_tensor: '//&
! 'expm failed.', 'KrylovOps_include.f90:1048', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * hessen(jj, jj - 1) * abs(prop(jj, 1)))
if(verbose >= 2) write(slog, *) 'beta', jj, hessen(jj + 1, jj), &
estimate, abs(prop(jj, 1))
if(abs(estimate) <= Cp%lanczos_tol) exit
if(abs(hessen(jj + 1, jj)) <= 10.0_rKind * numzero) exit
! Avoid increment to jj = nkrylov + 1
if(jj == nkrylov) exit
deallocate(hessencopy, prop)
end do
if(jj == nKrylov) write(slog, *) 'full two-site Krylov-Anroldi '//&
'failed to converge!', hessen(jj + 1, 1), estimate
prop(:, 1) = prop(:, 1) * beta0
call scale(0.0_rKind, Psi)
do ii = 1, jj
call gaxpy(Psi, prop(ii, 1), Krylov(ii))
call destroy(Krylov(ii))
end do
call destroy(Rs)
deallocate(Krylov, hessen, hessencopy, prop)
end subroutine krylov_arnoldi_full_local_tensor
"""
return
[docs]def krylov_arnoldi_full_local_tensorc():
"""
fortran-subroutine - November 2017 (dj)
Propagate 1 or 2 site Hamiltonian with LR-overlaps and Krylov-Arnoldi.
**Arguments**
Psi : TYPE(tensorc), inout
Tensor representing one or two sites of the MPS. Tensors
is evolved under the MPO plus left-right-overlap.
LL : , inout
Left overlap of MPS and MPO.
Wm : TYPE(sr_matrix_tensorc), inout
The MPO matrix for one or two sites represented as sparse
matrix.
RR : TYPE(tensorlistc), inout
Right overlap of MPS and MPO.
deltat : DT_TYPE, inout
Time step for the evolution. Complex for real time evolution
and real for imaginary time evolution.
leftmost : LOGICAL, inout
Flag if subroutine is acting on the leftmost site.
rightmost : LOGICAL, inout
Flag if subroutine is acting on the rightmost site.
Cp : TYPE(ConvParam), inout
Convergence parameters for the algorithm.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
This version uses the Krylov-Arnoldi algorithm for
general operators building the upper Hessenberg matrix in the
Krylov subspace.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_arnoldi_full_local_tensorc(Psi, LL, Wn, RR, deltat, &
leftmost, rightmost, Cp, pbc, errst)
type(tensorc), intent(inout) :: Psi
type(tensorlistc), intent(inout) :: LL, RR
type(sr_matrix_tensorc), intent(inout) :: Wn
complex(KIND=rKind), intent(in) :: deltat
logical, intent(in) :: leftmost, rightmost
type(ConvParam), intent(in) :: Cp
logical, intent(in) :: pbc
integer, intent(inout), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk
! number of Krylov vectors / Krylov basis
integer :: nkrylov
type(tensorc), pointer :: Krylov(:)
! Upper Hessenberg matrix / copy of it
complex(KIND=rKind), dimension(:, :), allocatable :: hessen, hessencopy
! initial norm
real(KIND=rKind) :: beta0
real(KIND=rKind) :: estimate
! vector for Krylov algorithm - modification carried out on r
type(tensorc) :: Rs, Tmp
! propagator
complex(KIND=rKind), allocatable :: prop(:, :)
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_arnoldi_full_local_MPO_TYPE:'//&
! 'Krylov does not have PBC.', 99, 'KrylovOps_include.f90:984', &
! errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
allocate(Krylov(nkrylov), hessen(nkrylov + 1, nkrylov + 1))
hessen = 0.0_rKind
call copy(Rs, Psi)
beta0 = sqrt(norm(Rs))
! First data point - initialize values
! ------------------------------------
jj = 1
call pointto(Krylov(jj), Rs)!, scalar=zone / beta0)
call scale(zone / beta0, Krylov(jj))
!call destroy(Rs)
call mpo_dot_psi(Rs, LL, Wn, RR, Krylov(jj), leftmost, rightmost)
hessen(jj + 1, jj) = sqrt(norm(Rs))
if(verbose >=2) write(slog, *) 'beta', jj, hessen(jj + 1, jj)
if(real(hessen(jj + 1, jj), KIND=rKind) < Cp%lanczos_tol) then
call scale(exp(deltat * hessen(1 , 1)), Psi)
call destroy(Rs)
call destroy(Krylov(1))
deallocate(Krylov, hessen)
return
end if
estimate = abs(hessen(2, 1))
! Iteration over vectors 2, ...
! -----------------------------
do jj = 2, nkrylov
if(verbose >= 2) write(slog, *) 'jj', jj, nkrylov, estimate
call pointto(Krylov(jj), Rs)!, scalar=zone / hessen(jj, jj -1))
call scale(zone / hessen(jj, jj -1), Krylov(jj))
!call destroy(Rs)
call mpo_dot_psi(Rs, LL, Wn, RR, Krylov(jj), leftmost, rightmost)
do kk = 1, (jj - 1)
hessen(kk, jj) = dot(Rs, Krylov(kk))
call gaxpy(Rs, - hessen(kk, jj), Krylov(kk))
end do
hessen(jj + 1, jj) = sqrt(norm(Rs))
! Exponentiate matrix
allocate(hessencopy(jj, jj), prop(jj, jj))
hessencopy = hessen(:jj, :jj)
call expm(prop, deltat, hessencopy, jj, errst=errst)
!if(prop_error('krylov_arnoldi_full_local_tensorc: '//&
! 'expm failed.', 'KrylovOps_include.f90:1048', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * hessen(jj, jj - 1) * abs(prop(jj, 1)))
if(verbose >= 2) write(slog, *) 'beta', jj, hessen(jj + 1, jj), &
estimate, abs(prop(jj, 1))
if(abs(estimate) <= Cp%lanczos_tol) exit
if(abs(hessen(jj + 1, jj)) <= 10.0_rKind * numzero) exit
! Avoid increment to jj = nkrylov + 1
if(jj == nkrylov) exit
deallocate(hessencopy, prop)
end do
if(jj == nKrylov) write(slog, *) 'full two-site Krylov-Anroldi '//&
'failed to converge!', hessen(jj + 1, 1), estimate
prop(:, 1) = prop(:, 1) * beta0
call scale(0.0_rKind, Psi)
do ii = 1, jj
call gaxpy(Psi, prop(ii, 1), Krylov(ii))
call destroy(Krylov(ii))
end do
call destroy(Rs)
deallocate(Krylov, hessen, hessencopy, prop)
end subroutine krylov_arnoldi_full_local_tensorc
"""
return
[docs]def krylov_arnoldi_full_local_qtensor():
"""
fortran-subroutine - November 2017 (dj)
Propagate 1 or 2 site Hamiltonian with LR-overlaps and Krylov-Arnoldi.
**Arguments**
Psi : TYPE(qtensorc), inout
Tensor representing one or two sites of the MPS. Tensors
is evolved under the MPO plus left-right-overlap.
LL : , inout
Left overlap of MPS and MPO.
Wm : TYPE(sr_matrix_qtensor), inout
The MPO matrix for one or two sites represented as sparse
matrix.
RR : TYPE(qtensorclist), inout
Right overlap of MPS and MPO.
deltat : DT_TYPE, inout
Time step for the evolution. Complex for real time evolution
and real for imaginary time evolution.
leftmost : LOGICAL, inout
Flag if subroutine is acting on the leftmost site.
rightmost : LOGICAL, inout
Flag if subroutine is acting on the rightmost site.
Cp : TYPE(ConvParam), inout
Convergence parameters for the algorithm.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
This version uses the Krylov-Arnoldi algorithm for
general operators building the upper Hessenberg matrix in the
Krylov subspace.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_arnoldi_full_local_qtensor(Psi, LL, Wn, RR, deltat, &
leftmost, rightmost, Cp, pbc, errst)
type(qtensorc), intent(inout) :: Psi
type(qtensorclist), intent(inout) :: LL, RR
type(sr_matrix_qtensor), intent(inout) :: Wn
complex(KIND=rKind), intent(in) :: deltat
logical, intent(in) :: leftmost, rightmost
type(ConvParam), intent(in) :: Cp
logical, intent(in) :: pbc
integer, intent(inout), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk
! number of Krylov vectors / Krylov basis
integer :: nkrylov
type(qtensorc), pointer :: Krylov(:)
! Upper Hessenberg matrix / copy of it
complex(KIND=rKind), dimension(:, :), allocatable :: hessen, hessencopy
! initial norm
real(KIND=rKind) :: beta0
real(KIND=rKind) :: estimate
! vector for Krylov algorithm - modification carried out on r
type(qtensorc) :: Rs, Tmp
! propagator
complex(KIND=rKind), allocatable :: prop(:, :)
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_arnoldi_full_local_MPO_TYPE:'//&
! 'Krylov does not have PBC.', 99, 'KrylovOps_include.f90:984', &
! errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
allocate(Krylov(nkrylov), hessen(nkrylov + 1, nkrylov + 1))
hessen = 0.0_rKind
call copy(Rs, Psi)
beta0 = sqrt(norm(Rs))
! First data point - initialize values
! ------------------------------------
jj = 1
call pointto(Krylov(jj), Rs)!, scalar=zone / beta0)
call scale(zone / beta0, Krylov(jj))
!call destroy(Rs)
call mpo_dot_psi(Rs, LL, Wn, RR, Krylov(jj), leftmost, rightmost)
hessen(jj + 1, jj) = sqrt(norm(Rs))
if(verbose >=2) write(slog, *) 'beta', jj, hessen(jj + 1, jj)
if(real(hessen(jj + 1, jj), KIND=rKind) < Cp%lanczos_tol) then
call scale(exp(deltat * hessen(1 , 1)), Psi)
call destroy(Rs)
call destroy(Krylov(1))
deallocate(Krylov, hessen)
return
end if
estimate = abs(hessen(2, 1))
! Iteration over vectors 2, ...
! -----------------------------
do jj = 2, nkrylov
if(verbose >= 2) write(slog, *) 'jj', jj, nkrylov, estimate
call pointto(Krylov(jj), Rs)!, scalar=zone / hessen(jj, jj -1))
call scale(zone / hessen(jj, jj -1), Krylov(jj))
!call destroy(Rs)
call mpo_dot_psi(Rs, LL, Wn, RR, Krylov(jj), leftmost, rightmost)
do kk = 1, (jj - 1)
hessen(kk, jj) = dot(Rs, Krylov(kk))
call gaxpy(Rs, - hessen(kk, jj), Krylov(kk))
end do
hessen(jj + 1, jj) = sqrt(norm(Rs))
! Exponentiate matrix
allocate(hessencopy(jj, jj), prop(jj, jj))
hessencopy = hessen(:jj, :jj)
call expm(prop, deltat, hessencopy, jj, errst=errst)
!if(prop_error('krylov_arnoldi_full_local_qtensor: '//&
! 'expm failed.', 'KrylovOps_include.f90:1048', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * hessen(jj, jj - 1) * abs(prop(jj, 1)))
if(verbose >= 2) write(slog, *) 'beta', jj, hessen(jj + 1, jj), &
estimate, abs(prop(jj, 1))
if(abs(estimate) <= Cp%lanczos_tol) exit
if(abs(hessen(jj + 1, jj)) <= 10.0_rKind * numzero) exit
! Avoid increment to jj = nkrylov + 1
if(jj == nkrylov) exit
deallocate(hessencopy, prop)
end do
if(jj == nKrylov) write(slog, *) 'full two-site Krylov-Anroldi '//&
'failed to converge!', hessen(jj + 1, 1), estimate
prop(:, 1) = prop(:, 1) * beta0
call scale(0.0_rKind, Psi)
do ii = 1, jj
call gaxpy(Psi, prop(ii, 1), Krylov(ii))
call destroy(Krylov(ii))
end do
call destroy(Rs)
deallocate(Krylov, hessen, hessencopy, prop)
end subroutine krylov_arnoldi_full_local_qtensor
"""
return
[docs]def krylov_arnoldi_full_local_qtensorc():
"""
fortran-subroutine - November 2017 (dj)
Propagate 1 or 2 site Hamiltonian with LR-overlaps and Krylov-Arnoldi.
**Arguments**
Psi : TYPE(qtensorc), inout
Tensor representing one or two sites of the MPS. Tensors
is evolved under the MPO plus left-right-overlap.
LL : , inout
Left overlap of MPS and MPO.
Wm : TYPE(sr_matrix_qtensorc), inout
The MPO matrix for one or two sites represented as sparse
matrix.
RR : TYPE(qtensorclist), inout
Right overlap of MPS and MPO.
deltat : DT_TYPE, inout
Time step for the evolution. Complex for real time evolution
and real for imaginary time evolution.
leftmost : LOGICAL, inout
Flag if subroutine is acting on the leftmost site.
rightmost : LOGICAL, inout
Flag if subroutine is acting on the rightmost site.
Cp : TYPE(ConvParam), inout
Convergence parameters for the algorithm.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
This version uses the Krylov-Arnoldi algorithm for
general operators building the upper Hessenberg matrix in the
Krylov subspace.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_arnoldi_full_local_qtensorc(Psi, LL, Wn, RR, deltat, &
leftmost, rightmost, Cp, pbc, errst)
type(qtensorc), intent(inout) :: Psi
type(qtensorclist), intent(inout) :: LL, RR
type(sr_matrix_qtensorc), intent(inout) :: Wn
complex(KIND=rKind), intent(in) :: deltat
logical, intent(in) :: leftmost, rightmost
type(ConvParam), intent(in) :: Cp
logical, intent(in) :: pbc
integer, intent(inout), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk
! number of Krylov vectors / Krylov basis
integer :: nkrylov
type(qtensorc), pointer :: Krylov(:)
! Upper Hessenberg matrix / copy of it
complex(KIND=rKind), dimension(:, :), allocatable :: hessen, hessencopy
! initial norm
real(KIND=rKind) :: beta0
real(KIND=rKind) :: estimate
! vector for Krylov algorithm - modification carried out on r
type(qtensorc) :: Rs, Tmp
! propagator
complex(KIND=rKind), allocatable :: prop(:, :)
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_arnoldi_full_local_MPO_TYPE:'//&
! 'Krylov does not have PBC.', 99, 'KrylovOps_include.f90:984', &
! errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
allocate(Krylov(nkrylov), hessen(nkrylov + 1, nkrylov + 1))
hessen = 0.0_rKind
call copy(Rs, Psi)
beta0 = sqrt(norm(Rs))
! First data point - initialize values
! ------------------------------------
jj = 1
call pointto(Krylov(jj), Rs)!, scalar=zone / beta0)
call scale(zone / beta0, Krylov(jj))
!call destroy(Rs)
call mpo_dot_psi(Rs, LL, Wn, RR, Krylov(jj), leftmost, rightmost)
hessen(jj + 1, jj) = sqrt(norm(Rs))
if(verbose >=2) write(slog, *) 'beta', jj, hessen(jj + 1, jj)
if(real(hessen(jj + 1, jj), KIND=rKind) < Cp%lanczos_tol) then
call scale(exp(deltat * hessen(1 , 1)), Psi)
call destroy(Rs)
call destroy(Krylov(1))
deallocate(Krylov, hessen)
return
end if
estimate = abs(hessen(2, 1))
! Iteration over vectors 2, ...
! -----------------------------
do jj = 2, nkrylov
if(verbose >= 2) write(slog, *) 'jj', jj, nkrylov, estimate
call pointto(Krylov(jj), Rs)!, scalar=zone / hessen(jj, jj -1))
call scale(zone / hessen(jj, jj -1), Krylov(jj))
!call destroy(Rs)
call mpo_dot_psi(Rs, LL, Wn, RR, Krylov(jj), leftmost, rightmost)
do kk = 1, (jj - 1)
hessen(kk, jj) = dot(Rs, Krylov(kk))
call gaxpy(Rs, - hessen(kk, jj), Krylov(kk))
end do
hessen(jj + 1, jj) = sqrt(norm(Rs))
! Exponentiate matrix
allocate(hessencopy(jj, jj), prop(jj, jj))
hessencopy = hessen(:jj, :jj)
call expm(prop, deltat, hessencopy, jj, errst=errst)
!if(prop_error('krylov_arnoldi_full_local_qtensorc: '//&
! 'expm failed.', 'KrylovOps_include.f90:1048', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * hessen(jj, jj - 1) * abs(prop(jj, 1)))
if(verbose >= 2) write(slog, *) 'beta', jj, hessen(jj + 1, jj), &
estimate, abs(prop(jj, 1))
if(abs(estimate) <= Cp%lanczos_tol) exit
if(abs(hessen(jj + 1, jj)) <= 10.0_rKind * numzero) exit
! Avoid increment to jj = nkrylov + 1
if(jj == nkrylov) exit
deallocate(hessencopy, prop)
end do
if(jj == nKrylov) write(slog, *) 'full two-site Krylov-Anroldi '//&
'failed to converge!', hessen(jj + 1, 1), estimate
prop(:, 1) = prop(:, 1) * beta0
call scale(0.0_rKind, Psi)
do ii = 1, jj
call gaxpy(Psi, prop(ii, 1), Krylov(ii))
call destroy(Krylov(ii))
end do
call destroy(Rs)
deallocate(Krylov, hessen, hessencopy, prop)
end subroutine krylov_arnoldi_full_local_qtensorc
"""
return
[docs]def krylov_local_tensor_complex():
"""
fortran-subroutine - August 2017 (dj, updated)
Propagate a local two-site Hamiltonian with the Krylov algorithm. This
version does not consider the left-right overlap.
**Arguments**
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
This version uses the two-site Hamiltonian in its tensor
representation. Scaling is chi^2 d^4 for each Krylov step, and an
additional d^4 kappa to build the two site Hamiltonian for this
subroutine. Using the MPO representation, might speed things up
to 2 chi^2 d^2 kappa. (kappa bond dimension of the MPO.)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_local_tensor_complex(Psi, Ham, deltat, Cp, &
pbc, errst)
type(tensorc), intent(inout) :: Psi
type(tensor), intent(inout) :: Ham
complex(KIND=rKind), intent(in) :: deltat
type(ConvParam), intent(in) :: Cp
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj!, kk, pp
! diagonal and 1st offdiagonal of tridiagonal upper hessenberg
real(KIND=rKind), allocatable :: beta(:), alpha(:)
! number of Krylov vectors and Krylov basis
integer :: nkrylov
type(tensorc), pointer :: Krylov(:)
! number of Krylov vectors used
integer :: top
! MPS at tensor w/ orthogonality center to be modified in
! iteration process
type(tensorc) :: Rs!, Tmp
! The exponential in the Krylov basis, including the previous step
complex(KIND=rKind), allocatable :: prop(:, :)
REAL(KIND=rKind) :: estimate!, tol
complex(KIND=rKind) :: overlap
!INTEGER :: k, n,p,q, chunk, i,j
!LOGICAL :: exitSwitch
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_local_tensor_'//&
! 'complex: Krylov does not have PBC.', 99, &
! 'KrylovOps_include.f90:1158', errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylov(nkrylov), alpha(nkrylov), beta(0:nkrylov))
call copy(Rs, Psi, errst=errst)
!if(prop_error('krylov_local_tensor : copy (1) '//&
! 'failed.', errst=errst)) return
beta(0) = sqrt(norm(Rs))
! First iteration
! ---------------
jj = 1
call pointto(Krylov(jj), Rs)!, scalar=zone / beta(jj - 1), errst=errst)
!!if(prop_error('krylov_local_tensor : copy (2) '//&
!! 'failed.', errst=errst)) return
call scale(zone / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo2_dot_psi(Rs, Ham, Krylov(jj), errst=errst)
!if(prop_error('krylov_local_tensor : mpo2_dot_psi (1) '//&
! 'failed.', errst=errst)) return
alpha(jj) = real(dot(Krylov(jj), Rs, errst=errst), KIND=rKind)
!if(prop_error('krylov_local_tensor : dot (1) '//&
! 'failed.', errst=errst)) then
! return
!end if
call gaxpy(Rs, -alpha(jj), Krylov(jj))
!if(prop_error('krylov_local_tensor : gaxpy (1) '//&
! 'failed.', errst=errst)) return
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs, errst=errst)
!if(prop_error('krylov_local_tensor : dot (2) '//&
! 'failed.', errst=errst)) return
call gaxpy(Rs, -overlap, Krylov(ii), errst=errst)
!if(prop_error('krylov_local_tensor : gaxpy (2) '//&
! 'failed.', errst=errst)) return
end do
beta(jj) = sqrt(norm(Rs))
if(verbose >= 2) write(slog, *) 'beta', jj, beta(jj)
if(beta(1) < Cp%lanczos_tol) then
call scale(exp(deltat * alpha(1)), Psi)
call destroy(Rs)
call destroy(Krylov(1))
deallocate(Krylov, alpha, beta)
return
end if
estimate = beta(1)
! Further iterations to build Krylov basis
! ----------------------------------------
do jj = 2, nkrylov
if(verbose >= 2) write(slog, *) 'jj', jj, nkrylov, estimate
call pointto(Krylov(jj), Rs)!, scalar=zone / beta(jj - 1), errst=errst)
!if(prop_error('krylov_local_tensor : copy (4) '//&
! 'failed.', errst=errst)) return
call scale(zone / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo2_dot_psi(Rs, Ham, Krylov(jj), errst=errst)
!if(prop_error('krylov_local_tensor : mpo2_dot_psi (2) '//&
! 'failed.', errst=errst)) return
alpha(jj) = real(dot(Krylov(jj), Rs, errst=errst), KIND=rKind)
!if(prop_error('krylov_local_tensor : dot (3) '//&
! 'failed.', errst=errst)) return
call gaxpy(Rs, -alpha(jj), Krylov(jj), errst=errst)
!if(prop_error('krylov_local_tensor : gaxpy (3) '//&
! 'failed.', errst=errst)) return
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs, errst=errst)
!if(prop_error('krylov_local_tensor : dot (4) '//&
! 'failed.', errst=errst)) return
call gaxpy(Rs, -overlap, Krylov(ii), errst=errst)
!if(prop_error('krylov_local_tensor : gaxpy (4) '//&
! 'failed.', errst=errst)) return
end do
beta(jj) = sqrt(norm(Rs))
! Exponentiate tridiagonal matrix
allocate(prop(jj, jj))
call expmh(prop, deltat, alpha(1:jj), beta(1:jj-1), jj, errst=errst)
!if(prop_error('krylov_local_tensor_complex: '//&
! 'expmh failed.', 'KrylovOps_include.f90:1267', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * beta(jj - 1) * abs(prop(jj,1)))
if(verbose >= 2) write(slog, *) 'beta', jj, beta(jj), estimate, &
abs(prop(jj, 1))
deallocate(prop)
if(abs(estimate) <= Cp%lanczos_tol) exit
if(abs(beta(jj)) <= 10.0_rKind * numzero) exit
end do
top = min(jj, nkrylov)
if(top == nkrylov) write(slog, *) 'two-site Krylov failed to '//&
'converge!', beta(ubound(beta, 1)), estimate
allocate(prop(1:top, 1:top))
call expmh(prop, deltat, alpha(1:top), beta(1:top - 1), top, &
errst=errst)
!if(prop_error('krylov_local_tensor_complex: '//&
! 'expmh failed.', 'KrylovOps_include.f90:1290', &
! errst=errst)) then
! return
!end if
call scale(0.0_rKind, Psi)
do ii = 1, top
call gaxpy(Psi, prop(ii,1), Krylov(ii), errst=errst)
!if(prop_error('krylov_local_tensor : gaxpy (5) '//&
! 'failed.', errst=errst)) return
call destroy(Krylov(ii))
end do
call destroy(Rs)
deallocate(Krylov, alpha, beta, prop)
end subroutine krylov_local_tensor_complex
"""
return
[docs]def krylov_local_tensorc_complex():
"""
fortran-subroutine - August 2017 (dj, updated)
Propagate a local two-site Hamiltonian with the Krylov algorithm. This
version does not consider the left-right overlap.
**Arguments**
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
This version uses the two-site Hamiltonian in its tensor
representation. Scaling is chi^2 d^4 for each Krylov step, and an
additional d^4 kappa to build the two site Hamiltonian for this
subroutine. Using the MPO representation, might speed things up
to 2 chi^2 d^2 kappa. (kappa bond dimension of the MPO.)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_local_tensorc_complex(Psi, Ham, deltat, Cp, &
pbc, errst)
type(tensorc), intent(inout) :: Psi
type(tensorc), intent(inout) :: Ham
complex(KIND=rKind), intent(in) :: deltat
type(ConvParam), intent(in) :: Cp
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj!, kk, pp
! diagonal and 1st offdiagonal of tridiagonal upper hessenberg
real(KIND=rKind), allocatable :: beta(:), alpha(:)
! number of Krylov vectors and Krylov basis
integer :: nkrylov
type(tensorc), pointer :: Krylov(:)
! number of Krylov vectors used
integer :: top
! MPS at tensor w/ orthogonality center to be modified in
! iteration process
type(tensorc) :: Rs!, Tmp
! The exponential in the Krylov basis, including the previous step
complex(KIND=rKind), allocatable :: prop(:, :)
REAL(KIND=rKind) :: estimate!, tol
complex(KIND=rKind) :: overlap
!INTEGER :: k, n,p,q, chunk, i,j
!LOGICAL :: exitSwitch
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_local_tensorc_'//&
! 'complex: Krylov does not have PBC.', 99, &
! 'KrylovOps_include.f90:1158', errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylov(nkrylov), alpha(nkrylov), beta(0:nkrylov))
call copy(Rs, Psi, errst=errst)
!if(prop_error('krylov_local_tensorc : copy (1) '//&
! 'failed.', errst=errst)) return
beta(0) = sqrt(norm(Rs))
! First iteration
! ---------------
jj = 1
call pointto(Krylov(jj), Rs)!, scalar=zone / beta(jj - 1), errst=errst)
!!if(prop_error('krylov_local_tensorc : copy (2) '//&
!! 'failed.', errst=errst)) return
call scale(zone / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo2_dot_psi(Rs, Ham, Krylov(jj), errst=errst)
!if(prop_error('krylov_local_tensorc : mpo2_dot_psi (1) '//&
! 'failed.', errst=errst)) return
alpha(jj) = real(dot(Krylov(jj), Rs, errst=errst), KIND=rKind)
!if(prop_error('krylov_local_tensorc : dot (1) '//&
! 'failed.', errst=errst)) then
! return
!end if
call gaxpy(Rs, -alpha(jj), Krylov(jj))
!if(prop_error('krylov_local_tensorc : gaxpy (1) '//&
! 'failed.', errst=errst)) return
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs, errst=errst)
!if(prop_error('krylov_local_tensorc : dot (2) '//&
! 'failed.', errst=errst)) return
call gaxpy(Rs, -overlap, Krylov(ii), errst=errst)
!if(prop_error('krylov_local_tensorc : gaxpy (2) '//&
! 'failed.', errst=errst)) return
end do
beta(jj) = sqrt(norm(Rs))
if(verbose >= 2) write(slog, *) 'beta', jj, beta(jj)
if(beta(1) < Cp%lanczos_tol) then
call scale(exp(deltat * alpha(1)), Psi)
call destroy(Rs)
call destroy(Krylov(1))
deallocate(Krylov, alpha, beta)
return
end if
estimate = beta(1)
! Further iterations to build Krylov basis
! ----------------------------------------
do jj = 2, nkrylov
if(verbose >= 2) write(slog, *) 'jj', jj, nkrylov, estimate
call pointto(Krylov(jj), Rs)!, scalar=zone / beta(jj - 1), errst=errst)
!if(prop_error('krylov_local_tensorc : copy (4) '//&
! 'failed.', errst=errst)) return
call scale(zone / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo2_dot_psi(Rs, Ham, Krylov(jj), errst=errst)
!if(prop_error('krylov_local_tensorc : mpo2_dot_psi (2) '//&
! 'failed.', errst=errst)) return
alpha(jj) = real(dot(Krylov(jj), Rs, errst=errst), KIND=rKind)
!if(prop_error('krylov_local_tensorc : dot (3) '//&
! 'failed.', errst=errst)) return
call gaxpy(Rs, -alpha(jj), Krylov(jj), errst=errst)
!if(prop_error('krylov_local_tensorc : gaxpy (3) '//&
! 'failed.', errst=errst)) return
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs, errst=errst)
!if(prop_error('krylov_local_tensorc : dot (4) '//&
! 'failed.', errst=errst)) return
call gaxpy(Rs, -overlap, Krylov(ii), errst=errst)
!if(prop_error('krylov_local_tensorc : gaxpy (4) '//&
! 'failed.', errst=errst)) return
end do
beta(jj) = sqrt(norm(Rs))
! Exponentiate tridiagonal matrix
allocate(prop(jj, jj))
call expmh(prop, deltat, alpha(1:jj), beta(1:jj-1), jj, errst=errst)
!if(prop_error('krylov_local_tensorc_complex: '//&
! 'expmh failed.', 'KrylovOps_include.f90:1267', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * beta(jj - 1) * abs(prop(jj,1)))
if(verbose >= 2) write(slog, *) 'beta', jj, beta(jj), estimate, &
abs(prop(jj, 1))
deallocate(prop)
if(abs(estimate) <= Cp%lanczos_tol) exit
if(abs(beta(jj)) <= 10.0_rKind * numzero) exit
end do
top = min(jj, nkrylov)
if(top == nkrylov) write(slog, *) 'two-site Krylov failed to '//&
'converge!', beta(ubound(beta, 1)), estimate
allocate(prop(1:top, 1:top))
call expmh(prop, deltat, alpha(1:top), beta(1:top - 1), top, &
errst=errst)
!if(prop_error('krylov_local_tensorc_complex: '//&
! 'expmh failed.', 'KrylovOps_include.f90:1290', &
! errst=errst)) then
! return
!end if
call scale(0.0_rKind, Psi)
do ii = 1, top
call gaxpy(Psi, prop(ii,1), Krylov(ii), errst=errst)
!if(prop_error('krylov_local_tensorc : gaxpy (5) '//&
! 'failed.', errst=errst)) return
call destroy(Krylov(ii))
end do
call destroy(Rs)
deallocate(Krylov, alpha, beta, prop)
end subroutine krylov_local_tensorc_complex
"""
return
[docs]def krylov_local_qtensor_complex():
"""
fortran-subroutine - August 2017 (dj, updated)
Propagate a local two-site Hamiltonian with the Krylov algorithm. This
version does not consider the left-right overlap.
**Arguments**
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
This version uses the two-site Hamiltonian in its tensor
representation. Scaling is chi^2 d^4 for each Krylov step, and an
additional d^4 kappa to build the two site Hamiltonian for this
subroutine. Using the MPO representation, might speed things up
to 2 chi^2 d^2 kappa. (kappa bond dimension of the MPO.)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_local_qtensor_complex(Psi, Ham, deltat, Cp, &
pbc, errst)
type(qtensorc), intent(inout) :: Psi
type(qtensor), intent(inout) :: Ham
complex(KIND=rKind), intent(in) :: deltat
type(ConvParam), intent(in) :: Cp
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj!, kk, pp
! diagonal and 1st offdiagonal of tridiagonal upper hessenberg
real(KIND=rKind), allocatable :: beta(:), alpha(:)
! number of Krylov vectors and Krylov basis
integer :: nkrylov
type(qtensorc), pointer :: Krylov(:)
! number of Krylov vectors used
integer :: top
! MPS at tensor w/ orthogonality center to be modified in
! iteration process
type(qtensorc) :: Rs!, Tmp
! The exponential in the Krylov basis, including the previous step
complex(KIND=rKind), allocatable :: prop(:, :)
REAL(KIND=rKind) :: estimate!, tol
complex(KIND=rKind) :: overlap
!INTEGER :: k, n,p,q, chunk, i,j
!LOGICAL :: exitSwitch
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_local_qtensor_'//&
! 'complex: Krylov does not have PBC.', 99, &
! 'KrylovOps_include.f90:1158', errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylov(nkrylov), alpha(nkrylov), beta(0:nkrylov))
call copy(Rs, Psi, errst=errst)
!if(prop_error('krylov_local_qtensor : copy (1) '//&
! 'failed.', errst=errst)) return
beta(0) = sqrt(norm(Rs))
! First iteration
! ---------------
jj = 1
call pointto(Krylov(jj), Rs)!, scalar=zone / beta(jj - 1), errst=errst)
!!if(prop_error('krylov_local_qtensor : copy (2) '//&
!! 'failed.', errst=errst)) return
call scale(zone / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo2_dot_psi(Rs, Ham, Krylov(jj), errst=errst)
!if(prop_error('krylov_local_qtensor : mpo2_dot_psi (1) '//&
! 'failed.', errst=errst)) return
alpha(jj) = real(dot(Krylov(jj), Rs, errst=errst), KIND=rKind)
!if(prop_error('krylov_local_qtensor : dot (1) '//&
! 'failed.', errst=errst)) then
! return
!end if
call gaxpy(Rs, -alpha(jj), Krylov(jj))
!if(prop_error('krylov_local_qtensor : gaxpy (1) '//&
! 'failed.', errst=errst)) return
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs, errst=errst)
!if(prop_error('krylov_local_qtensor : dot (2) '//&
! 'failed.', errst=errst)) return
call gaxpy(Rs, -overlap, Krylov(ii), errst=errst)
!if(prop_error('krylov_local_qtensor : gaxpy (2) '//&
! 'failed.', errst=errst)) return
end do
beta(jj) = sqrt(norm(Rs))
if(verbose >= 2) write(slog, *) 'beta', jj, beta(jj)
if(beta(1) < Cp%lanczos_tol) then
call scale(exp(deltat * alpha(1)), Psi)
call destroy(Rs)
call destroy(Krylov(1))
deallocate(Krylov, alpha, beta)
return
end if
estimate = beta(1)
! Further iterations to build Krylov basis
! ----------------------------------------
do jj = 2, nkrylov
if(verbose >= 2) write(slog, *) 'jj', jj, nkrylov, estimate
call pointto(Krylov(jj), Rs)!, scalar=zone / beta(jj - 1), errst=errst)
!if(prop_error('krylov_local_qtensor : copy (4) '//&
! 'failed.', errst=errst)) return
call scale(zone / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo2_dot_psi(Rs, Ham, Krylov(jj), errst=errst)
!if(prop_error('krylov_local_qtensor : mpo2_dot_psi (2) '//&
! 'failed.', errst=errst)) return
alpha(jj) = real(dot(Krylov(jj), Rs, errst=errst), KIND=rKind)
!if(prop_error('krylov_local_qtensor : dot (3) '//&
! 'failed.', errst=errst)) return
call gaxpy(Rs, -alpha(jj), Krylov(jj), errst=errst)
!if(prop_error('krylov_local_qtensor : gaxpy (3) '//&
! 'failed.', errst=errst)) return
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs, errst=errst)
!if(prop_error('krylov_local_qtensor : dot (4) '//&
! 'failed.', errst=errst)) return
call gaxpy(Rs, -overlap, Krylov(ii), errst=errst)
!if(prop_error('krylov_local_qtensor : gaxpy (4) '//&
! 'failed.', errst=errst)) return
end do
beta(jj) = sqrt(norm(Rs))
! Exponentiate tridiagonal matrix
allocate(prop(jj, jj))
call expmh(prop, deltat, alpha(1:jj), beta(1:jj-1), jj, errst=errst)
!if(prop_error('krylov_local_qtensor_complex: '//&
! 'expmh failed.', 'KrylovOps_include.f90:1267', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * beta(jj - 1) * abs(prop(jj,1)))
if(verbose >= 2) write(slog, *) 'beta', jj, beta(jj), estimate, &
abs(prop(jj, 1))
deallocate(prop)
if(abs(estimate) <= Cp%lanczos_tol) exit
if(abs(beta(jj)) <= 10.0_rKind * numzero) exit
end do
top = min(jj, nkrylov)
if(top == nkrylov) write(slog, *) 'two-site Krylov failed to '//&
'converge!', beta(ubound(beta, 1)), estimate
allocate(prop(1:top, 1:top))
call expmh(prop, deltat, alpha(1:top), beta(1:top - 1), top, &
errst=errst)
!if(prop_error('krylov_local_qtensor_complex: '//&
! 'expmh failed.', 'KrylovOps_include.f90:1290', &
! errst=errst)) then
! return
!end if
call scale(0.0_rKind, Psi)
do ii = 1, top
call gaxpy(Psi, prop(ii,1), Krylov(ii), errst=errst)
!if(prop_error('krylov_local_qtensor : gaxpy (5) '//&
! 'failed.', errst=errst)) return
call destroy(Krylov(ii))
end do
call destroy(Rs)
deallocate(Krylov, alpha, beta, prop)
end subroutine krylov_local_qtensor_complex
"""
return
[docs]def krylov_local_qtensorc_complex():
"""
fortran-subroutine - August 2017 (dj, updated)
Propagate a local two-site Hamiltonian with the Krylov algorithm. This
version does not consider the left-right overlap.
**Arguments**
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
This version uses the two-site Hamiltonian in its tensor
representation. Scaling is chi^2 d^4 for each Krylov step, and an
additional d^4 kappa to build the two site Hamiltonian for this
subroutine. Using the MPO representation, might speed things up
to 2 chi^2 d^2 kappa. (kappa bond dimension of the MPO.)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_local_qtensorc_complex(Psi, Ham, deltat, Cp, &
pbc, errst)
type(qtensorc), intent(inout) :: Psi
type(qtensorc), intent(inout) :: Ham
complex(KIND=rKind), intent(in) :: deltat
type(ConvParam), intent(in) :: Cp
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj!, kk, pp
! diagonal and 1st offdiagonal of tridiagonal upper hessenberg
real(KIND=rKind), allocatable :: beta(:), alpha(:)
! number of Krylov vectors and Krylov basis
integer :: nkrylov
type(qtensorc), pointer :: Krylov(:)
! number of Krylov vectors used
integer :: top
! MPS at tensor w/ orthogonality center to be modified in
! iteration process
type(qtensorc) :: Rs!, Tmp
! The exponential in the Krylov basis, including the previous step
complex(KIND=rKind), allocatable :: prop(:, :)
REAL(KIND=rKind) :: estimate!, tol
complex(KIND=rKind) :: overlap
!INTEGER :: k, n,p,q, chunk, i,j
!LOGICAL :: exitSwitch
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_local_qtensorc_'//&
! 'complex: Krylov does not have PBC.', 99, &
! 'KrylovOps_include.f90:1158', errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylov(nkrylov), alpha(nkrylov), beta(0:nkrylov))
call copy(Rs, Psi, errst=errst)
!if(prop_error('krylov_local_qtensorc : copy (1) '//&
! 'failed.', errst=errst)) return
beta(0) = sqrt(norm(Rs))
! First iteration
! ---------------
jj = 1
call pointto(Krylov(jj), Rs)!, scalar=zone / beta(jj - 1), errst=errst)
!!if(prop_error('krylov_local_qtensorc : copy (2) '//&
!! 'failed.', errst=errst)) return
call scale(zone / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo2_dot_psi(Rs, Ham, Krylov(jj), errst=errst)
!if(prop_error('krylov_local_qtensorc : mpo2_dot_psi (1) '//&
! 'failed.', errst=errst)) return
alpha(jj) = real(dot(Krylov(jj), Rs, errst=errst), KIND=rKind)
!if(prop_error('krylov_local_qtensorc : dot (1) '//&
! 'failed.', errst=errst)) then
! return
!end if
call gaxpy(Rs, -alpha(jj), Krylov(jj))
!if(prop_error('krylov_local_qtensorc : gaxpy (1) '//&
! 'failed.', errst=errst)) return
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs, errst=errst)
!if(prop_error('krylov_local_qtensorc : dot (2) '//&
! 'failed.', errst=errst)) return
call gaxpy(Rs, -overlap, Krylov(ii), errst=errst)
!if(prop_error('krylov_local_qtensorc : gaxpy (2) '//&
! 'failed.', errst=errst)) return
end do
beta(jj) = sqrt(norm(Rs))
if(verbose >= 2) write(slog, *) 'beta', jj, beta(jj)
if(beta(1) < Cp%lanczos_tol) then
call scale(exp(deltat * alpha(1)), Psi)
call destroy(Rs)
call destroy(Krylov(1))
deallocate(Krylov, alpha, beta)
return
end if
estimate = beta(1)
! Further iterations to build Krylov basis
! ----------------------------------------
do jj = 2, nkrylov
if(verbose >= 2) write(slog, *) 'jj', jj, nkrylov, estimate
call pointto(Krylov(jj), Rs)!, scalar=zone / beta(jj - 1), errst=errst)
!if(prop_error('krylov_local_qtensorc : copy (4) '//&
! 'failed.', errst=errst)) return
call scale(zone / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo2_dot_psi(Rs, Ham, Krylov(jj), errst=errst)
!if(prop_error('krylov_local_qtensorc : mpo2_dot_psi (2) '//&
! 'failed.', errst=errst)) return
alpha(jj) = real(dot(Krylov(jj), Rs, errst=errst), KIND=rKind)
!if(prop_error('krylov_local_qtensorc : dot (3) '//&
! 'failed.', errst=errst)) return
call gaxpy(Rs, -alpha(jj), Krylov(jj), errst=errst)
!if(prop_error('krylov_local_qtensorc : gaxpy (3) '//&
! 'failed.', errst=errst)) return
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs, errst=errst)
!if(prop_error('krylov_local_qtensorc : dot (4) '//&
! 'failed.', errst=errst)) return
call gaxpy(Rs, -overlap, Krylov(ii), errst=errst)
!if(prop_error('krylov_local_qtensorc : gaxpy (4) '//&
! 'failed.', errst=errst)) return
end do
beta(jj) = sqrt(norm(Rs))
! Exponentiate tridiagonal matrix
allocate(prop(jj, jj))
call expmh(prop, deltat, alpha(1:jj), beta(1:jj-1), jj, errst=errst)
!if(prop_error('krylov_local_qtensorc_complex: '//&
! 'expmh failed.', 'KrylovOps_include.f90:1267', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * beta(jj - 1) * abs(prop(jj,1)))
if(verbose >= 2) write(slog, *) 'beta', jj, beta(jj), estimate, &
abs(prop(jj, 1))
deallocate(prop)
if(abs(estimate) <= Cp%lanczos_tol) exit
if(abs(beta(jj)) <= 10.0_rKind * numzero) exit
end do
top = min(jj, nkrylov)
if(top == nkrylov) write(slog, *) 'two-site Krylov failed to '//&
'converge!', beta(ubound(beta, 1)), estimate
allocate(prop(1:top, 1:top))
call expmh(prop, deltat, alpha(1:top), beta(1:top - 1), top, &
errst=errst)
!if(prop_error('krylov_local_qtensorc_complex: '//&
! 'expmh failed.', 'KrylovOps_include.f90:1290', &
! errst=errst)) then
! return
!end if
call scale(0.0_rKind, Psi)
do ii = 1, top
call gaxpy(Psi, prop(ii,1), Krylov(ii), errst=errst)
!if(prop_error('krylov_local_qtensorc : gaxpy (5) '//&
! 'failed.', errst=errst)) return
call destroy(Krylov(ii))
end do
call destroy(Rs)
deallocate(Krylov, alpha, beta, prop)
end subroutine krylov_local_qtensorc_complex
"""
return
[docs]def krylov_local_tensor_real():
"""
fortran-subroutine - August 2017 (dj, updated)
Propagate a local two-site Hamiltonian with the Krylov algorithm. This
version does not consider the left-right overlap.
**Arguments**
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
This version uses the two-site Hamiltonian in its tensor
representation. Scaling is chi^2 d^4 for each Krylov step, and an
additional d^4 kappa to build the two site Hamiltonian for this
subroutine. Using the MPO representation, might speed things up
to 2 chi^2 d^2 kappa. (kappa bond dimension of the MPO.)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_local_tensor_real(Psi, Ham, deltat, Cp, &
pbc, errst)
type(tensor), intent(inout) :: Psi
type(tensor), intent(inout) :: Ham
real(KIND=rKind), intent(in) :: deltat
type(ConvParam), intent(in) :: Cp
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj!, kk, pp
! diagonal and 1st offdiagonal of tridiagonal upper hessenberg
real(KIND=rKind), allocatable :: beta(:), alpha(:)
! number of Krylov vectors and Krylov basis
integer :: nkrylov
type(tensor), pointer :: Krylov(:)
! number of Krylov vectors used
integer :: top
! MPS at tensor w/ orthogonality center to be modified in
! iteration process
type(tensor) :: Rs!, Tmp
! The exponential in the Krylov basis, including the previous step
real(KIND=rKind), allocatable :: prop(:, :)
REAL(KIND=rKind) :: estimate!, tol
real(KIND=rKind) :: overlap
!INTEGER :: k, n,p,q, chunk, i,j
!LOGICAL :: exitSwitch
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_local_tensor_'//&
! 'real: Krylov does not have PBC.', 99, &
! 'KrylovOps_include.f90:1158', errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylov(nkrylov), alpha(nkrylov), beta(0:nkrylov))
call copy(Rs, Psi, errst=errst)
!if(prop_error('krylov_local_tensor : copy (1) '//&
! 'failed.', errst=errst)) return
beta(0) = sqrt(norm(Rs))
! First iteration
! ---------------
jj = 1
call pointto(Krylov(jj), Rs)!, scalar=done / beta(jj - 1), errst=errst)
!!if(prop_error('krylov_local_tensor : copy (2) '//&
!! 'failed.', errst=errst)) return
call scale(done / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo2_dot_psi(Rs, Ham, Krylov(jj), errst=errst)
!if(prop_error('krylov_local_tensor : mpo2_dot_psi (1) '//&
! 'failed.', errst=errst)) return
alpha(jj) = real(dot(Krylov(jj), Rs, errst=errst), KIND=rKind)
!if(prop_error('krylov_local_tensor : dot (1) '//&
! 'failed.', errst=errst)) then
! return
!end if
call gaxpy(Rs, -alpha(jj), Krylov(jj))
!if(prop_error('krylov_local_tensor : gaxpy (1) '//&
! 'failed.', errst=errst)) return
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs, errst=errst)
!if(prop_error('krylov_local_tensor : dot (2) '//&
! 'failed.', errst=errst)) return
call gaxpy(Rs, -overlap, Krylov(ii), errst=errst)
!if(prop_error('krylov_local_tensor : gaxpy (2) '//&
! 'failed.', errst=errst)) return
end do
beta(jj) = sqrt(norm(Rs))
if(verbose >= 2) write(slog, *) 'beta', jj, beta(jj)
if(beta(1) < Cp%lanczos_tol) then
call scale(exp(deltat * alpha(1)), Psi)
call destroy(Rs)
call destroy(Krylov(1))
deallocate(Krylov, alpha, beta)
return
end if
estimate = beta(1)
! Further iterations to build Krylov basis
! ----------------------------------------
do jj = 2, nkrylov
if(verbose >= 2) write(slog, *) 'jj', jj, nkrylov, estimate
call pointto(Krylov(jj), Rs)!, scalar=done / beta(jj - 1), errst=errst)
!if(prop_error('krylov_local_tensor : copy (4) '//&
! 'failed.', errst=errst)) return
call scale(done / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo2_dot_psi(Rs, Ham, Krylov(jj), errst=errst)
!if(prop_error('krylov_local_tensor : mpo2_dot_psi (2) '//&
! 'failed.', errst=errst)) return
alpha(jj) = real(dot(Krylov(jj), Rs, errst=errst), KIND=rKind)
!if(prop_error('krylov_local_tensor : dot (3) '//&
! 'failed.', errst=errst)) return
call gaxpy(Rs, -alpha(jj), Krylov(jj), errst=errst)
!if(prop_error('krylov_local_tensor : gaxpy (3) '//&
! 'failed.', errst=errst)) return
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs, errst=errst)
!if(prop_error('krylov_local_tensor : dot (4) '//&
! 'failed.', errst=errst)) return
call gaxpy(Rs, -overlap, Krylov(ii), errst=errst)
!if(prop_error('krylov_local_tensor : gaxpy (4) '//&
! 'failed.', errst=errst)) return
end do
beta(jj) = sqrt(norm(Rs))
! Exponentiate tridiagonal matrix
allocate(prop(jj, jj))
call expmh(prop, deltat, alpha(1:jj), beta(1:jj-1), jj, errst=errst)
!if(prop_error('krylov_local_tensor_real: '//&
! 'expmh failed.', 'KrylovOps_include.f90:1267', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * beta(jj - 1) * abs(prop(jj,1)))
if(verbose >= 2) write(slog, *) 'beta', jj, beta(jj), estimate, &
abs(prop(jj, 1))
deallocate(prop)
if(abs(estimate) <= Cp%lanczos_tol) exit
if(abs(beta(jj)) <= 10.0_rKind * numzero) exit
end do
top = min(jj, nkrylov)
if(top == nkrylov) write(slog, *) 'two-site Krylov failed to '//&
'converge!', beta(ubound(beta, 1)), estimate
allocate(prop(1:top, 1:top))
call expmh(prop, deltat, alpha(1:top), beta(1:top - 1), top, &
errst=errst)
!if(prop_error('krylov_local_tensor_real: '//&
! 'expmh failed.', 'KrylovOps_include.f90:1290', &
! errst=errst)) then
! return
!end if
call scale(0.0_rKind, Psi)
do ii = 1, top
call gaxpy(Psi, prop(ii,1), Krylov(ii), errst=errst)
!if(prop_error('krylov_local_tensor : gaxpy (5) '//&
! 'failed.', errst=errst)) return
call destroy(Krylov(ii))
end do
call destroy(Rs)
deallocate(Krylov, alpha, beta, prop)
end subroutine krylov_local_tensor_real
"""
return
[docs]def krylov_local_tensorc_real():
"""
fortran-subroutine - August 2017 (dj, updated)
Propagate a local two-site Hamiltonian with the Krylov algorithm. This
version does not consider the left-right overlap.
**Arguments**
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
This version uses the two-site Hamiltonian in its tensor
representation. Scaling is chi^2 d^4 for each Krylov step, and an
additional d^4 kappa to build the two site Hamiltonian for this
subroutine. Using the MPO representation, might speed things up
to 2 chi^2 d^2 kappa. (kappa bond dimension of the MPO.)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_local_tensorc_real(Psi, Ham, deltat, Cp, &
pbc, errst)
type(tensorc), intent(inout) :: Psi
type(tensorc), intent(inout) :: Ham
real(KIND=rKind), intent(in) :: deltat
type(ConvParam), intent(in) :: Cp
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj!, kk, pp
! diagonal and 1st offdiagonal of tridiagonal upper hessenberg
real(KIND=rKind), allocatable :: beta(:), alpha(:)
! number of Krylov vectors and Krylov basis
integer :: nkrylov
type(tensorc), pointer :: Krylov(:)
! number of Krylov vectors used
integer :: top
! MPS at tensor w/ orthogonality center to be modified in
! iteration process
type(tensorc) :: Rs!, Tmp
! The exponential in the Krylov basis, including the previous step
real(KIND=rKind), allocatable :: prop(:, :)
REAL(KIND=rKind) :: estimate!, tol
complex(KIND=rKind) :: overlap
!INTEGER :: k, n,p,q, chunk, i,j
!LOGICAL :: exitSwitch
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_local_tensorc_'//&
! 'real: Krylov does not have PBC.', 99, &
! 'KrylovOps_include.f90:1158', errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylov(nkrylov), alpha(nkrylov), beta(0:nkrylov))
call copy(Rs, Psi, errst=errst)
!if(prop_error('krylov_local_tensorc : copy (1) '//&
! 'failed.', errst=errst)) return
beta(0) = sqrt(norm(Rs))
! First iteration
! ---------------
jj = 1
call pointto(Krylov(jj), Rs)!, scalar=zone / beta(jj - 1), errst=errst)
!!if(prop_error('krylov_local_tensorc : copy (2) '//&
!! 'failed.', errst=errst)) return
call scale(zone / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo2_dot_psi(Rs, Ham, Krylov(jj), errst=errst)
!if(prop_error('krylov_local_tensorc : mpo2_dot_psi (1) '//&
! 'failed.', errst=errst)) return
alpha(jj) = real(dot(Krylov(jj), Rs, errst=errst), KIND=rKind)
!if(prop_error('krylov_local_tensorc : dot (1) '//&
! 'failed.', errst=errst)) then
! return
!end if
call gaxpy(Rs, -alpha(jj), Krylov(jj))
!if(prop_error('krylov_local_tensorc : gaxpy (1) '//&
! 'failed.', errst=errst)) return
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs, errst=errst)
!if(prop_error('krylov_local_tensorc : dot (2) '//&
! 'failed.', errst=errst)) return
call gaxpy(Rs, -overlap, Krylov(ii), errst=errst)
!if(prop_error('krylov_local_tensorc : gaxpy (2) '//&
! 'failed.', errst=errst)) return
end do
beta(jj) = sqrt(norm(Rs))
if(verbose >= 2) write(slog, *) 'beta', jj, beta(jj)
if(beta(1) < Cp%lanczos_tol) then
call scale(exp(deltat * alpha(1)), Psi)
call destroy(Rs)
call destroy(Krylov(1))
deallocate(Krylov, alpha, beta)
return
end if
estimate = beta(1)
! Further iterations to build Krylov basis
! ----------------------------------------
do jj = 2, nkrylov
if(verbose >= 2) write(slog, *) 'jj', jj, nkrylov, estimate
call pointto(Krylov(jj), Rs)!, scalar=zone / beta(jj - 1), errst=errst)
!if(prop_error('krylov_local_tensorc : copy (4) '//&
! 'failed.', errst=errst)) return
call scale(zone / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo2_dot_psi(Rs, Ham, Krylov(jj), errst=errst)
!if(prop_error('krylov_local_tensorc : mpo2_dot_psi (2) '//&
! 'failed.', errst=errst)) return
alpha(jj) = real(dot(Krylov(jj), Rs, errst=errst), KIND=rKind)
!if(prop_error('krylov_local_tensorc : dot (3) '//&
! 'failed.', errst=errst)) return
call gaxpy(Rs, -alpha(jj), Krylov(jj), errst=errst)
!if(prop_error('krylov_local_tensorc : gaxpy (3) '//&
! 'failed.', errst=errst)) return
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs, errst=errst)
!if(prop_error('krylov_local_tensorc : dot (4) '//&
! 'failed.', errst=errst)) return
call gaxpy(Rs, -overlap, Krylov(ii), errst=errst)
!if(prop_error('krylov_local_tensorc : gaxpy (4) '//&
! 'failed.', errst=errst)) return
end do
beta(jj) = sqrt(norm(Rs))
! Exponentiate tridiagonal matrix
allocate(prop(jj, jj))
call expmh(prop, deltat, alpha(1:jj), beta(1:jj-1), jj, errst=errst)
!if(prop_error('krylov_local_tensorc_real: '//&
! 'expmh failed.', 'KrylovOps_include.f90:1267', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * beta(jj - 1) * abs(prop(jj,1)))
if(verbose >= 2) write(slog, *) 'beta', jj, beta(jj), estimate, &
abs(prop(jj, 1))
deallocate(prop)
if(abs(estimate) <= Cp%lanczos_tol) exit
if(abs(beta(jj)) <= 10.0_rKind * numzero) exit
end do
top = min(jj, nkrylov)
if(top == nkrylov) write(slog, *) 'two-site Krylov failed to '//&
'converge!', beta(ubound(beta, 1)), estimate
allocate(prop(1:top, 1:top))
call expmh(prop, deltat, alpha(1:top), beta(1:top - 1), top, &
errst=errst)
!if(prop_error('krylov_local_tensorc_real: '//&
! 'expmh failed.', 'KrylovOps_include.f90:1290', &
! errst=errst)) then
! return
!end if
call scale(0.0_rKind, Psi)
do ii = 1, top
call gaxpy(Psi, prop(ii,1), Krylov(ii), errst=errst)
!if(prop_error('krylov_local_tensorc : gaxpy (5) '//&
! 'failed.', errst=errst)) return
call destroy(Krylov(ii))
end do
call destroy(Rs)
deallocate(Krylov, alpha, beta, prop)
end subroutine krylov_local_tensorc_real
"""
return
[docs]def krylov_local_qtensor_real():
"""
fortran-subroutine - August 2017 (dj, updated)
Propagate a local two-site Hamiltonian with the Krylov algorithm. This
version does not consider the left-right overlap.
**Arguments**
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
This version uses the two-site Hamiltonian in its tensor
representation. Scaling is chi^2 d^4 for each Krylov step, and an
additional d^4 kappa to build the two site Hamiltonian for this
subroutine. Using the MPO representation, might speed things up
to 2 chi^2 d^2 kappa. (kappa bond dimension of the MPO.)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_local_qtensor_real(Psi, Ham, deltat, Cp, &
pbc, errst)
type(qtensor), intent(inout) :: Psi
type(qtensor), intent(inout) :: Ham
real(KIND=rKind), intent(in) :: deltat
type(ConvParam), intent(in) :: Cp
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj!, kk, pp
! diagonal and 1st offdiagonal of tridiagonal upper hessenberg
real(KIND=rKind), allocatable :: beta(:), alpha(:)
! number of Krylov vectors and Krylov basis
integer :: nkrylov
type(qtensor), pointer :: Krylov(:)
! number of Krylov vectors used
integer :: top
! MPS at tensor w/ orthogonality center to be modified in
! iteration process
type(qtensor) :: Rs!, Tmp
! The exponential in the Krylov basis, including the previous step
real(KIND=rKind), allocatable :: prop(:, :)
REAL(KIND=rKind) :: estimate!, tol
real(KIND=rKind) :: overlap
!INTEGER :: k, n,p,q, chunk, i,j
!LOGICAL :: exitSwitch
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_local_qtensor_'//&
! 'real: Krylov does not have PBC.', 99, &
! 'KrylovOps_include.f90:1158', errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylov(nkrylov), alpha(nkrylov), beta(0:nkrylov))
call copy(Rs, Psi, errst=errst)
!if(prop_error('krylov_local_qtensor : copy (1) '//&
! 'failed.', errst=errst)) return
beta(0) = sqrt(norm(Rs))
! First iteration
! ---------------
jj = 1
call pointto(Krylov(jj), Rs)!, scalar=done / beta(jj - 1), errst=errst)
!!if(prop_error('krylov_local_qtensor : copy (2) '//&
!! 'failed.', errst=errst)) return
call scale(done / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo2_dot_psi(Rs, Ham, Krylov(jj), errst=errst)
!if(prop_error('krylov_local_qtensor : mpo2_dot_psi (1) '//&
! 'failed.', errst=errst)) return
alpha(jj) = real(dot(Krylov(jj), Rs, errst=errst), KIND=rKind)
!if(prop_error('krylov_local_qtensor : dot (1) '//&
! 'failed.', errst=errst)) then
! return
!end if
call gaxpy(Rs, -alpha(jj), Krylov(jj))
!if(prop_error('krylov_local_qtensor : gaxpy (1) '//&
! 'failed.', errst=errst)) return
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs, errst=errst)
!if(prop_error('krylov_local_qtensor : dot (2) '//&
! 'failed.', errst=errst)) return
call gaxpy(Rs, -overlap, Krylov(ii), errst=errst)
!if(prop_error('krylov_local_qtensor : gaxpy (2) '//&
! 'failed.', errst=errst)) return
end do
beta(jj) = sqrt(norm(Rs))
if(verbose >= 2) write(slog, *) 'beta', jj, beta(jj)
if(beta(1) < Cp%lanczos_tol) then
call scale(exp(deltat * alpha(1)), Psi)
call destroy(Rs)
call destroy(Krylov(1))
deallocate(Krylov, alpha, beta)
return
end if
estimate = beta(1)
! Further iterations to build Krylov basis
! ----------------------------------------
do jj = 2, nkrylov
if(verbose >= 2) write(slog, *) 'jj', jj, nkrylov, estimate
call pointto(Krylov(jj), Rs)!, scalar=done / beta(jj - 1), errst=errst)
!if(prop_error('krylov_local_qtensor : copy (4) '//&
! 'failed.', errst=errst)) return
call scale(done / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo2_dot_psi(Rs, Ham, Krylov(jj), errst=errst)
!if(prop_error('krylov_local_qtensor : mpo2_dot_psi (2) '//&
! 'failed.', errst=errst)) return
alpha(jj) = real(dot(Krylov(jj), Rs, errst=errst), KIND=rKind)
!if(prop_error('krylov_local_qtensor : dot (3) '//&
! 'failed.', errst=errst)) return
call gaxpy(Rs, -alpha(jj), Krylov(jj), errst=errst)
!if(prop_error('krylov_local_qtensor : gaxpy (3) '//&
! 'failed.', errst=errst)) return
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs, errst=errst)
!if(prop_error('krylov_local_qtensor : dot (4) '//&
! 'failed.', errst=errst)) return
call gaxpy(Rs, -overlap, Krylov(ii), errst=errst)
!if(prop_error('krylov_local_qtensor : gaxpy (4) '//&
! 'failed.', errst=errst)) return
end do
beta(jj) = sqrt(norm(Rs))
! Exponentiate tridiagonal matrix
allocate(prop(jj, jj))
call expmh(prop, deltat, alpha(1:jj), beta(1:jj-1), jj, errst=errst)
!if(prop_error('krylov_local_qtensor_real: '//&
! 'expmh failed.', 'KrylovOps_include.f90:1267', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * beta(jj - 1) * abs(prop(jj,1)))
if(verbose >= 2) write(slog, *) 'beta', jj, beta(jj), estimate, &
abs(prop(jj, 1))
deallocate(prop)
if(abs(estimate) <= Cp%lanczos_tol) exit
if(abs(beta(jj)) <= 10.0_rKind * numzero) exit
end do
top = min(jj, nkrylov)
if(top == nkrylov) write(slog, *) 'two-site Krylov failed to '//&
'converge!', beta(ubound(beta, 1)), estimate
allocate(prop(1:top, 1:top))
call expmh(prop, deltat, alpha(1:top), beta(1:top - 1), top, &
errst=errst)
!if(prop_error('krylov_local_qtensor_real: '//&
! 'expmh failed.', 'KrylovOps_include.f90:1290', &
! errst=errst)) then
! return
!end if
call scale(0.0_rKind, Psi)
do ii = 1, top
call gaxpy(Psi, prop(ii,1), Krylov(ii), errst=errst)
!if(prop_error('krylov_local_qtensor : gaxpy (5) '//&
! 'failed.', errst=errst)) return
call destroy(Krylov(ii))
end do
call destroy(Rs)
deallocate(Krylov, alpha, beta, prop)
end subroutine krylov_local_qtensor_real
"""
return
[docs]def krylov_local_qtensorc_real():
"""
fortran-subroutine - August 2017 (dj, updated)
Propagate a local two-site Hamiltonian with the Krylov algorithm. This
version does not consider the left-right overlap.
**Arguments**
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
This version uses the two-site Hamiltonian in its tensor
representation. Scaling is chi^2 d^4 for each Krylov step, and an
additional d^4 kappa to build the two site Hamiltonian for this
subroutine. Using the MPO representation, might speed things up
to 2 chi^2 d^2 kappa. (kappa bond dimension of the MPO.)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_local_qtensorc_real(Psi, Ham, deltat, Cp, &
pbc, errst)
type(qtensorc), intent(inout) :: Psi
type(qtensorc), intent(inout) :: Ham
real(KIND=rKind), intent(in) :: deltat
type(ConvParam), intent(in) :: Cp
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj!, kk, pp
! diagonal and 1st offdiagonal of tridiagonal upper hessenberg
real(KIND=rKind), allocatable :: beta(:), alpha(:)
! number of Krylov vectors and Krylov basis
integer :: nkrylov
type(qtensorc), pointer :: Krylov(:)
! number of Krylov vectors used
integer :: top
! MPS at tensor w/ orthogonality center to be modified in
! iteration process
type(qtensorc) :: Rs!, Tmp
! The exponential in the Krylov basis, including the previous step
real(KIND=rKind), allocatable :: prop(:, :)
REAL(KIND=rKind) :: estimate!, tol
complex(KIND=rKind) :: overlap
!INTEGER :: k, n,p,q, chunk, i,j
!LOGICAL :: exitSwitch
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_local_qtensorc_'//&
! 'real: Krylov does not have PBC.', 99, &
! 'KrylovOps_include.f90:1158', errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
! Allocate space
allocate(Krylov(nkrylov), alpha(nkrylov), beta(0:nkrylov))
call copy(Rs, Psi, errst=errst)
!if(prop_error('krylov_local_qtensorc : copy (1) '//&
! 'failed.', errst=errst)) return
beta(0) = sqrt(norm(Rs))
! First iteration
! ---------------
jj = 1
call pointto(Krylov(jj), Rs)!, scalar=zone / beta(jj - 1), errst=errst)
!!if(prop_error('krylov_local_qtensorc : copy (2) '//&
!! 'failed.', errst=errst)) return
call scale(zone / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo2_dot_psi(Rs, Ham, Krylov(jj), errst=errst)
!if(prop_error('krylov_local_qtensorc : mpo2_dot_psi (1) '//&
! 'failed.', errst=errst)) return
alpha(jj) = real(dot(Krylov(jj), Rs, errst=errst), KIND=rKind)
!if(prop_error('krylov_local_qtensorc : dot (1) '//&
! 'failed.', errst=errst)) then
! return
!end if
call gaxpy(Rs, -alpha(jj), Krylov(jj))
!if(prop_error('krylov_local_qtensorc : gaxpy (1) '//&
! 'failed.', errst=errst)) return
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs, errst=errst)
!if(prop_error('krylov_local_qtensorc : dot (2) '//&
! 'failed.', errst=errst)) return
call gaxpy(Rs, -overlap, Krylov(ii), errst=errst)
!if(prop_error('krylov_local_qtensorc : gaxpy (2) '//&
! 'failed.', errst=errst)) return
end do
beta(jj) = sqrt(norm(Rs))
if(verbose >= 2) write(slog, *) 'beta', jj, beta(jj)
if(beta(1) < Cp%lanczos_tol) then
call scale(exp(deltat * alpha(1)), Psi)
call destroy(Rs)
call destroy(Krylov(1))
deallocate(Krylov, alpha, beta)
return
end if
estimate = beta(1)
! Further iterations to build Krylov basis
! ----------------------------------------
do jj = 2, nkrylov
if(verbose >= 2) write(slog, *) 'jj', jj, nkrylov, estimate
call pointto(Krylov(jj), Rs)!, scalar=zone / beta(jj - 1), errst=errst)
!if(prop_error('krylov_local_qtensorc : copy (4) '//&
! 'failed.', errst=errst)) return
call scale(zone / beta(jj - 1), Krylov(jj))
!call destroy(Rs)
call mpo2_dot_psi(Rs, Ham, Krylov(jj), errst=errst)
!if(prop_error('krylov_local_qtensorc : mpo2_dot_psi (2) '//&
! 'failed.', errst=errst)) return
alpha(jj) = real(dot(Krylov(jj), Rs, errst=errst), KIND=rKind)
!if(prop_error('krylov_local_qtensorc : dot (3) '//&
! 'failed.', errst=errst)) return
call gaxpy(Rs, -alpha(jj), Krylov(jj), errst=errst)
!if(prop_error('krylov_local_qtensorc : gaxpy (3) '//&
! 'failed.', errst=errst)) return
! Reorthogonalize
do ii = 1, jj
overlap = dot(Krylov(ii), Rs, errst=errst)
!if(prop_error('krylov_local_qtensorc : dot (4) '//&
! 'failed.', errst=errst)) return
call gaxpy(Rs, -overlap, Krylov(ii), errst=errst)
!if(prop_error('krylov_local_qtensorc : gaxpy (4) '//&
! 'failed.', errst=errst)) return
end do
beta(jj) = sqrt(norm(Rs))
! Exponentiate tridiagonal matrix
allocate(prop(jj, jj))
call expmh(prop, deltat, alpha(1:jj), beta(1:jj-1), jj, errst=errst)
!if(prop_error('krylov_local_qtensorc_real: '//&
! 'expmh failed.', 'KrylovOps_include.f90:1267', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * beta(jj - 1) * abs(prop(jj,1)))
if(verbose >= 2) write(slog, *) 'beta', jj, beta(jj), estimate, &
abs(prop(jj, 1))
deallocate(prop)
if(abs(estimate) <= Cp%lanczos_tol) exit
if(abs(beta(jj)) <= 10.0_rKind * numzero) exit
end do
top = min(jj, nkrylov)
if(top == nkrylov) write(slog, *) 'two-site Krylov failed to '//&
'converge!', beta(ubound(beta, 1)), estimate
allocate(prop(1:top, 1:top))
call expmh(prop, deltat, alpha(1:top), beta(1:top - 1), top, &
errst=errst)
!if(prop_error('krylov_local_qtensorc_real: '//&
! 'expmh failed.', 'KrylovOps_include.f90:1290', &
! errst=errst)) then
! return
!end if
call scale(0.0_rKind, Psi)
do ii = 1, top
call gaxpy(Psi, prop(ii,1), Krylov(ii), errst=errst)
!if(prop_error('krylov_local_qtensorc : gaxpy (5) '//&
! 'failed.', errst=errst)) return
call destroy(Krylov(ii))
end do
call destroy(Rs)
deallocate(Krylov, alpha, beta, prop)
end subroutine krylov_local_qtensorc_real
"""
return
[docs]def krylov_arnoldi_local_tensor_complex():
"""
fortran-subroutine - August 2017 (dj, updated)
Propagate an MPDO with the two-site Krylov-Arnoldi algorithm.
**Arguments**
Psi : TYPE(tensorc), inout
Represents the wave function of the tensor network. Has
to contain the orthogonality center.
Ham : TYPE(tensor), inout
The operator of which the exponential should be applied to
the state Psi. Must be two-site Hamiltonian.
deltat : COMPLEX, in
Time step for the time evolution.
Cp : TYPE(ConvParam), in
Contains the convergence parameters for the algorithm. Subroutine
accesses lanczos_tol and max_num_lanczos_iter.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
The Krylov-Arnoldi algorithm allows to propagate with
non-hermitian two-site operators.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_arnoldi_local_tensor_complex(Psi, Ham, deltat, &
Cp, pbc, errst)
type(tensorc), intent(inout) :: Psi
type(tensor), intent(inout) :: Ham
complex(KIND=rKind), intent(in) :: deltat
type(ConvParam), intent(in) :: Cp
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk
! number of krylov vectors / Krylov basis
integer :: nkrylov
type(tensorc), pointer :: Krylov(:)
! Upper Hessenberg matrix / upper Hessenberg weighted with dt
complex(KIND=rKind), dimension(:, :), allocatable :: hessen, hessendt
! intial norm
real(KIND=rKind) :: beta0
real(KIND=rKind) :: estimate
! vector for Krylov algorithm - modification carried out on r
type(tensorc) :: Rs, Tmp
! propagator
complex(KIND=rKind), allocatable :: prop(:, :)
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_arnoldi_local'//&
! '_tensor_complex: Krylov does not have PBC.', &
! 99, 'KrylovOps_include.f90:1405', errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
allocate(Krylov(nkrylov), hessen(nkrylov + 1, nkrylov + 1))
hessen = 0.0_rKind
call copy(Rs, Psi)
beta0 = sqrt(norm(Rs))
! First data point - initialize values
! ------------------------------------
jj = 1
call pointto(Krylov(jj), Rs)!, scalar=zone / beta0)
call scale(zone / beta0, Krylov(jj))
!call destroy(Rs)
call copy(Tmp, Krylov(jj))
call mpo2_dot_psi(Rs, Ham, Tmp)
call destroy(Tmp)
hessen(jj + 1, jj) = sqrt(norm(Rs))
if(verbose >= 2) write(slog, *) 'beta', jj, hessen(jj + 1, jj)
if(real(hessen(jj + 1, jj), KIND=rKind) < Cp%lanczos_tol) then
call scale(exp(deltat * hessen(1 , 1)), Psi)
call destroy(Rs)
call destroy(Krylov(1))
deallocate(Krylov, hessen)
return
end if
estimate = abs(hessen(2, 1))
! Iteration over vectors 2, ...
! -----------------------------
do jj = 2, nkrylov
if(verbose >= 2) write(slog, *) 'jj', jj, nkrylov, estimate
call copy(Krylov(jj), Rs, scalar=zone / hessen(jj, jj - 1))
call destroy(Rs)
call copy(Tmp, Krylov(jj))
call mpo2_dot_psi(Rs, Ham, Tmp)
call destroy(Tmp)
do kk = 1, jj - 1
hessen(kk, jj) = dot(Rs, Krylov(kk))
call gaxpy(Rs, - hessen(kk, jj), Krylov(kk))
end do
hessen(jj + 1, jj) = sqrt(norm(Rs))
! Exponentiate matrix
allocate(hessendt(jj, jj), prop(jj, jj))
hessendt = hessen
call expm(prop(1:jj, 1:jj), deltat, hessendt(1:jj, 1:jj), jj, &
errst=errst)
!if(prop_error('krylov_arnoldi_local_tensor_complex'//&
! ': expm failed.', 'KrylovOps_include.f90:1474', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * hessen(jj, jj - 1) * abs(prop(jj, 1)))
if(verbose >= 2) write(slog, *) 'beta', jj, hessen(jj + 1, jj), &
estimate, abs(prop(jj, 1))
if(abs(estimate) <= Cp%lanczos_tol) exit
if(abs(hessen(jj + 1, jj)) <= 10.0_rKind * numzero) exit
! Avoid increment to jj = nkrylov + 1
if(jj == nkrylov) exit
deallocate(hessendt, prop)
end do
! Get the final propagated state
! ------------------------------
if(jj == nKrylov) write(slog, *) 'two-site Krylov failed to '//&
'converge!', hessen(jj + 1, 1), estimate
prop(1:jj, 1) = prop(1:jj, 1) * beta0
call scale(0.0_rKind, Psi)
do ii = 1, jj
call gaxpy(Psi, prop(ii, 1), Krylov(ii))
call destroy(Krylov(ii))
end do
call destroy(Rs)
deallocate(Krylov, hessen, hessendt, prop)
end subroutine krylov_arnoldi_local_tensor_complex
"""
return
[docs]def krylov_arnoldi_local_tensorc_complex():
"""
fortran-subroutine - August 2017 (dj, updated)
Propagate an MPDO with the two-site Krylov-Arnoldi algorithm.
**Arguments**
Psi : TYPE(tensorc), inout
Represents the wave function of the tensor network. Has
to contain the orthogonality center.
Ham : TYPE(tensorc), inout
The operator of which the exponential should be applied to
the state Psi. Must be two-site Hamiltonian.
deltat : COMPLEX, in
Time step for the time evolution.
Cp : TYPE(ConvParam), in
Contains the convergence parameters for the algorithm. Subroutine
accesses lanczos_tol and max_num_lanczos_iter.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
The Krylov-Arnoldi algorithm allows to propagate with
non-hermitian two-site operators.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_arnoldi_local_tensorc_complex(Psi, Ham, deltat, &
Cp, pbc, errst)
type(tensorc), intent(inout) :: Psi
type(tensorc), intent(inout) :: Ham
complex(KIND=rKind), intent(in) :: deltat
type(ConvParam), intent(in) :: Cp
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk
! number of krylov vectors / Krylov basis
integer :: nkrylov
type(tensorc), pointer :: Krylov(:)
! Upper Hessenberg matrix / upper Hessenberg weighted with dt
complex(KIND=rKind), dimension(:, :), allocatable :: hessen, hessendt
! intial norm
real(KIND=rKind) :: beta0
real(KIND=rKind) :: estimate
! vector for Krylov algorithm - modification carried out on r
type(tensorc) :: Rs, Tmp
! propagator
complex(KIND=rKind), allocatable :: prop(:, :)
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_arnoldi_local'//&
! '_tensorc_complex: Krylov does not have PBC.', &
! 99, 'KrylovOps_include.f90:1405', errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
allocate(Krylov(nkrylov), hessen(nkrylov + 1, nkrylov + 1))
hessen = 0.0_rKind
call copy(Rs, Psi)
beta0 = sqrt(norm(Rs))
! First data point - initialize values
! ------------------------------------
jj = 1
call pointto(Krylov(jj), Rs)!, scalar=zone / beta0)
call scale(zone / beta0, Krylov(jj))
!call destroy(Rs)
call copy(Tmp, Krylov(jj))
call mpo2_dot_psi(Rs, Ham, Tmp)
call destroy(Tmp)
hessen(jj + 1, jj) = sqrt(norm(Rs))
if(verbose >= 2) write(slog, *) 'beta', jj, hessen(jj + 1, jj)
if(real(hessen(jj + 1, jj), KIND=rKind) < Cp%lanczos_tol) then
call scale(exp(deltat * hessen(1 , 1)), Psi)
call destroy(Rs)
call destroy(Krylov(1))
deallocate(Krylov, hessen)
return
end if
estimate = abs(hessen(2, 1))
! Iteration over vectors 2, ...
! -----------------------------
do jj = 2, nkrylov
if(verbose >= 2) write(slog, *) 'jj', jj, nkrylov, estimate
call copy(Krylov(jj), Rs, scalar=zone / hessen(jj, jj - 1))
call destroy(Rs)
call copy(Tmp, Krylov(jj))
call mpo2_dot_psi(Rs, Ham, Tmp)
call destroy(Tmp)
do kk = 1, jj - 1
hessen(kk, jj) = dot(Rs, Krylov(kk))
call gaxpy(Rs, - hessen(kk, jj), Krylov(kk))
end do
hessen(jj + 1, jj) = sqrt(norm(Rs))
! Exponentiate matrix
allocate(hessendt(jj, jj), prop(jj, jj))
hessendt = hessen
call expm(prop(1:jj, 1:jj), deltat, hessendt(1:jj, 1:jj), jj, &
errst=errst)
!if(prop_error('krylov_arnoldi_local_tensorc_complex'//&
! ': expm failed.', 'KrylovOps_include.f90:1474', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * hessen(jj, jj - 1) * abs(prop(jj, 1)))
if(verbose >= 2) write(slog, *) 'beta', jj, hessen(jj + 1, jj), &
estimate, abs(prop(jj, 1))
if(abs(estimate) <= Cp%lanczos_tol) exit
if(abs(hessen(jj + 1, jj)) <= 10.0_rKind * numzero) exit
! Avoid increment to jj = nkrylov + 1
if(jj == nkrylov) exit
deallocate(hessendt, prop)
end do
! Get the final propagated state
! ------------------------------
if(jj == nKrylov) write(slog, *) 'two-site Krylov failed to '//&
'converge!', hessen(jj + 1, 1), estimate
prop(1:jj, 1) = prop(1:jj, 1) * beta0
call scale(0.0_rKind, Psi)
do ii = 1, jj
call gaxpy(Psi, prop(ii, 1), Krylov(ii))
call destroy(Krylov(ii))
end do
call destroy(Rs)
deallocate(Krylov, hessen, hessendt, prop)
end subroutine krylov_arnoldi_local_tensorc_complex
"""
return
[docs]def krylov_arnoldi_local_qtensor_complex():
"""
fortran-subroutine - August 2017 (dj, updated)
Propagate an MPDO with the two-site Krylov-Arnoldi algorithm.
**Arguments**
Psi : TYPE(qtensorc), inout
Represents the wave function of the tensor network. Has
to contain the orthogonality center.
Ham : TYPE(qtensor), inout
The operator of which the exponential should be applied to
the state Psi. Must be two-site Hamiltonian.
deltat : COMPLEX, in
Time step for the time evolution.
Cp : TYPE(ConvParam), in
Contains the convergence parameters for the algorithm. Subroutine
accesses lanczos_tol and max_num_lanczos_iter.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
The Krylov-Arnoldi algorithm allows to propagate with
non-hermitian two-site operators.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_arnoldi_local_qtensor_complex(Psi, Ham, deltat, &
Cp, pbc, errst)
type(qtensorc), intent(inout) :: Psi
type(qtensor), intent(inout) :: Ham
complex(KIND=rKind), intent(in) :: deltat
type(ConvParam), intent(in) :: Cp
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk
! number of krylov vectors / Krylov basis
integer :: nkrylov
type(qtensorc), pointer :: Krylov(:)
! Upper Hessenberg matrix / upper Hessenberg weighted with dt
complex(KIND=rKind), dimension(:, :), allocatable :: hessen, hessendt
! intial norm
real(KIND=rKind) :: beta0
real(KIND=rKind) :: estimate
! vector for Krylov algorithm - modification carried out on r
type(qtensorc) :: Rs, Tmp
! propagator
complex(KIND=rKind), allocatable :: prop(:, :)
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_arnoldi_local'//&
! '_qtensor_complex: Krylov does not have PBC.', &
! 99, 'KrylovOps_include.f90:1405', errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
allocate(Krylov(nkrylov), hessen(nkrylov + 1, nkrylov + 1))
hessen = 0.0_rKind
call copy(Rs, Psi)
beta0 = sqrt(norm(Rs))
! First data point - initialize values
! ------------------------------------
jj = 1
call pointto(Krylov(jj), Rs)!, scalar=zone / beta0)
call scale(zone / beta0, Krylov(jj))
!call destroy(Rs)
call copy(Tmp, Krylov(jj))
call mpo2_dot_psi(Rs, Ham, Tmp)
call destroy(Tmp)
hessen(jj + 1, jj) = sqrt(norm(Rs))
if(verbose >= 2) write(slog, *) 'beta', jj, hessen(jj + 1, jj)
if(real(hessen(jj + 1, jj), KIND=rKind) < Cp%lanczos_tol) then
call scale(exp(deltat * hessen(1 , 1)), Psi)
call destroy(Rs)
call destroy(Krylov(1))
deallocate(Krylov, hessen)
return
end if
estimate = abs(hessen(2, 1))
! Iteration over vectors 2, ...
! -----------------------------
do jj = 2, nkrylov
if(verbose >= 2) write(slog, *) 'jj', jj, nkrylov, estimate
call copy(Krylov(jj), Rs, scalar=zone / hessen(jj, jj - 1))
call destroy(Rs)
call copy(Tmp, Krylov(jj))
call mpo2_dot_psi(Rs, Ham, Tmp)
call destroy(Tmp)
do kk = 1, jj - 1
hessen(kk, jj) = dot(Rs, Krylov(kk))
call gaxpy(Rs, - hessen(kk, jj), Krylov(kk))
end do
hessen(jj + 1, jj) = sqrt(norm(Rs))
! Exponentiate matrix
allocate(hessendt(jj, jj), prop(jj, jj))
hessendt = hessen
call expm(prop(1:jj, 1:jj), deltat, hessendt(1:jj, 1:jj), jj, &
errst=errst)
!if(prop_error('krylov_arnoldi_local_qtensor_complex'//&
! ': expm failed.', 'KrylovOps_include.f90:1474', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * hessen(jj, jj - 1) * abs(prop(jj, 1)))
if(verbose >= 2) write(slog, *) 'beta', jj, hessen(jj + 1, jj), &
estimate, abs(prop(jj, 1))
if(abs(estimate) <= Cp%lanczos_tol) exit
if(abs(hessen(jj + 1, jj)) <= 10.0_rKind * numzero) exit
! Avoid increment to jj = nkrylov + 1
if(jj == nkrylov) exit
deallocate(hessendt, prop)
end do
! Get the final propagated state
! ------------------------------
if(jj == nKrylov) write(slog, *) 'two-site Krylov failed to '//&
'converge!', hessen(jj + 1, 1), estimate
prop(1:jj, 1) = prop(1:jj, 1) * beta0
call scale(0.0_rKind, Psi)
do ii = 1, jj
call gaxpy(Psi, prop(ii, 1), Krylov(ii))
call destroy(Krylov(ii))
end do
call destroy(Rs)
deallocate(Krylov, hessen, hessendt, prop)
end subroutine krylov_arnoldi_local_qtensor_complex
"""
return
[docs]def krylov_arnoldi_local_qtensorc_complex():
"""
fortran-subroutine - August 2017 (dj, updated)
Propagate an MPDO with the two-site Krylov-Arnoldi algorithm.
**Arguments**
Psi : TYPE(qtensorc), inout
Represents the wave function of the tensor network. Has
to contain the orthogonality center.
Ham : TYPE(qtensorc), inout
The operator of which the exponential should be applied to
the state Psi. Must be two-site Hamiltonian.
deltat : COMPLEX, in
Time step for the time evolution.
Cp : TYPE(ConvParam), in
Contains the convergence parameters for the algorithm. Subroutine
accesses lanczos_tol and max_num_lanczos_iter.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Details**
The Krylov-Arnoldi algorithm allows to propagate with
non-hermitian two-site operators.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine krylov_arnoldi_local_qtensorc_complex(Psi, Ham, deltat, &
Cp, pbc, errst)
type(qtensorc), intent(inout) :: Psi
type(qtensorc), intent(inout) :: Ham
complex(KIND=rKind), intent(in) :: deltat
type(ConvParam), intent(in) :: Cp
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk
! number of krylov vectors / Krylov basis
integer :: nkrylov
type(qtensorc), pointer :: Krylov(:)
! Upper Hessenberg matrix / upper Hessenberg weighted with dt
complex(KIND=rKind), dimension(:, :), allocatable :: hessen, hessendt
! intial norm
real(KIND=rKind) :: beta0
real(KIND=rKind) :: estimate
! vector for Krylov algorithm - modification carried out on r
type(qtensorc) :: Rs, Tmp
! propagator
complex(KIND=rKind), allocatable :: prop(:, :)
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('krylov_arnoldi_local'//&
! '_qtensorc_complex: Krylov does not have PBC.', &
! 99, 'KrylovOps_include.f90:1405', errst=errst)
! return
!end if
! Short-cuts
nkrylov = Cp%max_num_lanczos_iter
allocate(Krylov(nkrylov), hessen(nkrylov + 1, nkrylov + 1))
hessen = 0.0_rKind
call copy(Rs, Psi)
beta0 = sqrt(norm(Rs))
! First data point - initialize values
! ------------------------------------
jj = 1
call pointto(Krylov(jj), Rs)!, scalar=zone / beta0)
call scale(zone / beta0, Krylov(jj))
!call destroy(Rs)
call copy(Tmp, Krylov(jj))
call mpo2_dot_psi(Rs, Ham, Tmp)
call destroy(Tmp)
hessen(jj + 1, jj) = sqrt(norm(Rs))
if(verbose >= 2) write(slog, *) 'beta', jj, hessen(jj + 1, jj)
if(real(hessen(jj + 1, jj), KIND=rKind) < Cp%lanczos_tol) then
call scale(exp(deltat * hessen(1 , 1)), Psi)
call destroy(Rs)
call destroy(Krylov(1))
deallocate(Krylov, hessen)
return
end if
estimate = abs(hessen(2, 1))
! Iteration over vectors 2, ...
! -----------------------------
do jj = 2, nkrylov
if(verbose >= 2) write(slog, *) 'jj', jj, nkrylov, estimate
call copy(Krylov(jj), Rs, scalar=zone / hessen(jj, jj - 1))
call destroy(Rs)
call copy(Tmp, Krylov(jj))
call mpo2_dot_psi(Rs, Ham, Tmp)
call destroy(Tmp)
do kk = 1, jj - 1
hessen(kk, jj) = dot(Rs, Krylov(kk))
call gaxpy(Rs, - hessen(kk, jj), Krylov(kk))
end do
hessen(jj + 1, jj) = sqrt(norm(Rs))
! Exponentiate matrix
allocate(hessendt(jj, jj), prop(jj, jj))
hessendt = hessen
call expm(prop(1:jj, 1:jj), deltat, hessendt(1:jj, 1:jj), jj, &
errst=errst)
!if(prop_error('krylov_arnoldi_local_qtensorc_complex'//&
! ': expm failed.', 'KrylovOps_include.f90:1474', &
! errst=errst)) then
! return
!end if
estimate = real(2.0_rKind * hessen(jj, jj - 1) * abs(prop(jj, 1)))
if(verbose >= 2) write(slog, *) 'beta', jj, hessen(jj + 1, jj), &
estimate, abs(prop(jj, 1))
if(abs(estimate) <= Cp%lanczos_tol) exit
if(abs(hessen(jj + 1, jj)) <= 10.0_rKind * numzero) exit
! Avoid increment to jj = nkrylov + 1
if(jj == nkrylov) exit
deallocate(hessendt, prop)
end do
! Get the final propagated state
! ------------------------------
if(jj == nKrylov) write(slog, *) 'two-site Krylov failed to '//&
'converge!', hessen(jj + 1, 1), estimate
prop(1:jj, 1) = prop(1:jj, 1) * beta0
call scale(0.0_rKind, Psi)
do ii = 1, jj
call gaxpy(Psi, prop(ii, 1), Krylov(ii))
call destroy(Krylov(ii))
end do
call destroy(Rs)
deallocate(Krylov, hessen, hessendt, prop)
end subroutine krylov_arnoldi_local_qtensorc_complex
"""
return