[ITensorMPS] Examples (#1408)

Author: mtfishman

Project.toml

@@ -1,7 +1,7 @@
 name = "ITensors"
 uuid = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 authors = ["Matthew Fishman <[email protected]>", "Miles Stoudenmire <[email protected]>"]
-version = "0.5.1"
+version = "0.5.2"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

docs/src/CodeTiming.md

@@ -7,93 +7,4 @@ If you are concerned about the performance of your code, a good place to start i
 
 ## Timing and benchmarking
 
-Julia has many nice timing tools available. Tools like [@time](https://docs.julialang.org/en/v1/base/base/#Base.@time) can be used to measure the time of specific lines of code. For microbenchmarking, we recommend the [BenchmarkTools](https://github.com/JuliaCI/BenchmarkTools.jl) package.
-
-## ProfileView.jl
-
-Julia has a built-in profiler in the standard module [Profile](https://docs.julialang.org/en/v1/manual/profile/). You can use it as follows:
-```julia
-julis> using ITensors
-
-julia> using Profile
-
-julia> @profile include(joinpath(pkgdir(ITensors), "examples", "dmrg", "1d_heisenberg_conserve_spin.jl"))
-After sweep 1 energy=-137.995732867390 maxlinkdim=10 maxerr=1.93E-02 time=0.862
-After sweep 2 energy=-138.801057557054 maxlinkdim=20 maxerr=3.37E-05 time=1.126
-After sweep 3 energy=-138.940075984826 maxlinkdim=91 maxerr=9.99E-11 time=1.880
-After sweep 4 energy=-138.940086063995 maxlinkdim=99 maxerr=1.00E-10 time=3.033
-After sweep 5 energy=-138.940086076330 maxlinkdim=95 maxerr=9.97E-11 time=2.824
-Final energy = -138.940086076330
-
-julia> Profile.print()
-Overhead ╎ [+additional indent] Count File:Line; Function
-=========================================================
-   ╎6576 @Base/client.jl:485; _start()
-   ╎ 6576 @Base/client.jl:302; exec_options(opts::Base.JLOptions)
-   ╎  6576 @Base/client.jl:372; run_main_repl(interactive::Bool, quiet::Bool, b...
-   ╎   6576 @Base/essentials.jl:706; invokelatest
-   ╎    6576 @Base/essentials.jl:708; #invokelatest#2
-   ╎     6576 @Base/client.jl:387; (::Base.var"#874#876"{Bool, Bool, Bool})(REPL...
-   ╎    ╎ 6576 .../stdlib/v1.6/REPL/src/REPL.jl:305; run_repl(repl::REPL.AbstractREPL, consumer::Any)
-   ╎    ╎  6576 ...stdlib/v1.6/REPL/src/REPL.jl:317; run_repl(repl::REPL.AbstractREPL, consumer::...
-[...]
-```
-We don't actually recommend printing it since the output is extremely long and difficult to read. Instead, it is best to use a visualization tool to visualize the profile, so you can see at a glance what parts of your code are taking up the most time.
-
-For a quick and easy way of seeing how fast different parts of your code are, we highly recommend [ProfileView.jl](https://github.com/timholy/ProfileView.jl). All you have to do is load `ProfileView` and append a simple macro `@profview`:
-```julia
-julia> using ITensors
-
-julia> using ProfileView
-
-julia> @profview include(joinpath(pkgdir(ITensors), "examples", "dmrg", "1d_heisenberg_conserve_spin.jl"));
-After sweep 1 energy=-137.995732867390 maxlinkdim=10 maxerr=1.93E-02 time=0.977
-After sweep 2 energy=-138.801057557054 maxlinkdim=20 maxerr=3.37E-05 time=1.252
-After sweep 3 energy=-138.940075984826 maxlinkdim=91 maxerr=9.99E-11 time=2.263
-After sweep 4 energy=-138.940086063995 maxlinkdim=99 maxerr=1.00E-10 time=2.938
-After sweep 5 energy=-138.940086076330 maxlinkdim=95 maxerr=9.97E-11 time=2.988
-Final energy = -138.940086076330
-```
-A window will pop up with a "flame graph", where the width of a bar corresponds to the percentage of time the function call took, and as you go up in the graph you see timings for more and more nested functions. The graph will still look quite complicated for this case, but at larger bond dimensions you should see that certain functions like matrix multiplications and decompositions will start to dominate.
-
-## TimerOutputs.jl
-
-For more focused timing of specific parts of your code, we recommend the package [TimerOutputs](https://github.com/KristofferC/TimerOutputs.jl). This can help measure lines of code that are called repeatedly, where you want to know the accumulated time. We have some timers defined internally in ITensors.jl for the `dmrg` function which can be called as follows:
-```julia
-julia> using ITensors
-
-julia> ITensors.TimerOutputs.enable_debug_timings(ITensors)
-timeit_debug_enabled (generic function with 1 method)
-
-julia> ITensors.TimerOutputs.reset_timer!(ITensors.NDTensors.timer)
- ──────────────────────────────────────────────────────────────────
-                           Time                   Allocations
-                   ──────────────────────   ───────────────────────
- Tot / % measured:     59.4μs / 0.00%              992B / 0.00%
-
- Section   ncalls     time   %tot     avg     alloc   %tot      avg
- ──────────────────────────────────────────────────────────────────
- ──────────────────────────────────────────────────────────────────
-
-julia> include("1d_heisenberg_conserve_spin.jl")
-After sweep 1 energy=-137.995732867390 maxlinkdim=10 maxerr=1.93E-02 time=0.597
-After sweep 2 energy=-138.801057557054 maxlinkdim=20 maxerr=3.37E-05 time=0.798
-After sweep 3 energy=-138.940075984826 maxlinkdim=91 maxerr=9.99E-11 time=1.285
-After sweep 4 energy=-138.940086063995 maxlinkdim=99 maxerr=1.00E-10 time=1.878
-After sweep 5 energy=-138.940086076330 maxlinkdim=95 maxerr=9.97E-11 time=1.936
-Final energy = -138.940086076330
-
-julia> ITensors.NDTensors.timer
- ────────────────────────────────────────────────────────────────────────────────
-                                         Time                   Allocations
-                                 ──────────────────────   ───────────────────────
-        Tot / % measured:             14.4s / 45.0%           9.59GiB / 95.1%
-
- Section                 ncalls     time   %tot     avg     alloc   %tot      avg
- ────────────────────────────────────────────────────────────────────────────────
- dmrg: eigsolve             990    4.72s  72.7%  4.76ms   7.27GiB  79.7%  7.52MiB
- dmrg: replacebond!         990    1.31s  20.2%  1.32ms   1.12GiB  12.2%  1.15MiB
- dmrg: position!            990    360ms  5.55%   364μs    571MiB  6.11%   590KiB
- dmrg: psi[b]*psi[b+1]      990    106ms  1.63%   107μs    182MiB  1.95%   188KiB
- ────────────────────────────────────────────────────────────────────────────────
-```
+Julia has many nice timing tools available. Tools like [@time](https://docs.julialang.org/en/v1/base/base/#Base.@time) and [TimerOutputs](https://github.com/KristofferC/TimerOutputs.jl) can be used to measure the time of specific lines of code. For microbenchmarking, we recommend the [BenchmarkTools](https://github.com/JuliaCI/BenchmarkTools.jl) package. For profiling your code, see the Julia documentation on [profiling](https://docs.julialang.org/en/v1/manual/profile/).

docs/src/DMRGObserver.md

@@ -6,11 +6,6 @@ for DMRG calculations such as measuring custom
 local observables at each step and stopping DMRG
 early if certain energy convergence conditions are met.
 
-In addition to the example code below, more detailed
-example code showing sample usage of DMRGObserver is included
-in the ITensor source, in the file `1d_ising_with_observer.jl`
-under the folder `examples/dmrg`.
-
 ## Sample Usage
 
 In the following example, we have already made a Hamiltonian MPO `H`

docs/src/Multithreading.md

@@ -150,144 +150,6 @@ Threaded contract:
 C_contract ≈ C_threaded_contract = true
 ```
 
-Here is a full example of making use of block sparse multithreading to speed up
-a DMRG calculation:
-```julia
-using ITensors
-using LinearAlgebra
-using Random
-
-include(joinpath(ITensors.examples_dir(), "src", "electronk.jl"))
-include(joinpath(ITensors.examples_dir(), "src", "hubbard.jl"))
-
-function main(; Nx::Int=6, Ny::Int=3, U::Float64=4.0, t::Float64=1.0,
-                maxdim::Int=3000, conserve_ky=true,
-                nsweeps=10, blas_num_threads=1, strided_num_threads=1,
-                threaded_blocksparse=true, outputlevel=1,
-                seed=1234)
-  Random.seed!(seed)
-  ITensors.Strided.set_num_threads(strided_num_threads)
-  BLAS.set_num_threads(blas_num_threads)
-  ITensors.enable_threaded_blocksparse(threaded_blocksparse)
-
-  if outputlevel > 0
-    @show Threads.nthreads()
-    @show Sys.CPU_THREADS
-    @show BLAS.get_num_threads()
-    @show ITensors.Strided.get_num_threads()
-    @show ITensors.using_threaded_blocksparse()
-    println()
-  end
-
-  N = Nx * Ny
-
-  maxdim = min.([100, 200, 400, 800, 2000, 3000, maxdim], maxdim)
-  cutoff = 1E-6
-  noise = [1E-6, 1E-7, 1E-8, 0.0]
-
-  sites = siteinds("ElecK", N; conserve_qns=true,
-                   conserve_ky, modulus_ky=Ny)
-
-  hubbard_ops = hubbard(Nx, Ny, t, U, ky=true)
-  H = MPO(hubbard_ops, sites)
-
-  # Number of structural nonzero elements in a bulk
-  # Hamiltonian MPO tensor
-  if outputlevel > 0
-    @show nnz(H[end÷2])
-    @show nnzblocks(H[end÷2])
-  end
-
-  # Create starting state with checkerboard
-  # pattern
-  state = map(CartesianIndices((Ny, Nx))) do I
-    return iseven(I[1]) ⊻ iseven(I[2]) ? "↓" : "↑"
-  end
-  display(state)
-
-  psi0 = randomMPS(sites, state; linkdims=10)
-
-  energy, psi = @time dmrg(H, psi0; nsweeps, maxdim, cutoff, noise, outputlevel)
-
-  if outputlevel > 0
-    @show Nx, Ny
-    @show t, U
-    @show flux(psi)
-    @show maxlinkdim(psi)
-    @show energy
-  end
-  return nothing
-end
-
-println("Without threaded block sparse:\n")
-main(; nsweeps=10, threaded_blocksparse=false)
-println()
-println("With threaded block sparse:\n")
-main(; nsweeps=10, threaded_blocksparse=true)
-println()
-```
-which when run on a laptop with 6 cores gives (after running once
-with 1 or 2 sweeps to trigger compilation):
-```
-Without threaded block sparse:
-
-Threads.nthreads() = 5
-Sys.CPU_THREADS = 6
-BLAS.get_num_threads() = 1
-ITensors.Strided.get_num_threads() = 1
-ITensors.using_threaded_blocksparse() = false
-
-splitblocks = true
-nnz(H[end ÷ 2]) = 67
-nnzblocks(H[end ÷ 2]) = 67
-After sweep 1 energy=-5.861157015737 maxlinkdim=78 time=0.633
-After sweep 2 energy=-12.397725587986 maxlinkdim=200 time=6.980
-After sweep 3 energy=-13.472412477115 maxlinkdim=400 time=14.257
-After sweep 4 energy=-13.627475223585 maxlinkdim=800 time=9.801
-After sweep 5 energy=-13.693628527487 maxlinkdim=2000 time=15.343
-After sweep 6 energy=-13.711928225391 maxlinkdim=3000 time=24.260
-After sweep 7 energy=-13.715575192226 maxlinkdim=3000 time=25.752
-After sweep 8 energy=-13.716394378223 maxlinkdim=3000 time=25.907
-After sweep 9 energy=-13.716535094341 maxlinkdim=3000 time=24.748
-After sweep 10 energy=-13.716556326664 maxlinkdim=3000 time=24.185
-171.903248 seconds (575.56 M allocations: 207.370 GiB, 9.37% gc time)
-(Nx, Ny) = (6, 3)
-(t, U) = (1.0, 4.0)
-flux(psi) = QN(("Ky",0,3),("Nf",18,-1),("Sz",0))
-maxlinkdim(psi) = 3000
-energy = -13.716556326663678
-
-With threaded block sparse:
-
-Threads.nthreads() = 5
-Sys.CPU_THREADS = 6
-BLAS.get_num_threads() = 1
-ITensors.Strided.get_num_threads() = 1
-ITensors.using_threaded_blocksparse() = true
-
-splitblocks = true
-nnz(H[end ÷ 2]) = 67
-nnzblocks(H[end ÷ 2]) = 67
-After sweep 1 energy=-5.861157015735 maxlinkdim=78 time=1.117
-After sweep 2 energy=-12.397725588011 maxlinkdim=200 time=6.587
-After sweep 3 energy=-13.472412477124 maxlinkdim=400 time=12.094
-After sweep 4 energy=-13.627475223588 maxlinkdim=800 time=8.760
-After sweep 5 energy=-13.693628527488 maxlinkdim=2000 time=12.447
-After sweep 6 energy=-13.711928225390 maxlinkdim=3000 time=17.401
-After sweep 7 energy=-13.715575192226 maxlinkdim=3000 time=18.863
-After sweep 8 energy=-13.716394378223 maxlinkdim=3000 time=19.258
-After sweep 9 energy=-13.716535094341 maxlinkdim=3000 time=19.627
-After sweep 10 energy=-13.716556326664 maxlinkdim=3000 time=18.446
-134.604491 seconds (1.69 G allocations: 300.213 GiB, 18.02% gc time)
-(Nx, Ny) = (6, 3)
-(t, U) = (1.0, 4.0)
-flux(psi) = QN(("Ky",0,3),("Nf",18,-1),("Sz",0))
-maxlinkdim(psi) = 3000
-energy = -13.71655632666368
-```
-Speedups should be greater for problems with tensors that are more
-sparse or have more blocks (for example larger system sizes).
-
 In addition, we plan to add more threading to other parts of the
 code beyond contraction (such as SVD) and improve composibility with
 other forms of threading like BLAS and Strided, so stay tuned!

docs/src/index.md

@@ -99,57 +99,6 @@ and associated "Codebase Release" for the version you have used. The current one
 }
 ```
 
-## Full Example Codes
-
-The ITensors.jl package contains a directory of examples, which we
-will continue to add to. You can find them online [here](https://github.com/ITensor/ITensors.jl/tree/main/examples).
-Additionally, once you have installed ITensors.jl you can find a local version
-of the examples in the directory `ITensors.examples_dir()`, and you can run them
-as follows from the Julia REPL:
-```julia
-julia> using ITensors
-
-julia> cd(ITensors.examples_dir())
-
-julia> readdir()
-7-element Array{String,1}:
- "basic_ops"
- "ctmrg"
- "dmrg"
- "gate_evolution"
- "krylov_methods"
- "src"
- "trg"
-
-julia> cd("dmrg")
-
-julia> readdir()
-8-element Array{String,1}:
- "1d_heisenberg.jl"
- "1d_heisenberg_conserve_spin.jl"
- "1d_hubbard_extended.jl"
- "1d_ising_with_observer.jl"
- "2d_heisenberg_conserve_spin.jl"
- "2d_hubbard_conserve_momentum.jl"
- "2d_hubbard_conserve_particles.jl"
- "input_files"
-
-julia> include("1d_heisenberg.jl")
-After sweep 1 energy=-138.837988775764 maxlinkdim=10 time=13.760
-After sweep 2 energy=-138.937408365962 maxlinkdim=20 time=0.249
-After sweep 3 energy=-138.940084788852 maxlinkdim=100 time=1.867
-After sweep 4 energy=-138.940086091070 maxlinkdim=100 time=3.824
-After sweep 5 energy=-138.940086113999 maxlinkdim=122 time=4.527
-Final energy = -138.940086113999
-
-julia> pwd()
-"[path_to_package_installation]/ITensors/examples/dmrg"
-```
-You can use your favorite text editor to view these examples. If you
-would like to modify them, either copy them into your own directory,
-or checkout ITensors.jl in development mode using the instructions
-in [Developing ITensors.jl](@ref).
-
 ## ITensor Code Samples
 
 ### Basic Overview
@@ -368,3 +317,5 @@ After sweep 4 energy=-138.940086009318 maxlinkdim=100 maxerr=1.05E-10 time=11.64
 After sweep 5 energy=-138.940086058840 maxlinkdim=96 maxerr=1.00E-10 time=12.771
 Final energy = -138.94008605883985
 ```
+You can find more examples of running `dmrg` and related algorithms [here](https://github.com/ITensor/ITensors.jl/tree/main/src/lib/ITensorMPS/examples).
+