From 40ce92c2ae695ad0aea82a4a47358a8c2e95116d Mon Sep 17 00:00:00 2001 From: Miles Date: Wed, 27 Mar 2024 15:34:23 +0900 Subject: [PATCH 01/26] Remove codes, now moved to ITensors.ITensorMPS --- src/ITensorTDVP.jl | 22 +- src/alternating_update.jl | 157 -------------- src/contract_mpo_mps.jl | 59 ------ src/defaults.jl | 50 ----- src/dmrg.jl | 58 ----- src/dmrg_x.jl | 16 -- src/linsolve.jl | 42 ---- src/projmpo_apply.jl | 90 -------- src/projmpo_mps2.jl | 47 ----- src/projmps2.jl | 123 ----------- src/solver_utils.jl | 69 ------ src/sweep_update.jl | 434 -------------------------------------- src/tdvp.jl | 136 ------------ src/tdvp_sweeps.jl | 17 -- src/tdvpinfo.jl | 7 - src/tdvporder.jl | 24 --- src/update_observer.jl | 6 - 17 files changed, 6 insertions(+), 1351 deletions(-) delete mode 100644 src/alternating_update.jl delete mode 100644 src/contract_mpo_mps.jl delete mode 100644 src/defaults.jl delete mode 100644 src/dmrg.jl delete mode 100644 src/dmrg_x.jl delete mode 100644 src/linsolve.jl delete mode 100644 src/projmpo_apply.jl delete mode 100644 src/projmpo_mps2.jl delete mode 100644 src/projmps2.jl delete mode 100644 src/solver_utils.jl delete mode 100644 src/sweep_update.jl delete mode 100644 src/tdvp.jl delete mode 100644 src/tdvp_sweeps.jl delete mode 100644 src/tdvpinfo.jl delete mode 100644 src/tdvporder.jl delete mode 100644 src/update_observer.jl diff --git a/src/ITensorTDVP.jl b/src/ITensorTDVP.jl index de12a47..4566a84 100644 --- a/src/ITensorTDVP.jl +++ b/src/ITensorTDVP.jl @@ -1,19 +1,9 @@ module ITensorTDVP -export TimeDependentSum, dmrg_x, linsolve, tdvp, to_vec -include("defaults.jl") -include("update_observer.jl") -include("solver_utils.jl") -include("tdvporder.jl") -include("tdvpinfo.jl") -include("sweep_update.jl") -include("alternating_update.jl") -include("tdvp.jl") -include("dmrg.jl") -include("dmrg_x.jl") -include("projmpo_apply.jl") -include("contract_mpo_mps.jl") -include("projmps2.jl") -include("projmpo_mps2.jl") -include("linsolve.jl") +using Reexport: @reexport +@reexport using ITensors.ITensorMPS: tdvp, dmrg_x, to_vec, TimeDependentSum, linsolve + +using ITensors: ITensorMPS +dmrg(args...; kwargs...) = ITensorMPS.alternate_dmrg(args...; kwargs...) + end diff --git a/src/alternating_update.jl b/src/alternating_update.jl deleted file mode 100644 index 4d5846f..0000000 --- a/src/alternating_update.jl +++ /dev/null @@ -1,157 +0,0 @@ -using ITensors: permute -using ITensors.ITensorMPS: - AbstractObserver, - MPO, - MPS, - ProjMPO, - ProjMPOSum, - check_hascommoninds, - checkdone!, - disk, - linkind, - maxlinkdim, - siteinds - -function _compute_nsweeps(t; time_step=default_time_step(t), nsweeps=default_nsweeps()) - if isinf(t) && isnothing(nsweeps) - nsweeps = 1 - elseif !isnothing(nsweeps) && time_step != t - error("Cannot specify both time_step and nsweeps in alternating_update") - elseif isfinite(time_step) && abs(time_step) > 0 && isnothing(nsweeps) - nsweeps = convert(Int, ceil(abs(t / time_step))) - if !(nsweeps * time_step ≈ t) - error("Time step $time_step not commensurate with total time t=$t") - end - end - return nsweeps -end - -function _extend_sweeps_param(param, nsweeps) - if param isa Number - eparam = fill(param, nsweeps) - else - length(param) == nsweeps && return param - eparam = Vector(undef, nsweeps) - eparam[1:length(param)] = param - eparam[(length(param) + 1):end] .= param[end] - end - return eparam -end - -function process_sweeps(; nsweeps, maxdim, mindim, cutoff, noise) - maxdim = _extend_sweeps_param(maxdim, nsweeps) - mindim = _extend_sweeps_param(mindim, nsweeps) - cutoff = _extend_sweeps_param(cutoff, nsweeps) - noise = _extend_sweeps_param(noise, nsweeps) - return (; maxdim, mindim, cutoff, noise) -end - -function alternating_update( - solver, - PH, - t::Number, - psi0::MPS; - nsweeps=default_nsweeps(), - checkdone=default_checkdone(), - write_when_maxdim_exceeds=default_write_when_maxdim_exceeds(), - nsite=default_nsite(), - reverse_step=default_reverse_step(), - time_start=default_time_start(), - time_step=default_time_step(t), - order=default_order(), - (observer!)=default_observer!(), - (step_observer!)=default_step_observer!(), - outputlevel=default_outputlevel(), - normalize=default_normalize(), - maxdim=default_maxdim(), - mindim=default_mindim(), - cutoff=default_cutoff(Float64), - noise=default_noise(), -) - nsweeps = _compute_nsweeps(t; time_step, nsweeps) - maxdim, mindim, cutoff, noise = process_sweeps(; nsweeps, maxdim, mindim, cutoff, noise) - forward_order = TDVPOrder(order, Base.Forward) - psi = copy(psi0) - # Keep track of the start of the current time step. - # Helpful for tracking the total time, for example - # when using time-dependent solvers. - # This will be passed as a keyword argument to the - # `solver`. - current_time = time_start - info = nothing - for sweep in 1:nsweeps - if !isnothing(write_when_maxdim_exceeds) && maxdim[sweep] > write_when_maxdim_exceeds - if outputlevel >= 2 - println( - "write_when_maxdim_exceeds = $write_when_maxdim_exceeds and maxdim(sweeps, sw) = $(maxdim(sweeps, sweep)), writing environment tensors to disk", - ) - end - PH = disk(PH) - end - sweep_time = @elapsed begin - psi, PH, info = sweep_update( - forward_order, - solver, - PH, - time_step, - psi; - nsite, - current_time, - reverse_step, - sweep, - observer!, - normalize, - maxdim=maxdim[sweep], - mindim=mindim[sweep], - cutoff=cutoff[sweep], - noise=noise[sweep], - ) - end - current_time += time_step - update!(step_observer!; psi, sweep, outputlevel, current_time) - if outputlevel >= 1 - print("After sweep ", sweep, ":") - print(" maxlinkdim=", maxlinkdim(psi)) - @printf(" maxerr=%.2E", info.maxtruncerr) - print(" current_time=", round(current_time; digits=3)) - print(" time=", round(sweep_time; digits=3)) - println() - flush(stdout) - end - isdone = false - if !isnothing(checkdone) - isdone = checkdone(; psi, sweep, outputlevel) - elseif observer! isa AbstractObserver - isdone = checkdone!(observer!; psi, sweep, outputlevel) - end - isdone && break - end - return psi -end - -# Convenience wrapper to not have to specify time step. -# Use a time step of `Inf` as a convention, since TDVP -# with an infinite time step corresponds to DMRG. -function alternating_update(solver, H, psi0::MPS; kwargs...) - return alternating_update(solver, H, ITensors.scalartype(psi0)(Inf), psi0; kwargs...) -end - -function alternating_update(solver, H::MPO, t::Number, psi0::MPS; kwargs...) - check_hascommoninds(siteinds, H, psi0) - check_hascommoninds(siteinds, H, psi0') - # Permute the indices to have a better memory layout - # and minimize permutations - H = permute(H, (linkind, siteinds, linkind)) - PH = ProjMPO(H) - return alternating_update(solver, PH, t, psi0; kwargs...) -end - -function alternating_update(solver, Hs::Vector{MPO}, t::Number, psi0::MPS; kwargs...) - for H in Hs - check_hascommoninds(siteinds, H, psi0) - check_hascommoninds(siteinds, H, psi0') - end - Hs .= ITensors.permute.(Hs, Ref((linkind, siteinds, linkind))) - PHs = ProjMPOSum(Hs) - return alternating_update(solver, PHs, t, psi0; kwargs...) -end diff --git a/src/contract_mpo_mps.jl b/src/contract_mpo_mps.jl deleted file mode 100644 index 8f2dac4..0000000 --- a/src/contract_mpo_mps.jl +++ /dev/null @@ -1,59 +0,0 @@ -using ITensors: - ITensors, - Index, - ITensor, - @Algorithm_str, - commoninds, - contract, - hasind, - linkinds, - replace_siteinds!, - sim, - siteinds - -function contractmpo_solver(; kwargs...) - function solver(PH, t, psi; kws...) - v = ITensor(true) - for j in (PH.lpos + 1):(PH.rpos - 1) - v *= PH.psi0[j] - end - Hpsi0 = contract(PH, v) - return Hpsi0, nothing - end - return solver -end - -function ITensors.contract( - ::Algorithm"fit", A::MPO, psi0::MPS; init_mps=psi0, nsweeps=1, kwargs... -)::MPS - n = length(A) - n != length(psi0) && - throw(DimensionMismatch("lengths of MPO ($n) and MPS ($(length(psi0))) do not match")) - if n == 1 - return MPS([A[1] * psi0[1]]) - end - any(i -> isempty(i), siteinds(commoninds, A, psi0)) && - error("In `contract(A::MPO, x::MPS)`, `A` and `x` must share a set of site indices") - # In case A and psi0 have the same link indices - A = sim(linkinds, A) - # Fix site and link inds of init_mps - init_mps = deepcopy(init_mps) - init_mps = sim(linkinds, init_mps) - Ai = siteinds(A) - ti = Vector{Index}(undef, n) - for j in 1:n - for i in Ai[j] - if !hasind(psi0[j], i) - ti[j] = i - break - end - end - end - replace_siteinds!(init_mps, ti) - reverse_step = false - PH = ProjMPOApply(psi0, A) - psi = alternating_update( - contractmpo_solver(; kwargs...), PH, init_mps; nsweeps, reverse_step, kwargs... - ) - return psi -end diff --git a/src/defaults.jl b/src/defaults.jl deleted file mode 100644 index 4a0f13b..0000000 --- a/src/defaults.jl +++ /dev/null @@ -1,50 +0,0 @@ -using ITensors: NoObserver -using KrylovKit: eigsolve, exponentiate - -default_nsweeps() = nothing -default_checkdone() = nothing -default_write_when_maxdim_exceeds() = nothing -default_nsite() = 2 -default_reverse_step() = true -default_time_start() = 0 -default_time_step(t) = t -default_order() = 2 -default_observer!() = NoObserver() -default_step_observer!() = NoObserver() -default_outputlevel() = 0 -default_normalize() = false -default_sweep() = 1 -default_current_time() = 0 - -# Truncation -default_maxdim() = typemax(Int) -default_mindim() = 1 -default_cutoff(type::Type{<:Number}) = eps(real(type)) -default_noise() = 0 - -# Solvers -default_tdvp_solver_backend() = "exponentiate" -default_ishermitian() = true -default_issymmetric() = true -default_solver_verbosity() = 0 - -# Customizable based on solver function -default_solver_outputlevel(::Function) = 0 - -default_solver_tol(::Function) = error("Not implemented") -default_solver_which_eigenvalue(::Function) = error("Not implemented") -default_solver_krylovdim(::Function) = error("Not implemented") -default_solver_verbosity(::Function) = error("Not implemented") - -## Solver-specific keyword argument defaults - -# dmrg/eigsolve -default_solver_tol(::typeof(eigsolve)) = 1e-14 -default_solver_krylovdim(::typeof(eigsolve)) = 3 -default_solver_maxiter(::typeof(eigsolve)) = 1 -default_solver_which_eigenvalue(::typeof(eigsolve)) = :SR - -# tdvp/exponentiate -default_solver_tol(::typeof(exponentiate)) = 1e-12 -default_solver_krylovdim(::typeof(exponentiate)) = 30 -default_solver_maxiter(::typeof(exponentiate)) = 100 diff --git a/src/dmrg.jl b/src/dmrg.jl deleted file mode 100644 index f736f88..0000000 --- a/src/dmrg.jl +++ /dev/null @@ -1,58 +0,0 @@ -function dmrg_solver( - f::typeof(eigsolve); - solver_which_eigenvalue, - ishermitian, - solver_tol, - solver_krylovdim, - solver_maxiter, - solver_verbosity, -) - function solver(H, t, psi0; current_time, outputlevel) - howmany = 1 - which = solver_which_eigenvalue - vals, vecs, info = f( - H, - psi0, - howmany, - which; - ishermitian=default_ishermitian(), - tol=solver_tol, - krylovdim=solver_krylovdim, - maxiter=solver_maxiter, - verbosity=solver_verbosity, - ) - psi = vecs[1] - return psi, info - end - return solver -end - -function dmrg( - H, - psi0::MPS; - ishermitian=default_ishermitian(), - solver_which_eigenvalue=default_solver_which_eigenvalue(eigsolve), - solver_tol=default_solver_tol(eigsolve), - solver_krylovdim=default_solver_krylovdim(eigsolve), - solver_maxiter=default_solver_maxiter(eigsolve), - solver_verbosity=default_solver_verbosity(), - kwargs..., -) - reverse_step = false - psi = alternating_update( - dmrg_solver( - eigsolve; - solver_which_eigenvalue, - ishermitian, - solver_tol, - solver_krylovdim, - solver_maxiter, - solver_verbosity, - ), - H, - psi0; - reverse_step, - kwargs..., - ) - return psi -end diff --git a/src/dmrg_x.jl b/src/dmrg_x.jl deleted file mode 100644 index 7a82deb..0000000 --- a/src/dmrg_x.jl +++ /dev/null @@ -1,16 +0,0 @@ -using ITensors: MPS, array, contract, dag, uniqueind, onehot -using LinearAlgebra: eigen - -function dmrg_x_solver(PH, t, psi0; current_time, outputlevel) - H = contract(PH, ITensor(true)) - D, U = eigen(H; ishermitian=true) - u = uniqueind(U, H) - max_overlap, max_ind = findmax(abs, array(psi0 * dag(U))) - U_max = U * dag(onehot(eltype(U), u => max_ind)) - return U_max, nothing -end - -function dmrg_x(PH, psi0::MPS; reverse_step=false, kwargs...) - psi = alternating_update(dmrg_x_solver, PH, psi0; reverse_step, kwargs...) - return psi -end diff --git a/src/linsolve.jl b/src/linsolve.jl deleted file mode 100644 index 88dd403..0000000 --- a/src/linsolve.jl +++ /dev/null @@ -1,42 +0,0 @@ -using ITensors: MPO, MPS -using KrylovKit: KrylovKit, linsolve - -""" -Compute a solution x to the linear system: - -(a₀ + a₁ * A)*x = b - -using starting guess x₀. Leaving a₀, a₁ -set to their default values solves the -system A*x = b. - -To adjust the balance between accuracy of solution -and speed of the algorithm, it is recommed to first try -adjusting the solver keyword arguments as descibed below. - -Keyword arguments: - - `nsweeps`, `cutoff`, `maxdim`, etc. (like for other MPO/MPS solvers). - - `solver_kwargs=(;)` - a `NamedTuple` containing keyword arguments that will get forwarded to the local solver, - in this case `KrylovKit.linsolve` which is a GMRES linear solver. For example: - ```juli - linsolve(A, b, x; maxdim=100, cutoff=1e-8, nsweeps=10, solver_kwargs=(; ishermitian=true, tol=1e-6, maxiter=20, krylovdim=30)) - ``` - See `KrylovKit.jl` documentation for more details on available keyword arguments. -""" -function KrylovKit.linsolve( - A::MPO, - b::MPS, - x₀::MPS, - a₀::Number=false, - a₁::Number=true; - solver_kwargs=(;), - tdvp_kwargs..., -) - function linsolve_solver(P::ProjMPO_MPS2, t, x₀; current_time, outputlevel) - b = dag(only(proj_mps(P))) - x, info = linsolve(P, b, x₀, a₀, a₁; solver_kwargs...) - return x, nothing - end - P = ProjMPO_MPS2(A, b) - return alternating_update(linsolve_solver, P, x₀; reverse_step=false, tdvp_kwargs...) -end diff --git a/src/projmpo_apply.jl b/src/projmpo_apply.jl deleted file mode 100644 index 0493e48..0000000 --- a/src/projmpo_apply.jl +++ /dev/null @@ -1,90 +0,0 @@ -using ITensors: ITensor -using ITensors.ITensorMPS: ITensorMPS, AbstractProjMPO, MPO, MPS - -""" -A ProjMPOApply represents the application of an -MPO `H` onto an MPS `psi0` but "projected" by -the basis of a different MPS `psi` (which -could be an approximation to H|psi>). - -As an implementation of the AbstractProjMPO -type, it supports multiple `nsite` values for -one- and two-site algorithms. - -``` - *--*--*- -*--*--*--*--*--* -``` -""" -mutable struct ProjMPOApply <: AbstractProjMPO - lpos::Int - rpos::Int - nsite::Int - psi0::MPS - H::MPO - LR::Vector{ITensor} -end - -function ProjMPOApply(psi0::MPS, H::MPO) - return ProjMPOApply(0, length(H) + 1, 2, psi0, H, Vector{ITensor}(undef, length(H))) -end - -function Base.copy(P::ProjMPOApply) - return ProjMPOApply(P.lpos, P.rpos, P.nsite, copy(P.psi0), copy(P.H), copy(P.LR)) -end - -function ITensorMPS.set_nsite!(P::ProjMPOApply, nsite) - P.nsite = nsite - return P -end - -function ITensorMPS.makeL!(P::ProjMPOApply, psi::MPS, k::Int) - # Save the last `L` that is made to help with caching - # for DiskProjMPO - ll = P.lpos - if ll ≥ k - # Special case when nothing has to be done. - # Still need to change the position if lproj is - # being moved backward. - P.lpos = k - return nothing - end - # Make sure ll is at least 0 for the generic logic below - ll = max(ll, 0) - L = lproj(P) - while ll < k - L = L * P.psi0[ll + 1] * P.H[ll + 1] * dag(psi[ll + 1]) - P.LR[ll + 1] = L - ll += 1 - end - # Needed when moving lproj backward. - P.lpos = k - return P -end - -function ITensorMPS.makeR!(P::ProjMPOApply, psi::MPS, k::Int) - # Save the last `R` that is made to help with caching - # for DiskProjMPO - rl = P.rpos - if rl ≤ k - # Special case when nothing has to be done. - # Still need to change the position if rproj is - # being moved backward. - P.rpos = k - return nothing - end - N = length(P.H) - # Make sure rl is no bigger than `N + 1` for the generic logic below - rl = min(rl, N + 1) - R = rproj(P) - while rl > k - R = R * P.psi0[rl - 1] * P.H[rl - 1] * dag(psi[rl - 1]) - P.LR[rl - 1] = R - rl -= 1 - end - P.rpos = k - return P -end diff --git a/src/projmpo_mps2.jl b/src/projmpo_mps2.jl deleted file mode 100644 index bb9801d..0000000 --- a/src/projmpo_mps2.jl +++ /dev/null @@ -1,47 +0,0 @@ -using ITensors: contract -using ITensors.ITensorMPS: AbstractProjMPO, ProjMPO, makeL!, makeR!, nsite, set_nsite! - -mutable struct ProjMPO_MPS2 <: AbstractProjMPO - PH::ProjMPO - Ms::Vector{ProjMPS2} -end - -function ProjMPO_MPS2(H::MPO, M::MPS) - return ProjMPO_MPS2(ProjMPO(H), [ProjMPS2(M)]) -end - -function ProjMPO_MPS2(H::MPO, Mv::Vector{MPS}) - return ProjMPO_MPS2(ProjMPO(H), [ProjMPS2(m) for m in Mv]) -end - -Base.copy(P::ProjMPO_MPS2) = ProjMPO_MPS2(copy(P.PH), copy(P.Ms)) - -ITensorMPS.nsite(P::ProjMPO_MPS2) = nsite(P.PH) - -function ITensorMPS.set_nsite!(P::ProjMPO_MPS2, nsite) - set_nsite!(P.PH, nsite) - for m in P.Ms - set_nsite!(m, nsite) - end - return P -end - -function ITensorMPS.makeL!(P::ProjMPO_MPS2, psi::MPS, k::Int) - makeL!(P.PH, psi, k) - for m in P.Ms - makeL!(m, psi, k) - end - return P -end - -function ITensorMPS.makeR!(P::ProjMPO_MPS2, psi::MPS, k::Int) - makeR!(P.PH, psi, k) - for m in P.Ms - makeR!(m, psi, k) - end - return P -end - -ITensors.contract(P::ProjMPO_MPS2, v::ITensor) = contract(P.PH, v) - -proj_mps(P::ProjMPO_MPS2) = [proj_mps(m) for m in P.Ms] diff --git a/src/projmps2.jl b/src/projmps2.jl deleted file mode 100644 index c324529..0000000 --- a/src/projmps2.jl +++ /dev/null @@ -1,123 +0,0 @@ -using ITensors: ITensors, ITensor, dag, dim, prime -using ITensors.ITensorMPS: ITensorMPS, AbstractProjMPO, OneITensor, lproj, rproj, site_range - -""" -Holds the following data where psi -is the MPS being optimized and M is the -MPS held constant by the ProjMPS. -``` - o--o--o--o--o--o--o--o--o--o--o -``` -""" -mutable struct ProjMPS2 <: AbstractProjMPO - lpos::Int - rpos::Int - nsite::Int - M::MPS - LR::Vector{ITensor} -end - -function ProjMPS2(M::MPS) - return ProjMPS2(0, length(M) + 1, 2, M, Vector{ITensor}(undef, length(M))) -end - -Base.length(P::ProjMPS2) = length(P.M) - -function Base.copy(P::ProjMPS2) - return ProjMPS2(P.lpos, P.rpos, P.nsite, copy(P.M), copy(P.LR)) -end - -function ITensorMPS.set_nsite!(P::ProjMPS2, nsite) - P.nsite = nsite - return P -end - -function ITensorMPS.makeL!(P::ProjMPS2, psi::MPS, k::Int) - # Save the last `L` that is made to help with caching - # for DiskProjMPO - ll = P.lpos - if ll ≥ k - # Special case when nothing has to be done. - # Still need to change the position if lproj is - # being moved backward. - P.lpos = k - return nothing - end - # Make sure ll is at least 0 for the generic logic below - ll = max(ll, 0) - L = lproj(P) - while ll < k - L = L * psi[ll + 1] * dag(prime(P.M[ll + 1], "Link")) - P.LR[ll + 1] = L - ll += 1 - end - # Needed when moving lproj backward. - P.lpos = k - return P -end - -function ITensorMPS.makeR!(P::ProjMPS2, psi::MPS, k::Int) - # Save the last `R` that is made to help with caching - # for DiskProjMPO - rl = P.rpos - if rl ≤ k - # Special case when nothing has to be done. - # Still need to change the position if rproj is - # being moved backward. - P.rpos = k - return nothing - end - N = length(P.M) - # Make sure rl is no bigger than `N + 1` for the generic logic below - rl = min(rl, N + 1) - R = rproj(P) - while rl > k - R = R * psi[rl - 1] * dag(prime(P.M[rl - 1], "Link")) - P.LR[rl - 1] = R - rl -= 1 - end - P.rpos = k - return P -end - -function ITensors.contract(P::ProjMPS2, v::ITensor) - itensor_map = Union{ITensor,OneITensor}[lproj(P)] - append!(itensor_map, [prime(t, "Link") for t in P.M[site_range(P)]]) - push!(itensor_map, rproj(P)) - - # Reverse the contraction order of the map if - # the first tensor is a scalar (for example we - # are at the left edge of the system) - if dim(first(itensor_map)) == 1 - reverse!(itensor_map) - end - - # Apply the map - Mv = v - for it in itensor_map - Mv *= it - end - return Mv -end - -function proj_mps(P::ProjMPS2) - itensor_map = Union{ITensor,OneITensor}[lproj(P)] - append!(itensor_map, [dag(prime(t, "Link")) for t in P.M[site_range(P)]]) - push!(itensor_map, rproj(P)) - - # Reverse the contraction order of the map if - # the first tensor is a scalar (for example we - # are at the left edge of the system) - if dim(first(itensor_map)) == 1 - reverse!(itensor_map) - end - - # Apply the map - m = ITensor(true) - for it in itensor_map - m *= it - end - return m -end diff --git a/src/solver_utils.jl b/src/solver_utils.jl deleted file mode 100644 index f0a7942..0000000 --- a/src/solver_utils.jl +++ /dev/null @@ -1,69 +0,0 @@ -using ITensors: ITensor, ProjMPOSum, apply, array, inds, itensor, permute - -# Utilities for making it easier -# to define solvers (like ODE solvers) -# for TDVP - -""" - to_vec(x) -Transform `x` into a `Vector`, and return the vector, and a closure which inverts the -transformation. - -Modeled after FiniteDifferences.to_vec: - -https://github.com/JuliaDiff/FiniteDifferences.jl/blob/main/src/to_vec.jl -""" -to_vec(x) = error("Not implemented") - -function to_vec(x::ITensor) - function to_itensor(x_vec) - return itensor(x_vec, inds(x)) - end - return vec(array(x)), to_itensor -end - -# Represents a time-dependent sum of terms: -# -# H(t) = f[1](t) * H0[1] + f[2](t) * H0[2] + … -# -struct TimeDependentSum{S,T} - f::Vector{S} - H0::T -end -TimeDependentSum(f::Vector, H0::ProjMPOSum) = TimeDependentSum(f, ITensors.terms(H0)) -Base.length(H::TimeDependentSum) = length(H.f) - -function Base.:*(c::Number, H::TimeDependentSum) - return TimeDependentSum([t -> c * fₙ(t) for fₙ in H.f], H.H0) -end -Base.:*(H::TimeDependentSum, c::Number) = c * H - -# Calling a `TimeDependentOpSum` at a certain time like: -# -# H(t) -# -# Returns a `ScaledSum` at that time. -(H::TimeDependentSum)(t::Number) = ScaledSum([fₙ(t) for fₙ in H.f], H.H0) - -# Represents the sum of scaled terms: -# -# H = coefficient[1] * H[1] + coefficient * H[2] + … -# -struct ScaledSum{S,T} - coefficients::Vector{S} - H::T -end -Base.length(H::ScaledSum) = length(H.coefficients) - -# Apply the scaled sum of terms: -# -# H(ψ₀) = coefficient[1] * H[1](ψ₀) + coefficient[2] * H[2](ψ₀) + … -# -# onto ψ₀. -function (H::ScaledSum)(ψ₀) - ψ = ITensor(inds(ψ₀)) - for n in 1:length(H) - ψ += H.coefficients[n] * apply(H.H[n], ψ₀) - end - return permute(ψ, inds(ψ₀)) -end diff --git a/src/sweep_update.jl b/src/sweep_update.jl deleted file mode 100644 index 3679e7a..0000000 --- a/src/sweep_update.jl +++ /dev/null @@ -1,434 +0,0 @@ -using ITensors: uniqueinds -using ITensors.ITensorMPS: - ITensorMPS, MPS, isortho, orthocenter, orthogonalize!, position!, replacebond!, set_nsite! -using LinearAlgebra: norm, normalize!, svd -using Observers: update! -using Printf: @printf - -function sweep_update( - order::TDVPOrder, solver, PH, time_step::Number, psi::MPS; current_time=0, kwargs... -) - order_orderings = orderings(order) - order_sub_time_steps = eltype(time_step).(sub_time_steps(order)) - order_sub_time_steps *= time_step - info = nothing - for substep in 1:length(order_sub_time_steps) - psi, PH, info = sub_sweep_update( - order_orderings[substep], - solver, - PH, - order_sub_time_steps[substep], - psi; - current_time, - kwargs..., - ) - current_time += order_sub_time_steps[substep] - end - return psi, PH, info -end - -isforward(direction::Base.ForwardOrdering) = true -isforward(direction::Base.ReverseOrdering) = false -isreverse(direction) = !isforward(direction) - -function sweep_bonds(direction::Base.ForwardOrdering, n::Int; ncenter::Int) - return 1:(n - ncenter + 1) -end - -function sweep_bonds(direction::Base.ReverseOrdering, n::Int; ncenter::Int) - return reverse(sweep_bonds(Base.Forward, n; ncenter)) -end - -is_forward_done(direction::Base.ForwardOrdering, b, n; ncenter) = (b + ncenter - 1 == n) -is_forward_done(direction::Base.ReverseOrdering, b, n; ncenter) = false -is_reverse_done(direction::Base.ForwardOrdering, b, n; ncenter) = false -is_reverse_done(direction::Base.ReverseOrdering, b, n; ncenter) = (b == 1) -function is_half_sweep_done(direction, b, n; ncenter) - return is_forward_done(direction, b, n; ncenter) || - is_reverse_done(direction, b, n; ncenter) -end - -function sub_sweep_update( - direction::Base.Ordering, - solver, - PH, - time_step::Number, - psi::MPS; - which_decomp=nothing, - svd_alg=nothing, - sweep=default_sweep(), - current_time=default_current_time(), - nsite=default_nsite(), - reverse_step=default_reverse_step(), - normalize=default_normalize(), - (observer!)=default_observer!(), - outputlevel=default_outputlevel(), - maxdim=default_maxdim(), - mindim=default_mindim(), - cutoff=default_cutoff(time_step), - noise=default_noise(), -) - PH = copy(PH) - psi = copy(psi) - if length(psi) == 1 - error( - "`tdvp`, `dmrg`, `linsolve`, etc. currently does not support system sizes of 1. You can diagonalize the MPO tensor directly with tools like `LinearAlgebra.eigen`, `KrylovKit.exponentiate`, etc.", - ) - end - N = length(psi) - set_nsite!(PH, nsite) - if isforward(direction) - if !isortho(psi) || orthocenter(psi) != 1 - orthogonalize!(psi, 1) - end - @assert isortho(psi) && orthocenter(psi) == 1 - position!(PH, psi, 1) - elseif isreverse(direction) - if !isortho(psi) || orthocenter(psi) != N - nsite + 1 - orthogonalize!(psi, N - nsite + 1) - end - @assert(isortho(psi) && (orthocenter(psi) == N - nsite + 1)) - position!(PH, psi, N - nsite + 1) - end - maxtruncerr = 0.0 - info = nothing - for b in sweep_bonds(direction, N; ncenter=nsite) - current_time, maxtruncerr, spec, info = region_update!( - solver, - PH, - psi, - b; - nsite, - reverse_step, - current_time, - outputlevel, - time_step, - normalize, - direction, - noise, - which_decomp, - svd_alg, - cutoff, - maxdim, - mindim, - maxtruncerr, - ) - if outputlevel >= 2 - if nsite == 1 - @printf("Sweep %d, direction %s, bond (%d,) \n", sweep, direction, b) - elseif nsite == 2 - @printf("Sweep %d, direction %s, bond (%d,%d) \n", sweep, direction, b, b + 1) - end - print(" Truncated using") - @printf(" cutoff=%.1E", cutoff) - @printf(" maxdim=%.1E", maxdim) - print(" mindim=", mindim) - print(" current_time=", round(current_time; digits=3)) - println() - if spec != nothing - @printf( - " Trunc. err=%.2E, bond dimension %d\n", spec.truncerr, dim(linkind(psi, b)) - ) - end - flush(stdout) - end - update!( - observer!; - psi, - bond=b, - sweep, - half_sweep=isforward(direction) ? 1 : 2, - spec, - outputlevel, - half_sweep_is_done=is_half_sweep_done(direction, b, N; ncenter=nsite), - current_time, - info, - ) - end - # Just to be sure: - normalize && normalize!(psi) - return psi, PH, TDVPInfo(maxtruncerr) -end - -function region_update!( - solver, - PH, - psi, - b; - nsite, - reverse_step, - current_time, - outputlevel, - time_step, - normalize, - direction, - noise, - which_decomp, - svd_alg, - cutoff, - maxdim, - mindim, - maxtruncerr, -) - return region_update!( - Val(nsite), - Val(reverse_step), - solver, - PH, - psi, - b; - current_time, - outputlevel, - time_step, - normalize, - direction, - noise, - which_decomp, - svd_alg, - cutoff, - maxdim, - mindim, - maxtruncerr, - ) -end - -function region_update!( - nsite_val::Val{1}, - reverse_step_val::Val{false}, - solver, - PH, - psi, - b; - current_time, - outputlevel, - time_step, - normalize, - direction, - noise, - which_decomp, - svd_alg, - cutoff, - maxdim, - mindim, - maxtruncerr, -) - N = length(psi) - nsite = 1 - # Do 'forwards' evolution step - set_nsite!(PH, nsite) - position!(PH, psi, b) - phi1 = psi[b] - phi1, info = solver(PH, time_step, phi1; current_time, outputlevel) - current_time += time_step - normalize && (phi1 /= norm(phi1)) - spec = nothing - psi[b] = phi1 - if !is_half_sweep_done(direction, b, N; ncenter=nsite) - # Move ortho center - Δ = (isforward(direction) ? +1 : -1) - orthogonalize!(psi, b + Δ) - end - return current_time, maxtruncerr, spec, info -end - -function region_update!( - nsite_val::Val{1}, - reverse_step_val::Val{true}, - solver, - PH, - psi, - b; - current_time, - outputlevel, - time_step, - normalize, - direction, - noise, - which_decomp, - svd_alg, - cutoff, - maxdim, - mindim, - maxtruncerr, -) - N = length(psi) - nsite = 1 - # Do 'forwards' evolution step - set_nsite!(PH, nsite) - position!(PH, psi, b) - phi1 = psi[b] - phi1, info = solver(PH, time_step, phi1; current_time, outputlevel) - current_time += time_step - normalize && (phi1 /= norm(phi1)) - spec = nothing - psi[b] = phi1 - if !is_half_sweep_done(direction, b, N; ncenter=nsite) - # Do backwards evolution step - b1 = (isforward(direction) ? b + 1 : b) - Δ = (isforward(direction) ? +1 : -1) - uinds = uniqueinds(phi1, psi[b + Δ]) - U, S, V = svd(phi1, uinds) - psi[b] = U - phi0 = S * V - if isforward(direction) - ITensorMPS.setleftlim!(psi, b) - elseif isreverse(direction) - ITensorMPS.setrightlim!(psi, b) - end - set_nsite!(PH, nsite - 1) - position!(PH, psi, b1) - phi0, info = solver(PH, -time_step, phi0; current_time, outputlevel) - current_time -= time_step - normalize && (phi0 ./= norm(phi0)) - psi[b + Δ] = phi0 * psi[b + Δ] - if isforward(direction) - ITensorMPS.setrightlim!(psi, b + Δ + 1) - elseif isreverse(direction) - ITensorMPS.setleftlim!(psi, b + Δ - 1) - end - set_nsite!(PH, nsite) - end - return current_time, maxtruncerr, spec, info -end - -function region_update!( - nsite_val::Val{2}, - reverse_step_val::Val{false}, - solver, - PH, - psi, - b; - current_time, - outputlevel, - time_step, - normalize, - direction, - noise, - which_decomp, - svd_alg, - cutoff, - maxdim, - mindim, - maxtruncerr, -) - N = length(psi) - nsite = 2 - # Do 'forwards' evolution step - set_nsite!(PH, nsite) - position!(PH, psi, b) - phi1 = psi[b] * psi[b + 1] - phi1, info = solver(PH, time_step, phi1; current_time, outputlevel) - current_time += time_step - normalize && (phi1 /= norm(phi1)) - spec = nothing - ortho = isforward(direction) ? "left" : "right" - drho = nothing - if noise > 0.0 && isforward(direction) - drho = noise * noiseterm(PH, phi, ortho) - end - spec = replacebond!( - psi, - b, - phi1; - maxdim, - mindim, - cutoff, - eigen_perturbation=drho, - ortho=ortho, - normalize, - which_decomp, - svd_alg, - ) - maxtruncerr = max(maxtruncerr, spec.truncerr) - return current_time, maxtruncerr, spec, info -end - -function region_update!( - nsite_val::Val{2}, - reverse_step_val::Val{true}, - solver, - PH, - psi, - b; - current_time, - outputlevel, - time_step, - normalize, - direction, - noise, - which_decomp, - svd_alg, - cutoff, - maxdim, - mindim, - maxtruncerr, -) - N = length(psi) - nsite = 2 - # Do 'forwards' evolution step - set_nsite!(PH, nsite) - position!(PH, psi, b) - phi1 = psi[b] * psi[b + 1] - phi1, info = solver(PH, time_step, phi1; current_time, outputlevel) - current_time += time_step - normalize && (phi1 /= norm(phi1)) - spec = nothing - ortho = isforward(direction) ? "left" : "right" - drho = nothing - if noise > 0.0 && isforward(direction) - drho = noise * noiseterm(PH, phi, ortho) - end - spec = replacebond!( - psi, - b, - phi1; - maxdim, - mindim, - cutoff, - eigen_perturbation=drho, - ortho=ortho, - normalize, - which_decomp, - svd_alg, - ) - maxtruncerr = max(maxtruncerr, spec.truncerr) - if !is_half_sweep_done(direction, b, N; ncenter=nsite) - # Do backwards evolution step - b1 = (isforward(direction) ? b + 1 : b) - Δ = (isforward(direction) ? +1 : -1) - phi0 = psi[b1] - set_nsite!(PH, nsite - 1) - position!(PH, psi, b1) - phi0, info = solver(PH, -time_step, phi0; current_time, outputlevel) - current_time -= time_step - normalize && (phi0 ./= norm(phi0)) - psi[b1] = phi0 - set_nsite!(PH, nsite) - end - return current_time, maxtruncerr, spec, info -end - -function region_update!( - ::Val{nsite}, - ::Val{reverse_step}, - solver, - PH, - psi, - b; - current_time, - outputlevel, - time_step, - normalize, - direction, - noise, - which_decomp, - svd_alg, - cutoff, - maxdim, - mindim, - maxtruncerr, -) where {nsite,reverse_step} - return error( - "`tdvp`, `dmrg`, `linsolve`, etc. with `nsite=$nsite` and `reverse_step=$reverse_step` not implemented.", - ) -end diff --git a/src/tdvp.jl b/src/tdvp.jl deleted file mode 100644 index 7c40795..0000000 --- a/src/tdvp.jl +++ /dev/null @@ -1,136 +0,0 @@ -using ITensors: Algorithm, MPO, MPS, @Algorithm_str - -# Select solver function -solver_function(solver_backend::String) = solver_function(Algorithm(solver_backend)) -solver_function(::Algorithm"exponentiate") = exponentiate -function solver_function(solver_backend::Algorithm) - return error( - "solver_backend=$(String(solver_backend)) not recognized (only \"exponentiate\" is supported)", - ) -end - -# Kept for backwards compatibility -function solver_function(::Algorithm"applyexp") - println( - "Warning: the `solver_backend` option `\"applyexp\"` in `tdvp` has been removed. `\"exponentiate\"` will be used instead. To remove this warning, don't specify the `solver_backend` keyword argument.", - ) - return solver_function(Algorithm"exponentiate"()) -end - -function tdvp_solver( - f::typeof(exponentiate); - ishermitian, - issymmetric, - solver_tol, - solver_krylovdim, - solver_maxiter, - solver_outputlevel, -) - function solver(H, t, psi0; current_time, outputlevel) - psi, info = f( - H, - t, - psi0; - ishermitian, - issymmetric, - tol=solver_tol, - krylovdim=solver_krylovdim, - maxiter=solver_maxiter, - verbosity=solver_outputlevel, - eager=true, - ) - return psi, info - end - return solver -end - -function tdvp( - H, - t::Number, - psi0::MPS; - ishermitian=default_ishermitian(), - issymmetric=default_issymmetric(), - solver_backend=default_tdvp_solver_backend(), - solver_function=solver_function(solver_backend), - solver_tol=default_solver_tol(solver_function), - solver_krylovdim=default_solver_krylovdim(solver_function), - solver_maxiter=default_solver_maxiter(solver_function), - solver_outputlevel=default_solver_outputlevel(solver_function), - kwargs..., -) - return tdvp( - tdvp_solver( - solver_function; - ishermitian, - issymmetric, - solver_tol, - solver_krylovdim, - solver_maxiter, - solver_outputlevel, - ), - H, - t, - psi0; - kwargs..., - ) -end - -function tdvp(t::Number, H, psi0::MPS; kwargs...) - return tdvp(H, t, psi0; kwargs...) -end - -function tdvp(H, psi0::MPS, t::Number; kwargs...) - return tdvp(H, t, psi0; kwargs...) -end - -""" - tdvp(H::MPO,psi0::MPS,t::Number; kwargs...) - tdvp(H::MPO,psi0::MPS,t::Number; kwargs...) - -Use the time dependent variational principle (TDVP) algorithm -to compute `exp(t*H)*psi0` using an efficient algorithm based -on alternating optimization of the MPS tensors and local Krylov -exponentiation of H. - -Returns: -* `psi::MPS` - time-evolved MPS - -Optional keyword arguments: -* `outputlevel::Int = 1` - larger outputlevel values resulting in printing more information and 0 means no output -* `observer` - object implementing the [Observer](@ref observer) interface which can perform measurements and stop early -* `write_when_maxdim_exceeds::Int` - when the allowed maxdim exceeds this value, begin saving tensors to disk to free memory in large calculations -""" -function tdvp(solver, H::MPO, t::Number, psi0::MPS; kwargs...) - return alternating_update(solver, H, t, psi0; kwargs...) -end - -function tdvp(solver, t::Number, H, psi0::MPS; kwargs...) - return tdvp(solver, H, t, psi0; kwargs...) -end - -function tdvp(solver, H, psi0::MPS, t::Number; kwargs...) - return tdvp(solver, H, t, psi0; kwargs...) -end - -""" - tdvp(Hs::Vector{MPO},psi0::MPS,t::Number; kwargs...) - tdvp(Hs::Vector{MPO},psi0::MPS,t::Number, sweeps::Sweeps; kwargs...) - -Use the time dependent variational principle (TDVP) algorithm -to compute `exp(t*H)*psi0` using an efficient algorithm based -on alternating optimization of the MPS tensors and local Krylov -exponentiation of H. - -This version of `tdvp` accepts a representation of H as a -Vector of MPOs, Hs = [H1,H2,H3,...] such that H is defined -as H = H1+H2+H3+... -Note that this sum of MPOs is not actually computed; rather -the set of MPOs [H1,H2,H3,..] is efficiently looped over at -each step of the algorithm when optimizing the MPS. - -Returns: -* `psi::MPS` - time-evolved MPS -""" -function tdvp(solver, Hs::Vector{MPO}, t::Number, psi0::MPS; kwargs...) - return alternating_update(solver, Hs, t, psi0; kwargs...) -end diff --git a/src/tdvp_sweeps.jl b/src/tdvp_sweeps.jl deleted file mode 100644 index d39ea92..0000000 --- a/src/tdvp_sweeps.jl +++ /dev/null @@ -1,17 +0,0 @@ -function process_sweeps(s::Sweeps) - return (; - nsweeps=s.nsweep, maxdim=s.maxdim, mindim=s.mindim, cutoff=s.cutoff, noise=s.noise - ) -end - -function tdvp(H, t::Number, psi0::MPS, sweeps::Sweeps; kwargs...) - return tdvp(H, t, psi0; process_sweeps(sweeps)..., kwargs...) -end - -function tdvp(solver, H, t::Number, psi0::MPS, sweeps::Sweeps; kwargs...) - return tdvp(solver, H, t, psi0; process_sweeps(sweeps)..., kwargs...) -end - -function dmrg(H, psi0::MPS, sweeps::Sweeps; kwargs...) - return dmrg(H, psi0; process_sweeps(sweeps)..., kwargs...) -end diff --git a/src/tdvpinfo.jl b/src/tdvpinfo.jl deleted file mode 100644 index fb4c128..0000000 --- a/src/tdvpinfo.jl +++ /dev/null @@ -1,7 +0,0 @@ -""" -#fields -- `maxtruncerr::Float64`: the maximum tuncation error -""" -struct TDVPInfo - maxtruncerr::Float64 -end diff --git a/src/tdvporder.jl b/src/tdvporder.jl deleted file mode 100644 index e37e00e..0000000 --- a/src/tdvporder.jl +++ /dev/null @@ -1,24 +0,0 @@ -struct TDVPOrder{order,direction} end - -TDVPOrder(order::Int, direction::Base.Ordering) = TDVPOrder{order,direction}() - -orderings(::TDVPOrder) = error("Not implemented") -sub_time_steps(::TDVPOrder) = error("Not implemented") - -function orderings(::TDVPOrder{1,direction}) where {direction} - return [direction, Base.ReverseOrdering(direction)] -end -sub_time_steps(::TDVPOrder{1}) = [1, 0] - -function orderings(::TDVPOrder{2,direction}) where {direction} - return [direction, Base.ReverseOrdering(direction)] -end -sub_time_steps(::TDVPOrder{2}) = [1 / 2, 1 / 2] - -function orderings(::TDVPOrder{4,direction}) where {direction} - return [direction, Base.ReverseOrdering(direction)] -end -function sub_time_steps(::TDVPOrder{4}) - s = 1 / (2 - 2^(1 / 3)) - return [s / 2, s / 2, (1 - 2 * s) / 2, (1 - 2 * s) / 2, s / 2, s / 2] -end diff --git a/src/update_observer.jl b/src/update_observer.jl deleted file mode 100644 index 3e75f9a..0000000 --- a/src/update_observer.jl +++ /dev/null @@ -1,6 +0,0 @@ -using ITensors: AbstractObserver, measure! -using Observers: Observers - -function Observers.update!(observer::AbstractObserver; kwargs...) - return measure!(observer; kwargs...) -end From b7ee90ba5e56026681f33abdc7a123aa0d272911 Mon Sep 17 00:00:00 2001 From: Miles Date: Wed, 27 Mar 2024 16:05:32 +0900 Subject: [PATCH 02/26] Update tests --- test/test_dmrg.jl | 4 ++-- test/test_dmrg_x.jl | 4 +++- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/test/test_dmrg.jl b/test/test_dmrg.jl index 159b73b..4414621 100644 --- a/test/test_dmrg.jl +++ b/test/test_dmrg.jl @@ -1,5 +1,5 @@ @eval module $(gensym()) -using ITensors: ITensors, MPO, OpSum, inner, randomMPS, siteinds +using ITensors: ITensorMPS, ITensors, MPO, OpSum, inner, randomMPS, siteinds using ITensorTDVP: ITensorTDVP using Test: @test, @testset @testset "DMRG (eltype=$elt, nsite=$nsite, conserve_qns=$conserve_qns)" for elt in ( @@ -24,7 +24,7 @@ using Test: @test, @testset psi = ITensorTDVP.dmrg( H, psi; nsweeps, maxdim, cutoff, nsite, solver_krylovdim=3, solver_maxiter=1 ) - e2, psi2 = ITensors.dmrg(H, psi; nsweeps, maxdim, cutoff, outputlevel=0) + e2, psi2 = ITensorMPS.dmrg(H, psi; nsweeps, maxdim, cutoff, outputlevel=0) @test ITensors.scalartype(psi2) == elt @test e2 isa real(elt) @test inner(psi', H, psi) ≈ inner(psi2', H, psi2) rtol = √(eps(real(elt))) * 10 diff --git a/test/test_dmrg_x.jl b/test/test_dmrg_x.jl index 8e6d7f4..4e21cb8 100644 --- a/test/test_dmrg_x.jl +++ b/test/test_dmrg_x.jl @@ -1,8 +1,10 @@ @eval module $(gensym()) -using ITensors: ITensors, MPO, MPS, OpSum, ProjMPO, inner, siteinds +using ITensors: ITensors, inner, OpSum, siteinds +using ITensors.ITensorMPS: MPO, MPS, ProjMPO using ITensorTDVP: dmrg_x using Random: Random using Test: @test, @testset + @testset "DMRG-X (eltype=$elt, conserve_qns=$conserve_qns)" for elt in ( Float32, Float64, Complex{Float32}, Complex{Float64} ), From 86afc5daba5a443e24e062cb7c94a970fcdc4874 Mon Sep 17 00:00:00 2001 From: Miles Date: Wed, 27 Mar 2024 16:12:17 +0900 Subject: [PATCH 03/26] Add Reexport to dependencies --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 29971f7..02dd134 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Observers = "338f10d5-c7f1-4033-a7d1-f9dec39bcaa0" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" [compat] From 344d485b8ca39a230d7db921890ec57ba21e62e5 Mon Sep 17 00:00:00 2001 From: Miles Date: Wed, 27 Mar 2024 16:20:20 +0900 Subject: [PATCH 04/26] Update test_tdvp to avoid warning --- test/test_tdvp.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/test/test_tdvp.jl b/test/test_tdvp.jl index 46b27da..d230dd6 100644 --- a/test/test_tdvp.jl +++ b/test/test_tdvp.jl @@ -41,8 +41,6 @@ using Test: @test, @testset @test ITensors.scalartype(ψ1) == complex(elt) @test ψ1 ≈ tdvp(-time_step, H, ψ0; nsweeps=1, cutoff, nsite=1) @test ψ1 ≈ tdvp(H, ψ0, -time_step; nsweeps=1, cutoff, nsite=1) - #Different backend solvers, default solver_backend = "exponentiate" - @test ψ1 ≈ tdvp(H, ψ0, -time_step; nsweeps=1, cutoff, nsite=1, solver_backend="applyexp") @test norm(ψ1) ≈ 1 rtol = eps(real(elt)) * 10^3 ## Should lose fidelity: #@test abs(inner(ψ0,ψ1)) < 0.9 From 4e5d5d9f6e0c65a09af72ae0856a473deb04b7a8 Mon Sep 17 00:00:00 2001 From: Miles Date: Wed, 27 Mar 2024 16:25:25 +0900 Subject: [PATCH 05/26] Update example code --- examples/05_utils.jl | 33 ++++++++++++++++----------------- 1 file changed, 16 insertions(+), 17 deletions(-) diff --git a/examples/05_utils.jl b/examples/05_utils.jl index cba0f75..89ff174 100644 --- a/examples/05_utils.jl +++ b/examples/05_utils.jl @@ -1,5 +1,4 @@ -using ITensors: MPS, maxlinkdim -using ITensorTDVP: ITensorTDVP +using ITensors.ITensorMPS: ITensorMPS, MPS, maxlinkdim using Observers: observer, update! using Printf: @printf @@ -12,22 +11,22 @@ function tdvp_nonuniform_timesteps( time_start=0.0, order=2, (step_observer!)=observer(), - maxdim=ITensorTDVP.default_maxdim(), - mindim=ITensorTDVP.default_mindim(), - cutoff=ITensorTDVP.default_cutoff(), - noise=ITensorTDVP.default_noise(), - outputlevel=ITensorTDVP.default_outputlevel(), + maxdim=ITensorMPS.default_maxdim(), + mindim=ITensorMPS.default_mindim(), + cutoff=ITensorMPS.default_cutoff(), + noise=ITensorMPS.default_noise(), + outputlevel=ITensorMPS.default_outputlevel(), kwargs..., ) nsweeps = length(time_steps) - maxdim, mindim, cutoff, noise = ITensorTDVP.process_sweeps(; + maxdim, mindim, cutoff, noise = ITensorMPS.process_sweeps(; nsweeps, maxdim, mindim, cutoff, noise ) - tdvp_order = ITensorTDVP.TDVPOrder(order, Base.Forward) + tdvp_order = ITensorMPS.TDVPOrder(order, Base.Forward) current_time = time_start for sw in 1:nsweeps sw_time = @elapsed begin - psi, PH, info = ITensorTDVP.sweep_update( + psi, PH, info = ITensorMPS.sweep_update( tdvp_order, solver, PH, @@ -61,16 +60,16 @@ end function tdvp_nonuniform_timesteps( H, psi::MPS; - ishermitian=ITensorTDVP.default_ishermitian(), - issymmetric=ITensorTDVP.default_issymmetric(), - solver_tol=ITensorTDVP.default_solver_tol(exponentiate), - solver_krylovdim=ITensorTDVP.default_solver_krylovdim(exponentiate), - solver_maxiter=ITensorTDVP.default_solver_maxiter(exponentiate), - solver_outputlevel=ITensorTDVP.default_solver_outputlevel(exponentiate), + ishermitian=ITensorMPS.default_ishermitian(), + issymmetric=ITensorMPS.default_issymmetric(), + solver_tol=ITensorMPS.default_solver_tol(exponentiate), + solver_krylovdim=ITensorMPS.default_solver_krylovdim(exponentiate), + solver_maxiter=ITensorMPS.default_solver_maxiter(exponentiate), + solver_outputlevel=ITensorMPS.default_solver_outputlevel(exponentiate), kwargs..., ) return tdvp_nonuniform_timesteps( - ITensorTDVP.tdvp_solver( + ITensorMPS.tdvp_solver( exponentiate; ishermitian, issymmetric, From 7ddf3474d0453d2b00dcd24f43d6cbaa64a6539d Mon Sep 17 00:00:00 2001 From: Miles Date: Fri, 29 Mar 2024 14:15:43 +0900 Subject: [PATCH 06/26] Update src/ITensorTDVP.jl Co-authored-by: Matt Fishman --- src/ITensorTDVP.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ITensorTDVP.jl b/src/ITensorTDVP.jl index 4566a84..9f9b255 100644 --- a/src/ITensorTDVP.jl +++ b/src/ITensorTDVP.jl @@ -3,7 +3,7 @@ module ITensorTDVP using Reexport: @reexport @reexport using ITensors.ITensorMPS: tdvp, dmrg_x, to_vec, TimeDependentSum, linsolve -using ITensors: ITensorMPS +using ITensors.ITensorMPS: ITensorMPS dmrg(args...; kwargs...) = ITensorMPS.alternate_dmrg(args...; kwargs...) end From f4d56fc436842e352a1c5e32e5bcac0e8f5ae241 Mon Sep 17 00:00:00 2001 From: Miles Date: Fri, 29 Mar 2024 14:15:58 +0900 Subject: [PATCH 07/26] Update test/test_dmrg.jl Co-authored-by: Matt Fishman --- test/test_dmrg.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/test_dmrg.jl b/test/test_dmrg.jl index 4414621..3127642 100644 --- a/test/test_dmrg.jl +++ b/test/test_dmrg.jl @@ -1,5 +1,6 @@ @eval module $(gensym()) -using ITensors: ITensorMPS, ITensors, MPO, OpSum, inner, randomMPS, siteinds +using ITensors: ITensors, inner +using ITensors.ITensorMPS: ITensorMPS, OpSum, MPO, randomMPS, siteinds using ITensorTDVP: ITensorTDVP using Test: @test, @testset @testset "DMRG (eltype=$elt, nsite=$nsite, conserve_qns=$conserve_qns)" for elt in ( From a8eeacc53fc60e5d596b7b3d5488d4d3f49d2cca Mon Sep 17 00:00:00 2001 From: Miles Date: Fri, 29 Mar 2024 14:16:06 +0900 Subject: [PATCH 08/26] Update test/test_dmrg_x.jl Co-authored-by: Matt Fishman --- test/test_dmrg_x.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_dmrg_x.jl b/test/test_dmrg_x.jl index 4e21cb8..4cac5d6 100644 --- a/test/test_dmrg_x.jl +++ b/test/test_dmrg_x.jl @@ -1,6 +1,6 @@ @eval module $(gensym()) -using ITensors: ITensors, inner, OpSum, siteinds -using ITensors.ITensorMPS: MPO, MPS, ProjMPO +using ITensors: ITensors, inner +using ITensors.ITensorMPS: MPO, MPS, OpSum, ProjMPO, siteinds using ITensorTDVP: dmrg_x using Random: Random using Test: @test, @testset From eadd062e17f278ba5dc2a43f5f48d13e363e9355 Mon Sep 17 00:00:00 2001 From: Miles Date: Fri, 29 Mar 2024 14:18:27 +0900 Subject: [PATCH 09/26] Restore test of applyexp solver backend --- test/test_tdvp.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/test_tdvp.jl b/test/test_tdvp.jl index d230dd6..02784a2 100644 --- a/test/test_tdvp.jl +++ b/test/test_tdvp.jl @@ -41,6 +41,9 @@ using Test: @test, @testset @test ITensors.scalartype(ψ1) == complex(elt) @test ψ1 ≈ tdvp(-time_step, H, ψ0; nsweeps=1, cutoff, nsite=1) @test ψ1 ≈ tdvp(H, ψ0, -time_step; nsweeps=1, cutoff, nsite=1) + #Different backend solvers, default solver_backend = "exponentiate" + #May give a warning but it is intended + @test ψ1 ≈ tdvp(H, ψ0, -time_step; nsweeps=1, cutoff, nsite=1, solver_backend="applyexp") @test norm(ψ1) ≈ 1 rtol = eps(real(elt)) * 10^3 ## Should lose fidelity: #@test abs(inner(ψ0,ψ1)) < 0.9 From 80b4da794275536c346d600ae9ac4f0c6638ae26 Mon Sep 17 00:00:00 2001 From: Matt Fishman Date: Sat, 30 Mar 2024 12:26:13 -0400 Subject: [PATCH 10/26] Bump ITensors compat to v0.3.59 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 02dd134..29ffa9b 100644 --- a/Project.toml +++ b/Project.toml @@ -13,7 +13,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" [compat] -ITensors = "0.3.58" +ITensors = "0.3.59" KrylovKit = "0.6" Observers = "0.2" TimerOutputs = "0.5" From 91cd565cbcda340be27a634982c1879d161bb883 Mon Sep 17 00:00:00 2001 From: Miles Date: Sun, 31 Mar 2024 11:51:35 -0400 Subject: [PATCH 11/26] Reflect rename to ITensorMPS.alternating_update_dmrg --- src/ITensorTDVP.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ITensorTDVP.jl b/src/ITensorTDVP.jl index 9f9b255..2047e64 100644 --- a/src/ITensorTDVP.jl +++ b/src/ITensorTDVP.jl @@ -4,6 +4,6 @@ using Reexport: @reexport @reexport using ITensors.ITensorMPS: tdvp, dmrg_x, to_vec, TimeDependentSum, linsolve using ITensors.ITensorMPS: ITensorMPS -dmrg(args...; kwargs...) = ITensorMPS.alternate_dmrg(args...; kwargs...) +dmrg(args...; kwargs...) = ITensorMPS.alternating_update_dmrg(args...; kwargs...) end From 261ff90a22d0aafad5d93d1604d088b1a7c39929 Mon Sep 17 00:00:00 2001 From: Miles Date: Sun, 31 Mar 2024 11:54:48 -0400 Subject: [PATCH 12/26] Remove unused dependencies --- Project.toml | 8 -------- 1 file changed, 8 deletions(-) diff --git a/Project.toml b/Project.toml index 29ffa9b..0cdfae0 100644 --- a/Project.toml +++ b/Project.toml @@ -5,18 +5,10 @@ version = "0.2.1" [deps] ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" -KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -Observers = "338f10d5-c7f1-4033-a7d1-f9dec39bcaa0" -Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" -TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" [compat] ITensors = "0.3.59" -KrylovKit = "0.6" -Observers = "0.2" -TimerOutputs = "0.5" julia = "1.6" [extras] From 204a765d77cbcf71ce22d53dac332e6dad6a718b Mon Sep 17 00:00:00 2001 From: Miles Date: Sun, 31 Mar 2024 11:56:47 -0400 Subject: [PATCH 13/26] Add compat entry for Reexport --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 0cdfae0..bfa7400 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69" [compat] ITensors = "0.3.59" +Reexport = "1" julia = "1.6" [extras] From 010a696bfaa7d75386e612a5bf00b1ad926178e1 Mon Sep 17 00:00:00 2001 From: Matt Fishman Date: Sun, 31 Mar 2024 17:11:43 -0400 Subject: [PATCH 14/26] Bump ITensors compat to v0.3.60 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index bfa7400..3d029d2 100644 --- a/Project.toml +++ b/Project.toml @@ -8,7 +8,7 @@ ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" [compat] -ITensors = "0.3.59" +ITensors = "0.3.60" Reexport = "1" julia = "1.6" From 8042972cc289347e049c0c9ef10261e74b1c5007 Mon Sep 17 00:00:00 2001 From: Miles Date: Sun, 31 Mar 2024 22:19:07 -0400 Subject: [PATCH 15/26] Remove example that was using code internals --- examples/05_tdvp_nonuniform_timesteps.jl | 52 --------------- examples/05_utils.jl | 85 ------------------------ 2 files changed, 137 deletions(-) delete mode 100644 examples/05_tdvp_nonuniform_timesteps.jl delete mode 100644 examples/05_utils.jl diff --git a/examples/05_tdvp_nonuniform_timesteps.jl b/examples/05_tdvp_nonuniform_timesteps.jl deleted file mode 100644 index 0a76206..0000000 --- a/examples/05_tdvp_nonuniform_timesteps.jl +++ /dev/null @@ -1,52 +0,0 @@ -using ITensors: MPO, MPS, OpSum, ProjMPO, siteinds -using KrylovKit: exponentiate -using Observers: observer - -include("05_utils.jl") - -function main() - function heisenberg(N) - os = OpSum() - for j in 1:(N - 1) - os += 0.5, "S+", j, "S-", j + 1 - os += 0.5, "S-", j, "S+", j + 1 - os += "Sz", j, "Sz", j + 1 - end - return os - end - - N = 10 - cutoff = 1e-12 - outputlevel = 1 - nsteps = 10 - time_steps = [n ≤ 2 ? -0.2im : -0.1im for n in 1:nsteps] - - obs = observer("times" => (; current_time) -> current_time, "psis" => (; psi) -> psi) - - s = siteinds("S=1/2", N; conserve_qns=true) - H = MPO(heisenberg(N), s) - - psi0 = MPS(s, n -> isodd(n) ? "Up" : "Dn") - psi = tdvp_nonuniform_timesteps( - ProjMPO(H), psi0; time_steps, cutoff, outputlevel, (step_observer!)=obs - ) - - times = obs.times - psis = obs.psis - - println("\nResults") - println("=======") - print("step = ", 0) - print(", time = ", zero(ComplexF64)) - print(", ⟨Sᶻ⟩ = ", round(expect(psi0, "Sz"; sites=N ÷ 2); digits=3)) - println() - for n in 1:length(times) - print("step = ", n) - print(", time = ", round(times[n]; digits=3)) - print(", ⟨Sᶻ⟩ = ", round(expect(psis[n], "Sz"; sites=N ÷ 2); digits=3)) - println() - end - return nothing -end - -main() diff --git a/examples/05_utils.jl b/examples/05_utils.jl deleted file mode 100644 index 89ff174..0000000 --- a/examples/05_utils.jl +++ /dev/null @@ -1,85 +0,0 @@ -using ITensors.ITensorMPS: ITensorMPS, MPS, maxlinkdim -using Observers: observer, update! -using Printf: @printf - -function tdvp_nonuniform_timesteps( - solver, - PH, - psi::MPS; - time_steps, - reverse_step=true, - time_start=0.0, - order=2, - (step_observer!)=observer(), - maxdim=ITensorMPS.default_maxdim(), - mindim=ITensorMPS.default_mindim(), - cutoff=ITensorMPS.default_cutoff(), - noise=ITensorMPS.default_noise(), - outputlevel=ITensorMPS.default_outputlevel(), - kwargs..., -) - nsweeps = length(time_steps) - maxdim, mindim, cutoff, noise = ITensorMPS.process_sweeps(; - nsweeps, maxdim, mindim, cutoff, noise - ) - tdvp_order = ITensorMPS.TDVPOrder(order, Base.Forward) - current_time = time_start - for sw in 1:nsweeps - sw_time = @elapsed begin - psi, PH, info = ITensorMPS.sweep_update( - tdvp_order, - solver, - PH, - time_steps[sw], - psi; - kwargs..., - current_time, - reverse_step, - sweep=sw, - maxdim=maxdim[sw], - mindim=mindim[sw], - cutoff=cutoff[sw], - noise=noise[sw], - ) - end - current_time += time_steps[sw] - update!(step_observer!; psi, sweep=sw, outputlevel, current_time) - if outputlevel ≥ 1 - print("After sweep ", sw, ":") - print(" maxlinkdim=", maxlinkdim(psi)) - @printf(" maxerr=%.2E", info.maxtruncerr) - print(" current_time=", round(current_time; digits=3)) - print(" time=", round(sw_time; digits=3)) - println() - flush(stdout) - end - end - return psi -end - -function tdvp_nonuniform_timesteps( - H, - psi::MPS; - ishermitian=ITensorMPS.default_ishermitian(), - issymmetric=ITensorMPS.default_issymmetric(), - solver_tol=ITensorMPS.default_solver_tol(exponentiate), - solver_krylovdim=ITensorMPS.default_solver_krylovdim(exponentiate), - solver_maxiter=ITensorMPS.default_solver_maxiter(exponentiate), - solver_outputlevel=ITensorMPS.default_solver_outputlevel(exponentiate), - kwargs..., -) - return tdvp_nonuniform_timesteps( - ITensorMPS.tdvp_solver( - exponentiate; - ishermitian, - issymmetric, - solver_tol, - solver_krylovdim, - solver_maxiter, - solver_outputlevel, - ), - H, - psi; - kwargs..., - ) -end From d7a587b15af661556478252d0841b70fda19737a Mon Sep 17 00:00:00 2001 From: Miles Date: Sun, 31 Mar 2024 22:26:39 -0400 Subject: [PATCH 16/26] Restore solver_utils.jl and associated exports --- src/ITensorTDVP.jl | 3 ++ src/solver_utils.jl | 69 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 72 insertions(+) create mode 100644 src/solver_utils.jl diff --git a/src/ITensorTDVP.jl b/src/ITensorTDVP.jl index 2047e64..f73661f 100644 --- a/src/ITensorTDVP.jl +++ b/src/ITensorTDVP.jl @@ -6,4 +6,7 @@ using Reexport: @reexport using ITensors.ITensorMPS: ITensorMPS dmrg(args...; kwargs...) = ITensorMPS.alternating_update_dmrg(args...; kwargs...) +include("solver_utils.jl") +export TimeDependentSum, to_vec + end diff --git a/src/solver_utils.jl b/src/solver_utils.jl new file mode 100644 index 0000000..6372e56 --- /dev/null +++ b/src/solver_utils.jl @@ -0,0 +1,69 @@ +using ITensors: ITensor, apply, array, inds, itensor, permute + +# Utilities for making it easier +# to define solvers (like ODE solvers) +# for TDVP + +""" + to_vec(x) +Transform `x` into a `Vector`, and return the vector, and a closure which inverts the +transformation. + +Modeled after FiniteDifferences.to_vec: + +https://github.com/JuliaDiff/FiniteDifferences.jl/blob/main/src/to_vec.jl +""" +to_vec(x) = error("Not implemented") + +function to_vec(x::ITensor) + function to_itensor(x_vec) + return itensor(x_vec, inds(x)) + end + return vec(array(x)), to_itensor +end + +# Represents a time-dependent sum of terms: +# +# H(t) = f[1](t) * H0[1] + f[2](t) * H0[2] + … +# +struct TimeDependentSum{S,T} + f::Vector{S} + H0::T +end +TimeDependentSum(f::Vector, H0::ProjMPOSum) = TimeDependentSum(f, ITensors.terms(H0)) +Base.length(H::TimeDependentSum) = length(H.f) + +function Base.:*(c::Number, H::TimeDependentSum) + return TimeDependentSum([t -> c * fₙ(t) for fₙ in H.f], H.H0) +end +Base.:*(H::TimeDependentSum, c::Number) = c * H + +# Calling a `TimeDependentOpSum` at a certain time like: +# +# H(t) +# +# Returns a `ScaledSum` at that time. +(H::TimeDependentSum)(t::Number) = ScaledSum([fₙ(t) for fₙ in H.f], H.H0) + +# Represents the sum of scaled terms: +# +# H = coefficient[1] * H[1] + coefficient * H[2] + … +# +struct ScaledSum{S,T} + coefficients::Vector{S} + H::T +end +Base.length(H::ScaledSum) = length(H.coefficients) + +# Apply the scaled sum of terms: +# +# H(ψ₀) = coefficient[1] * H[1](ψ₀) + coefficient[2] * H[2](ψ₀) + … +# +# onto ψ₀. +function (H::ScaledSum)(ψ₀) + ψ = ITensor(inds(ψ₀)) + for n in 1:length(H) + ψ += H.coefficients[n] * apply(H.H[n], ψ₀) + end + return permute(ψ, inds(ψ₀)) +end From 62d08ebc3adad2a5d191708c6fdabe77333cb3a4 Mon Sep 17 00:00:00 2001 From: Miles Date: Sun, 31 Mar 2024 23:50:25 -0400 Subject: [PATCH 17/26] Call new internal names for ITensorTDVP functions inside ITensorMPS --- Project.toml | 3 +-- src/ITensorTDVP.jl | 21 +++++++++++++++++---- src/solver_utils.jl | 2 +- test/test_examples.jl | 6 +----- 4 files changed, 20 insertions(+), 12 deletions(-) diff --git a/Project.toml b/Project.toml index 3d029d2..494f556 100644 --- a/Project.toml +++ b/Project.toml @@ -5,11 +5,10 @@ version = "0.2.1" [deps] ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" -Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" [compat] ITensors = "0.3.60" -Reexport = "1" julia = "1.6" [extras] diff --git a/src/ITensorTDVP.jl b/src/ITensorTDVP.jl index f73661f..1f20635 100644 --- a/src/ITensorTDVP.jl +++ b/src/ITensorTDVP.jl @@ -1,10 +1,23 @@ module ITensorTDVP -using Reexport: @reexport -@reexport using ITensors.ITensorMPS: tdvp, dmrg_x, to_vec, TimeDependentSum, linsolve - +using ITensors: ITensors, MPO +using ITensors.NDTensors: @Algorithm_str using ITensors.ITensorMPS: ITensorMPS -dmrg(args...; kwargs...) = ITensorMPS.alternating_update_dmrg(args...; kwargs...) +using KrylovKit: KrylovKit + +tdvp(args...; kwargs...) = ITensorMPS.itensortdvp_tdvp(args...; kwargs...) + +function KrylovKit.linsolve(A::MPO, args...; kwargs...) + return ITensorMPS.itensortdvp_linsolve(A, args...; kwargs...) +end + +dmrg(args...; kwargs...) = ITensorMPS.itensortdvp_dmrg(args...; kwargs...) + +dmrg_x(args...; kwargs...) = ITensorMPS.itensortdvp_dmrg_x(args...; kwargs...) + +function ITensors.contract(alg::Algorithm"fit", args...; kwargs...) + return ITensorMPS.itensortdvp_contract(alg, args...; kwargs...) +end include("solver_utils.jl") export TimeDependentSum, to_vec diff --git a/src/solver_utils.jl b/src/solver_utils.jl index 6372e56..ae2ab1b 100644 --- a/src/solver_utils.jl +++ b/src/solver_utils.jl @@ -1,4 +1,4 @@ -using ITensors: ITensor, apply, array, inds, itensor, permute +using ITensors: ITensor, apply, array, inds, itensor, permute, ProjMPOSum # Utilities for making it easier # to define solvers (like ODE solvers) diff --git a/test/test_examples.jl b/test/test_examples.jl index 025b888..88f0ac7 100644 --- a/test/test_examples.jl +++ b/test/test_examples.jl @@ -4,11 +4,7 @@ using ITensorTDVP: ITensorTDVP using Test: @testset @testset "Run examples" begin examples_files = [ - "01_tdvp.jl", - "02_dmrg-x.jl", - "03_tdvp_time_dependent.jl", - "04_tdvp_observers.jl", - "05_tdvp_nonuniform_timesteps.jl", + "01_tdvp.jl", "02_dmrg-x.jl", "03_tdvp_time_dependent.jl", "04_tdvp_observers.jl" ] examples_path = joinpath(pkgdir(ITensorTDVP), "examples") @testset "Running example file $f" for f in examples_files From 9358bf359a80952a8138eb834c4933b62d8afcf3 Mon Sep 17 00:00:00 2001 From: Matt Fishman Date: Mon, 1 Apr 2024 09:22:52 -0400 Subject: [PATCH 18/26] Bump to v0.3 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 494f556..b815b4e 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ITensorTDVP" uuid = "25707e16-a4db-4a07-99d9-4d67b7af0342" authors = ["Matthew Fishman and contributors"] -version = "0.2.1" +version = "0.3" [deps] ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" From 58ef5d7ab5c99766d99966d363abe214846d8910 Mon Sep 17 00:00:00 2001 From: Miles Date: Mon, 1 Apr 2024 13:11:02 -0400 Subject: [PATCH 19/26] More specific wrappers for ITensorTDVP interface --- src/ITensorTDVP.jl | 38 ++++++++++++++++++++++++++++++-------- 1 file changed, 30 insertions(+), 8 deletions(-) diff --git a/src/ITensorTDVP.jl b/src/ITensorTDVP.jl index 1f20635..d1dba32 100644 --- a/src/ITensorTDVP.jl +++ b/src/ITensorTDVP.jl @@ -1,22 +1,44 @@ module ITensorTDVP -using ITensors: ITensors, MPO +using ITensors: ITensors, MPO, MPS using ITensors.NDTensors: @Algorithm_str using ITensors.ITensorMPS: ITensorMPS using KrylovKit: KrylovKit -tdvp(args...; kwargs...) = ITensorMPS.itensortdvp_tdvp(args...; kwargs...) +function tdvp(H_generic, t::Number, psi0::MPS; kwargs...) + return ITensorMPS.itensortdvp_tdvp(H_generic, t, psi0; kwargs...) +end +function tdvp(t::Number, H_generic, psi0::MPS; kwargs...) + return ITensorMPS.itensortdvp_tdvp(t, H_generic, psi0; kwargs...) +end +function tdvp(H_generic, psi0::MPS, t::Number; kwargs...) + return ITensorMPS.itensortdvp_tdvp(H_generic, psi0, t; kwargs...) +end +function tdvp(solver, H::MPO, t::Number, psi0::MPS; kwargs...) + return ITensorMPS.itensortdvp_tdvp(solver, H, t, psi0; kwargs...) +end +function tdvp(solver, t::Number, H::MPO, psi0::MPS; kwargs...) + return ITensorMPS.itensortdvp_tdvp(solver, t, H, psi0; kwargs...) +end +function tdvp(solver, psi0::MPS, t::Number, H::MPO; kwargs...) + return ITensorMPS.itensortdvp_tdvp(solver, psi0, t, H; kwargs...) +end +function tdvp(solver, Hs::Vector{MPO}, t::Number, psi0::MPS; kwargs...) + return ITensorMPS.itensortdvp_tdvp(solver, Hs, t, psi0; kwargs...) +end -function KrylovKit.linsolve(A::MPO, args...; kwargs...) - return ITensorMPS.itensortdvp_linsolve(A, args...; kwargs...) +function KrylovKit.linsolve( + A::MPO, b::MPS, x0::MPS, a0::Number=false, a1::Number=true; kwargs... +) + return ITensorMPS.itensortdvp_linsolve(A, b, x0, a0, a1; kwargs...) end -dmrg(args...; kwargs...) = ITensorMPS.itensortdvp_dmrg(args...; kwargs...) +dmrg(H, psi0::MPS; kwargs...) = ITensorMPS.itensortdvp_dmrg(H, psi0; kwargs...) -dmrg_x(args...; kwargs...) = ITensorMPS.itensortdvp_dmrg_x(args...; kwargs...) +dmrg_x(H, psi0::MPS; kwargs...) = ITensorMPS.itensortdvp_dmrg_x(H, psi0; kwargs...) -function ITensors.contract(alg::Algorithm"fit", args...; kwargs...) - return ITensorMPS.itensortdvp_contract(alg, args...; kwargs...) +function ITensors.contract(alg::Algorithm"fit", A::MPO, psi::MPS; kwargs...) + return ITensorMPS.itensortdvp_contract(alg, A, psi; kwargs...) end include("solver_utils.jl") From 34e4b19b9858e36dc145bdd9fc34685031073a08 Mon Sep 17 00:00:00 2001 From: Miles Date: Mon, 1 Apr 2024 13:57:38 -0400 Subject: [PATCH 20/26] Update naming conventions --- src/ITensorTDVP.jl | 40 ++++++++++++++++++++++------------------ 1 file changed, 22 insertions(+), 18 deletions(-) diff --git a/src/ITensorTDVP.jl b/src/ITensorTDVP.jl index d1dba32..5f7d883 100644 --- a/src/ITensorTDVP.jl +++ b/src/ITensorTDVP.jl @@ -1,30 +1,30 @@ module ITensorTDVP -using ITensors: ITensors, MPO, MPS +using ITensors: ITensors +using ITensors.ITensorMPS: ITensorMPS, MPO, MPS using ITensors.NDTensors: @Algorithm_str -using ITensors.ITensorMPS: ITensorMPS using KrylovKit: KrylovKit -function tdvp(H_generic, t::Number, psi0::MPS; kwargs...) - return ITensorMPS.itensortdvp_tdvp(H_generic, t, psi0; kwargs...) +function tdvp(operator, t::Number, initial_state::MPS; kwargs...) + return ITensorMPS.itensortdvp_tdvp(operator, t, initial_state; kwargs...) end -function tdvp(t::Number, H_generic, psi0::MPS; kwargs...) - return ITensorMPS.itensortdvp_tdvp(t, H_generic, psi0; kwargs...) +function tdvp(t::Number, operator, initial_state::MPS; kwargs...) + return ITensorMPS.itensortdvp_tdvp(t, operator, initial_state; kwargs...) end -function tdvp(H_generic, psi0::MPS, t::Number; kwargs...) - return ITensorMPS.itensortdvp_tdvp(H_generic, psi0, t; kwargs...) +function tdvp(operator, initial_state::MPS, t::Number; kwargs...) + return ITensorMPS.itensortdvp_tdvp(operator, initial_state, t; kwargs...) end -function tdvp(solver, H::MPO, t::Number, psi0::MPS; kwargs...) - return ITensorMPS.itensortdvp_tdvp(solver, H, t, psi0; kwargs...) +function tdvp(solver, operator::MPO, t::Number, initial_state::MPS; kwargs...) + return ITensorMPS.itensortdvp_tdvp(solver, operator, t, initial_state; kwargs...) end -function tdvp(solver, t::Number, H::MPO, psi0::MPS; kwargs...) - return ITensorMPS.itensortdvp_tdvp(solver, t, H, psi0; kwargs...) +function tdvp(solver, t::Number, operator::MPO, initial_state::MPS; kwargs...) + return ITensorMPS.itensortdvp_tdvp(solver, t, operator, initial_state; kwargs...) end -function tdvp(solver, psi0::MPS, t::Number, H::MPO; kwargs...) - return ITensorMPS.itensortdvp_tdvp(solver, psi0, t, H; kwargs...) +function tdvp(solver, initial_state::MPS, t::Number, operator::MPO; kwargs...) + return ITensorMPS.itensortdvp_tdvp(solver, initial_state, t, operator; kwargs...) end -function tdvp(solver, Hs::Vector{MPO}, t::Number, psi0::MPS; kwargs...) - return ITensorMPS.itensortdvp_tdvp(solver, Hs, t, psi0; kwargs...) +function tdvp(solver, operators::Vector{MPO}, t::Number, initial_state::MPS; kwargs...) + return ITensorMPS.itensortdvp_tdvp(solver, operators, t, initial_state; kwargs...) end function KrylovKit.linsolve( @@ -33,9 +33,13 @@ function KrylovKit.linsolve( return ITensorMPS.itensortdvp_linsolve(A, b, x0, a0, a1; kwargs...) end -dmrg(H, psi0::MPS; kwargs...) = ITensorMPS.itensortdvp_dmrg(H, psi0; kwargs...) +function dmrg(operator, initial_state::MPS; kwargs...) + return ITensorMPS.itensortdvp_dmrg(operator, initial_state; kwargs...) +end -dmrg_x(H, psi0::MPS; kwargs...) = ITensorMPS.itensortdvp_dmrg_x(H, psi0; kwargs...) +function dmrg_x(operator, initial_state::MPS; kwargs...) + return ITensorMPS.itensortdvp_dmrg_x(operator, initial_state; kwargs...) +end function ITensors.contract(alg::Algorithm"fit", A::MPO, psi::MPS; kwargs...) return ITensorMPS.itensortdvp_contract(alg, A, psi; kwargs...) From e8251299e070ee505a18d39a4c036dfd32bda209 Mon Sep 17 00:00:00 2001 From: Miles Date: Tue, 2 Apr 2024 22:21:03 -0400 Subject: [PATCH 21/26] Update function argument names --- src/ITensorTDVP.jl | 40 ++++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/src/ITensorTDVP.jl b/src/ITensorTDVP.jl index 5f7d883..81ebd3f 100644 --- a/src/ITensorTDVP.jl +++ b/src/ITensorTDVP.jl @@ -5,26 +5,26 @@ using ITensors.ITensorMPS: ITensorMPS, MPO, MPS using ITensors.NDTensors: @Algorithm_str using KrylovKit: KrylovKit -function tdvp(operator, t::Number, initial_state::MPS; kwargs...) - return ITensorMPS.itensortdvp_tdvp(operator, t, initial_state; kwargs...) +function tdvp(operator, t::Number, init_state::MPS; kwargs...) + return ITensorMPS.itensortdvp_tdvp(operator, t, init_state; kwargs...) end -function tdvp(t::Number, operator, initial_state::MPS; kwargs...) - return ITensorMPS.itensortdvp_tdvp(t, operator, initial_state; kwargs...) +function tdvp(t::Number, operator, init_state::MPS; kwargs...) + return ITensorMPS.itensortdvp_tdvp(t, operator, init_state; kwargs...) end -function tdvp(operator, initial_state::MPS, t::Number; kwargs...) - return ITensorMPS.itensortdvp_tdvp(operator, initial_state, t; kwargs...) +function tdvp(operator, init_state::MPS, t::Number; kwargs...) + return ITensorMPS.itensortdvp_tdvp(operator, init_state, t; kwargs...) end -function tdvp(solver, operator::MPO, t::Number, initial_state::MPS; kwargs...) - return ITensorMPS.itensortdvp_tdvp(solver, operator, t, initial_state; kwargs...) +function tdvp(solver, operator::MPO, t::Number, init_state::MPS; kwargs...) + return ITensorMPS.itensortdvp_tdvp(solver, operator, t, init_state; kwargs...) end -function tdvp(solver, t::Number, operator::MPO, initial_state::MPS; kwargs...) - return ITensorMPS.itensortdvp_tdvp(solver, t, operator, initial_state; kwargs...) +function tdvp(solver, t::Number, operator::MPO, init_state::MPS; kwargs...) + return ITensorMPS.itensortdvp_tdvp(solver, t, operator, init_state; kwargs...) end -function tdvp(solver, initial_state::MPS, t::Number, operator::MPO; kwargs...) - return ITensorMPS.itensortdvp_tdvp(solver, initial_state, t, operator; kwargs...) +function tdvp(solver, init_state::MPS, t::Number, operator::MPO; kwargs...) + return ITensorMPS.itensortdvp_tdvp(solver, init_state, t, operator; kwargs...) end -function tdvp(solver, operators::Vector{MPO}, t::Number, initial_state::MPS; kwargs...) - return ITensorMPS.itensortdvp_tdvp(solver, operators, t, initial_state; kwargs...) +function tdvp(solver, operators::Vector{MPO}, t::Number, init_state::MPS; kwargs...) + return ITensorMPS.itensortdvp_tdvp(solver, operators, t, init_state; kwargs...) end function KrylovKit.linsolve( @@ -33,16 +33,16 @@ function KrylovKit.linsolve( return ITensorMPS.itensortdvp_linsolve(A, b, x0, a0, a1; kwargs...) end -function dmrg(operator, initial_state::MPS; kwargs...) - return ITensorMPS.itensortdvp_dmrg(operator, initial_state; kwargs...) +function dmrg(operator, init_state::MPS; kwargs...) + return ITensorMPS.itensortdvp_dmrg(operator, init_state; kwargs...) end -function dmrg_x(operator, initial_state::MPS; kwargs...) - return ITensorMPS.itensortdvp_dmrg_x(operator, initial_state; kwargs...) +function dmrg_x(operator, init_state::MPS; kwargs...) + return ITensorMPS.itensortdvp_dmrg_x(operator, init_state; kwargs...) end -function ITensors.contract(alg::Algorithm"fit", A::MPO, psi::MPS; kwargs...) - return ITensorMPS.itensortdvp_contract(alg, A, psi; kwargs...) +function ITensors.contract(alg::Algorithm"fit", tn1::MPO, tn2::MPS; kwargs...) + return ITensorMPS.itensortdvp_contract(alg, tn1, tn2; kwargs...) end include("solver_utils.jl") From 37777528773231069ee730cde2ecd302410c9132 Mon Sep 17 00:00:00 2001 From: Miles Date: Thu, 25 Apr 2024 17:41:52 -0400 Subject: [PATCH 22/26] Add Reexport as dependency --- Project.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Project.toml b/Project.toml index f979f9c..36019d3 100644 --- a/Project.toml +++ b/Project.toml @@ -6,10 +6,12 @@ version = "0.3" [deps] ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" [compat] ITensors = "0.3.58, 0.4" KrylovKit = "0.6, 0.7" +Reexport = "1" julia = "1.6" [extras] From 2e1f1ad7f576ca6a0e8376502d25c040ba5281c4 Mon Sep 17 00:00:00 2001 From: Miles Date: Thu, 25 Apr 2024 17:43:43 -0400 Subject: [PATCH 23/26] Change most functions to be reexported from ITensorMPS --- src/ITensorTDVP.jl | 43 +++---------------------------------------- 1 file changed, 3 insertions(+), 40 deletions(-) diff --git a/src/ITensorTDVP.jl b/src/ITensorTDVP.jl index 81ebd3f..43aed20 100644 --- a/src/ITensorTDVP.jl +++ b/src/ITensorTDVP.jl @@ -1,50 +1,13 @@ module ITensorTDVP -using ITensors: ITensors -using ITensors.ITensorMPS: ITensorMPS, MPO, MPS -using ITensors.NDTensors: @Algorithm_str -using KrylovKit: KrylovKit - -function tdvp(operator, t::Number, init_state::MPS; kwargs...) - return ITensorMPS.itensortdvp_tdvp(operator, t, init_state; kwargs...) -end -function tdvp(t::Number, operator, init_state::MPS; kwargs...) - return ITensorMPS.itensortdvp_tdvp(t, operator, init_state; kwargs...) -end -function tdvp(operator, init_state::MPS, t::Number; kwargs...) - return ITensorMPS.itensortdvp_tdvp(operator, init_state, t; kwargs...) -end -function tdvp(solver, operator::MPO, t::Number, init_state::MPS; kwargs...) - return ITensorMPS.itensortdvp_tdvp(solver, operator, t, init_state; kwargs...) -end -function tdvp(solver, t::Number, operator::MPO, init_state::MPS; kwargs...) - return ITensorMPS.itensortdvp_tdvp(solver, t, operator, init_state; kwargs...) -end -function tdvp(solver, init_state::MPS, t::Number, operator::MPO; kwargs...) - return ITensorMPS.itensortdvp_tdvp(solver, init_state, t, operator; kwargs...) -end -function tdvp(solver, operators::Vector{MPO}, t::Number, init_state::MPS; kwargs...) - return ITensorMPS.itensortdvp_tdvp(solver, operators, t, init_state; kwargs...) -end - -function KrylovKit.linsolve( - A::MPO, b::MPS, x0::MPS, a0::Number=false, a1::Number=true; kwargs... -) - return ITensorMPS.itensortdvp_linsolve(A, b, x0, a0, a1; kwargs...) -end +using Reexport: @reexport +@reexport using ITensors.ITensorMPS: contract, dmrg_x, tdvp +using ITensors.ITensorMPS: ITensorMPS, MPS function dmrg(operator, init_state::MPS; kwargs...) return ITensorMPS.itensortdvp_dmrg(operator, init_state; kwargs...) end -function dmrg_x(operator, init_state::MPS; kwargs...) - return ITensorMPS.itensortdvp_dmrg_x(operator, init_state; kwargs...) -end - -function ITensors.contract(alg::Algorithm"fit", tn1::MPO, tn2::MPS; kwargs...) - return ITensorMPS.itensortdvp_contract(alg, tn1, tn2; kwargs...) -end - include("solver_utils.jl") export TimeDependentSum, to_vec From f78b291bf0305795568e09dd8d5375924f686340 Mon Sep 17 00:00:00 2001 From: Miles Date: Thu, 25 Apr 2024 17:52:11 -0400 Subject: [PATCH 24/26] Reexport linsolve from KrylovKit --- src/ITensorTDVP.jl | 3 ++- test/test_linsolve.jl | 3 +-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/ITensorTDVP.jl b/src/ITensorTDVP.jl index 43aed20..9340400 100644 --- a/src/ITensorTDVP.jl +++ b/src/ITensorTDVP.jl @@ -2,8 +2,9 @@ module ITensorTDVP using Reexport: @reexport @reexport using ITensors.ITensorMPS: contract, dmrg_x, tdvp - +@reexport using KrylovKit: linsolve using ITensors.ITensorMPS: ITensorMPS, MPS + function dmrg(operator, init_state::MPS; kwargs...) return ITensorMPS.itensortdvp_dmrg(operator, init_state; kwargs...) end diff --git a/test/test_linsolve.jl b/test/test_linsolve.jl index 42d0902..870eb4d 100644 --- a/test/test_linsolve.jl +++ b/test/test_linsolve.jl @@ -1,7 +1,6 @@ @eval module $(gensym()) using ITensors: ITensors, MPO, OpSum, apply, randomMPS, siteinds -using ITensorTDVP: ITensorTDVP, dmrg -using KrylovKit: linsolve +using ITensorTDVP: ITensorTDVP, dmrg, linsolve using LinearAlgebra: norm using Test: @test, @testset using Random: Random From 0094d7eb713ba1853a9954901b76685e34706054 Mon Sep 17 00:00:00 2001 From: Miles Date: Sat, 27 Apr 2024 21:19:11 -0400 Subject: [PATCH 25/26] Add test of exports --- test/test_exports.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) create mode 100644 test/test_exports.jl diff --git a/test/test_exports.jl b/test/test_exports.jl new file mode 100644 index 0000000..2a3f526 --- /dev/null +++ b/test/test_exports.jl @@ -0,0 +1,10 @@ +@eval module $(gensym()) +using ITensorTDVP: ITensorTDVP +using Test: @test, @testset +@testset "Test exports" begin + itensortdvp_exports = [ + :ITensorTDVP, :TimeDependentSum, :dmrg_x, :linsolve, :tdvp, :to_vec + ] + @test issetequal(names(ITensorTDVP), itensortdvp_exports) +end +end From 1bcded0dbc13896ec9c2537aedfef73804a3bc38 Mon Sep 17 00:00:00 2001 From: Miles Date: Sat, 27 Apr 2024 21:46:01 -0400 Subject: [PATCH 26/26] Restore ITensorMPS wrappers and remove Reexport dependency --- Project.toml | 2 -- src/ITensorTDVP.jl | 48 ++++++++++++++++++++++++++++++++++++++----- test/test_linsolve.jl | 3 ++- 3 files changed, 45 insertions(+), 8 deletions(-) diff --git a/Project.toml b/Project.toml index 36019d3..f979f9c 100644 --- a/Project.toml +++ b/Project.toml @@ -6,12 +6,10 @@ version = "0.3" [deps] ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" -Reexport = "189a3867-3050-52da-a836-e630ba90ab69" [compat] ITensors = "0.3.58, 0.4" KrylovKit = "0.6, 0.7" -Reexport = "1" julia = "1.6" [extras] diff --git a/src/ITensorTDVP.jl b/src/ITensorTDVP.jl index 9340400..c6fc238 100644 --- a/src/ITensorTDVP.jl +++ b/src/ITensorTDVP.jl @@ -1,15 +1,53 @@ module ITensorTDVP -using Reexport: @reexport -@reexport using ITensors.ITensorMPS: contract, dmrg_x, tdvp -@reexport using KrylovKit: linsolve -using ITensors.ITensorMPS: ITensorMPS, MPS +export TimeDependentSum, dmrg_x, linsolve, tdvp, to_vec + +using ITensors: ITensors +using ITensors.ITensorMPS: ITensorMPS, MPO, MPS +using ITensors.NDTensors: @Algorithm_str +using KrylovKit: KrylovKit +function tdvp(operator, t::Number, init_state::MPS; kwargs...) + return ITensorMPS.itensortdvp_tdvp(operator, t, init_state; kwargs...) +end +function tdvp(t::Number, operator, init_state::MPS; kwargs...) + return ITensorMPS.itensortdvp_tdvp(t, operator, init_state; kwargs...) +end +function tdvp(operator, init_state::MPS, t::Number; kwargs...) + return ITensorMPS.itensortdvp_tdvp(operator, init_state, t; kwargs...) +end +function tdvp(solver, operator::MPO, t::Number, init_state::MPS; kwargs...) + return ITensorMPS.itensortdvp_tdvp(solver, operator, t, init_state; kwargs...) +end +function tdvp(solver, t::Number, operator::MPO, init_state::MPS; kwargs...) + return ITensorMPS.itensortdvp_tdvp(solver, t, operator, init_state; kwargs...) +end +function tdvp(solver, init_state::MPS, t::Number, operator::MPO; kwargs...) + return ITensorMPS.itensortdvp_tdvp(solver, init_state, t, operator; kwargs...) +end +function tdvp(solver, operators::Vector{MPO}, t::Number, init_state::MPS; kwargs...) + return ITensorMPS.itensortdvp_tdvp(solver, operators, t, init_state; kwargs...) +end + +function KrylovKit.linsolve( + A::MPO, b::MPS, x0::MPS, a0::Number=false, a1::Number=true; kwargs... +) + return ITensorMPS.itensortdvp_linsolve(A, b, x0, a0, a1; kwargs...) +end + +using ITensors.ITensorMPS: ITensorMPS, MPS function dmrg(operator, init_state::MPS; kwargs...) return ITensorMPS.itensortdvp_dmrg(operator, init_state; kwargs...) end +function dmrg_x(operator, init_state::MPS; kwargs...) + return ITensorMPS.itensortdvp_dmrg_x(operator, init_state; kwargs...) +end + +function ITensors.contract(alg::Algorithm"fit", tn1::MPO, tn2::MPS; kwargs...) + return ITensorMPS.itensortdvp_contract(alg, tn1, tn2; kwargs...) +end + include("solver_utils.jl") -export TimeDependentSum, to_vec end diff --git a/test/test_linsolve.jl b/test/test_linsolve.jl index 870eb4d..824ef8b 100644 --- a/test/test_linsolve.jl +++ b/test/test_linsolve.jl @@ -1,7 +1,8 @@ @eval module $(gensym()) using ITensors: ITensors, MPO, OpSum, apply, randomMPS, siteinds -using ITensorTDVP: ITensorTDVP, dmrg, linsolve +using ITensorTDVP: ITensorTDVP, dmrg using LinearAlgebra: norm +using KrylovKit: linsolve using Test: @test, @testset using Random: Random @testset "linsolve (eltype=$elt, conserve_qns=$conserve_qns)" for elt in (