src/FiniteTLanczos.jl

@@ -32,11 +32,17 @@ export icgs,
 
 #ftlm_methods.jl
 export FTLMOptions, TraceEstimator,
-       full_ed,
+       full_ed_ftrace, full_ed_operator_average,
        simple_ftlm, 
        replaced_ftlm,
-       orthogonalized_ftlm
-
+       orthogonalized_ftlm,
+       evaluate_ftrace
+
+export full_ed_operator_average, 
+       simple_ftlm_operator_average, 
+       replaced_ftlm_operator_average,
+       orthogonalized_ftlm_operator_average,
+       evaluate_operator_average
 
 include("solver.jl")
 include("utilities.jl")

src/evaluate.jl

 
 """
-    evaluate(trace_estimator::TraceEstimator, args...)
+    evaluate_ftrace(trace_estimator::TraceEstimator, args...)
 
-Evaluate the trace of a matrix functions by a certain `trace_estimator`, where
-`evaluate(trace_estimator, M, expr)` is usd to estimate ``Tr[expr(M)]`` 
-and `evaluate(trace_estimator, M, Op; expr)` is used to estimate ``Tr[O expr(M)]``,
-where `M` and `O` are matrix-like instances and `M` is required to be hermitian.
+Evaluate the trace of a matrix functions ``Tr[expr(M)]`` by a certain 
+`trace_estimator`, where `M` is a matrix-like instances and is required 
+to be hermitian.
 
 The `trace_estimator` is chosen from the Finite Temperature Lanczos Methods(FTLM) family,
 which includes four methods:
-- `full_ed`: The exact diagonalization method
-- `simple_ftlm`: The naive FTLM, which is also konwn as the Huchinson's method
-- `replaced_ftlm`: Replaced FTLM.
-- `orthogonalized_ftlm`: orthogonalized FTLM
+- `full_ed_ftrace`: The exact diagonalization method
+- `simple_ftlm_ftrace`: The naive FTLM, which is also konwn as the Huchinson's method
+- `replaced_ftlm_ftrace`: Replaced FTLM.
+- `orthogonalized_ftlm_ftrace`: orthogonalized FTLM
 
 
 Note:
@@ -20,4 +19,27 @@ Note:
    `size`, `length`, `getindex`, `iterate` and matrix-vector products are required.
 2. `M` is not explicity symmetrized in the following methods. 
 """
-function evaluate end
+function evaluate_ftrace end
+
+
+"""
+    evaluate_operator_average(trace_estimator::TraceEstimator, args...)
+
+Evaluate ``Tr[O expr(M)]``, the thermal average for a operator ``O``, 
+whose density operator is ``expr(M)``, where `M` and `O` are matrix-like 
+instances and `M` is required to be hermitian.
+
+The `trace_estimator` is chosen from the Finite Temperature Lanczos Methods(FTLM) family,
+which includes four methods:
+- `full_ed_operator_average`: The exact diagonalization method
+- `simple_ftlm_operator_average`: The naive FTLM, which is also konwn as the Huchinson's method
+- `replaced_ftlm_operator_average`: Replaced FTLM.
+- `orthogonalized_ftlm_operator_average`: orthogonalized FTLM
+
+
+Note:
+1. A matrix-like instance `M` behaves as a matrix, whose `eltype`, 
+   `size`, `length`, `getindex`, `iterate` and matrix-vector products are required.
+2. `M` is not explicity symmetrized in the following methods. 
+"""
+function evaluate_operator_average end

src/ftlm_methods.jl

@@ -26,19 +26,19 @@ Struct to store the options when appling FTLM.
 end
 
 """ 
-    full_ed(M, expr::Function [; options::FTLMOptions = FTLMOptions()])
-    full_ed(M, expr::Tuple{Function,N}[; options::FTLMOptions = FTLMOptions()])
+    full_ed_ftrace(M, expr::Function [; options::FTLMOptions = FTLMOptions()])
+    full_ed_ftrace(M, expr::Tuple{Function,N}[; options::FTLMOptions = FTLMOptions()])
 
 Calculate ``Tr[expr(M)]`` by full exact diagonalization method.
 """
-function full_ed(M, expr::Function; options::FTLMOptions=FTLMOptions())
+function full_ed_ftrace(M, expr::Function; options::FTLMOptions=FTLMOptions())
     solver = solver_function(options.processor)
     M = solver(M)
     vals, _ = eigensolver(M)
     return reduce(+, map(expr , vals))
 end
 
-function full_ed(M, expr::Vararg{Function,N}; options::FTLMOptions=FTLMOptions()) where {N}
+function full_ed_ftrace(M, expr::Vararg{Function,N}; options::FTLMOptions=FTLMOptions()) where {N}
     solver = solver_function(options.processor)
     M = solver(M)
     vals, _ = eigensolver(M)
@@ -53,17 +53,18 @@ function _sum_expr_w1_w2(expr::Function, vals, w1, w2)
 end
 
 """
-    simple_ftlm(M, expr::Function [; options::FTLMOptions = FTLMOptions()])
-    simple_ftlm(M, expr::Tuple{Function,N}[; options::FTLMOptions = FTLMOptions()])
+    simple_ftlm_ftrace(M, expr::Function [; options::FTLMOptions = FTLMOptions()])
+    simple_ftlm_ftrace(M, expr::Tuple{Function,N}[; options::FTLMOptions = FTLMOptions()])
 
 Calculate ``Tr[expr(M)]`` by FTLM, which is also konwn as the Huchinson's method.
 """
-function simple_ftlm(M, expr::Function; options::FTLMOptions = FTLMOptions())
+function simple_ftlm_ftrace(M, expr::Function; options::FTLMOptions = FTLMOptions())
     @unpack distr, Nr, Nk, processor = options
     solver = solver_function(processor)
     M = solver(M)
 
     Ns = size(M, 1)
+    Nk = min(Nk, Ns)
     res = 0.
     for r = 1: Nr
         init_vector = random_unit_vector(Ns, distr) |> solver
@@ -74,7 +75,7 @@ function simple_ftlm(M, expr::Function; options::FTLMOptions = FTLMOptions())
     return res * factor
 end
 
-function simple_ftlm(M, expr::Vararg{Function,N}; options::FTLMOptions = FTLMOptions()) where {N}
+function simple_ftlm_ftrace(M, expr::Vararg{Function,N}; options::FTLMOptions = FTLMOptions()) where {N}
     @unpack distr, Nr, Nk, processor = options
     solver = solver_function(processor)
     M = solver(M)
@@ -92,12 +93,12 @@ end
 
 
 """
-    replaced_ftlm(M, expr::Function [; options::FTLMOptions = FTLMOptions()])
-    replaced_ftlm(M, expr::Tuple{Function,N}[; options::FTLMOptions = FTLMOptions()])
+    replaced_ftlm_ftrace(M, expr::Function [; options::FTLMOptions = FTLMOptions()])
+    replaced_ftlm_ftrace(M, expr::Tuple{Function,N}[; options::FTLMOptions = FTLMOptions()])
 
 Calculate ``Tr[expr(M)]`` by replaced FTLM.
 """
-function replaced_ftlm(M, expr::Function; options::FTLMOptions = FTLMOptions())
+function replaced_ftlm_ftrace(M, expr::Function; options::FTLMOptions = FTLMOptions())
     @unpack distr, Nr, Nk, Ne, processor, which = options
     solver = solver_function(processor)
     M = solver(M)
@@ -126,7 +127,7 @@ function replaced_ftlm(M, expr::Function; options::FTLMOptions = FTLMOptions())
     return res * factor
 end
 
-function replaced_ftlm(M, expr::Vararg{Function,N}; options::FTLMOptions = FTLMOptions()) where {N}
+function replaced_ftlm_ftrace(M, expr::Vararg{Function,N}; options::FTLMOptions = FTLMOptions()) where {N}
     @unpack distr, Nr, Nk, Ne, processor, which = options
     solver = solver_function(processor)
     M = solver(M)
@@ -157,12 +158,12 @@ end
 
 
 """
-    orthogonalized_ftlm(M, expr::Function [; options::FTLMOptions = FTLMOptions()])
-    orthogonalized_ftlm(M, expr::Tuple{Function,N}[; options::FTLMOptions = FTLMOptions()])
+    orthogonalized_ftlm_ftrace(M, expr::Function [; options::FTLMOptions = FTLMOptions()])
+    orthogonalized_ftlm_ftrace(M, expr::Tuple{Function,N}[; options::FTLMOptions = FTLMOptions()])
 
 Calculate ``Tr[expr(M)]`` by orthogonalized FTLM.
 """
-function orthogonalized_ftlm(M, expr::Function; options::FTLMOptions=FTLMOptions())
+function orthogonalized_ftlm_ftrace(M, expr::Function; options::FTLMOptions=FTLMOptions())
     @unpack distr, Nr, Nk, Ne, processor, which = options
     solver = solver_function(processor)
     M = solver(M)
@@ -187,7 +188,7 @@ function orthogonalized_ftlm(M, expr::Function; options::FTLMOptions=FTLMOptions
     return res1 + res2 * factor
 end
 
-function orthogonalized_ftlm(M, expr::Vararg{Function,N}; options::FTLMOptions=FTLMOptions()) where {N}
+function orthogonalized_ftlm_ftrace(M, expr::Vararg{Function,N}; options::FTLMOptions=FTLMOptions()) where {N}
     @unpack distr, Nr, Nk, Ne, processor, which = options
     solver = solver_function(processor)
     M = solver(M)
@@ -223,7 +224,7 @@ Composite struct of a method function and its keyword arguments
     options::To
 end
 
-function evaluate(method::TraceEstimator, M, expr::Vararg{Function,N}) where {N}
+function evaluate_ftrace(method::TraceEstimator, M, expr::Vararg{Function,N}) where {N}
     @unpack estimator, options = method
     return estimator(M, expr...; options)
 end

src/solver.jl

@@ -31,8 +31,22 @@ end
 
 cpu_solver(M::AbstractArray{T,N}) where {T,N} = convert(Array{T,N}, M)
 cpu_solver(expr, M::AbstractArray{T,N}) where {T,N} = expr(convert(Array{T,N}, M))
+function cpu_solver(expr, M::AbstractArray...)
+    M1 = Vector{Array}(undef, length(M))
+    for i in eachindex(M)
+        M1[i] = convert(Array, M[i])
+    end
+    return expr(M1...)
+end
+
 
 gpu_solver(M::AbstractArray{T,N}) where {T,N} = convert(CuArray{T,N}, M)
 gpu_solver(expr, M::AbstractArray{T,N}) where {T,N} = expr(convert(CuArray{T,N}, M))
-
+function gpu_solver(expr, M::AbstractArray) 
+    M1 = Vector{CuArray}(undef, length(M))
+    for i in eachindex(M)
+        M1[i] = convert(CuArray, M[i])
+    end
+    return expr(M1...)
+end
 

src/thermal_average.jl

 function _sum_diag_vecs_Op(Λ, vecs, Op)
+    vecs2 = Op * vecs
+    ave = ein"ik, ki -> i"(conj(vecs), vecs2)
+    return reduce(+, map(*, Λ, ave))
+end
+
+function _sum_diag_vecs_Op(Λ, vecs, Op::AbstractArray)
     ave = ein"ik, kl, li -> i"(conj(vecs), Op, vecs)
     return reduce(+, map(*, Λ, ave))
 end
 
-function _sum_diag_w1_v0_Op_vecs(Λ, w1,v0, Op, vecs)
+
+function _sum_diag_w1_v0_Op_vecs(Λ, w1, v0, Op)
+    vecs2 = Op * vecs
+    w2 = ein"i,ij -> j"(conj(v0), vecs2)
+    return reduce(+,map(*, Λ, w1, w2))
+end
+
+function _sum_diag_w1_v0_Op_vecs(Λ, w1, v0, Op::AbstractArray, vecs)
     w2 = ein"i,ik,kj -> j"(conj(v0), Op, vecs)
     return reduce(+,map(*, Λ, w1, w2))
 end
 
+
+
 """
-    full_ed(M, Op::AbstractArray; expr::Function[, options::FTLMOptions=FTLMOptions()])
+    full_ed_operator_average(M, Op; expr::Function[, options::FTLMOptions=FTLMOptions()])
 
 Calculate ``Tr[Op expr(M)]`` by full exact diagonalization method.
 """
-function full_ed(M, Op::AbstractArray; expr::Function, options::FTLMOptions=FTLMOptions())
+function full_ed_operator_average(M, Op; expr::Function, options::FTLMOptions=FTLMOptions())
     solver = solver_function(options.processor)
     M = solver(M)
 
@@ -26,11 +41,11 @@ end
 
 
 """
-    simple_ftlm(M, Op::AbstractArray; expr::Function[, options::FTLMOptions=FTLMOptions()])
+    simple_ftlm_operator_average(M, Op; expr::Function[, options::FTLMOptions=FTLMOptions()])
 
 Calculate ``Tr[Op expr(M)]`` by simple FTLM.
 """
-function simple_FTLM(M, Op::AbstractArray; expr::Function,options::FTLMOptions = FTLMOptions())
+function simple_FTLM_operator_average(M, Op; expr::Function,options::FTLMOptions = FTLMOptions())
     @unpack distr, Nr, Nk, processor = options
     solver = solver_function(processor)
     M = solver(M)
@@ -53,11 +68,11 @@ end
 
 
 """
-    replaced_ftlm(M, Op::AbstractArray; expr::Function[, options::FTLMOptions=FTLMOptions()])
+    replaced_ftlm_operator_average(M, Op; expr::Function[, options::FTLMOptions=FTLMOptions()])
 
 Calculate ``Tr[Op expr(M)]`` by replaced FTLM.
 """
-function replaced_ftlm(M, Op::AbstractArray; expr::Function,options::FTLMOptions = FTLMOptions())
+function replaced_ftlm_operator_average(M, Op; expr::Function,options::FTLMOptions = FTLMOptions())
     @unpack distr, Nr, Nk, Ne, processor, which = options
     solver = solver_function(processor)
     M = solver(M)
@@ -93,11 +108,11 @@ end
 
 
 """
-    orthogonalized_ftlm(M, Op::AbstractArray; expr::Function[, options::FTLMOptions=FTLMOptions())
+    orthogonalized_ftlm_operator_average(M, Op; expr::Function[, options::FTLMOptions=FTLMOptions())
 
 Calculate ``Tr[Op expr(M)]`` by orthogonalized FTLM.
 """
-function orthogonalized_ftlm(M, Op::AbstractArray; expr::Function,options::FTLMOptions = FTLMOptions())
+function orthogonalized_ftlm_operator_average(M, Op; expr::Function,options::FTLMOptions = FTLMOptions())
     @unpack distr, Nr, Nk, Ne, processor, which = options
     solver = solver_function(processor)
     M = solver(M)
@@ -130,18 +145,18 @@ function orthogonalized_ftlm(M, Op::AbstractArray; expr::Function,options::FTLMO
 end
 
 
-function evaluate(method::TraceEstimator, M, Op::AbstractArray; expr::Function)
+function evaluate_operator_average(method::TraceEstimator, M, Op; expr::Function)
     @unpack estimator, options = method
     return estimator(M, Op; expr, options)
 end
 
 
 """
-    thermal_average(M, β, Op::AbstractArray; trace_estimator::TraceEstimator)
+    thermal_average(M, β, Op; trace_estimator::TraceEstimator)
 
 Calculate the thermal average ``⟨Op⟩ = \\frac{1}{Z} Tr(Op e^{-βH})`` of operator ``Op``.
 """
-function thermal_average(M, β, Op::AbstractArray; trace_estimator::TraceEstimator)
+function thermal_average(M, β, Op; trace_estimator::TraceEstimator)
     @unpack processor = trace_estimator.options
     processor == CPU ? x0 = rand(size(M,1)) : x0 = CUDA.rand(Float64, size(M,1))
     e0, _, _ = eigsolve(M, x0, 1, :SR, ishermitian = true)

src/thermaldynamic.jl

@@ -20,7 +20,7 @@ function partitian(M, β; trace_estimator::TraceEstimator)
     e0, _, _ = eigsolve(M, x0, 1, :SR, ishermitian = true)
     e1 = e0[1]
     exprZ = e -> _boltzmann_factor(e, e1, β)
-    return evaluate(estimator, M, exprZ)[1]
+    return evaluate_ftrace(estimator, M, exprZ)[1]
 end
 
 
@@ -38,7 +38,7 @@ function free_energy(M, β; trace_estimator::TraceEstimator)
     e0, _, _ = eigsolve(M, x0, 1, :SR, ishermitian = true)
     e1 = e0[1]
     exprZ = e -> _boltzmann_factor(e, e1, β)
-    res = evaluate(estimator, M, exprZ)[1]
+    res = evaluate_ftrace(estimator, M, exprZ)[1]
     return -log(res)/β + e1
 end
 
@@ -57,7 +57,7 @@ function energy(M, β; trace_estimator::TraceEstimator)
     e1 = e0[1]
     exprZ = e -> _boltzmann_factor(e, e1, β)
     exprE = e -> e * _boltzmann_factor(e, e1, β)
-    Z, E = evaluate(estimator, M, exprZ, exprE)
+    Z, E = evaluate_ftrace(estimator, M, exprZ, exprE)
     return E/Z
 end
 
@@ -78,7 +78,7 @@ function specific_heat(M, β; trace_estimator::TraceEstimator)
     exprE = e -> e * _boltzmann_factor(e, e1, β)
     exprCv = e -> e^2 * _boltzmann_factor(e, e1, β)
    
-    Z, E, Cv = evaluate(estimator, M, exprZ, exprE, exprCv)
+    Z, E, Cv = evaluate_ftrace(estimator, M, exprZ, exprE, exprCv)
     E = E/Z
     return β^2*(Cv/Z - E^2)
 end
@@ -98,7 +98,7 @@ function entropy(M, β; trace_estimator::TraceEstimator)
     e1 = e0[1]
     exprZ = e -> _boltzmann_factor(e, e1, β)
     exprE = e -> e * _boltzmann_factor(e, e1, β)
-    Z, E = evaluate(estimator, M, exprZ, exprE)
+    Z, E = evaluate_ftrace(estimator, M, exprZ, exprE)
     F = -log(Z)/β + e1
     E = E/Z
     return β*(E - F)