Merge pull request #38 [Fix] Type compatibility

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]

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)

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
 

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,30 +118,32 @@ 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
         print("Calculating spectrum for bosonic systems...\n")
         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,22 +190,24 @@ 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
         print("Calculating spectrum for fermionic systems...\n")
         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
 

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...

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

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)