diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl index a1dfa39c..0703adaa 100644 --- a/src/solvers/eigsolve.jl +++ b/src/solvers/eigsolve.jl @@ -1,38 +1,39 @@ -function eigsolve_solver(; - solver_which_eigenvalue=:SR, #TODO: settle on pattern to pass solver kwargs +function eigsolve_updater( + init; + psi_ref!, + PH_ref!, + outputlevel, + which_sweep, + region_updates, + which_region_update, + region_kwargs, + updater_kwargs, + ) + + default_updater_kwargs=(; + solver_which_eigenvalue=:SR, ishermitian=true, - solver_tol=1e-14, - solver_krylovdim=3, - solver_maxiter=1, - solver_verbosity=0, -) - function solver( - init; - psi_ref!, - PH_ref!, - normalize, - region, - sweep_regions, - sweep_step, - sweep_kwargs..., - # slurp solver_kwargs? #TODO: homogenize how the solver kwargs are passed + tol=1e-14, + krylovdim=3, + maxiter=1, + outputlevel=0, + eager=false, ) - howmany = 1 - which = solver_which_eigenvalue - vals, vecs, info = eigsolve( - PH_ref![], + updater_kwargs=merge(default_updater_kwargs,updater_kwargs) #last collection has precedence + howmany=1 + which=updater_kwargs[:solver_which_eigenvalue] + vals, vecs, info = KrylovKit.eigsolve(PH_ref![], init, howmany, which; - ishermitian, - tol=solver_tol, - krylovdim=solver_krylovdim, - maxiter=solver_maxiter, - verbosity=solver_verbosity, + ishermitian=updater_kwargs[:ishermitian], + tol=updater_kwargs[:tol], + krylovdim=updater_kwargs[:krylovdim], + maxiter=updater_kwargs[:maxiter], + verbosity=updater_kwargs[:outputlevel], + eager=updater_kwargs[:eager], ) - phi = vecs[1] - return phi, (; solver_info=info, energies=vals) - end - return solver + return vecs[1], (; info, energies=vals) end + diff --git a/src/solvers/exponentiate.jl b/src/solvers/exponentiate.jl index e70edc80..cf08c4d0 100644 --- a/src/solvers/exponentiate.jl +++ b/src/solvers/exponentiate.jl @@ -3,10 +3,10 @@ psi_ref!, PH_ref!, outputlevel, - which_sweep, # keep, change name + which_sweep, region_updates, which_region_update, - region_kwargs, # region_kwargs (timestep for solver) + region_kwargs, updater_kwargs, ) default_updater_kwargs=(; diff --git a/src/treetensornetworks/solvers/dmrg.jl b/src/treetensornetworks/solvers/dmrg.jl index 8f24fd63..be4c3039 100644 --- a/src/treetensornetworks/solvers/dmrg.jl +++ b/src/treetensornetworks/solvers/dmrg.jl @@ -1,34 +1,59 @@ """ Overload of `ITensors.dmrg`. """ + + +function dmrg_sweep( + nsite::Int, + graph::AbstractGraph; + root_vertex=default_root_vertex(graph), +) + return tdvp_sweep(2,nsite,Inf,graph;root_vertex,reverse_step=false) +end + + +function default_sweep_regions(nsite, graph::AbstractGraph; kwargs...) ###move this to a different file, algorithmic level idea + return vcat( + [ + half_sweep( + direction(half), + graph, + make_region; + nsite, + region_args=(; substep=half), + kwargs..., + ) for half in 1:2 + ]..., + ) +end + function dmrg( + updater, H, init::AbstractTTN; - solver_which_eigenvalue=:SR, - ishermitian=true, - solver_tol=1e-14, - solver_krylovdim=3, - solver_maxiter=1, - solver_verbosity=0, + nsweeps, #it makes sense to require this to be defined + nsite=2, + (sweep_observer!)=observer(), + root_vertex=default_root_vertex(init), + updater_kwargs=NamedTuple(;), kwargs..., -) - return alternating_update( - eigsolve_solver(; - solver_which_eigenvalue, - ishermitian, - solver_tol, - solver_krylovdim, - solver_maxiter, - solver_verbosity, - ), - H, - init; - reverse_step=false, - kwargs..., ) + region_updates = dmrg_sweep(nsite,init; root_vertex) + + psi = alternating_update( + updater, H, init; nsweeps, sweep_observer!, region_updates, updater_kwargs, kwargs... + ) + return psi +end + + +function dmrg(H, init::AbstractTTN; updater=eigsolve_updater, kwargs...) + return dmrg(updater, H, init; kwargs...) end """ Overload of `KrylovKit.eigsolve`. """ eigsolve(H, init::AbstractTTN; kwargs...) = dmrg(H, init; kwargs...) + + diff --git a/test/test_treetensornetworks/test_solvers/test_dmrg.jl b/test/test_treetensornetworks/test_solvers/test_dmrg.jl index 939b4fee..1e523b28 100644 --- a/test/test_treetensornetworks/test_solvers/test_dmrg.jl +++ b/test/test_treetensornetworks/test_solvers/test_dmrg.jl @@ -30,21 +30,23 @@ using Observers psi_mps = MPS([psi[v] for v in 1:nv(psi)]) e2, psi2 = dmrg(H_mpo, psi_mps; nsweeps, maxdim, outputlevel=0) - psi = dmrg(H, psi; nsweeps, maxdim, cutoff, nsite, solver_krylovdim=3, solver_maxiter=1) + psi = dmrg(H, psi; nsweeps, maxdim, cutoff, nsite, updater_kwargs=(;krylovdim=3, maxiter=1)) @test inner(psi', H, psi) ≈ inner(psi2', H_mpo, psi2) # Alias for `ITensorNetworks.dmrg` psi = eigsolve( - H, psi; nsweeps, maxdim, cutoff, nsite, solver_krylovdim=3, solver_maxiter=1 + H, psi; nsweeps, maxdim, cutoff, nsite, updater_kwargs=(;krylovdim=3, maxiter=1) ) @test inner(psi', H, psi) ≈ inner(psi2', H_mpo, psi2) - # Test custom sweep regions + # Test custom sweep regions #BROKEN, ToDo: Make proper custom sweep regions for test + #= orig_E = inner(psi', H, psi) sweep_regions = [[1], [2], [3], [3], [2], [1]] psi = dmrg(H, psi; nsweeps, maxdim, cutoff, sweep_regions) new_E = inner(psi', H, psi) @test new_E ≈ orig_E + =# end @testset "Observers" begin @@ -67,20 +69,20 @@ end # # Make observers # - sweep(; sweep, kw...) = sweep + sweep(; which_sweep, kw...) = which_sweep sweep_observer! = observer(sweep) - region(; region, kw...) = region + region(; which_region_update, region_updates, kw...) = first(region_updates[which_region_update]) energy(; energies, kw...) = energies[1] - step_observer! = observer(region, sweep, energy) + region_observer! = observer(region, sweep, energy) - psi = dmrg(H, psi; nsweeps, maxdim, cutoff, sweep_observer!, step_observer!) + psi = dmrg(H, psi; nsweeps, maxdim, cutoff, sweep_observer!, region_observer!) # # Test out certain values # - @test step_observer![9, :region] == [2, 1] - @test step_observer![30, :energy] < -4.25 + @test region_observer![9, :region] == [2, 1] + @test region_observer![30, :energy] < -4.25 end @testset "Regression test: Arrays of Parameters" begin @@ -126,7 +128,7 @@ end sweeps = Sweeps(nsweeps) # number of sweeps is 5 maxdim!(sweeps, 10, 20, 40, 100) # gradually increase states kept cutoff!(sweeps, cutoff) - psi = dmrg(H, psi; nsweeps, maxdim, cutoff, nsite, solver_krylovdim=3, solver_maxiter=1) + psi = dmrg(H, psi; nsweeps, maxdim, cutoff, nsite, updater_kwargs=(;krylovdim=3, maxiter=1)) # Compare to `ITensors.MPO` version of `dmrg` linear_order = [4, 1, 2, 5, 3, 6]