diff --git a/docs/src/extensions/CUDA.md b/docs/src/extensions/CUDA.md index d9e2cb8e..3ee68c85 100644 --- a/docs/src/extensions/CUDA.md +++ b/docs/src/extensions/CUDA.md @@ -40,8 +40,8 @@ kT = 0.5 N = 5 tier = 3 -tlist = 0f0:1f-1:10f0 # same as 0:0.1:10 but in the type of `Float32` -ωlist = -10f0:1f0:10f0 # same as -10:1:10 but in the type of `Float32` +tlist = 0:0.1:10 +ωlist = -10:1:10 σm = [0 1; 0 0] σz = [1 0; 0 -1] diff --git a/src/HeomBase.jl b/src/HeomBase.jl index 1a341f0f..b82b6f06 100644 --- a/src/HeomBase.jl +++ b/src/HeomBase.jl @@ -21,6 +21,23 @@ function HandleMatrixType(M::AbstractMatrix, dim::Int=0, MatrixName::String="") end end +function _HandleFloatType(ElType::Type{T}, V::StepRangeLen) where T <: Number + if real(ElType) == Float32 + return StepRangeLen(Float32(V.ref), Float32(V.step), Int32(V.len), Int64(V.offset)) + else + return StepRangeLen(Float64(V.ref), Float64(V.step), Int64(V.len), Int64(V.offset)) + end +end + +function _HandleFloatType(ElType::Type{T}, V::Any) where T <: Number + FType = real(ElType) + if eltype(V) == FType + return V + else + convert.(FType, V) + end +end + function _get_pkg_version(pkg_name::String) D = Pkg.dependencies() for uuid in keys(D) diff --git a/src/HierarchicalEOM.jl b/src/HierarchicalEOM.jl index 5d59ee5e..1e53fe13 100644 --- a/src/HierarchicalEOM.jl +++ b/src/HierarchicalEOM.jl @@ -54,7 +54,7 @@ module HierarchicalEOM import SparseArrays: sparse, spzeros, sparsevec, reshape, SparseVector, SparseMatrixCSC, AbstractSparseMatrix import ProgressMeter: Progress, next! import FastExpm: fastExpm - import ..HeomBase: PROGBAR_OPTIONS, HandleMatrixType + import ..HeomBase: PROGBAR_OPTIONS, HandleMatrixType, _HandleFloatType # for solving time evolution import SciMLOperators: MatrixOperator @@ -85,7 +85,7 @@ module HierarchicalEOM import SparseArrays: sparse, sparsevec, SparseMatrixCSC import LinearSolve: LinearProblem, init, solve!, UMFPACKFactorization import ProgressMeter: Progress, next! - import ..HeomBase: PROGBAR_OPTIONS, HandleMatrixType + import ..HeomBase: PROGBAR_OPTIONS, HandleMatrixType, _HandleFloatType export spectrum diff --git a/src/Spectrum.jl b/src/Spectrum.jl index 9cf83c11..501f478a 100644 --- a/src/Spectrum.jl +++ b/src/Spectrum.jl @@ -1,5 +1,5 @@ @doc raw""" - spectrum(M, ρ, op, ω_list; solver, verbose, filename, SOLVEROptions...) + spectrum(M, ρ, op, ωlist; solver, verbose, filename, SOLVEROptions...) Calculate spectrum for the system. # To calculate spectrum for bosonic systems (usually known as power spectrum): @@ -22,7 +22,7 @@ remember to set the parameters: - `M::AbstractHEOMLSMatrix` : the matrix given from HEOM model. - `ρ` : the system density matrix or the auxiliary density operators. - `op` : The annihilation operator acting on the system. -- `ω_list::AbstractVector` : the specific frequency points to solve. +- `ωlist::AbstractVector` : the specific frequency points to solve. - `solver` : solver in package `LinearSolve.jl`. Default to `UMFPACKFactorization()`. - `verbose::Bool` : To display verbose output and progress bar during the process or not. Defaults to `true`. - `filename::String` : If filename was specified, the value of spectrum for each ω will be saved into the file "filename.txt" during the solving process. @@ -31,13 +31,13 @@ remember to set the parameters: For more details about solvers and extra options, please refer to [`LinearSolve.jl`](http://linearsolve.sciml.ai/stable/) # Returns -- `spec::AbstractVector` : the spectrum list corresponds to the specified `ω_list` +- `spec::AbstractVector` : the spectrum list corresponds to the specified `ωlist` """ function spectrum( M::AbstractHEOMLSMatrix, ρ, op, - ω_list::AbstractVector; + ωlist::AbstractVector; solver=UMFPACKFactorization(), verbose::Bool = true, filename::String = "", @@ -71,14 +71,14 @@ function spectrum( # check parity and calculate spectrum if (typeof(M.parity) == OddParity) - return _density_of_states(M, ados_vec, _op, ω_list; + return _density_of_states(M, ados_vec, _op, ωlist; solver = solver, verbose = verbose, filename = filename, SOLVEROptions... ) else - return _power_spectrum(M, ados_vec, _op, ω_list; + return _power_spectrum(M, ados_vec, _op, ωlist; solver = solver, verbose = verbose, filename = filename, @@ -91,7 +91,7 @@ end M::AbstractHEOMLSMatrix, ados_vec::AbstractVector, op, - ω_list::AbstractVector; + ωlist::AbstractVector; solver=UMFPACKFactorization(), verbose::Bool = true, filename::String = "", @@ -118,7 +118,9 @@ end a_dagger = kron(I_heom, spre(op')) X = _HandleVectorType(typeof(M.data), a_normal * ados_vec) - Length = length(ω_list) + i = convert(eltype(M), 1im) + ωList = _HandleFloatType(eltype(M), ωlist) + Length = length(ωList) Sω = Vector{Float64}(undef, Length) if verbose @@ -126,22 +128,22 @@ end flush(stdout) prog = Progress(Length; desc="Progress : ", PROGBAR_OPTIONS...) end - Iω = 1im * ω_list[1] * I_total + Iω = i * ωList[1] * I_total cache = init(LinearProblem(M.data - Iω, X), solver, SOLVEROptions...) sol = solve!(cache) - @inbounds for (i, ω) in enumerate(ω_list) - if i > 1 - Iω = 1im * ω * I_total + @inbounds for (j, ω) in enumerate(ωList) + if j > 1 + Iω = i * ω * I_total cache.A = M.data - Iω sol = solve!(cache) end # trace over the Hilbert space of system (expectation value) - Sω[i] = -1 * real(I_dual_vec * (a_dagger * _HandleVectorType(sol.u, false))[1:(M.sup_dim)]) + Sω[j] = -1 * real(I_dual_vec * (a_dagger * _HandleVectorType(sol.u, false))[1:(M.sup_dim)]) if SAVE open(FILENAME, "a") do file - write(file, "$(Sω[i]),\n") + write(file, "$(Sω[j]),\n") end end @@ -160,7 +162,7 @@ end M::AbstractHEOMLSMatrix, ados_vec::AbstractVector, op, - ω_list::AbstractVector; + ωlist::AbstractVector; solver=UMFPACKFactorization(), verbose::Bool = true, filename::String = "", @@ -188,7 +190,9 @@ end X_m = _HandleVectorType(typeof(M.data), d_normal * ados_vec) X_p = _HandleVectorType(typeof(M.data), d_dagger * ados_vec) - Length = length(ω_list) + i = convert(eltype(M), 1im) + ωList = _HandleFloatType(eltype(M), ωlist) + Length = length(ωList) Aω = Vector{Float64}(undef, Length) if verbose @@ -196,14 +200,14 @@ end flush(stdout) prog = Progress(Length; desc="Progress : ", PROGBAR_OPTIONS...) end - Iω = 1im * ω_list[1] * I_total + Iω = i * ωList[1] * I_total cache_m = init(LinearProblem(M.data - Iω, X_m), solver, SOLVEROptions...) cache_p = init(LinearProblem(M.data + Iω, X_p), solver, SOLVEROptions...) sol_m = solve!(cache_m) sol_p = solve!(cache_p) - @inbounds for (i, ω) in enumerate(ω_list) - if i > 1 - Iω = 1im * ω * I_total + @inbounds for (j, ω) in enumerate(ωList) + if j > 1 + Iω = i * ω * I_total cache_m.A = M.data - Iω sol_m = solve!(cache_m) @@ -215,14 +219,14 @@ end Cω_p = d_normal * _HandleVectorType(sol_p.u, false) # trace over the Hilbert space of system (expectation value) - Aω[i] = -1 * ( + Aω[j] = -1 * ( real(I_dual_vec * Cω_p[1:(M.sup_dim)]) + real(I_dual_vec * Cω_m[1:(M.sup_dim)]) ) if SAVE open(FILENAME, "a") do file - write(file, "$(Aω[i]),\n") + write(file, "$(Aω[j]),\n") end end diff --git a/src/SteadyState.jl b/src/SteadyState.jl index 29d2821c..78450219 100644 --- a/src/SteadyState.jl +++ b/src/SteadyState.jl @@ -151,7 +151,7 @@ For more details about solvers and extra options, please refer to [`Differential end sol = solve( SteadyStateProblem(prob), - DynamicSS(solver; abstol = abstol, reltol = reltol); + DynamicSS(solver; abstol = _HandleFloatType(eltype(M), abstol), reltol = _HandleFloatType(eltype(M), reltol)); maxiters = maxiters, save_everystep = save_everystep, SOLVEROptions... diff --git a/src/evolution.jl b/src/evolution.jl index 335800de..95ac99c2 100644 --- a/src/evolution.jl +++ b/src/evolution.jl @@ -259,23 +259,24 @@ For more details about solvers and extra options, please refer to [`Differential end end + Tlist = _HandleFloatType(eltype(M), tlist) ADOs_list::Vector{ADOs} = [ados] if SAVE jldopen(FILENAME, "a") do file - file[string(tlist[1])] = ados + file[string(Tlist[1])] = ados end end # problem: dρ/dt = L * ρ(0) L = MatrixOperator(M.data) - prob = ODEProblem(L, _HandleVectorType(typeof(M.data), ados.data), (tlist[1], tlist[end])) + prob = ODEProblem(L, _HandleVectorType(typeof(M.data), ados.data), (Tlist[1], Tlist[end])) # setup integrator integrator = init( prob, solver; - reltol = reltol, - abstol = abstol, + reltol = _HandleFloatType(eltype(M), reltol), + abstol = _HandleFloatType(eltype(M), abstol), maxiters = maxiters, save_everystep = save_everystep, SOLVEROptions... @@ -285,10 +286,10 @@ For more details about solvers and extra options, please refer to [`Differential if verbose print("Solving time evolution for auxiliary density operators...\n") flush(stdout) - prog = Progress(length(tlist); start=1, desc="Progress : ", PROGBAR_OPTIONS...) + prog = Progress(length(Tlist); start=1, desc="Progress : ", PROGBAR_OPTIONS...) end idx = 1 - dt_list = diff(tlist) + dt_list = diff(Tlist) for dt in dt_list idx += 1 step!(integrator, dt, true) @@ -299,7 +300,7 @@ For more details about solvers and extra options, please refer to [`Differential if SAVE jldopen(FILENAME, "a") do file - file[string(tlist[idx])] = ados + file[string(Tlist[idx])] = ados end end if verbose @@ -432,28 +433,29 @@ For more details about solvers and extra options, please refer to [`Differential end end + Tlist = _HandleFloatType(eltype(M), tlist) ADOs_list::Vector{ADOs} = [ados] if SAVE jldopen(FILENAME, "a") do file - file[string(tlist[1])] = ados + file[string(Tlist[1])] = ados end end - Ht = H(param, tlist[1]) - _Ht = HandleMatrixType(Ht, M.dim, "H (Hamiltonian) at t=$(tlist[1])") + Ht = H(param, Tlist[1]) + _Ht = HandleMatrixType(Ht, M.dim, "H (Hamiltonian) at t=$(Tlist[1])") Lt = kron(sparse(I, M.N, M.N), - 1im * (spre(_Ht) - spost(_Ht))) L = MatrixOperator(M.data + Lt, update_func! = _update_L!) # problem: dρ/dt = L(t) * ρ(0) ## M.dim will check whether the returned time-dependent Hamiltonian has the correct dimension - prob = ODEProblem(L, _HandleVectorType(typeof(M.data), ados.data), (tlist[1], tlist[end]), (M, H, param)) + prob = ODEProblem(L, _HandleVectorType(typeof(M.data), ados.data), (Tlist[1], Tlist[end]), (M, H, param)) # setup integrator integrator = init( prob, solver; - reltol = reltol, - abstol = abstol, + reltol = _HandleFloatType(eltype(M), reltol), + abstol = _HandleFloatType(eltype(M), abstol), maxiters = maxiters, save_everystep = save_everystep, SOLVEROptions... @@ -463,10 +465,10 @@ For more details about solvers and extra options, please refer to [`Differential if verbose print("Solving time evolution for auxiliary density operators with time-dependent Hamiltonian...\n") flush(stdout) - prog = Progress(length(tlist); start=1, desc="Progress : ", PROGBAR_OPTIONS...) + prog = Progress(length(Tlist); start=1, desc="Progress : ", PROGBAR_OPTIONS...) end idx = 1 - dt_list = diff(tlist) + dt_list = diff(Tlist) for dt in dt_list idx += 1 step!(integrator, dt, true) @@ -477,7 +479,7 @@ For more details about solvers and extra options, please refer to [`Differential if SAVE jldopen(FILENAME, "a") do file - file[string(tlist[idx])] = ados + file[string(Tlist[idx])] = ados end end if verbose diff --git a/test/CUDA/CUDAExt.jl b/test/CUDA/CUDAExt.jl index 41acd0f0..38b19c1c 100644 --- a/test/CUDA/CUDAExt.jl +++ b/test/CUDA/CUDAExt.jl @@ -39,28 +39,28 @@ Fbath = Fermion_Lorentz_Pade(Qf, λ, μ, W, kT, N) ## Schrodinger HEOMLS L_cpu = M_S(Hsys; verbose=false) L_gpu = cu(L_cpu) -ados_cpu = evolution(L_cpu, ρ0, [0f0, 10f0]; verbose=false) -ados_gpu = evolution(L_gpu, ρ0, [0f0, 10f0]; verbose=false) +ados_cpu = evolution(L_cpu, ρ0, [0, 10]; verbose=false) +ados_gpu = evolution(L_gpu, ρ0, [0, 10]; verbose=false) @test _is_Matrix_approx(getRho(ados_cpu[end]), getRho(ados_cpu[end])) ## Boson HEOMLS L_cpu = M_Boson(Hsys, tier, Bbath; verbose=false) L_gpu = cu(L_cpu) -ados_cpu = evolution(L_cpu, ρ0, [0f0, 10f0]; verbose=false) -ados_gpu = evolution(L_gpu, ρ0, [0f0, 10f0]; verbose=false) +ados_cpu = evolution(L_cpu, ρ0, [0, 10]; verbose=false) +ados_gpu = evolution(L_gpu, ρ0, [0, 10]; verbose=false) @test _is_Matrix_approx(getRho(ados_cpu[end]), getRho(ados_cpu[end])) ## Boson Fermion HEOMLS L_cpu = M_Fermion(Hsys, tier, Fbath; verbose=false) L_gpu = cu(L_cpu) -ados_cpu = evolution(L_cpu, ρ0, [0f0, 10f0]; verbose=false) -ados_gpu = evolution(L_gpu, ρ0, [0f0, 10f0]; verbose=false) +ados_cpu = evolution(L_cpu, ρ0, [0, 10]; verbose=false) +ados_gpu = evolution(L_gpu, ρ0, [0, 10]; verbose=false) @test _is_Matrix_approx(getRho(ados_cpu[end]), getRho(ados_cpu[end])) ## Boson Fermion HEOMLS L_cpu = M_Boson_Fermion(Hsys, tier, tier, Bbath, Fbath; verbose=false) L_gpu = cu(L_cpu) -tlist = 0f0:1f0:10f0 +tlist = 0:1:10 ados_cpu = evolution(L_cpu, ρ0, tlist; verbose=false) ados_gpu = evolution(L_gpu, ρ0, tlist; verbose=false) for i in 1:length(tlist) @@ -87,7 +87,7 @@ bath_dn = Fermion_Lorentz_Pade(d_dn, Γ, μ, W, kT, N) bath_list = [bath_up, bath_dn] ## solve density of states -ωlist = -10:1:10 +ωlist = -5:0.5:5 L_cpu = M_Fermion(Hsys, tier, bath_list; verbose=false) L_odd_cpu = M_Fermion(Hsys, tier, bath_list, ODD; verbose=false) L_odd_gpu = cu(L_odd_cpu)